X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2Fq_measure.F;h=5077b40adb25ad7a1595d1d43ff39007b9146494;hb=34d3ad3987785642be58fb2f26557d3314215577;hp=8f12dc1087770e12890dd7b84ef460f925ab5001;hpb=f690e8b70bab14132839afebf080d4a28363b226;p=unres.git diff --git a/source/unres/src_MD-M/q_measure.F b/source/unres/src_MD-M/q_measure.F index 8f12dc1..5077b40 100644 --- a/source/unres/src_MD-M/q_measure.F +++ b/source/unres/src_MD-M/q_measure.F @@ -14,6 +14,10 @@ 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 @@ -26,6 +30,7 @@ & (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( @@ -34,10 +39,15 @@ & (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 @@ -53,6 +63,7 @@ & (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( @@ -62,14 +73,16 @@ 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------------------------------------------------------------------- @@ -90,6 +103,11 @@ 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 @@ -115,7 +133,6 @@ c------------------------------------------------------------------- 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( @@ -133,6 +150,10 @@ c------------------------------------------------------------------- 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 @@ -178,6 +199,9 @@ c------------------------------------------------------------------- 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 @@ -230,262 +254,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---------------------------------------------------------------------------