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 integer i,j,jl,k,l,il,kl,nl,np,seg1,seg2,seg3,seg4,secseg
11 double precision dist,qm
12 double precision qq,qqij,qqijCM,dij,d0ij,dijCM,d0ijCM
13 logical lprn /.false./
24 write (iout,*) "seg1",seg1," seg2",seg2," seg3",seg3," seg4",seg4,
32 if (itype(il).ne.10) then
37 if (itype(jl).ne.10) then
42 qqijCM = qcontrib(il,jl,ilnres,jlnres)
45 write (iout,*) "qqijCM",qqijCM
51 write (iout,*) "nl",nl," qq",qq
56 if((seg3-il).lt.3) then
63 if (itype(il).ne.10) then
68 if (itype(jl).ne.10) then
73 qqijCM = qcontrib(il,jl,ilnres,jlnres)
76 write (iout,*) "qqijCM",qqijCM
86 dqwol(j,i)=dqwol(j,i)/nl
87 dxqwol(j,i)=dxqwol(j,i)/nl
92 c-------------------------------------------------------------------
93 subroutine qwol_num(seg1,seg2,flag,seg3,seg4)
94 implicit real*8 (a-h,o-z)
96 include 'COMMON.IOUNITS'
97 include 'COMMON.CHAIN'
98 include 'COMMON.INTERACT'
101 integer seg1,seg2,seg3,seg4
103 double precision qwolan(3,0:maxres),cdummy(3,0:maxres2),
104 & qwolxan(3,0:maxres),q1,q2
105 double precision delta /1.0d-7/
106 write (iout,*) "seg1",seg1," seg2",seg2," seg3",seg3," seg4",seg4
107 write(iout,*) "dQ/dc backbone "
109 write(iout,'(i5,3e15.5)') i, (dqwol(j,i),j=1,3)
111 write(iout,*) "dQ/dX side chain "
113 write(iout,'(i5,3e15.5)') i,(dxqwol(j,i),j=1,3)
119 q1=qwolynes(seg1,seg2,flag,seg3,seg4)
120 c(j,i)=cdummy(j,i)+delta
121 q2=qwolynes(seg1,seg2,flag,seg3,seg4)
122 qwolan(j,i)=0.5d0*(q2-q1)/delta
124 c write (iout,*) "i",i," j",j," q1",q1," a2",q2
129 cdummy(j,i+nres)=c(j,i+nres)
130 c(j,i+nres)=c(j,i+nres)-delta
131 q1=qwolynes(seg1,seg2,flag,seg3,seg4)
132 c(j,i+nres)=cdummy(j,i+nres)+delta
133 q2=qwolynes(seg1,seg2,flag,seg3,seg4)
134 qwolxan(j,i)=0.5d0*(q2-q1)/delta
135 c(j,i+nres)=cdummy(j,i+nres)
138 write(iout,*) "Numerical Q cartesian gradients backbone: "
140 write(iout,'(i5,3e15.5)') i, (qwolan(j,i),j=1,3)
142 write(iout,*) "Numerical Q cartesian gradients side-chain: "
144 write(iout,'(i5,3e15.5)') i, (qwolxan(j,i),j=1,3)
148 c------------------------------------------------------------------------
150 c MD with umbrella_sampling using Wolyne's distance measure as a constraint
151 implicit real*8 (a-h,o-z)
153 include 'COMMON.CONTROL'
157 include 'COMMON.LANGEVIN'
160 include 'COMMON.LANGEVIN.lang0.5diag'
162 include 'COMMON.LANGEVIN.lang0'
165 include 'COMMON.CHAIN'
166 include 'COMMON.DERIV'
168 include 'COMMON.LOCAL'
169 include 'COMMON.INTERACT'
170 include 'COMMON.IOUNITS'
171 include 'COMMON.NAMES'
172 include 'COMMON.TIME1'
173 double precision uzap1,uzap2,hm1,hm2,hmnum
174 double precision ucdelan,dUcartan(3,0:MAXRES)
175 & ,dUxcartan(3,0:MAXRES),cdummy(3,0:MAXRES),
176 & duconst(3,0:MAXRES),duxconst(3,0:MAXRES)
177 integer kstart,kend,lstart,lend,idummy
178 double precision delta /1.0d-7/
189 qfrag(i)=qwolynes(ifrag(1,i,iset),ifrag(2,i,iset),.true.
191 Uconst=Uconst+wfrag(i,iset)*harmonic(qfrag(i),qinfrag(i,iset))
192 c Calculating the derivatives of Constraint energy with respect to Q
193 Ucdfrag=wfrag(i,iset)*harmonicprim(qfrag(i),qinfrag(i,iset))
194 c Calculating the derivatives of Q with respect to cartesian coordinates
197 duconst(j,ii)=dUconst(j,ii)+ucdfrag*dqwol(j,ii)
198 dUxconst(j,ii)=dUxconst(j,ii)+ucdfrag*dxqwol(j,ii)
201 c write (iout,*) "Calling qwol_num"
202 c call qwol_num(ifrag(1,i),ifrag(2,i),.true.,idummy,idummy)
205 kstart=ifrag(1,ipair(1,i,iset),iset)
206 kend=ifrag(2,ipair(1,i,iset),iset)
207 lstart=ifrag(1,ipair(2,i,iset),iset)
208 lend=ifrag(2,ipair(2,i,iset),iset)
209 qpair(i)=qwolynes(kstart,kend,.false.,lstart,lend)
210 Uconst=Uconst+wpair(i,iset)*harmonic(qpair(i),qinpair(i,iset))
212 Ucdpair=wpair(i,iset)*harmonicprim(qpair(i),qinpair(i,iset))
216 duconst(j,ii)=dUconst(j,ii)+ucdpair*dqwol(j,ii)
217 dUxconst(j,ii)=dUxconst(j,ii)+ucdpair*dxqwol(j,ii)
221 c write(iout,*) "Uconst inside subroutine ", Uconst
222 c Transforming the gradients from Cs to dCs for the backbone
226 dudconst(k,i)=dudconst(k,i)+duconst(k,j)+duxconst(k,j)
230 c Transforming the gradients from Cs to dCs for the side chains
233 dudxconst(j,i)=duxconst(j,i)
236 c write(iout,*) "dU/dc backbone "
238 c write(iout,'(i5,3e15.5)') ii, (duconst(j,ii),j=1,3)
240 c write(iout,*) "dU/dX side chain "
242 c write(iout,'(i5,3e15.5)') ii,(duxconst(j,ii),j=1,3)
244 c write(iout,*) "dU/ddc backbone "
246 c write(iout,'(i5,3e15.5)') ii, (dudconst(j,ii),j=1,3)
248 c write(iout,*) "dU/ddX side chain "
250 c write(iout,'(i5,3e15.5)') ii,(dudxconst(j,ii),j=1,3)
252 c Calculating numerical gradients of dUconst/ddc and dUconst/ddx
256 c-----------------------------------------------------------------------
257 subroutine dEconstrQ_num
258 c Calculating numerical dUconst/ddc and dUconst/ddx
259 implicit real*8 (a-h,o-z)
261 include 'COMMON.CONTROL'
265 include 'COMMON.LANGEVIN'
268 include 'COMMON.LANGEVIN.lang0.5diag'
270 include 'COMMON.LANGEVIN.lang0'
273 include 'COMMON.CHAIN'
274 include 'COMMON.DERIV'
276 include 'COMMON.LOCAL'
277 include 'COMMON.INTERACT'
278 include 'COMMON.IOUNITS'
279 include 'COMMON.NAMES'
280 include 'COMMON.TIME1'
281 double precision uzap1,uzap2
282 double precision dUcartan(3,0:MAXRES)
283 & ,dUxcartan(3,0:MAXRES),cdummy(3,0:MAXRES)
284 integer kstart,kend,lstart,lend,idummy
285 double precision delta /1.0d-7/
291 dc(j,i)=dc(j,i)+delta
295 qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),
296 & .true.,idummy,idummy)
297 uzap2=uzap2+wfrag(ii,iset)*harmonic(qfrag(ii),
301 kstart=ifrag(1,ipair(1,ii,iset),iset)
302 kend=ifrag(2,ipair(1,ii,iset),iset)
303 lstart=ifrag(1,ipair(2,ii,iset),iset)
304 lend=ifrag(2,ipair(2,ii,iset),iset)
305 qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
306 uzap2=uzap2+wpair(ii,iset)*
307 & harmonic(qpair(ii),qinpair(ii,iset))
313 qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),
314 & .true.,idummy,idummy)
315 uzap1=uzap1+wfrag(ii,iset)*
316 & harmonic(qfrag(ii),qinfrag(ii,iset))
319 kstart=ifrag(1,ipair(1,ii,iset),iset)
320 kend=ifrag(2,ipair(1,ii,iset),iset)
321 lstart=ifrag(1,ipair(2,ii,iset),iset)
322 lend=ifrag(2,ipair(2,ii,iset),iset)
323 qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
324 uzap1=uzap1+wpair(ii,iset)*
325 & harmonic(qpair(ii),qinpair(ii,iset))
327 ducartan(j,i)=(uzap2-uzap1)/(delta)
330 c Calculating numerical gradients for dU/ddx
336 cdummy(j,i)=dc(j,i+nres)
337 dc(j,i+nres)=dc(j,i+nres)+delta
341 qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),
342 & .true.,idummy,idummy)
343 uzap2=uzap2+wfrag(ii,iset)*
344 & harmonic(qfrag(ii),qinfrag(ii,iset))
347 kstart=ifrag(1,ipair(1,ii,iset),iset)
348 kend=ifrag(2,ipair(1,ii,iset),iset)
349 lstart=ifrag(1,ipair(2,ii,iset),iset)
350 lend=ifrag(2,ipair(2,ii,iset),iset)
351 qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
352 uzap2=uzap2+wpair(ii,iset)*
353 & harmonic(qpair(ii),qinpair(ii,iset))
355 dc(j,i+nres)=cdummy(j,i)
359 qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),
360 & .true.,idummy,idummy)
361 uzap1=uzap1+wfrag(ii,iset)*
362 & harmonic(qfrag(ii),qinfrag(ii,iset))
365 kstart=ifrag(1,ipair(1,ii,iset),iset)
366 kend=ifrag(2,ipair(1,ii,iset),iset)
367 lstart=ifrag(1,ipair(2,ii,iset),iset)
368 lend=ifrag(2,ipair(2,ii,iset),iset)
369 qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
370 uzap1=uzap1+wpair(ii,iset)*
371 & harmonic(qpair(ii),qinpair(ii,iset))
373 duxcartan(j,i)=(uzap2-uzap1)/(delta)
376 write(iout,*) "Numerical dUconst/ddc backbone "
378 write(iout,'(i5,3e15.5)') ii,(dUcartan(j,ii),j=1,3)
380 write(iout,*) "Numerical dUconst/ddx side-chain "
382 write(iout,'(i5,3e15.5)') ii,(dUxcartan(j,ii),j=1,3)
386 c---------------------------------------------------------------------------
387 double precision function qcontrib(il,jl,il1,jl1)
390 include 'COMMON.IOUNITS'
391 include 'COMMON.CHAIN'
392 include 'COMMON.INTERACT'
394 integer i,j,k,il,jl,il1,jl1,nd
395 double precision dist
397 double precision dij1,dij2,dij3,dij4,d0ij1,d0ij2,d0ij3,d0ij4,fac,
398 & fac1,ddave,ssij,ddqij
399 logical lprn /.false./
400 d0ij1=dsqrt((cref(1,jl)-cref(1,il))**2+
401 & (cref(2,jl)-cref(2,il))**2+
402 & (cref(3,jl)-cref(3,il))**2)
404 ddave=(dij1-d0ij1)**2
407 d0ij2=dsqrt((cref(1,jl1)-cref(1,il))**2+
408 & (cref(2,jl1)-cref(2,il))**2+
409 & (cref(3,jl1)-cref(3,il))**2)
411 ddave=ddave+(dij2-d0ij2)**2
415 d0ij3=dsqrt((cref(1,jl)-cref(1,il1))**2+
416 & (cref(2,jl)-cref(2,il1))**2+
417 & (cref(3,jl)-cref(3,il1))**2)
419 ddave=ddave+(dij3-d0ij3)**2
422 if (il1.ne.il .and. jl1.ne.jl) then
423 d0ij4=dsqrt((cref(1,jl1)-cref(1,il1))**2+
424 & (cref(2,jl1)-cref(2,il1))**2+
425 & (cref(3,jl1)-cref(3,il1))**2)
427 ddave=ddave+(dij4-d0ij4)**2
432 write (iout,*) "il",il," jl",jl,
433 & " itype",itype(il),itype(jl)," nd",nd
434 write (iout,*)"d0ij",d0ij1,d0ij2,d0ij3,d0ij4,
435 & " dij",dij1,dij2,dij3,dij4," ddave",ddave
438 c ssij = (0.25d0*d0ij1)**2
439 if (il.ne.il1 .and. jl.ne.jl1) then
440 ssij = 16.0d0/(d0ij1*d0ij4)
442 ssij = 16.0d0/(d0ij1*d0ij1)
444 qcontrib = dexp(-0.5d0*ddave*ssij)
446 fac1 = qcontrib*ssij/nd
447 fac = fac1*(dij1-d0ij1)/dij1
449 ddqij = (c(k,il)-c(k,jl))*fac
450 dqwol(k,il)=dqwol(k,il)+ddqij
451 dqwol(k,jl)=dqwol(k,jl)-ddqij
454 fac = fac1*(dij2-d0ij2)/dij2
456 ddqij = (c(k,il)-c(k,jl1))*fac
457 dqwol(k,il)=dqwol(k,il)+ddqij
458 dxqwol(k,jl)=dxqwol(k,jl)-ddqij
462 fac = fac1*(dij3-d0ij3)/dij3
464 ddqij = (c(k,il1)-c(k,jl))*fac
465 dxqwol(k,il)=dxqwol(k,il)+ddqij
466 dqwol(k,jl)=dqwol(k,jl)-ddqij
469 if (il1.ne.il .and. jl1.ne.jl) then
470 fac = fac1*(dij4-d0ij4)/dij4
472 ddqij = (c(k,il1)-c(k,jl1))*fac
473 dxqwol(k,il)=dxqwol(k,il)+ddqij
474 dxqwol(k,jl)=dxqwol(k,jl)-ddqij