X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M-newcorr%2Fq_measure3.F;fp=source%2Funres%2Fsrc_MD-M-newcorr%2Fq_measure3.F;h=f0a030e65d638d5d2cdcebb5d7a94964c31d2211;hb=d101c97dea752458d76055fdbae49c26fff03c1f;hp=0000000000000000000000000000000000000000;hpb=325eda160c9ad2982501e091ca40606a29043712;p=unres.git diff --git a/source/unres/src_MD-M-newcorr/q_measure3.F b/source/unres/src_MD-M-newcorr/q_measure3.F new file mode 100644 index 0000000..f0a030e --- /dev/null +++ b/source/unres/src_MD-M-newcorr/q_measure3.F @@ -0,0 +1,529 @@ + 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,iset),ifrag(2,i,iset),.true.,idummy,idummy) + enddo +c stop + 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' + include 'COMMON.LOCAL' + integer i,j,k,il,jl,il1,jl1,nd,itl,jtl + double precision dist + external dist + double precision dij,dij1,d0ij,d0ij1,om1,om2,om12,om10,om20,om120 + & ,fac,fac1,ddave,ssij,ddqij,d0ii1,d0jj1,rij,eom1,eom2,eom12 + double precision u(3),v(3),er(3),er0(3),dcosom1(3),dcosom2(3), + & aux1,aux2 + double precision scalar + external scalar + logical lprn /.false./ + if (lprn) write (iout,*) "il",il," jl",jl," il1",il1," jl1",jl1 + 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) + dij1=dist(il1,jl1) + do i=1,3 + er(i)=(c(i,jl1)-c(i,il1))/dij1 + enddo + do i=1,3 + er0(i)=cref(i,jl1)-cref(i,il1) + enddo + d0ij1=dsqrt(scalar(er0,er0)) + do i=1,3 + er0(i)=er0(i)/d0ij1 + enddo + if (il.ne.il1 .or. jl.ne.jl1) then + ddave=0.5d0*((dij-d0ij)**2+(dij1-d0ij1)**2) + nd=2 + else + ddave=(dij-d0ij)**2 + nd=1 + endif + if (il.ne.il1) then + do i=1,3 + u(i)=cref(i,il1)-cref(i,il) + enddo + d0ii1=dsqrt(scalar(u,u)) + do i=1,3 + u(i)=u(i)/d0ii1 + enddo + if (lprn) then + write (iout,*) "u",(u(i),i=1,3) + write (iout,*) "er0",(er0(i),i=1,3) + om10=scalar(er0,u) + om1=scalar(er,dc_norm(1,il1)) + write (iout,*) "om10",om10," om1",om1 + endif + else + om1=0.0d0 + om10=0.0d0 + endif + if (jl.ne.jl1) then + do i=1,3 + v(i)=cref(i,jl1)-cref(i,jl) + enddo + d0jj1=dsqrt(scalar(v,v)) + do i=1,3 + v(i)=v(i)/d0jj1 + enddo + if (lprn) then + write (iout,*) "v",(v(i),i=1,3) + write (iout,*) "er0",(er0(i),i=1,3) + om20=scalar(er,v) + om2=scalar(er,dc_norm(1,jl1)) + write (iout,*) "om20",om20," om2",om2 + endif + else + om2=0.0d0 + om20=0.0d0 + endif + if (il.ne.il1 .and. jl.ne.jl1) then + om120=scalar(u,v) + om12=scalar(dc_norm(1,il1),dc_norm(1,jl1)) + else + om12=0.0d0 + om120=0.0d0 + endif + if (lprn) then + write (iout,*) "il",il," jl",jl,itype(il),itype(jl) + write (iout,*)"d0ij",d0ij," om10",om10," om20",om20, + & " om120",om120, + & " dij",dij," om1",om1," om2",om2," om12",om12 + call flush(iout) + endif + ssij = 16.0d0/(d0ij*d0ij) + qcontrib = dexp(-0.5d0*(ddave*ssij+((om1-om10)**2 + & +(om2-om20)**2+(om12-om120)**2))) + if (lprn) write (iout,*) "ssij",ssij," qcontrib",qcontrib +c qcontrib = dexp(-0.5d0*(ddave*ssij)+(om1-om10)**2+(om2-om20)**2) +c qcontrib = dexp(-0.5d0*(ddave*ssij)) +c Compute gradient - radial component + fac1 = qcontrib*ssij/nd + fac = fac1*(dij-d0ij)/dij + 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 (il1.ne.il .or. jl1.ne.jl) then + fac = fac1*(dij1-d0ij1)/dij1 + do k=1,3 + ddqij = (c(k,il1)-c(k,jl1))*fac + if (il1.ne.il) then + dxqwol(k,il)=dxqwol(k,il)+ddqij + else + dqwol(k,il)=dqwol(k,il)+ddqij + endif + if (jl1.ne.jl) then + dxqwol(k,jl)=dxqwol(k,jl)-ddqij + else + dqwol(k,jl)=dqwol(k,jl)-ddqij + endif + enddo + endif +c return +c Orientational contributions + rij=1.0d0/dij1 + eom1=qcontrib*(om1-om10) + eom2=qcontrib*(om2-om20) + eom12=qcontrib*(om12-om120) + do k=1,3 + dcosom1(k)=rij*(dc_norm(k,il1)-om1*er(k)) + dcosom2(k)=rij*(dc_norm(k,jl1)-om2*er(k)) + enddo + do k=1,3 + ddqij=eom1*dcosom1(k)+eom2*dcosom2(k) + aux1=(eom12*(dc_norm(k,jl1)-om12*dc_norm(k,il1)) + & +eom1*(er(k)-om1*dc_norm(k,il1)))*vbld_inv(il1) + aux2=(eom12*(dc_norm(k,il1)-om12*dc_norm(k,jl1)) + & +eom2*(er(k)-om2*dc_norm(k,jl1)))*vbld_inv(jl1) + dqwol(k,il)=dqwol(k,il)-ddqij-aux1 + dqwol(k,jl)=dqwol(k,jl)+ddqij-aux2 + dxqwol(k,il)=dxqwol(k,il)-ddqij+aux1 +c & +(eom12*(dc_norm(k,jl1)-om12*dc_norm(k,il1)) +c & +eom1*(er(k)-om1*dc_norm(k,il1)))*vbld_inv(il1) + dxqwol(k,jl)=dxqwol(k,jl)+ddqij+aux2 +c & +(eom12*(dc_norm(k,il1)-om12*dc_norm(k,jl1)) +c & +eom2*(er(k)-om2*dc_norm(k,jl1)))*vbld_inv(jl1) + enddo + return + end