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.
& +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)
& +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
energia(19)=esccor
energia(20)=edihcnstr
energia(21)=evdw_t
+ energia(24)=ethetacnstr
c detecting NaNQ
#ifdef ISNAN
#ifdef AIX
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,
& 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)'/
& '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
& 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)'/
& '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
common /srutu/icall
integer icant
external icant
- logical energy_dec /.false./
+ logical energy_dec /.true./
c print *,'Entering EGB nnt=',nnt,' nct=',nct,' expon=',expon
evdw=0.0D0
evdw_t=0.0d0
include 'COMMON.DERIV'
include 'COMMON.VAR'
include 'COMMON.INTERACT'
+ include 'COMMON.CONTROL'
dimension ggg(3)
ehpb=0.0D0
cd print *,'edis: nhpb=',nhpb,' fbr=',fbr
endif
C 24/11/03 AL: SS bridges handled separately because of introducing a specific
C distance and angle dependent SS bond potential.
- if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and.
- & iabs(itype(jjj)).eq.1) then
+C if (ii.gt.nres .and. iabs(itype(iii)).eq.1 .and.
+C & iabs(itype(jjj)).eq.1) then
+C call ssbond_ene(iii,jjj,eij)
+C ehpb=ehpb+2*eij
+C else
+ 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
- else
+ 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 print *,"TUTU"
+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 !constr_dist.eq.11
+ 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,*) "beta nmr",
+c & dd,2*forcon(i)*gnmr1(dd,dhpb(i),dhpb1(i))
+ else !dhpb(i).gt.0.00
+
C Calculate the distance between the two points and its difference from the
C target distance.
dd=dist(ii,jj)
C Evaluate gradient.
C
fac=waga*rdis/dd
+ endif !dhpb(i).gt.0
+ endif
cd print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd,
cd & ' waga=',waga,' fac=',fac
do j=1,3
ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j)
enddo
endif
+ 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)
+C Calculate the contribution to energy.
+ ehpb=ehpb+waga*rdis*rdis
+c write (iout,*) "alpha reg",dd,waga*rdis*rdis
+C
+C Evaluate gradient.
+C
+ fac=waga*rdis/dd
+ endif
+ endif
+ do j=1,3
+ ggg(j)=fac*(c(j,jj)-c(j,ii))
+ enddo
+cd print '(i3,3(1pe14.5))',i,(ggg(j),j=1,3)
+C If this is a SC-SC distance, we need to calculate the contributions to the
+C Cartesian gradient in the SC vectors (ghpbx).
+ if (iii.lt.ii) then
+ do j=1,3
+ ghpbx(j,iii)=ghpbx(j,iii)-ggg(j)
+ ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j)
+ enddo
+ endif
do j=iii,jjj-1
do k=1,3
ghpbc(k,j)=ghpbc(k,j)+ggg(k)
enddo
endif
enddo
- ehpb=0.5D0*ehpb
+ if (constr_dist.ne.11) ehpb=0.5D0*ehpb
return
end
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.
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
c 1215 continue
enddo
C Ufff.... We've done all this!!!
+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
C---------------------------------------------------------------------------
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.
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),
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
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)
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