X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fwham_calc.F90;h=8a25261317b12589fa8c1e40df54b3038784baa2;hb=4db3a5f0b00eb4e0dc343054d4560d5e76f548f5;hp=81093c83d95c3b480f37bf613b16c389a8101951;hpb=f0d780276769530461648dd4c4ad83650d31cf7c;p=unres4.git diff --git a/source/wham/wham_calc.F90 b/source/wham/wham_calc.F90 index 81093c8..8a25261 100644 --- a/source/wham/wham_calc.F90 +++ b/source/wham/wham_calc.F90 @@ -79,7 +79,7 @@ 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,potEmin_t_all(maxT_h,Max_Parm)!,entmin_p,entmax_p ! integer :: histent_p(0:2000) @@ -98,7 +98,8 @@ hrmsrgy(0:MaxBinRgy,0:MaxBinRms,maxT_h),& potEmin,ent,& hfin_ent(0:MaxHdim),vmax,aux,fi_min(MaxR,maxT_h,nParmSet), & - potEmin_all(maxT_h,Max_Parm),potEmin_min,entfac_min + 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 @@ -107,7 +108,9 @@ eello_turn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor, & 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 + 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 @@ -373,11 +376,20 @@ 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 @@ -387,7 +399,7 @@ ! +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 @@ -397,7 +409,7 @@ ! +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 @@ -406,7 +418,7 @@ 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 & @@ -417,14 +429,16 @@ +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 + +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 & @@ -435,7 +449,11 @@ +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 + +wpepbase*epepbase+wscpho*escpho+wpeppho*epeppho+wcatnucl*ecation_nucl& + +elipbond+elipang+eliplj+elipelec+wcat_tran*ecat_prottran+ecation_protang& + +wliptran*eliptran + + #endif @@ -455,6 +473,28 @@ 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 @@ -469,6 +509,7 @@ etot,potEmin,etot-potEmin,v(i,kk,ib,iparm) #endif enddo ! kk + endif enddo ! ib enddo ! iparm enddo ! i @@ -804,12 +845,20 @@ 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+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 & @@ -821,13 +870,16 @@ *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 + +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+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 & @@ -839,7 +891,11 @@ *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 + +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) @@ -1091,7 +1147,15 @@ 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 @@ -1178,7 +1242,7 @@ 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 & @@ -1190,7 +1254,12 @@ *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 + +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+ & @@ -1211,7 +1280,7 @@ 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 & @@ -1223,7 +1292,12 @@ *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 + +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+ & @@ -1242,36 +1316,36 @@ +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 +! 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 @@ -1306,6 +1380,825 @@ 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)