+++ /dev/null
- subroutine monte_carlo
-C Does Boltzmann and entropic sampling without energy minimization
- 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,nbins
- integer conf_comp
- double precision RandOrPert
- double precision varia(maxvar),elowest,elowest1,
- & ehighest,ehighest1,eold
- double precision przes(3),obr(3,3)
- double precision varold(maxvar)
- logical non_conv
- integer moves1(-1:MaxMoveType+1,0:MaxProcs-1),
- & moves_acc1(-1:MaxMoveType+1,0:MaxProcs-1)
-#ifdef MPL
- double precision etot_temp,etot_all(0:MaxProcs)
- external d_vadd,d_vmin,d_vmax
- double precision entropy1(-max_ene:max_ene),
- & nhist1(-max_ene:max_ene)
- integer nbond_move1(maxres*(MaxProcs+1)),
- & nbond_acc1(maxres*(MaxProcs+1)),itemp(2)
-#endif
- double precision var_lowest(maxvar)
- double precision energia(0:n_ene),energia_ave(0:n_ene)
-C
- write(iout,'(a,i8,2x,a,f10.5)')
- & 'pool_read_freq=',pool_read_freq,' pool_fraction=',pool_fraction
- open (istat,file=statname)
- WhatsUp=0
- indminn=-max_ene
- indmaxx=max_ene
- facee=1.0D0/(maxacc*delte)
-C Number of bins in energy histogram
- nbins=e_up/delte-1
- write (iout,*) 'NBINS=',nbins
- 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)=0.0D0
- 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
- elowest=1.0D+10
- ehighest=-1.0D+10
-C----------------------------------------------------------------------------
-C Entropy-sampling simulations with continually updated entropy;
-C set NSWEEP=1 for Boltzmann sampling.
-C Loop thru simulations
-C----------------------------------------------------------------------------
- DO ISWEEP=1,NSWEEP
-C
-C Initialize the IFINISH array.
-C
-#ifdef MPL
- do i=1,nctasks
- ifinish(i)=0
- enddo
-#endif
-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
- nbond_acc(i)=0
- enddo
-C Initialize total and accepted number of moves of various kind.
- do i=-1,MaxMoveType
- moves(i)=0
- moves_acc(i)=0
- enddo
-C Total number of energy evaluations.
- neneval=0
- nfun=0
-C----------------------------------------------------------------------------
-C Take a conformation from the pool
-C----------------------------------------------------------------------------
- rewind(istat)
- write (iout,*) 'emin=',emin,' emax=',emax
- 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 if (isweep.gt.1) then
- if (eold.lt.emax) then
- do i=1,nvar
- varia(i)=varold(i)
- enddo
- else
- do i=1,nvar
- varia(i)=var_lowest(i)
- enddo
- endif
- call var_to_geom(nvar,varia)
- endif
-C----------------------------------------------------------------------------
-C Compute and print initial energies.
-C----------------------------------------------------------------------------
- nsave=0
- 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 geom_to_var(nvar,varia)
- call etotal(energia(0))
- etot = energia(0)
- call enerprint(energia(0))
- if (refstr) then
- call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup),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,'(i10,16(1pe14.5))') 0,
- & (energia(print_order(i)),i=1,nprint_ene),
- & etot,rms,frac,co
- else
- write (istat,'(i10,14(1pe14.5))') 0,
- & (energia(print_order(i)),i=1,nprint_ene),etot
- endif
-c close(istat)
- neneval=neneval+1
- if (.not. ent_read) then
-C Initialize the entropy array
-#ifdef MPL
-C Collect total energies from other processors.
- etot_temp=etot
- etot_all(0)=etot
- call mp_gather(etot_temp,etot_all,8,MasterID,cgGroupID)
- if (MyID.eq.MasterID) then
-C Get the lowest and the highest energy.
- print *,'MASTER: etot_temp: ',(etot_all(i),i=0,nprocs-1),
- & ' emin=',emin,' emax=',emax
- emin=1.0D10
- emax=-1.0D10
- do i=0,nprocs
- if (emin.gt.etot_all(i)) emin=etot_all(i)
- if (emax.lt.etot_all(i)) emax=etot_all(i)
- enddo
- emax=emin+e_up
- endif ! MyID.eq.MasterID
- etot_all(1)=emin
- etot_all(2)=emax
- print *,'Processor',MyID,' calls MP_BCAST to send/recv etot_all'
- call mp_bcast(etot_all(1),16,MasterID,cgGroupID)
- print *,'Processor',MyID,' MP_BCAST to send/recv etot_all ended'
- if (MyID.ne.MasterID) then
- print *,'Processor:',MyID,etot_all(1),etot_all(2),
- & etot_all(1),etot_all(2)
- emin=etot_all(1)
- emax=etot_all(2)
- endif ! MyID.ne.MasterID
- write (iout,*) 'After MP_GATHER etot_temp=',
- & etot_temp,' emin=',emin
-#else
- emin=etot
- emax=emin+e_up
- indminn=0
- indmin=0
-#endif
- IF (MULTICAN) THEN
-C Multicanonical sampling - start from Boltzmann distribution
- do i=-max_ene,max_ene
- entropy(i)=(emin+i*delte)*betbol
- enddo
- ELSE
-C Entropic sampling - start from uniform distribution of the density of states
- do i=-max_ene,max_ene
- entropy(i)=0.0D0
- enddo
- ENDIF ! MULTICAN
- write (iout,'(/a)') 'Initial entropy'
- do i=indminn,indmaxx
- write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
- enddo
- if (isweep.eq.1) then
- emax=emin+e_up
- indminn=0
- indmin=0
- indmaxx=indminn+nbins
- indmax=indmaxx
- endif ! isweep.eq.1
- endif ! .not. 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=.true.
-#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
- call entropia(eold,sold,indeold)
-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
-C Main loop.
-c----------------------------------------------------------------------------
-C Zero out average energies
- do i=0,n_ene
- energia_ave(i)=0.0d0
- enddo
-C Initialize energy histogram
- do i=-max_ene,max_ene
- nhist(i)=0.0D0
- enddo ! i
-C Zero out iteration counter.
- it=0
- do j=1,nvar
- varold(j)=varia(j)
- enddo
-C Begin MC iteration loop.
- do while (not_done)
- it=it+1
-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
- call var_to_geom(nvar,varia)
-C Rebuild the chain.
- call chainbuild
- MoveType=0
- nbond=0
- lprint=.true.
-C Decide whether to take a conformation from the pool or generate/perturb one
-C randomly
- from_pool=ran_number(0.0D0,1.0D0)
- if (npool.gt.0 .and. from_pool.lt.pool_fraction) then
-C Throw a dice to choose the conformation from the pool
- ii=iran_num(1,npool)
- do i=1,nvar
- varia(i)=xpool(i,ii)
- enddo
- call var_to_geom(nvar,varia)
- call chainbuild
-cd call intout
-cd write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
- if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it)
- & write (iout,'(a,i3,a,f10.5)')
- & 'Try conformation',ii,' from the pool energy=',epool(ii)
- MoveType=-1
- moves(-1)=moves(-1)+1
- else
-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 .and. (it/print_freq)*print_freq.eq.it)
- & write (iout,'(a)') 'Perturbation-generated conformation.'
- call perturb(error,lprint,MoveType,nbond,0.1D0)
- 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 .and. (it/print_freq)*print_freq.eq.it)
- & 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)
- endif ! pool
-Cd write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
- ngen=ngen+1
- if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it)
- & write (iout,'(a,i5,a,i10,a,i10)')
- & 'Processor',MyId,' trial move',ntrial,' total generated:',ngen
- if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it)
- & write (*,'(a,i5,a,i10,a,i10)')
- & 'Processor',MyId,' trial move',ntrial,' total generated:',ngen
- call etotal(energia(0))
- etot = energia(0)
- neneval=neneval+1
-cd call enerprint(energia(0))
-cd write(iout,*)'it=',it,' etot=',etot
- if (etot-elowest.gt.overlap_cut) then
- if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it)
- & write (iout,'(a,i5,a,1pe14.5)') 'Iteration',it,
- & ' Overlap detected in the current conf.; energy is',etot
- accepted=.false.
- noverlap=noverlap+1
- if (noverlap.gt.maxoverlap) then
- write (iout,'(a)') 'Too many overlapping confs.'
- goto 20
- endif
- else
-C--------------------------------------------------------------------------
-C... Acceptance test
-C--------------------------------------------------------------------------
- accepted=.false.
- if (WhatsUp.eq.0)
- & call accept_mc(it,etot,eold,scur,sold,varia,varold,accepted)
- if (accepted) then
- nacc=nacc+1
- nacc_tot=nacc_tot+1
- if (elowest.gt.etot) then
- elowest=etot
- do i=1,nvar
- var_lowest(i)=varia(i)
- enddo
- endif
- 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 Compare with reference structure.
- if (refstr) then
- call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup),
- & nsup,przes,obr,non_conv)
- rms=dsqrt(rms)
- call contact(.false.,ncont,icont,co)
- frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
- endif ! refstr
-C
-C Periodically save average energies and confs.
-C
- do i=0,n_ene
- energia_ave(i)=energia_ave(i)+energia(i)
- enddo
- moves(MaxMoveType+1)=nmove
- moves_acc(MaxMoveType+1)=nacc
- IF ((it/save_frequency)*save_frequency.eq.it) THEN
- do i=0,n_ene
- energia_ave(i)=energia_ave(i)/save_frequency
- enddo
- etot_ave=energia_ave(0)
-C#ifdef AIX
-C open (istat,file=statname,position='append')
-C#else
-C open (istat,file=statname,access='append')
-Cendif
- if (print_mc.gt.0)
- & write (iout,'(80(1h*)/20x,a,i20)')
- & 'Iteration #',it
- if (refstr .and. print_mc.gt.0) then
- 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
- endif
- if (print_stat) then
- if (refstr) then
- write (istat,'(i10,10(1pe14.5))') it,
- & (energia_ave(print_order(i)),i=1,nprint_ene),
- & etot_ave,rms_ave,frac_ave
- else
- write (istat,'(i10,10(1pe14.5))') it,
- & (energia_ave(print_order(i)),i=1,nprint_ene),
- & etot_ave
- endif
- endif
-c close(istat)
- if (print_mc.gt.0)
- & call statprint(nacc,nfun,iretcode,etot,elowest)
-C Print internal coordinates.
- if (print_int) call briefout(nacc,etot)
- do i=0,n_ene
- energia_ave(i)=0.0d0
- enddo
- ENDIF ! ( (it/save_frequency)*save_frequency.eq.it)
-C Update histogram
- inde=icialosc((etot-emin)/delte)
- nhist(inde)=nhist(inde)+1.0D0
-#ifdef MPL
- if ( (it/message_frequency)*message_frequency.eq.it
- & .and. (MyID.ne.MasterID) ) then
- call recv_stop_sig(Kwita)
- call send_MCM_info(message_frequency)
- endif
-#endif
-C Store the accepted conf. and its energy.
- eold=etot
- sold=scur
- do i=1,nvar
- varold(i)=varia(i)
- enddo
-#ifdef MPL
- if (Kwita.eq.0) call recv_stop_sig(kwita)
-#endif
- endif ! accepted
- endif ! overlap
-#ifdef MPL
- if (MyID.eq.MasterID .and.
- & (it/message_frequency)*message_frequency.eq.it) then
- call receive_MC_info
- if (nacc_tot.ge.maxacc) accepted=.true.
- endif
-#endif
-C if ((ntrial.gt.maxtrial_iter
-C & .or. (it/pool_read_freq)*pool_read_freq.eq.it)
-C & .and. npool.gt.0) then
-C Take a conformation from the pool
-C ii=iran_num(1,npool)
-C do i=1,nvar
-C varold(i)=xpool(i,ii)
-C enddo
-C if (ntrial.gt.maxtrial_iter)
-C & write (iout,*) 'Iteration',it,' max. # of trials exceeded.'
-C write (iout,*)
-C & 'Take conformation',ii,' from the pool energy=',epool(ii)
-C if (print_mc.gt.2)
-C & write (iout,'(10f8.3)') (rad2deg*varold(i),i=1,nvar)
-C ntrial=0
-C eold=epool(ii)
-C call entropia(eold,sold,indeold)
-C accepted=.true.
-C endif ! (ntrial.gt.maxtrial_iter .and. npool.gt.0)
- 30 continue
- enddo ! accepted
-#ifdef MPL
- if (MyID.eq.MasterID .and.
- & (it/message_frequency)*message_frequency.eq.it) then
- call receive_MC_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
- if (MyID.ne.MasterID) call send_MCM_info(-1)
- 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)
- if (MyID.ne.MasterID) call send_MCM_info(-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)
- if (MyID.ne.MasterID) call send_MCM_info(-1)
- 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)
- if (MyID.ne.MasterID) call send_MCM_info(-1)
-#endif
- 21 continue
- write (iout,'(/a)') 'Energy histogram'
- do i=-100,100
- write (iout,'(i5,2f20.5)') i,emin+i*delte,nhist(i)
- enddo
-#ifdef MPL
-C Wait until every processor has sent complete MC info.
- if (MyID.eq.MasterID) then
- not_done=.true.
- do while (not_done)
-C write (*,*) 'The IFINISH array:'
-C write (*,*) (ifinish(i),i=1,nctasks)
- not_done=.false.
- do i=2,nctasks
- not_done=not_done.or.(ifinish(i).ge.0)
- enddo
- if (not_done) call receive_MC_info
- enddo
- endif
-C Make collective histogram from the work of all processors.
- msglen=(2*max_ene+1)*8
- print *,
- & 'Processor',MyID,' calls MP_REDUCE to send/receive histograms',
- & ' msglen=',msglen
- call mp_reduce(nhist,nhist1,msglen,MasterID,d_vadd,
- & cgGroupID)
- print *,'Processor',MyID,' MP_REDUCE accomplished for histogr.'
- do i=-max_ene,max_ene
- nhist(i)=nhist1(i)
- enddo
-C Collect min. and max. energy
- print *,
- &'Processor',MyID,' calls MP_REDUCE to send/receive energy borders'
- call mp_reduce(elowest,elowest1,8,MasterID,d_vmin,cgGroupID)
- call mp_reduce(ehighest,ehighest1,8,MasterID,d_vmax,cgGroupID)
- print *,'Processor',MyID,' MP_REDUCE accomplished for energies.'
- IF (MyID.eq.MasterID) THEN
- elowest=elowest1
- ehighest=ehighest1
-#endif
- write (iout,'(a,i10)') '# of accepted confs:',nacc_tot
- write (iout,'(a,f10.5,a,f10.5)') 'Lowest energy:',elowest,
- & ' Highest energy',ehighest
- indmin=icialosc((elowest-emin)/delte)
- imdmax=icialosc((ehighest-emin)/delte)
- if (indmin.lt.indminn) then
- emax=emin+indmin*delte+e_up
- indmaxx=indmin+nbins
- indminn=indmin
- endif
- if (.not.ent_read) ent_read=.true.
- write(iout,*)'indminn=',indminn,' indmaxx=',indmaxx
-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
- 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)') 'Energy histogram'
- do i=indminn,indmaxx
- write (iout,'(i5,2f20.5)') i,emin+i*delte,nhist(i)
- enddo
- write (iout,'(/a)') 'Entropy'
- do i=indminn,indmaxx
- write (iout,'(i5,2f20.5)') i,emin+i*delte,entropy(i)
- enddo
-C-----------------------------------------------------------------
-C... End of energy histogram construction
-C-----------------------------------------------------------------
-#ifdef MPL
- ELSE
- if (.not. ent_read) ent_read=.true.
- ENDIF ! MyID .eq. MaterID
- if (MyID.eq.MasterID) then
- itemp(1)=indminn
- itemp(2)=indmaxx
- endif
- print *,'before mp_bcast processor',MyID,' indminn=',indminn,
- & ' indmaxx=',indmaxx,' itemp=',itemp(1),itemp(2)
- call mp_bcast(itemp(1),8,MasterID,cgGroupID)
- call mp_bcast(emax,8,MasterID,cgGroupID)
- print *,'after mp_bcast processor',MyID,' indminn=',indminn,
- & ' indmaxx=',indmaxx,' itemp=',itemp(1),itemp(2)
- if (MyID .ne. MasterID) then
- indminn=itemp(1)
- indmaxx=itemp(2)
- endif
- msglen=(indmaxx-indminn+1)*8
- print *,'processor',MyID,' calling mp_bcast msglen=',msglen,
- & ' indminn=',indminn,' indmaxx=',indmaxx,' isweep=',isweep
- call mp_bcast(entropy(indminn),msglen,MasterID,cgGroupID)
- IF(MyID.eq.MasterID .and. .not. ovrtim() .and. WhatsUp.ge.0)THEN
- 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,*) '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 ! MyID.eq.MasterID
- print *,'Processor',MyID,' calls MP_GATHER'
- call mp_gather(nbond_move,nbond_move1,4*Nbm,MasterID,
- & cgGroupID)
- call mp_gather(nbond_acc,nbond_acc1,4*Nbm,MasterID,
- & cgGroupID)
- print *,'Processor',MyID,' MP_GATHER call accomplished'
- if (MyID.eq.MasterID) then
-
- write (iout,'(/80(1h*)/20x,a)') 'Summary run statistics:'
- call statprint(nacc_tot,nfun,iretcode,etot,elowest)
- write (iout,'(a)')
- & 'Statistics of multiple-bond motions. Total motions:'
- write (iout,'(8i10)') (nbond_move(i),i=1,Nbm)
- write (iout,'(a)') 'Accepted motions:'
- write (iout,'(8i10)') (nbond_acc(i),i=1,Nbm)
-
- write (iout,'(a)')
- & 'Statistics of multi-bond moves of respective processors:'
- do iproc=1,Nprocs-1
- do i=1,Nbm
- ind=iproc*nbm+i
- nbond_move(i)=nbond_move(i)+nbond_move1(ind)
- nbond_acc(i)=nbond_acc(i)+nbond_acc1(ind)
- enddo
- enddo
- do iproc=0,NProcs-1
- write (iout,*) 'Processor',iproc,' nbond_move:',
- & (nbond_move1(iproc*nbm+i),i=1,Nbm),
- & ' nbond_acc:',(nbond_acc1(iproc*nbm+i),i=1,Nbm)
- enddo
- endif
- call mp_gather(moves,moves1,4*(MaxMoveType+3),MasterID,
- & cgGroupID)
- call mp_gather(moves_acc,moves_acc1,4*(MaxMoveType+3),
- & MasterID,cgGroupID)
- if (MyID.eq.MasterID) then
- do iproc=1,Nprocs-1
- do i=-1,MaxMoveType+1
- moves(i)=moves(i)+moves1(i,iproc)
- moves_acc(i)=moves_acc(i)+moves_acc1(i,iproc)
- enddo
- enddo
- nmove=0
- do i=0,MaxMoveType+1
- nmove=nmove+moves(i)
- enddo
- do iproc=0,NProcs-1
- write (iout,*) 'Processor',iproc,' moves',
- & (moves1(i,iproc),i=0,MaxMoveType+1),
- & ' moves_acc:',(moves_acc1(i,iproc),i=0,MaxMoveType+1)
- enddo
- endif
-#else
- 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)
-#endif
- write (iout,'(/80(1h*)/20x,a)') 'Summary run statistics:'
- call statprint(nacc_tot,nfun,iretcode,etot,elowest)
- write (iout,'(a)')
- & 'Statistics of multiple-bond motions. Total motions:'
- write (iout,'(8i10)') (nbond_move(i),i=1,Nbm)
- write (iout,'(a)') 'Accepted motions:'
- write (iout,'(8i10)') (nbond_acc(i),i=1,Nbm)
- if (ovrtim() .or. WhatsUp.lt.0) return
-
-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 accept_mc(it,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)
- 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 .and. (it/print_freq)*print_freq.eq.it)
- & write (iout,'(a)') 'Conformation rejected as too high in energy'
- return
- endif
-C Else evaluate the entropy of the conf and compare it with that of the previous
-C one.
- call entropia(ecur,scur,indecur)
-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.0 .and. (it/print_freq)*print_freq.eq.it) then
- write(iout,*)'it=',it,'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 .and. (it/print_freq)*print_freq.eq.it)
- & write (iout,'(a)') 'Conformation accepted.'
- else
- accepted=.false.
- if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it)
- & write (iout,'(a)') 'Conformation rejected.'
- endif
- endif
- return
- end
-c--------------------------------------------------------------------------
- integer function icialosc(x)
- double precision x
- if (x.lt.0.0D0) then
- icialosc=dint(x)-1
- else
- icialosc=dint(x)
- endif
- return
- end
-c--------------------------------------------------------------------------
- subroutine entropia(ecur,scur,indecur)
- implicit real*8 (a-h,o-z)
- include 'DIMENSIONS'
- include 'COMMON.MCM'
- include 'COMMON.MCE'
- include 'COMMON.IOUNITS'
- double precision ecur,scur
- integer indecur
- indecur=icialosc((ecur-emin)/delte)
- if (iabs(indecur).gt.max_ene) then
- if ((it/print_freq)*it.eq.it) write (iout,'(a,2i5)')
- & 'Accepting: Index out of range:',indecur
- scur=1000.0D0
- else if (indecur.ge.indmaxx) then
- scur=entropy(indecur)
- if (print_mc.gt.0 .and. (it/print_freq)*it.eq.it)
- & 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
- return
- end
projects / unres.git / blobdiff