X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fsrc%2Fwham_calc1.F;h=ce5de230423b57da17f79d3b5a10a1a1747d1399;hb=aa037b39c46474f358b050551261d7c5ac3d2357;hp=5783bd4cdd5c403d911f4b9e6e2e1d80d1e34523;hpb=24bbf0f3f06452c69bf4b69b4a9ba620e0260d7b;p=unres.git diff --git a/source/wham/src/wham_calc1.F b/source/wham/src/wham_calc1.F index 5783bd4..ce5de23 100644 --- a/source/wham/src/wham_calc1.F +++ b/source/wham/src/wham_calc1.F @@ -38,6 +38,7 @@ c parameter (MaxHdim=200000) parameter (MaxPoint=MaxStr, & MaxPointProc=MaxStr_Proc) double precision finorm_max,potfac,entmin,entmax,expfac,vf + double precision entfac_min,entfac_min_t parameter (finorm_max=1.0d0) integer islice integer i,ii,j,jj,k,kk,l,m,ind,iter,t,tmax,ient,ientmax,iln @@ -48,7 +49,8 @@ c parameter (MaxHdim=200000) double precision energia(0:max_ene) #ifdef MPI integer tmax_t,upindE_p - double precision fi_p(MaxR,MaxT_h,Max_Parm) + double precision fi_p(MaxR,MaxT_h,Max_Parm), + & fi_p_min(MaxR,MaxT_h,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), @@ -71,6 +73,7 @@ c parameter (MaxHdim=200000) & sumEprim(0:NGridT,Max_Parm),sumEbis(0:NGridT,Max_Parm),betaT, & weight,econstr double precision fi(MaxR,maxT_h,Max_Parm), + & fi_min(MaxR,maxT_h,Max_Parm), & dd,dd1,dd2,hh,dmin,denom,finorm,avefi,pom, & hfin(0:MaxHdim,maxT_h),histE(0:maxindE), & hrmsrgy(0:MaxBinRgy,0:MaxBinRms,maxT_h), @@ -81,7 +84,9 @@ c parameter (MaxHdim=200000) & eplus,eminus,logfac,tanhT,tt double precision etot,evdw,evdw_t,evdw2,ees,evdw1,ebe,etors, & escloc,ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3, - & eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor + & eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor, + & ehomology_constr,edfadis,edfator,edfanei,edfabet + integer ind_point(maxpoint),upindE,indE character*16 plik @@ -91,7 +96,6 @@ c parameter (MaxHdim=200000) character*128 nazwa integer ilen external ilen - write(licz2,'(bz,i2.2)') islice nbin1 = 1.0d0/delta write (iout,'(//80(1h-)/"Solving WHAM equations for slice", @@ -217,8 +221,8 @@ c parameter (MaxHdim=200000) 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) + write (iout,'(2i5,22f8.2)') i,iparm, + & (enetb(k,i,iparm),k=1,22) #endif call restore_parm(iparm) #ifdef DEBUG @@ -303,32 +307,40 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft estr=enetb(18,i,iparm) esccor=enetb(19,i,iparm) edihcnstr=enetb(20,i,iparm) + ehomology_constr=enetb(22,i,iparm) + edfadis=enetb(23,i,iparm) + edfator=enetb(24,i,iparm) + edfanei=enetb(25,i,iparm) + edfabet=enetb(26,i,iparm) #ifdef DEBUG - write (iout,'(3i5,6f5.2,14f12.3)') i,ib,iparm,(ft(l),l=1,6), + write (iout,'(3i5,6f5.2,16f12.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 + & etors,etors_d,eello_turn3,eello_turn4,esccor, + & ehomology_constr,edfadis,edfator,edfanei,edfabet #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 + & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 & +ft(2)*wturn3*eello_turn3 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor - & +wbond*estr + & +wbond*estr+wdfa_dist*edfadis + & +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet #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*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor - & +wbond*estr + & +wbond*estr+wdfa_dist*edfadis + & +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet #endif #ifdef DEBUG write (iout,*) i,iparm,1.0d0/(beta_h(ib,iparm)*1.987D-3), @@ -344,6 +356,24 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft call enerprint(energia(0),fT) endif #endif +#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 @@ -358,6 +388,8 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft & etot,v(i,kk,ib,iparm) #endif enddo ! kk + + endif enddo ! ib enddo ! iparm enddo ! i @@ -368,6 +400,7 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft ! Compute new free-energy values corresponding to the righ-hand side of the ! equation and their derivatives. write (iout,*) "------------------------fi" + entfac_min=1.0d10 #ifdef MPI do t=1,scount(me1) #else @@ -393,10 +426,67 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft 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 enddo +c#ifdef MPI +c write (iout,*) "entfac_min before AllReduce",entfac_min +c call MPI_AllReduce(entfac_min,entfac_min_t,1, +c & MPI_DOUBLE_PRECISION,MPI_MAX,WHAM_COMM,IERROR) +c entfac_min=entfac_min_t +c write (iout,*) "entfac_min after AllReduce",entfac_min +c#endif +c#ifdef MPI +c do t=1,scount(me) +c entfac(t)=entfac(t)-entfac_min +c enddo +c#else +c do t=1,ntot(islice) +c entfac(t)=entfac(t)-entfac_min +c enddo +c#endif + do iparm=1,nParmSet + do iib=1,nT_h(iparm) + do ii=1,nR(iib,iparm) +#ifdef MPI + fi_p_min(ii,iib,iparm)=-1.0d10 + do t=1,scount(me) + aux=v(t,ii,iib,iparm)+entfac(t) + if (aux.gt.fi_p_min(ii,iib,iparm)) + & fi_p_min(ii,iib,iparm)=aux + enddo +#else + do t=1,ntot(islice) + 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_min before AllReduce" + do i=1,nParmSet + do j=1,nT_h(i) + write (iout,*) (i,j,k,fi_p_min(k,j,i),k=1,nR(j,i)) + 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 +#endif do iparm=1,nParmSet do iib=1,nT_h(iparm) do ii=1,nR(iib,iparm) @@ -404,17 +494,18 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft fi_p(ii,iib,iparm)=0.0d0 do t=1,scount(me) fi_p(ii,iib,iparm)=fi_p(ii,iib,iparm) - & +dexp(v(t,ii,iib,iparm)+entfac(t)) + & +dexp(v(t,ii,iib,iparm)+entfac(t)-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 @@ -432,10 +523,12 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft 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 call MPI_Reduce(fi_p(1,1,1),fi(1,1,1),maxR*MaxT_h*nParmSet, & MPI_DOUBLE_PRECISION, & MPI_SUM,Master,WHAM_COMM,IERROR) @@ -455,7 +548,7 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft do iparm=1,nParmSet do ib=1,nT_h(iparm) do i=1,nR(ib,iparm) - fi(i,ib,iparm)=-dlog(fi(i,ib,iparm)) + fi(i,ib,iparm)=-dlog(fi(i,ib,iparm))-fi_min(i,ib,iparm) avefi=avefi+fi(i,ib,iparm) enddo enddo @@ -516,7 +609,7 @@ c write (9,'(3i5,f10.5)') i,(iparm,potE(i,iparm),iparm=1,nParmSet) do iparm=1,nParmSet #ifdef DEBUG write (iout,'(2i5,21f8.2)') i,iparm, - & (enetb(k,i,iparm),k=1,21) + & (enetb(k,i,iparm),k=1,22) #endif call restore_parm(iparm) #ifdef DEBUG @@ -601,39 +694,54 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft estr=enetb(18,i,iparm) esccor=enetb(19,i,iparm) edihcnstr=enetb(20,i,iparm) + ehomology_constr=enetb(22,i,iparm) + if (homol_nset.gt.1) + & ehomology_constr=waga_homology(ihset)*ehomology_constr + edfadis=enetb(23,i,iparm) + edfator=enetb(24,i,iparm) + edfanei=enetb(25,i,iparm) + edfabet=enetb(26,i,iparm) #ifdef DEBUG - write (iout,'(3i5,6f5.2,14f12.3)') i,ib,iparm,(ft(l),l=1,6), + write (iout,'(3i5,6f5.2,17f12.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 + & etors,etors_d,eello_turn3,eello_turn4,esccor,edihcnstr, + & ehomology_constr+wdfa_dist*edfadis+wdfa_tor*edfator+ + & wdfa_nei*edfanei+wdfa_beta*edfabet #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 + & +wstrain*ehpb+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 & +ft(2)*wturn3*eello_turn3 & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor - & +wbond*estr + & +wbond*estr+ehomology_constr+wdfa_dist*edfadis + & +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet #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*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor - & +wbond*estr + & +wbond*estr+ehomology_constr+wdfa_dist*edfadis + & +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet #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) +#ifdef DEBUG + write (iout,*) "i",i," ib",ib, + & " temp",1.0d0/(1.987d-3*beta_h(ib,iparm))," etot",etot, + & " entfac",entfac(i)," ecorr",etot-entfac(i)/beta_h(ib,iparm) etot=etot-entfac(i)/beta_h(ib,iparm) +#endif if(etot.lt.potEmin_all(ib,iparm)) potEmin_all(ib,iparm)=etot -c write (iout,*) "efree",etot," potEmin",potEmin_all(ib,iparm) +#ifdef DEBUG + write (iout,*) "efree",etot," potEmin",potEmin_all(ib,iparm) +#endif enddo ! ib enddo ! iparm enddo ! i @@ -674,7 +782,6 @@ C Determine the minimum energes for all parameter sets and temperatures enddo write (iout,*) "potEmin_min",potEmin_min #endif -#undef DEBUG #ifdef MPI do t=0,tmax @@ -847,6 +954,13 @@ c write (iout,*) "me1",me1," scount",scount(me1) estr=enetb(18,t,iparm) esccor=enetb(19,t,iparm) edihcnstr=enetb(20,t,iparm) + ehomology_constr=enetb(22,t,iparm) + if (homol_nset.gt.1) + & ehomology_constr=waga_homology(ihset)*ehomology_constr + edfadis=enetb(23,t,iparm) + edfator=enetb(24,t,iparm) + edfanei=enetb(25,t,iparm) + edfabet=enetb(26,t,iparm) do k=0,nGridT betaT=startGridT+k*delta_T temper=betaT @@ -950,12 +1064,12 @@ c write (iout,*) ib," PotEmin",potEmin 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*eturn6+ft(2)*wel_loc*eel_loc & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor - & +wbond*estr + & +wbond*estr+ehomology_constr eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees & +ftprim(1)*wtor*etors+ & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+ @@ -973,12 +1087,12 @@ c write (iout,*) ib," PotEmin",potEmin 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*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor - & +wbond*estr + & +wbond*estr+ehomology_constr eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1) & +ftprim(1)*wtor*etors+ & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+ @@ -995,6 +1109,11 @@ c write (iout,*) ib," PotEmin",potEmin #endif weight=dexp(-betaT*(etot-potEmin)+entfac(t)) #ifdef DEBUG + write (iout,'(3i5,6f5.2,17f12.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, + & ehomology_constr+wdfa_dist*edfadis+wdfa_tor*edfator+ + & wdfa_nei*edfanei+wdfa_beta*edfabet write (iout,*) "iparm",iparm," t",t," temper",temper, & " etot",etot," entfac",entfac(t), & " efree",etot-entfac(t)/betaT," potEmin",potEmin,