2 c MD with umbrella_sampling using Wolyne's distance measure as a constraint
9 include 'COMMON.CONTROL'
12 include 'COMMON.QRESTR'
14 include 'COMMON.LANGEVIN'
17 include 'COMMON.LANGEVIN.lang0.5diag'
19 include 'COMMON.LANGEVIN.lang0'
22 include 'COMMON.CHAIN'
23 include 'COMMON.DERIV'
25 include 'COMMON.LOCAL'
26 include 'COMMON.INTERACT'
27 include 'COMMON.IOUNITS'
28 include 'COMMON.NAMES'
29 include 'COMMON.TIME1'
30 double precision uzap1,uzap2,hm1,hm2,hmnum
31 double precision ucdelan,dUcartan(3,0:MAXRES)
32 & ,dUxcartan(3,0:MAXRES),cdummy(3,0:MAXRES),
33 & duconst(3,0:MAXRES),duxconst(3,0:MAXRES)
34 integer kstart,kend,lstart,lend,idummy
35 double precision delta /1.0d-7/
37 double precision qwolynes,harmonic,harmonicprim
38 double precision ePMF,ePMF_q
48 write (iout,*) "Econstrq adaptive",adaptive," iset",iset
52 qfrag(i)=qwolynes(ifrag(1,i,iset),ifrag(2,i,iset),.true.
55 write (iout,*) "fragment",i," qfrag",qfrag(i)," iset",iset,
56 & " ifrag",ifrag(1,i,iset),ifrag(2,i,iset)," qinfrag",
59 Uconst=Uconst+wfrag(i,iset)*harmonic(qfrag(i),qinfrag(i,iset))
60 c Calculating the derivatives of Constraint energy with respect to Q
61 Ucdfrag=wfrag(i,iset)*harmonicprim(qfrag(i),
64 write (iout,*) "Uconst",Uconst," Ucdfrag",Ucdfrag
66 if (adaptive.and.i.eq.1) then
67 c write (iout,*) "Econstrq adative"
68 call PMF_energy(qfrag(1),ePMF,ePMF_q)
70 Ucdfrag=Ucdfrag+ePMF_q
72 write (iout,*) "PMF Uconst",Uconst," Ucdfrag",Ucdfrag
75 c hm1=harmonic(qfrag(i,iset),qinfrag(i,iset))
76 c hm2=harmonic(qfrag(i,iset)+delta,qinfrag(i,iset))
77 c hmnum=(hm2-hm1)/delta
78 c write(iout,*) "harmonicprim frag",harmonicprim(qfrag(i,iset),
80 c write(iout,*) "harmonicnum frag", hmnum
81 c Calculating the derivatives of Q with respect to cartesian coordinates
82 call qwolynes_prim(ifrag(1,i,iset),ifrag(2,i,iset),.true.
85 write(iout,*) "dqwol "
87 write(iout,'(i5,3e15.5)') ii,(dqwol(j,ii),j=1,3)
89 write(iout,*) "dxqwol "
91 write(iout,'(i5,3e15.5)') ii,(dxqwol(j,ii),j=1,3)
94 c The gradients of Uconst in Cs
97 duconst(j,ii)=dUconst(j,ii)+ucdfrag*dqwol(j,ii)
98 dUxconst(j,ii)=dUxconst(j,ii)+ucdfrag*dxqwol(j,ii)
103 kstart=ifrag(1,ipair(1,i,iset),iset)
104 kend=ifrag(2,ipair(1,i,iset),iset)
105 lstart=ifrag(1,ipair(2,i,iset),iset)
106 lend=ifrag(2,ipair(2,i,iset),iset)
107 qpair(i)=qwolynes(kstart,kend,.false.,lstart,lend)
108 Uconst=Uconst+wpair(i,iset)*harmonic(qpair(i),qinpair(i,iset))
110 Ucdpair=wpair(i,iset)*harmonicprim(qpair(i),qinpair(i,iset))
111 c hm1=harmonic(qpair(i),qinpair(i,iset))
112 c hm2=harmonic(qpair(i)+delta,qinpair(i,iset))
113 c hmnum=(hm2-hm1)/delta
114 c write(iout,*) "harmonicprim pair ",harmonicprim(qpair(i),
116 c write(iout,*) "harmonicnum pair ", hmnum
118 call qwolynes_prim(kstart,kend,.false.
120 c write(iout,*) "dqwol "
122 c write(iout,'(i5,3e15.5)') ii,(dqwol(j,ii),j=1,3)
124 c write(iout,*) "dxqwol "
126 c write(iout,'(i5,3e15.5)') ii,(dxqwol(j,ii),j=1,3)
128 c The gradients of Uconst in Cs
131 duconst(j,ii)=dUconst(j,ii)+ucdpair*dqwol(j,ii)
132 dUxconst(j,ii)=dUxconst(j,ii)+ucdpair*dxqwol(j,ii)
136 c write(iout,*) "Uconst inside subroutine ", Uconst
137 c Transforming the gradients from Cs to dCs for the backbone
141 dudconst(k,i)=dudconst(k,i)+duconst(k,j)+duxconst(k,j)
145 c Transforming the gradients from Cs to dCs for the side chains
148 dudxconst(j,i)=duxconst(j,i)
152 write(iout,*) "dU/ddc backbone "
154 write(iout,'(i5,3e15.5)') ii, (dudconst(j,ii),j=1,3)
156 write(iout,*) "dU/ddX side chain "
158 write(iout,'(i5,3e15.5)') ii,(duxconst(j,ii),j=1,3)