--- /dev/null
+ 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
projects / unres4.git / blobdiff