[e40e335] | 1 | ! ****************************************************************
|
---|
| 2 | ! Trial version implementing the semi-local conformational update
|
---|
| 3 | ! BGS (Biased Gaussian Steps). This file presently contains the
|
---|
[32289cd] | 4 | ! functions initlund, bgsposs and bgs.
|
---|
[e40e335] | 5 | !
|
---|
| 6 | ! Copyright 2007 Frank Eisenmenger, U.H.E. Hansmann,
|
---|
| 7 | ! Jan H. Meinke, Sandipan Mohanty
|
---|
| 8 | ! ****************************************************************
|
---|
| 9 |
|
---|
| 10 |
|
---|
[32289cd] | 11 | ! Checks if it is possible to perform a BGS update starting at the
|
---|
[e40e335] | 12 | ! variable indexed ipos. Calls: none.
|
---|
| 13 | logical function bgsposs(ips)
|
---|
| 14 | include 'INCL.H'
|
---|
| 15 | include 'incl_lund.h'
|
---|
[32289cd] | 16 | integer jv, ips, iaa, nursvr, nnonfx, i
|
---|
| 17 |
|
---|
[e40e335] | 18 | logical ians
|
---|
| 19 |
|
---|
| 20 | jv=idvr(ips)
|
---|
| 21 | iaa=nursvr(jv)
|
---|
| 22 | ians=.true.
|
---|
| 23 | ! print *,'evaluating bgs possibility for ',ips,nmvr(jv)
|
---|
[32289cd] | 24 | if (nmvr(jv).ne.'phi') then
|
---|
[e40e335] | 25 | ! print *,'bgs not possible because variable name is ',nmvr(jv)
|
---|
| 26 | ians=.false.
|
---|
[32289cd] | 27 | else if (iaa.gt.(irsml2(mlvr(jv))-3)) then
|
---|
[e40e335] | 28 | ! print *,'bgs impossible, residue too close to end'
|
---|
| 29 | ! print *,'iaa = ',iaa,' end = ',irsml2(mlvr(jv))
|
---|
| 30 | ians=.false.
|
---|
[32289cd] | 31 | else
|
---|
| 32 | nnonfx=0
|
---|
[e40e335] | 33 | do i=iaa,iaa+3
|
---|
[32289cd] | 34 | if (iphi(i).gt.0) then
|
---|
[e40e335] | 35 | if (.not.fxvr(iphi(i))) then
|
---|
| 36 | nnonfx=nnonfx+1
|
---|
| 37 | endif
|
---|
| 38 | endif
|
---|
[32289cd] | 39 | if (.not.fxvr(ipsi(i))) then
|
---|
[e40e335] | 40 | nnonfx=nnonfx+1
|
---|
| 41 | endif
|
---|
| 42 | enddo
|
---|
| 43 | if (nnonfx.lt.6) then
|
---|
| 44 | ! print *,iaa,'bgs impossible because ndof = ',nnonfx
|
---|
| 45 | ians=.false.
|
---|
| 46 | endif
|
---|
| 47 | endif
|
---|
| 48 | if (ians) then
|
---|
| 49 | ! print *,'bgs is possible for angle ',ips,jv,nmvr(jv)
|
---|
| 50 | endif
|
---|
| 51 | bgsposs=ians
|
---|
| 52 | return
|
---|
| 53 | end
|
---|
| 54 |
|
---|
| 55 | ! Biased Gaussian Steps. Implements a semi-local conformational update
|
---|
[32289cd] | 56 | ! which modifies the protein backbone locally in a certain range of
|
---|
| 57 | ! amino acids. The 'down-stream' parts of the molecule outside the
|
---|
[e40e335] | 58 | ! region of update get small rigid body translations and rotations.
|
---|
| 59 | !
|
---|
[32289cd] | 60 | ! Use the update sparingly. It is rather local, and is not of great
|
---|
[e40e335] | 61 | ! value if we need big changes in the conformation. It is recommended
|
---|
[32289cd] | 62 | ! that this update be used to refine a structure around a low energy
|
---|
[e40e335] | 63 | ! candidate structure. Even at low energies, if you always
|
---|
[32289cd] | 64 | ! perform BGS, the chances of coming out of that minimum are small.
|
---|
| 65 | ! So, there is a probability bgsprob, which decides whether BGS or the
|
---|
| 66 | ! normal single angle update is used.
|
---|
[e40e335] | 67 | !
|
---|
| 68 | ! Calls: energy, dummy (function provided as argument), addang, (rand)
|
---|
| 69 | !
|
---|
| 70 |
|
---|
| 71 | integer function bgs(eol1,dummy)
|
---|
| 72 | include 'INCL.H'
|
---|
| 73 | include 'incl_lund.h'
|
---|
[32289cd] | 74 | double precision grnd, xiv, bv, ab, rv, dv, a, sum, p, r1, r2
|
---|
| 75 | double precision ppsi, wfw, addang, enw, energy, wbw, rd, delta
|
---|
| 76 | double precision dummy, eol1
|
---|
| 77 |
|
---|
| 78 | integer ivar, i, nph, jv, ia, nursvr, icurraa, j, k, l
|
---|
| 79 |
|
---|
[e40e335] | 80 | external dummy
|
---|
| 81 | dimension xiv(8,3),bv(8,3),rv(3,3),dv(3,8,3)
|
---|
| 82 | dimension ab(8), A(8,8),p(8),ppsi(8)
|
---|
| 83 | double precision ovr(mxvr)
|
---|
| 84 | ! Initialize
|
---|
| 85 | ! print *,'using BGS on angle ',nmvr(idvr(ivar))
|
---|
[32289cd] | 86 | if (bgsnvar.eq.0) then
|
---|
[e40e335] | 87 | bgs=0
|
---|
| 88 | goto 171
|
---|
| 89 | endif
|
---|
| 90 | ivar=1+grnd()*bgsnvar
|
---|
| 91 | do i=1,8
|
---|
| 92 | iph(i)=-50000
|
---|
| 93 | dph(i)=0
|
---|
| 94 | enddo
|
---|
| 95 | nph=0
|
---|
| 96 | jv=idvr(ivar)
|
---|
| 97 | ia=nursvr(jv)
|
---|
| 98 | ! Get BGS matrices based on coordinates of atoms in 4 amino acids
|
---|
| 99 | do i=1,4
|
---|
| 100 | icurraa=ia+i-1
|
---|
[32289cd] | 101 | if (iphi(icurraa).gt.0.and..not.fxvr(iphi(icurraa))) then
|
---|
[e40e335] | 102 | nph=nph+1
|
---|
| 103 | xiv(nph,1)=xat(iCa(icurraa))
|
---|
| 104 | xiv(nph,2)=yat(iCa(icurraa))
|
---|
| 105 | xiv(nph,3)=zat(iCa(icurraa))
|
---|
| 106 | bv(nph,1)=xiv(nph,1)-xat(iN(icurraa))
|
---|
| 107 | bv(nph,2)=xiv(nph,2)-yat(iN(icurraa))
|
---|
| 108 | bv(nph,3)=xiv(nph,3)-zat(iN(icurraa))
|
---|
| 109 | ab(nph)=bv(nph,1)*bv(nph,1)+bv(nph,2)*bv(nph,2)
|
---|
[bd2278d] | 110 | & +bv(nph,3)*bv(nph,3)
|
---|
[e40e335] | 111 | iph(nph)=iphi(icurraa)
|
---|
[32289cd] | 112 | endif
|
---|
[e40e335] | 113 | if (.not.fxvr(ipsi(icurraa))) then
|
---|
| 114 | nph=nph+1
|
---|
| 115 | xiv(nph,1)=xat(iC(icurraa))
|
---|
| 116 | xiv(nph,2)=yat(iC(icurraa))
|
---|
| 117 | xiv(nph,3)=zat(iC(icurraa))
|
---|
| 118 | bv(nph,1)=xiv(nph,1)-xat(iCa(icurraa))
|
---|
| 119 | bv(nph,2)=xiv(nph,2)-yat(iCa(icurraa))
|
---|
| 120 | bv(nph,3)=xiv(nph,3)-zat(iCa(icurraa))
|
---|
| 121 | ab(nph)=bv(nph,1)*bv(nph,1)+bv(nph,2)*bv(nph,2)
|
---|
[bd2278d] | 122 | & +bv(nph,3)*bv(nph,3)
|
---|
[e40e335] | 123 | iph(nph)=ipsi(icurraa)
|
---|
| 124 | endif
|
---|
| 125 | enddo
|
---|
| 126 | rv(1,1)=xat(iCa(ia+3))
|
---|
| 127 | rv(1,2)=yat(iCa(ia+3))
|
---|
| 128 | rv(1,3)=zat(iCa(ia+3))
|
---|
| 129 | rv(2,1)=xat(iC(ia+3))
|
---|
| 130 | rv(2,2)=yat(iC(ia+3))
|
---|
| 131 | rv(2,3)=zat(iC(ia+3))
|
---|
| 132 | rv(3,1)=xat(iC(ia+3)+1)
|
---|
| 133 | rv(3,2)=yat(iC(ia+3)+1)
|
---|
| 134 | rv(3,3)=zat(iC(ia+3)+1)
|
---|
| 135 |
|
---|
| 136 | do i=1,3
|
---|
| 137 | do j=1,nph
|
---|
| 138 | dv(i,j,1)=(1.0/ab(j))*(bv(j,2)*(rv(i,3)-xiv(j,3))-
|
---|
[bd2278d] | 139 | & bv(j,3)*(rv(i,2)-xiv(j,2)))
|
---|
[e40e335] | 140 | dv(i,j,2)=(-1.0/ab(j))*(bv(j,1)*(rv(i,3)-xiv(j,3))-
|
---|
[bd2278d] | 141 | & bv(j,3)*(rv(i,1)-xiv(j,1)))
|
---|
[e40e335] | 142 | dv(i,j,3)=(1.0/ab(j))*(bv(j,1)*(rv(i,2)-xiv(j,2))-
|
---|
[bd2278d] | 143 | & bv(j,2)*(rv(i,1)-xiv(j,1)))
|
---|
[e40e335] | 144 | enddo
|
---|
| 145 | enddo
|
---|
[32289cd] | 146 | do i=1,nph
|
---|
[e40e335] | 147 | do j=i,nph
|
---|
| 148 | A(i,j)=0
|
---|
| 149 | do k=1,3
|
---|
| 150 | do l=1,3
|
---|
| 151 | A(i,j)=A(i,j)+dv(k,i,l)*(dv(k,j,l))
|
---|
| 152 | enddo
|
---|
| 153 | enddo
|
---|
| 154 | A(i,j)=bbgs*A(i,j)
|
---|
| 155 | if (i.eq.j) then
|
---|
| 156 | A(i,j)=A(i,j)+1
|
---|
| 157 | endif
|
---|
| 158 | A(i,j)=0.5*abgs*A(i,j)
|
---|
| 159 | enddo
|
---|
| 160 | enddo
|
---|
| 161 | do i=1,nph
|
---|
| 162 | do j=i,nph
|
---|
| 163 | sum=A(i,j)
|
---|
| 164 | do k=i-1,1,-1
|
---|
| 165 | sum=sum-A(i,k)*A(j,k)
|
---|
| 166 | enddo
|
---|
[32289cd] | 167 | if (i.eq.j) then
|
---|
| 168 | p(i)=sqrt(sum)
|
---|
| 169 | else
|
---|
[e40e335] | 170 | A(j,i)=sum/p(i)
|
---|
| 171 | endif
|
---|
| 172 | enddo
|
---|
| 173 | enddo
|
---|
| 174 | ! Generate 8 Gaussian distributed small random angles
|
---|
| 175 | do i=1,8,2
|
---|
| 176 | r1=grnd()
|
---|
| 177 | ! In the rare event that this random number is 0, just take the next
|
---|
| 178 | if (r1.le.0) r1=grnd()
|
---|
| 179 | r1=sqrt(-log(r1))
|
---|
| 180 | r2=grnd()
|
---|
| 181 | ppsi(i)=r1*cos(pi2*r2)
|
---|
| 182 | ppsi(i+1)=r1*sin(pi2*r2)
|
---|
| 183 | enddo
|
---|
| 184 | do i=1,nph
|
---|
| 185 | dph(i)=0
|
---|
| 186 | enddo
|
---|
[32289cd] | 187 | ! Solve lower triangular matrix to get dphi proposals
|
---|
[e40e335] | 188 | do i=nph,1,-1
|
---|
| 189 | sum=ppsi(i)
|
---|
| 190 | do k=i+1,nph
|
---|
| 191 | sum=sum-A(k,i)*dph(k)
|
---|
| 192 | enddo
|
---|
| 193 | dph(i)=sum/p(i)
|
---|
| 194 | enddo
|
---|
| 195 | ! Calculate intrinsic (non-Boltzmann) weight for forward process
|
---|
| 196 | ! print *,'calculating intrinsic weight for forward process'
|
---|
| 197 | sum=0
|
---|
| 198 | do i=1,nph
|
---|
| 199 | sum=sum+ppsi(i)*ppsi(i)
|
---|
[32289cd] | 200 | enddo
|
---|
[e40e335] | 201 | wfw=exp(-sum)
|
---|
| 202 | do i=1,nph
|
---|
| 203 | wfw=wfw*p(i)
|
---|
| 204 | enddo
|
---|
| 205 |
|
---|
| 206 | ! Reconstruct chain and calculate new energy
|
---|
| 207 | ! print *,'going to assign changes to the chain'
|
---|
| 208 | ovr = vlvr
|
---|
| 209 | do i=1,nph
|
---|
| 210 | vlvr(iph(i))=addang(vlvr(iph(i)),dph(i))
|
---|
| 211 | enddo
|
---|
| 212 | enw = energy()
|
---|
| 213 | ! Calculate weight for reverse process for detail balance
|
---|
| 214 | ! print *,'proceeding to calculate weight for the reverse process'
|
---|
| 215 | nph=0
|
---|
| 216 | do i=1,4
|
---|
| 217 | icurraa=ia+i-1
|
---|
[32289cd] | 218 | if (iphi(icurraa).gt.0.and..not.fxvr(iphi(icurraa))) then
|
---|
[e40e335] | 219 | nph=nph+1
|
---|
| 220 | xiv(nph,1)=xat(iCa(icurraa))
|
---|
| 221 | xiv(nph,2)=yat(iCa(icurraa))
|
---|
| 222 | xiv(nph,3)=zat(iCa(icurraa))
|
---|
| 223 | bv(nph,1)=xiv(nph,1)-xat(iN(icurraa))
|
---|
| 224 | bv(nph,2)=xiv(nph,2)-yat(iN(icurraa))
|
---|
| 225 | bv(nph,3)=xiv(nph,3)-zat(iN(icurraa))
|
---|
| 226 | ab(nph)=bv(nph,1)*bv(nph,1)+bv(nph,2)*bv(nph,2)
|
---|
[bd2278d] | 227 | & +bv(nph,3)*bv(nph,3)
|
---|
[e40e335] | 228 | iph(nph)=iphi(icurraa)
|
---|
[32289cd] | 229 | endif
|
---|
| 230 | if (.not.fxvr(ipsi(icurraa))) then
|
---|
[e40e335] | 231 | nph=nph+1
|
---|
| 232 | xiv(nph,1)=xat(iC(icurraa))
|
---|
| 233 | xiv(nph,2)=yat(iC(icurraa))
|
---|
| 234 | xiv(nph,3)=zat(iC(icurraa))
|
---|
| 235 | bv(nph,1)=xiv(nph,1)-xat(iCa(icurraa))
|
---|
| 236 | bv(nph,2)=xiv(nph,2)-yat(iCa(icurraa))
|
---|
| 237 | bv(nph,3)=xiv(nph,3)-zat(iCa(icurraa))
|
---|
| 238 | ab(nph)=bv(nph,1)*bv(nph,1)+bv(nph,2)*bv(nph,2)
|
---|
[bd2278d] | 239 | & +bv(nph,3)*bv(nph,3)
|
---|
[e40e335] | 240 | iph(nph)=ipsi(icurraa)
|
---|
| 241 | endif
|
---|
| 242 | enddo
|
---|
| 243 | rv(1,1)=xat(iCa(ia+3))
|
---|
| 244 | rv(1,2)=yat(iCa(ia+3))
|
---|
| 245 | rv(1,3)=zat(iCa(ia+3))
|
---|
| 246 | rv(2,1)=xat(iC(ia+3))
|
---|
| 247 | rv(2,2)=yat(iC(ia+3))
|
---|
| 248 | rv(2,3)=zat(iC(ia+3))
|
---|
| 249 | rv(3,1)=xat(iC(ia+3)+1)
|
---|
| 250 | rv(3,2)=yat(iC(ia+3)+1)
|
---|
| 251 | rv(3,3)=zat(iC(ia+3)+1)
|
---|
| 252 |
|
---|
| 253 | do i=1,3
|
---|
| 254 | do j=1,nph
|
---|
| 255 | dv(i,j,1)=(1.0/ab(j))*(bv(j,2)*(rv(i,3)-xiv(j,3))-
|
---|
[bd2278d] | 256 | & bv(j,3)*(rv(i,2)-xiv(j,2)))
|
---|
[e40e335] | 257 | dv(i,j,2)=(-1.0/ab(j))*(bv(j,1)*(rv(i,3)-xiv(j,3))-
|
---|
[bd2278d] | 258 | & bv(j,3)*(rv(i,1)-xiv(j,1)))
|
---|
[e40e335] | 259 | dv(i,j,3)=(1.0/ab(j))*(bv(j,1)*(rv(i,2)-xiv(j,2))-
|
---|
[bd2278d] | 260 | & bv(j,2)*(rv(i,1)-xiv(j,1)))
|
---|
[e40e335] | 261 | enddo
|
---|
| 262 | enddo
|
---|
[32289cd] | 263 | do i=1,nph
|
---|
[e40e335] | 264 | do j=i,nph
|
---|
| 265 | A(i,j)=0
|
---|
| 266 | do k=1,3
|
---|
| 267 | do l=1,3
|
---|
| 268 | A(i,j)=A(i,j)+dv(k,i,l)*(dv(k,j,l))
|
---|
| 269 | enddo
|
---|
| 270 | enddo
|
---|
| 271 | A(i,j)=bbgs*A(i,j)
|
---|
[32289cd] | 272 | if (i.eq.j) then
|
---|
[e40e335] | 273 | A(i,j)=A(i,j)+1
|
---|
| 274 | endif
|
---|
| 275 | A(i,j)=0.5*abgs*A(i,j)
|
---|
| 276 | enddo
|
---|
| 277 | enddo
|
---|
| 278 | do i=1,nph
|
---|
| 279 | do j=i,nph
|
---|
| 280 | sum=A(i,j)
|
---|
| 281 | do k=i-1,1,-1
|
---|
| 282 | sum=sum-A(i,k)*A(j,k)
|
---|
| 283 | enddo
|
---|
[32289cd] | 284 | if (i.eq.j) then
|
---|
| 285 | p(i)=sqrt(sum)
|
---|
| 286 | else
|
---|
[e40e335] | 287 | A(j,i)=sum/p(i)
|
---|
| 288 | endif
|
---|
| 289 | enddo
|
---|
| 290 | enddo
|
---|
| 291 | do i=1,nph
|
---|
| 292 | ppsi(i)=p(i)*dph(i)
|
---|
| 293 | do j=i+1,nph
|
---|
| 294 | ppsi(i)=ppsi(i)+A(j,i)*dph(j)
|
---|
[32289cd] | 295 | enddo
|
---|
[e40e335] | 296 | enddo
|
---|
| 297 | sum=0
|
---|
| 298 | do i=1,nph
|
---|
| 299 | sum=sum+ppsi(i)*ppsi(i)
|
---|
| 300 | enddo
|
---|
[32289cd] | 301 | wbw=exp(-sum)
|
---|
[e40e335] | 302 | do i=1,nph
|
---|
| 303 | wbw=wbw*p(i)
|
---|
| 304 | enddo
|
---|
| 305 | ! Acceptance condition (includes additional weight)
|
---|
| 306 | rd=grnd()
|
---|
| 307 | ! print *,'generated selection random number ',rd
|
---|
| 308 | ! print *,'wfw/wbw = ',wfw/wbw
|
---|
| 309 | rd=-log(rd*wfw/wbw)
|
---|
| 310 | ! print *,'modified rd = ',rd
|
---|
| 311 | ! print *,'before calculating energy change'
|
---|
| 312 | delta = dummy(enw)-dummy(eol1)
|
---|
| 313 | ! print *,'delta = ',delta
|
---|
| 314 | ! call outpdb(0,'after.pdb')
|
---|
| 315 | ! print *,'after outpdb for after.pdb'
|
---|
[32289cd] | 316 | ! do i=1,nph
|
---|
[e40e335] | 317 | ! print *,'BGS>',i,iph(i),vlvr(iph(i)),dph(i)
|
---|
| 318 | ! enddo
|
---|
[32289cd] | 319 | if (rd.ge.delta) then
|
---|
[e40e335] | 320 | ! accept
|
---|
| 321 | eol1=enw
|
---|
| 322 | bgs=1
|
---|
| 323 | ! print *,'BGS move accepted'
|
---|
[32289cd] | 324 | else
|
---|
[e40e335] | 325 | ! reject
|
---|
| 326 | vlvr = ovr
|
---|
| 327 | ! enw=energy()
|
---|
| 328 | ! if (abs(enw-eol1).gt.0.000001) then
|
---|
[38d77eb] | 329 | ! write (logString, *) 'rejected bgs move: energy change :',eol1,enw
|
---|
[e40e335] | 330 | ! endif
|
---|
[38d77eb] | 331 | ! write (logString, *) 'rejected bgs move: ',eol1,enw,wfw,wbw,rd
|
---|
[e40e335] | 332 | bgs=0
|
---|
| 333 | endif
|
---|
| 334 | 171 continue
|
---|
| 335 | return
|
---|
| 336 | end
|
---|
| 337 |
|
---|