--- /dev/null
+ subroutine entmcm
+C Does modified entropic sampling in the space of minima.
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.IOUNITS'
+#ifdef MPL
+ include 'COMMON.INFO'
+#endif
+ include 'COMMON.GEO'
+ include 'COMMON.CHAIN'
+ include 'COMMON.MCM'
+ include 'COMMON.MCE'
+ include 'COMMON.CONTACTS'
+ include 'COMMON.CONTROL'
+ include 'COMMON.VAR'
+ include 'COMMON.INTERACT'
+ include 'COMMON.THREAD'
+ include 'COMMON.NAMES'
+ logical accepted,not_done,over,ovrtim,error,lprint
+ integer MoveType,nbond
+ integer conf_comp
+ double precision RandOrPert
+ double precision varia(maxvar),elowest,ehighest,eold
+ double precision przes(3),obr(3,3)
+ double precision varold(maxvar)
+ logical non_conv
+ double precision energia(0:n_ene),energia_ave(0:n_ene)
+C
+cd write (iout,*) 'print_mc=',print_mc
+ WhatsUp=0
+ maxtrial_iter=50
+c---------------------------------------------------------------------------
+C Initialize counters.
+c---------------------------------------------------------------------------
+C Total number of generated confs.
+ ngen=0
+C Total number of moves. In general this won't be equal to the number of
+C attempted moves, because we may want to reject some "bad" confs just by
+C overlap check.
+ nmove=0
+C Total number of shift (nbond_move(1)), spike, crankshaft, three-bond,...
+C motions.
+ do i=1,nres
+ nbond_move(i)=0
+ enddo
+C Initialize total and accepted number of moves of various kind.
+ do i=0,MaxMoveType
+ moves(i)=0
+ moves_acc(i)=0
+ enddo
+C Total number of energy evaluations.
+ neneval=0
+ nfun=0
+ indminn=-max_ene
+ indmaxx=max_ene
+ delte=0.5D0
+ facee=1.0D0/(maxacc*delte)
+ conste=dlog(facee)
+C Read entropy from previous simulations.
+ if (ent_read) then
+ read (ientin,*) indminn,indmaxx,emin,emax
+ print *,'indminn=',indminn,' indmaxx=',indmaxx,' emin=',emin,
+ & ' emax=',emax
+ do i=-max_ene,max_ene
+ entropy(i)=(emin+i*delte)*betbol
+ enddo
+ read (ientin,*) (ijunk,ejunk,entropy(i),i=indminn,indmaxx)
+ indmin=indminn
+ indmax=indmaxx
+ write (iout,*) 'indminn=',indminn,' indmaxx=',indmaxx,
+ & ' emin=',emin,' emax=',emax
+ write (iout,'(/a)') 'Initial entropy'
+ do i=indminn,indmaxx
+ write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
+ enddo
+ endif ! ent_read
+C Read the pool of conformations
+ call read_pool
+C----------------------------------------------------------------------------
+C Entropy-sampling simulations with continually updated entropy
+C Loop thru simulations
+C----------------------------------------------------------------------------
+ DO ISWEEP=1,NSWEEP
+C----------------------------------------------------------------------------
+C Take a conformation from the pool
+C----------------------------------------------------------------------------
+ if (npool.gt.0) then
+ ii=iran_num(1,npool)
+ do i=1,nvar
+ varia(i)=xpool(i,ii)
+ enddo
+ write (iout,*) 'Took conformation',ii,' from the pool energy=',
+ & epool(ii)
+ call var_to_geom(nvar,varia)
+C Print internal coordinates of the initial conformation
+ call intout
+ else
+ call gen_rand_conf(1,*20)
+ endif
+C----------------------------------------------------------------------------
+C Compute and print initial energies.
+C----------------------------------------------------------------------------
+ nsave=0
+#ifdef MPL
+ if (MyID.eq.MasterID) then
+ do i=1,nctasks
+ nsave_part(i)=0
+ enddo
+ endif
+#endif
+ Kwita=0
+ WhatsUp=0
+ write (iout,'(/80(1h*)/a,i2/80(1h*)/)') 'MCE iteration #',isweep
+ write (iout,'(/80(1h*)/a)') 'Initial energies:'
+ call chainbuild
+ call etotal(energia(0))
+ etot = energia(0)
+ call enerprint(energia(0))
+C Minimize the energy of the first conformation.
+ if (minim) then
+ call geom_to_var(nvar,varia)
+ call minimize(etot,varia,iretcode,nfun)
+ call etotal(energia(0))
+ etot = energia(0)
+ write (iout,'(/80(1h*)/a/80(1h*))')
+ & 'Results of the first energy minimization:'
+ call enerprint(energia(0))
+ endif
+ if (refstr) then
+ kkk=1
+ call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,kkk),
+ &nsup,przes,
+ & obr,non_conv)
+ rms=dsqrt(rms)
+ call contact(.false.,ncont,icont,co)
+ frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
+ write (iout,'(a,f8.3,a,f8.3,a,f8.3)')
+ & 'RMS deviation from the reference structure:',rms,
+ & ' % of native contacts:',frac*100,' contact order:',co
+ write (istat,'(i5,11(1pe14.5))') 0,
+ & (energia(print_order(i)),i=1,nprint_ene),etot,rms,frac,co
+ else
+ write (istat,'(i5,9(1pe14.5))') 0,
+ & (energia(print_order(i)),i=1,nprint_ene),etot
+ endif
+ close(istat)
+ neneval=neneval+nfun+1
+ if (.not. ent_read) then
+C Initialize the entropy array
+ do i=-max_ene,max_ene
+ emin=etot
+C Uncomment the line below for actual entropic sampling (start with uniform
+C energy distribution).
+c entropy(i)=0.0D0
+C Uncomment the line below for multicanonical sampling (start with Boltzmann
+C distribution).
+ entropy(i)=(emin+i*delte)*betbol
+ enddo
+ emax=10000000.0D0
+ emin=etot
+ write (iout,'(/a)') 'Initial entropy'
+ do i=indminn,indmaxx
+ write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
+ enddo
+ endif ! ent_read
+#ifdef MPL
+ call recv_stop_sig(Kwita)
+ if (whatsup.eq.1) then
+ call send_stop_sig(-2)
+ not_done=.false.
+ else if (whatsup.le.-2) then
+ not_done=.false.
+ else if (whatsup.eq.2) then
+ not_done=.false.
+ else
+ not_done=.true.
+ endif
+#else
+ not_done = (iretcode.ne.11)
+#endif
+ write (iout,'(/80(1h*)/20x,a/80(1h*))')
+ & 'Enter Monte Carlo procedure.'
+ close(igeom)
+ call briefout(0,etot)
+ do i=1,nvar
+ varold(i)=varia(i)
+ enddo
+ eold=etot
+ indeold=(eold-emin)/delte
+ deix=eold-(emin+indeold*delte)
+ dent=entropy(indeold+1)-entropy(indeold)
+cd write (iout,*) 'indeold=',indeold,' deix=',deix,' dent=',dent
+cd write (*,*) 'Processor',MyID,' indeold=',indeold,' deix=',deix,
+cd & ' dent=',dent
+ sold=entropy(indeold)+(dent/delte)*deix
+ elowest=etot
+ write (iout,*) 'eold=',eold,' sold=',sold,' elowest=',etot
+ write (*,*) 'Processor',MyID,' eold=',eold,' sold=',sold,
+ & ' elowest=',etot
+ if (minim) call zapis(varia,etot)
+ nminima(1)=1.0D0
+C NACC is the counter for the accepted conformations of a given processor
+ nacc=0
+C NACC_TOT counts the total number of accepted conformations
+ nacc_tot=0
+#ifdef MPL
+ if (MyID.eq.MasterID) then
+ call receive_MCM_info
+ else
+ call send_MCM_info(2)
+ endif
+#endif
+ do iene=indminn,indmaxx
+ nhist(iene)=0.0D0
+ enddo
+ do i=2,maxsave
+ nminima(i)=0.0D0
+ enddo
+C Main loop.
+c----------------------------------------------------------------------------
+ elowest=1.0D10
+ ehighest=-1.0D10
+ it=0
+ do while (not_done)
+ it=it+1
+ if (print_mc.gt.0) write (iout,'(80(1h*)/20x,a,i7)')
+ & 'Beginning iteration #',it
+C Initialize local counter.
+ ntrial=0 ! # of generated non-overlapping confs.
+ noverlap=0 ! # of overlapping confs.
+ accepted=.false.
+ do while (.not. accepted .and. WhatsUp.eq.0 .and. Kwita.eq.0)
+ ntrial=ntrial+1
+C Retrieve the angles of previously accepted conformation
+ do j=1,nvar
+ varia(j)=varold(j)
+ enddo
+cd write (iout,'(a)') 'Old variables:'
+cd write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
+ call var_to_geom(nvar,varia)
+C Rebuild the chain.
+ call chainbuild
+ MoveType=0
+ nbond=0
+ lprint=.true.
+C Decide whether to generate a random conformation or perturb the old one
+ RandOrPert=ran_number(0.0D0,1.0D0)
+ if (RandOrPert.gt.RanFract) then
+ if (print_mc.gt.0)
+ & write (iout,'(a)') 'Perturbation-generated conformation.'
+ call perturb(error,lprint,MoveType,nbond,1.0D0)
+ if (error) goto 20
+ if (MoveType.lt.1 .or. MoveType.gt.MaxMoveType) then
+ write (iout,'(/a,i7,a/)') 'Error - unknown MoveType=',
+ & MoveType,' returned from PERTURB.'
+ goto 20
+ endif
+ call chainbuild
+ else
+ MoveType=0
+ moves(0)=moves(0)+1
+ nstart_grow=iran_num(3,nres)
+ if (print_mc.gt.0)
+ & write (iout,'(2a,i3)') 'Random-generated conformation',
+ & ' - chain regrown from residue',nstart_grow
+ call gen_rand_conf(nstart_grow,*30)
+ endif
+ call geom_to_var(nvar,varia)
+cd write (iout,'(a)') 'New variables:'
+cd write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
+ ngen=ngen+1
+ if (print_mc.gt.0) write (iout,'(a,i5,a,i10,a,i10)')
+ & 'Processor',MyId,' trial move',ntrial,' total generated:',ngen
+ if (print_mc.gt.0) write (*,'(a,i5,a,i10,a,i10)')
+ & 'Processor',MyId,' trial move',ntrial,' total generated:',ngen
+ call etotal(energia(0))
+ etot = energia(0)
+c call enerprint(energia(0))
+c write (iout,'(2(a,1pe14.5))') 'Etot=',Etot,' Elowest=',Elowest
+ if (etot-elowest.gt.overlap_cut) then
+ write (iout,'(a,i5,a,1pe14.5)') 'Iteration',it,
+ & ' Overlap detected in the current conf.; energy is',etot
+ neneval=neneval+1
+ accepted=.false.
+ noverlap=noverlap+1
+ if (noverlap.gt.maxoverlap) then
+ write (iout,'(a)') 'Too many overlapping confs.'
+ goto 20
+ endif
+ else
+ if (minim) then
+ call minimize(etot,varia,iretcode,nfun)
+cd write (iout,'(a)') 'Variables after minimization:'
+cd write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
+ call etotal(energia(0))
+ etot = energia(0)
+ neneval=neneval+nfun+1
+ endif
+ if (print_mc.gt.2) then
+ write (iout,'(a)') 'Total energies of trial conf:'
+ call enerprint(energia(0))
+ else if (print_mc.eq.1) then
+ write (iout,'(a,i6,a,1pe16.6)')
+ & 'Trial conformation:',ngen,' energy:',etot
+ endif
+C--------------------------------------------------------------------------
+C... Acceptance test
+C--------------------------------------------------------------------------
+ accepted=.false.
+ if (WhatsUp.eq.0)
+ & call accepting(etot,eold,scur,sold,varia,varold,
+ & accepted)
+ if (accepted) then
+ nacc=nacc+1
+ nacc_tot=nacc_tot+1
+ if (elowest.gt.etot) elowest=etot
+ if (ehighest.lt.etot) ehighest=etot
+ moves_acc(MoveType)=moves_acc(MoveType)+1
+ if (MoveType.eq.1) then
+ nbond_acc(nbond)=nbond_acc(nbond)+1
+ endif
+C Check against conformation repetitions.
+ irep=conf_comp(varia,etot)
+#if defined(AIX) || defined(PGI)
+ open (istat,file=statname,position='append')
+#else
+ open (istat,file=statname,access='append')
+#endif
+ if (refstr) then
+ kkk=1
+ call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,kkk),
+ & nsup,
+ & przes,obr,non_conv)
+ rms=dsqrt(rms)
+ call contact(.false.,ncont,icont,co)
+ frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
+ if (print_mc.gt.0)
+ & write (iout,'(a,f8.3,a,f8.3,a,f8.3)')
+ & 'RMS deviation from the reference structure:',rms,
+ & ' % of native contacts:',frac*100,' contact order:',co
+ if (print_stat)
+ & write (istat,'(i5,11(1pe14.5))') it,
+ & (energia(print_order(i)),i=1,nprint_ene),etot,
+ & rms,frac,co
+ elseif (print_stat) then
+ write (istat,'(i5,10(1pe14.5))') it,
+ & (energia(print_order(i)),i=1,nprint_ene),etot
+ endif
+ close(istat)
+ if (print_mc.gt.1)
+ & call statprint(nacc,nfun,iretcode,etot,elowest)
+C Print internal coordinates.
+ if (print_int) call briefout(nacc,etot)
+#ifdef MPL
+ if (MyID.ne.MasterID) then
+ call recv_stop_sig(Kwita)
+cd print *,'Processor:',MyID,' STOP=',Kwita
+ if (irep.eq.0) then
+ call send_MCM_info(2)
+ else
+ call send_MCM_info(1)
+ endif
+ endif
+#endif
+C Store the accepted conf. and its energy.
+ eold=etot
+ sold=scur
+ do i=1,nvar
+ varold(i)=varia(i)
+ enddo
+ if (irep.eq.0) then
+ irep=nsave+1
+cd write (iout,*) 'Accepted conformation:'
+cd write (iout,*) (rad2deg*varia(i),i=1,nphi)
+ if (minim) call zapis(varia,etot)
+ do i=1,n_ene
+ ener(i,nsave)=energia(i)
+ enddo
+ ener(n_ene+1,nsave)=etot
+ ener(n_ene+2,nsave)=frac
+ endif
+ nminima(irep)=nminima(irep)+1.0D0
+c print *,'irep=',irep,' nminima=',nminima(irep)
+#ifdef MPL
+ if (Kwita.eq.0) call recv_stop_sig(kwita)
+#endif
+ endif ! accepted
+ endif ! overlap
+#ifdef MPL
+ if (MyID.eq.MasterID) then
+ call receive_MCM_info
+ if (nacc_tot.ge.maxacc) accepted=.true.
+ endif
+#endif
+ if (ntrial.gt.maxtrial_iter .and. npool.gt.0) then
+C Take a conformation from the pool
+ ii=iran_num(1,npool)
+ do i=1,nvar
+ varia(i)=xpool(i,ii)
+ enddo
+ write (iout,*) 'Iteration',it,' max. # of trials exceeded.'
+ write (iout,*)
+ & 'Take conformation',ii,' from the pool energy=',epool(ii)
+ if (print_mc.gt.2)
+ & write (iout,'(10f8.3)') (rad2deg*varia(i),i=1,nvar)
+ ntrial=0
+ endif ! (ntrial.gt.maxtrial_iter .and. npool.gt.0)
+ 30 continue
+ enddo ! accepted
+#ifdef MPL
+ if (MyID.eq.MasterID) then
+ call receive_MCM_info
+ endif
+ if (Kwita.eq.0) call recv_stop_sig(kwita)
+#endif
+ if (ovrtim()) WhatsUp=-1
+cd write (iout,*) 'WhatsUp=',WhatsUp,' Kwita=',Kwita
+ not_done = (nacc_tot.lt.maxacc) .and. (WhatsUp.eq.0)
+ & .and. (Kwita.eq.0)
+cd write (iout,*) 'not_done=',not_done
+#ifdef MPL
+ if (Kwita.lt.0) then
+ print *,'Processor',MyID,
+ & ' has received STOP signal =',Kwita,' in EntSamp.'
+cd print *,'not_done=',not_done
+ if (Kwita.lt.-1) WhatsUp=Kwita
+ else if (nacc_tot.ge.maxacc) then
+ print *,'Processor',MyID,' calls send_stop_sig,',
+ & ' because a sufficient # of confs. have been collected.'
+cd print *,'not_done=',not_done
+ call send_stop_sig(-1)
+ else if (WhatsUp.eq.-1) then
+ print *,'Processor',MyID,
+ & ' calls send_stop_sig because of timeout.'
+cd print *,'not_done=',not_done
+ call send_stop_sig(-2)
+ endif
+#endif
+ enddo ! not_done
+
+C-----------------------------------------------------------------
+C... Construct energy histogram & update entropy
+C-----------------------------------------------------------------
+ go to 21
+ 20 WhatsUp=-3
+#ifdef MPL
+ write (iout,*) 'Processor',MyID,
+ & ' is broadcasting ERROR-STOP signal.'
+ write (*,*) 'Processor',MyID,
+ & ' is broadcasting ERROR-STOP signal.'
+ call send_stop_sig(-3)
+#endif
+ 21 continue
+#ifdef MPL
+ if (MyID.eq.MasterID) then
+c call receive_MCM_results
+ call receive_energies
+#endif
+ do i=1,nsave
+ if (esave(i).lt.elowest) elowest=esave(i)
+ if (esave(i).gt.ehighest) ehighest=esave(i)
+ enddo
+ write (iout,'(a,i10)') '# of accepted confs:',nacc_tot
+ write (iout,'(a,f10.5,a,f10.5)') 'Lowest energy:',elowest,
+ & ' Highest energy',ehighest
+ if (isweep.eq.1 .and. .not.ent_read) then
+ emin=elowest
+ emax=ehighest
+ write (iout,*) 'EMAX=',emax
+ indminn=0
+ indmaxx=(ehighest-emin)/delte
+ indmin=indminn
+ indmax=indmaxx
+ do i=-max_ene,max_ene
+ entropy(i)=(emin+i*delte)*betbol
+ enddo
+ ent_read=.true.
+ else
+ indmin=(elowest-emin)/delte
+ indmax=(ehighest-emin)/delte
+ if (indmin.lt.indminn) indminn=indmin
+ if (indmax.gt.indmaxx) indmaxx=indmax
+ endif
+ write(iout,*)'indminn=',indminn,' indmaxx=',indmaxx
+C Construct energy histogram
+ do i=1,nsave
+ inde=(esave(i)-emin)/delte
+ nhist(inde)=nhist(inde)+nminima(i)
+ enddo
+C Update entropy (density of states)
+ do i=indmin,indmax
+ if (nhist(i).gt.0) then
+ entropy(i)=entropy(i)+dlog(nhist(i)+0.0D0)
+ endif
+ enddo
+Cd do i=indmaxx+1
+Cd entropy(i)=1.0D+10
+Cd enddo
+ write (iout,'(/80(1h*)/a,i2/80(1h*)/)')
+ & 'End of macroiteration',isweep
+ write (iout,'(a,f10.5,a,f10.5)') 'Elowest=',elowest,
+ & ' Ehighest=',ehighest
+ write (iout,'(a)') 'Frequecies of minima'
+ do i=1,nsave
+ write (iout,'(i5,f5.0,f10.5)') i,nminima(i),esave(i)
+ enddo
+ write (iout,'(/a)') 'Energy histogram'
+ do i=indminn,indmaxx
+ write (iout,'(i5,2f10.5)') i,emin+i*delte,nhist(i)
+ enddo
+ write (iout,'(/a)') 'Entropy'
+ do i=indminn,indmaxx
+ write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
+ enddo
+C-----------------------------------------------------------------
+C... End of energy histogram construction
+C-----------------------------------------------------------------
+#ifdef MPL
+ entropy(-max_ene-4)=dfloat(indminn)
+ entropy(-max_ene-3)=dfloat(indmaxx)
+ entropy(-max_ene-2)=emin
+ entropy(-max_ene-1)=emax
+ call send_MCM_update
+cd print *,entname,ientout
+ open (ientout,file=entname,status='unknown')
+ write (ientout,'(2i5,2e25.17)') indminn,indmaxx,emin,emax
+ do i=indminn,indmaxx
+ write (ientout,'(i5,f10.5,f20.15)') i,emin+i*delte,entropy(i)
+ enddo
+ close(ientout)
+ else
+ write (iout,'(a)') 'Frequecies of minima'
+ do i=1,nsave
+ write (iout,'(i5,f5.0,f10.5)') i,nminima(i),esave(i)
+ enddo
+c call send_MCM_results
+ call send_energies
+ call receive_MCM_update
+ indminn=entropy(-max_ene-4)
+ indmaxx=entropy(-max_ene-3)
+ emin=entropy(-max_ene-2)
+ emax=entropy(-max_ene-1)
+ write (iout,*) 'Received from master:'
+ write (iout,*) 'indminn=',indminn,' indmaxx=',indmaxx,
+ & ' emin=',emin,' emax=',emax
+ write (iout,'(/a)') 'Entropy'
+ do i=indminn,indmaxx
+ write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
+ enddo
+ endif
+ if (WhatsUp.lt.-1) return
+#else
+ if (ovrtim() .or. WhatsUp.lt.0) return
+#endif
+
+ write (iout,'(/80(1h*)/20x,a)') 'Summary run statistics:'
+ call statprint(nacc,nfun,iretcode,etot,elowest)
+ write (iout,'(a)')
+ & 'Statistics of multiple-bond motions. Total motions:'
+ write (iout,'(16i5)') (nbond_move(i),i=1,Nbm)
+ write (iout,'(a)') 'Accepted motions:'
+ write (iout,'(16i5)') (nbond_acc(i),i=1,Nbm)
+ write (iout,'(a,i10)') 'Number of chain regrowths:',nregrow
+ write (iout,'(a,i10)') 'Accepted chain regrowths:',nregrow_acc
+
+C---------------------------------------------------------------------------
+ ENDDO ! ISWEEP
+C---------------------------------------------------------------------------
+
+ runtime=tcpu()
+
+ if (isweep.eq.nsweep .and. it.ge.maxacc)
+ &write (iout,'(/80(1h*)/20x,a/80(1h*)/)') 'All iterations done.'
+ return
+ end
+c------------------------------------------------------------------------------
+ subroutine accepting(ecur,eold,scur,sold,x,xold,accepted)
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.MCM'
+ include 'COMMON.MCE'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.VAR'
+#ifdef MPL
+ include 'COMMON.INFO'
+#endif
+ include 'COMMON.GEO'
+ double precision ecur,eold,xx,ran_number,bol
+ double precision x(maxvar),xold(maxvar)
+ double precision tole /1.0D-1/, tola /5.0D0/
+ logical accepted
+C Check if the conformation is similar.
+cd write (iout,*) 'Enter ACCEPTING'
+cd write (iout,*) 'Old PHI angles:'
+cd write (iout,*) (rad2deg*xold(i),i=1,nphi)
+cd write (iout,*) 'Current angles'
+cd write (iout,*) (rad2deg*x(i),i=1,nphi)
+cd ddif=dif_ang(nphi,x,xold)
+cd write (iout,*) 'Angle norm:',ddif
+cd write (iout,*) 'ecur=',ecur,' emax=',emax
+ if (ecur.gt.emax) then
+ accepted=.false.
+ if (print_mc.gt.0)
+ & write (iout,'(a)') 'Conformation rejected as too high in energy'
+ return
+ else if (dabs(ecur-eold).lt.tole .and.
+ & dif_ang(nphi,x,xold).lt.tola) then
+ accepted=.false.
+ if (print_mc.gt.0)
+ & write (iout,'(a)') 'Conformation rejected as too similar'
+ return
+ endif
+C Else evaluate the entropy of the conf and compare it with that of the previous
+C one.
+ indecur=(ecur-emin)/delte
+ if (iabs(indecur).gt.max_ene) then
+ write (iout,'(a,2i5)')
+ & 'Accepting: Index out of range:',indecur
+ scur=1000.0D0
+ else if (indecur.eq.indmaxx) then
+ scur=entropy(indecur)
+ if (print_mc.gt.0) write (iout,*)'Energy boundary reached',
+ & indmaxx,indecur,entropy(indecur)
+ else
+ deix=ecur-(emin+indecur*delte)
+ dent=entropy(indecur+1)-entropy(indecur)
+ scur=entropy(indecur)+(dent/delte)*deix
+ endif
+cd print *,'Processor',MyID,' ecur=',ecur,' indecur=',indecur,
+cd & ' scur=',scur,' eold=',eold,' sold=',sold
+cd print *,'deix=',deix,' dent=',dent,' delte=',delte
+ if (print_mc.gt.1) then
+ write(iout,*)'ecur=',ecur,' indecur=',indecur,' scur=',scur
+ write(iout,*)'eold=',eold,' sold=',sold
+ endif
+ if (scur.le.sold) then
+ accepted=.true.
+ else
+C Else carry out acceptance test
+ xx=ran_number(0.0D0,1.0D0)
+ xxh=scur-sold
+ if (xxh.gt.50.0D0) then
+ bol=0.0D0
+ else
+ bol=exp(-xxh)
+ endif
+ if (bol.gt.xx) then
+ accepted=.true.
+ if (print_mc.gt.0) write (iout,'(a)')
+ & 'Conformation accepted.'
+ else
+ accepted=.false.
+ if (print_mc.gt.0) write (iout,'(a)')
+ & 'Conformation rejected.'
+ endif
+ endif
+ return
+ end
+c-----------------------------------------------------------------------------
+ subroutine read_pool
+ implicit real*8 (a-h,o-z)
+ include 'DIMENSIONS'
+ include 'COMMON.IOUNITS'
+ include 'COMMON.GEO'
+ include 'COMMON.MCM'
+ include 'COMMON.MCE'
+ include 'COMMON.VAR'
+ double precision varia(maxvar)
+ print '(a)','Call READ_POOL'
+ do npool=1,max_pool
+ print *,'i=',i
+ read (intin,'(i5,f10.5)',end=10,err=10) iconf,epool(npool)
+ if (epool(npool).eq.0.0D0) goto 10
+ call read_angles(intin,*10)
+ call geom_to_var(nvar,xpool(1,npool))
+ enddo
+ goto 11
+ 10 npool=npool-1
+ 11 write (iout,'(a,i5)') 'Number of pool conformations:',npool
+ if (print_mc.gt.2) then
+ do i=1,npool
+ write (iout,'(a,i5,a,1pe14.5)') 'Pool conformation',i,' energy',
+ & epool(i)
+ write (iout,'(10f8.3)') (rad2deg*xpool(j,i),j=1,nvar)
+ enddo
+ endif ! (print_mc.gt.2)
+ return
+ end
projects / unres.git / blobdiff