module probability !----------------------------------------------------------------------------- use clust_data implicit none !----------------------------------------------------------------------------- ! ! !----------------------------------------------------------------------------- contains !---------------------------------------------------------------------------- ! probabl.f90 !---------------------------------------------------------------------------- subroutine probabl(ib,nlist,ncon,*) ! construct the conformational ensembles at REMD temperatures ! implicit none ! include "DIMENSIONS" ! include "sizesclu.dat" use io_units use io_clust, only: daread_ccoords use geometry_data, only: nres,c use energy_data, only: nss,ihpb,jhpb,max_ene use geometry, only: int_from_cart1 use energy, only: etotal,rescale_weights use energy, only: rescale_mode,enerprint,weights #ifdef MPI use MPI_data include "mpif.h" ! include "COMMON.MPI" integer :: ierror,errcode !,status(MPI_STATUS_SIZE) #endif ! include "COMMON.IOUNITS" ! include "COMMON.FREE" ! include "COMMON.FFIELD" ! include "COMMON.INTERACT" ! include "COMMON.SBRIDGE" ! include "COMMON.CHAIN" ! include "COMMON.CLUSTER" real(kind=4) :: csingle(3,2*nres) real(kind=8) :: fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl,& eprim,ebis,temper,kfac=2.4d0,T0=300.0d0 real(kind=8) :: 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 integer :: i,ii,ik,iproc,iscor,j,k,l,ib,nlist,ncon real(kind=8) :: qfree,sumprob,eini,efree,rmsdev character(len=80) :: bxname character(len=2) :: licz1 character(len=5) :: ctemper ! integer ilen ! external ilen real(kind=4) :: Fdimless(maxconf),Fdimless_(maxconf) real(kind=8) :: totfree_(0:maxconf),entfac_(maxconf) real(kind=8) :: energia(0:max_ene) integer,dimension(0:nprocs) :: scount_ do i=1,ncon list_conf(i)=i enddo ! do i=1,ncon ! write (iout,*) i,list_conf(i) ! enddo #ifdef MPI write (iout,*) me," indstart",indstart(me)," indend",indend(me) call daread_ccoords(indstart(me),indend(me)) #endif ! write (iout,*) "ncon",ncon temper=1.0d0/(beta_h(ib)*1.987D-3) !elwrite(iout,*)"rescale_mode", rescale_mode ! write (iout,*) "ib",ib," beta_h",beta_h(ib)," temper",temper ! 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 ! EL start old rescale------------------------------- ! 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 ! 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 !el write (iout,*) 1.0d0/(beta_h(ib)*1.987D-3) !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 !EL end old rescele---------------------- ! call rescale_weights(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 ! write (iout,*) "i",i," ii",ii ! 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) enddo enddo do k=1,3 c(k,nres+1)=c(k,1) c(k,nres+nres)=c(k,nres) enddo !el do j=1,2*nres ! do k=1,3 !write(iout,*)"c, k, j",k,j,c(k,j) ! enddo !el enddo 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.) ! call etotal(energia(0),fT) call etotal(energia) totfree(i)=energia(0) ! write (iout,'(8f10.5)') ((c(l,k),l=1,3),k=1,nres) ! write (iout,'(8f10.5)') ((c(l,k+nres),l=1,3),k=nnt,nct) ! call enerprint(energia) ! call pdbout(totfree(i),16,i) #ifdef DEBUG write (iout,*) i," energia",(energia(j),j=0,21) write (iout,*) "etot", etot ! write (iout,*) "ft(6)", ft(6) #endif do k=1,max_ene enetb(k,i)=energia(k) enddo endif !el evdw=enetb(1,i) ! write (iout,*) evdw etot=energia(0) #ifdef SCP14 !el evdw2_14=enetb(17,i) evdw2_14=enetb(18,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) estr=enetb(17,i) ! esccor=enetb(19,i) esccor=enetb(21,i) ! edihcnstr=enetb(20,i) edihcnstr=enetb(19,i) !#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 !#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) totfree(i)=etot #ifdef DEBUG write (iout,*) "fdim calc", i,ii,ib,& 1.0d0/(1.987d-3*beta_h(ib)),totfree(i),& entfac(ii),Fdimless(i) #endif Fdimless_(i)=Fdimless(i) totfree_(i)=totfree(i) call enerprint(energia(0)) !el enddo ! i do i=1,maxconf entfac_(i)=entfac(i) enddo do i=0,nprocs scount_(i)=scount(i) enddo #ifdef MPI call MPI_Gatherv(Fdimless_(1),scount_(me),& MPI_REAL,Fdimless(1),& scount(0),idispl(0),MPI_REAL,Master,& MPI_COMM_WORLD, IERROR) call MPI_Gatherv(totfree_(1),scount(me),& MPI_DOUBLE_PRECISION,totfree(1),& scount_(0),idispl(0),MPI_DOUBLE_PRECISION,Master,& MPI_COMM_WORLD, IERROR) call MPI_Gatherv(entfac_(indstart(me)+1),scount_(me),& MPI_DOUBLE_PRECISION,entfac(1),& scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master,& MPI_COMM_WORLD, IERROR) if (me.eq.Master) then #endif #ifdef DEBUG write (iout,*) "The FDIMLESS array before sorting" do i=1,ncon write (iout,*) i,fdimless(i) enddo #endif call mysort1(ncon,Fdimless,list_conf) #ifdef DEBUG write (iout,*) "The FDIMLESS array after sorting" do i=1,ncon write (iout,*) i,list_conf(i),fdimless(i) enddo #endif do i=1,ncon totfree(i)=fdimless(i) enddo qfree=0.0d0 do i=1,ncon qfree=qfree+exp(-fdimless(i)+fdimless(1)) ! write (iout,*) "fdimless", fdimless(i) enddo write (iout,*) "qfree",qfree !d nlist=1 sumprob=0.0 write (iout,*) "ncon", ncon,maxstr_proc do i=1,min0(ncon,maxstr_proc)-1 sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree #ifdef DEBUG write (iout,*) i,ib,beta_h(ib),& 1.0d0/(1.987d-3*beta_h(ib)),list_conf(i),& totfree(list_conf(i)),& -entfac(list_conf(i)),fdimless(i),sumprob #endif if (sumprob.gt.prob_limit) goto 122 ! if (sumprob.gt.1.00d0) goto 122 nlist=nlist+1 enddo 122 continue #ifdef MPI endif call MPI_Bcast(nlist, 1, MPI_INTEGER, Master, MPI_COMM_WORLD, & IERROR) call MPI_Bcast(list_conf,nlist,MPI_INTEGER,Master,MPI_COMM_WORLD,& IERROR) ! do iproc=0,nprocs ! write (iout,*) "iproc",iproc," indstart",indstart(iproc), ! & " indend",indend(iproc) ! enddo #endif !write(iout,*)"koniec probabl" return end subroutine probabl !-------------------------------------------------- subroutine mysort1(n, x, ipermut) ! implicit none integer :: i,j,imax,ipm,n real(kind=4) :: x(n) integer :: ipermut(n) real(kind=4) :: xtemp do i=1,n xtemp=x(i) imax=i do j=i+1,n if (x(j).lt.xtemp) then imax=j xtemp=x(j) endif enddo x(imax)=x(i) x(i)=xtemp ipm=ipermut(imax) ipermut(imax)=ipermut(i) ipermut(i)=ipm enddo return end subroutine mysort1 !----------------------------------------------------------------------------- !----------------------------------------------------------------------------- end module probability