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=09d4fb10f3c623621964b33c806fe05678a95beb;hb=75c81b9dbe2f7e18e73a9d61ee02980790335164;hp=ec0aef92c147b72de51eb82ba118958ee2a6d49f;hpb=42b30ef61fe3a8077aee8292bd43e4de4e6b8522;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 ec0aef9..09d4fb1 100644 --- a/source/unres/src_MD-M/energy_p_new_barrier.F +++ b/source/unres/src_MD-M/energy_p_new_barrier.F @@ -193,6 +193,7 @@ C call ebend(ebe,ethetacnstr) else ebe=0 + ethetacnstr=0 endif c print *,"Processor",myrank," computed UB" C @@ -301,7 +302,9 @@ C energia(17)=estr energia(20)=Uconst+Uconst_back energia(21)=esccor -C energia(22)= +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" @@ -393,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 @@ -971,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, @@ -1461,9 +1467,11 @@ c write(iout,*) "PRZED ZWYKLE", evdwij c write(iout,*) "PO ZWYKLE", evdwij evdw=evdw+evdwij +c write(iout,*) "DISULFIDY:", i,j,evdwij if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)') & 'evdw',i,j,evdwij,' ss' C triple bond artifac removal +C MODIFIED j+1 to j+2 TO AVOID EBERGY BARRIER FOR X-Cys-Cys-X situations do k=j+1,iend(i,iint) C search over all next residues if (dyn_ss_mask(k)) then @@ -1473,6 +1481,7 @@ C write(iout,*) 'k=',k c write(iout,*) "PRZED TRI", evdwij evdwij_przed_tri=evdwij call triple_ssbond_ene(i,j,k,evdwij) +c write(iout,*) "TRISULFIDY:", i,j,k,evdwij c if(evdwij_przed_tri.ne.evdwij) then c write (iout,*) "TRI:", evdwij, evdwij_przed_tri c endif @@ -4405,7 +4414,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. @@ -4422,6 +4431,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 @@ -4525,18 +4535,26 @@ C Derivatives of the "mean" values in gamma1 and gamma2. & E_theta,E_tc) endif etheta=etheta+ethetai + if (energy_dec) write (iout,'(a6,i5,0pf7.3)') + & 'ebend',i,ethetai + if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*E_tc*dthetg1 + 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(itori) + thetiii=theta(itheta) difi=pinorm(thetiii-theta_constr0(i)) if (difi.gt.theta_drange(i)) then difi=difi-theta_drange(i) - ethetcnstr=ethetcnstr+0.25d0*for_thet_constr(i)*difi**4 + 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) - ethetcnstr=ethetcnstr+0.25d0*for_thet_constr(i)*difi**4 + 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 @@ -4550,12 +4568,7 @@ C Derivatives of the "mean" values in gamma1 and gamma2. & gloc(itheta+nphi-2,icg) endif enddo - if (energy_dec) write (iout,'(a6,i5,0pf7.3)') - & 'ebend',i,ethetai - if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*E_tc*dthetg1 - 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 + C Ufff.... We've done all this!!! return end @@ -4669,7 +4682,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. @@ -4688,6 +4701,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), @@ -4862,19 +4876,33 @@ C print *,"tu",cosph1(k),sinph1(k),cosph2(k),sinph2(k) enddo enddo 10 continue -C now we have the theta_constrains +c lprn1=.true. +C print *,ethetai + if (lprn1) + & write (iout,'(i2,3f8.1,9h ethetai ,f10.5)') + & i,theta(i)*rad2deg,phii*rad2deg, + & phii1*rad2deg,ethetai +c lprn1=.false. + etheta=etheta+ethetai + if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*dephii + 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(itori) + thetiii=theta(itheta) difi=pinorm(thetiii-theta_constr0(i)) if (difi.gt.theta_drange(i)) then difi=difi-theta_drange(i) - ethetcnstr=ethetcnstr+0.25d0*for_thet_constr(i)*difi**4 + 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) - ethetcnstr=ethetcnstr+0.25d0*for_thet_constr(i)*difi**4 + 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 @@ -4889,18 +4917,6 @@ C now we have the theta_constrains endif enddo -c lprn1=.true. -C print *,ethetai - if (lprn1) - & write (iout,'(i2,3f8.1,9h ethetai ,f10.5)') - & i,theta(i)*rad2deg,phii*rad2deg, - & phii1*rad2deg,ethetai -c lprn1=.false. - etheta=etheta+ethetai - if (i.gt.3) gloc(i-3,icg)=gloc(i-3,icg)+wang*dephii - 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 return end #endif