2 c MD with umbrella_sampling using Wolyne's distance measure as a constraint
3 implicit real*8 (a-h,o-z)
9 include 'COMMON.CONTROL'
13 include 'COMMON.LANGEVIN'
15 include 'COMMON.LANGEVIN.lang0'
17 include 'COMMON.CHAIN'
18 include 'COMMON.DERIV'
20 include 'COMMON.LOCAL'
21 include 'COMMON.INTERACT'
22 include 'COMMON.IOUNITS'
23 include 'COMMON.NAMES'
24 include 'COMMON.TIME1'
25 double precision uzap1,uzap2,hm1,hm2,hmnum
26 double precision ucdelan,dUcartan(3,0:MAXRES)
27 & ,dUxcartan(3,0:MAXRES),cdummy(3,0:MAXRES),
28 & duconst(3,0:MAXRES),duxconst(3,0:MAXRES)
29 integer kstart,kend,lstart,lend,idummy
30 double precision delta /1.0d-7/
40 write (iout,*) "Econstrq adaptive",adaptive," iset",iset
44 qfrag(i)=qwolynes(ifrag(1,i,iset),ifrag(2,i,iset),.true.
47 write (iout,*) "fragment",i," qfrag",qfrag(i)," iset",iset,
48 & " ifrag",ifrag(1,i,iset),ifrag(2,i,iset)," qinfrag",
51 Uconst=Uconst+wfrag(i,iset)*harmonic(qfrag(i),qinfrag(i,iset))
52 c Calculating the derivatives of Constraint energy with respect to Q
53 Ucdfrag=wfrag(i,iset)*harmonicprim(qfrag(i),
56 write (iout,*) "Uconst",Uconst," Ucdfrag",Ucdfrag
58 if (adaptive.and.i.eq.1) then
59 c write (iout,*) "Econstrq adative"
60 call PMF_energy(qfrag(1),ePMF,ePMF_q)
62 Ucdfrag=Ucdfrag+ePMF_q
64 write (iout,*) "PMF Uconst",Uconst," Ucdfrag",Ucdfrag
67 c hm1=harmonic(qfrag(i,iset),qinfrag(i,iset))
68 c hm2=harmonic(qfrag(i,iset)+delta,qinfrag(i,iset))
69 c hmnum=(hm2-hm1)/delta
70 c write(iout,*) "harmonicprim frag",harmonicprim(qfrag(i,iset),
72 c write(iout,*) "harmonicnum frag", hmnum
73 c Calculating the derivatives of Q with respect to cartesian coordinates
74 call qwolynes_prim(ifrag(1,i,iset),ifrag(2,i,iset),.true.
77 write(iout,*) "dqwol "
79 write(iout,'(i5,3e15.5)') ii,(dqwol(j,ii),j=1,3)
81 write(iout,*) "dxqwol "
83 write(iout,'(i5,3e15.5)') ii,(dxqwol(j,ii),j=1,3)
86 c Calculating numerical gradients of dU/dQi and dQi/dxi
87 c call qwol_num(ifrag(1,i,iset),ifrag(2,i,iset),.true.
89 c The gradients of Uconst in Cs
92 duconst(j,ii)=dUconst(j,ii)+ucdfrag*dqwol(j,ii)
93 dUxconst(j,ii)=dUxconst(j,ii)+ucdfrag*dxqwol(j,ii)
98 kstart=ifrag(1,ipair(1,i,iset),iset)
99 kend=ifrag(2,ipair(1,i,iset),iset)
100 lstart=ifrag(1,ipair(2,i,iset),iset)
101 lend=ifrag(2,ipair(2,i,iset),iset)
102 qpair(i)=qwolynes(kstart,kend,.false.,lstart,lend)
103 Uconst=Uconst+wpair(i,iset)*harmonic(qpair(i),qinpair(i,iset))
105 Ucdpair=wpair(i,iset)*harmonicprim(qpair(i),qinpair(i,iset))
106 c hm1=harmonic(qpair(i),qinpair(i,iset))
107 c hm2=harmonic(qpair(i)+delta,qinpair(i,iset))
108 c hmnum=(hm2-hm1)/delta
109 c write(iout,*) "harmonicprim pair ",harmonicprim(qpair(i),
111 c write(iout,*) "harmonicnum pair ", hmnum
113 call qwolynes_prim(kstart,kend,.false.
115 c write(iout,*) "dqwol "
117 c write(iout,'(i5,3e15.5)') ii,(dqwol(j,ii),j=1,3)
119 c write(iout,*) "dxqwol "
121 c write(iout,'(i5,3e15.5)') ii,(dxqwol(j,ii),j=1,3)
123 c Calculating numerical gradients
124 c call qwol_num(kstart,kend,.false.
126 c The gradients of Uconst in Cs
129 duconst(j,ii)=dUconst(j,ii)+ucdpair*dqwol(j,ii)
130 dUxconst(j,ii)=dUxconst(j,ii)+ucdpair*dxqwol(j,ii)
134 c write(iout,*) "Uconst inside subroutine ", Uconst
135 c Transforming the gradients from Cs to dCs for the backbone
139 dudconst(k,i)=dudconst(k,i)+duconst(k,j)+duxconst(k,j)
143 c Transforming the gradients from Cs to dCs for the side chains
146 dudxconst(j,i)=duxconst(j,i)
150 write(iout,*) "dU/ddc backbone "
152 write(iout,'(i5,3e15.5)') ii, (dudconst(j,ii),j=1,3)
154 write(iout,*) "dU/ddX side chain "
156 write(iout,'(i5,3e15.5)') ii,(duxconst(j,ii),j=1,3)
159 c Calculating numerical gradients of dUconst/ddc and dUconst/ddx