subroutine make_ensembles(islice,*) ! construct the conformational ensembles at REMD temperatures implicit none include "DIMENSIONS" include "DIMENSIONS.ZSCOPT" include "DIMENSIONS.FREE" #ifdef MPI include "mpif.h" include "COMMON.MPI" integer ierror,errcode,status(MPI_STATUS_SIZE) #endif include "COMMON.IOUNITS" include "COMMON.CONTROL" include "COMMON.HOMOLOGY" include "COMMON.FREE" include "COMMON.ENERGIES" include "COMMON.FFIELD" include "COMMON.INTERACT" include "COMMON.SBRIDGE" include "COMMON.CHAIN" include "COMMON.PROTFILES" include "COMMON.PROT" 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,evdw_t,evdw2,ees,evdw1,ebe,etors, & escloc,eliptran,esaxs, & ehpb,ecorr,ecorr5,ecorr6,eello_turn4,eello_turn3, & eturn6,eel_loc,edihcnstr,etors_d,estr,evdw2_14,esccor,tt, & ehomology_constr,edfadis,edfator,edfanei,edfabet integer i,ii,ik,iproc,iscor,j,k,l,ib,iparm,iprot,nlist double precision qfree,sumprob,eini,efree,rmsdev character*80 bxname character*2 licz1,licz2 character*3 licz3,licz4 character*5 ctemper integer ilen external ilen real*4 Fdimless(MaxStr),Fdimless_(MaxStr) double precision enepot(MaxStr) integer iperm(MaxStr) integer islice #ifdef MPI if (me.eq.Master) then #endif write (licz2,'(bz,i2.2)') islice if (nslice.eq.1) then if (.not.separate_parset) then bxname = prefix(:ilen(prefix))//".bx" else write (licz3,'(bz,i3.3)') myparm bxname = prefix(:ilen(prefix))//"_par"//licz3//".bx" endif else if (.not.separate_parset) then bxname = prefix(:ilen(prefix))//"_slice_"//licz2//".bx" else write (licz3,'(bz,i3.3)') myparm bxname = prefix(:ilen(prefix))//"par_"//licz3// & "_slice_"//licz2//".bx" endif endif open (ientout,file=bxname,status="unknown", & form="unformatted",access="direct",recl=lenrec1) #ifdef MPI endif #endif do iparm=1,nParmSet if (iparm.ne.iparmprint) exit call restore_parm(iparm) do ib=1,nT_h(iparm) #ifdef DEBUG write (iout,*) "iparm",iparm," ib",ib #endif temper=1.0d0/(beta_h(ib,iparm)*1.987D-3) c quot=1.0d0/(T0*beta_h(ib,iparm)*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 if (rescale_mode.eq.1) then quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3) #if defined(FUNCTH) tt=1.0d0/(beta_h(ib,iparm)*1.987D-3) ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/320.0d0 #elif defined(FUNCT) ft(6)=quot #else ft(6)=1.0d0 #endif 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 else if (rescale_mode.eq.2) then quot=1.0d0/(T0*beta_h(ib,iparm)*1.987D-3) #if defined(FUNCTH) tt=1.0d0/(beta_h(ib,iparm)*1.987D-3) ft(6)=(320.0d0+80.0d0*dtanh((tt-320.0d0)/80.0d0))/3200.d0 #elif defined(FUNCT) ft(6)=quot #else ft(6)=1.0d0 #endif 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,iparm)*1.987D-3),ft else if (rescale_mode.eq.0) then do l=1,5 fT(l)=0.0d0 enddo else write (iout,*) & "Error in MAKE_ENSEMBLE: Wrong RESCALE_MODE:",rescale_mode call flush(iout) return1 endif #ifdef MPI do i=1,scount(me1) #else do i=1,ntot(islice) #endif evdw=enetb(1,i,iparm) evdw_t=enetb(21,i,iparm) #ifdef SCP14 evdw2_14=enetb(17,i,iparm) evdw2=enetb(2,i,iparm)+evdw2_14 #else evdw2=enetb(2,i,iparm) evdw2_14=0.0d0 #endif #ifdef SPLITELE ees=enetb(3,i,iparm) evdw1=enetb(16,i,iparm) #else ees=enetb(3,i,iparm) evdw1=0.0d0 #endif ecorr=enetb(4,i,iparm) ecorr5=enetb(5,i,iparm) ecorr6=enetb(6,i,iparm) eel_loc=enetb(7,i,iparm) eello_turn3=enetb(8,i,iparm) eello_turn4=enetb(9,i,iparm) eturn6=enetb(10,i,iparm) ebe=enetb(11,i,iparm) escloc=enetb(12,i,iparm) etors=enetb(13,i,iparm) etors_d=enetb(14,i,iparm) ehpb=enetb(15,i,iparm) estr=enetb(18,i,iparm) esccor=enetb(19,i,iparm) edihcnstr=enetb(20,i,iparm) eliptran=enetb(22,i,iparm) esaxs=enetb(26,i,iparm) ehomology_constr=enetb(27,i,iparm) edfadis=enetb(28,i,iparm) edfator=enetb(29,i,iparm) edfanei=enetb(30,i,iparm) edfabet=enetb(31,i,iparm) if (homol_nset.gt.1) & ehomology_constr=waga_homology(homol_nset)*ehomology_constr #ifdef SPLITELE if (shield_mode.gt.0) then etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*wscp*evdw2 & +ft(1)*welec*ees & +ft(1)*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+wliptran*eliptran+wsaxs*esaxs & +ehomology_constr & +wdfa_dist*edfadis & +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet else 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+wliptran*eliptran+wsaxs*esaxs & +ehomology_constr & +wdfa_dist*edfadis & +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet endif #else if (shield_mode.gt.0) then etot=ft(1)*wsc*(evdw+ft(6)*evdw_t)+ft(1)*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+wliptran*eliptran+wsaxs*esaxs & +ehomology_constr & +wdfa_dist*edfadis & +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet 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+wliptran*eliptran+wsaxs*esaxs & +ehomology_constr & +wdfa_dist*edfadis & +wdfa_tor*edfator+wdfa_nei*edfanei+wdfa_beta*edfabet endif #endif #ifdef MPI Fdimless_(i)= & beta_h(ib,iparm)*etot-entfac(i) potE(i,iparm)=etot #ifdef DEBUG write (iout,*) i,indstart(me)+i-1,ib, & 1.0d0/(1.987d-3*beta_h(ib,iparm)),potE(i,iparm), & -entfac(i),Fdimless(i) #endif #else Fdimless(i)=beta_h(ib,iparm)*etot-entfac(i) potE(i,iparm)=etot #endif enddo ! i #ifdef MPI call MPI_Gatherv(Fdimless_(1),scount(me), & MPI_REAL,Fdimless(1), & scount(0),idispl(0),MPI_REAL,Master, & WHAM_COMM, IERROR) #ifdef DEBUG call MPI_Gatherv(potE(1,iparm),scount(me), & MPI_DOUBLE_PRECISION,potE(1,iparm), & scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master, & WHAM_COMM, IERROR) call MPI_Gatherv(entfac(1),scount(me), & MPI_DOUBLE_PRECISION,entfac(1), & scount(0),idispl(0),MPI_DOUBLE_PRECISION,Master, & WHAM_COMM, IERROR) #endif if (me.eq.Master) then #ifdef DEBUG write (iout,*) "The FDIMLESS array before sorting" do i=1,ntot(islice) write (iout,*) i,fdimless(i) enddo #endif #endif do i=1,ntot(islice) iperm(i)=i enddo call mysort1(ntot(islice),Fdimless,iperm) #ifdef DEBUG write (iout,*) "The FDIMLESS array after sorting" do i=1,ntot(islice) write (iout,*) i,iperm(i),fdimless(i) enddo #endif qfree=0.0d0 do i=1,ntot(islice) qfree=qfree+exp(-fdimless(i)+fdimless(1)) enddo c write (iout,*) "qfree",qfree nlist=1 sumprob=0.0 do i=1,min0(ntot(islice),ensembles) sumprob=sumprob+exp(-fdimless(i)+fdimless(1))/qfree #ifdef DEBUG write (iout,*) i,ib,beta_h(ib,iparm), & 1.0d0/(1.987d-3*beta_h(ib,iparm)),iperm(i), & potE(iperm(i),iparm), & -entfac(iperm(i)),fdimless(i),sumprob #endif if (sumprob.gt.0.99d0) goto 122 nlist=nlist+1 enddo 122 continue #ifdef MPI endif call MPI_Bcast(nlist, 1, MPI_INTEGER, Master, WHAM_COMM, & IERROR) call MPI_Bcast(iperm,nlist,MPI_INTEGER,Master,WHAM_COMM, & IERROR) do i=1,nlist ii=iperm(i) iproc=0 do while (ii.lt.indstart(iproc).or.ii.gt.indend(iproc)) iproc=iproc+1 enddo if (iproc.ge.nprocs) then write (iout,*) "Fatal error: processor out of range",iproc call flush(iout) if (bxfile) then close (ientout) else close (ientout,status="delete") endif return1 endif ik=ii-indstart(iproc)+1 if (iproc.ne.Master) then if (me.eq.iproc) then #ifdef DEBUG write (iout,*) "i",i," ii",ii," iproc",iproc," ik",ik, & " energy",potE(ik,iparm) #endif call MPI_Send(potE(ik,iparm),1,MPI_DOUBLE_PRECISION, & Master,i,WHAM_COMM,IERROR) else if (me.eq.Master) then call MPI_Recv(enepot(i),1,MPI_DOUBLE_PRECISION,iproc,i, & WHAM_COMM,STATUS,IERROR) endif else if (me.eq.Master) then enepot(i)=potE(ik,iparm) endif enddo #else do i=1,nlist enepot(i)=potE(iperm(i),iparm) enddo #endif #ifdef MPI if (me.eq.Master) then #endif write(licz3,'(bz,i3.3)') iparm write(licz2,'(bz,i2.2)') islice if (temper.lt.100.0d0) then write(ctemper,'(f3.0)') temper else if (temper.lt.1000.0) then write (ctemper,'(f4.0)') temper else write (ctemper,'(f5.0)') temper endif if (nparmset.eq.1) then if (separate_parset) then write(licz4,'(bz,i3.3)') myparm pdbname=prefix(:ilen(prefix))//"_par"//licz4 else pdbname=prefix(:ilen(prefix)) endif else pdbname=prefix(:ilen(prefix))//"_parm_"//licz3 endif if (nslice.eq.1) then pdbname=pdbname(:ilen(pdbname))//"_T_"// & ctemper(:ilen(ctemper))//"pdb" else pdbname=pdbname(:ilen(pdbname))//"_slice_"//licz2//"_T_"// & ctemper(:ilen(ctemper))//"pdb" endif open(ipdb,file=pdbname) write (iout,*) "Before reading nlist",nlist do i=1,nlist read (ientout,rec=iperm(i)) & ((csingle(l,k),l=1,3),k=1,nres), & ((csingle(l,k+nres),l=1,3),k=nnt,nct), & nss,(ihpb(k),jhpb(k),k=1,nss), & eini,efree,rmsdev,iscor do j=1,2*nres do k=1,3 c(k,j)=csingle(k,j) enddo enddo eini=fdimless(i) call pdbout(iperm(i),temper,eini,enepot(i),efree,rmsdev) enddo #ifdef MPI endif #endif enddo ! ib enddo ! iparm if (bxfile) then close(ientout) else close(ientout,status="delete") 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