From: Cezary Czaplewski Date: Sat, 5 Sep 2020 09:03:58 +0000 (+0200) Subject: cluster average over multiple temperatures X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?p=unres.git;a=commitdiff_plain;h=ae8d9c0b861e24dfda2917c3a7f495eed21c227d cluster average over multiple temperatures --- diff --git a/source/cluster/wham/src-M/COMMON.CONTROL b/source/cluster/wham/src-M/COMMON.CONTROL index bb11a5e..6894685 100644 --- a/source/cluster/wham/src-M/COMMON.CONTROL +++ b/source/cluster/wham/src-M/COMMON.CONTROL @@ -13,14 +13,15 @@ & with_dihed_constr,with_theta_constr,out1file, & print_homology_restraints, & print_contact_map,print_homology_models, - & read2sigma,l_homo,read_homol_frag,dist1cut + & read2sigma,l_homo,read_homol_frag,dist1cut,cumul_prob common /cntrl/ betaT,iscode,indpdb,refstr,pdbref,outpdb,outmol2, & punch_dist,print_dist,caonly,lside,lprint_cart,lprint_int, & from_cart,from_bx,from_cx, with_dihed_constr,with_theta_constr, & efree,iopt,nstart,nend,symetr, & tor_mode,shield_mode, & constr_dist,out1file, - & constr_homology,homol_nset,read2sigma,read_homol_frag,dist1cut + & constr_homology,homol_nset,read2sigma,read_homol_frag,dist1cut, + & cumul_prob common /homol/ waga_homology(10), & waga_dist,waga_angle,waga_theta,waga_d,dist_cut,dist2_cut, & iset,ihset,l_homo(max_template,maxdim), diff --git a/source/cluster/wham/src-M/main_clust.F b/source/cluster/wham/src-M/main_clust.F index 3c7e430..5b3cfe8 100644 --- a/source/cluster/wham/src-M/main_clust.F +++ b/source/cluster/wham/src-M/main_clust.F @@ -41,6 +41,7 @@ C double precision hrtime,mintime,sectime logical eof external ilen + integer nTend #ifdef MPI call MPI_Init( IERROR ) call MPI_Comm_rank( MPI_COMM_WORLD, me, IERROR ) @@ -123,7 +124,12 @@ c call flush(iout) write (iout,*) 'from read_coords: ncon',ncon write (iout,*) "nT",nT - do iT=1,nT + if (cumul_prob) then + nTend=1 + else + nTend=nT + endif + do iT=1,nTend write (iout,*) "iT",iT #ifdef MPI call work_partition(.true.,ncon) diff --git a/source/cluster/wham/src-M/probabl.F b/source/cluster/wham/src-M/probabl.F index 22db88f..198c5f5 100644 --- a/source/cluster/wham/src-M/probabl.F +++ b/source/cluster/wham/src-M/probabl.F @@ -1,4 +1,4 @@ - subroutine probabl(ib,nlist,ncon,*) + subroutine probabl(ibb,nlist,ncon,*) ! construct the conformational ensembles at REMD temperatures implicit none include "DIMENSIONS" @@ -8,6 +8,7 @@ include "COMMON.MPI" integer ierror,errcode,status(MPI_STATUS_SIZE) #endif + include "COMMON.CONTROL" include "COMMON.IOUNITS" include "COMMON.FREE" include "COMMON.FFIELD" @@ -21,9 +22,11 @@ double precision etot,evdw,evdw2,ees,evdw1,ebe,etors,escloc, & ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3, & eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor, - & evdw_t,esaxs - integer i,ii,ik,iproc,iscor,j,k,l,ib,nlist,ncon - double precision qfree,sumprob,eini,efree,rmsdev + & evdw_t,esaxs,eliptran,ehomology_constr + integer i,ii,ik,iproc,iscor,j,k,l,ib,ibb,ib_start,ib_end,nlist, + & ncon + double precision qfree,sumprob,eini,rmsdev + real*4 probab(maxconf),totprob(maxconf),aux character*80 bxname character*2 licz1 character*5 ctemper @@ -32,82 +35,86 @@ real*4 Fdimless(maxconf), Fdimless_buf(maxconf) double precision energia(0:max_ene), totfree_buf(0:maxconf), & entfac_buf(maxconf) + if (cumul_prob) then + write (iout,*) + & "Probabilities will be summed over all temperatures" + ib_start=1 + ib_end=nT + else + ib_start=ibb + ib_end=ibb + endif + probab=0.0 + totprob=0.0 + Fdimless=0.0d0 + Fdimless_buf=0.0d0 + totfree=0.0d0 + totfree_buf=0.0d0 do i=1,ncon list_conf(i)=i enddo -c do i=1,ncon -c write (iout,*) i,list_conf(i) -c enddo -#ifdef MPI - write (iout,*) me," indstart",indstart(me)," indend",indend(me) - call daread_ccoords(indstart(me),indend(me)) -#endif -C write (iout,*) "ncon",ncon -C call flush(iout) - temper=1.0d0/(beta_h(ib)*1.987D-3) -c write (iout,*) "ib",ib," beta_h",beta_h(ib)," temper",temper -c quot=1.0d0/(T0*beta_h(ib)*1.987D-3) -c quotl=1.0d0 -c kfacl=1.0d0 -c do l=1,5 -c quotl1=quotl -c quotl=quotl*quot -c kfacl=kfacl*kfac -c fT(l)=kfacl/(kfacl-1.0d0+quotl) -c enddo -C#define DEBUG - if (rescale_mode.eq.1) then - quot=1.0d0/(T0*beta_h(ib)*1.987D-3) - 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 + ib=1 + if (rescale_mode.eq.1) then + quot=1.0d0/(T0*beta_h(ib)*1.987D-3) + 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) - 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) + 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 + 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 + fT(6)=1.0d0 + ftprim(6)=0.0d0 + ftbis(6)=0.0d0 #endif - else if (rescale_mode.eq.2) then - quot=1.0d0/(T0*beta_h(ib)*1.987D-3) - quotl=1.0d0 - do l=1,5 - quotl=quotl*quot - fT(l)=1.12692801104297249644d0/ - & dlog(dexp(quotl)+dexp(-quotl)) - enddo -c write (iout,*) 1.0d0/(beta_h(ib)*1.987D-3),ft -c call flush(iout) + else if (rescale_mode.eq.2) then + quot=1.0d0/(T0*beta_h(ib)*1.987D-3) + quotl=1.0d0 + do l=1,5 + quotl=quotl*quot + fT(l)=1.12692801104297249644d0/ + & dlog(dexp(quotl)+dexp(-quotl)) + enddo +c write (iout,*) 1.0d0/(beta_h(ib)*1.987D-3),ft +c call flush(iout) #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) + 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 + 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 + fT(6)=1.0d0 + ftprim(6)=0.0d0 + ftbis(6)=0.0d0 +#endif + endif +c do i=1,ncon +c write (iout,*) i,list_conf(i) +c enddo +#ifdef MPI + write (iout,*) me," indstart",indstart(me)," indend",indend(me) + call daread_ccoords(indstart(me),indend(me)) #endif - endif +C write (iout,*) "ncon",ncon +C call flush(iout) #ifdef MPI do i=1,scount(me) @@ -116,151 +123,263 @@ c call flush(iout) do i=1,ncon ii=i #endif -C write (iout,*) "i",i," ii",ii,"ib",ib,scount(me) -c call flush(iout) - if (ib.eq.1) then - do j=1,nres - do k=1,3 - c(k,j)=allcart(k,j,i) - c(k,j+nres)=allcart(k,j+nres,i) -C write(iout,*) "coord",i,j,k,allcart(k,j,i),c(k,j), -C & c(k,j+nres),allcart(k,j+nres,i) - enddo - enddo -C write(iout,*) "out of j loop" -C call flush(iout) + do j=1,nres do k=1,3 - c(k,nres+1)=c(k,1) - c(k,nres+nres)=c(k,nres) + c(k,j)=allcart(k,j,i) + c(k,j+nres)=allcart(k,j+nres,i) +C write(iout,*) "coord",i,j,k,allcart(k,j,i),c(k,j), +C & c(k,j+nres),allcart(k,j+nres,i) enddo -C write(iout,*) "after nres+nres",nss_all(i) -C call flush(iout) - nss=nss_all(i) - do j=1,nss - ihpb(j)=ihpb_all(j,i) - jhpb(j)=jhpb_all(j,i) - enddo - call int_from_cart1(.false.) -C write(iout,*) "before etotal" -C call flush(iout) - call etotal(energia(0),fT) - totfree(i)=energia(0) - totfree_buf(i)=totfree(i) + enddo +C write(iout,*) "out of j loop" +C call flush(iout) + do k=1,3 + c(k,nres+1)=c(k,1) + c(k,nres+nres)=c(k,nres) + enddo +C write(iout,*) "after nres+nres",nss_all(i) +C call flush(iout) + nss=nss_all(i) + do j=1,nss + ihpb(j)=ihpb_all(j,i) + jhpb(j)=jhpb_all(j,i) + enddo + call int_from_cart1(.false.) +C write(iout,*) "before etotal" +C call flush(iout) + call etotal(energia(0),fT) c write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres) c write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct) c call enerprint(energia(0),fT) c call pdbout(totfree(i),16,i) #ifdef DEBUG - write (iout,*) i," energia",(energia(j),j=0,max_ene) - write (iout,*) "etot", etot - write (iout,*) "ft(6)", ft(6) + write (iout,*) i," energia",(energia(j),j=0,max_ene) + write (iout,*) "etot", etot + write (iout,*) "ft(6)", ft(6) #endif - do k=1,max_ene - enetb(k,i)=energia(k) - enddo - endif - evdw=enetb(1,i) -c write (iout,*) evdw - etot=energia(0) + do k=1,max_ene + enetb(k,i)=energia(k) + enddo + enddo ! i +c#endif +C write (iout,*) "i",i," ii",ii,"ib",ib,scount(me) +c call flush(iout) + do ib=ib_start,ib_end + temper=1.0d0/(beta_h(ib)*1.987D-3) +#ifdef MPI + do i=1,scount(me) + ii=i+indstart(me)-1 +#else + do i=1,ncon + ii=i +#endif + evdw=enetb(1,i) +c write (iout,*) evdw + etot=energia(0) #ifdef SCP14 - evdw2_14=enetb(17,i) - evdw2=enetb(2,i)+evdw2_14 + evdw2_14=enetb(17,i) + evdw2=enetb(2,i)+evdw2_14 +#else + evdw2=enetb(2,i) + evdw2_14=0.0d0 +#endif +#ifdef SPLITELE + ees=enetb(3,i) + evdw1=enetb(16,i) +#else + ees=enetb(3,i) + evdw1=0.0d0 +#endif + ecorr=enetb(4,i) + ecorr5=enetb(5,i) + ecorr6=enetb(6,i) + eel_loc=enetb(7,i) + eello_turn3=enetb(8,i) + eello_turn4=enetb(9,i) + eturn6=enetb(10,i) + ebe=enetb(11,i) + escloc=enetb(12,i) + etors=enetb(13,i) + etors_d=enetb(14,i) + ehpb=enetb(15,i) + estr=enetb(18,i) + esccor=enetb(19,i) + edihcnstr=enetb(20,i) + evdw_t=enetb(21,i) + eliptran=enetb(22,i) + ehomology_constr=enetb(24,i) + esaxs=enetb(25,i) + if (rescale_mode.eq.1) then + quot=1.0d0/(T0*beta_h(ib)*1.987D-3) + 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) + 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 - evdw2=enetb(2,i) - evdw2_14=0.0d0 + fT(6)=1.0d0 + ftprim(6)=0.0d0 + ftbis(6)=0.0d0 #endif + + else if (rescale_mode.eq.2) then + quot=1.0d0/(T0*beta_h(ib)*1.987D-3) + quotl=1.0d0 + do l=1,5 + quotl=quotl*quot + fT(l)=1.12692801104297249644d0/ + & dlog(dexp(quotl)+dexp(-quotl)) + enddo +c write (iout,*) 1.0d0/(beta_h(ib)*1.987D-3),ft +c call flush(iout) +#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 + endif +c write (iout,*) 1.0d0/(beta_h(ib)*1.987D-3),ft +c call flush(iout) +c call etotal(energia(0),fT) +c#ifdef CHUJ #ifdef SPLITELE - ees=enetb(3,i) - evdw1=enetb(16,i) + if (shield_mode.gt.0) then + etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2 + & +welec*ft(1)*ees + & +ft(1)*wvdwpp*evdw1 + & +wang*ebe+wtor*ft(1)*etors+wscloc*escloc + & +wstrain*ehpb+wcorr*ft(3)*ecorr+wcorr5*ft(4)*ecorr5 + & +wcorr6*ft(5)*ecorr6+wturn4*ft(3)*eello_turn4 + & +wturn3*ft(2)*eello_turn3+wturn6*ft(5)*eturn6 + & +wel_loc*ft(2)*eel_loc+edihcnstr+wtor_d*ft(2)*etors_d + & +wbond*estr+wsccor*ft(1)*esccor+!ethetacnstr + & +wliptran*eliptran+wsaxs*esaxs + else + etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+welec*ft(1)*ees + & +wvdwpp*evdw1 + & +wang*ebe+wtor*ft(1)*etors+wscloc*escloc + & +wstrain*ehpb+wcorr*ft(3)*ecorr+wcorr5*ft(4)*ecorr5 + & +wcorr6*ft(5)*ecorr6+wturn4*ft(3)*eello_turn4 + & +wturn3*ft(2)*eello_turn3+wturn6*ft(5)*eturn6 + & +wel_loc*ft(2)*eel_loc+edihcnstr+wtor_d*ft(2)*etors_d + & +wbond*estr+wsccor*ft(1)*esccor+ehomology_constr + & +wliptran*eliptran+wsaxs*esaxs + endif #else - ees=enetb(3,i) - evdw1=0.0d0 + if (shield_mode.gt.0) then + etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2 + & +welec*ft(1)*(ees+evdw1) + & +wang*ebe+wtor*ft(1)*etors+wscloc*escloc + & +wstrain*ehpb+wcorr*ft(3)*ecorr+wcorr5*ft(4)*ecorr5 + & +wcorr6*ft(5)*ecorr6+wturn4*ft(3)*eello_turn4 + & +wturn3*ft(2)*eello_turn3+wturn6*ft(5)*eturn6 + & +wel_loc*ft(2)*eel_loc+edihcnstr+wtor_d*ft(2)*etors_d + & +wbond*estr+wsccor*ft(1)*esccor+ehomology_constr + & +wliptran*eliptran+wsaxs*esaxs + else + etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2 + & +welec*ft(1)*(ees+evdw1) + & +wang*ebe+wtor*ft(1)*etors+wscloc*escloc + & +wstrain*ehpb+wcorr*ft(3)*ecorr+wcorr5*ft(4)*ecorr5 + & +wcorr6*ft(5)*ecorr6+wturn4*ft(3)*eello_turn4 + & +wturn3*ft(2)*eello_turn3+wturn6*ft(5)*eturn6 + & +wel_loc*ft(2)*eel_loc+edihcnstr+wtor_d*ft(2)*etors_d + & +wbond*estr+wsccor*ft(1)*esccor+!ethetacnstr + & +wliptran*eliptran+wsaxs*esaxs + endif #endif - ecorr=enetb(4,i) - ecorr5=enetb(5,i) - ecorr6=enetb(6,i) - eel_loc=enetb(7,i) - eello_turn3=enetb(8,i) - eello_turn4=enetb(9,i) - eturn6=enetb(10,i) - ebe=enetb(11,i) - escloc=enetb(12,i) - etors=enetb(13,i) - etors_d=enetb(14,i) - ehpb=enetb(15,i) - estr=enetb(18,i) - esccor=enetb(19,i) - edihcnstr=enetb(20,i) - esaxs=enetb(25,i) -c#ifdef SPLITELE -c etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ -c &ft(1)*welec*ees+wvdwpp*evdw1 -c & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc -c & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 -c & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 -c & +ft(2)*wturn3*eello_turn3 -c & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc -c & +edihcnstr+ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor -c & +wbond*estr +#ifdef DEBUG + write (iout,*) "etot2", etot + write (iout,*) "evdw", wsc, evdw,evdw_t + write (iout,*) "evdw2", wscp, evdw2 + write (iout,*) "welec", ft(1),welec,ees + write (iout,*) "evdw1", wvdwpp,evdw1 + write (iout,*) "ebe", ebe,wang +#endif c#else -c etot=wsc*(evdw+ft(6)*evdw_t)+wscp*evdw2+ft(1)*welec*(ees+evdw1) -c & +wang*ebe+ft(1)*wtor*etors+wscloc*escloc -c & +wstrain*ehpb+nss*ebr+ft(3)*wcorr*ecorr+ft(4)*wcorr5*ecorr5 -c & +ft(5)*wcorr6*ecorr6+ft(3)*wturn4*eello_turn4 -c & +ft(2)*wturn3*eello_turn3 -c & +ft(5)*wturn6*eturn6+ft(2)*wel_loc*eel_loc+edihcnstr -c & +ft(2)*wtor_d*etors_d+ft(1)*wsccor*esccor -c & +wbond*estr +c etot=energia(0) c#endif #ifdef DEBUG - write (iout,*) "etot2", etot - write (iout,*) "evdw", wsc, evdw,evdw_t - write (iout,*) "evdw2", wscp, evdw2 - write (iout,*) "welec", ft(1),welec,ees - write (iout,*) "evdw1", wvdwpp,evdw1 - write (iout,*) "ebe", ebe,wang -#endif - Fdimless(i)=beta_h(ib)*etot+entfac(ii) - Fdimless_buf(i)=Fdimless(i) - totfree(i)=etot - totfree_buf(i)=totfree(i) + write (iout,'(2i5,30f10.2)') i,ii,etot, + & (enetb(ijk,i),ijk=1,max_ene) +#endif + Fdimless(i)=beta_h(ib)*etot+entfac(ii) #ifdef DEBUG - write (iout,*) "fdim calc", i,ii,ib, - & 1.0d0/(1.987d-3*beta_h(ib)),totfree(i), - & entfac(ii),Fdimless(i) + write (iout,*) "fdim calc", i,ii,ib, + & 1.0d0/(1.987d-3*beta_h(ib)),totfree(i), + & entfac(ii),Fdimless(i) #endif - enddo ! i - - do ijk=1,maxconf - entfac_buf(ijk)=entfac(ijk) - Fdimless_buf(ijk)=Fdimless(ijk) - enddo - do ijk=0,maxconf - totfree_buf(ijk)=totfree(ijk) - enddo - - -c scount_buf=scount(me) -c scount_buf2=scount(0) - -c entfac_buf(indstart(me)+1)=entfac(indstart(me)+1) + Fdimless_buf(i)=Fdimless(i) + totfree_buf(i)=totfree(i) + enddo ! i #ifdef MPI c WRITE (iout,*) "Wchodze do call MPI_Gatherv1 (Propabl)" - call MPI_Gatherv(Fdimless_buf(1),scount(me), - & MPI_REAL,Fdimless(1), - & scount(0),idispl(0),MPI_REAL,Master, - & MPI_COMM_WORLD, IERROR) + call MPI_Gatherv(Fdimless_buf(1),scount(me), + & MPI_REAL,Fdimless(1), + & scount(0),idispl(0),MPI_REAL,Master, + & MPI_COMM_WORLD, IERROR) c WRITE (iout,*) "Wchodze do call MPI_Gatherv2 (Propabl)" - call MPI_Gatherv(totfree_buf(1),scount(me), - & MPI_DOUBLE_PRECISION,totfree(1), - & scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master, - & MPI_COMM_WORLD, IERROR) + call MPI_Gatherv(totfree_buf(1),scount(me), + & MPI_DOUBLE_PRECISION,totfree(1), + & scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master, + & MPI_COMM_WORLD, IERROR) c WRITE (iout,*) "Wchodze do call MPI_Gatherv3 (Propabl)" +c WRITE (iout,*) "Wychodze z call MPI_Gatherv (Propabl)" + aux=Fdimless(1) + do i=1,ncon + if (Fdimless(i).lt.aux) aux=Fdimless(i) + enddo + sumprob=0.0d0 + do i=1,ncon + probab(i)=exp(aux-fdimless(i)) + sumprob=sumprob+probab(i) + enddo + do i=1,ncon + probab(i)=probab(i)/sumprob + totprob(i)=totprob(i)+probab(i) +c write (iout,*) ib,i,probab(i),totprob(i) + enddo + enddo ! ib + sumprob=0.0d0 + do i=1,ncon + sumprob=sumprob+totprob(i) + enddo + write (iout,*) "sumprob",sumprob + do i=1,ncon + totprob(i)=totprob(i)/sumprob + enddo + do i=1,ncon +c Fdimless(i)=-dlog(Fdimless(i)) + Fdimless(i)=-alog(totprob(i)) + enddo call MPI_Gatherv(entfac_buf(indstart(me)+1),scount(me), - & MPI_DOUBLE_PRECISION,entfac(1), - & scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master, - & MPI_COMM_WORLD, IERROR) + & MPI_DOUBLE_PRECISION,entfac(1), + & scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master, + & MPI_COMM_WORLD, IERROR) c WRITE (iout,*) "Wychodze z call MPI_Gatherv (Propabl)" if (me.eq.Master) then c WRITE (iout,*) "me.eq.Master" @@ -319,9 +438,16 @@ c & " indend",indend(iproc) c enddo write (iout,*) "nlist",nlist #endif - write (iout,'(80(1h-)/i4,a,f6.3,a,f5.1/80(1h-))') + if (cumul_prob) then + write (iout,'(80(1h-)/i4,a,f6.3,a,10f6.1/80(1h-))') & nlist," conformations found within",prob_limit, - & " probablity cut_off at temperature ",1.0d0/(1.987d-3*beta_h(ib)) + & " probablity cut_off in the temperature range ", + & (1.0d0/(1.987d-3*beta_h(ib)),ib=1,nT) + else + write (iout,'(80(1h-)/i4,a,f6.3,a,f6.1/80(1h-))') + & nlist," conformations found within",prob_limit, + &" probablity cut_off at temperature ",1.0d0/(1.987d-3*beta_h(ibb)) + endif return end !-------------------------------------------------- diff --git a/source/cluster/wham/src-M/readrtns.F b/source/cluster/wham/src-M/readrtns.F index 05d6f22..5a42f4b 100644 --- a/source/cluster/wham/src-M/readrtns.F +++ b/source/cluster/wham/src-M/readrtns.F @@ -110,7 +110,9 @@ C long axis of side chain lside = index(controlcard,"SIDE").gt.0 efree = index(controlcard,"EFREE").gt.0 call readi(controlcard,'NTEMP',nT,1) - write (iout,*) "nT",nT + cumul_prob=nt.lt.0 + nT=iabs(nT) + write (iout,*) "nT",nT," cumul_prob",cumul_prob call multreada(controlcard,'TEMPER',beta_h,nT,300.0d0) write (iout,*) "nT",nT write (iout,*) 'beta_h',(beta_h(i),i=1,nT)