X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fwham_calc.F90;h=81093c83d95c3b480f37bf613b16c389a8101951;hb=f0d780276769530461648dd4c4ad83650d31cf7c;hp=d1ed2f7b8b5f02df0019c64b7bce9643a9f0f37a;hpb=db2040af09e73f3fc871158afe11ee242ec6cfb1;p=unres4.git diff --git a/source/wham/wham_calc.F90 b/source/wham/wham_calc.F90 index d1ed2f7..81093c8 100644 --- a/source/wham/wham_calc.F90 +++ b/source/wham/wham_calc.F90 @@ -71,7 +71,8 @@ !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,& @@ -80,7 +81,7 @@ hfin_ent_p(0:MaxHdim),histE_p(0:maxindE),sumH,& hrmsrgy_p(0:MaxBinRgy,0:MaxBinRms,maxT_h) 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 @@ -96,14 +97,18 @@ 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 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 + integer :: ind_point(maxpoint),upindE,indE character(len=16) :: plik @@ -125,6 +130,12 @@ 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 @@ -204,8 +215,8 @@ 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,& @@ -214,7 +225,7 @@ 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 @@ -245,16 +256,16 @@ #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 ! @@ -345,6 +356,24 @@ 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) + #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,& @@ -382,7 +411,15 @@ +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 + + #else etot=wsc*evdw+wscp*evdw2 & +welec*(ees+evdw1) & @@ -392,7 +429,15 @@ +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 + + #endif #ifdef DEBUG @@ -459,6 +504,7 @@ 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 @@ -467,20 +513,79 @@ 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 @@ -498,6 +603,13 @@ 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,& @@ -521,7 +633,7 @@ 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 @@ -571,6 +683,211 @@ 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) + +#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*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 + +#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*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 + +#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 @@ -757,6 +1074,24 @@ 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) + do k=0,nGridT betaT=startGridT+k*delta_T @@ -848,20 +1183,30 @@ +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 + 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) & @@ -871,20 +1216,31 @@ +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 + 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 + 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 @@ -924,6 +1280,7 @@ 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)) @@ -936,6 +1293,7 @@ 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)) @@ -1154,6 +1512,21 @@ 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)