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' include 'COMMON.MD' integer i,j,jl,k,l,il,kl,nl,np,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 qq = 0.0d0 nl=0 do i=0,nres do j=1,3 dqwol(j,i)=0.0d0 dxqwol(j,i)=0.0d0 enddo enddo if (lprn) then write (iout,*) "seg1",seg1," seg2",seg2," seg3",seg3," seg4",seg4, & " flag",flag call flush(iout) endif if (flag) then do il=seg1+nsep,seg2 do jl=seg1,il-nsep nl=nl+1 if (itype(il).ne.10) then ilnres=il+nres else ilnres=il endif if (itype(jl).ne.10) then jlnres=jl+nres else jlnres=jl endif qqijCM = qcontrib(il,jl,ilnres,jlnres) qq = qq+qqijCM if (lprn) then write (iout,*) "qqijCM",qqijCM call flush(iout) endif enddo enddo if (lprn) then write (iout,*) "nl",nl," qq",qq call flush(iout) endif 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 if (itype(il).ne.10) then ilnres=il+nres else ilnres=il endif if (itype(jl).ne.10) then jlnres=jl+nres else jlnres=jl endif qqijCM = qcontrib(il,jl,ilnres,jlnres) qq = qq+qqijCM if (lprn) then write (iout,*) "qqijCM",qqijCM call flush(iout) endif enddo enddo endif qq = qq/nl qwolynes=1.0d0-qq 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' include 'COMMON.MD' 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-7/ write (iout,*) "seg1",seg1," seg2",seg2," seg3",seg3," seg4",seg4 write(iout,*) "dQ/dc backbone " do i=0,nres write(iout,'(i5,3e15.5)') i, (dqwol(j,i),j=1,3) enddo write(iout,*) "dQ/dX side chain " do i=1,nres write(iout,'(i5,3e15.5)') i,(dxqwol(j,i),j=1,3) enddo do i=1,nres do j=1,3 cdummy(j,i)=c(j,i) c(j,i)=c(j,i)-delta q1=qwolynes(seg1,seg2,flag,seg3,seg4) c(j,i)=cdummy(j,i)+delta q2=qwolynes(seg1,seg2,flag,seg3,seg4) qwolan(j,i)=0.5d0*(q2-q1)/delta c(j,i)=cdummy(j,i) c write (iout,*) "i",i," j",j," q1",q1," a2",q2 enddo enddo do i=1,nres do j=1,3 cdummy(j,i+nres)=c(j,i+nres) c(j,i+nres)=c(j,i+nres)-delta q1=qwolynes(seg1,seg2,flag,seg3,seg4) c(j,i+nres)=cdummy(j,i+nres)+delta q2=qwolynes(seg1,seg2,flag,seg3,seg4) qwolxan(j,i)=0.5d0*(q2-q1)/delta c(j,i+nres)=cdummy(j,i+nres) enddo enddo write(iout,*) "Numerical Q cartesian gradients backbone: " do i=0,nres write(iout,'(i5,3e15.5)') i, (qwolan(j,i),j=1,3) enddo write(iout,*) "Numerical Q cartesian gradients side-chain: " do i=0,nres write(iout,'(i5,3e15.5)') i, (qwolxan(j,i),j=1,3) 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 Calculating the derivatives of Q with respect to cartesian coordinates 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 c write (iout,*) "Calling qwol_num" c call qwol_num(ifrag(1,i),ifrag(2,i),.true.,idummy,idummy) 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 Calculating dQ/dXi 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/dc backbone " c do ii=0,nres c write(iout,'(i5,3e15.5)') ii, (duconst(j,ii),j=1,3) c enddo c write(iout,*) "dU/dX side chain " c do ii=1,nres c write(iout,'(i5,3e15.5)') ii,(duxconst(j,ii),j=1,3) c 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,(dudxconst(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 do j=1,3 duxcartan(j,i)=0.0d0 enddo 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 write(iout,*) "Numerical dUconst/ddx side-chain " do ii=1,nres write(iout,'(i5,3e15.5)') ii,(dUxcartan(j,ii),j=1,3) enddo return end c--------------------------------------------------------------------------- double precision function qcontrib(il,jl,il1,jl1) implicit none include 'DIMENSIONS' include 'COMMON.IOUNITS' include 'COMMON.CHAIN' include 'COMMON.INTERACT' include 'COMMON.MD' integer i,j,k,il,jl,il1,jl1,nd double precision dist external dist double precision dij1,dij2,dij3,dij4,d0ij1,d0ij2,d0ij3,d0ij4,fac, & fac1,ddave,ssij,ddqij logical lprn /.false./ d0ij1=dsqrt((cref(1,jl)-cref(1,il))**2+ & (cref(2,jl)-cref(2,il))**2+ & (cref(3,jl)-cref(3,il))**2) dij1=dist(il,jl) ddave=(dij1-d0ij1)**2 nd=1 if (jl1.ne.jl) then d0ij2=dsqrt((cref(1,jl1)-cref(1,il))**2+ & (cref(2,jl1)-cref(2,il))**2+ & (cref(3,jl1)-cref(3,il))**2) dij2=dist(il,jl1) ddave=ddave+(dij2-d0ij2)**2 nd=nd+1 endif if (il1.ne.il) then d0ij3=dsqrt((cref(1,jl)-cref(1,il1))**2+ & (cref(2,jl)-cref(2,il1))**2+ & (cref(3,jl)-cref(3,il1))**2) dij3=dist(il1,jl) ddave=ddave+(dij3-d0ij3)**2 nd=nd+1 endif if (il1.ne.il .and. jl1.ne.jl) then d0ij4=dsqrt((cref(1,jl1)-cref(1,il1))**2+ & (cref(2,jl1)-cref(2,il1))**2+ & (cref(3,jl1)-cref(3,il1))**2) dij4=dist(il1,jl1) ddave=ddave+(dij4-d0ij4)**2 nd=nd+1 endif ddave=ddave/nd if (lprn) then write (iout,*) "il",il," jl",jl, & " itype",itype(il),itype(jl)," nd",nd write (iout,*)"d0ij",d0ij1,d0ij2,d0ij3,d0ij4, & " dij",dij1,dij2,dij3,dij4," ddave",ddave call flush(iout) endif c ssij = (0.25d0*d0ij1)**2 if (il.ne.il1 .and. jl.ne.jl1) then ssij = 16.0d0/(d0ij1*d0ij4) else ssij = 16.0d0/(d0ij1*d0ij1) endif qcontrib = dexp(-0.5d0*ddave*ssij) c Compute gradient fac1 = qcontrib*ssij/nd fac = fac1*(dij1-d0ij1)/dij1 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 (jl1.ne.jl) then fac = fac1*(dij2-d0ij2)/dij2 do k=1,3 ddqij = (c(k,il)-c(k,jl1))*fac dqwol(k,il)=dqwol(k,il)+ddqij dxqwol(k,jl)=dxqwol(k,jl)-ddqij enddo endif if (il1.ne.il) then fac = fac1*(dij3-d0ij3)/dij3 do k=1,3 ddqij = (c(k,il1)-c(k,jl))*fac dxqwol(k,il)=dxqwol(k,il)+ddqij dqwol(k,jl)=dqwol(k,jl)-ddqij enddo endif if (il1.ne.il .and. jl1.ne.jl) then fac = fac1*(dij4-d0ij4)/dij4 do k=1,3 ddqij = (c(k,il1)-c(k,jl1))*fac dxqwol(k,il)=dxqwol(k,il)+ddqij dxqwol(k,jl)=dxqwol(k,jl)-ddqij enddo endif return end