parameter (NGridT=400)
integer MaxBinRms,MaxBinRgy
parameter (MaxBinRms=100,MaxBinRgy=100)
- integer MaxHdim
-c parameter (MaxHdim=200000)
- parameter (MaxHdim=200)
+c integer MaxHdim
+c parameter (MaxHdim=200)
integer maxinde
parameter (maxinde=200)
#ifdef MPI
include "COMMON.SBRIDGE"
include "COMMON.PROT"
include "COMMON.ENEPS"
+ include "COMMON.SHIELD"
integer MaxPoint,MaxPointProc
parameter (MaxPoint=MaxStr,
& MaxPointProc=MaxStr_Proc)
double precision energia(0:max_ene)
#ifdef MPI
integer tmax_t,upindE_p
- double precision fi_p(MaxR,MaxT_h,Max_Parm)
+ double precision fi_p(MaxR,MaxT_h,Max_Parm),
+ & fimax_p(MaxR,MaxT_h,Max_Parm)
double precision sumW_p(0:nGridT,Max_Parm),
& sumE_p(0:nGridT,Max_Parm),sumEsq_p(0:nGridT,Max_Parm),
& sumQ_p(MaxQ1,0:nGridT,Max_Parm),
& hrmsrgy_p(0:MaxBinRgy,0:MaxBinRms,maxT_h)
double precision rgymin_t,rmsmin_t,rgymax_t,rmsmax_t
double precision potEmin_t,entmin_p,entmax_p
+ double precision ePMF,ePMF_q
+ double precision weimax_(0:ngridT)
integer histent_p(0:2000)
logical lprint /.true./
#endif
& sumEprim(0:NGridT,Max_Parm),sumEbis(0:NGridT,Max_Parm),betaT,
& weight,econstr
double precision fi(MaxR,maxT_h,Max_Parm),
+ & fimax(MaxR,maxT_h,Max_Parm),
& dd,dd1,dd2,hh,dmin,denom,finorm,avefi,pom,
& hfin(0:MaxHdim,maxT_h),histE(0:maxindE),
& hrmsrgy(0:MaxBinRgy,0:MaxBinRms,maxT_h),
& potEmin,ent,
- & hfin_ent(0:MaxHdim),vmax,aux
+ & hfin_ent(0:MaxHdim),vmax,aux,weimax(0:nGridT,Max_Parm)
double precision fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl,
& eprim,ebis,temper,kfac/2.4d0/,T0/300.0d0/,startGridT/200.0d0/,
& eplus,eminus,logfac,tanhT,tt
double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,
& escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
- & eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor
+ & eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,
+ & eliptran
integer ind_point(maxpoint),upindE,indE
character*16 plik
do t=0,MaxN
htot(t)=0
enddo
+C#define DEBUG
#ifdef MPI
do i=1,scount(me1)
#else
do i=1,ntot(islice)
#endif
+c write (iout,*) "i",i," potE",(potE(i,j),j=1,nParmset)
do j=1,nParmSet
if (potE(i,j).le.potEmin) potEmin=potE(i,j)
enddo
endif
c write (iout,*) "i",i," j",j," q",q(j,i)," ind_point",
c & ind_point(i)
- call flush(iout)
+c call flush(iout)
if (ind_point(i).lt.0 .or. ind_point(i).gt.MaxHdim) then
write (iout,*) "Error - index exceeds range for point",i,
& " q=",q(j,i)," ind",ind_point(i)
call flush(iout)
#endif
enddo ! i
- call flush(iout)
+ write (iout,*) "potEmin before reduce",potEmin
nbin=nbin1**nQ-1
write (iout,'(a)') "Numbers of counts in Q bins"
do t=0,tmax
& MPI_MIN,WHAM_COMM,IERROR)
call MPI_AllReduce(rgymax,rgymax_t,1,MPI_DOUBLE_PRECISION,
& MPI_MAX,WHAM_COMM,IERROR)
- potEmin=potEmin_t/2
+c potEmin=potEmin_t/2
+ potEmin=potEmin_t
rgymin=rgymin_t
rgymax=rgymax_t
rmsmin=rmsmin_t
do iparm=1,nParmSet
#ifdef DEBUG
write (iout,'(2i5,21f8.2)') i,iparm,
- & (enetb(k,i,iparm),k=1,21)
+ & (enetb(k,i,iparm),k=1,22)
#endif
call restore_parm(iparm)
#ifdef DEBUG
estr=enetb(18,i,iparm)
esccor=enetb(19,i,iparm)
edihcnstr=enetb(20,i,iparm)
+ eliptran=enetb(22,i,iparm)
+
#ifdef DEBUG
write (iout,'(3i5,6f5.2,14f12.3)') i,ib,iparm,(ft(l),l=1,6),
& evdw+evdw_t,evdw2,ees,evdw1,ecorr,eel_loc,estr,ebe,escloc,
#endif
#ifdef SPLITELE
+ if (shield_mode.gt.0) then
+ etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
+ & +ft(1)*welec*ees
+ & +ft(1)*wvdwpp*evdw1
+ & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
+! & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
+ & +ft(2)*wturn3*eello_turn3
+ & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
+ & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
+ & +wbond*estr+wliptran*eliptran
+ else
etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
& +wvdwpp*evdw1
& +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
- & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+c & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
& +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
& +ft(2)*wturn3*eello_turn3
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
& +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
- & +wbond*estr
+ & +wbond*estr+wliptran*eliptran
+ endif
#else
+ if (shield_mode.gt.0) then
+ etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
+ & +ft(1)*welec*(ees+evdw1)
+ & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
+c & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
+ & +ft(2)*wturn3*eello_turn3
+ & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
+ & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
+ & +wbond*estr+wliptran*eliptran
+ else
etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
& +ft(1)*welec*(ees+evdw1)
& +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
- & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+c & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
& +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
& +ft(2)*wturn3*eello_turn3
& +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
& +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
- & +wbond*estr
+ & +wbond*estr+wliptran*eliptran
+ endif
+
#endif
#ifdef DEBUG
write (iout,*) i,iparm,1.0d0/(beta_h(ib,iparm)*1.987D-3),
& etot,potEmin
#endif
+#define DEBUG
#ifdef DEBUG
if (iparm.eq.1 .and. ib.eq.1) then
write (iout,*)"Conformation",i
call enerprint(energia(0),fT)
endif
#endif
+#undef DEBUG
do kk=1,nR(ib,iparm)
Econstr=0.0d0
do j=1,nQ
Econstr=Econstr+Kh(j,kk,ib,iparm)
& *(dd-q0(j,kk,ib,iparm))**2
enddo
+c Adaptive potential contribution
+ if (adaptive) then
+ call PMF_energy(q(1,i),ib,kk,iparm,ePMF,ePMF_q)
+ Econstr=Econstr+ePMF
+ endif
v(i,kk,ib,iparm)=
& -beta_h(ib,iparm)*(etot-potEmin+Econstr)
#ifdef DEBUG
write (iout,*) t,"vmax",vmax," denom",denom,"entfac",entfac(t)
#endif
enddo
+
+ do iparm=1,nParmSet
+ do iib=1,nT_h(iparm)
+ do ii=1,nR(iib,iparm)
+#ifdef MPI
+ fimax_p(ii,iib,iparm)=v(1,ii,iib,iparm)+entfac(1)
+ do t=2,scount(me)
+ if(v(t,ii,iib,iparm)+entfac(t).gt.fimax_p(ii,iib,iparm))
+ & fimax_p(ii,iib,iparm)=v(t,ii,iib,iparm)+entfac(t)
+ enddo
+#else
+ fimax(ii,iib,iparm)=v(1,ii,iib,iparm)+entfac(1)
+ do t=2,ntot(islice)
+ if(v(t,ii,iib,iparm)+entfac(t).gt.fimax(ii,iib,iparm))
+ & fimax(ii,iib,iparm)=v(t,ii,iib,iparm)+entfac(t)
+ enddo
+#endif
+ enddo ! ii
+ enddo ! iib
+ enddo ! iparm
+#ifdef MPI
+ call MPI_AllReduce(fimax_p(1,1,1),fimax(1,1,1),
+ & maxR*MaxT_h*nParmSet,MPI_DOUBLE_PRECISION,
+ & MPI_MAX,WHAM_COMM,IERROR)
+#endif
do iparm=1,nParmSet
do iib=1,nT_h(iparm)
do ii=1,nR(iib,iparm)
fi_p(ii,iib,iparm)=0.0d0
do t=1,scount(me)
fi_p(ii,iib,iparm)=fi_p(ii,iib,iparm)
- & +dexp(v(t,ii,iib,iparm)+entfac(t))
+ & +dexp(v(t,ii,iib,iparm)+entfac(t)-fimax(ii,iib,iparm))
#ifdef DEBUG
write (iout,'(4i5,3e15.5)') t,ii,iib,iparm,
& v(t,ii,iib,iparm),entfac(t),fi_p(ii,iib,iparm)
fi(ii,iib,iparm)=0.0d0
do t=1,ntot(islice)
fi(ii,iib,iparm)=fi(ii,iib,iparm)
- & +dexp(v(t,ii,iib,iparm)+entfac(t))
+ & +dexp(v(t,ii,iib,iparm)+entfac(t)-fimax(ii,iib,iparm))
enddo
#endif
enddo ! ii
enddo
enddo
#endif
- write (iout,*) "REDUCE size",maxR,MaxT_h,nParmSet,
- & maxR*MaxT_h*nParmSet
- write (iout,*) "MPI_COMM_WORLD",MPI_COMM_WORLD,
- & " WHAM_COMM",WHAM_COMM
+c write (iout,*) "REDUCE size",maxR,MaxT_h,nParmSet,
+c & maxR*MaxT_h*nParmSet
+c write (iout,*) "MPI_COMM_WORLD",MPI_COMM_WORLD,
+c & " WHAM_COMM",WHAM_COMM
call MPI_Reduce(fi_p(1,1,1),fi(1,1,1),maxR*MaxT_h*nParmSet,
& MPI_DOUBLE_PRECISION,
& MPI_SUM,Master,WHAM_COMM,IERROR)
do iparm=1,nParmSet
do ib=1,nT_h(iparm)
do i=1,nR(ib,iparm)
- fi(i,ib,iparm)=-dlog(fi(i,ib,iparm))
+ fi(i,ib,iparm)=-dlog(fi(i,ib,iparm))-fimax(i,ib,iparm)
avefi=avefi+fi(i,ib,iparm)
enddo
enddo
write(iout,'(2f15.4)') entmin+i,histent(i)
enddo
#endif
-
-#ifdef MPI
-c write (iout,*) "me1",me1," scount",scount(me1)
-
do iparm=1,nParmSet
+ call restore_parm(iparm)
+c
+C Histograms
+c
#ifdef MPI
do ib=1,nT_h(iparm)
do t=0,tmax
enddo
enddo
enddo
-
+#ifdef MPI
do t=1,scount(me1)
#else
do t=1,ntot(islice)
#endif
- ind = ind_point(t)
+ indE = aint(potE(t,iparm)-aint(potEmin))
+ if (indE.ge.0 .and. indE.le.maxinde) then
+ if (indE.gt.upindE_p) upindE_p=indE
+ histE_p(indE)=histE_p(indE)+dexp(-entfac(t))
+ endif
#ifdef MPI
- hfin_ent_p(ind)=hfin_ent_p(ind)+dexp(entfac(t))
+ do ib=1,nT_h(iparm)
+ expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
+ hfin_p(ind,ib)=hfin_p(ind,ib)+
+ & dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
+ if (rmsrgymap) then
+ indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
+ indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
+ hrmsrgy_p(indrgy,indrms,ib)=
+ & hrmsrgy_p(indrgy,indrms,ib)+expfac
+ endif
+ enddo
#else
- hfin_ent(ind)=hfin_ent(ind)+dexp(entfac(t))
+ do ib=1,nT_h(iparm)
+ expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
+ hfin(ind,ib)=hfin(ind,ib)+
+ & dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
+ if (rmsrgymap) then
+ indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
+ indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
+ hrmsrgy(indrgy,indrms,ib)=
+ & hrmsrgy(indrgy,indrms,ib)+expfac
+ endif
+ enddo
#endif
-c write (iout,'(2i5,20f8.2)') t,t,(enetb(k,t,iparm),k=1,18)
- call restore_parm(iparm)
- evdw=enetb(21,t,iparm)
- evdw_t=enetb(1,t,iparm)
-#ifdef SCP14
- evdw2_14=enetb(17,t,iparm)
- evdw2=enetb(2,t,iparm)+evdw2_14
+ enddo ! t
+c
+c Thermo and ensemble averages
+c
+ do k=0,nGridT
+ betaT=startGridT+k*delta_T
+ call temp_scalfac(betaT,ft,ftprim,ftbis,*10)
+c write (iout,*) "ftprim",ftprim
+c write (iout,*) "ftbis",ftbis
+ betaT=1.0d0/(1.987D-3*betaT)
+c 7/10/18 AL Determine the max Botzmann weights for each temerature
+ call sum_ene(1,iparm,ft,etot)
+ weimax(k,iparm)=-betaT*(etot-potEmin)+entfac(1)
+c write (iout,*) "k",k," t",1," weight",weimax(k,iparm)
+#ifdef MPI
+ do t=2,scount(me1)
#else
- evdw2=enetb(2,t,iparm)
- evdw2_14=0.0d0
+ do t=2,ntot(islice)
#endif
-#ifdef SPLITELE
- ees=enetb(3,t,iparm)
- evdw1=enetb(16,t,iparm)
-#else
- ees=enetb(3,t,iparm)
- evdw1=0.0d0
+ call sum_ene(t,iparm,ft,etot)
+ weight=-betaT*(etot-potEmin)+entfac(t)
+c write (iout,*) "k",k," t",t," weight",weight
+ if (weight.gt.weimax(k,iparm)) weimax(k,iparm)=weight
+ enddo
+#ifdef MPI
+ enddo
+#ifdef DEBUG
+ write (iout,*) "weimax before REDUCE"
+ write (iout,*) (weimax(k,iparm),k=0,ngridt)
#endif
- ecorr=enetb(4,t,iparm)
- ecorr5=enetb(5,t,iparm)
- ecorr6=enetb(6,t,iparm)
- eel_loc=enetb(7,t,iparm)
- eello_turn3=enetb(8,t,iparm)
- eello_turn4=enetb(9,t,iparm)
- eturn6=enetb(10,t,iparm)
- ebe=enetb(11,t,iparm)
- escloc=enetb(12,t,iparm)
- etors=enetb(13,t,iparm)
- etors_d=enetb(14,t,iparm)
- ehpb=enetb(15,t,iparm)
- estr=enetb(18,t,iparm)
- esccor=enetb(19,t,iparm)
- edihcnstr=enetb(20,t,iparm)
- edihcnstr=0.0d0
- do k=0,nGridT
- betaT=startGridT+k*delta_T
- temper=betaT
-c fT=T0/betaT
-c ft=2*T0/(T0+betaT)
- if (rescale_mode.eq.1) then
- quot=betaT/T0
- quotl=1.0d0
- kfacl=1.0d0
- do l=1,5
- quotl1=quotl
- quotl=quotl*quot
- kfacl=kfacl*kfac
- denom=kfacl-1.0d0+quotl
- fT(l)=kfacl/denom
- ftprim(l)=-l*ft(l)*quotl1/(T0*denom)
- ftbis(l)=l*kfacl*quotl1*
- & (2*l*quotl-(l-1)*denom)/(quot*t0*t0*denom**3)
- enddo
-#if defined(FUNCTH)
- ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
- & 320.0d0
- ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
- ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
- & /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
-#elif defined(FUNCT)
- fT(6)=betaT/T0
- ftprim(6)=1.0d0/T0
- ftbis(6)=0.0d0
-#else
- fT(6)=1.0d0
- ftprim(6)=0.0d0
- ftbis(6)=0.0d0
+ do k=0,nGridT
+ weimax_(k)=weimax(k,iparm)
+ enddo
+ call MPI_Allreduce(weimax_(0),weimax(0,iparm),nGridT+1,
+ & MPI_DOUBLE_PRECISION,MPI_MAX,WHAM_COMM,IERROR)
+#ifdef DEBUG
+ write (iout,*) "weimax"
+ write (iout,*) (weimax(k,iparm),k=0,ngridt)
#endif
- else if (rescale_mode.eq.2) then
- quot=betaT/T0
- quotl=1.0d0
- do l=1,5
- quotl1=quotl
- quotl=quotl*quot
- eplus=dexp(quotl)
- eminus=dexp(-quotl)
- logfac=1.0d0/dlog(eplus+eminus)
- tanhT=(eplus-eminus)/(eplus+eminus)
- fT(l)=1.12692801104297249644d0*logfac
- ftprim(l)=-l*quotl1*ft(l)*tanhT*logfac/T0
- ftbis(l)=(l-1)*ftprim(l)/(quot*T0)-
- & 2*l*quotl1/T0*logfac*
- & (2*l*quotl1*ft(l)/(T0*(eplus+eminus)**2)
- & +ftprim(l)*tanhT)
- enddo
-#if defined(FUNCTH)
- ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
- & 320.0d0
- ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
- ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
- & /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
-#elif defined(FUNCT)
- fT(6)=betaT/T0
- ftprim(6)=1.0d0/T0
- ftbis(6)=0.0d0
+ do k=0,nGridT
+ temper=startGridT+k*delta_T
+ betaT=1.0d0/(1.987D-3*temper)
+ call temp_scalfac(temper,ft,ftprim,ftbis,*10)
+ do t=1,scount(me1)
#else
- fT(6)=1.0d0
- ftprim(6)=0.0d0
- ftbis(6)=0.0d0
+ do t=1,ntot(islice)
#endif
- else if (rescale_mode.eq.0) then
- do l=1,5
- fT(l)=1.0d0
- ftprim(l)=0.0d0
- enddo
- else
- write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
- & rescale_mode
- call flush(iout)
- return1
- endif
-c write (iout,*) "ftprim",ftprim
-c write (iout,*) "ftbis",ftbis
- betaT=1.0d0/(1.987D-3*betaT)
-#ifdef SPLITELE
- etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
- & +wvdwpp*evdw1
- & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
- & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
- & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
- & +ft(2)*wturn3*eello_turn3
- & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
- & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
- & +wbond*estr
- eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees
- & +ftprim(1)*wtor*etors+
- & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
- & ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
- & ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
- & ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
- & ftprim(1)*wsccor*esccor
- ebis=ftbis(1)*welec*ees+ftbis(1)*wtor*etors+
- & ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
- & ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
- & ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
- & ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
- & ftbis(1)*wsccor*esccor
+ ind = ind_point(t)
+#ifdef MPI
+ hfin_ent_p(ind)=hfin_ent_p(ind)+dexp(entfac(t))
#else
- etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
- & +ft(1)*welec*(ees+evdw1)
- & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
- & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
- & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
- & +ft(2)*wturn3*eello_turn3
- & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
- & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
- & +wbond*estr
- eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1)
- & +ftprim(1)*wtor*etors+
- & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
- & ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
- & ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
- & ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
- & ftprim(1)*wsccor*esccor
- ebis=ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+
- & ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
- & ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
- & ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
- & ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
- & ftprim(1)*wsccor*esccor
+ hfin_ent(ind)=hfin_ent(ind)+dexp(entfac(t))
#endif
- weight=dexp(-betaT*(etot-potEmin)+entfac(t))
+c write (iout,'(2i5,20f8.2)') t,t,(enetb(k,t,iparm),k=1,18)
+c call restore_parm(iparm)
+ call sum_ene_deriv(t,iparm,ft,ftprim,ftbis,etot,eprim,ebis)
+ weight=dexp(-betaT*(etot-potEmin)+entfac(t)-weimax(k,iparm))
#ifdef DEBUG
write (iout,*) "iparm",iparm," t",t," betaT",betaT,
- & " etot",etot," entfac",entfac(t),
+ & " etot",etot," entfac",entfac(t)," boltz",
+ & -betaT*(etot-potEmin)+entfac(t)," weimax",weimax(k,iparm),
& " weight",weight," ebis",ebis
#endif
etot=etot-temper*eprim
& +etot*q(j,t)*weight
enddo
#endif
- enddo
- indE = aint(potE(t,iparm)-aint(potEmin))
- if (indE.ge.0 .and. indE.le.maxinde) then
- if (indE.gt.upindE_p) upindE_p=indE
- histE_p(indE)=histE_p(indE)+dexp(-entfac(t))
- endif
-#ifdef MPI
- do ib=1,nT_h(iparm)
- expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
- hfin_p(ind,ib)=hfin_p(ind,ib)+
- & dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
- if (rmsrgymap) then
- indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
- indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
- hrmsrgy_p(indrgy,indrms,ib)=
- & hrmsrgy_p(indrgy,indrms,ib)+expfac
- endif
- enddo
-#else
- do ib=1,nT_h(iparm)
- expfac=dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
- hfin(ind,ib)=hfin(ind,ib)+
- & dexp(-beta_h(ib,iparm)*(etot-potEmin)+entfac(t))
- if (rmsrgymap) then
- indrgy=dint((q(nQ+2,t)-rgymin)/deltrgy)
- indrms=dint((q(nQ+1,t)-rmsmin)/deltrms)
- hrmsrgy(indrgy,indrms,ib)=
- & hrmsrgy(indrgy,indrms,ib)+expfac
- endif
- enddo
-#endif
- enddo ! t
+ enddo ! t
+ enddo ! k
do ib=1,nT_h(iparm)
if (histout) call MPI_Reduce(hfin_p(0,ib),hfin(0,ib),nbin,
& MPI_DOUBLE_PRECISION,MPI_SUM,Master,WHAM_COMM,IERROR)
sumEQ(j,i,iparm)=sumEQ(j,i,iparm)/sumW(i,iparm)
& -sumQ(j,i,iparm)*sumE(i,iparm)
enddo
- sumW(i,iparm)=-dlog(sumW(i,iparm))*(1.987D-3*
+ sumW(i,iparm)=(-dlog(sumW(i,iparm))-weimax(i,iparm))*(1.987D-3*
& (startGridT+i*delta_T))+potEmin
write (iout,'(f7.1,2f15.5,$)') startGridT+i*delta_T,
& sumW(i,iparm),sumE(i,iparm)
#endif
return
-
+C#undef DEBUG
+ 10 return1
+ end
+c------------------------------------------------------------------------
+ subroutine temp_scalfac(betaT,ft,ftprim,ftbis,*)
+ implicit none
+ include "DIMENSIONS"
+ include "DIMENSIONS.FREE"
+ include "COMMON.CONTROL"
+ include "COMMON.FREE"
+ include "COMMON.IOUNITS"
+ double precision fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl,
+ & kfac/2.4d0/,T0/300.0d0/,startGridT/200.0d0/,
+ & logfac,tanhT,betaT,denom,eplus,eminus
+ integer l
+ if (rescale_mode.eq.1) then
+ quot=betaT/T0
+ quotl=1.0d0
+ kfacl=1.0d0
+ do l=1,5
+ quotl1=quotl
+ quotl=quotl*quot
+ kfacl=kfacl*kfac
+ denom=kfacl-1.0d0+quotl
+ fT(l)=kfacl/denom
+ ftprim(l)=-l*ft(l)*quotl1/(T0*denom)
+ ftbis(l)=l*kfacl*quotl1*
+ & (2*l*quotl-(l-1)*denom)/(quot*t0*t0*denom**3)
+ enddo
+#if defined(FUNCTH)
+ ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
+ & 320.0d0
+ ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
+ ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
+ & /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
+#elif defined(FUNCT)
+ fT(6)=betaT/T0
+ ftprim(6)=1.0d0/T0
+ ftbis(6)=0.0d0
+#else
+ fT(6)=1.0d0
+ ftprim(6)=0.0d0
+ ftbis(6)=0.0d0
+#endif
+ else if (rescale_mode.eq.2) then
+ quot=betaT/T0
+ quotl=1.0d0
+ do l=1,5
+ quotl1=quotl
+ quotl=quotl*quot
+ eplus=dexp(quotl)
+ eminus=dexp(-quotl)
+ logfac=1.0d0/dlog(eplus+eminus)
+ tanhT=(eplus-eminus)/(eplus+eminus)
+ fT(l)=1.12692801104297249644d0*logfac
+ ftprim(l)=-l*quotl1*ft(l)*tanhT*logfac/T0
+ ftbis(l)=(l-1)*ftprim(l)/(quot*T0)-
+ & 2*l*quotl1/T0*logfac*
+ & (2*l*quotl1*ft(l)/(T0*(eplus+eminus)**2)
+ & +ftprim(l)*tanhT)
+ enddo
+#if defined(FUNCTH)
+ ft(6)=(320.0d0+80.0d0*dtanh((betaT-320.0d0)/80.0d0))/
+ & 320.0d0
+ ftprim(6)=1.0d0/(320.0d0*dcosh((betaT-320.0d0)/80.0d0)**2)
+ ftbis(6)=-2.0d0*dtanh((betaT-320.0d0)/80.0d0)
+ & /(320.0d0*80.0d0*dcosh((betaT-320.0d0)/80.0d0)**3)
+#elif defined(FUNCT)
+ fT(6)=betaT/T0
+ ftprim(6)=1.0d0/T0
+ ftbis(6)=0.0d0
+#else
+ fT(6)=1.0d0
+ ftprim(6)=0.0d0
+ ftbis(6)=0.0d0
+#endif
+ else if (rescale_mode.eq.0) then
+ do l=1,5
+ fT(l)=1.0d0
+ ftprim(l)=0.0d0
+ enddo
+ else
+ write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE",
+ & rescale_mode
+ call flush(iout)
+ return1
+ endif
+ return
+ end
+c--------------------------------------------------------------------
+ subroutine sum_ene(t,iparm,ft,etot)
+ implicit none
+ include 'DIMENSIONS'
+ include 'DIMENSIONS.ZSCOPT'
+ include 'DIMENSIONS.FREE'
+ include 'COMMON.CONTROL'
+ include 'COMMON.FFIELD'
+ include "COMMON.SBRIDGE"
+ include "COMMON.ENERGIES"
+ include "COMMON.IOUNITS"
+ integer t,iparm
+ double precision fT(6)
+ double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,
+ & escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
+ & eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,
+ & eliptran
+ evdw=enetb(21,t,iparm)
+ evdw_t=enetb(1,t,iparm)
+#ifdef SCP14
+ evdw2_14=enetb(17,t,iparm)
+ evdw2=enetb(2,t,iparm)+evdw2_14
+#else
+ evdw2=enetb(2,t,iparm)
+ evdw2_14=0.0d0
+#endif
+#ifdef SPLITELE
+ ees=enetb(3,t,iparm)
+ evdw1=enetb(16,t,iparm)
+#else
+ ees=enetb(3,t,iparm)
+ evdw1=0.0d0
+#endif
+ ecorr=enetb(4,t,iparm)
+ ecorr5=enetb(5,t,iparm)
+ ecorr6=enetb(6,t,iparm)
+ eel_loc=enetb(7,t,iparm)
+ eello_turn3=enetb(8,t,iparm)
+ eello_turn4=enetb(9,t,iparm)
+ eturn6=enetb(10,t,iparm)
+ ebe=enetb(11,t,iparm)
+ escloc=enetb(12,t,iparm)
+ etors=enetb(13,t,iparm)
+ etors_d=enetb(14,t,iparm)
+ ehpb=enetb(15,t,iparm)
+ estr=enetb(18,t,iparm)
+ esccor=enetb(19,t,iparm)
+ edihcnstr=enetb(20,t,iparm)
+ edihcnstr=0.0d0
+#ifdef SPLITELE
+ if (shield_mode.gt.0) then
+ etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
+ & +ft(1)*welec*ees
+ & +ft(1)*wvdwpp*evdw1
+ & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
+c & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
+ & +ft(2)*wturn3*eello_turn3
+ & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
+ & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
+ & +wbond*estr+wliptran*eliptran
+ else
+ etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
+ & +wvdwpp*evdw1
+ & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
+c & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
+ & +ft(2)*wturn3*eello_turn3
+ & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
+ & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
+ & +wbond*estr+wliptran*eliptran
+ endif
+#else
+ if (shield_mode.gt.0) then
+ etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
+ & +ft(1)*welec*(ees+evdw1)
+ & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
+c & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
+ & +ft(2)*wturn3*eello_turn3
+ & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
+ & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
+ & +wbond*estr+wliptran*eliptran
+ else
+ etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
+ & +ft(1)*welec*(ees+evdw1)
+ & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
+c & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
+ & +ft(2)*wturn3*eello_turn3
+ & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
+ & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
+ & +wbond*estr+wliptran*eliptran
+ endif
+#endif
+ return
+ end
+c--------------------------------------------------------------------
+ subroutine sum_ene_deriv(t,iparm,ft,ftprim,ftbis,etot,eprim,ebis)
+ implicit none
+ include 'DIMENSIONS'
+ include 'DIMENSIONS.ZSCOPT'
+ include 'DIMENSIONS.FREE'
+ include 'COMMON.CONTROL'
+ include 'COMMON.FFIELD'
+ include "COMMON.SBRIDGE"
+ include 'COMMON.ENERGIES'
+ include "COMMON.IOUNITS"
+ integer t,iparm
+ double precision fT(6),fTprim(6),fTbis(6),
+ & eprim,ebis,temper
+ double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,
+ & escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,
+ & eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,
+ & eliptran
+ evdw=enetb(21,t,iparm)
+ evdw_t=enetb(1,t,iparm)
+#ifdef SCP14
+ evdw2_14=enetb(17,t,iparm)
+ evdw2=enetb(2,t,iparm)+evdw2_14
+#else
+ evdw2=enetb(2,t,iparm)
+ evdw2_14=0.0d0
+#endif
+#ifdef SPLITELE
+ ees=enetb(3,t,iparm)
+ evdw1=enetb(16,t,iparm)
+#else
+ ees=enetb(3,t,iparm)
+ evdw1=0.0d0
+#endif
+ ecorr=enetb(4,t,iparm)
+ ecorr5=enetb(5,t,iparm)
+ ecorr6=enetb(6,t,iparm)
+ eel_loc=enetb(7,t,iparm)
+ eello_turn3=enetb(8,t,iparm)
+ eello_turn4=enetb(9,t,iparm)
+ eturn6=enetb(10,t,iparm)
+ ebe=enetb(11,t,iparm)
+ escloc=enetb(12,t,iparm)
+ etors=enetb(13,t,iparm)
+ etors_d=enetb(14,t,iparm)
+ ehpb=enetb(15,t,iparm)
+ estr=enetb(18,t,iparm)
+ esccor=enetb(19,t,iparm)
+ edihcnstr=enetb(20,t,iparm)
+ edihcnstr=0.0d0
+#ifdef SPLITELE
+ if (shield_mode.gt.0) then
+ etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
+ & +ft(1)*welec*ees
+ & +ft(1)*wvdwpp*evdw1
+ & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
+c & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
+ & +ft(2)*wturn3*eello_turn3
+ & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
+ & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
+ & +wbond*estr+wliptran*eliptran
+ eprim=ftprim(1)*(ft(6)*evdw_t+evdw)
+C & +ftprim(6)*evdw_t
+ & +ftprim(1)*wscp*evdw2
+ & +ftprim(1)*welec*ees
+ & +ftprim(1)*wvdwpp*evdw1
+ & +ftprim(1)*wtor*etors+
+ & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
+ & ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
+ & ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
+ & ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
+ & ftprim(1)*wsccor*esccor
+ ebis=ftbis(1)*wsc*(evdw+ft(6)*evdw_t)
+ & +ftbis(1)*wscp*evdw2+
+ & ftbis(1)*welec*ees
+ & +ftbis(1)*wvdwpp*evdw
+ & +ftbis(1)*wtor*etors+
+ & ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
+ & ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
+ & ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
+ & ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
+ & ftbis(1)*wsccor*esccor
+ else
+ etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees
+ & +wvdwpp*evdw1
+ & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
+c & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
+ & +ft(2)*wturn3*eello_turn3
+ & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc
+ & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
+ & +wbond*estr+wliptran*eliptran
+ eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees
+ & +ftprim(1)*wtor*etors+
+ & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
+ & ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
+ & ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
+ & ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
+ & ftprim(1)*wsccor*esccor
+ ebis=ftbis(1)*welec*ees+ftbis(1)*wtor*etors+
+ & ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
+ & ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
+ & ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
+ & ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
+ & ftbis(1)*wsccor*esccor
+ endif
+#else
+ if (shield_mode.gt.0) then
+ etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2
+ & +ft(1)*welec*(ees+evdw1)
+ & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
+c & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
+ & +ft(2)*wturn3*eello_turn3
+ & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
+ & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
+ & +wbond*estr+wliptran*eliptran
+ eprim=ftprim(1)*(evdw+ft(6)*evdw_t)
+ & +ftprim(1)*welec*(ees+evdw1)
+ & +ftprim(1)*wtor*etors+
+ & ftprim(1)*wscp*evdw2+
+ & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
+ & ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
+ & ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
+ & ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
+ & ftprim(1)*wsccor*esccor
+ ebis= ftbis(1)*(evdw+ft(6)*evdw_t)
+ & +ftbis(1)*wscp*evdw2
+ & +ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+
+ & ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
+ & ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
+ & ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
+ & ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
+ & ftprim(1)*wsccor*esccor
+ else
+ etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2
+ & +ft(1)*welec*(ees+evdw1)
+ & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc
+c & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5
+ & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4
+ & +ft(2)*wturn3*eello_turn3
+ & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr
+ & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor
+ & +wbond*estr+wliptran*eliptran
+ eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1)
+ & +ftprim(1)*wtor*etors+
+ & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+
+ & ftprim(5)*wcorr6*ecorr6+ftprim(3)*wturn4*eello_turn4+
+ & ftprim(2)*wturn3*eello_turn3+ftprim(5)*wturn6*eturn6+
+ & ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+
+ & ftprim(1)*wsccor*esccor
+ ebis=ftbis(1)*welec*(ees+evdw1)+ftbis(1)*wtor*etors+
+ & ftbis(3)*wcorr*ecorr+ftbis(4)*wcorr5*ecorr5+
+ & ftbis(5)*wcorr6*ecorr6+ftbis(3)*wturn4*eello_turn4+
+ & ftbis(2)*wturn3*eello_turn3+ftbis(5)*wturn6*eturn6+
+ & ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
+ & ftprim(1)*wsccor*esccor
+ endif
+#endif
+ return
end