--- /dev/null
+ double precision function qwolynes(seg1,seg2,flag,seg3,seg4)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.CHAIN'
+ include 'COMMON.INTERACT'
+ include 'COMMON.VAR'
+ integer i,j,jl,k,l,il,kl,nl,np,ip,kp,seg1,seg2,seg3,seg4,
+ & secseg
+ integer nsep /3/
+ double precision dist,qm
+ double precision qq,qqij,qqijCM,dij,d0ij,dijCM,d0ijCM
+ logical lprn /.false./
+ logical flag
+ double precision sigm,x
+ sigm(x)=0.25d0*x
+ qq = 0.0d0
+ nl=0
+ if(flag) then
+ do il=seg1+nsep,seg2
+ do jl=seg1,il-nsep
+ nl=nl+1
+ d0ij=dsqrt((cref(1,jl)-cref(1,il))**2+
+ & (cref(2,jl)-cref(2,il))**2+
+ & (cref(3,jl)-cref(3,il))**2)
+ dij=dist(il,jl)
+ qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2)
+ if (itype(il).ne.10 .or. itype(jl).ne.10) then
+ nl=nl+1
+ d0ijCM=dsqrt(
+ & (cref(1,jl+nres)-cref(1,il+nres))**2+
+ & (cref(2,jl+nres)-cref(2,il+nres))**2+
+ & (cref(3,jl+nres)-cref(3,il+nres))**2)
+ dijCM=dist(il+nres,jl+nres)
+ qqijCM = dexp(-0.5d0*((dijCM-d0ijCM)/(sigm(d0ijCM)))**2)
+ endif
+ qq = qq+qqij+qqijCM
+ enddo
+ enddo
+ qq = qq/nl
+ else
+ do il=seg1,seg2
+ if((seg3-il).lt.3) then
+ secseg=il+3
+ else
+ secseg=seg3
+ endif
+ do jl=secseg,seg4
+ nl=nl+1
+ d0ij=dsqrt((cref(1,jl)-cref(1,il))**2+
+ & (cref(2,jl)-cref(2,il))**2+
+ & (cref(3,jl)-cref(3,il))**2)
+ dij=dist(il,jl)
+ qqij = dexp(-0.5d0*((dij-d0ij)/(sigm(d0ij)))**2)
+ if (itype(il).ne.10 .or. itype(jl).ne.10) then
+ nl=nl+1
+ d0ijCM=dsqrt(
+ & (cref(1,jl+nres)-cref(1,il+nres))**2+
+ & (cref(2,jl+nres)-cref(2,il+nres))**2+
+ & (cref(3,jl+nres)-cref(3,il+nres))**2)
+ dijCM=dist(il+nres,jl+nres)
+ qqijCM = dexp(-0.5d0*((dijCM-d0ijCM)/(sigm(d0ijCM)))**2)
+ endif
+ qq = qq+qqij+qqijCM
+ enddo
+ enddo
+ qq = qq/nl
+ endif
+ qwolynes=1.0d0-qq
+ return
+ end
+c-------------------------------------------------------------------
+ subroutine qwolynes_prim(seg1,seg2,flag,seg3,seg4)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.CHAIN'
+ include 'COMMON.INTERACT'
+ include 'COMMON.VAR'
+ include 'COMMON.MD'
+ integer i,j,jl,k,l,il,nl,seg1,seg2,seg3,seg4,
+ & secseg
+ integer nsep /3/
+ double precision dist
+ double precision dij,d0ij,dijCM,d0ijCM
+ logical lprn /.false./
+ logical flag
+ double precision sigm,x,sim,dd0,fac,ddqij
+ sigm(x)=0.25d0*x
+
+ do i=0,nres
+ do j=1,3
+ dqwol(j,i)=0.0d0
+ dxqwol(j,i)=0.0d0
+ enddo
+ enddo
+ nl=0
+ if(flag) then
+ do il=seg1+nsep,seg2
+ do jl=seg1,il-nsep
+ nl=nl+1
+ d0ij=dsqrt((cref(1,jl)-cref(1,il))**2+
+ & (cref(2,jl)-cref(2,il))**2+
+ & (cref(3,jl)-cref(3,il))**2)
+ dij=dist(il,jl)
+ sim = 1.0d0/sigm(d0ij)
+ sim = sim*sim
+ dd0 = dij-d0ij
+ fac = dd0*sim/dij*dexp(-0.5d0*dd0*dd0*sim)
+ do k=1,3
+ ddqij = (c(k,il)-c(k,jl))*fac
+ 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(
+ & (cref(1,jl+nres)-cref(1,il+nres))**2+
+ & (cref(2,jl+nres)-cref(2,il+nres))**2+
+ & (cref(3,jl+nres)-cref(3,il+nres))**2)
+ dijCM=dist(il+nres,jl+nres)
+ sim = 1.0d0/sigm(d0ijCM)
+ sim = sim*sim
+ dd0=dijCM-d0ijCM
+ fac=dd0*sim/dijCM*dexp(-0.5d0*dd0*dd0*sim)
+ do k=1,3
+ ddqij = (c(k,il+nres)-c(k,jl+nres))*fac
+ dxqwol(k,il)=dxqwol(k,il)+ddqij
+ dxqwol(k,jl)=dxqwol(k,jl)-ddqij
+ enddo
+ endif
+ enddo
+ enddo
+ else
+ do il=seg1,seg2
+ if((seg3-il).lt.3) then
+ secseg=il+3
+ else
+ secseg=seg3
+ endif
+ do jl=secseg,seg4
+ nl=nl+1
+ d0ij=dsqrt((cref(1,jl)-cref(1,il))**2+
+ & (cref(2,jl)-cref(2,il))**2+
+ & (cref(3,jl)-cref(3,il))**2)
+ dij=dist(il,jl)
+ sim = 1.0d0/sigm(d0ij)
+ sim = sim*sim
+ dd0 = dij-d0ij
+ fac = dd0*sim/dij*dexp(-0.5d0*dd0*dd0*sim)
+ do k=1,3
+ ddqij = (c(k,il)-c(k,jl))*fac
+ 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(
+ & (cref(1,jl+nres)-cref(1,il+nres))**2+
+ & (cref(2,jl+nres)-cref(2,il+nres))**2+
+ & (cref(3,jl+nres)-cref(3,il+nres))**2)
+ dijCM=dist(il+nres,jl+nres)
+ sim = 1.0d0/sigm(d0ijCM)
+ sim=sim*sim
+ dd0 = dijCM-d0ijCM
+ fac = dd0*sim/dijCM*dexp(-0.5d0*dd0*dd0*sim)
+ do k=1,3
+ ddqij = (c(k,il+nres)-c(k,jl+nres))*fac
+ dxqwol(k,il)=dxqwol(k,il)+ddqij
+ dxqwol(k,jl)=dxqwol(k,jl)-ddqij
+ enddo
+ endif
+ enddo
+ enddo
+ endif
+ do i=0,nres
+ do j=1,3
+ dqwol(j,i)=dqwol(j,i)/nl
+ dxqwol(j,i)=dxqwol(j,i)/nl
+ enddo
+ enddo
+ return
+ end
+c-------------------------------------------------------------------
+ subroutine qwol_num(seg1,seg2,flag,seg3,seg4)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.CHAIN'
+ include 'COMMON.INTERACT'
+ include 'COMMON.VAR'
+ integer seg1,seg2,seg3,seg4
+ logical flag
+ double precision qwolan(3,0:maxres),cdummy(3,0:maxres2),
+ & qwolxan(3,0:maxres),q1,q2
+ double precision delta /1.0d-10/
+ do i=0,nres
+ do j=1,3
+ q1=qwolynes(seg1,seg2,flag,seg3,seg4)
+ cdummy(j,i)=c(j,i)
+ c(j,i)=c(j,i)+delta
+ q2=qwolynes(seg1,seg2,flag,seg3,seg4)
+ qwolan(j,i)=(q2-q1)/delta
+ c(j,i)=cdummy(j,i)
+ enddo
+ enddo
+ do i=0,nres
+ do j=1,3
+ q1=qwolynes(seg1,seg2,flag,seg3,seg4)
+ cdummy(j,i+nres)=c(j,i+nres)
+ c(j,i+nres)=c(j,i+nres)+delta
+ q2=qwolynes(seg1,seg2,flag,seg3,seg4)
+ qwolxan(j,i)=(q2-q1)/delta
+ c(j,i+nres)=cdummy(j,i+nres)
+ enddo
+ enddo
+c write(iout,*) "Numerical Q carteisan gradients backbone: "
+c do i=0,nct
+c write(iout,'(i5,3e15.5)') i, (qwolan(j,i),j=1,3)
+c enddo
+c write(iout,*) "Numerical Q carteisan gradients side-chain: "
+c do i=0,nct
+c write(iout,'(i5,3e15.5)') i, (qwolxan(j,i),j=1,3)
+c enddo
+ 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