2 c MD with umbrella_sampling using Wolyne's distance measure as a constraint
3 implicit real*8 (a-h,o-z)
5 include 'COMMON.CONTROL'
9 include 'COMMON.LANGEVIN'
11 include 'COMMON.LANGEVIN.lang0'
13 include 'COMMON.CHAIN'
14 include 'COMMON.DERIV'
16 include 'COMMON.LOCAL'
17 include 'COMMON.INTERACT'
18 include 'COMMON.IOUNITS'
19 include 'COMMON.NAMES'
20 include 'COMMON.TIME1'
21 double precision uzap1,uzap2,hm1,hm2,hmnum
22 double precision ucdelan,dUcartan(3,0:MAXRES)
23 & ,dUxcartan(3,0:MAXRES),cdummy(3,0:MAXRES),
24 & duconst(3,0:MAXRES),duxconst(3,0:MAXRES)
25 integer kstart,kend,lstart,lend,idummy
26 double precision delta /1.0d-7/
37 write (iout,*) "econstrq: nfrag",nfrag," iset",iset,
38 & " ifrag",(ifrag(1,i,iset),ifrag(2,i,iset),i=1,nfrag)
40 qfrag(i)=qwolynes(ifrag(1,i,iset),ifrag(2,i,iset),.true.
42 write (iout,*) "fragment",i," qfrag",qfrag(i)
44 Uconst=Uconst+wfrag(i,iset)*harmonic(qfrag(i),qinfrag(i,iset))
45 c Calculating the derivatives of Constraint energy with respect to Q
46 Ucdfrag=wfrag(i,iset)*harmonicprim(qfrag(i),
48 c hm1=harmonic(qfrag(i,iset),qinfrag(i,iset))
49 c hm2=harmonic(qfrag(i,iset)+delta,qinfrag(i,iset))
50 c hmnum=(hm2-hm1)/delta
51 c write(iout,*) "harmonicprim frag",harmonicprim(qfrag(i,iset),
53 c write(iout,*) "harmonicnum frag", hmnum
54 c Calculating the derivatives of Q with respect to cartesian coordinates
55 call qwolynes_prim(ifrag(1,i,iset),ifrag(2,i,iset),.true.
57 c write(iout,*) "dqwol "
59 c write(iout,'(i5,3e15.5)') ii,(dqwol(j,ii),j=1,3)
61 c write(iout,*) "dxqwol "
63 c write(iout,'(i5,3e15.5)') ii,(dxqwol(j,ii),j=1,3)
65 c Calculating numerical gradients of dU/dQi and dQi/dxi
66 c call qwol_num(ifrag(1,i,iset),ifrag(2,i,iset),.true.
68 c The gradients of Uconst in Cs
71 duconst(j,ii)=dUconst(j,ii)+ucdfrag*dqwol(j,ii)
72 dUxconst(j,ii)=dUxconst(j,ii)+ucdfrag*dxqwol(j,ii)
77 kstart=ifrag(1,ipair(1,i,iset),iset)
78 kend=ifrag(2,ipair(1,i,iset),iset)
79 lstart=ifrag(1,ipair(2,i,iset),iset)
80 lend=ifrag(2,ipair(2,i,iset),iset)
81 qpair(i)=qwolynes(kstart,kend,.false.,lstart,lend)
82 Uconst=Uconst+wpair(i,iset)*harmonic(qpair(i),qinpair(i,iset))
84 Ucdpair=wpair(i,iset)*harmonicprim(qpair(i),qinpair(i,iset))
85 c hm1=harmonic(qpair(i),qinpair(i,iset))
86 c hm2=harmonic(qpair(i)+delta,qinpair(i,iset))
87 c hmnum=(hm2-hm1)/delta
88 c write(iout,*) "harmonicprim pair ",harmonicprim(qpair(i),
90 c write(iout,*) "harmonicnum pair ", hmnum
92 call qwolynes_prim(kstart,kend,.false.
94 c write(iout,*) "dqwol "
96 c write(iout,'(i5,3e15.5)') ii,(dqwol(j,ii),j=1,3)
98 c write(iout,*) "dxqwol "
100 c write(iout,'(i5,3e15.5)') ii,(dxqwol(j,ii),j=1,3)
102 c Calculating numerical gradients
103 c call qwol_num(kstart,kend,.false.
105 c The gradients of Uconst in Cs
108 duconst(j,ii)=dUconst(j,ii)+ucdpair*dqwol(j,ii)
109 dUxconst(j,ii)=dUxconst(j,ii)+ucdpair*dxqwol(j,ii)
113 c write(iout,*) "Uconst inside subroutine ", Uconst
114 c Transforming the gradients from Cs to dCs for the backbone
118 dudconst(k,i)=dudconst(k,i)+duconst(k,j)+duxconst(k,j)
122 c Transforming the gradients from Cs to dCs for the side chains
125 dudxconst(j,i)=duxconst(j,i)
128 c write(iout,*) "dU/ddc backbone "
130 c write(iout,'(i5,3e15.5)') ii, (dudconst(j,ii),j=1,3)
132 c write(iout,*) "dU/ddX side chain "
134 c write(iout,'(i5,3e15.5)') ii,(duxconst(j,ii),j=1,3)
136 c Calculating numerical gradients of dUconst/ddc and dUconst/ddx