X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fsrc-M%2Fwham_calc1.F;h=35a9278dd42558fe96bbcf1f107d37934055f145;hb=34d3ad3987785642be58fb2f26557d3314215577;hp=15d6716a051c997b7d10060e36198cb27c42554a;hpb=f690e8b70bab14132839afebf080d4a28363b226;p=unres.git diff --git a/source/wham/src-M/wham_calc1.F b/source/wham/src-M/wham_calc1.F index 15d6716..35a9278 100644 --- a/source/wham/src-M/wham_calc1.F +++ b/source/wham/src-M/wham_calc1.F @@ -18,9 +18,8 @@ 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 @@ -36,6 +35,7 @@ c parameter (MaxHdim=200000) include "COMMON.SBRIDGE" include "COMMON.PROT" include "COMMON.ENEPS" + include "COMMON.SHIELD" integer MaxPoint,MaxPointProc parameter (MaxPoint=MaxStr, & MaxPointProc=MaxStr_Proc) @@ -50,7 +50,8 @@ c parameter (MaxHdim=200000) 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), @@ -63,6 +64,8 @@ c parameter (MaxHdim=200000) & 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 @@ -74,17 +77,19 @@ c parameter (MaxHdim=200000) & 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 @@ -112,11 +117,13 @@ c parameter (MaxHdim=200000) 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 @@ -134,7 +141,7 @@ c parameter (MaxHdim=200000) 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) @@ -154,8 +161,8 @@ c & 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 @@ -194,7 +201,8 @@ c & ind_point(i) & 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 @@ -227,7 +235,7 @@ c write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet) 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 @@ -312,6 +320,8 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft 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, @@ -319,30 +329,61 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft #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 @@ -353,6 +394,7 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft call enerprint(energia(0),fT) endif #endif +#undef DEBUG do kk=1,nR(ib,iparm) Econstr=0.0d0 do j=1,nQ @@ -360,6 +402,11 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft 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 @@ -406,6 +453,31 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft 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) @@ -413,7 +485,7 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft 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) @@ -423,7 +495,7 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft 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 @@ -441,10 +513,10 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft 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) @@ -464,7 +536,7 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft 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 @@ -607,12 +679,12 @@ c ient=-dlog(entfac(t))-entmin 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 @@ -642,182 +714,102 @@ c write (iout,*) "me1",me1," scount",scount(me1) 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 @@ -844,38 +836,8 @@ c write (iout,*) "ftbis",ftbis & +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) @@ -1094,7 +1056,7 @@ c write (iout,*) "ftbis",ftbis 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) @@ -1190,5 +1152,359 @@ c write (iout,*) "ftbis",ftbis #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