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.
etot=wsc*(evdw+fact(6)*evdw_t)+wscp*evdw2+welec*fact(1)*ees
& +wvdwpp*evdw1
& +wang*ebe+wtor*fact(1)*etors+wscloc*escloc
- & +wstrain*ehpb+nss*ebr+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5
+ & +wstrain*ehpb+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5
& +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)
& +wang*ebe+wtor*fact(1)*etors+wscloc*escloc
- & +wstrain*ehpb+nss*ebr+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5
+ & +wstrain*ehpb+wcorr*fact(3)*ecorr+wcorr5*fact(4)*ecorr5
& +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
IF (dyn_ss_mask(i).and.dyn_ss_mask(j)) THEN
call dyn_ssbond_ene(i,j,evdwij)
evdw=evdw+evdwij
- write (iout,'(a6,2i5,0pf7.3,a3)')
- & 'evdw',i,j,evdwij,' ss'
+C write (iout,'(a6,2i5,0pf7.3,a3,2f10.3)')
+C & 'evdw',i,j,evdwij,' ss',evdw,evdw_t
C triple bond artifac removal
do k=j+1,iend(i,iint)
C search over all next residues
C call the energy function that removes the artifical triple disulfide
C bond the soubroutine is located in ssMD.F
evdw=evdw+evdwij
-C if (energy_dec) write (iout,'(a6,2i5,0pf7.3,a3)')
-C & 'evdw',i,j,evdwij,'tss'
+C write (iout,'(a6,2i5,0pf7.3,a3,2f10.3)')
+C & 'evdw',i,j,evdwij,'tss',evdw,evdw_t
endif!dyn_ss_mask(k)
enddo! k
ELSE
C Calculate angular part of the gradient.
call sc_grad
endif
+C write(iout,*) "partial sum", evdw, evdw_t
ENDIF ! dyn_ss
enddo ! j
enddo ! iint
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
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
- else
-C Calculate the distance between the two points and its difference from the
-C target distance.
- dd=dist(ii,jj)
- rdis=dd-dhpb(i)
+ 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))
+ 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
+ dd=dist(ii,jj)
+ 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,*) "beta reg",dd,waga*rdis*rdis
+C
+C Evaluate gradient.
+C
+ fac=waga*rdis/dd
+ 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
+ do j=1,3
+ ghpbx(j,iii)=ghpbx(j,iii)-ggg(j)
+ ghpbx(j,jjj)=ghpbx(j,jjj)+ggg(j)
+ enddo
+ do k=1,3
+ ghpbc(k,jjj)=ghpbc(k,jjj)+ggg(k)
+ ghpbc(k,iii)=ghpbc(k,iii)-ggg(k)
+ enddo
+ 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)
+ waga=forcon(i)
C Calculate the contribution to energy.
- ehpb=ehpb+waga*rdis*rdis
+ ehpb=ehpb+waga*rdis*rdis
+c write (iout,*) "alpha reg",dd,waga*rdis*rdis
C
C Evaluate gradient.
C
- fac=waga*rdis/dd
-cd print *,'i=',i,' ii=',ii,' jj=',jj,' dhpb=',dhpb(i),' dd=',dd,
-cd & ' waga=',waga,' fac=',fac
+ fac=waga*rdis/dd
+ endif
+ endif
+
do j=1,3
ggg(j)=fac*(c(j,jj)-c(j,ii))
enddo
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
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
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)
+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
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
esccor=esccor+v1ij*cosphi+v2ij*sinphi
gloci=gloci+j*(v2ij*cosphi-v1ij*sinphi)
enddo
+C write (iout,*)"EBACK_SC_COR",esccor,i
c write (iout,*) "EBACK_SC_COR",i,v1ij*cosphi+v2ij*sinphi,intertyp,
c & nterm_sccor(isccori,isccori1),isccori,isccori1
c gloc_sc(intertyp,i-3,icg)=gloc_sc(intertyp,i-3,icg)+wsccor*gloci
projects / unres.git / blobdiff