X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-M%2Fenergy_p_new_barrier.F;h=9414f1cd19507191665a8b05e510413c18352696;hb=3a74889c3f6514588d40ec326a78f33b3663f5be;hp=662a7967e51b00f3828fcd11d1072fa178de4309;hpb=c9dddd7cdc2cbf55f6fa8df9404eea7ecfe59605;p=unres.git diff --git a/source/unres/src_MD-M/energy_p_new_barrier.F b/source/unres/src_MD-M/energy_p_new_barrier.F index 662a796..9414f1c 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -190,9 +190,10 @@ C C Calculate the virtual-bond-angle energy. C if (wang.gt.0d0) then - call ebend(ebe) + call ebend(ebe,ethetacnstr) else ebe=0 + ethetacnstr=0 endif c print *,"Processor",myrank," computed UB" C @@ -301,6 +302,9 @@ C energia(17)=estr energia(20)=Uconst+Uconst_back energia(21)=esccor +C energia(22)=eliptrans (the energy for lipid transfere implemented in lipid branch) +C energia(23)= ... (energy for AFM, steered molecular dynamics) + energia(24)=ethetacnstr c Here are the energies showed per procesor if the are more processors c per molecule then we sum it up in sum_energy subroutine c print *," Processor",myrank," calls SUM_ENERGY" @@ -392,20 +396,22 @@ cMS$ATTRIBUTES C :: proc_proc estr=energia(17) Uconst=energia(20) esccor=energia(21) + ethetacnstr=energia(24) + #ifdef SPLITELE etot=wsc*evdw+wscp*evdw2+welec*ees+wvdwpp*evdw1 & +wang*ebe+wtor*etors+wscloc*escloc & +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5 & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3 & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d - & +wbond*estr+Uconst+wsccor*esccor + & +wbond*estr+Uconst+wsccor*esccor+ethetacnstr #else etot=wsc*evdw+wscp*evdw2+welec*(ees+evdw1) & +wang*ebe+wtor*etors+wscloc*escloc & +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5 & +wcorr6*ecorr6+wturn4*eello_turn4+wturn3*eello_turn3 & +wturn6*eturn6+wel_loc*eel_loc+edihcnstr+wtor_d*etors_d - & +wbond*estr+Uconst+wsccor*esccor + & +wbond*estr+Uconst+wsccor*esccor+ethetacnstr #endif energia(0)=etot c detecting NaNQ @@ -970,6 +976,7 @@ C------------------------------------------------------------------------ estr=energia(17) Uconst=energia(20) esccor=energia(21) + ethetacnstr=energia(24) #ifdef SPLITELE write (iout,10) evdw,wsc,evdw2,wscp,ees,welec,evdw1,wvdwpp, & estr,wbond,ebe,wang, @@ -977,7 +984,8 @@ C------------------------------------------------------------------------ & ecorr,wcorr, & ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3, & eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccor, - & edihcnstr,ebr*nss, + & edihcnstr, + & ethetacnstr,ebr*nss, & Uconst,etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ @@ -1000,6 +1008,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)'/ & 'UCONST= ',1pE16.6,' (Constraint energy)'/ & 'ETOT= ',1pE16.6,' (total)') @@ -1010,6 +1019,7 @@ C------------------------------------------------------------------------ & ecorr,wcorr, & ecorr5,wcorr5,ecorr6,wcorr6,eel_loc,wel_loc,eello_turn3,wturn3, & eello_turn4,wturn4,eello_turn6,wturn6,esccor,wsccro,edihcnstr, + & ethetacnstr, & ebr*nss,Uconst,etot 10 format (/'Virtual-chain energies:'// & 'EVDW= ',1pE16.6,' WEIGHT=',1pD16.6,' (SC-SC)'/ @@ -1031,6 +1041,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)'/ & 'UCONST=',1pE16.6,' (Constraint energy)'/ & 'ETOT= ',1pE16.6,' (total)') @@ -4400,7 +4411,7 @@ c 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. @@ -4417,6 +4428,7 @@ C include 'COMMON.NAMES' include 'COMMON.FFIELD' include 'COMMON.CONTROL' + 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 @@ -4526,6 +4538,34 @@ C Derivatives of the "mean" values in gamma1 and gamma2. if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*E_tc*dthetg2 gloc(nphi+i-2,icg)=wang*(E_theta+E_tc*dthett) enddo + ethetacnstr=0.0d0 +C print *,ithetaconstr_start,ithetaconstr_end,"TU" + do i=ithetaconstr_start,ithetaconstr_end + 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=ethetcnstr+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=ethetcnstr+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 + if (energy_dec) then + write (iout,'(a6,2i5,4f8.3,2e14.5)') "ethetc", + & i,itheta,rad2deg*thetiii, + & rad2deg*theta_constr0(i), rad2deg*theta_drange(i), + & rad2deg*difi,0.25d0*for_thet_constr(i)*difi**4, + & gloc(itheta+nphi-2,icg) + endif + enddo + C Ufff.... We've done all this!!! return end @@ -4639,7 +4679,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. @@ -4658,6 +4698,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), @@ -4844,6 +4885,35 @@ c lprn1=.false. if (i.lt.nres) gloc(i-2,icg)=gloc(i-2,icg)+wang*dephii1 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=ithetaconstr_start,ithetaconstr_end + 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 + if (energy_dec) then + write (iout,'(a6,2i5,4f8.3,2e14.5)') "ethetc", + & i,itheta,rad2deg*thetiii, + & rad2deg*theta_constr0(i), rad2deg*theta_drange(i), + & rad2deg*difi,0.25d0*for_thet_constr(i)*difi**4, + & gloc(itheta+nphi-2,icg) + endif + enddo + return end #endif @@ -5663,12 +5733,12 @@ c write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg) difi=phii-phi0(i) if (difi.gt.drange(i)) then difi=difi-drange(i) - edihcnstr=edihcnstr+0.25d0*ftors*difi**4 - gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3 + edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4 + gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3 else if (difi.lt.-drange(i)) then difi=difi+drange(i) - edihcnstr=edihcnstr+0.25d0*ftors*difi**4 - gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3 + edihcnstr=edihcnstr+0.25d0*ftors(i)**difi**4 + gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3 endif ! write (iout,'(2i5,2f8.3,2e14.5)') i,itori,rad2deg*phii, ! & rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg) @@ -5767,12 +5837,12 @@ c do i=1,ndih_constr difi=pinorm(phii-phi0(i)) if (difi.gt.drange(i)) then difi=difi-drange(i) - edihcnstr=edihcnstr+0.25d0*ftors*difi**4 - gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3 + edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4 + gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3 else if (difi.lt.-drange(i)) then difi=difi+drange(i) - edihcnstr=edihcnstr+0.25d0*ftors*difi**4 - gloc(itori-3,icg)=gloc(itori-3,icg)+ftors*difi**3 + edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4 + gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3 else difi=0.0 endif @@ -5780,7 +5850,7 @@ c do i=1,ndih_constr write (iout,'(a6,2i5,4f8.3,2e14.5)') "edihc", & i,itori,rad2deg*phii, & rad2deg*phi0(i), rad2deg*drange(i), - & rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg) + & rad2deg*difi,0.25d0*ftors(i)*difi**4,gloc(itori-3,icg) endif enddo cd write (iout,*) 'edihcnstr',edihcnstr