X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fsrc%2Fwham_calc1.F;h=5783bd4cdd5c403d911f4b9e6e2e1d80d1e34523;hb=24bbf0f3f06452c69bf4b69b4a9ba620e0260d7b;hp=379e93ba8e715c8851f03a2b20a6855c5ab05e84;hpb=fbc0f12bd8022b2c766eb8c350e7bb0451816c1d;p=unres.git diff --git a/source/wham/src/wham_calc1.F b/source/wham/src/wham_calc1.F index 379e93b..5783bd4 100644 --- a/source/wham/src/wham_calc1.F +++ b/source/wham/src/wham_calc1.F @@ -14,10 +14,8 @@ include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" include "DIMENSIONS.FREE" - integer nGridT - parameter (NGridT=400) integer MaxBinRms,MaxBinRgy - parameter (MaxBinRms=1000,MaxBinRgy=1000) + parameter (MaxBinRms=100,MaxBinRgy=100) integer MaxHdim c parameter (MaxHdim=200000) parameter (MaxHdim=200) @@ -51,22 +49,21 @@ c parameter (MaxHdim=200000) #ifdef MPI integer tmax_t,upindE_p double precision fi_p(MaxR,MaxT_h,Max_Parm) - double precision sumW_p(0:nGridT,Max_Parm), - & sumE_p(0:nGridT,Max_Parm),sumEsq_p(0:nGridT,Max_Parm), - & sumQ_p(MaxQ1,0:nGridT,Max_Parm), - & sumQsq_p(MaxQ1,0:nGridT,Max_Parm), - & sumEQ_p(MaxQ1,0:nGridT,Max_Parm), - & sumEprim_p(MaxQ1,0:nGridT,Max_Parm), - & sumEbis_p(0:nGridT,Max_Parm) + double precision sumW_p(0:Max_GridT,Max_Parm), + & sumE_p(0:Max_GridT,Max_Parm),sumEsq_p(0:Max_GridT,Max_Parm), + & sumQ_p(MaxQ1,0:Max_GridT,Max_Parm), + & sumQsq_p(MaxQ1,0:Max_GridT,Max_Parm), + & sumEQ_p(MaxQ1,0:Max_GridT,Max_Parm), + & sumEprim_p(MaxQ1,0:Max_GridT,Max_Parm), + & sumEbis_p(0:Max_GridT,Max_Parm) double precision hfin_p(0:MaxHdim,maxT_h), & hfin_ent_p(0:MaxHdim),histE_p(0:maxindE),sumH, & hrmsrgy_p(0:MaxBinRgy,0:MaxBinRms,maxT_h) double precision rgymin_t,rmsmin_t,rgymax_t,rmsmax_t - double precision potEmin_t,entmin_p,entmax_p + double precision potEmin_t_all(maxT_h,Max_Parm),entmin_p,entmax_p integer histent_p(0:2000) logical lprint /.true./ #endif - double precision delta_T /1.0d0/ double precision rgymin,rmsmin,rgymax,rmsmax double precision sumW(0:NGridT,Max_Parm),sumE(0:NGridT,Max_Parm), & sumEsq(0:NGridT,Max_Parm),sumQ(MaxQ1,0:NGridT,Max_Parm), @@ -77,10 +74,10 @@ c parameter (MaxHdim=200000) & dd,dd1,dd2,hh,dmin,denom,finorm,avefi,pom, & hfin(0:MaxHdim,maxT_h),histE(0:maxindE), & hrmsrgy(0:MaxBinRgy,0:MaxBinRms,maxT_h), - & potEmin,ent, + & potEmin_all(maxT_h,Max_Parm),potEmin,potEmin_min,ent, & hfin_ent(0:MaxHdim),vmax,aux double precision fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl, - & eprim,ebis,temper,kfac/2.4d0/,T0/300.0d0/,startGridT/200.0d0/, + & eprim,ebis,temper,kfac/2.4d0/,T0/300.0d0/, & eplus,eminus,logfac,tanhT,tt double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors, & escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3, @@ -104,7 +101,11 @@ c parameter (MaxHdim=200000) call flush(iout) dmin=0.0d0 tmax=0 - potEmin=1.0d10 + do i=1,nParmset + do j=1,nT_h(i) + potEmin_all(j,i)=1.0d10 + enddo + enddo rgymin=1.0d10 rmsmin=1.0d10 rgymax=0.0d0 @@ -117,9 +118,6 @@ c parameter (MaxHdim=200000) #else do i=1,ntot(islice) #endif - do j=1,nParmSet - if (potE(i,j).le.potEmin) potEmin=potE(i,j) - enddo if (q(nQ+1,i).lt.rmsmin) rmsmin=q(nQ+1,i) if (q(nQ+1,i).gt.rmsmax) rmsmax=q(nQ+1,i) if (q(nQ+2,i).lt.rgymin) rgymin=q(nQ+2,i) @@ -132,9 +130,6 @@ c parameter (MaxHdim=200000) else ind_point(i)=ind_point(i)+nbin1**(j-1)*ind endif -c write (iout,*) "i",i," j",j," q",q(j,i)," ind_point", -c & ind_point(i) - call flush(iout) if (ind_point(i).lt.0 .or. ind_point(i).gt.MaxHdim) then write (iout,*) "Error - index exceeds range for point",i, & " q=",q(j,i)," ind",ind_point(i) @@ -184,8 +179,6 @@ c & ind_point(i) 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(rmsmin,rmsmin_t,1,MPI_DOUBLE_PRECISION, & MPI_MIN,WHAM_COMM,IERROR) call MPI_AllReduce(rmsmax,rmsmax_t,1,MPI_DOUBLE_PRECISION, @@ -194,12 +187,10 @@ c & ind_point(i) & MPI_MIN,WHAM_COMM,IERROR) call MPI_AllReduce(rgymax,rgymax_t,1,MPI_DOUBLE_PRECISION, & MPI_MAX,WHAM_COMM,IERROR) - potEmin=potEmin_t/2 rgymin=rgymin_t rgymax=rgymax_t rmsmin=rmsmin_t rmsmax=rmsmax_t - write (iout,*) "potEmin",potEmin #endif rmsmin=deltrms*dint(rmsmin/deltrms) rmsmax=deltrms*dint(rmsmax/deltrms) @@ -341,7 +332,7 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft #endif #ifdef DEBUG write (iout,*) i,iparm,1.0d0/(beta_h(ib,iparm)*1.987D-3), - & etot,potEmin + & etot #endif #ifdef DEBUG if (iparm.eq.1 .and. ib.eq.1) then @@ -361,10 +352,10 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft & *(dd-q0(j,kk,ib,iparm))**2 enddo v(i,kk,ib,iparm)= - & -beta_h(ib,iparm)*(etot-potEmin+Econstr) + & -beta_h(ib,iparm)*(etot+Econstr) #ifdef DEBUG write (iout,'(4i5,4e15.5)') i,kk,ib,iparm, - & etot,potEmin,etot-potEmin,v(i,kk,ib,iparm) + & etot,v(i,kk,ib,iparm) #endif enddo ! kk enddo ! ib @@ -434,7 +425,7 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft #ifdef DEBUG write (iout,*) "fi before MPI_Reduce me",me,' master',master do iparm=1,nParmSet - do ib=1,nT_h(nparmset) + do ib=1,nT_h(iparm) 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)) @@ -514,6 +505,177 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft 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,21) +#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_mode.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_mode.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_mode.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) + return1 + endif + evdw=enetb(1,i,iparm) + evdw_t=enetb(21,i,iparm) +#ifdef SCP14 + evdw2_14=enetb(17,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) + eturn6=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) + esccor=enetb(19,i,iparm) + edihcnstr=enetb(20,i,iparm) +#ifdef DEBUG + write (iout,'(3i5,6f5.2,14f12.3)') i,ib,iparm,(ft(l),l=1,6), + & evdw+evdw_t,evdw2,ees,evdw1,ecorr,eel_loc,estr,ebe,escloc, + & etors,etors_d,eello_turn3,eello_turn4,esccor,edihcnstr +#endif + +#ifdef SPLITELE + etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees + & +wvdwpp*evdw1 + & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc + & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 + & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 + & +ft(2)*wturn3*eello_turn3 + & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc + & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor + & +wbond*estr +#else + etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 + & +ft(1)*welec*(ees+evdw1) + & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc + & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 + & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 + & +ft(2)*wturn3*eello_turn3 + & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr + & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor + & +wbond*estr +#endif +c write (iout,*) "i",i," ib",ib, +c & " temp",1.0d0/(1.987d-3*beta_h(ib,iparm))," etot",etot, +c & " entfac",entfac(i) + etot=etot-entfac(i)/beta_h(ib,iparm) + if(etot.lt.potEmin_all(ib,iparm)) potEmin_all(ib,iparm)=etot +c write (iout,*) "efree",etot," potEmin",potEmin_all(ib,iparm) + 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,PotEmin_all(j,i) + enddo + enddo + write (iout,*) "potEmin_min",potEmin_min +#endif +#undef DEBUG + #ifdef MPI do t=0,tmax hfin_ent_p(t)=0.0d0 @@ -653,7 +815,6 @@ c write (iout,*) "me1",me1," scount",scount(me1) #else hfin_ent(ind)=hfin_ent(ind)+dexp(entfac(t)) #endif -c write (iout,'(2i5,20f8.2)') t,t,(enetb(k,t,iparm),k=1,18) call restore_parm(iparm) evdw=enetb(21,t,iparm) evdw_t=enetb(1,t,iparm) @@ -686,7 +847,6 @@ c write (iout,'(2i5,20f8.2)') t,t,(enetb(k,t,iparm),k=1,18) estr=enetb(18,t,iparm) esccor=enetb(19,t,iparm) edihcnstr=enetb(20,t,iparm) - edihcnstr=0.0d0 do k=0,nGridT betaT=startGridT+k*delta_T temper=betaT @@ -767,6 +927,25 @@ c ft=2*T0/(T0+betaT) c write (iout,*) "ftprim",ftprim c write (iout,*) "ftbis",ftbis betaT=1.0d0/(1.987D-3*betaT) + if (betaT.ge.beta_h(1,iparm)) then + potEmin=potEmin_all(1,iparm) +c write(iout,*) "first",temper,potEmin + else if (betaT.lt.beta_h(nT_h(iparm),iparm)) then + potEmin=potEmin_all(nT_h(iparm),iparm) +c write (iout,*) "last",temper,potEmin + 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) +c write (iout,*) "l",l, +c & betaT,1.0d0/(1.987D-3*beta_h(l,iparm)), +c & 1.0d0/(1.987D-3*beta_h(l+1,iparm)),temper,potEmin + exit + endif + enddo + endif +c write (iout,*) ib," PotEmin",potEmin #ifdef SPLITELE etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees & +wvdwpp*evdw1 @@ -816,8 +995,10 @@ c write (iout,*) "ftbis",ftbis #endif weight=dexp(-betaT*(etot-potEmin)+entfac(t)) #ifdef DEBUG - write (iout,*) "iparm",iparm," t",t," betaT",betaT, + write (iout,*) "iparm",iparm," t",t," temper",temper, & " etot",etot," entfac",entfac(t), + & " efree",etot-entfac(t)/betaT," potEmin",potEmin, + & " boltz",-betaT*(etot-potEmin)+entfac(t), & " weight",weight," ebis",ebis #endif etot=etot-temper*eprim @@ -852,6 +1033,7 @@ c write (iout,*) "ftbis",ftbis 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)) @@ -864,6 +1046,7 @@ c write (iout,*) "ftbis",ftbis 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)) @@ -1082,6 +1265,20 @@ c write (iout,*) "ftbis",ftbis 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)