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=49b4d93d38312bdd6e5e9c6bb845431396c9b0f5;hp=6e19387514503b4677ce038146e76accb3b7fe71;hpb=c45326004a77063de99e919c46db1a8bba1d3eba;p=unres.git diff --git a/source/wham/src-M/energy_p_new.F b/source/wham/src-M/energy_p_new.F index 6e19387..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 @@ -2947,10 +2951,13 @@ C include 'COMMON.DERIV' include 'COMMON.VAR' include 'COMMON.INTERACT' + include 'COMMON.CONTROL' + include 'COMMON.IOUNITS' dimension ggg(3) ehpb=0.0D0 cd print *,'edis: nhpb=',nhpb,' fbr=',fbr cd print *,'link_start=',link_start,' link_end=',link_end +C write(iout,*) link_end, "link_end" if (link_end.eq.0) return do i=link_start,link_end C If ihpb(i) and jhpb(i) > NRES, this is a SC-SC distance, otherwise a @@ -2969,21 +2976,29 @@ C 24/11/03 AL: SS bridges handled separately because of introducing a specific C distance and angle dependent SS bond potential. C if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and. C & iabs(itype(jjj)).eq.1) then - +C write(iout,*) constr_dist,"const" if (.not.dyn_ss .and. i.le.nss) then if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and. & iabs(itype(jjj)).eq.1) then call ssbond_ene(iii,jjj,eij) ehpb=ehpb+2*eij - endif + endif !ii.gt.neres else if (ii.gt.nres .and. jj.gt.nres) then c Restraints from contact prediction dd=dist(ii,jj) if (constr_dist.eq.11) then +C ehpb=ehpb+fordepth(i)**4.0d0 +C & *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)) ehpb=ehpb+fordepth(i)**4.0d0 & *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)) fac=fordepth(i)**4.0d0 & *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd +C write (iout,'(a6,2i5,3f8.3)') "edisl",ii,jj, +C & ehpb,fordepth(i),dd +C write(iout,*) ehpb,"atu?" +C ehpb,"tu?" +C fac=fordepth(i)**4.0d0 +C & *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd else if (dhpb1(i).gt.0.0d0) then ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) @@ -3002,8 +3017,8 @@ C C Evaluate gradient. C fac=waga*rdis/dd - endif - endif + endif !end dhpb1(i).gt.0 + endif !end const_dist=11 do j=1,3 ggg(j)=fac*(c(j,jj)-c(j,ii)) enddo @@ -3015,7 +3030,29 @@ C ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k) ghpbc(k,iii)=ghpbc(k,iii)-ggg(k) enddo - else + else !ii.gt.nres +C write(iout,*) "before" + dd=dist(ii,jj) +C write(iout,*) "after",dd + if (constr_dist.eq.11) then + ehpb=ehpb+fordepth(i)**4.0d0 + & *rlornmr1(dd,dhpb(i),dhpb1(i),forcon(i)) + fac=fordepth(i)**4.0d0 + & *rlornmr1prim(dd,dhpb(i),dhpb1(i),forcon(i))/dd +C ehpb=ehpb+fordepth(i)**4*rlornmr1(dd,dhpb(i),dhpb1(i)) +C fac=fordepth(i)**4*rlornmr1prim(dd,dhpb(i),dhpb1(i))/dd +C print *,ehpb,"tu?" +C write(iout,*) ehpb,"btu?", +C & dd,dhpb(i),dhpb1(i),fordepth(i),forcon(i) +C write (iout,'(a6,2i5,3f8.3)') "edisl",ii,jj, +C & ehpb,fordepth(i),dd + else + if (dhpb1(i).gt.0.0d0) then + ehpb=ehpb+2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) + fac=forcon(i)*gnmr1prim(dd,dhpb(i),dhpb1(i))/dd +c write (iout,*) "alph nmr", +c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i)) + else rdis=dd-dhpb(i) C Get the force constant corresponding to this distance. waga=forcon(i) @@ -3232,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. @@ -3249,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 @@ -3367,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 @@ -3480,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. @@ -3500,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), @@ -3679,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 @@ -4445,15 +4539,16 @@ 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) +C write (iout,'(a6,2i5,2f8.3,2e14.5)') "edih", +C & i,itori,rad2deg*phii, +C & rad2deg*difi,0.25d0*ftors(i)*difi**4,gloc(itori-3,icg) enddo ! write (iout,*) 'edihcnstr',edihcnstr return @@ -4540,21 +4635,24 @@ c write (iout,*) 'i=',i,' gloc=',gloc(i-3,icg) edihi=0.0d0 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 - edihi=0.25d0*ftors*difi**4 + edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4 + gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3 + edihi=0.25d0*ftors(i)*difi**4 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 - edihi=0.25d0*ftors*difi**4 + edihcnstr=edihcnstr+0.25d0*ftors(i)*difi**4 + gloc(itori-3,icg)=gloc(itori-3,icg)+ftors(i)*difi**3 + edihi=0.25d0*ftors(i)*difi**4 else difi=0.0d0 endif + write (iout,'(a6,2i5,2f8.3,2e14.5)') "edih", + & i,itori,rad2deg*phii, + & rad2deg*difi,0.25d0*ftors(i)*difi**4 c write (iout,'(2i5,4f10.5,e15.5)') i,itori,phii,phi0(i),difi, c & drange(i),edihi ! write (iout,'(2i5,2f8.3,2e14.5)') i,itori,rad2deg*phii, -! & rad2deg*difi,0.25d0*ftors*difi**4,gloc(itori-3,icg) +! & rad2deg*difi,0.25d0*ftors(i)*difi**4,gloc(itori-3,icg) enddo ! write (iout,*) 'edihcnstr',edihcnstr return