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 include 'COMMON.QRESTR'
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'
159 include 'COMMON.LANGEVIN.lang0'
161 include 'COMMON.CHAIN'
162 include 'COMMON.DERIV'
164 include 'COMMON.LOCAL'
165 include 'COMMON.INTERACT'
166 include 'COMMON.IOUNITS'
167 include 'COMMON.NAMES'
168 include 'COMMON.TIME1'
169 double precision uzap1,uzap2,hm1,hm2,hmnum
170 double precision ucdelan,dUcartan(3,0:MAXRES)
171 & ,dUxcartan(3,0:MAXRES),cdummy(3,0:MAXRES),
172 & duconst(3,0:MAXRES),duxconst(3,0:MAXRES)
173 integer kstart,kend,lstart,lend,idummy
174 double precision delta /1.0d-7/
185 qfrag(i)=qwolynes(ifrag(1,i,iset),ifrag(2,i,iset),.true.
187 Uconst=Uconst+wfrag(i,iset)*harmonic(qfrag(i),qinfrag(i,iset))
188 c Calculating the derivatives of Constraint energy with respect to Q
189 Ucdfrag=wfrag(i,iset)*harmonicprim(qfrag(i),qinfrag(i,iset))
190 c Calculating the derivatives of Q with respect to cartesian coordinates
193 duconst(j,ii)=dUconst(j,ii)+ucdfrag*dqwol(j,ii)
194 dUxconst(j,ii)=dUxconst(j,ii)+ucdfrag*dxqwol(j,ii)
197 c write (iout,*) "Calling qwol_num"
198 c call qwol_num(ifrag(1,i),ifrag(2,i),.true.,idummy,idummy)
201 kstart=ifrag(1,ipair(1,i,iset),iset)
202 kend=ifrag(2,ipair(1,i,iset),iset)
203 lstart=ifrag(1,ipair(2,i,iset),iset)
204 lend=ifrag(2,ipair(2,i,iset),iset)
205 qpair(i)=qwolynes(kstart,kend,.false.,lstart,lend)
206 Uconst=Uconst+wpair(i,iset)*harmonic(qpair(i),qinpair(i,iset))
208 Ucdpair=wpair(i,iset)*harmonicprim(qpair(i),qinpair(i,iset))
212 duconst(j,ii)=dUconst(j,ii)+ucdpair*dqwol(j,ii)
213 dUxconst(j,ii)=dUxconst(j,ii)+ucdpair*dxqwol(j,ii)
217 c write(iout,*) "Uconst inside subroutine ", Uconst
218 c Transforming the gradients from Cs to dCs for the backbone
222 dudconst(k,i)=dudconst(k,i)+duconst(k,j)+duxconst(k,j)
226 c Transforming the gradients from Cs to dCs for the side chains
229 dudxconst(j,i)=duxconst(j,i)
232 c write(iout,*) "dU/dc backbone "
234 c write(iout,'(i5,3e15.5)') ii, (duconst(j,ii),j=1,3)
236 c write(iout,*) "dU/dX side chain "
238 c write(iout,'(i5,3e15.5)') ii,(duxconst(j,ii),j=1,3)
240 c write(iout,*) "dU/ddc backbone "
242 c write(iout,'(i5,3e15.5)') ii, (dudconst(j,ii),j=1,3)
244 c write(iout,*) "dU/ddX side chain "
246 c write(iout,'(i5,3e15.5)') ii,(dudxconst(j,ii),j=1,3)
248 c Calculating numerical gradients of dUconst/ddc and dUconst/ddx
253 c-----------------------------------------------------------------------
254 subroutine dEconstrQ_num
255 c Calculating numerical dUconst/ddc and dUconst/ddx
256 implicit real*8 (a-h,o-z)
258 include 'COMMON.CONTROL'
262 include 'COMMON.LANGEVIN'
264 include 'COMMON.LANGEVIN.lang0'
266 include 'COMMON.CHAIN'
267 include 'COMMON.DERIV'
269 include 'COMMON.LOCAL'
270 include 'COMMON.INTERACT'
271 include 'COMMON.IOUNITS'
272 include 'COMMON.NAMES'
273 include 'COMMON.TIME1'
274 double precision uzap1,uzap2
275 double precision dUcartan(3,0:MAXRES)
276 & ,dUxcartan(3,0:MAXRES),cdummy(3,0:MAXRES)
277 integer kstart,kend,lstart,lend,idummy
278 double precision delta /1.0d-7/
284 dc(j,i)=dc(j,i)+delta
288 qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),
289 & .true.,idummy,idummy)
290 uzap2=uzap2+wfrag(ii,iset)*harmonic(qfrag(ii),
294 kstart=ifrag(1,ipair(1,ii,iset),iset)
295 kend=ifrag(2,ipair(1,ii,iset),iset)
296 lstart=ifrag(1,ipair(2,ii,iset),iset)
297 lend=ifrag(2,ipair(2,ii,iset),iset)
298 qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
299 uzap2=uzap2+wpair(ii,iset)*
300 & harmonic(qpair(ii),qinpair(ii,iset))
306 qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),
307 & .true.,idummy,idummy)
308 uzap1=uzap1+wfrag(ii,iset)*
309 & harmonic(qfrag(ii),qinfrag(ii,iset))
312 kstart=ifrag(1,ipair(1,ii,iset),iset)
313 kend=ifrag(2,ipair(1,ii,iset),iset)
314 lstart=ifrag(1,ipair(2,ii,iset),iset)
315 lend=ifrag(2,ipair(2,ii,iset),iset)
316 qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
317 uzap1=uzap1+wpair(ii,iset)*
318 & harmonic(qpair(ii),qinpair(ii,iset))
320 ducartan(j,i)=(uzap2-uzap1)/(delta)
323 c Calculating numerical gradients for dU/ddx
329 cdummy(j,i)=dc(j,i+nres)
330 dc(j,i+nres)=dc(j,i+nres)+delta
334 qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),
335 & .true.,idummy,idummy)
336 uzap2=uzap2+wfrag(ii,iset)*
337 & harmonic(qfrag(ii),qinfrag(ii,iset))
340 kstart=ifrag(1,ipair(1,ii,iset),iset)
341 kend=ifrag(2,ipair(1,ii,iset),iset)
342 lstart=ifrag(1,ipair(2,ii,iset),iset)
343 lend=ifrag(2,ipair(2,ii,iset),iset)
344 qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
345 uzap2=uzap2+wpair(ii,iset)*
346 & harmonic(qpair(ii),qinpair(ii,iset))
348 dc(j,i+nres)=cdummy(j,i)
352 qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),
353 & .true.,idummy,idummy)
354 uzap1=uzap1+wfrag(ii,iset)*
355 & harmonic(qfrag(ii),qinfrag(ii,iset))
358 kstart=ifrag(1,ipair(1,ii,iset),iset)
359 kend=ifrag(2,ipair(1,ii,iset),iset)
360 lstart=ifrag(1,ipair(2,ii,iset),iset)
361 lend=ifrag(2,ipair(2,ii,iset),iset)
362 qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
363 uzap1=uzap1+wpair(ii,iset)*
364 & harmonic(qpair(ii),qinpair(ii,iset))
366 duxcartan(j,i)=(uzap2-uzap1)/(delta)
369 write(iout,*) "Numerical dUconst/ddc backbone "
371 write(iout,'(i5,3e15.5)') ii,(dUcartan(j,ii),j=1,3)
373 write(iout,*) "Numerical dUconst/ddx side-chain "
375 write(iout,'(i5,3e15.5)') ii,(dUxcartan(j,ii),j=1,3)
380 c---------------------------------------------------------------------------
381 double precision function qcontrib(il,jl,il1,jl1)
384 include 'COMMON.IOUNITS'
385 include 'COMMON.CHAIN'
386 include 'COMMON.INTERACT'
387 include 'COMMON.QRESTR'
388 integer i,j,k,il,jl,il1,jl1,nd
389 double precision dist
391 double precision dij1,dij2,dij3,dij4,d0ij1,d0ij2,d0ij3,d0ij4,fac,
392 & fac1,ddave,ssij,ddqij
393 logical lprn /.false./
394 d0ij1=dsqrt((cref(1,jl)-cref(1,il))**2+
395 & (cref(2,jl)-cref(2,il))**2+
396 & (cref(3,jl)-cref(3,il))**2)
398 ddave=(dij1-d0ij1)**2
401 d0ij2=dsqrt((cref(1,jl1)-cref(1,il))**2+
402 & (cref(2,jl1)-cref(2,il))**2+
403 & (cref(3,jl1)-cref(3,il))**2)
405 ddave=ddave+(dij2-d0ij2)**2
409 d0ij3=dsqrt((cref(1,jl)-cref(1,il1))**2+
410 & (cref(2,jl)-cref(2,il1))**2+
411 & (cref(3,jl)-cref(3,il1))**2)
413 ddave=ddave+(dij3-d0ij3)**2
416 if (il1.ne.il .and. jl1.ne.jl) then
417 d0ij4=dsqrt((cref(1,jl1)-cref(1,il1))**2+
418 & (cref(2,jl1)-cref(2,il1))**2+
419 & (cref(3,jl1)-cref(3,il1))**2)
421 ddave=ddave+(dij4-d0ij4)**2
426 write (iout,*) "il",il," jl",jl,
427 & " itype",itype(il),itype(jl)," nd",nd
428 write (iout,*)"d0ij",d0ij1,d0ij2,d0ij3,d0ij4,
429 & " dij",dij1,dij2,dij3,dij4," ddave",ddave
432 c ssij = (0.25d0*d0ij1)**2
433 if (il.ne.il1 .and. jl.ne.jl1) then
434 ssij = 16.0d0/(d0ij1*d0ij4)
436 ssij = 16.0d0/(d0ij1*d0ij1)
438 qcontrib = dexp(-0.5d0*ddave*ssij)
440 fac1 = qcontrib*ssij/nd
441 fac = fac1*(dij1-d0ij1)/dij1
443 ddqij = (c(k,il)-c(k,jl))*fac
444 dqwol(k,il)=dqwol(k,il)+ddqij
445 dqwol(k,jl)=dqwol(k,jl)-ddqij
448 fac = fac1*(dij2-d0ij2)/dij2
450 ddqij = (c(k,il)-c(k,jl1))*fac
451 dqwol(k,il)=dqwol(k,il)+ddqij
452 dxqwol(k,jl)=dxqwol(k,jl)-ddqij
456 fac = fac1*(dij3-d0ij3)/dij3
458 ddqij = (c(k,il1)-c(k,jl))*fac
459 dxqwol(k,il)=dxqwol(k,il)+ddqij
460 dqwol(k,jl)=dqwol(k,jl)-ddqij
463 if (il1.ne.il .and. jl1.ne.jl) then
464 fac = fac1*(dij4-d0ij4)/dij4
466 ddqij = (c(k,il1)-c(k,jl1))*fac
467 dxqwol(k,il)=dxqwol(k,il)+ddqij
468 dxqwol(k,jl)=dxqwol(k,jl)-ddqij