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'
8 include 'COMMON.QRESTR'
10 include 'COMMON.LANGEVIN'
13 include 'COMMON.LANGEVIN.lang0.5diag'
15 include 'COMMON.LANGEVIN.lang0'
18 include 'COMMON.CHAIN'
19 include 'COMMON.DERIV'
21 include 'COMMON.LOCAL'
22 include 'COMMON.INTERACT'
23 include 'COMMON.IOUNITS'
24 include 'COMMON.NAMES'
25 include 'COMMON.TIME1'
26 double precision uzap1,uzap2,hm1,hm2,hmnum
27 double precision ucdelan,dUcartan(3,0:MAXRES)
28 & ,dUxcartan(3,0:MAXRES),cdummy(3,0:MAXRES),
29 & duconst(3,0:MAXRES),duxconst(3,0:MAXRES)
30 integer kstart,kend,lstart,lend,idummy
31 double precision delta /1.0d-7/
42 write (iout,*) "econstrq: nfrag",nfrag," iset",iset,
43 & " ifrag",(ifrag(1,i,iset),ifrag(2,i,iset),i=1,nfrag)
45 qfrag(i)=qwolynes(ifrag(1,i,iset),ifrag(2,i,iset),.true.
47 write (iout,*) "fragment",i," qfrag",qfrag(i)
49 Uconst=Uconst+wfrag(i,iset)*harmonic(qfrag(i),qinfrag(i,iset))
50 c Calculating the derivatives of Constraint energy with respect to Q
51 Ucdfrag=wfrag(i,iset)*harmonicprim(qfrag(i),
53 c hm1=harmonic(qfrag(i,iset),qinfrag(i,iset))
54 c hm2=harmonic(qfrag(i,iset)+delta,qinfrag(i,iset))
55 c hmnum=(hm2-hm1)/delta
56 c write(iout,*) "harmonicprim frag",harmonicprim(qfrag(i,iset),
58 c write(iout,*) "harmonicnum frag", hmnum
59 c Calculating the derivatives of Q with respect to cartesian coordinates
60 call qwolynes_prim(ifrag(1,i,iset),ifrag(2,i,iset),.true.
62 c write(iout,*) "dqwol "
64 c write(iout,'(i5,3e15.5)') ii,(dqwol(j,ii),j=1,3)
66 c write(iout,*) "dxqwol "
68 c write(iout,'(i5,3e15.5)') ii,(dxqwol(j,ii),j=1,3)
70 c The gradients of Uconst in Cs
73 duconst(j,ii)=dUconst(j,ii)+ucdfrag*dqwol(j,ii)
74 dUxconst(j,ii)=dUxconst(j,ii)+ucdfrag*dxqwol(j,ii)
79 kstart=ifrag(1,ipair(1,i,iset),iset)
80 kend=ifrag(2,ipair(1,i,iset),iset)
81 lstart=ifrag(1,ipair(2,i,iset),iset)
82 lend=ifrag(2,ipair(2,i,iset),iset)
83 qpair(i)=qwolynes(kstart,kend,.false.,lstart,lend)
84 Uconst=Uconst+wpair(i,iset)*harmonic(qpair(i),qinpair(i,iset))
86 Ucdpair=wpair(i,iset)*harmonicprim(qpair(i),qinpair(i,iset))
87 c hm1=harmonic(qpair(i),qinpair(i,iset))
88 c hm2=harmonic(qpair(i)+delta,qinpair(i,iset))
89 c hmnum=(hm2-hm1)/delta
90 c write(iout,*) "harmonicprim pair ",harmonicprim(qpair(i),
92 c write(iout,*) "harmonicnum pair ", hmnum
94 call qwolynes_prim(kstart,kend,.false.
96 c write(iout,*) "dqwol "
98 c write(iout,'(i5,3e15.5)') ii,(dqwol(j,ii),j=1,3)
100 c write(iout,*) "dxqwol "
102 c write(iout,'(i5,3e15.5)') ii,(dxqwol(j,ii),j=1,3)
104 c The gradients of Uconst in Cs
107 duconst(j,ii)=dUconst(j,ii)+ucdpair*dqwol(j,ii)
108 dUxconst(j,ii)=dUxconst(j,ii)+ucdpair*dxqwol(j,ii)
112 c write(iout,*) "Uconst inside subroutine ", Uconst
113 c Transforming the gradients from Cs to dCs for the backbone
117 dudconst(k,i)=dudconst(k,i)+duconst(k,j)+duxconst(k,j)
121 c Transforming the gradients from Cs to dCs for the side chains
124 dudxconst(j,i)=duxconst(j,i)
127 c write(iout,*) "dU/ddc backbone "
129 c write(iout,'(i5,3e15.5)') ii, (dudconst(j,ii),j=1,3)
131 c write(iout,*) "dU/ddX side chain "
133 c write(iout,'(i5,3e15.5)') ii,(duxconst(j,ii),j=1,3)