1 double precision function qwolynes(seg1,seg2,flag,seg3,seg4)
2 implicit real*8 (a-h,o-z)
4 include 'COMMON.IOUNITS'
6 include 'COMMON.INTERACT'
8 integer i,j,jl,k,l,il,kl,nl,np,ip,kp,seg1,seg2,seg3,seg4,
11 double precision dist,qm
12 double precision qq,qqij,qqijCM,dij,d0ij,dijCM,d0ijCM
13 logical lprn /.false./
15 double precision sigm,x
24 d0ij=dsqrt((cref(1,jl,kkk)-cref(1,il,kkk))**2+
25 & (cref(2,jl,kkk)-cref(2,il,kkk))**2+
26 & (cref(3,jl,kkk)-cref(3,il,kkk))**2)
28 qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2)
29 if (itype(il).ne.10 .or. itype(jl).ne.10) then
32 & (cref(1,jl+nres,kkk)-cref(1,il+nres,kkk))**2+
33 & (cref(2,jl+nres,kkk)-cref(2,il+nres,kkk))**2+
34 & (cref(3,jl+nres,kkk)-cref(3,il+nres,kkk))**2)
35 dijCM=dist(il+nres,jl+nres)
36 qqijCM = dexp(-0.5d0*((dijCM-d0ijCM)/(sigm(d0ijCM)))**2)
44 if((seg3-il).lt.3) then
51 d0ij=dsqrt((cref(1,jl,kkk)-cref(1,il,kkk))**2+
52 & (cref(2,jl,kkk)-cref(2,il,kkk))**2+
53 & (cref(3,jl,kkk)-cref(3,il,kkk))**2)
55 qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2)
56 if (itype(il).ne.10 .or. itype(jl).ne.10) then
59 & (cref(1,jl+nres,kkk)-cref(1,il+nres,kkk))**2+
60 & (cref(2,jl+nres,kkk)-cref(2,il+nres,kkk))**2+
61 & (cref(3,jl+nres,kkk)-cref(3,il+nres,kkk))**2)
62 dijCM=dist(il+nres,jl+nres)
63 qqijCM = dexp(-0.5d0*((dijCM-d0ijCM)/(sigm(d0ijCM)))**2)
70 if (qqmax.le.qq) qqmax=qq
75 c-------------------------------------------------------------------
76 subroutine qwolynes_prim(seg1,seg2,flag,seg3,seg4)
77 implicit real*8 (a-h,o-z)
79 include 'COMMON.IOUNITS'
80 include 'COMMON.CHAIN'
81 include 'COMMON.INTERACT'
84 integer i,j,jl,k,l,il,nl,seg1,seg2,seg3,seg4,
88 double precision dij,d0ij,dijCM,d0ijCM
89 logical lprn /.false./
91 double precision sigm,x,sim,dd0,fac,ddqij
105 d0ij=dsqrt((cref(1,jl,kkk)-cref(1,il,kkk))**2+
106 & (cref(2,jl,kkk)-cref(2,il,kkk))**2+
107 & (cref(3,jl,kkk)-cref(3,il,kkk))**2)
109 sim = 1.0d0/sigm(d0ij)
112 fac = dd0*sim/dij*dexp(-0.5d0*dd0*dd0*sim)
114 ddqij = (c(k,il)-c(k,jl))*fac
115 dqwol(k,il)=dqwol(k,il)+ddqij
116 dqwol(k,jl)=dqwol(k,jl)-ddqij
119 if (itype(il).ne.10 .or. itype(jl).ne.10) then
122 & (cref(1,jl+nres,kkk)-cref(1,il+nres,kkk))**2+
123 & (cref(2,jl+nres,kkk)-cref(2,il+nres,kkk))**2+
124 & (cref(3,jl+nres,kkk)-cref(3,il+nres,kkk))**2)
125 dijCM=dist(il+nres,jl+nres)
126 sim = 1.0d0/sigm(d0ijCM)
129 fac=dd0*sim/dijCM*dexp(-0.5d0*dd0*dd0*sim)
131 ddqij = (c(k,il+nres)-c(k,jl+nres))*fac
132 dxqwol(k,il)=dxqwol(k,il)+ddqij
133 dxqwol(k,jl)=dxqwol(k,jl)-ddqij
140 if((seg3-il).lt.3) then
147 d0ij=dsqrt((cref(1,jl,kkk)-cref(1,il,kkk))**2+
148 & (cref(2,jl,kkk)-cref(2,il,kkk))**2+
149 & (cref(3,jl,kkk)-cref(3,il,kkk))**2)
151 sim = 1.0d0/sigm(d0ij)
154 fac = dd0*sim/dij*dexp(-0.5d0*dd0*dd0*sim)
156 ddqij = (c(k,il)-c(k,jl))*fac
157 dqwol(k,il)=dqwol(k,il)+ddqij
158 dqwol(k,jl)=dqwol(k,jl)-ddqij
160 if (itype(il).ne.10 .or. itype(jl).ne.10) then
163 & (cref(1,jl+nres,kkk)-cref(1,il+nres,kkk))**2+
164 & (cref(2,jl+nres,kkk)-cref(2,il+nres,kkk))**2+
165 & (cref(3,jl+nres,kkk)-cref(3,il+nres,kkk))**2)
166 dijCM=dist(il+nres,jl+nres)
167 sim = 1.0d0/sigm(d0ijCM)
170 fac = dd0*sim/dijCM*dexp(-0.5d0*dd0*dd0*sim)
172 ddqij = (c(k,il+nres)-c(k,jl+nres))*fac
173 dxqwol(k,il)=dxqwol(k,il)+ddqij
174 dxqwol(k,jl)=dxqwol(k,jl)-ddqij
183 dqwol(j,i)=dqwol(j,i)/nl
184 dxqwol(j,i)=dxqwol(j,i)/nl
189 c-------------------------------------------------------------------
190 subroutine qwol_num(seg1,seg2,flag,seg3,seg4)
191 implicit real*8 (a-h,o-z)
193 include 'COMMON.IOUNITS'
194 include 'COMMON.CHAIN'
195 include 'COMMON.INTERACT'
197 integer seg1,seg2,seg3,seg4
199 double precision qwolan(3,0:maxres),cdummy(3,0:maxres2),
200 & qwolxan(3,0:maxres),q1,q2
201 double precision delta /1.0d-10/
204 q1=qwolynes(seg1,seg2,flag,seg3,seg4)
207 q2=qwolynes(seg1,seg2,flag,seg3,seg4)
208 qwolan(j,i)=(q2-q1)/delta
214 q1=qwolynes(seg1,seg2,flag,seg3,seg4)
215 cdummy(j,i+nres)=c(j,i+nres)
216 c(j,i+nres)=c(j,i+nres)+delta
217 q2=qwolynes(seg1,seg2,flag,seg3,seg4)
218 qwolxan(j,i)=(q2-q1)/delta
219 c(j,i+nres)=cdummy(j,i+nres)
222 c write(iout,*) "Numerical Q carteisan gradients backbone: "
224 c write(iout,'(i5,3e15.5)') i, (qwolan(j,i),j=1,3)
226 c write(iout,*) "Numerical Q carteisan gradients side-chain: "
228 c write(iout,'(i5,3e15.5)') i, (qwolxan(j,i),j=1,3)
232 c------------------------------------------------------------------------
234 c MD with umbrella_sampling using Wolyne's distance measure as a constraint
235 implicit real*8 (a-h,o-z)
237 include 'COMMON.CONTROL'
241 include 'COMMON.LANGEVIN'
244 include 'COMMON.LANGEVIN.lang0.5diag'
246 include 'COMMON.LANGEVIN.lang0'
249 include 'COMMON.CHAIN'
250 include 'COMMON.DERIV'
252 include 'COMMON.LOCAL'
253 include 'COMMON.INTERACT'
254 include 'COMMON.IOUNITS'
255 include 'COMMON.NAMES'
256 include 'COMMON.TIME1'
257 double precision uzap1,uzap2,hm1,hm2,hmnum
258 double precision ucdelan,dUcartan(3,0:MAXRES)
259 & ,dUxcartan(3,0:MAXRES),cdummy(3,0:MAXRES),
260 & duconst(3,0:MAXRES),duxconst(3,0:MAXRES)
261 integer kstart,kend,lstart,lend,idummy
262 double precision delta /1.0d-7/
273 qfrag(i)=qwolynes(ifrag(1,i,iset),ifrag(2,i,iset),.true.
275 Uconst=Uconst+wfrag(i,iset)*harmonic(qfrag(i),qinfrag(i,iset))
276 c Calculating the derivatives of Constraint energy with respect to Q
277 Ucdfrag=wfrag(i,iset)*harmonicprim(qfrag(i),
279 c hm1=harmonic(qfrag(i,iset),qinfrag(i,iset))
280 c hm2=harmonic(qfrag(i,iset)+delta,qinfrag(i,iset))
281 c hmnum=(hm2-hm1)/delta
282 c write(iout,*) "harmonicprim frag",harmonicprim(qfrag(i,iset),
284 c write(iout,*) "harmonicnum frag", hmnum
285 c Calculating the derivatives of Q with respect to cartesian coordinates
286 call qwolynes_prim(ifrag(1,i,iset),ifrag(2,i,iset),.true.
288 c write(iout,*) "dqwol "
290 c write(iout,'(i5,3e15.5)') ii,(dqwol(j,ii),j=1,3)
292 c write(iout,*) "dxqwol "
294 c write(iout,'(i5,3e15.5)') ii,(dxqwol(j,ii),j=1,3)
296 c Calculating numerical gradients of dU/dQi and dQi/dxi
297 c call qwol_num(ifrag(1,i,iset),ifrag(2,i,iset),.true.
299 c The gradients of Uconst in Cs
302 duconst(j,ii)=dUconst(j,ii)+ucdfrag*dqwol(j,ii)
303 dUxconst(j,ii)=dUxconst(j,ii)+ucdfrag*dxqwol(j,ii)
308 kstart=ifrag(1,ipair(1,i,iset),iset)
309 kend=ifrag(2,ipair(1,i,iset),iset)
310 lstart=ifrag(1,ipair(2,i,iset),iset)
311 lend=ifrag(2,ipair(2,i,iset),iset)
312 qpair(i)=qwolynes(kstart,kend,.false.,lstart,lend)
313 Uconst=Uconst+wpair(i,iset)*harmonic(qpair(i),qinpair(i,iset))
315 Ucdpair=wpair(i,iset)*harmonicprim(qpair(i),qinpair(i,iset))
316 c hm1=harmonic(qpair(i),qinpair(i,iset))
317 c hm2=harmonic(qpair(i)+delta,qinpair(i,iset))
318 c hmnum=(hm2-hm1)/delta
319 c write(iout,*) "harmonicprim pair ",harmonicprim(qpair(i),
321 c write(iout,*) "harmonicnum pair ", hmnum
323 call qwolynes_prim(kstart,kend,.false.
325 c write(iout,*) "dqwol "
327 c write(iout,'(i5,3e15.5)') ii,(dqwol(j,ii),j=1,3)
329 c write(iout,*) "dxqwol "
331 c write(iout,'(i5,3e15.5)') ii,(dxqwol(j,ii),j=1,3)
333 c Calculating numerical gradients
334 c call qwol_num(kstart,kend,.false.
336 c The gradients of Uconst in Cs
339 duconst(j,ii)=dUconst(j,ii)+ucdpair*dqwol(j,ii)
340 dUxconst(j,ii)=dUxconst(j,ii)+ucdpair*dxqwol(j,ii)
344 c write(iout,*) "Uconst inside subroutine ", Uconst
345 c Transforming the gradients from Cs to dCs for the backbone
349 dudconst(k,i)=dudconst(k,i)+duconst(k,j)+duxconst(k,j)
353 c Transforming the gradients from Cs to dCs for the side chains
356 dudxconst(j,i)=duxconst(j,i)
359 c write(iout,*) "dU/ddc backbone "
361 c write(iout,'(i5,3e15.5)') ii, (dudconst(j,ii),j=1,3)
363 c write(iout,*) "dU/ddX side chain "
365 c write(iout,'(i5,3e15.5)') ii,(duxconst(j,ii),j=1,3)
367 c Calculating numerical gradients of dUconst/ddc and dUconst/ddx
371 c-----------------------------------------------------------------------
372 subroutine dEconstrQ_num
373 c Calculating numerical dUconst/ddc and dUconst/ddx
374 implicit real*8 (a-h,o-z)
376 include 'COMMON.CONTROL'
380 include 'COMMON.LANGEVIN'
383 include 'COMMON.LANGEVIN.lang0.5diag'
385 include 'COMMON.LANGEVIN.lang0'
388 include 'COMMON.CHAIN'
389 include 'COMMON.DERIV'
391 include 'COMMON.LOCAL'
392 include 'COMMON.INTERACT'
393 include 'COMMON.IOUNITS'
394 include 'COMMON.NAMES'
395 include 'COMMON.TIME1'
396 double precision uzap1,uzap2
397 double precision dUcartan(3,0:MAXRES)
398 & ,dUxcartan(3,0:MAXRES),cdummy(3,0:MAXRES)
399 integer kstart,kend,lstart,lend,idummy
400 double precision delta /1.0d-7/
406 dc(j,i)=dc(j,i)+delta
410 qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),.true.
412 uzap2=uzap2+wfrag(ii,iset)*harmonic(qfrag(ii),
416 kstart=ifrag(1,ipair(1,ii,iset),iset)
417 kend=ifrag(2,ipair(1,ii,iset),iset)
418 lstart=ifrag(1,ipair(2,ii,iset),iset)
419 lend=ifrag(2,ipair(2,ii,iset),iset)
420 qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
421 uzap2=uzap2+wpair(ii,iset)*harmonic(qpair(ii),
428 qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),.true.
430 uzap1=uzap1+wfrag(ii,iset)*harmonic(qfrag(ii),
434 kstart=ifrag(1,ipair(1,ii,iset),iset)
435 kend=ifrag(2,ipair(1,ii,iset),iset)
436 lstart=ifrag(1,ipair(2,ii,iset),iset)
437 lend=ifrag(2,ipair(2,ii,iset),iset)
438 qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
439 uzap1=uzap1+wpair(ii,iset)*harmonic(qpair(ii),
442 ducartan(j,i)=(uzap2-uzap1)/(delta)
445 c Calculating numerical gradients for dU/ddx
449 cdummy(j,i)=dc(j,i+nres)
450 dc(j,i+nres)=dc(j,i+nres)+delta
454 qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),.true.
456 uzap2=uzap2+wfrag(ii,iset)*harmonic(qfrag(ii),
460 kstart=ifrag(1,ipair(1,ii,iset),iset)
461 kend=ifrag(2,ipair(1,ii,iset),iset)
462 lstart=ifrag(1,ipair(2,ii,iset),iset)
463 lend=ifrag(2,ipair(2,ii,iset),iset)
464 qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
465 uzap2=uzap2+wpair(ii,iset)*harmonic(qpair(ii),
468 dc(j,i+nres)=cdummy(j,i)
472 qfrag(ii)=qwolynes(ifrag(1,ii,iset),
473 & ifrag(2,ii,iset),.true.,idummy,idummy)
474 uzap1=uzap1+wfrag(ii,iset)*harmonic(qfrag(ii),
478 kstart=ifrag(1,ipair(1,ii,iset),iset)
479 kend=ifrag(2,ipair(1,ii,iset),iset)
480 lstart=ifrag(1,ipair(2,ii,iset),iset)
481 lend=ifrag(2,ipair(2,ii,iset),iset)
482 qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
483 uzap1=uzap1+wpair(ii,iset)*harmonic(qpair(ii),
486 duxcartan(j,i)=(uzap2-uzap1)/(delta)
489 write(iout,*) "Numerical dUconst/ddc backbone "
491 write(iout,'(i5,3e15.5)') ii,(dUcartan(j,ii),j=1,3)
493 c write(iout,*) "Numerical dUconst/ddx side-chain "
495 c write(iout,'(i5,3e15.5)') ii,(dUxcartan(j,ii),j=1,3)
499 c---------------------------------------------------------------------------