X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=blobdiff_plain;f=source%2Fwham%2Fsrc%2Fwham_calc1.F;fp=source%2Fwham%2Fsrc%2Fwham_calc1.F;h=314644809179e868ebbcdb1026be418d3fb8c4eb;hp=61b95048055fb2208776175fa330b22d67e38111;hb=9453fc761eb545fcb727824c94d012dbf3931951;hpb=6f521277aa2a382d409f5189957283b0998b0d07 diff --git a/source/wham/src/wham_calc1.F b/source/wham/src/wham_calc1.F index 61b9504..3146448 100644 --- a/source/wham/src/wham_calc1.F +++ b/source/wham/src/wham_calc1.F @@ -92,6 +92,7 @@ c parameter (MaxHdim=200000) character*2 licz2 character*3 licz3 character*128 nazwa + character*30 frm_write integer ilen external ilen write(licz2,'(bz,i2.2)') islice @@ -531,8 +532,8 @@ c#endif do i=1,nR(ib,iparm) fi(i,ib,iparm)=fi(i,ib,iparm)-avefi enddo - write (iout,'(8f10.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm)) - write (iout,'(8f10.5)') (f(i,ib,iparm),i=1,nR(ib,iparm)) + write (iout,'(6f15.5)') (fi(i,ib,iparm),i=1,nR(ib,iparm)) + write (iout,'(6f15.5)') (f(i,ib,iparm),i=1,nR(ib,iparm)) enddo enddo @@ -732,7 +733,8 @@ C Determine the minimum energes for all parameter sets and temperatures do i=1,nParmSet write (iout,*) "Parameter set",i do j=1,nT_h(i) - write (iout,*) j,PotEmin_all(j,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 @@ -990,24 +992,38 @@ 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 + potEmin=potEmin_all(1,iparm)+ + & (potEmin_all(1,iparm)-potEmin_all(2,iparm))/ + & (1.0/beta_h(1,iparm)-1.0/beta_h(2,iparm))* + & (1.0/betaT-1.0/beta_h(1,iparm)) +#ifdef DEBUG + write(iout,*) "first",temper,potEmin +#endif + else if (betaT.le.beta_h(nT_h(iparm),iparm)) then + potEmin=potEmin_all(nT_h(iparm),iparm)+ + &(potEmin_all(nT_h(iparm),iparm)-potEmin_all(nT_h(iparm)-1,iparm))/ + &(1.0/beta_h(nT_h(iparm),iparm)-1.0/beta_h(nT_h(iparm)-1,iparm))* + &(1.0/betaT-1.0/beta_h(nt_h(iparm),iparm)) +#ifdef DEBUG + write (iout,*) "last",temper,potEmin +#endif 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 +#ifdef DEBUG + write (iout,*) "l",l, + & betaT,1.0d0/(1.987D-3*beta_h(l,iparm)), + & 1.0d0/(1.987D-3*beta_h(l+1,iparm)),temper,potEmin +#endif exit endif enddo endif -c write (iout,*) ib," PotEmin",potEmin +#ifdef DEBUG + write (iout,*) "k",k," potEmin",potEmin +#endif #ifdef SPLITELE etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*ees & +wvdwpp*evdw1 @@ -1341,6 +1357,9 @@ c write (iout,*) ib," PotEmin",potEmin endif enddo endif + +c write (iout,*) "i",i," potEmin",potEmin + sumE(i,iparm)=sumE(i,iparm)/sumW(i,iparm) sumEbis(i,iparm)=(startGridT+i*delta_T)*sumEbis(i,iparm)/ & sumW(i,iparm) @@ -1363,8 +1382,10 @@ c write (iout,*) ib," PotEmin",potEmin write (iout,*) write (34,'(f7.1,2f15.5,$)') startGridT+i*delta_T, & sumW(i,iparm),sumE(i,iparm) - write (34,'(f10.5,$)') (sumQ(j,i,iparm),j=1,nQ+2) - write (34,'(e15.5,$)') sumEsq(i,iparm)-sumEbis(i,iparm), + write(frm_write,'( "(",i3,"e15.5,$)" )' ) nQ+2 + write (34,frm_write) (sumQ(j,i,iparm),j=1,nQ+2) + write(frm_write,'( "(",i3,"e15.5,$)" )' ) (nQ+2)*2+1 + write (34,frm_write) sumEsq(i,iparm)-sumEbis(i,iparm), & (sumQsq(j,i,iparm),j=1,nQ+2),(sumEQ(j,i,iparm),j=1,nQ+2) write (34,*) enddo