subroutine probabl(ibb,nlist,ncon,*) ! construct the conformational ensembles at REMD temperatures implicit none include "DIMENSIONS" include "sizesclu.dat" #ifdef MPI include "mpif.h" include "COMMON.MPI" integer ierror,errcode,status(MPI_STATUS_SIZE) #endif include "COMMON.CONTROL" include "COMMON.IOUNITS" include "COMMON.FREE" include "COMMON.FFIELD" include "COMMON.INTERACT" include "COMMON.SBRIDGE" include "COMMON.CHAIN" include "COMMON.CLUSTER" real*4 csingle(3,maxres2) double precision fT(6),fTprim(6),fTbis(6),quot,quotl1,quotl,kfacl, & eprim,ebis,temper,kfac/2.4d0/,T0/300.0d0/ 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,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 integer ilen,ijk external ilen 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 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) #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 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 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) #ifdef MPI do i=1,scount(me) ii=i+indstart(me)-1 #else do i=1,ncon ii=i #endif 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 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) #endif 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 #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 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 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 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 #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=energia(0) c#endif #ifdef DEBUG 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) #endif 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) 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) 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) c WRITE (iout,*) "Wychodze z call MPI_Gatherv (Propabl)" if (me.eq.Master) then c WRITE (iout,*) "me.eq.Master" #endif #ifdef DEBUG write (iout,*) "The FDIMLESS array before sorting" do i=1,ncon c write (iout,*) i,fdimless(i) enddo #endif c WRITE (iout,*) "Wchodze do call mysort1" call mysort1(ncon,Fdimless,list_conf) c WRITE (iout,*) "Wychodze z call mysort1" #ifdef DEBUG write (iout,*) "The FDIMLESS array after sorting" do i=1,ncon write (iout,*) i,list_conf(i),fdimless(i) enddo #endif c WRITE (iout,*) "Wchodze do petli i=1,ncon totfree(i)=fdimless(i)" do i=1,ncon totfree(i)=fdimless(i) enddo qfree=0.0d0 do i=1,ncon qfree=qfree+exp(-fdimless(i)+fdimless(1)) c write (iout,*) "fdimless", fdimless(i) enddo c write (iout,*) "qfree",qfree 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,*) "tu szukaj ponizej 7" 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 c 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) c do iproc=0,nprocs c write (iout,*) "iproc",iproc," indstart",indstart(iproc), c & " indend",indend(iproc) c enddo write (iout,*) "nlist",nlist #endif 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 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 !-------------------------------------------------- subroutine mysort1(n, x, ipermut) implicit none integer i,j,imax,ipm,n real x(n) integer ipermut(n) real 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