+++ /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'
- 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
projects / unres.git / blobdiff