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'
83 include 'COMMON.QRESTR'
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'
239 include 'COMMON.QRESTR'
241 include 'COMMON.LANGEVIN'
243 include 'COMMON.LANGEVIN.lang0'
245 include 'COMMON.CHAIN'
246 include 'COMMON.DERIV'
248 include 'COMMON.LOCAL'
249 include 'COMMON.INTERACT'
250 include 'COMMON.IOUNITS'
251 include 'COMMON.NAMES'
252 include 'COMMON.TIME1'
253 double precision uzap1,uzap2,hm1,hm2,hmnum
254 double precision ucdelan,dUcartan(3,0:MAXRES)
255 & ,dUxcartan(3,0:MAXRES),cdummy(3,0:MAXRES),
256 & duconst(3,0:MAXRES),duxconst(3,0:MAXRES)
257 integer kstart,kend,lstart,lend,idummy
258 double precision delta /1.0d-7/
269 qfrag(i)=qwolynes(ifrag(1,i,iset),ifrag(2,i,iset),.true.
271 Uconst=Uconst+wfrag(i,iset)*harmonic(qfrag(i),qinfrag(i,iset))
272 c Calculating the derivatives of Constraint energy with respect to Q
273 Ucdfrag=wfrag(i,iset)*harmonicprim(qfrag(i),
275 c hm1=harmonic(qfrag(i,iset),qinfrag(i,iset))
276 c hm2=harmonic(qfrag(i,iset)+delta,qinfrag(i,iset))
277 c hmnum=(hm2-hm1)/delta
278 c write(iout,*) "harmonicprim frag",harmonicprim(qfrag(i,iset),
280 c write(iout,*) "harmonicnum frag", hmnum
281 c Calculating the derivatives of Q with respect to cartesian coordinates
282 call qwolynes_prim(ifrag(1,i,iset),ifrag(2,i,iset),.true.
284 c write(iout,*) "dqwol "
286 c write(iout,'(i5,3e15.5)') ii,(dqwol(j,ii),j=1,3)
288 c write(iout,*) "dxqwol "
290 c write(iout,'(i5,3e15.5)') ii,(dxqwol(j,ii),j=1,3)
292 c Calculating numerical gradients of dU/dQi and dQi/dxi
293 c call qwol_num(ifrag(1,i,iset),ifrag(2,i,iset),.true.
295 c The gradients of Uconst in Cs
298 duconst(j,ii)=dUconst(j,ii)+ucdfrag*dqwol(j,ii)
299 dUxconst(j,ii)=dUxconst(j,ii)+ucdfrag*dxqwol(j,ii)
304 kstart=ifrag(1,ipair(1,i,iset),iset)
305 kend=ifrag(2,ipair(1,i,iset),iset)
306 lstart=ifrag(1,ipair(2,i,iset),iset)
307 lend=ifrag(2,ipair(2,i,iset),iset)
308 qpair(i)=qwolynes(kstart,kend,.false.,lstart,lend)
309 Uconst=Uconst+wpair(i,iset)*harmonic(qpair(i),qinpair(i,iset))
311 Ucdpair=wpair(i,iset)*harmonicprim(qpair(i),qinpair(i,iset))
312 c hm1=harmonic(qpair(i),qinpair(i,iset))
313 c hm2=harmonic(qpair(i)+delta,qinpair(i,iset))
314 c hmnum=(hm2-hm1)/delta
315 c write(iout,*) "harmonicprim pair ",harmonicprim(qpair(i),
317 c write(iout,*) "harmonicnum pair ", hmnum
319 call qwolynes_prim(kstart,kend,.false.
321 c write(iout,*) "dqwol "
323 c write(iout,'(i5,3e15.5)') ii,(dqwol(j,ii),j=1,3)
325 c write(iout,*) "dxqwol "
327 c write(iout,'(i5,3e15.5)') ii,(dxqwol(j,ii),j=1,3)
329 c Calculating numerical gradients
330 c call qwol_num(kstart,kend,.false.
332 c The gradients of Uconst in Cs
335 duconst(j,ii)=dUconst(j,ii)+ucdpair*dqwol(j,ii)
336 dUxconst(j,ii)=dUxconst(j,ii)+ucdpair*dxqwol(j,ii)
340 c write(iout,*) "Uconst inside subroutine ", Uconst
341 c Transforming the gradients from Cs to dCs for the backbone
345 dudconst(k,i)=dudconst(k,i)+duconst(k,j)+duxconst(k,j)
349 c Transforming the gradients from Cs to dCs for the side chains
352 dudxconst(j,i)=duxconst(j,i)
355 c write(iout,*) "dU/ddc backbone "
357 c write(iout,'(i5,3e15.5)') ii, (dudconst(j,ii),j=1,3)
359 c write(iout,*) "dU/ddX side chain "
361 c write(iout,'(i5,3e15.5)') ii,(duxconst(j,ii),j=1,3)
363 c Calculating numerical gradients of dUconst/ddc and dUconst/ddx
368 c-----------------------------------------------------------------------
369 subroutine dEconstrQ_num
370 c Calculating numerical dUconst/ddc and dUconst/ddx
371 implicit real*8 (a-h,o-z)
373 include 'COMMON.CONTROL'
375 include 'COMMON.QRESTR'
377 include 'COMMON.LANGEVIN'
379 include 'COMMON.LANGEVIN.lang0'
381 include 'COMMON.CHAIN'
382 include 'COMMON.DERIV'
384 include 'COMMON.LOCAL'
385 include 'COMMON.INTERACT'
386 include 'COMMON.IOUNITS'
387 include 'COMMON.NAMES'
388 include 'COMMON.TIME1'
389 double precision uzap1,uzap2
390 double precision dUcartan(3,0:MAXRES)
391 & ,dUxcartan(3,0:MAXRES),cdummy(3,0:MAXRES)
392 integer kstart,kend,lstart,lend,idummy
393 double precision delta /1.0d-7/
399 dc(j,i)=dc(j,i)+delta
403 qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),.true.
405 uzap2=uzap2+wfrag(ii,iset)*harmonic(qfrag(ii),
409 kstart=ifrag(1,ipair(1,ii,iset),iset)
410 kend=ifrag(2,ipair(1,ii,iset),iset)
411 lstart=ifrag(1,ipair(2,ii,iset),iset)
412 lend=ifrag(2,ipair(2,ii,iset),iset)
413 qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
414 uzap2=uzap2+wpair(ii,iset)*harmonic(qpair(ii),
421 qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),.true.
423 uzap1=uzap1+wfrag(ii,iset)*harmonic(qfrag(ii),
427 kstart=ifrag(1,ipair(1,ii,iset),iset)
428 kend=ifrag(2,ipair(1,ii,iset),iset)
429 lstart=ifrag(1,ipair(2,ii,iset),iset)
430 lend=ifrag(2,ipair(2,ii,iset),iset)
431 qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
432 uzap1=uzap1+wpair(ii,iset)*harmonic(qpair(ii),
435 ducartan(j,i)=(uzap2-uzap1)/(delta)
438 c Calculating numerical gradients for dU/ddx
442 cdummy(j,i)=dc(j,i+nres)
443 dc(j,i+nres)=dc(j,i+nres)+delta
447 qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),.true.
449 uzap2=uzap2+wfrag(ii,iset)*harmonic(qfrag(ii),
453 kstart=ifrag(1,ipair(1,ii,iset),iset)
454 kend=ifrag(2,ipair(1,ii,iset),iset)
455 lstart=ifrag(1,ipair(2,ii,iset),iset)
456 lend=ifrag(2,ipair(2,ii,iset),iset)
457 qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
458 uzap2=uzap2+wpair(ii,iset)*harmonic(qpair(ii),
461 dc(j,i+nres)=cdummy(j,i)
465 qfrag(ii)=qwolynes(ifrag(1,ii,iset),
466 & ifrag(2,ii,iset),.true.,idummy,idummy)
467 uzap1=uzap1+wfrag(ii,iset)*harmonic(qfrag(ii),
471 kstart=ifrag(1,ipair(1,ii,iset),iset)
472 kend=ifrag(2,ipair(1,ii,iset),iset)
473 lstart=ifrag(1,ipair(2,ii,iset),iset)
474 lend=ifrag(2,ipair(2,ii,iset),iset)
475 qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
476 uzap1=uzap1+wpair(ii,iset)*harmonic(qpair(ii),
479 duxcartan(j,i)=(uzap2-uzap1)/(delta)
482 write(iout,*) "Numerical dUconst/ddc backbone "
484 write(iout,'(i5,3e15.5)') ii,(dUcartan(j,ii),j=1,3)
486 c write(iout,*) "Numerical dUconst/ddx side-chain "
488 c write(iout,'(i5,3e15.5)') ii,(dUxcartan(j,ii),j=1,3)
492 c---------------------------------------------------------------------------