X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD%2Fsrc%2Fmc.F;fp=source%2Funres%2Fsrc_MD%2Fsrc%2Fmc.F;h=0000000000000000000000000000000000000000;hb=0a11a2c4ccee14ed99ae44f2565b270ba8d4bbb6;hp=0f39d488cda4604ada2a5f4b1dcf2eb33244e4ea;hpb=5eb407964903815242c59de10960f42761139e10;p=unres.git diff --git a/source/unres/src_MD/src/mc.F b/source/unres/src_MD/src/mc.F deleted file mode 100644 index 0f39d48..0000000 --- a/source/unres/src_MD/src/mc.F +++ /dev/null @@ -1,819 +0,0 @@ - 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