!el real(kind=8) :: energia(0:max_ene)
#ifdef MPI
integer :: tmax_t,upindE_p
- real(kind=8) :: fi_p(MaxR,MaxT_h,nParmSet) !(MaxR,MaxT_h,Max_Parm)
+ real(kind=8) :: fi_p(MaxR,MaxT_h,nParmSet), &
+ fi_p_min(MaxR,MaxT_h,nParmSet) !(MaxR,MaxT_h,Max_Parm)
real(kind=8),dimension(0:nGridT,nParmSet) :: sumW_p,sumE_p,&
sumEbis_p,sumEsq_p !(0:nGridT,Max_Parm)
real(kind=8),dimension(MaxQ1,0:nGridT,nParmSet) :: sumQ_p,&
sumQsq_p,sumEQ_p,sumEprim_p !(MaxQ1,0:nGridT,Max_Parm)
real(kind=8) :: hfin_p(0:MaxHdim,maxT_h),&
hfin_ent_p(0:MaxHdim),histE_p(0:maxindE),sumH,&
- hrmsrgy_p(0:MaxBinRgy,0:MaxBinRms,maxT_h)
+ hrmsrgy_p(0:MaxBinRgy,0:MaxBinRms,maxT_h),weimax_(0:ngridT)
real(kind=8) :: rgymin_t,rmsmin_t,rgymax_t,rmsmax_t
- real(kind=8) :: potEmin_t!,entmin_p,entmax_p
+ real(kind=8) :: potEmin_t,potEmin_t_all(maxT_h,Max_Parm)!,entmin_p,entmax_p
! integer :: histent_p(0:2000)
logical :: lprint=.true.
#endif
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,fi_min(MaxR,maxT_h,nParmSet), &
+ potEmin_all(maxT_h,Max_Parm),potEmin_min,entfac_min,&
+ weimax(0:nGridT,Max_Parm)
real(kind=8) :: 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
real(kind=8) :: etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors,&
escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3,&
eello_turn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor, &
- ecationcation,ecation_prot
+ ecationcation,ecation_prot, evdwpp,eespp ,evdwpsb,eelpsb, &
+ evdwsb, eelsb, estr_nucl,ebe_nucl,esbloc,etors_nucl,etors_d_nucl,&
+ ecorr_nucl, ecorr3_nucl,epeppho, escpho, epepbase,escbase,ecation_nucl,&
+ elipbond,elipang,eliplj,elipelec,ecat_prottran,ecation_protang,eliptran,&
+ ehomology_constr
+
integer :: ind_point(maxpoint),upindE,indE
character(len=16) :: plik
dmin=0.0d0
tmax=0
potEmin=1.0d10
+ do i=1,nParmset
+ do j=1,nT_h(i)
+ potEmin_all(j,i)=1.0d15
+ enddo
+ enddo
+ entfac_min=1.0d10
rgymin=1.0d10
rmsmin=1.0d10
rgymax=0.0d0
call MPI_AllReduce(tmax,tmax_t,1,MPI_INTEGER,MPI_MAX,&
WHAM_COMM,IERROR)
tmax=tmax_t
- call MPI_AllReduce(potEmin,potEmin_t,1,MPI_DOUBLE_PRECISION,&
- MPI_MIN,WHAM_COMM,IERROR)
+! call MPI_AllReduce(potEmin,potEmin_t,1,MPI_DOUBLE_PRECISION,&
+! MPI_MIN,WHAM_COMM,IERROR) !????
call MPI_AllReduce(rmsmin,rmsmin_t,1,MPI_DOUBLE_PRECISION,&
MPI_MIN,WHAM_COMM,IERROR)
call MPI_AllReduce(rmsmax,rmsmax_t,1,MPI_DOUBLE_PRECISION,&
MPI_MIN,WHAM_COMM,IERROR)
call MPI_AllReduce(rgymax,rgymax_t,1,MPI_DOUBLE_PRECISION,&
MPI_MAX,WHAM_COMM,IERROR)
- potEmin=potEmin_t !/2 try now??
+! potEmin=potEmin_t !/2 try now??
rgymin=rgymin_t
rgymax=rgymax_t
rmsmin=rmsmin_t
#endif
! write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet)
do iparm=1,nParmSet
-!#ifdef DEBUG
+#ifdef DEBUG
write (iout,'(2i5,21f8.2)') i,iparm,&
(enetb(k,i,iparm),k=1,21)
-!#endif
+#endif
call restore_parm(iparm)
-!#ifdef DEBUG
+#ifdef DEBUG
write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc,&
wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc,&
wtor_d,wsccor,wbond,wcatcat
-!#endif
+#endif
do ib=1,nT_h(iparm)
!el old rascale weights
!
edihcnstr=enetb(19,i,iparm)
ecationcation=enetb(41,i,iparm)
ecation_prot=enetb(42,i,iparm)
+ evdwpp = enetb(26,i,iparm)
+ eespp = enetb(27,i,iparm)
+ evdwpsb = enetb(28,i,iparm)
+ eelpsb = enetb(29,i,iparm)
+ evdwsb = enetb(30,i,iparm)
+ eelsb = enetb(31,i,iparm)
+ estr_nucl = enetb(32,i,iparm)
+ ebe_nucl = enetb(33,i,iparm)
+ esbloc = enetb(34,i,iparm)
+ etors_nucl = enetb(35,i,iparm)
+ etors_d_nucl = enetb(36,i,iparm)
+ ecorr_nucl = enetb(37,i,iparm)
+ ecorr3_nucl = enetb(38,i,iparm)
+ epeppho= enetb(49,i,iparm)
+ escpho= enetb(48,i,iparm)
+ epepbase= enetb(47,i,iparm)
+ escbase= enetb(46,i,iparm)
+ ecation_nucl= enetb(50,i,iparm)
+ elipbond=enetb(52,i,iparm)
+ elipang=enetb(53,i,iparm)
+ eliplj=enetb(54,i,iparm)
+ elipelec=enetb(55,i,iparm)
+ ecat_prottran=enetb(56,i,iparm)
+ ecation_protang=enetb(57,i,iparm)
+ eliptran=enetb(22,i,iparm)
+ ehomology_constr=enetb(51,i,iparm)
#ifdef DEBUG
write (iout,'(3i5,6f5.2,15f12.3)') i,ib,iparm,(ft(l),l=1,6),&
evdw+evdw_t,evdw2,ees,evdw1,ecorr,eel_loc,estr,ebe,escloc,&
- etors,etors_d,eello_turn3,eello_turn4,esccor,ecationcation
+ etors,etors_d,eello_turn3,eello_turn4,esccor,ecationcation,&
+ ehomology_constr
#endif
!#ifdef SPLITELE
! +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 &
+! +ft(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc &
! +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
! +wbond*estr
!#else
! +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(5)*wturn6*eello_turn6+ft(2)*wel_loc*eel_loc+edihcnstr &
! +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
! +wbond*estr
!#endif
etot=wsc*evdw+wscp*evdw2+welec*ees &
+wvdwpp*evdw1 &
+wang*ebe+wtor*etors+wscloc*escloc &
- +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5 &
+ +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5 &
+wcorr6*ecorr6+wturn4*eello_turn4 &
+wturn3*eello_turn3 &
+wturn6*eello_turn6+wel_loc*eel_loc &
+edihcnstr+wtor_d*etors_d+wsccor*esccor &
- +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation
+ +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
+ +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
+ +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
+ +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
+ +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl&
+ +wscbase*escbase&
+ +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho+wcatnucl*ecation_nucl&
+ +elipbond+elipang+eliplj+elipelec+wcat_tran*ecat_prottran+ecation_protang&
+ +wliptran*eliptran
+
+
#else
etot=wsc*evdw+wscp*evdw2 &
+welec*(ees+evdw1) &
+wang*ebe+wtor*etors+wscloc*escloc &
- +wstrain*ehpb+nss*ebr+wcorr*ecorr+wcorr5*ecorr5 &
+ +wstrain*ehpb+wcorr*ecorr+wcorr5*ecorr5 &
+wcorr6*ecorr6+wturn4*eello_turn4 &
+wturn3*eello_turn3 &
+wturn6*eello_turn6+wel_loc*eel_loc+edihcnstr &
+wtor_d*etors_d+wsccor*esccor &
- +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation
+ +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
+ +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
+ +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
+ +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
+ +wtor_d_nucl*etors_d_nucl+wcorr_nucl*ecorr_nucl+wcorr3_nucl*ecorr3_nucl&
+ +wscbase*escbase&
+ +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho+wcatnucl*ecation_nucl&
+ +elipbond+elipang+eliplj+elipelec+wcat_tran*ecat_prottran+ecation_protang&
+ +wliptran*eliptran
+
+
+
+
#endif
#ifdef DEBUG
call enerprint(energia(0))
endif
#endif
+ write (iout,*) "homol_nset",homol_nset,nR(ib,iparm)
+ write (iout,*) "waga_homology", waga_homology
+#ifdef DEBUG
+ write (iout,*) "homol_nset",homol_nset,nR(ib,iparm)
+#endif
+ if (homol_nset.gt.1) then
+
+ do kk=1,nR(ib,iparm)
+ Econstr=waga_homology(kk)*ehomology_constr
+ v(i,kk,ib,iparm)= &
+ -beta_h(ib,iparm)*(etot+Econstr)
+#ifdef DEBUG
+ write (iout,'(4i5,4e15.5)') i,kk,ib,iparm, &
+ etot,Econstr,v(i,kk,ib,iparm)
+#endif
+ enddo ! kk
+
+ else
+
+ etot=etot+ehomology_constr
+
+
do kk=1,nR(ib,iparm)
Econstr=0.0d0
do j=1,nQ
etot,potEmin,etot-potEmin,v(i,kk,ib,iparm)
#endif
enddo ! kk
+ endif
enddo ! ib
enddo ! iparm
enddo ! i
enddo
enddo
entfac(t)=-dlog(denom)-vmax
+ if (entfac(t).lt.entfac_min) entfac_min=entfac(t)
#ifdef DEBUG
write (iout,*) t,"vmax",vmax," denom",denom,"entfac",entfac(t)
#endif
do iib=1,nT_h(iparm)
do ii=1,nR(iib,iparm)
#ifdef MPI
+ fi_p_min(ii,iib,iparm)=-1.0d5
+ do t=1,scount(me)
+! fi_p(ii,iib,iparm)=fi_p(ii,iib,iparm) &
+! +dexp(v(t,ii,iib,iparm)+entfac(t))
+ aux=v(t,ii,iib,iparm)+entfac(t)
+ if (aux.gt.fi_p_min(ii,iib,iparm)) fi_p_min(ii,iib,iparm)=aux
+
+!#define DEBUG
+#ifdef DEBUG
+ write (iout,'(4i5,4e15.5)') t,ii,iib,iparm,&
+ v(t,ii,iib,iparm),entfac(t),fi_p(ii,iib,iparm),fi_p_min(ii,iib,iparm)
+#endif
+!#undef DEBUG
+ enddo
+#else
+! 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))
+ aux=v(t,ii,iib,iparm)+entfac(t)
+ if (aux.gt.fi_min(ii,iib,iparm))
+ & fi_min(ii,iib,iparm)=aux
+
+ enddo
+#endif
+ enddo ! ii
+ enddo ! iib
+ enddo ! iparm
+
+#ifdef MPI
+#ifdef DEBUG
+ write (iout,*) "fi before MPI_Reduce me",me,' master',master
+ do iparm=1,nParmSet
+ do ib=1,nT_h(nparmset)
+ write (iout,*) "iparm",iparm," ib",ib
+ write (iout,*) "beta=",beta_h(ib,iparm)
+ write (iout,'(8e15.5)') (fi_p(i,ib,iparm),i=1,nR(ib,iparm))
+ write (iout,'(8e15.5)') (fi_p_min(i,ib,iparm),i=1,nR(ib,iparm))
+ enddo
+ enddo
+#endif
+ call MPI_AllReduce(fi_p_min,fi_min,MaxR*MaxT_h*nParmSet, &
+ MPI_DOUBLE_PRECISION,MPI_MAX,WHAM_COMM,IERROR)
+
+#ifdef DEBUG
+ write (iout,*) "fi_min after AllReduce"
+ do i=1,nParmSet
+ do j=1,nT_h(i)
+ write (iout,*) (i,j,k,fi_min(k,j,i),k=1,nR(j,i))
+ enddo
+ enddo
+#endif
+!#undef DEBUG
+#endif
+ do iparm=1,nParmSet
+ do iib=1,nT_h(iparm)
+ do ii=1,nR(iib,iparm)
+#ifdef MPI
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)-fi_min(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)
+ write (iout,'(4i5,4e15.5)') t,ii,iib,iparm, &
+ v(t,ii,iib,iparm),entfac(t),fi_min(ii,iib,iparm), &
+ fi_p(ii,iib,iparm)
#endif
enddo
#else
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)-fi_min(ii,iib,iparm))
enddo
#endif
enddo ! ii
enddo
enddo
#endif
+#ifdef DEBUG
+ write (iout,*) "REDUCE size",maxR,MaxT_h,nParmSet, &
+ maxR*MaxT_h*nParmSet
+ write (iout,*) "MPI_COMM_WORLD",MPI_COMM_WORLD, &
+ " WHAM_COMM",WHAM_COMM
+#endif
+
write (iout,*) "REDUCE size",maxR,MaxT_h,nParmSet,&
maxR*MaxT_h*nParmSet
write (iout,*) "MPI_COMM_WORLD",MPI_COMM_WORLD,&
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))-fi_min(i,ib,iparm)
avefi=avefi+fi(i,ib,iparm)
enddo
enddo
20 continue
! Now, put together the histograms from all simulations, in order to get the
! unbiased total histogram.
+!C Determine the minimum free energies
+#ifdef MPI
+ do i=1,scount(me1)
+#else
+ do i=1,ntot(islice)
+#endif
+!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,26)
+#endif
+ call restore_parm(iparm)
+#ifdef DEBUG
+ write (iout,*) wsc,wscp,welec,wvdwpp,wang,wtor,wscloc, &
+ wcorr,wcorr5,wcorr6,wturn4,wturn3,wturn6,wel_loc, &
+ wtor_d,wsccor,wbond
+#endif
+ do ib=1,nT_h(iparm)
+ if (rescale_modeW.eq.1) then
+ quot=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
+ quotl=1.0d0
+ kfacl=1.0d0
+ do l=1,5
+ quotl1=quotl
+ quotl=quotl*quot
+ kfacl=kfacl*kfac
+ fT(l)=kfacl/(kfacl-1.0d0+quotl)
+ enddo
+#if defined(FUNCTH)
+ tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
+ ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
+#elif defined(FUNCT)
+ ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
+#else
+ ft(6)=1.0d0
+#endif
+ else if (rescale_modeW.eq.2) then
+ quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3)
+ quotl=1.0d0
+ do l=1,5
+ quotl=quotl*quot
+ fT(l)=1.12692801104297249644d0/ &
+ dlog(dexp(quotl)+dexp(-quotl))
+ enddo
+#if defined(FUNCTH)
+ tt = 1.0d0/(beta_h(ib,iparm)*1.987D-3)
+ ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0
+#elif defined(FUNCT)
+ ft(6)=1.0d0/(beta_h(ib,iparm)*1.987D-3*T0)
+#else
+ ft(6)=1.0d0
+#endif
+!c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft
+ else if (rescale_modeW.eq.0) then
+ do l=1,6
+ fT(l)=1.0d0
+ enddo
+ else
+ write (iout,*) "Error in WHAM_CALC: wrong RESCALE_MODE", &
+ rescale_mode
+ call flush(iout)
+ return 1
+ endif
+ evdw=enetb(1,i,iparm)
+! evdw_t=enetb(21,i,iparm)
+ evdw_t=enetb(20,i,iparm)
+#ifdef SCP14
+! evdw2_14=enetb(17,i,iparm)
+ evdw2_14=enetb(18,i,iparm)
+ evdw2=enetb(2,i,iparm)+evdw2_14
+#else
+ evdw2=enetb(2,i,iparm)
+ evdw2_14=0.0d0
+#endif
+#ifdef SPLITELE
+ ees=enetb(3,i,iparm)
+ evdw1=enetb(16,i,iparm)
+#else
+ ees=enetb(3,i,iparm)
+ evdw1=0.0d0
+#endif
+ ecorr=enetb(4,i,iparm)
+ ecorr5=enetb(5,i,iparm)
+ ecorr6=enetb(6,i,iparm)
+ eel_loc=enetb(7,i,iparm)
+ eello_turn3=enetb(8,i,iparm)
+ eello_turn4=enetb(9,i,iparm)
+ eello_turn6=enetb(10,i,iparm)
+ ebe=enetb(11,i,iparm)
+ escloc=enetb(12,i,iparm)
+ etors=enetb(13,i,iparm)
+ etors_d=enetb(14,i,iparm)
+ ehpb=enetb(15,i,iparm)
+! estr=enetb(18,i,iparm)
+ estr=enetb(17,i,iparm)
+! esccor=enetb(19,i,iparm)
+ esccor=enetb(21,i,iparm)
+! edihcnstr=enetb(20,i,iparm)
+ edihcnstr=enetb(19,i,iparm)
+! ehomology_constr=enetb(22,i,iparm)
+! esaxs=enetb(26,i,iparm)
+ ecationcation=enetb(41,i,iparm)
+ ecation_prot=enetb(42,i,iparm)
+ evdwpp = enetb(26,i,iparm)
+ eespp = enetb(27,i,iparm)
+ evdwpsb = enetb(28,i,iparm)
+ eelpsb = enetb(29,i,iparm)
+ evdwsb = enetb(30,i,iparm)
+ eelsb = enetb(31,i,iparm)
+ estr_nucl = enetb(32,i,iparm)
+ ebe_nucl = enetb(33,i,iparm)
+ esbloc = enetb(34,i,iparm)
+ etors_nucl = enetb(35,i,iparm)
+ etors_d_nucl = enetb(36,i,iparm)
+ ecorr_nucl = enetb(37,i,iparm)
+ ecorr3_nucl = enetb(38,i,iparm)
+ epeppho= enetb(49,i,iparm)
+ escpho= enetb(48,i,iparm)
+ epepbase= enetb(47,i,iparm)
+ escbase= enetb(46,i,iparm)
+ ecation_nucl=enetb(50,i,iparm)
+ elipbond=enetb(52,i,iparm)
+ elipang=enetb(53,i,iparm)
+ eliplj=enetb(54,i,iparm)
+ elipelec=enetb(55,i,iparm)
+ ecat_prottran=enetb(56,i,iparm)
+ ecation_protang=enetb(57,i,iparm)
+ eliptran=enetb(22,i,iparm)
+ ehomology_constr=enetb(51,i,iparm)
+#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+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*eello_turn6+ft(2)*wel_loc*eel_loc &
+ +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
+ +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
+ +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
+ +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
+ +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
+ *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
+ +wcorr3_nucl*ecorr3_nucl&
+ +wscbase*escbase&
+ +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho+ wcatnucl*ecation_nucl&
+ +elipbond+elipang+eliplj+elipelec+wcat_tran*ecat_prottran+ecation_protang&
+ +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+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*eello_turn6+ft(2)*wel_loc*eel_loc+edihcnstr &
+ +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
+ +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
+ +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
+ +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
+ +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
+ *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
+ +wcorr3_nucl*ecorr3_nucl&
+ +wscbase*escbase&
+ +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho+ wcatnucl*ecation_nucl&
+ +elipbond+elipang+eliplj+elipelec+wcat_tran*ecat_prottran+ecation_protang&
+ +wliptran*eliptran
+
+
+
+#endif
+ write (iout,*) "WTF,",etot,potEmin_all(ib,iparm),entfac(i)/beta_h(ib,iparm)
+ etot=etot-entfac(i)/beta_h(ib,iparm)
+ if(etot.lt.potEmin_all(ib,iparm)) potEmin_all(ib,iparm)=etot
+
+ enddo ! ib
+ enddo ! iparm
+ enddo ! i
+#ifdef DEBUG
+ write (iout,*) "The potEmin array before reduction"
+ do i=1,nParmSet
+ write (iout,*) "Parameter set",i
+ do j=1,nT_h(i)
+ write (iout,*) j,PotEmin_all(j,i)
+ enddo
+ enddo
+ write (iout,*) "potEmin_min",potEmin_min
+#endif
+#ifdef MPI
+!C Determine the minimum energes for all parameter sets and temperatures
+ call MPI_AllReduce(potEmin_all(1,1),potEmin_t_all(1,1), &
+ maxT_h*nParmSet,MPI_DOUBLE_PRECISION,MPI_MIN,WHAM_COMM,IERROR)
+ do i=1,nParmSet
+ do j=1,nT_h(i)
+ potEmin_all(j,i)=potEmin_t_all(j,i)
+ enddo
+ enddo
+#endif
+ potEmin_min=potEmin_all(1,1)
+ do i=1,nParmSet
+ do j=1,nT_h(i)
+ if (potEmin_all(j,i).lt.potEmin_min) &
+ potEmin_min=potEmin_all(j,i)
+ enddo
+ enddo
+#ifdef DEBUG
+ write (iout,*) "The potEmin array"
+ do i=1,nParmSet
+ write (iout,*) "Parameter set",i
+ do j=1,nT_h(i)
+ write (iout,*) j,1.0d0/(1.987d-3*beta_h(j,i)), &
+ PotEmin_all(j,i)
+ enddo
+ enddo
+ write (iout,*) "potEmin_min",potEmin_min
+#endif
+
#ifdef MPI
do t=0,tmax
hfin_ent_p(t)=0.0d0
edihcnstr=0.0d0
ecationcation=enetb(41,t,iparm)
ecation_prot=enetb(42,t,iparm)
+ evdwpp = enetb(26,t,iparm)
+ eespp = enetb(27,t,iparm)
+ evdwpsb = enetb(28,t,iparm)
+ eelpsb = enetb(29,t,iparm)
+ evdwsb = enetb(30,t,iparm)
+ eelsb = enetb(31,t,iparm)
+ estr_nucl = enetb(32,t,iparm)
+ ebe_nucl = enetb(33,t,iparm)
+ esbloc = enetb(34,t,iparm)
+ etors_nucl = enetb(35,t,iparm)
+ etors_d_nucl = enetb(36,t,iparm)
+ ecorr_nucl = enetb(37,t,iparm)
+ ecorr3_nucl = enetb(38,t,iparm)
+ epeppho= enetb(49,t,iparm)
+ escpho= enetb(48,t,iparm)
+ epepbase= enetb(47,t,iparm)
+ escbase= enetb(46,t,iparm)
+ ecation_nucl=enetb(50,t,iparm)
+ elipbond=enetb(52,t,iparm)
+ elipang=enetb(53,t,iparm)
+ eliplj=enetb(54,t,iparm)
+ elipelec=enetb(55,t,iparm)
+ ecat_prottran=enetb(56,t,iparm)
+ ecation_protang=enetb(57,t,iparm)
+ eliptran=enetb(22,t,iparm)
+ ehomology_constr=enetb(51,t,iparm)
do k=0,nGridT
betaT=startGridT+k*delta_T
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 &
+ +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*eello_turn6+ft(2)*wel_loc*eel_loc &
+edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
- +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation
+ +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
+ +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
+ +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
+ +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
+ *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
+ +wcorr3_nucl*ecorr3_nucl&
+ +wscbase*escbase&
+ +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho &
+ + wcatnucl*ecation_nucl&
+ +elipbond+elipang+eliplj+elipelec+wcat_tran*ecat_prottran+ecation_protang&
+ +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*eello_turn6+ &
ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+ &
- ftprim(1)*wsccor*esccor
+ ftprim(1)*wsccor*esccor +ftprim(1)*wtor_nucl*etors_nucl&
+ +wtor_d_nucl*ftprim(2)*etors_d_nucl+ftprim(1)*wpepbase*epepbase
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*eello_turn6+ &
ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+ &
- ftbis(1)*wsccor*esccor
+ ftbis(1)*wsccor*esccor+ftbis(1)*wtor_nucl*etors_nucl&
+ +wtor_d_nucl*ftbis(2)*etors_d_nucl+ftbis(1)*wpepbase*epepbase
#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 &
+ +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*eello_turn6+ft(2)*wel_loc*eel_loc+edihcnstr &
+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
- +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation
+ +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
+ +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
+ +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
+ +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
+ *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
+ +wcorr3_nucl*ecorr3_nucl&
+ +wscbase*escbase&
+ +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho&
+ + wcatnucl*ecation_nucl&
+ +elipbond+elipang+eliplj+elipelec+wcat_tran*ecat_prottran+ecation_protang&
+ +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*eello_turn6+ &
ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+ &
- ftprim(1)*wsccor*esccor
+ ftprim(1)*wsccor*esccor+ftprim(1)*wtor_nucl*etors_nucl&
+ +wtor_d_nucl*ftprim(2)*etors_d_nucl+ftprim(1)*wpepbase*epepbase
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*eello_turn6+ &
ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+ &
- ftprim(1)*wsccor*esccor
-#endif
- weight=dexp(-betaT*(etot-potEmin)+entfac(t))
-#ifdef DEBUG
- write (iout,*) "iparm",iparm," t",t," betaT",betaT,&
- " etot",etot," entfac",entfac(t),&
- " weight",weight," ebis",ebis
-#endif
- etot=etot-temper*eprim
-#ifdef MPI
- sumW_p(k,iparm)=sumW_p(k,iparm)+weight
- sumE_p(k,iparm)=sumE_p(k,iparm)+etot*weight
- sumEbis_p(k,iparm)=sumEbis_p(k,iparm)+ebis*weight
- sumEsq_p(k,iparm)=sumEsq_p(k,iparm)+etot**2*weight
- do j=1,nQ+2
- sumQ_p(j,k,iparm)=sumQ_p(j,k,iparm)+q(j,t)*weight
- sumQsq_p(j,k,iparm)=sumQsq_p(j,k,iparm)+q(j,t)**2*weight
- sumEQ_p(j,k,iparm)=sumEQ_p(j,k,iparm) &
- +etot*q(j,t)*weight
- enddo
-#else
- sumW(k,iparm)=sumW(k,iparm)+weight
- sumE(k,iparm)=sumE(k,iparm)+etot*weight
- sumEbis(k,iparm)=sumEbis(k,iparm)+ebis*weight
- sumEsq(k,iparm)=sumEsq(k,iparm)+etot**2*weight
- do j=1,nQ+2
- sumQ(j,k,iparm)=sumQ(j,k,iparm)+q(j,t)*weight
- sumQsq(j,k,iparm)=sumQsq(j,k,iparm)+q(j,t)**2*weight
- sumEQ(j,k,iparm)=sumEQ(j,k,iparm) &
- +etot*q(j,t)*weight
- enddo
+ ftbis(1)*wsccor*esccor+ftbis(1)*wtor_nucl*etors_nucl&
+ +wtor_d_nucl*ftbis(2)*etors_d_nucl*ftbis(1)*wpepbase*epepbase
+
#endif
+! weight=dexp(-betaT*(etot-potEmin)+entfac(t))
+!#ifdef DEBUG
+! write (iout,*) "iparm",iparm," t",t," betaT",betaT,&
+! " etot",etot," entfac",entfac(t),&
+! " weight",weight," ebis",ebis
+!#endif
+! etot=etot-temper*eprim
+!#ifdef MPI
+! sumW_p(k,iparm)=sumW_p(k,iparm)+weight
+! sumE_p(k,iparm)=sumE_p(k,iparm)+etot*weight
+! sumEbis_p(k,iparm)=sumEbis_p(k,iparm)+ebis*weight
+! sumEsq_p(k,iparm)=sumEsq_p(k,iparm)+etot**2*weight
+! do j=1,nQ+2
+! sumQ_p(j,k,iparm)=sumQ_p(j,k,iparm)+q(j,t)*weight
+! sumQsq_p(j,k,iparm)=sumQsq_p(j,k,iparm)+q(j,t)**2*weight
+! sumEQ_p(j,k,iparm)=sumEQ_p(j,k,iparm) &
+! +etot*q(j,t)*weight
+! enddo
+!#else
+! sumW(k,iparm)=sumW(k,iparm)+weight
+! sumE(k,iparm)=sumE(k,iparm)+etot*weight
+! sumEbis(k,iparm)=sumEbis(k,iparm)+ebis*weight
+! sumEsq(k,iparm)=sumEsq(k,iparm)+etot**2*weight
+! do j=1,nQ+2
+! sumQ(j,k,iparm)=sumQ(j,k,iparm)+q(j,t)*weight
+! sumQsq(j,k,iparm)=sumQsq(j,k,iparm)+q(j,t)**2*weight
+! sumEQ(j,k,iparm)=sumEQ(j,k,iparm) &
+! +etot*q(j,t)*weight
+! enddo
+!#endif
enddo
indE = aint(potE(t,iparm)-aint(potEmin))
if (indE.ge.0 .and. indE.le.maxinde) then
endif
#ifdef MPI
do ib=1,nT_h(iparm)
+ potEmin=potEmin_all(ib,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))
enddo
#else
do ib=1,nT_h(iparm)
+ potEmin=potEmin_all(ib,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))
enddo
#endif
enddo ! t
+!
+! Thermo and ensemble averages
+!
+ do k=0,nGridT
+ betaT=startGridT+k*delta_T
+! call temp_scalfac(betaT,ft,ftprim,ftbis,*10)
+ 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
+
+! write (iout,*) "ftprim",ftprim
+! write (iout,*) "ftbis",ftbis
+ betaT=1.0d0/(1.987D-3*betaT)
+! 7/10/18 AL Determine the max Botzmann weights for each temerature
+! call sum_ene(1,iparm,ft,etot)
+ evdw=enetb(1,1,iparm)
+ evdw_t=enetb(20,1,iparm)
+#ifdef SCP14
+ evdw2_14=enetb(18,1,iparm)
+ evdw2=enetb(2,1,iparm)+evdw2_14
+#else
+ evdw2=enetb(2,1,iparm)
+ evdw2_14=0.0d0
+#endif
+#ifdef SPLITELE
+ ees=enetb(3,1,iparm)
+ evdw1=enetb(16,1,iparm)
+#else
+ ees=enetb(3,1,iparm)
+ evdw1=0.0d0
+#endif
+ ecorr=enetb(4,1,iparm)
+ ecorr5=enetb(5,1,iparm)
+ ecorr6=enetb(6,1,iparm)
+ eel_loc=enetb(7,1,iparm)
+ eello_turn3=enetb(8,1,iparm)
+ eello_turn4=enetb(9,1,iparm)
+ eello_turn6=enetb(10,1,iparm)
+ ebe=enetb(11,1,iparm)
+ escloc=enetb(12,1,iparm)
+ etors=enetb(13,1,iparm)
+ etors_d=enetb(14,1,iparm)
+ ehpb=enetb(15,1,iparm)
+ estr=enetb(17,1,iparm)
+ esccor=enetb(21,1,iparm)
+ edihcnstr=enetb(19,1,iparm)
+! eliptran=enetb(22,1,iparm)
+! esaxs=enetb(26,1,iparm)
+ ecationcation=enetb(41,1,iparm)
+ ecation_prot=enetb(42,1,iparm)
+ evdwpp = enetb(26,1,iparm)
+ eespp = enetb(27,1,iparm)
+ evdwpsb = enetb(28,1,iparm)
+ eelpsb = enetb(29,1,iparm)
+ evdwsb = enetb(30,1,iparm)
+ eelsb = enetb(31,1,iparm)
+ estr_nucl = enetb(32,1,iparm)
+ ebe_nucl = enetb(33,1,iparm)
+ esbloc = enetb(34,1,iparm)
+ etors_nucl = enetb(35,1,iparm)
+ etors_d_nucl = enetb(36,1,iparm)
+ ecorr_nucl = enetb(37,1,iparm)
+ ecorr3_nucl = enetb(38,1,iparm)
+ epeppho= enetb(49,1,iparm)
+ escpho= enetb(48,1,iparm)
+ epepbase= enetb(47,1,iparm)
+ escbase= enetb(46,1,iparm)
+ ecation_nucl= enetb(50,1,iparm)
+ ehomology_constr=enetb(51,1,iparm)
+
+#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*eello_turn6+ft(2)*wel_loc*eel_loc &
+ +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
+ +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
+ +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
+ +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
+ +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
+ *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
+ +wcorr3_nucl*ecorr3_nucl&
+ +wscbase*escbase&
+ +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho &
+ + wcatnucl*ecation_nucl
+
+
+ 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
+ +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*eello_turn6+ft(2)*wel_loc*eel_loc &
+ +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
+ +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
+ +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
+ +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
+ +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
+ *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
+ +wcorr3_nucl*ecorr3_nucl&
+ +wscbase*escbase&
+ +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho &
+ + wcatnucl*ecation_nucl
+
+ 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 &
+! & +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*eello_turn6+ft(2)*wel_loc*eel_loc+edihcnstr &
+ +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
+ +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
+ +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
+ +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
+ +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
+ *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
+ +wcorr3_nucl*ecorr3_nucl&
+ +wscbase*escbase&
+ +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho &
+ + wcatnucl*ecation_nucl
+
+ 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
+ +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*eello_turn6+ft(2)*wel_loc*eel_loc+edihcnstr &
+ +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
+ +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
+ +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
+ +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
+ +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
+ *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
+ +wcorr3_nucl*ecorr3_nucl&
+ +wscbase*escbase&
+ +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho &
+ + wcatnucl*ecation_nucl
+
+ endif
+#endif
+
+ weimax(k,iparm)=-betaT*(etot-potEmin)+entfac(1)
+! write (iout,*) "k",k," t",1," weight",weimax(k,iparm)
+#ifdef MPI
+ do t=2,scount(me1)
+#else
+ do t=2,ntot(islice)
+#endif
+! call sum_ene(t,iparm,ft,etot)
+ evdw=enetb(1,t,iparm)
+ evdw_t=enetb(20,t,iparm)
+#ifdef SCP14
+ evdw2_14=enetb(18,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)
+ eello_turn6=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(17,t,iparm)
+ esccor=enetb(21,t,iparm)
+ edihcnstr=enetb(19,t,iparm)
+! eliptran=enetb(22,t,iparm)
+! esaxs=enetb(26,t,iparm)
+ ecationcation=enetb(41,t,iparm)
+ ecation_prot=enetb(42,t,iparm)
+ evdwpp = enetb(26,t,iparm)
+ eespp = enetb(27,t,iparm)
+ evdwpsb = enetb(28,t,iparm)
+ eelpsb = enetb(29,t,iparm)
+ evdwsb = enetb(30,t,iparm)
+ eelsb = enetb(31,t,iparm)
+ estr_nucl = enetb(32,t,iparm)
+ ebe_nucl = enetb(33,t,iparm)
+ esbloc = enetb(34,t,iparm)
+ etors_nucl = enetb(35,t,iparm)
+ etors_d_nucl = enetb(36,t,iparm)
+ ecorr_nucl = enetb(37,t,iparm)
+ ecorr3_nucl = enetb(38,t,iparm)
+ epeppho= enetb(49,t,iparm)
+ escpho= enetb(48,t,iparm)
+ epepbase= enetb(47,t,iparm)
+ escbase= enetb(46,t,iparm)
+ ecation_nucl=enetb(50,t,iparm)
+ ehomology_constr=enetb(51,t,iparm)
+
+#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*eello_turn6+ft(2)*wel_loc*eel_loc &
+ +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
+ +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
+ +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
+ +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
+ +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
+ *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
+ +wcorr3_nucl*ecorr3_nucl&
+ +wscbase*escbase&
+ +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho &
+ + wcatnucl*ecation_nucl
+
+ 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
+ +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*eello_turn6+ft(2)*wel_loc*eel_loc &
+ +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
+ +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
+ +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
+ +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
+ +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
+ *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
+ +wcorr3_nucl*ecorr3_nucl&
+ +wscbase*escbase&
+ +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho &
+ + wcatnucl*ecation_nucl
+ 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 &
+! & +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*eello_turn6+ft(2)*wel_loc*eel_loc+edihcnstr &
+ +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
+ +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
+ +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
+ +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
+ +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
+ *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
+ +wcorr3_nucl*ecorr3_nucl&
+ +wscbase*escbase&
+ +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho &
+ + wcatnucl*ecation_nucl
+ 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
+ +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*eello_turn6+ft(2)*wel_loc*eel_loc+edihcnstr &
+ +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
+ +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
+ +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
+ +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
+ +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
+ *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
+ +wcorr3_nucl*ecorr3_nucl&
+ +wscbase*escbase&
+ +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho &
+ + wcatnucl*ecation_nucl
+ endif
+#endif
+
+ weight=-betaT*(etot-potEmin)+entfac(t)
+! 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
+ 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
+ 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
+ do k=0,nGridT
+ temper=startGridT+k*delta_T
+ betaT=1.0d0/(1.987D-3*temper)
+! call temp_scalfac(temper,ft,ftprim,ftbis,*10)
+ if (rescale_mode.eq.1) then
+ quot=temper/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((temper-320.0d0)/80.0d0))/ &
+ 320.0d0
+ ftprim(6)=1.0d0/(320.0d0*dcosh((temper-320.0d0)/80.0d0)**2)
+ ftbis(6)=-2.0d0*dtanh((temper-320.0d0)/80.0d0) &
+ /(320.0d0*80.0d0*dcosh((temper-320.0d0)/80.0d0)**3)
+#elif defined(FUNCT)
+ fT(6)=temper/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=temper/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((temper-320.0d0)/80.0d0))/ &
+ 320.0d0
+ ftprim(6)=1.0d0/(320.0d0*dcosh((temper-320.0d0)/80.0d0)**2)
+ ftbis(6)=-2.0d0*dtanh((temper-320.0d0)/80.0d0) &
+ /(320.0d0*80.0d0*dcosh((temper-320.0d0)/80.0d0)**3)
+#elif defined(FUNCT)
+ fT(6)=temper/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
+
+
+
+ do t=1,scount(me1)
+#else
+ do t=1,ntot(islice)
+#endif
+ ind = ind_point(t)
+#ifdef MPI
+ hfin_ent_p(ind)=hfin_ent_p(ind)+dexp(entfac(t))
+#else
+ hfin_ent(ind)=hfin_ent(ind)+dexp(entfac(t))
+#endif
+! write (iout,'(2i5,20f8.2)') t,t,(enetb(k,t,iparm),k=1,18)
+! call restore_parm(iparm)
+! call sum_ene_deriv(t,iparm,ft,ftprim,ftbis,etot,eprim,ebis)
+ evdw=enetb(1,t,iparm)
+ evdw_t=enetb(20,t,iparm)
+#ifdef SCP14
+ evdw2_14=enetb(18,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)
+ eello_turn6=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(17,t,iparm)
+ esccor=enetb(21,t,iparm)
+ edihcnstr=enetb(19,t,iparm)
+! eliptran=enetb(22,t,iparm)
+! esaxs=enetb(26,t,iparm)
+! ehomology_constr=enetb(27,t,iparm)
+! edfadis=enetb(28,t,iparm)
+! edfator=enetb(29,t,iparm)
+! edfanei=enetb(30,t,iparm)
+! edfabet=enetb(31,t,iparm)
+ ecationcation=enetb(41,i,iparm)
+ ecation_prot=enetb(42,i,iparm)
+ evdwpp = enetb(26,i,iparm)
+ eespp = enetb(27,i,iparm)
+ evdwpsb = enetb(28,i,iparm)
+ eelpsb = enetb(29,i,iparm)
+ evdwsb = enetb(30,i,iparm)
+ eelsb = enetb(31,i,iparm)
+ estr_nucl = enetb(32,i,iparm)
+ ebe_nucl = enetb(33,i,iparm)
+ esbloc = enetb(34,i,iparm)
+ etors_nucl = enetb(35,i,iparm)
+ etors_d_nucl = enetb(36,i,iparm)
+ ecorr_nucl = enetb(37,i,iparm)
+ ecorr3_nucl = enetb(38,i,iparm)
+ epeppho= enetb(49,i,iparm)
+ escpho= enetb(48,i,iparm)
+ epepbase= enetb(47,i,iparm)
+ escbase= enetb(46,i,iparm)
+ ecation_nucl= enetb(50,i,iparm)
+ ehomology_constr=enetb(51,i,iparm)
+
+#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*eello_turn6+ft(2)*wel_loc*eel_loc &
+ +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
+ +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation &
+ +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
+ +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
+ +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
+ *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
+ +wcorr3_nucl*ecorr3_nucl&
+ +wscbase*escbase&
+ +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho &
+ + wcatnucl*ecation_nucl
+
+! eprim=ftprim(1)*(ft(6)*evdw_t+evdw)
+!! & +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*eello_turn6+
+! & 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*eello_turn6+
+! & ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
+! & ftbis(1)*wsccor*esccor
+
+ 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*eello_turn6+ &
+ ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+ &
+ ftprim(1)*wsccor*esccor+ftprim(1)*wtor_nucl*etors_nucl&
+ +wtor_d_nucl*ftprim(2)*etors_d_nucl+ftprim(1)*wpepbase*epepbase
+ 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*eello_turn6+ &
+ ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+ &
+ ftbis(1)*wsccor*esccor+ftbis(1)*wtor_nucl*etors_nucl&
+ +wtor_d_nucl*ftbis(2)*etors_d_nucl*ftbis(1)*wpepbase*epepbase
+
+
+
+ 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*eello_turn6+ft(2)*wel_loc*eel_loc &
+ +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
+ +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
+ +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
+ +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
+ +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
+ *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
+ +wcorr3_nucl*ecorr3_nucl&
+ +wscbase*escbase&
+ +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho&
+ + wcatnucl*ecation_nucl
+
+
+ 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*eello_turn6+ &
+ ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+ &
+ ftprim(1)*wsccor*esccor+ftprim(1)*wtor_nucl*etors_nucl&
+ +wtor_d_nucl*ftprim(2)*etors_d_nucl+ftprim(1)*wpepbase*epepbase
+ 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*eello_turn6+ &
+ ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+ &
+ ftbis(1)*wsccor*esccor+ftbis(1)*wtor_nucl*etors_nucl&
+ +wtor_d_nucl*ftbis(2)*etors_d_nucl*ftbis(1)*wpepbase*epepbase
+
+
+
+! 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*eello_turn6+
+! & 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*eello_turn6+
+! & 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*eello_turn6+ft(2)*wel_loc*eel_loc+edihcnstr &
+ +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
+ +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
+ +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
+ +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
+ +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
+ *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
+ +wcorr3_nucl*ecorr3_nucl&
+ +wscbase*escbase&
+ +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho&
+ + wcatnucl*ecation_nucl
+
+
+ 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*eello_turn6+ &
+ ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+ &
+ ftprim(1)*wsccor*esccor+ftprim(1)*wtor_nucl*etors_nucl&
+ +wtor_d_nucl*ftprim(2)*etors_d_nucl+ftprim(1)*wpepbase*epepbase
+ 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*eello_turn6+ &
+ ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+ &
+ ftbis(1)*wsccor*esccor+ftbis(1)*wtor_nucl*etors_nucl&
+ +wtor_d_nucl*ftbis(2)*etors_d_nucl*ftbis(1)*wpepbase*epepbase
+
+
+
+
+! 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*eello_turn6+
+! & 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*eello_turn6+
+! & 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*eello_turn6+ft(2)*wel_loc*eel_loc+edihcnstr &
+ +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor &
+ +wbond*estr+wcatprot*ecation_prot+wcatcat*ecationcation&
+ +wbond_nucl*estr_nucl+wang_nucl*ebe_nucl&
+ +wvdwpp_nucl*evdwpp+welpp*eespp+wvdwpsb*evdwpsb+welpsb*eelpsb&
+ +wvdwsb*evdwsb+welsb*eelsb+wsbloc*esbloc+wtor_nucl*etors_nucl&
+ *ft(1)+wtor_d_nucl*ft(2)*etors_d_nucl+wcorr_nucl*ecorr_nucl &
+ +wcorr3_nucl*ecorr3_nucl&
+ +wscbase*escbase&
+ +ft(1)*wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho&
+ + wcatnucl*ecation_nucl
+
+
+ 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*eello_turn6+ &
+ ftprim(2)*wel_loc*eel_loc+ftprim(2)*wtor_d*etors_d+ &
+ ftprim(1)*wsccor*esccor+ftprim(1)*wtor_nucl*etors_nucl&
+ +wtor_d_nucl*ftprim(2)*etors_d_nucl+ftprim(1)*wpepbase*epepbase
+ 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*eello_turn6+ &
+ ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+ &
+ ftbis(1)*wsccor*esccor+ftbis(1)*wtor_nucl*etors_nucl&
+ +wtor_d_nucl*ftbis(2)*etors_d_nucl*ftbis(1)*wpepbase*epepbase
+
+
+
+
+
+! 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*eello_turn6+
+! & 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*eello_turn6+
+! & ftbis(2)*wel_loc*eel_loc+ftbis(2)*wtor_d*etors_d+
+! & ftprim(1)*wsccor*esccor
+ endif
+#endif
+
+ 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)," boltz", &
+ -betaT*(etot-potEmin)+entfac(t)," weimax",weimax(k,iparm), &
+ " weight",weight," ebis",ebis
+#endif
+ etot=etot-temper*eprim
+#ifdef MPI
+ sumW_p(k,iparm)=sumW_p(k,iparm)+weight
+ sumE_p(k,iparm)=sumE_p(k,iparm)+etot*weight
+ sumEbis_p(k,iparm)=sumEbis_p(k,iparm)+ebis*weight
+ sumEsq_p(k,iparm)=sumEsq_p(k,iparm)+etot**2*weight
+ do j=1,nQ+2
+ sumQ_p(j,k,iparm)=sumQ_p(j,k,iparm)+q(j,t)*weight
+ sumQsq_p(j,k,iparm)=sumQsq_p(j,k,iparm)+q(j,t)**2*weight
+ sumEQ_p(j,k,iparm)=sumEQ_p(j,k,iparm) &
+ +etot*q(j,t)*weight
+ enddo
+#else
+ sumW(k,iparm)=sumW(k,iparm)+weight
+ sumE(k,iparm)=sumE(k,iparm)+etot*weight
+ sumEbis(k,iparm)=sumEbis(k,iparm)+ebis*weight
+ sumEsq(k,iparm)=sumEsq(k,iparm)+etot**2*weight
+ do j=1,nQ+2
+ sumQ(j,k,iparm)=sumQ(j,k,iparm)+q(j,t)*weight
+ sumQsq(j,k,iparm)=sumQsq(j,k,iparm)+q(j,t)**2*weight
+ sumEQ(j,k,iparm)=sumEQ(j,k,iparm) &
+ +etot*q(j,t)*weight
+ enddo
+#endif
+ 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)
write (iout,'(a,i3)') "Parameter set",iparm
endif
do i=0,NGridT
+ betaT=1.0d0/(1.987D-3*(startGridT+i*delta_T))
+ if (betaT.ge.beta_h(1,iparm)) then
+ potEmin=potEmin_all(1,iparm)
+ else if (betaT.lt.beta_h(nT_h(iparm),iparm)) then
+ potEmin=potEmin_all(nT_h(iparm),iparm)
+ else
+ do l=1,nT_h(iparm)-1
+ if (betaT.le.beta_h(l,iparm) .and. &
+ betaT.gt.beta_h(l+1,iparm)) then
+ potEmin=potEmin_all(l,iparm)
+ exit
+ endif
+ enddo
+ endif
+
sumE(i,iparm)=sumE(i,iparm)/sumW(i,iparm)
sumEbis(i,iparm)=(startGridT+i*delta_T)*sumEbis(i,iparm)/ &
sumW(i,iparm)
projects / unres4.git / blobdiff