X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Fwham%2Fsrc%2Fwham_calc1.F;h=41f42bc41334f38200099d11ad01cfc751e0e007;hb=ff639b60bc8a1c011c8fed81bc1f82e4b35f4e2a;hp=ec46232f3bb8ac1b0162c222e23862ef3b097d3f;hpb=357148eb01c1ea7e5cfd362d34f1e1ab61c8f2a8;p=unres.git diff --git a/source/wham/src/wham_calc1.F b/source/wham/src/wham_calc1.F index ec46232..41f42bc 100644 --- a/source/wham/src/wham_calc1.F +++ b/source/wham/src/wham_calc1.F @@ -313,7 +313,7 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft 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, & ehomology_constr,edfadis,edfator,edfanei,edfabet @@ -328,7 +328,7 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft & +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+ehomology_constr+wdfa_dist*edfadis + & +wbond*estr+wdfa_dist*edfadis & +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet #else etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 @@ -339,7 +339,7 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft & +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+ehomology_constr+wdfa_dist*edfadis + & +wbond*estr+wdfa_dist*edfadis & +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet #endif #ifdef DEBUG @@ -356,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 @@ -370,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 @@ -541,8 +561,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 @@ -674,12 +694,15 @@ 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, & ehomology_constr+wdfa_dist*edfadis+wdfa_tor*edfator+ @@ -709,12 +732,16 @@ c write (iout,*) 1.0d0/(beta_h(ib,iparm)*1.987D-3),ft & +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) +#endif 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) +#ifdef DEBUG + write (iout,*) "efree",etot," potEmin",potEmin_all(ib,iparm) +#endif enddo ! ib enddo ! iparm enddo ! i @@ -750,7 +777,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 @@ -928,6 +956,13 @@ c write (iout,*) "me1",me1," scount",scount(me1) 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) +c do k=0,nGridT do k=0,nGridT betaT=startGridT+k*delta_T temper=betaT @@ -1009,24 +1044,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 @@ -1037,6 +1086,8 @@ c write (iout,*) ib," PotEmin",potEmin & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor & +wbond*estr+ehomology_constr + & +wdfa_dist*edfadis + & +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet eprim=ftprim(6)*evdw_t+ftprim(1)*welec*ees & +ftprim(1)*wtor*etors+ & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+ @@ -1060,6 +1111,8 @@ c write (iout,*) ib," PotEmin",potEmin & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor & +wbond*estr+ehomology_constr + & +wdfa_dist*edfadis + & +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet eprim=ftprim(6)*evdw_t+ftprim(1)*welec*(ees+evdw1) & +ftprim(1)*wtor*etors+ & ftprim(3)*wcorr*ecorr+ftprim(4)*wcorr5*ecorr5+ @@ -1076,6 +1129,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, @@ -1345,6 +1403,7 @@ c write (iout,*) ib," PotEmin",potEmin else write (iout,'(a,i3)') "Parameter set",iparm endif +c do i=0,NGridT do i=0,NGridT betaT=1.0d0/(1.987D-3*(startGridT+i*delta_T)) if (betaT.ge.beta_h(1,iparm)) then @@ -1360,6 +1419,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)