X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fsrc-M%2Fenergy_p_new.F;h=44e5cba1dbc6e633cce1abf56814388f34e960f4;hb=7baaf36efbc265e029ec93a4bb1adfadfd55d7bf;hp=33e670a803d1746763e2f7ed284336100f58f76f;hpb=c65d296f3cdc480d3737eae335a510332b4f0dc8;p=unres.git diff --git a/source/wham/src-M/energy_p_new.F b/source/wham/src-M/energy_p_new.F index 33e670a..44e5cba 100644 --- a/source/wham/src-M/energy_p_new.F +++ b/source/wham/src-M/energy_p_new.F @@ -67,7 +67,7 @@ cd print *,'EHPB exitted succesfully.' C C Calculate the virtual-bond-angle energy. C - call ebend(ebe) + call ebend(ebe,ethetacnstr) cd print *,'Bend energy finished.' C C Calculate the SC local energy. @@ -111,7 +111,7 @@ c write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t & +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4 & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6 & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d - & +wbond*estr+wsccor*fact(1)*esccor + & +wbond*estr+wsccor*fact(1)*esccor+ethetacnstr #else etot=wsc*(evdw+fact(6)*evdw_t)+wscp*evdw2 & +welec*fact(1)*(ees+evdw1) @@ -120,7 +120,7 @@ c write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t & +wcorr6*fact(5)*ecorr6+wturn4*fact(3)*eello_turn4 & +wturn3*fact(2)*eello_turn3+wturn6*fact(5)*eturn6 & +wel_loc*fact(2)*eel_loc+edihcnstr+wtor_d*fact(2)*etors_d - & +wbond*estr+wsccor*fact(1)*esccor + & +wbond*estr+wsccor*fact(1)*esccor+ethetacnstr #endif energia(0)=etot energia(1)=evdw @@ -154,6 +154,7 @@ c write (iout,*) "ft(6)",fact(6)," evdw",evdw," evdw_t",evdw_t energia(19)=esccor energia(20)=edihcnstr energia(21)=evdw_t + energia(24)=ethetacnstr c detecting NaNQ #ifdef ISNAN #ifdef AIX @@ -270,6 +271,7 @@ C------------------------------------------------------------------------ esccor=energia(19) edihcnstr=energia(20) estr=energia(18) + ethetacnstr=energia(24) #ifdef SPLITELE write (iout,10) evdw,wsc,evdw2,wscp,ees,welec*fact(1),evdw1, & wvdwpp, @@ -278,7 +280,7 @@ C------------------------------------------------------------------------ & ecorr,wcorr*fact(3),ecorr5,wcorr5*fact(4),ecorr6,wcorr6*fact(5), & eel_loc,wel_loc*fact(2),eello_turn3,wturn3*fact(2), & eello_turn4,wturn4*fact(3),eello_turn6,wturn6*fact(5), - & esccor,wsccor*fact(1),edihcnstr,ebr*nss,etot + & esccor,wsccor*fact(1),edihcnstr,ethetacnstr,ebr*nss,etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ @@ -300,6 +302,7 @@ C------------------------------------------------------------------------ & 'ETURN6=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 6th order)'/ & 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/ & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/ + & 'ETHETC= ',1pE16.6,' (valence angle constraints)'/ & 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ & 'ETOT= ',1pE16.6,' (total)') #else @@ -309,7 +312,7 @@ C------------------------------------------------------------------------ & ecorr6,wcorr6*fact(5),eel_loc,wel_loc*fact(2), & eello_turn3,wturn3*fact(2),eello_turn4,wturn4*fact(3), & eello_turn6,wturn6*fact(5),esccor*fact(1),wsccor, - & edihcnstr,ebr*nss,etot + & edihcnstr,ethetacnstr,ebr*nss,etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ & 'EVDW2= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-p)'/ @@ -330,6 +333,7 @@ C------------------------------------------------------------------------ & 'ETURN6=',1pE16.6,' WEIGHT=',1pD16.6,' (turns, 6th order)'/ & 'ESCCOR=',1pE16.6,' WEIGHT=',1pD16.6,' (backbone-rotamer corr)'/ & 'EDIHC= ',1pE16.6,' (dihedral angle constraints)'/ + & 'ETHETC= ',1pE16.6,' (valence angle constraints)'/ & 'ESS= ',1pE16.6,' (disulfide-bridge intrinsic energy)'/ & 'ETOT= ',1pE16.6,' (total)') #endif @@ -3265,7 +3269,7 @@ c & AKSC(j,iti),abond0(j,iti),u(j),j=1,nbi) end #ifdef CRYST_THETA C-------------------------------------------------------------------------- - subroutine ebend(etheta) + subroutine ebend(etheta,ethetacnstr) C C Evaluate the virtual-bond-angle energy given the virtual-bond dihedral C angles gamma and its derivatives in consecutive thetas and gammas. @@ -3282,6 +3286,7 @@ C include 'COMMON.IOUNITS' include 'COMMON.NAMES' include 'COMMON.FFIELD' + include 'COMMON.TORCNSTR' common /calcthet/ term1,term2,termm,diffak,ratak, & ak,aktc,termpre,termexp,sigc,sig0i,time11,time12,sigcsq, & delthe0,sig0inv,sigtc,sigsqtc,delthec,it @@ -3400,6 +3405,33 @@ c & rad2deg*phii,rad2deg*phii1,ethetai gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett) c 1215 continue enddo + ethetacnstr=0.0d0 +C print *,ithetaconstr_start,ithetaconstr_end,"TU" + do i=1,ntheta_constr + itheta=itheta_constr(i) + thetiii=theta(itheta) + difi=pinorm(thetiii-theta_constr0(i)) + if (difi.gt.theta_drange(i)) then + difi=difi-theta_drange(i) + ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4 + gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) + & +for_thet_constr(i)*difi**3 + else if (difi.lt.-drange(i)) then + difi=difi+drange(i) + ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4 + gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) + & +for_thet_constr(i)*difi**3 + else + difi=0.0 + endif +C if (energy_dec) then +C write (iout,'(a6,2i5,4f8.3,2e14.5)') "ethetc", +C & i,itheta,rad2deg*thetiii, +C & rad2deg*theta_constr0(i), rad2deg*theta_drange(i), +C & rad2deg*difi,0.25d0*for_thet_constr(i)*difi**4, +C & gloc(itheta+nphi-2,icg) + endif + enddo C Ufff.... We've done all this!!! return end @@ -3513,7 +3545,7 @@ C "Thank you" to MAPLE (probably spared one day of hand-differentiation). end #else C-------------------------------------------------------------------------- - subroutine ebend(etheta) + subroutine ebend(etheta,ethetacnstr) C C Evaluate the virtual-bond-angle energy given the virtual-bond dihedral C angles gamma and its derivatives in consecutive thetas and gammas. @@ -3533,6 +3565,7 @@ C include 'COMMON.NAMES' include 'COMMON.FFIELD' include 'COMMON.CONTROL' + include 'COMMON.TORCNSTR' double precision coskt(mmaxtheterm),sinkt(mmaxtheterm), & cosph1(maxsingle),sinph1(maxsingle),cosph2(maxsingle), & sinph2(maxsingle),cosph1ph2(maxdouble,maxdouble), @@ -3712,6 +3745,34 @@ c call flush(iout) c gloc(nphi+i-2,icg)=wang*dethetai gloc(nphi+i-2,icg)=gloc(nphi+i-2,icg)+wang*dethetai enddo +C now constrains + ethetacnstr=0.0d0 +C print *,ithetaconstr_start,ithetaconstr_end,"TU" + do i=1,ntheta_constr + itheta=itheta_constr(i) + thetiii=theta(itheta) + difi=pinorm(thetiii-theta_constr0(i)) + if (difi.gt.theta_drange(i)) then + difi=difi-theta_drange(i) + ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4 + gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) + & +for_thet_constr(i)*difi**3 + else if (difi.lt.-drange(i)) then + difi=difi+drange(i) + ethetacnstr=ethetacnstr+0.25d0*for_thet_constr(i)*difi**4 + gloc(itheta+nphi-2,icg)=gloc(itheta+nphi-2,icg) + & +for_thet_constr(i)*difi**3 + else + difi=0.0 + endif +C if (energy_dec) then +C write (iout,'(a6,2i5,4f8.3,2e14.5)') "ethetc", +C & i,itheta,rad2deg*thetiii, +C & rad2deg*theta_constr0(i), rad2deg*theta_drange(i), +C & rad2deg*difi,0.25d0*for_thet_constr(i)*difi**4, +C & gloc(itheta+nphi-2,icg) +C endif + enddo return end #endif