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'
9 include 'COMMON.QRESTR'
10 integer i,j,jl,k,l,il,kl,nl,np,seg1,seg2,seg3,seg4,secseg
12 double precision dist,qm
13 double precision qq,qqij,qqijCM,dij,d0ij,dijCM,d0ijCM
14 logical lprn /.false./
25 write (iout,*) "seg1",seg1," seg2",seg2," seg3",seg3," seg4",seg4,
33 if (itype(il).ne.10) then
38 if (itype(jl).ne.10) then
43 qqijCM = qcontrib(il,jl,ilnres,jlnres)
46 write (iout,*) "qqijCM",qqijCM
52 write (iout,*) "nl",nl," qq",qq
57 if((seg3-il).lt.3) then
64 if (itype(il).ne.10) then
69 if (itype(jl).ne.10) then
74 qqijCM = qcontrib(il,jl,ilnres,jlnres)
77 write (iout,*) "qqijCM",qqijCM
87 dqwol(j,i)=dqwol(j,i)/nl
88 dxqwol(j,i)=dxqwol(j,i)/nl
93 c-------------------------------------------------------------------
94 subroutine qwol_num(seg1,seg2,flag,seg3,seg4)
95 implicit real*8 (a-h,o-z)
97 include 'COMMON.IOUNITS'
98 include 'COMMON.CHAIN'
99 include 'COMMON.INTERACT'
102 include 'COMMON.QRESTR'
103 integer seg1,seg2,seg3,seg4
105 double precision qwolan(3,0:maxres),cdummy(3,0:maxres2),
106 & qwolxan(3,0:maxres),q1,q2
107 double precision delta /1.0d-7/
108 write (iout,*) "seg1",seg1," seg2",seg2," seg3",seg3," seg4",seg4
109 write(iout,*) "dQ/dc backbone "
111 write(iout,'(i5,3e15.5)') i, (dqwol(j,i),j=1,3)
113 write(iout,*) "dQ/dX side chain "
115 write(iout,'(i5,3e15.5)') i,(dxqwol(j,i),j=1,3)
121 q1=qwolynes(seg1,seg2,flag,seg3,seg4)
122 c(j,i)=cdummy(j,i)+delta
123 q2=qwolynes(seg1,seg2,flag,seg3,seg4)
124 qwolan(j,i)=0.5d0*(q2-q1)/delta
126 c write (iout,*) "i",i," j",j," q1",q1," a2",q2
131 cdummy(j,i+nres)=c(j,i+nres)
132 c(j,i+nres)=c(j,i+nres)-delta
133 q1=qwolynes(seg1,seg2,flag,seg3,seg4)
134 c(j,i+nres)=cdummy(j,i+nres)+delta
135 q2=qwolynes(seg1,seg2,flag,seg3,seg4)
136 qwolxan(j,i)=0.5d0*(q2-q1)/delta
137 c(j,i+nres)=cdummy(j,i+nres)
140 write(iout,*) "Numerical Q cartesian gradients backbone: "
142 write(iout,'(i5,3e15.5)') i, (qwolan(j,i),j=1,3)
144 write(iout,*) "Numerical Q cartesian gradients side-chain: "
146 write(iout,'(i5,3e15.5)') i, (qwolxan(j,i),j=1,3)
150 c------------------------------------------------------------------------
152 c MD with umbrella_sampling using Wolyne's distance measure as a constraint
153 implicit real*8 (a-h,o-z)
155 include 'COMMON.CONTROL'
159 include 'COMMON.LANGEVIN'
161 include 'COMMON.LANGEVIN.lang0'
163 include 'COMMON.CHAIN'
164 include 'COMMON.DERIV'
166 include 'COMMON.LOCAL'
167 include 'COMMON.INTERACT'
168 include 'COMMON.IOUNITS'
169 include 'COMMON.NAMES'
170 include 'COMMON.TIME1'
171 double precision uzap1,uzap2,hm1,hm2,hmnum
172 double precision ucdelan,dUcartan(3,0:MAXRES)
173 & ,dUxcartan(3,0:MAXRES),cdummy(3,0:MAXRES),
174 & duconst(3,0:MAXRES),duxconst(3,0:MAXRES)
175 integer kstart,kend,lstart,lend,idummy
176 double precision delta /1.0d-7/
187 qfrag(i)=qwolynes(ifrag(1,i,iset),ifrag(2,i,iset),.true.
189 Uconst=Uconst+wfrag(i,iset)*harmonic(qfrag(i),qinfrag(i,iset))
190 c Calculating the derivatives of Constraint energy with respect to Q
191 Ucdfrag=wfrag(i,iset)*harmonicprim(qfrag(i),qinfrag(i,iset))
192 c Calculating the derivatives of Q with respect to cartesian coordinates
195 duconst(j,ii)=dUconst(j,ii)+ucdfrag*dqwol(j,ii)
196 dUxconst(j,ii)=dUxconst(j,ii)+ucdfrag*dxqwol(j,ii)
199 c write (iout,*) "Calling qwol_num"
200 c call qwol_num(ifrag(1,i,iset),ifrag(2,i,iset),.true.,idummy,idummy)
204 kstart=ifrag(1,ipair(1,i,iset),iset)
205 kend=ifrag(2,ipair(1,i,iset),iset)
206 lstart=ifrag(1,ipair(2,i,iset),iset)
207 lend=ifrag(2,ipair(2,i,iset),iset)
208 qpair(i)=qwolynes(kstart,kend,.false.,lstart,lend)
209 Uconst=Uconst+wpair(i,iset)*harmonic(qpair(i),qinpair(i,iset))
211 Ucdpair=wpair(i,iset)*harmonicprim(qpair(i),qinpair(i,iset))
215 duconst(j,ii)=dUconst(j,ii)+ucdpair*dqwol(j,ii)
216 dUxconst(j,ii)=dUxconst(j,ii)+ucdpair*dxqwol(j,ii)
220 c write(iout,*) "Uconst inside subroutine ", Uconst
221 c Transforming the gradients from Cs to dCs for the backbone
225 dudconst(k,i)=dudconst(k,i)+duconst(k,j)+duxconst(k,j)
229 c Transforming the gradients from Cs to dCs for the side chains
232 dudxconst(j,i)=duxconst(j,i)
235 c write(iout,*) "dU/dc backbone "
237 c write(iout,'(i5,3e15.5)') ii, (duconst(j,ii),j=1,3)
239 c write(iout,*) "dU/dX side chain "
241 c write(iout,'(i5,3e15.5)') ii,(duxconst(j,ii),j=1,3)
243 c write(iout,*) "dU/ddc backbone "
245 c write(iout,'(i5,3e15.5)') ii, (dudconst(j,ii),j=1,3)
247 c write(iout,*) "dU/ddX side chain "
249 c write(iout,'(i5,3e15.5)') ii,(dudxconst(j,ii),j=1,3)
251 c Calculating numerical gradients of dUconst/ddc and dUconst/ddx
255 c-----------------------------------------------------------------------
256 subroutine dEconstrQ_num
257 c Calculating numerical dUconst/ddc and dUconst/ddx
258 implicit real*8 (a-h,o-z)
260 include 'COMMON.CONTROL'
264 include 'COMMON.LANGEVIN'
266 include 'COMMON.LANGEVIN.lang0'
268 include 'COMMON.CHAIN'
269 include 'COMMON.DERIV'
271 include 'COMMON.LOCAL'
272 include 'COMMON.INTERACT'
273 include 'COMMON.IOUNITS'
274 include 'COMMON.NAMES'
275 include 'COMMON.TIME1'
276 double precision uzap1,uzap2
277 double precision dUcartan(3,0:MAXRES)
278 & ,dUxcartan(3,0:MAXRES),cdummy(3,0:MAXRES)
279 integer kstart,kend,lstart,lend,idummy
280 double precision delta /1.0d-7/
286 dc(j,i)=dc(j,i)+delta
290 qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),
291 & .true.,idummy,idummy)
292 uzap2=uzap2+wfrag(ii,iset)*
293 & harmonic(qfrag(ii),qinfrag(ii,iset))
296 kstart=ifrag(1,ipair(1,ii,iset),iset)
297 kend=ifrag(2,ipair(1,ii,iset),iset)
298 lstart=ifrag(1,ipair(2,ii,iset),iset)
299 lend=ifrag(2,ipair(2,ii,iset),iset)
300 qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
301 uzap2=uzap2+wpair(ii,iset)*
302 & harmonic(qpair(ii),qinpair(ii,iset))
308 qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),
309 & .true.,idummy,idummy)
310 uzap1=uzap1+wfrag(ii,iset)*
311 & harmonic(qfrag(ii),qinfrag(ii,iset))
314 kstart=ifrag(1,ipair(1,ii,iset),iset)
315 kend=ifrag(2,ipair(1,ii,iset),iset)
316 lstart=ifrag(1,ipair(2,ii,iset),iset)
317 lend=ifrag(2,ipair(2,ii,iset),iset)
318 qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
319 uzap1=uzap1+wpair(ii,iset)*
320 & harmonic(qpair(ii),qinpair(ii,iset))
322 ducartan(j,i)=(uzap2-uzap1)/(delta)
325 c Calculating numerical gradients for dU/ddx
331 cdummy(j,i)=dc(j,i+nres)
332 dc(j,i+nres)=dc(j,i+nres)+delta
336 qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),
337 & .true.,idummy,idummy)
338 uzap2=uzap2+wfrag(ii,iset)*
339 & harmonic(qfrag(ii),qinfrag(ii,iset))
342 kstart=ifrag(1,ipair(1,ii,iset),iset)
343 kend=ifrag(2,ipair(1,ii,iset),iset)
344 lstart=ifrag(1,ipair(2,ii,iset),iset)
345 lend=ifrag(2,ipair(2,ii,iset),iset)
346 qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
347 uzap2=uzap2+wpair(ii,iset)*
348 & harmonic(qpair(ii),qinpair(ii,iset))
350 dc(j,i+nres)=cdummy(j,i)
354 qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),
355 & .true.,idummy,idummy)
356 uzap1=uzap1+wfrag(ii,iset)*
357 & harmonic(qfrag(ii),qinfrag(ii,iset))
360 kstart=ifrag(1,ipair(1,ii,iset),iset)
361 kend=ifrag(2,ipair(1,ii,iset),iset)
362 lstart=ifrag(1,ipair(2,ii,iset),iset)
363 lend=ifrag(2,ipair(2,ii,iset),iset)
364 qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
365 uzap1=uzap1+wpair(ii,iset)*
366 & harmonic(qpair(ii),qinpair(ii,iset))
368 duxcartan(j,i)=(uzap2-uzap1)/(delta)
371 write(iout,*) "Numerical dUconst/ddc backbone "
373 write(iout,'(i5,3e15.5)') ii,(dUcartan(j,ii),j=1,3)
375 write(iout,*) "Numerical dUconst/ddx side-chain "
377 write(iout,'(i5,3e15.5)') ii,(dUxcartan(j,ii),j=1,3)
381 c---------------------------------------------------------------------------
382 double precision function qcontrib(il,jl,il1,jl1)
385 include 'COMMON.IOUNITS'
386 include 'COMMON.CHAIN'
387 include 'COMMON.INTERACT'
389 include 'COMMON.LOCAL'
390 integer i,j,k,il,jl,il1,jl1,nd,itl,jtl
391 double precision dist
393 double precision dij,dij1,d0ij,d0ij1,om1,om2,om12,om10,om20,om120
394 & ,fac,fac1,ddave,ssij,ddqij,d0ii1,d0jj1,rij,eom1,eom2,eom12
395 double precision u(3),v(3),er(3),er0(3),dcosom1(3),dcosom2(3),
397 double precision scalar
399 logical lprn /.false./
400 if (lprn) write (iout,*) "il",il," jl",jl," il1",il1," jl1",jl1
401 d0ij=dsqrt((cref(1,jl)-cref(1,il))**2+
402 & (cref(2,jl)-cref(2,il))**2+
403 & (cref(3,jl)-cref(3,il))**2)
407 er(i)=(c(i,jl1)-c(i,il1))/dij1
410 er0(i)=cref(i,jl1)-cref(i,il1)
412 d0ij1=dsqrt(scalar(er0,er0))
416 if (il.ne.il1 .or. jl.ne.jl1) then
417 ddave=0.5d0*((dij-d0ij)**2+(dij1-d0ij1)**2)
425 u(i)=cref(i,il1)-cref(i,il)
427 d0ii1=dsqrt(scalar(u,u))
432 write (iout,*) "u",(u(i),i=1,3)
433 write (iout,*) "er0",(er0(i),i=1,3)
435 om1=scalar(er,dc_norm(1,il1))
436 write (iout,*) "om10",om10," om1",om1
444 v(i)=cref(i,jl1)-cref(i,jl)
446 d0jj1=dsqrt(scalar(v,v))
451 write (iout,*) "v",(v(i),i=1,3)
452 write (iout,*) "er0",(er0(i),i=1,3)
454 om2=scalar(er,dc_norm(1,jl1))
455 write (iout,*) "om20",om20," om2",om2
461 if (il.ne.il1 .and. jl.ne.jl1) then
463 om12=scalar(dc_norm(1,il1),dc_norm(1,jl1))
469 write (iout,*) "il",il," jl",jl,itype(il),itype(jl)
470 write (iout,*)"d0ij",d0ij," om10",om10," om20",om20,
472 & " dij",dij," om1",om1," om2",om2," om12",om12
475 ssij = 16.0d0/(d0ij*d0ij)
476 qcontrib = dexp(-0.5d0*(ddave*ssij+((om1-om10)**2
477 & +(om2-om20)**2+(om12-om120)**2)))
478 if (lprn) write (iout,*) "ssij",ssij," qcontrib",qcontrib
479 c qcontrib = dexp(-0.5d0*(ddave*ssij)+(om1-om10)**2+(om2-om20)**2)
480 c qcontrib = dexp(-0.5d0*(ddave*ssij))
481 c Compute gradient - radial component
482 fac1 = qcontrib*ssij/nd
483 fac = fac1*(dij-d0ij)/dij
485 ddqij = (c(k,il)-c(k,jl))*fac
486 dqwol(k,il)=dqwol(k,il)+ddqij
487 dqwol(k,jl)=dqwol(k,jl)-ddqij
489 if (il1.ne.il .or. jl1.ne.jl) then
490 fac = fac1*(dij1-d0ij1)/dij1
492 ddqij = (c(k,il1)-c(k,jl1))*fac
494 dxqwol(k,il)=dxqwol(k,il)+ddqij
496 dqwol(k,il)=dqwol(k,il)+ddqij
499 dxqwol(k,jl)=dxqwol(k,jl)-ddqij
501 dqwol(k,jl)=dqwol(k,jl)-ddqij
506 c Orientational contributions
508 eom1=qcontrib*(om1-om10)
509 eom2=qcontrib*(om2-om20)
510 eom12=qcontrib*(om12-om120)
512 dcosom1(k)=rij*(dc_norm(k,il1)-om1*er(k))
513 dcosom2(k)=rij*(dc_norm(k,jl1)-om2*er(k))
516 ddqij=eom1*dcosom1(k)+eom2*dcosom2(k)
517 aux1=(eom12*(dc_norm(k,jl1)-om12*dc_norm(k,il1))
518 & +eom1*(er(k)-om1*dc_norm(k,il1)))*vbld_inv(il1)
519 aux2=(eom12*(dc_norm(k,il1)-om12*dc_norm(k,jl1))
520 & +eom2*(er(k)-om2*dc_norm(k,jl1)))*vbld_inv(jl1)
521 dqwol(k,il)=dqwol(k,il)-ddqij-aux1
522 dqwol(k,jl)=dqwol(k,jl)+ddqij-aux2
523 dxqwol(k,il)=dxqwol(k,il)-ddqij+aux1
524 c & +(eom12*(dc_norm(k,jl1)-om12*dc_norm(k,il1))
525 c & +eom1*(er(k)-om1*dc_norm(k,il1)))*vbld_inv(il1)
526 dxqwol(k,jl)=dxqwol(k,jl)+ddqij+aux2
527 c & +(eom12*(dc_norm(k,il1)-om12*dc_norm(k,jl1))
528 c & +eom2*(er(k)-om2*dc_norm(k,jl1)))*vbld_inv(jl1)