logical flag
double precision sigm,x
sigm(x)=0.25d0*x
+#ifdef DEBUG
+ write (iout,*) "qwolynes: nperm",nperm," flag",flag,
+ & " seg1",seg1," seg2",seg2," nsep",nsep
+#endif
do kkk=1,nperm
qq = 0.0d0
nl=0
& (cref(3,jl,kkk)-cref(3,il,kkk))**2)
dij=dist(il,jl)
qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2)
+ qq = qq+qqij
if (itype(il).ne.10 .or. itype(jl).ne.10) then
nl=nl+1
d0ijCM=dsqrt(
& (cref(3,jl+nres,kkk)-cref(3,il+nres,kkk))**2)
dijCM=dist(il+nres,jl+nres)
qqijCM = dexp(-0.5d0*((dijCM-d0ijCM)/(sigm(d0ijCM)))**2)
+ qq = qq+qqijCM
endif
- qq = qq+qqij+qqijCM
+c write (iout,*) "il",il,itype(il)," jl",jl,itype(jl),
+c & " qqiij",qqij," qqijCM",qqijCM
enddo
enddo
+#ifdef DEBUG
+ write (iout,*) "qwolynes: nl",nl
+#endif
qq = qq/nl
else
do il=seg1,seg2
& (cref(3,jl,kkk)-cref(3,il,kkk))**2)
dij=dist(il,jl)
qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2)
+ qq = qq+qqij
if (itype(il).ne.10 .or. itype(jl).ne.10) then
nl=nl+1
d0ijCM=dsqrt(
dijCM=dist(il+nres,jl+nres)
qqijCM = dexp(-0.5d0*((dijCM-d0ijCM)/(sigm(d0ijCM)))**2)
endif
- qq = qq+qqij+qqijCM
+c write (iout,*) "il",il,itype(il)," jl",jl,itype(jl),
+c & " qqiij",qqij," qqijCM",qqijCM
+ qq = qq+qqijCM
enddo
enddo
qq = qq/nl
endif
- if (qqmax.le.qq) qqmax=qq
enddo
- qwolynes=1.0d0-qqmax
+c write (iout,*) "qq",qq
+ qwolynes=1.0d0-qq
return
end
c-------------------------------------------------------------------
logical flag
double precision sigm,x,sim,dd0,fac,ddqij
sigm(x)=0.25d0*x
+#ifdef DEBUG
+ write (iout,*) "qwolynes: flag",flag," seg1 seg1",seg1,seg2,
+ & " nsep",nsep
+ write (iout,*) "nperm",nperm
+#endif
do kkk=1,nperm
do i=0,nres
do j=1,3
dqwol(k,il)=dqwol(k,il)+ddqij
dqwol(k,jl)=dqwol(k,jl)-ddqij
enddo
-
if (itype(il).ne.10 .or. itype(jl).ne.10) then
nl=nl+1
d0ijCM=dsqrt(
dxqwol(k,jl)=dxqwol(k,jl)-ddqij
enddo
endif
+#ifdef DEBUG
+ write (iout,*) "prim il",il,itype(il)," jl",jl,itype(jl),
+ & " dqwol",(dqwol(k,il),k=1,3)," dxqwol",(dxqwol(k,il),k=1,3)
+#endif
enddo
enddo
else
enddo
endif
enddo
+#ifdef DEBUG
+ write (iout,*) "qwolynes: nl",nl
+#endif
do i=0,nres
do j=1,3
dqwol(j,i)=dqwol(j,i)/nl
return
end
c------------------------------------------------------------------------
- subroutine EconstrQ
-c MD with umbrella_sampling using Wolyne's distance measure as a constraint
- implicit real*8 (a-h,o-z)
- include 'DIMENSIONS'
- include 'COMMON.CONTROL'
- include 'COMMON.VAR'
- include 'COMMON.MD'
-#ifndef LANG0
- include 'COMMON.LANGEVIN'
-#else
- include 'COMMON.LANGEVIN.lang0'
-#endif
- include 'COMMON.CHAIN'
- include 'COMMON.DERIV'
- include 'COMMON.GEO'
- include 'COMMON.LOCAL'
- include 'COMMON.INTERACT'
- include 'COMMON.IOUNITS'
- include 'COMMON.NAMES'
- include 'COMMON.TIME1'
- double precision uzap1,uzap2,hm1,hm2,hmnum
- double precision ucdelan,dUcartan(3,0:MAXRES)
- & ,dUxcartan(3,0:MAXRES),cdummy(3,0:MAXRES),
- & duconst(3,0:MAXRES),duxconst(3,0:MAXRES)
- integer kstart,kend,lstart,lend,idummy
- double precision delta /1.0d-7/
- do i=0,nres
- do j=1,3
- duconst(j,i)=0.0d0
- dudconst(j,i)=0.0d0
- duxconst(j,i)=0.0d0
- dudxconst(j,i)=0.0d0
- enddo
- enddo
- Uconst=0.0d0
- do i=1,nfrag
- qfrag(i)=qwolynes(ifrag(1,i,iset),ifrag(2,i,iset),.true.
- & ,idummy,idummy)
- Uconst=Uconst+wfrag(i,iset)*harmonic(qfrag(i),qinfrag(i,iset))
-c Calculating the derivatives of Constraint energy with respect to Q
- Ucdfrag=wfrag(i,iset)*harmonicprim(qfrag(i),
- & qinfrag(i,iset))
-c hm1=harmonic(qfrag(i,iset),qinfrag(i,iset))
-c hm2=harmonic(qfrag(i,iset)+delta,qinfrag(i,iset))
-c hmnum=(hm2-hm1)/delta
-c write(iout,*) "harmonicprim frag",harmonicprim(qfrag(i,iset),
-c & qinfrag(i,iset))
-c write(iout,*) "harmonicnum frag", hmnum
-c Calculating the derivatives of Q with respect to cartesian coordinates
- call qwolynes_prim(ifrag(1,i,iset),ifrag(2,i,iset),.true.
- & ,idummy,idummy)
-c write(iout,*) "dqwol "
-c do ii=1,nres
-c write(iout,'(i5,3e15.5)') ii,(dqwol(j,ii),j=1,3)
-c enddo
-c write(iout,*) "dxqwol "
-c do ii=1,nres
-c write(iout,'(i5,3e15.5)') ii,(dxqwol(j,ii),j=1,3)
-c enddo
-c Calculating numerical gradients of dU/dQi and dQi/dxi
-c call qwol_num(ifrag(1,i,iset),ifrag(2,i,iset),.true.
-c & ,idummy,idummy)
-c The gradients of Uconst in Cs
- do ii=0,nres
- do j=1,3
- duconst(j,ii)=dUconst(j,ii)+ucdfrag*dqwol(j,ii)
- dUxconst(j,ii)=dUxconst(j,ii)+ucdfrag*dxqwol(j,ii)
- enddo
- enddo
- enddo
- do i=1,npair
- kstart=ifrag(1,ipair(1,i,iset),iset)
- kend=ifrag(2,ipair(1,i,iset),iset)
- lstart=ifrag(1,ipair(2,i,iset),iset)
- lend=ifrag(2,ipair(2,i,iset),iset)
- qpair(i)=qwolynes(kstart,kend,.false.,lstart,lend)
- Uconst=Uconst+wpair(i,iset)*harmonic(qpair(i),qinpair(i,iset))
-c Calculating dU/dQ
- Ucdpair=wpair(i,iset)*harmonicprim(qpair(i),qinpair(i,iset))
-c hm1=harmonic(qpair(i),qinpair(i,iset))
-c hm2=harmonic(qpair(i)+delta,qinpair(i,iset))
-c hmnum=(hm2-hm1)/delta
-c write(iout,*) "harmonicprim pair ",harmonicprim(qpair(i),
-c & qinpair(i,iset))
-c write(iout,*) "harmonicnum pair ", hmnum
-c Calculating dQ/dXi
- call qwolynes_prim(kstart,kend,.false.
- & ,lstart,lend)
-c write(iout,*) "dqwol "
-c do ii=1,nres
-c write(iout,'(i5,3e15.5)') ii,(dqwol(j,ii),j=1,3)
-c enddo
-c write(iout,*) "dxqwol "
-c do ii=1,nres
-c write(iout,'(i5,3e15.5)') ii,(dxqwol(j,ii),j=1,3)
-c enddo
-c Calculating numerical gradients
-c call qwol_num(kstart,kend,.false.
-c & ,lstart,lend)
-c The gradients of Uconst in Cs
- do ii=0,nres
- do j=1,3
- duconst(j,ii)=dUconst(j,ii)+ucdpair*dqwol(j,ii)
- dUxconst(j,ii)=dUxconst(j,ii)+ucdpair*dxqwol(j,ii)
- enddo
- enddo
- enddo
-c write(iout,*) "Uconst inside subroutine ", Uconst
-c Transforming the gradients from Cs to dCs for the backbone
- do i=0,nres
- do j=i+1,nres
- do k=1,3
- dudconst(k,i)=dudconst(k,i)+duconst(k,j)+duxconst(k,j)
- enddo
- enddo
- enddo
-c Transforming the gradients from Cs to dCs for the side chains
- do i=1,nres
- do j=1,3
- dudxconst(j,i)=duxconst(j,i)
- enddo
- enddo
-c write(iout,*) "dU/ddc backbone "
-c do ii=0,nres
-c write(iout,'(i5,3e15.5)') ii, (dudconst(j,ii),j=1,3)
-c enddo
-c write(iout,*) "dU/ddX side chain "
-c do ii=1,nres
-c write(iout,'(i5,3e15.5)') ii,(duxconst(j,ii),j=1,3)
-c enddo
-c Calculating numerical gradients of dUconst/ddc and dUconst/ddx
-c call dEconstrQ_num
- return
- end
-c-----------------------------------------------------------------------
- subroutine dEconstrQ_num
-c Calculating numerical dUconst/ddc and dUconst/ddx
- implicit real*8 (a-h,o-z)
- include 'DIMENSIONS'
- include 'COMMON.CONTROL'
- include 'COMMON.VAR'
- include 'COMMON.MD'
-#ifndef LANG0
- include 'COMMON.LANGEVIN'
-#else
- include 'COMMON.LANGEVIN.lang0'
-#endif
- include 'COMMON.CHAIN'
- include 'COMMON.DERIV'
- include 'COMMON.GEO'
- include 'COMMON.LOCAL'
- include 'COMMON.INTERACT'
- include 'COMMON.IOUNITS'
- include 'COMMON.NAMES'
- include 'COMMON.TIME1'
- double precision uzap1,uzap2
- double precision dUcartan(3,0:MAXRES)
- & ,dUxcartan(3,0:MAXRES),cdummy(3,0:MAXRES)
- integer kstart,kend,lstart,lend,idummy
- double precision delta /1.0d-7/
-c For the backbone
- do i=0,nres-1
- do j=1,3
- dUcartan(j,i)=0.0d0
- cdummy(j,i)=dc(j,i)
- dc(j,i)=dc(j,i)+delta
- call chainbuild_cart
- uzap2=0.0d0
- do ii=1,nfrag
- qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),.true.
- & ,idummy,idummy)
- uzap2=uzap2+wfrag(ii,iset)*harmonic(qfrag(ii),
- & qinfrag(ii,iset))
- enddo
- do ii=1,npair
- kstart=ifrag(1,ipair(1,ii,iset),iset)
- kend=ifrag(2,ipair(1,ii,iset),iset)
- lstart=ifrag(1,ipair(2,ii,iset),iset)
- lend=ifrag(2,ipair(2,ii,iset),iset)
- qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
- uzap2=uzap2+wpair(ii,iset)*harmonic(qpair(ii),
- & qinpair(ii,iset))
- enddo
- dc(j,i)=cdummy(j,i)
- call chainbuild_cart
- uzap1=0.0d0
- do ii=1,nfrag
- qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),.true.
- & ,idummy,idummy)
- uzap1=uzap1+wfrag(ii,iset)*harmonic(qfrag(ii),
- & qinfrag(ii,iset))
- enddo
- do ii=1,npair
- kstart=ifrag(1,ipair(1,ii,iset),iset)
- kend=ifrag(2,ipair(1,ii,iset),iset)
- lstart=ifrag(1,ipair(2,ii,iset),iset)
- lend=ifrag(2,ipair(2,ii,iset),iset)
- qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
- uzap1=uzap1+wpair(ii,iset)*harmonic(qpair(ii),
- & qinpair(ii,iset))
- enddo
- ducartan(j,i)=(uzap2-uzap1)/(delta)
- enddo
- enddo
-c Calculating numerical gradients for dU/ddx
- do i=0,nres-1
- duxcartan(j,i)=0.0d0
- do j=1,3
- cdummy(j,i)=dc(j,i+nres)
- dc(j,i+nres)=dc(j,i+nres)+delta
- call chainbuild_cart
- uzap2=0.0d0
- do ii=1,nfrag
- qfrag(ii)=qwolynes(ifrag(1,ii,iset),ifrag(2,ii,iset),.true.
- & ,idummy,idummy)
- uzap2=uzap2+wfrag(ii,iset)*harmonic(qfrag(ii),
- & qinfrag(ii,iset))
- enddo
- do ii=1,npair
- kstart=ifrag(1,ipair(1,ii,iset),iset)
- kend=ifrag(2,ipair(1,ii,iset),iset)
- lstart=ifrag(1,ipair(2,ii,iset),iset)
- lend=ifrag(2,ipair(2,ii,iset),iset)
- qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
- uzap2=uzap2+wpair(ii,iset)*harmonic(qpair(ii),
- & qinpair(ii,iset))
- enddo
- dc(j,i+nres)=cdummy(j,i)
- call chainbuild_cart
- uzap1=0.0d0
- do ii=1,nfrag
- qfrag(ii)=qwolynes(ifrag(1,ii,iset),
- & ifrag(2,ii,iset),.true.,idummy,idummy)
- uzap1=uzap1+wfrag(ii,iset)*harmonic(qfrag(ii),
- & qinfrag(ii,iset))
- enddo
- do ii=1,npair
- kstart=ifrag(1,ipair(1,ii,iset),iset)
- kend=ifrag(2,ipair(1,ii,iset),iset)
- lstart=ifrag(1,ipair(2,ii,iset),iset)
- lend=ifrag(2,ipair(2,ii,iset),iset)
- qpair(ii)=qwolynes(kstart,kend,.false.,lstart,lend)
- uzap1=uzap1+wpair(ii,iset)*harmonic(qpair(ii),
- & qinpair(ii,iset))
- enddo
- duxcartan(j,i)=(uzap2-uzap1)/(delta)
- enddo
- enddo
- write(iout,*) "Numerical dUconst/ddc backbone "
- do ii=0,nres
- write(iout,'(i5,3e15.5)') ii,(dUcartan(j,ii),j=1,3)
- enddo
-c write(iout,*) "Numerical dUconst/ddx side-chain "
-c do ii=1,nres
-c write(iout,'(i5,3e15.5)') ii,(dUxcartan(j,ii),j=1,3)
-c enddo
- return
- end
-c---------------------------------------------------------------------------
projects / unres.git / blobdiff