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