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