X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2FMCM_MD.f90;fp=source%2Funres%2FMCM_MD.f90;h=0000000000000000000000000000000000000000;hb=2611e37f82b576a1366c2d78fce87c1a55852037;hp=afb31bb4a99a86a7c2d5108f3720ad5d810f8582;hpb=a1e9d5a6b704c90ebd10f86a655780212a880dce;p=unres4.git diff --git a/source/unres/MCM_MD.f90 b/source/unres/MCM_MD.f90 deleted file mode 100644 index afb31bb..0000000 --- a/source/unres/MCM_MD.f90 +++ /dev/null @@ -1,3514 +0,0 @@ - module mcm_md -!----------------------------------------------------------------------------- - use io_units - use names - use math - use geometry_data, only: nres,nvar,rad2deg - use random, only: iran_num,ran_number - use MD_data - use MCM_data - use geometry - use energy - - implicit none -!----------------------------------------------------------------------------- -! Max. number of move types in MCM -! integer,parameter :: maxmovetype=4 -!----------------------------------------------------------------------------- -! Max. number of conformations in Master's cache array - integer,parameter :: max_cache=10 -!----------------------------------------------------------------------------- -! Max. number of stored confs. in MC/MCM simulation -! integer,parameter :: maxsave=20 -!----------------------------------------------------------------------------- -! Number of threads in deformation - integer,parameter :: max_thread=4, max_thread2=2*max_thread -!----------------------------------------------------------------------------- -! Number of structures to compare at t=0 - integer,parameter :: max_threadss=8,max_threadss2=2*max_threadss -!----------------------------------------------------------------------------- -! Max. number of conformations in the pool - integer,parameter :: max_pool=10 -!----------------------------------------------------------------------------- -!----------------------------------------------------------------------------- -! commom.cache -! common /cache/ - integer :: ncache -! integer,dimension(max_cache) :: CachSrc nie używane -! integer,dimension(max_cache) :: isent,iused -! logical :: cache_update -! real(kind=8),dimension(max_cache) :: ecache -! real(kind=8),dimension(:,:),allocatable :: xcache !(maxvar,max_cache) -!----------------------------------------------------------------------------- -! common.mce -! common /mce/ -! real(kind=8) :: emin,emax - real(kind=8),dimension(:),allocatable :: entropy !(-max_ene-4:max_ene) - real(kind=8),dimension(:),allocatable :: nhist !(-max_ene:max_ene) - real(kind=8),dimension(:),allocatable :: nminima !(maxsave) -! logical :: ent_read - logical :: multican - integer :: indminn,indmaxx -! common /pool/ - integer :: npool -! real(kind=8) :: pool_fraction - real(kind=8),dimension(:,:),allocatable :: xpool !(maxvar,max_pool) - real(kind=8),dimension(:),allocatable :: epool !(max_pool) -! common /mce_counters/ -!------------------------------------------------------------------------------ -!... Following COMMON block contains variables controlling motion. -!------------------------------------------------------------------------------ -! common /move/ -! real(kind=8),dimension(0:MaxMoveType) :: sumpro_type !(0:MaxMoveType) - real(kind=8),dimension(:),allocatable :: sumpro_bond !(0:maxres) - integer :: koniecl,Nbm,MaxSideMove!,nmove - integer,dimension(:),allocatable :: nbond_move,nbond_acc !(maxres) -! integer,dimension(-1:MaxMoveType+1) :: moves,moves_acc !(-1:MaxMoveType+1) -! common /accept_stats/ -! integer :: nacc_tot - integer,dimension(:),allocatable :: nacc_part !(0:MaxProcs) !el nie uzywane??? -! common /windows/ -! integer :: nwindow -! integer,dimension(:),allocatable :: winstart,winend,winlen !(maxres) -! common /moveID/ -! character(len=16),dimension(-1:MaxMoveType+1) :: MovTypID !(-1:MaxMoveType+1) -!------------------------------------------------------------------------------ -!... koniecl - the number of bonds to be considered "end bonds" subjected to -!... end moves; -!... Nbm - The maximum length of N-bond segment to be moved; -!... MaxSideMove - maximum number of side chains subjected to local moves -!... simultaneously; -!... nmove - the current number of attempted moves; -!... nbond_move(*) array that stores the total numbers of 2-bond,3-bond,... -!... moves; -!... nendmove - number of endmoves; -!... nbackmove - number of backbone moves; -!... nsidemove - number of local side chain moves; -!... sumpro_type(*) - array that stores the lower and upper boundary of the -!... random-number range that determines the type of move -!... (N-bond, backbone or side chain); -!... sumpro_bond(*) - array that stores the probabilities to perform bond -!... moves of consecutive segment length. -!... winstart(*) - the starting position of the perturbation window; -!... winend(*) - the end position of the perturbation window; -!... winlen(*) - length of the perturbation window; -!... nwindow - the number of perturbation windows (0 - entire chain). -!----------------------------------------------------------------------------- -! -! -!----------------------------------------------------------------------------- - contains -!----------------------------------------------------------------------------- -! compare_s1.F -!----------------------------------------------------------------------------- - subroutine compare_s1(n_thr,num_thread_save,energyx,x,icomp,enetbss,& - coordss,rms_d,modif,iprint) - -! This subroutine compares the new conformation, whose variables are in X -! with the previously accumulated conformations whose energies and variables -! are stored in ENETBSS and COORDSS, respectively. The meaning of other -! variables is as follows: -! -! N_THR - on input the previous # of accumulated confs, on output the current -! # of accumulated confs. -! N_REPEAT - an array that indicates how many times the structure has already -! been used to start the reversed-reversing procedure. Addition of -! a new structure replacement of a structure with a similar, but -! lower-energy structure resets the respective entry in N_REPEAT to zero -! I9 - output unit -! ENERGYX,X - the energy and variables of the new conformations. -! ICOMP - comparison result: -! 0 - the new structure is similar to one of the previous ones and does -! not have a remarkably lower energy and is therefore rejected; -! 1 - the new structure is different and is added to the set, because -! there is still room in the COORDSS and ENETBSS arrays; -! 2 - the new structure is different, but higher in energy than any -! previous one and is therefore rejected -! 3 - there is no more room in the COORDSS and ENETBSS arrays, but -! the new structure is lower in energy than at least the highest- -! energy previous structure and therefore replaces it. -! 9 - the new structure is similar to a number of previous structures, -! but has a remarkably lower energy than any of them; therefore -! replaces all these structures; -! MODIF - a logical variable that shows whether to include the new structure -! in the set of accumulated structures - -! implicit real*8 (a-h,o-z) - use geometry_data - use regularize_, only:fitsq -! include 'DIMENSIONS' -! include 'COMMON.GEO' -! include 'COMMON.VAR' -!rc include 'COMMON.DEFORM' -! include 'COMMON.IOUNITS' -!el#ifdef UNRES -!el use geometry_data !include 'COMMON.CHAIN' -!el#endif - - real(kind=8),dimension(6*nres) :: x,x1 !(maxvar) (maxvar=6*maxres) - real(kind=8) :: przes(3),obrot(3,3) - integer :: list(max_thread) - logical :: non_conv,modif - real(kind=8) :: enetbss(max_threadss) - real(kind=8) :: coordss(6*nres,max_threadss) - -!!! local variables - el - integer :: n_thr,num_thread_save,icomp,minimize_s_flag,iprint - real(kind=8) :: energyx,energyy,rms_d - integer :: nlist,k,kk,j,i,iresult - real(kind=8) :: enex_jp,roznica - - nlist=0 -#ifdef UNRES - call var_to_geom(nvar,x) - call chainbuild - do k=1,2*nres - do kk=1,3 - cref(kk,k,1)=c(kk,k) - enddo - enddo -#endif -! write(iout,*)'*ene=',energyx - j=0 - enex_jp=-1.0d+99 - do i=1,n_thr - do k=1,nvar - x1(k)=coordss(k,i) - enddo - if (iprint.gt.3) then - write (iout,*) 'Compare_ss, i=',i - write (iout,*) 'New structure Energy:',energyx - write (iout,'(10f8.3)') (rad2deg*x(k),k=1,nvar) - write (iout,*) 'Template structure Energy:',enetbss(i) - write (iout,'(10f8.3)') (rad2deg*x1(k),k=1,nvar) - endif - -#ifdef UNRES - call var_to_geom(nvar,x1) - call chainbuild -!d write(iout,*)'C and CREF' -!d write(iout,'(i5,3f10.5,5x,3f10.5)')(k,(c(j,k),j=1,3), -!d & (cref(j,k),j=1,3),k=1,nres) - call fitsq(roznica,c(1,1),cref(1,1,1),nres,przes,obrot,non_conv) - if (non_conv) then - print *,'Problems in FITSQ!!!' - print *,'X' - print '(10f8.3)',(x(k),k=1,nvar) - print *,'X1' - print '(10f8.3)',(x1(k),k=1,nvar) - print *,'C and CREF' - print '(i5,3f10.5,5x,3f10.5)',(k,(c(j,k),j=1,3),& - (cref(j,k,1),j=1,3),k=1,nres) - endif - roznica=dsqrt(dabs(roznica)) - iresult = 1 - if (roznica.lt.rms_d) iresult = 0 -#else - energyy=enetbss(i) -!el call cmprs(x,x1,roznica,energyx,energyy,iresult) -#endif - if (iprint.gt.1) write(iout,'(i5,f10.6,$)') i,roznica -! print '(i5,f8.3)',i,roznica - if(iresult.eq.0) then - nlist = nlist + 1 - list(nlist)=i - if (iprint.gt.1) write(iout,*) - if(energyx.ge.enetbss(i)) then - if (iprint.gt.1) & - write(iout,*)'s*>> structure rejected - same as nr ',i, & - ' RMS',roznica - minimize_s_flag=0 - icomp=0 - go to 1106 - endif - endif - if(energyx.lt.enetbss(i).and.enex_jp.lt.enetbss(i))then - j=i - enex_jp=enetbss(i) - endif - enddo - if (iprint.gt.1) write(iout,*) - if(nlist.gt.0) then - if (modif) then - if (iprint.gt.1) & - write(iout,'(a,i3,$)')'s*>> structure accepted1 - repl nr ',& - list(1) - else - if (iprint.gt.1) & - write(iout,'(a,i3)') & - 's*>> structure accepted1 - would repl nr ',list(1) - endif - icomp=9 - if (.not. modif) goto 1106 - j=list(1) - enetbss(j)=energyx - do i=1,nvar - coordss(i,j)=x(i) - enddo - do j=2,nlist - if (iprint.gt.1) write(iout,'(i3,$)')list(j) - do kk=list(j)+1,nlist - enetbss(kk-1)=enetbss(kk) - do i=1,nvar - coordss(i,kk-1)=coordss(i,kk) - enddo - enddo - enddo - if (iprint.gt.1) write(iout,*) - go to 1106 - endif - if(n_thr.lt.num_thread_save) then - icomp=1 - if (modif) then - if (iprint.gt.1) & - write(iout,*)'s*>> structure accepted - add with nr ',n_thr+1 - else - if (iprint.gt.1) & - write(iout,*)'s*>> structure accepted - would add with nr ',& - n_thr+1 - goto 1106 - endif - n_thr=n_thr+1 - enetbss(n_thr)=energyx - do i=1,nvar - coordss(i,n_thr)=x(i) - enddo - else - if(j.eq.0) then - if (iprint.gt.1) & - write(iout,*)'s*>> structure rejected - too high energy' - icomp=2 - go to 1106 - end if - icomp=3 - if (modif) then - if (iprint.gt.1) & - write(iout,*)'s*>> structure accepted - repl nr ',j - else - if (iprint.gt.1) & - write(iout,*)'s*>> structure accepted - would repl nr ',j - goto 1106 - endif - enetbss(j)=energyx - do i=1,nvar - coordss(i,j)=x(i) - enddo - end if - -1106 continue - return - end subroutine compare_s1 -!----------------------------------------------------------------------------- -! entmcm.F -!----------------------------------------------------------------------------- - subroutine entmcm - - use energy_data - use geometry_data - use MPI_data, only:WhatsUp,MyID - use compare_data, only: ener - use control_data, only: minim,refstr - use io_base - use regularize_, only:fitsq - use control, only: tcpu,ovrtim - use compare - use minimm, only:minimize -! Does modified entropic sampling in the space of minima. -! implicit real*8 (a-h,o-z) -! include 'DIMENSIONS' -! include 'COMMON.IOUNITS' -#ifdef MPL - use MPI_data !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,error,lprint !,ovrtim - integer :: MoveType,nbond -! integer :: conf_comp - real(kind=8) :: RandOrPert - real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres) - real(kind=8) :: elowest,ehighest,eold - real(kind=8) :: przes(3),obr(3,3) - real(kind=8),dimension(6*nres) :: varold !(maxvar) (maxvar=6*maxres) - logical :: non_conv - real(kind=8),dimension(0:n_ene) :: energia,energia_ave - -!!! local variables -el - integer :: i,ii,kkk,it,j,nacc,nfun,ijunk,indmin,indmax,& - ISWEEP,Kwita,iretcode,indeold,iene,noverlap,& - irep,nstart_grow,inde - real(kind=8) :: facee,conste,ejunk,etot,rms,co,frac,& - deix,dent,sold,scur,runtime -! - -! if(.not.allocated(varsave)) allocate(varsave(maxvar,maxsave)) !(maxvar,maxsave) -!d write (iout,*) 'print_mc=',print_mc - WhatsUp=0 - maxtrial_iter=50 -!--------------------------------------------------------------------------- -! Initialize counters. -!--------------------------------------------------------------------------- -! Total number of generated confs. - ngen=0 -! Total number of moves. In general this won't be equal to the number of -! attempted moves, because we may want to reject some "bad" confs just by -! overlap check. - nmove=0 -! Total number of shift (nbond_move(1)), spike, crankshaft, three-bond,... -! motions. -!el allocate(nbond_move(nres)) !(maxres) - - do i=1,nres - nbond_move(i)=0 - enddo -! Initialize total and accepted number of moves of various kind. - do i=0,MaxMoveType - moves(i)=0 - moves_acc(i)=0 - enddo -! 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) -! 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 -! Read the pool of conformations - call read_pool -!---------------------------------------------------------------------------- -! Entropy-sampling simulations with continually updated entropy -! Loop thru simulations -!---------------------------------------------------------------------------- - DO ISWEEP=1,NSWEEP -!---------------------------------------------------------------------------- -! Take a conformation from the pool -!---------------------------------------------------------------------------- - 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) -! Print internal coordinates of the initial conformation - call intout - else - call gen_rand_conf(1,*20) - endif -!---------------------------------------------------------------------------- -! Compute and print initial energies. -!---------------------------------------------------------------------------- - nsave=0 -#ifdef MPL - allocate(nsave_part(nctasks)) - 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) - etot = energia(0) - call enerprint(energia) -! 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) - etot = energia(0) - write (iout,'(/80(1h*)/a/80(1h*))') & - 'Results of the first energy minimization:' - call enerprint(energia) - 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 -! Initialize the entropy array - do i=-max_ene,max_ene - emin=etot -! Uncomment the line below for actual entropic sampling (start with uniform -! energy distribution). -! entropy(i)=0.0D0 -! Uncomment the line below for multicanonical sampling (start with Boltzmann -! 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) -!d write (iout,*) 'indeold=',indeold,' deix=',deix,' dent=',dent -!d write (*,*) 'Processor',MyID,' indeold=',indeold,' deix=',deix, -!d & ' 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 -! NACC is the counter for the accepted conformations of a given processor - nacc=0 -! 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 -! Main loop. -!---------------------------------------------------------------------------- - 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 -! 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 -! Retrieve the angles of previously accepted conformation - do j=1,nvar - varia(j)=varold(j) - enddo -!d write (iout,'(a)') 'Old variables:' -!d write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar) - call var_to_geom(nvar,varia) -! Rebuild the chain. - call chainbuild - MoveType=0 - nbond=0 - lprint=.true. -! 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) -!d write (iout,'(a)') 'New variables:' -!d 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) - etot = energia(0) -! call enerprint(energia(0)) -! 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) -!d write (iout,'(a)') 'Variables after minimization:' -!d write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar) - call etotal(energia) - 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) - else if (print_mc.eq.1) then - write (iout,'(a,i6,a,1pe16.6)') & - 'Trial conformation:',ngen,' energy:',etot - endif -!-------------------------------------------------------------------------- -!... Acceptance test -!-------------------------------------------------------------------------- - 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 -! 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) -! Print internal coordinates. - if (print_int) call briefout(nacc,etot) -#ifdef MPL - if (MyID.ne.MasterID) then - call recv_stop_sig(Kwita) -!d print *,'Processor:',MyID,' STOP=',Kwita - if (irep.eq.0) then - call send_MCM_info(2) - else - call send_MCM_info(1) - endif - endif -#endif -! 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 -!d write (iout,*) 'Accepted conformation:' -!d 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 -! 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 -! 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 -!d write (iout,*) 'WhatsUp=',WhatsUp,' Kwita=',Kwita - not_done = (nacc_tot.lt.maxacc) .and. (WhatsUp.eq.0) & - .and. (Kwita.eq.0) -!d write (iout,*) 'not_done=',not_done -#ifdef MPL - if (Kwita.lt.0) then - print *,'Processor',MyID,& - ' has received STOP signal =',Kwita,' in EntSamp.' -!d 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.' -!d 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.' -!d print *,'not_done=',not_done - call send_stop_sig(-2) - endif -#endif - enddo ! not_done - -!----------------------------------------------------------------- -!... Construct energy histogram & update entropy -!----------------------------------------------------------------- - 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 -! 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 -! Construct energy histogram - do i=1,nsave - inde=(esave(i)-emin)/delte - nhist(inde)=nhist(inde)+nminima(i) - enddo -! 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 -!d do i=indmaxx+1 -!d entropy(i)=1.0D+10 -!d 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 -!----------------------------------------------------------------- -!... End of energy histogram construction -!----------------------------------------------------------------- -#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 -!d 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 -! 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) -!el write (iout,'(a,i10)') 'Number of chain regrowths:',nregrow -!el write (iout,'(a,i10)') 'Accepted chain regrowths:',nregrow_acc - -!--------------------------------------------------------------------------- - ENDDO ! ISWEEP -!--------------------------------------------------------------------------- - - runtime=tcpu() - - if (isweep.eq.nsweep .and. it.ge.maxacc) & - write (iout,'(/80(1h*)/20x,a/80(1h*)/)') 'All iterations done.' - return - end subroutine entmcm -!----------------------------------------------------------------------------- - subroutine accepting(ecur,eold,scur,sold,x,xold,accepted) - - use geometry_data, only: nphi - use energy_data, only: max_ene -! implicit real*8 (a-h,o-z) -! include 'DIMENSIONS' -! include 'COMMON.MCM' -! include 'COMMON.MCE' -! include 'COMMON.IOUNITS' -! include 'COMMON.VAR' -#ifdef MPL - use MPI_data !include 'COMMON.INFO' -#endif -! include 'COMMON.GEO' - real(kind=8) :: ecur,eold,xx,bol !,ran_number - real(kind=8),dimension(6*nres) :: x,xold !(maxvar) (maxvar=6*maxres) - real(kind=8) :: tole=1.0D-1, tola=5.0D0 - logical :: accepted - -!!! local variables - el - integer :: indecur - real(kind=8) :: scur,sold,xxh,deix,dent - -! Check if the conformation is similar. -!d write (iout,*) 'Enter ACCEPTING' -!d write (iout,*) 'Old PHI angles:' -!d write (iout,*) (rad2deg*xold(i),i=1,nphi) -!d write (iout,*) 'Current angles' -!d write (iout,*) (rad2deg*x(i),i=1,nphi) -!d ddif=dif_ang(nphi,x,xold) -!d write (iout,*) 'Angle norm:',ddif -!d 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 -! Else evaluate the entropy of the conf and compare it with that of the previous -! 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 -!d print *,'Processor',MyID,' ecur=',ecur,' indecur=',indecur, -!d & ' scur=',scur,' eold=',eold,' sold=',sold -!d 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 -! 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 subroutine accepting -!----------------------------------------------------------------------------- - subroutine read_pool - - use io_base, only:read_angles -! implicit real*8 (a-h,o-z) -! include 'DIMENSIONS' -! include 'COMMON.IOUNITS' -! include 'COMMON.GEO' -! include 'COMMON.MCM' -! include 'COMMON.MCE' -! include 'COMMON.VAR' - real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres) - -!!! local variables - el - integer :: j,i,iconf - - 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 subroutine read_pool -!----------------------------------------------------------------------------- -! mc.F -!----------------------------------------------------------------------------- - subroutine monte_carlo - - use energy_data - use geometry_data - use MPI_data, only:ifinish,nctasks,WhatsUp,MyID - use control_data, only:refstr,MaxProcs - use io_base - use control, only:tcpu,ovrtim - use regularize_, only:fitsq - use compare -! use control -! Does Boltzmann and entropic sampling without energy minimization -! implicit real*8 (a-h,o-z) -! include 'DIMENSIONS' -! include 'COMMON.IOUNITS' -#ifdef MPL - use MPI_data !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,error,lprint !,ovrtim - integer :: MoveType,nbond,nbins -! integer :: conf_comp - real(kind=8) :: RandOrPert - real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres) - real(kind=8) :: elowest,elowest1,ehighest,ehighest1,eold - real(kind=8) :: przes(3),obr(3,3) - real(kind=8),dimension(6*nres) :: varold !(maxvar) (maxvar=6*maxres) - logical :: non_conv - integer,dimension(-1:MaxMoveType+1,0:MaxProcs-1) :: moves1,moves_acc1 !(-1:MaxMoveType+1,0:MaxProcs-1) -#ifdef MPL - real(kind=8) :: etot_temp,etot_all(0:MaxProcs) - external d_vadd,d_vmin,d_vmax - real(kind=8),dimension(-max_ene:max_ene) :: entropy1,nhist1 - integer,dimension(nres*(MaxProcs+1)) :: nbond_move1,nbond_acc1 - integer,dimension(2) :: itemp -#endif - real(kind=8),dimension(6*nres) :: var_lowest !(maxvar) (maxvar=6*maxres) - real(kind=8),dimension(0:n_ene) :: energia,energia_ave -! -!!! local variables - el - integer :: i,j,it,ii,iproc,nacc,ISWEEP,nfun,indmax,indmin,ijunk,& - Kwita,indeold,imdmax,inde,iretcode,nstart_grow,noverlap - real(kind=8) :: facee,conste,ejunk,etot,sold,frac,runtime,& - frac_ave,rms_ave,etot_ave,scur,from_pool,co,rms - - 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) -! Number of bins in energy histogram - nbins=e_up/delte-1 - write (iout,*) 'NBINS=',nbins - conste=dlog(facee) -! 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 -! Read the pool of conformations - call read_pool - elowest=1.0D+10 - ehighest=-1.0D+10 -!---------------------------------------------------------------------------- -! Entropy-sampling simulations with continually updated entropy; -! set NSWEEP=1 for Boltzmann sampling. -! Loop thru simulations -!---------------------------------------------------------------------------- - allocate(ifinish(nctasks)) - DO ISWEEP=1,NSWEEP -! -! Initialize the IFINISH array. -! -#ifdef MPL - do i=1,nctasks - ifinish(i)=0 - enddo -#endif -!--------------------------------------------------------------------------- -! Initialize counters. -!--------------------------------------------------------------------------- -! Total number of generated confs. - ngen=0 -! Total number of moves. In general this won't be equal to the number of -! attempted moves, because we may want to reject some "bad" confs just by -! overlap check. - nmove=0 -! Total number of shift (nbond_move(1)), spike, crankshaft, three-bond,... -! motions. -!el allocate(nbond_move(nres)) !(maxres) -!el allocate(nbond_acc(nres)) !(maxres) - - do i=1,nres - nbond_move(i)=0 - nbond_acc(i)=0 - enddo -! Initialize total and accepted number of moves of various kind. - do i=-1,MaxMoveType - moves(i)=0 - moves_acc(i)=0 - enddo -! Total number of energy evaluations. - neneval=0 - nfun=0 -!---------------------------------------------------------------------------- -! Take a conformation from the pool -!---------------------------------------------------------------------------- - 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) -! 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 -!---------------------------------------------------------------------------- -! Compute and print initial energies. -!---------------------------------------------------------------------------- - 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) - etot = energia(0) - call enerprint(energia) - if (refstr) then - call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,1),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 -! close(istat) - neneval=neneval+1 - if (.not. ent_read) then -! Initialize the entropy array -#ifdef MPL -! 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 -! 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 -! Multicanonical sampling - start from Boltzmann distribution - do i=-max_ene,max_ene - entropy(i)=(emin+i*delte)*betbol - enddo - ELSE -! 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) -! NACC is the counter for the accepted conformations of a given processor - nacc=0 -! NACC_TOT counts the total number of accepted conformations - nacc_tot=0 -! Main loop. -!---------------------------------------------------------------------------- -! Zero out average energies - do i=0,n_ene - energia_ave(i)=0.0d0 - enddo -! Initialize energy histogram - do i=-max_ene,max_ene - nhist(i)=0.0D0 - enddo ! i -! Zero out iteration counter. - it=0 - do j=1,nvar - varold(j)=varia(j) - enddo -! Begin MC iteration loop. - do while (not_done) - it=it+1 -! 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 -! Retrieve the angles of previously accepted conformation - do j=1,nvar - varia(j)=varold(j) - enddo - call var_to_geom(nvar,varia) -! Rebuild the chain. - call chainbuild - MoveType=0 - nbond=0 - lprint=.true. -! Decide whether to take a conformation from the pool or generate/perturb one -! randomly - from_pool=ran_number(0.0D0,1.0D0) - if (npool.gt.0 .and. from_pool.lt.pool_fraction) then -! 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 -!d call intout -!d 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 -! 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 -!d 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) - etot = energia(0) - neneval=neneval+1 -!d call enerprint(energia(0)) -!d 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 -!-------------------------------------------------------------------------- -!... Acceptance test -!-------------------------------------------------------------------------- - 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 -! Compare with reference structure. - if (refstr) then - call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,1),& - 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 -! -! Periodically save average energies and confs. -! - 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) -!#ifdef AIX -! open (istat,file=statname,position='append') -!#else -! open (istat,file=statname,access='append') -!endif - 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 -! close(istat) - if (print_mc.gt.0) & - call statprint(nacc,nfun,iretcode,etot,elowest) -! 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) -! 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 -! 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 -! if ((ntrial.gt.maxtrial_iter -! & .or. (it/pool_read_freq)*pool_read_freq.eq.it) -! & .and. npool.gt.0) then -! Take a conformation from the pool -! ii=iran_num(1,npool) -! do i=1,nvar -! varold(i)=xpool(i,ii) -! enddo -! if (ntrial.gt.maxtrial_iter) -! & 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*varold(i),i=1,nvar) -! ntrial=0 -! eold=epool(ii) -! call entropia(eold,sold,indeold) -! accepted=.true. -! 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 -!d write (iout,*) 'WhatsUp=',WhatsUp,' Kwita=',Kwita - not_done = (nacc_tot.lt.maxacc) .and. (WhatsUp.eq.0) & - .and. (Kwita.eq.0) -!d write (iout,*) 'not_done=',not_done -#ifdef MPL - if (Kwita.lt.0) then - print *,'Processor',MyID,& - ' has received STOP signal =',Kwita,' in EntSamp.' -!d 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.' -!d 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.' -!d print *,'not_done=',not_done - call send_stop_sig(-2) - if (MyID.ne.MasterID) call send_MCM_info(-1) - endif -#endif - enddo ! not_done - -!----------------------------------------------------------------- -!... Construct energy histogram & update entropy -!----------------------------------------------------------------- - 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 -! Wait until every processor has sent complete MC info. - if (MyID.eq.MasterID) then - not_done=.true. - do while (not_done) -! write (*,*) 'The IFINISH array:' -! 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 -! 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 -! 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 -! 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 -!----------------------------------------------------------------- -!... End of energy histogram construction -!----------------------------------------------------------------- -#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 - -!--------------------------------------------------------------------------- - ENDDO ! ISWEEP -!--------------------------------------------------------------------------- - - runtime=tcpu() - - if (isweep.eq.nsweep .and. it.ge.maxacc) & - write (iout,'(/80(1h*)/20x,a/80(1h*)/)') 'All iterations done.' - return - end subroutine monte_carlo -!----------------------------------------------------------------------------- - 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 - use MPI_data !include 'COMMON.INFO' -#endif -! include 'COMMON.GEO' - real(kind=8) :: ecur,eold,xx,bol - real(kind=8),dimension(6*nres) :: x,xold !(maxvar) (maxvar=6*maxres) - logical :: accepted - -!el local variables - integer :: it,indecur - real(kind=8) :: scur,sold,xxh -! Check if the conformation is similar. -!d write (iout,*) 'Enter ACCEPTING' -!d write (iout,*) 'Old PHI angles:' -!d write (iout,*) (rad2deg*xold(i),i=1,nphi) -!d write (iout,*) 'Current angles' -!d write (iout,*) (rad2deg*x(i),i=1,nphi) -!d ddif=dif_ang(nphi,x,xold) -!d write (iout,*) 'Angle norm:',ddif -!d 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 -! Else evaluate the entropy of the conf and compare it with that of the previous -! one. - call entropia(ecur,scur,indecur) -!d print *,'Processor',MyID,' ecur=',ecur,' indecur=',indecur, -!d & ' scur=',scur,' eold=',eold,' sold=',sold -!d 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 -! 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 subroutine accept_mc -!----------------------------------------------------------------------------- - integer function icialosc(x) - - real(kind=8) :: x - if (x.lt.0.0D0) then - icialosc=dint(x)-1 - else - icialosc=dint(x) - endif - return - end function icialosc -!----------------------------------------------------------------------------- - subroutine entropia(ecur,scur,indecur) - - use energy_data, only: max_ene -! implicit real*8 (a-h,o-z) -! include 'DIMENSIONS' -! include 'COMMON.MCM' -! include 'COMMON.MCE' -! include 'COMMON.IOUNITS' - real(kind=8) :: ecur,scur,deix,dent - integer :: indecur,it !???el - - 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 subroutine entropia -!----------------------------------------------------------------------------- -! mcm.F -!----------------------------------------------------------------------------- - subroutine mcm_setup - - use energy_data -! implicit real*8 (a-h,o-z) -! include 'DIMENSIONS' -! include 'COMMON.IOUNITS' -! include 'COMMON.MCM' -! include 'COMMON.CONTROL' -! include 'COMMON.INTERACT' -! include 'COMMON.NAMES' -! include 'COMMON.CHAIN' -! include 'COMMON.VAR' -! -!!! local variables - el - integer :: i,i1,i2,it1,it2,ngly,mmm,maxwinlen - -! Set up variables used in MC/MCM. -! -! allocate(sumpro_bond(0:nres)) !(0:maxres) - - write (iout,'(80(1h*)/20x,a/80(1h*))') 'MCM control parameters:' - write (iout,'(5(a,i7))') 'Maxacc:',maxacc,' MaxTrial:',MaxTrial,& - ' MaxRepm:',MaxRepm,' MaxGen:',MaxGen,' MaxOverlap:',MaxOverlap - write (iout,'(4(a,f8.1)/2(a,i3))') & - 'Tmin:',Tmin,' Tmax:',Tmax,' TstepH:',TstepH,& - ' TstepC:',TstepC,'NstepH:',NstepH,' NstepC:',NstepC - if (nwindow.gt.0) then - write (iout,'(a)') 'Perturbation windows:' - do i=1,nwindow - i1=winstart(i) - i2=winend(i) - it1=itype(i1) - it2=itype(i2) - write (iout,'(a,i3,a,i3,a,i3)') restyp(it1),i1,restyp(it2),i2,& - ' length',winlen(i) - enddo - endif -! Rbolt=8.3143D-3*2.388459D-01 kcal/(mol*K) - RBol=1.9858D-3 -! Number of "end bonds". - koniecl=0 -! koniecl=nphi - print *,'koniecl=',koniecl - write (iout,'(a)') 'Probabilities of move types:' - write (*,'(a)') 'Probabilities of move types:' - do i=1,MaxMoveType - write (iout,'(a,f10.5)') MovTypID(i),& - sumpro_type(i)-sumpro_type(i-1) - write (*,'(a,f10.5)') MovTypID(i),& - sumpro_type(i)-sumpro_type(i-1) - enddo - write (iout,*) -! Maximum length of N-bond segment to be moved -! nbm=nres-1-(2*koniecl-1) - if (nwindow.gt.0) then - maxwinlen=winlen(1) - do i=2,nwindow - if (winlen(i).gt.maxwinlen) maxwinlen=winlen(i) - enddo - nbm=min0(maxwinlen,6) - write (iout,'(a,i3,a,i3)') 'Nbm=',Nbm,' Maxwinlen=',Maxwinlen - else - nbm=min0(6,nres-2) - endif - sumpro_bond(0)=0.0D0 - sumpro_bond(1)=0.0D0 - do i=2,nbm - sumpro_bond(i)=sumpro_bond(i-1)+1.0D0/dfloat(i) - enddo - write (iout,'(a)') 'The SumPro_Bond array:' - write (iout,'(8f10.5)') (sumpro_bond(i),i=1,nbm) - write (*,'(a)') 'The SumPro_Bond array:' - write (*,'(8f10.5)') (sumpro_bond(i),i=1,nbm) -! Maximum number of side chains moved simultaneously -! print *,'nnt=',nnt,' nct=',nct - ngly=0 - do i=nnt,nct - if (itype(i).eq.10) ngly=ngly+1 - enddo - mmm=nct-nnt-ngly+1 - if (mmm.gt.0) then - MaxSideMove=min0((nct-nnt+1)/2,mmm) - endif -! print *,'MaxSideMove=',MaxSideMove -! Max. number of generated confs (not used at present). - maxgen=10000 -! Set initial temperature - Tcur=Tmin - betbol=1.0D0/(Rbol*Tcur) - write (iout,'(a,f8.1,a,f10.5)') 'Initial temperature:',Tcur,& - ' BetBol:',betbol - write (iout,*) 'RanFract=',ranfract - return - end subroutine mcm_setup -!----------------------------------------------------------------------------- -#ifndef MPI - subroutine do_mcm(i_orig) - - use geometry_data - use energy_data - use MPI_data, only:Whatsup - use control_data, only:refstr,minim,iprint - use io_base - use control, only:tcpu,ovrtim - use regularize_, only:fitsq - use compare - use minimm, only:minimize -! Monte-Carlo-with-Minimization calculations - serial code. -! implicit real*8 (a-h,o-z) -! include 'DIMENSIONS' -! include 'COMMON.IOUNITS' -! include 'COMMON.GEO' -! include 'COMMON.CHAIN' -! include 'COMMON.MCM' -! include 'COMMON.CONTACTS' -! include 'COMMON.CONTROL' -! include 'COMMON.VAR' -! include 'COMMON.INTERACT' -! include 'COMMON.CACHE' -!rc include 'COMMON.DEFORM' -!rc include 'COMMON.DEFORM1' -! include 'COMMON.NAMES' - logical :: accepted,over,error,lprint,not_done,my_conf,& - enelower,non_conv !,ovrtim - integer :: MoveType,nbond !,conf_comp - integer,dimension(max_cache) :: ifeed - real(kind=8),dimension(6*nres) :: varia,varold !(maxvar) (maxvar=6*maxres) - real(kind=8) :: elowest,eold - real(kind=8) :: przes(3),obr(3,3) - real(kind=8) :: energia(0:n_ene) - real(kind=8) :: coord1(6*nres,max_thread2),enetb1(max_threadss) !el -!!! local variables - el - integer :: i,nf,nacc,it,nout,j,i_orig,nfun,Kwita,iretcode,& - noverlap,nstart_grow,irepet,n_thr,ii - real(kind=8) :: etot,frac,rms,co,RandOrPert,& - rms_deform,runtime -!--------------------------------------------------------------------------- -! Initialize counters. -!--------------------------------------------------------------------------- -! Total number of generated confs. - ngen=0 -! Total number of moves. In general this won't be equal to the number of -! attempted moves, because we may want to reject some "bad" confs just by -! overlap check. - nmove=0 -! Total number of temperature jumps. - ntherm=0 -! Total number of shift (nbond_move(1)), spike, crankshaft, three-bond,... -! motions. -! if(.not.allocated(varsave)) allocate(varsave(maxvar,maxsave)) !(maxvar,maxsave) -! allocate(nbond_move(nres)) !(maxres) - - ncache=0 - do i=1,nres - nbond_move(i)=0 - enddo -! Initialize total and accepted number of moves of various kind. - do i=0,MaxMoveType - moves(i)=0 - moves_acc(i)=0 - enddo -! Total number of energy evaluations. - neneval=0 - nfun=0 - nsave=0 - - write (iout,*) 'RanFract=',RanFract - - WhatsUp=0 - Kwita=0 - -!---------------------------------------------------------------------------- -! Compute and print initial energies. -!---------------------------------------------------------------------------- - call intout - write (iout,'(/80(1h*)/a)') 'Initial energies:' - call chainbuild - nf=0 - - call etotal(energia) - etot = energia(0) -! Minimize the energy of the first conformation. - if (minim) then - call geom_to_var(nvar,varia) -! write (iout,*) 'The VARIA array' -! write (iout,'(8f10.4)') (rad2deg*varia(i),i=1,nvar) - call minimize(etot,varia,iretcode,nfun) - call var_to_geom(nvar,varia) - call chainbuild - write (iout,*) 'etot from MINIMIZE:',etot -! write (iout,*) 'Tha VARIA array' -! write (iout,'(8f10.4)') (rad2deg*varia(i),i=1,nvar) - - call etotal(energia) - etot=energia(0) - call enerprint(energia) - endif - if (refstr) then - call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,1),nsup,przes,& !el cref(1,nstart_sup) - 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 - if (print_stat) & - write (istat,'(i5,17(1pe14.5))') 0,& - (energia(print_order(i)),i=1,nprint_ene),& - etot,rms,frac,co - else - if (print_stat) write (istat,'(i5,16(1pe14.5))') 0,& - (energia(print_order(i)),i=1,nprint_ene),etot - endif - if (print_stat) close(istat) - neneval=neneval+nfun+1 - write (iout,'(/80(1h*)/20x,a/80(1h*))') & - 'Enter Monte Carlo procedure.' - if (print_int) then - close(igeom) - call briefout(0,etot) - endif - eold=etot - do i=1,nvar - varold(i)=varia(i) - enddo - elowest=etot - call zapis(varia,etot) - nacc=0 ! total # of accepted confs of the current processor. - nacc_tot=0 ! total # of accepted confs of all processors. - - not_done = (iretcode.ne.11) - -!---------------------------------------------------------------------------- -! Main loop. -!---------------------------------------------------------------------------- - it=0 - nout=0 - do while (not_done) - it=it+1 - write (iout,'(80(1h*)/20x,a,i7)') & - 'Beginning iteration #',it -! Initialize local counter. - ntrial=0 ! # of generated non-overlapping confs. - accepted=.false. - do while (.not. accepted) - -! Retrieve the angles of previously accepted conformation - noverlap=0 ! # of overlapping confs. - do j=1,nvar - varia(j)=varold(j) - enddo - call var_to_geom(nvar,varia) -! Rebuild the chain. - call chainbuild -! Heat up the system, if necessary. - call heat(over) -! If temperature cannot be further increased, stop. - if (over) goto 20 - MoveType=0 - nbond=0 - lprint=.true. -!d write (iout,'(a)') 'Old variables:' -!d write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar) -! 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) -!d write (iout,'(a)') 'New variables:' -!d write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar) - ngen=ngen+1 - - call etotal(energia) - etot=energia(0) -! call enerprint(energia(0)) -! write (iout,'(2(a,1pe14.5))') 'Etot=',Etot,' Elowest=',Elowest - if (etot-elowest.gt.overlap_cut) then - if(iprint.gt.1.or.etot.lt.1d20) & - write (iout,'(a,1pe14.5)') & - '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) -!d write (iout,*) 'etot from MINIMIZE:',etot -!d write (iout,'(a)') 'Variables after minimization:' -!d write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar) - - call etotal(energia) - etot = energia(0) - neneval=neneval+nfun+2 - endif -! call enerprint(energia(0)) - write (iout,'(a,i6,a,1pe16.6)') 'Conformation:',ngen,& - ' energy:',etot -!-------------------------------------------------------------------------- -!... Do Metropolis test -!-------------------------------------------------------------------------- - accepted=.false. - my_conf=.false. - - if (WhatsUp.eq.0 .and. Kwita.eq.0) then - call metropolis(nvar,varia,varold,etot,eold,accepted,& - my_conf,EneLower) - endif - write (iout,*) 'My_Conf=',My_Conf,' EneLower=',EneLower - if (accepted) then - - nacc=nacc+1 - nacc_tot=nacc_tot+1 - if (elowest.gt.etot) elowest=etot - moves_acc(MoveType)=moves_acc(MoveType)+1 - if (MoveType.eq.1) then - nbond_acc(nbond)=nbond_acc(nbond)+1 - endif -! Check against conformation repetitions. - irepet=conf_comp(varia,etot) - if (print_stat) then -#if defined(AIX) || defined(PGI) - open (istat,file=statname,position='append') -#else - open (istat,file=statname,access='append') -#endif - endif - call statprint(nacc,nfun,iretcode,etot,elowest) - if (refstr) then - call var_to_geom(nvar,varia) - call chainbuild - call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,1),& !el 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)') & - 'RMS deviation from the reference structure:',rms,& - ' % of native contacts:',frac*100,' contact order',co - endif ! refstr - if (My_Conf) then - nout=nout+1 - write (iout,*) 'Writing new conformation',nout - if (refstr) then - write (istat,'(i5,16(1pe14.5))') nout,& - (energia(print_order(i)),i=1,nprint_ene),& - etot,rms,frac - else - if (print_stat) & - write (istat,'(i5,17(1pe14.5))') nout,& - (energia(print_order(i)),i=1,nprint_ene),etot - endif ! refstr - if (print_stat) close(istat) -! Print internal coordinates. - if (print_int) call briefout(nout,etot) -! Accumulate the newly accepted conf in the coord1 array, if it is different -! from all confs that are already there. - call compare_s1(n_thr,max_thread2,etot,varia,ii,& - enetb1,coord1,rms_deform,.true.,iprint) - write (iout,*) 'After compare_ss: n_thr=',n_thr - if (ii.eq.1 .or. ii.eq.3) then - write (iout,'(8f10.4)') & - (rad2deg*coord1(i,n_thr),i=1,nvar) - endif - else - write (iout,*) 'Conformation from cache, not written.' - endif ! My_Conf - - if (nrepm.gt.maxrepm) then - write (iout,'(a)') 'Too many conformation repetitions.' - goto 20 - endif -! Store the accepted conf. and its energy. - eold=etot - do i=1,nvar - varold(i)=varia(i) - enddo - if (irepet.eq.0) call zapis(varia,etot) -! Lower the temperature, if necessary. - call cool - - else - - ntrial=ntrial+1 - endif ! accepted - endif ! overlap - - 30 continue - enddo ! accepted -! Check for time limit. - if (ovrtim()) WhatsUp=-1 - not_done = (nacc_tot.lt.maxacc) .and. (WhatsUp.eq.0) & - .and. (Kwita.eq.0) - - enddo ! not_done - goto 21 - 20 WhatsUp=-3 - - 21 continue - runtime=tcpu() - 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) - if (it.ge.maxacc) & - write (iout,'(/80(1h*)/20x,a/80(1h*)/)') 'All iterations done.' - !(maxvar) (maxvar=6*maxres) - return - end subroutine do_mcm -#endif -!----------------------------------------------------------------------------- -#ifdef MPI - subroutine do_mcm(i_orig) - -! Monte-Carlo-with-Minimization calculations - parallel code. - use MPI_data - use control_data, only:refstr!,tag - use io_base, only:intout,briefout - use control, only:ovrtim,tcpu - use compare, only:contact,contact_fract - use minimm, only:minimize - use regularize_, only:fitsq -! use contact_, only:contact -! use minim -! implicit real*8 (a-h,o-z) -! include 'DIMENSIONS' - include 'mpif.h' -! include 'COMMON.IOUNITS' -! include 'COMMON.GEO' -! include 'COMMON.CHAIN' -! include 'COMMON.MCM' -! include 'COMMON.CONTACTS' -! include 'COMMON.CONTROL' -! include 'COMMON.VAR' -! include 'COMMON.INTERACT' -! include 'COMMON.INFO' -! include 'COMMON.CACHE' -!rc include 'COMMON.DEFORM' -!rc include 'COMMON.DEFORM1' -!rc include 'COMMON.DEFORM2' -! include 'COMMON.MINIM' -! include 'COMMON.NAMES' - logical :: accepted,over,error,lprint,not_done,similar,& - enelower,non_conv,flag,finish !,ovrtim - integer :: MoveType,nbond !,conf_comp - real(kind=8),dimension(6*nres) :: x1,varold1,varold,varia !(maxvar) (maxvar=6*maxres) - real(kind=8) :: elowest,eold - real(kind=8) :: przes(3),obr(3,3) - integer :: iparentx(max_threadss2) - integer :: iparentx1(max_threadss2) - integer :: imtasks(150),imtasks_n - real(kind=8) :: energia(0:n_ene) - -!el local variables - integer :: nfun,nodenum,i_orig,i,nf,nacc,it,nout,j,kkk,is,& - Kwita,iretcode,noverlap,nstart_grow,ierr,iitt,& - ii_grnum_d,ii_ennum_d,ii_hesnum_d,i_grnum_d,i_ennum_d,& - i_hesnum_d,i_minimiz,irepet - real(kind=8) :: etot,frac,eneglobal,RandOrPert,eold1,co,& - runtime,rms - -! if(.not.allocated(varsave)) allocate(varsave(maxvar,maxsave)) !(maxvar,maxsave) - print *,'Master entered DO_MCM' - nodenum = nprocs - - finish=.false. - imtasks_n=0 - do i=1,nodenum-1 - imtasks(i)=0 - enddo -!--------------------------------------------------------------------------- -! Initialize counters. -!--------------------------------------------------------------------------- -! Total number of generated confs. - ngen=0 -! Total number of moves. In general this won`t be equal to the number of -! attempted moves, because we may want to reject some "bad" confs just by -! overlap check. - nmove=0 -! Total number of temperature jumps. - ntherm=0 -! Total number of shift (nbond_move(1)), spike, crankshaft, three-bond,... -! motions. - allocate(nbond_move(nres)) !(maxres) - - ncache=0 - do i=1,nres - nbond_move(i)=0 - enddo -! Initialize total and accepted number of moves of various kind. - do i=0,MaxMoveType - moves(i)=0 - moves_acc(i)=0 - enddo -! Total number of energy evaluations. - neneval=0 - nfun=0 - nsave=0 -! write (iout,*) 'RanFract=',RanFract - WhatsUp=0 - Kwita=0 -!---------------------------------------------------------------------------- -! Compute and print initial energies. -!---------------------------------------------------------------------------- - call intout - write (iout,'(/80(1h*)/a)') 'Initial energies:' - call chainbuild - nf=0 - call etotal(energia) - etot = energia(0) - call enerprint(energia) -! Request energy computation from slave processors. - call geom_to_var(nvar,varia) -! write (iout,*) 'The VARIA array' -! write (iout,'(8f10.4)') (rad2deg*varia(i),i=1,nvar) - call minimize(etot,varia,iretcode,nfun) - call var_to_geom(nvar,varia) - call chainbuild - write (iout,*) 'etot from MINIMIZE:',etot -! write (iout,*) 'Tha VARIA array' -! write (iout,'(8f10.4)') (rad2deg*varia(i),i=1,nvar) - neneval=0 - eneglobal=1.0d99 - if (print_mc .gt. 0) write (iout,'(/80(1h*)/20x,a/80(1h*))') & - 'Enter Monte Carlo procedure.' - if (print_mc .gt. 0) write (iout,'(i5,1pe14.5)' ) i_orig,etot - eold=etot - do i=1,nvar - varold(i)=varia(i) - enddo - elowest=etot - call zapis(varia,etot) -! diagnostics - call var_to_geom(nvar,varia) - call chainbuild - call etotal(energia) - if (print_mc.gt.0) write (iout,*) 'Initial energy:',etot -! end diagnostics - nacc=0 ! total # of accepted confs of the current processor. - nacc_tot=0 ! total # of accepted confs of all processors. - not_done=.true. -!---------------------------------------------------------------------------- -! Main loop. -!---------------------------------------------------------------------------- - it=0 - nout=0 - LOOP1:do while (not_done) - it=it+1 - if (print_mc.gt.0) write (iout,'(80(1h*)/20x,a,i7)') & - 'Beginning iteration #',it -! Initialize local counter. - ntrial=0 ! # of generated non-overlapping confs. - noverlap=0 ! # of overlapping confs. - accepted=.false. - LOOP2:do while (.not. accepted) - - LOOP3:do while (imtasks_n.lt.nodenum-1.and..not.finish) - do i=1,nodenum-1 - if(imtasks(i).eq.0) then - is=i - exit - endif - enddo -! Retrieve the angles of previously accepted conformation - do j=1,nvar - varia(j)=varold(j) - enddo - call var_to_geom(nvar,varia) -! Rebuild the chain. - call chainbuild -! Heat up the system, if necessary. - call heat(over) -! If temperature cannot be further increased, stop. - if (over) then - finish=.true. - endif - MoveType=0 - nbond=0 -! write (iout,'(a)') 'Old variables:' -! write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar) -! 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) -! print *,'after perturb',error,finish - if (error) finish = .true. - if (MoveType.lt.1 .or. MoveType.gt.MaxMoveType) then - write (iout,'(/a,i7,a/)') 'Error - unknown MoveType=',& - MoveType,' returned from PERTURB.' - finish=.true. - write (*,'(/a,i7,a/)') 'Error - unknown MoveType=',& - MoveType,' returned from PERTURB.' - 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) - ngen=ngen+1 -! print *,'finish=',finish - if (etot-elowest.gt.overlap_cut) then - if (print_mc.gt.1) write (iout,'(a,1pe14.5)') & - 'Overlap detected in the current conf.; energy is',etot - if(iprint.gt.1.or.etot.lt.1d19) print *,& - 'Overlap detected in the current conf.; energy is',etot - neneval=neneval+1 - accepted=.false. - noverlap=noverlap+1 - if (noverlap.gt.maxoverlap) then - write (iout,*) 'Too many overlapping confs.',& - ' etot, elowest, overlap_cut', etot, elowest, overlap_cut - finish=.true. - endif - else if (.not. finish) then -! Distribute tasks to processors -! print *,'Master sending order' - call MPI_SEND(12, 1, MPI_INTEGER, is, tag,& - CG_COMM, ierr) -! write (iout,*) '12: tag=',tag -! print *,'Master sent order to processor',is - call MPI_SEND(it, 1, MPI_INTEGER, is, tag,& - CG_COMM, ierr) -! write (iout,*) 'it: tag=',tag - call MPI_SEND(eold, 1, MPI_DOUBLE_PRECISION, is, tag,& - CG_COMM, ierr) -! write (iout,*) 'eold: tag=',tag - call MPI_SEND(varia(1), nvar, MPI_DOUBLE_PRECISION, & - is, tag,& - CG_COMM, ierr) -! write (iout,*) 'varia: tag=',tag - call MPI_SEND(varold(1), nvar, MPI_DOUBLE_PRECISION, & - is, tag,& - CG_COMM, ierr) -! write (iout,*) 'varold: tag=',tag -#ifdef AIX - call flush_(iout) -#else - call flush(iout) -#endif - imtasks(is)=1 - imtasks_n=imtasks_n+1 -! End distribution - endif ! overlap - enddo LOOP3 - - flag = .false. - LOOP_RECV:do while(.not.flag) - do is=1, nodenum-1 - call MPI_IPROBE(is,tag,CG_COMM,flag,status,ierr) - if(flag) then - call MPI_RECV(iitt, 1, MPI_INTEGER, is, tag,& - CG_COMM, status, ierr) - call MPI_RECV(eold1, 1, MPI_DOUBLE_PRECISION, is, tag,& - CG_COMM, status, ierr) - call MPI_RECV(etot, 1, MPI_DOUBLE_PRECISION, is, tag,& - CG_COMM, status, ierr) - call MPI_RECV(varia(1), nvar, MPI_DOUBLE_PRECISION,is,tag,& - CG_COMM, status, ierr) - call MPI_RECV(varold1(1), nvar, MPI_DOUBLE_PRECISION, is, & - tag, CG_COMM, status, ierr) - call MPI_RECV(ii_grnum_d, 1, MPI_INTEGER, is, tag,& - CG_COMM, status, ierr) - call MPI_RECV(ii_ennum_d, 1, MPI_INTEGER, is, tag,& - CG_COMM, status, ierr) - call MPI_RECV(ii_hesnum_d, 1, MPI_INTEGER, is, tag,& - CG_COMM, status, ierr) - i_grnum_d=i_grnum_d+ii_grnum_d - i_ennum_d=i_ennum_d+ii_ennum_d - neneval = neneval+ii_ennum_d - i_hesnum_d=i_hesnum_d+ii_hesnum_d - i_minimiz=i_minimiz+1 - imtasks(is)=0 - imtasks_n=imtasks_n-1 - exit - endif - enddo - enddo LOOP_RECV - - if(print_mc.gt.0) write (iout,'(a,i6,a,i6,a,i6,a,1pe16.6)') & - 'From Worker #',is,' iitt',iitt,& - ' Conformation:',ngen,' energy:',etot -!-------------------------------------------------------------------------- -!... Do Metropolis test -!-------------------------------------------------------------------------- - call metropolis(nvar,varia,varold1,etot,eold1,accepted,& - similar,EneLower) - if(iitt.ne.it.and..not.similar) then - call metropolis(nvar,varia,varold,etot,eold,accepted,& - similar,EneLower) - accepted=enelower - endif - if(etot.lt.eneglobal)eneglobal=etot -! if(mod(it,100).eq.0) - write(iout,*)'CHUJOJEB ',neneval,eneglobal - if (accepted) then -! Write the accepted conformation. - nout=nout+1 - if (refstr) then - call var_to_geom(nvar,varia) - call chainbuild - 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 - endif ! refstr - if (print_mc.gt.0) & - write (iout,*) 'Writing new conformation',nout - if (print_stat) then - call var_to_geom(nvar,varia) -#if defined(AIX) || defined(PGI) - open (istat,file=statname,position='append') -#else - open (istat,file=statname,access='append') -#endif - if (refstr) then - write (istat,'(i5,16(1pe14.5))') nout,& - (energia(print_order(i)),i=1,nprint_ene),& - etot,rms,frac - else - write (istat,'(i5,16(1pe14.5))') nout,& - (energia(print_order(i)),i=1,nprint_ene),etot - endif ! refstr - close(istat) - endif ! print_stat -! Print internal coordinates. - if (print_int) call briefout(nout,etot) - nacc=nacc+1 - nacc_tot=nacc_tot+1 - if (elowest.gt.etot) elowest=etot - moves_acc(MoveType)=moves_acc(MoveType)+1 - if (MoveType.eq.1) then - nbond_acc(nbond)=nbond_acc(nbond)+1 - endif -! Check against conformation repetitions. - irepet=conf_comp(varia,etot) - if (nrepm.gt.maxrepm) then - if (print_mc.gt.0) & - write (iout,'(a)') 'Too many conformation repetitions.' - finish=.true. - endif -! Store the accepted conf. and its energy. - eold=etot - do i=1,nvar - varold(i)=varia(i) - enddo - if (irepet.eq.0) call zapis(varia,etot) -! Lower the temperature, if necessary. - call cool - else - ntrial=ntrial+1 - endif ! accepted - 30 continue - if(finish.and.imtasks_n.eq.0)exit LOOP2 - enddo LOOP2 ! accepted -! Check for time limit. - not_done = (it.lt.max_mcm_it) .and. (nacc_tot.lt.maxacc) - if(.not.not_done .or. finish) then - if(imtasks_n.gt.0) then - not_done=.true. - else - not_done=.false. - endif - finish=.true. - endif - enddo LOOP1 ! not_done - runtime=tcpu() - if (print_mc.gt.0) then - 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) - if (it.ge.maxacc) & - write (iout,'(/80(1h*)/20x,a/80(1h*)/)') 'All iterations done.' - endif -#ifdef AIX - call flush_(iout) -#else - call flush(iout) -#endif - do is=1,nodenum-1 - call MPI_SEND(999, 1, MPI_INTEGER, is, tag,& - CG_COMM, ierr) - enddo - return - end subroutine do_mcm -!----------------------------------------------------------------------------- - subroutine execute_slave(nodeinfo,iprint) - - use MPI_data - use minimm, only:minimize -! use minim -! implicit real*8 (a-h,o-z) -! include 'DIMENSIONS' - include 'mpif.h' -! include 'COMMON.TIME1' -! include 'COMMON.IOUNITS' -!rc include 'COMMON.DEFORM' -!rc include 'COMMON.DEFORM1' -!rc include 'COMMON.DEFORM2' -! include 'COMMON.LOCAL' -! include 'COMMON.VAR' -! include 'COMMON.INFO' -! include 'COMMON.MINIM' - character(len=10) :: nodeinfo - real(kind=8),dimension(6*nres) :: x,x1 !(maxvar) (maxvar=6*maxres) - integer :: nfun,iprint,i_switch,ierr,i_grnum_d,i_ennum_d,& - i_hesnum_d,iitt,iretcode,iminrep - real(kind=8) :: ener,energyx - - nodeinfo='chujwdupe' -! print *,'Processor:',MyID,' Entering execute_slave' - tag=0 -! call MPI_SEND(nodeinfo, 10, MPI_CHARACTER, 0, tag, -! & CG_COMM, ierr) - -1001 call MPI_RECV(i_switch, 1, MPI_INTEGER, 0, tag,& - CG_COMM, status, ierr) -! write(iout,*)'12: tag=',tag - if(iprint.ge.2)print *, MyID,' recv order ',i_switch - if (i_switch.eq.12) then - i_grnum_d=0 - i_ennum_d=0 - i_hesnum_d=0 - call MPI_RECV(iitt, 1, MPI_INTEGER, 0, tag,& - CG_COMM, status, ierr) -! write(iout,*)'12: tag=',tag - call MPI_RECV(ener, 1, MPI_DOUBLE_PRECISION, 0, tag,& - CG_COMM, status, ierr) -! write(iout,*)'ener: tag=',tag - call MPI_RECV(x(1), nvar, MPI_DOUBLE_PRECISION, 0, tag,& - CG_COMM, status, ierr) -! write(iout,*)'x: tag=',tag - call MPI_RECV(x1(1), nvar, MPI_DOUBLE_PRECISION, 0, tag,& - CG_COMM, status, ierr) -! write(iout,*)'x1: tag=',tag -#ifdef AIX - call flush_(iout) -#else - call flush(iout) -#endif -! print *,'calling minimize' - call minimize(energyx,x,iretcode,nfun) - if(iprint.gt.0) & - write(iout,100)'minimized energy = ',energyx,& - ' # funeval:',nfun,' iret ',iretcode - write(*,100)'minimized energy = ',energyx,& - ' # funeval:',nfun,' iret ',iretcode - 100 format(a20,f10.5,a12,i5,a6,i2) - if(iretcode.eq.10) then - do iminrep=2,3 - if(iprint.gt.1) & - write(iout,*)' ... not converged - trying again ',iminrep - call minimize(energyx,x,iretcode,nfun) - if(iprint.gt.1) & - write(iout,*)'minimized energy = ',energyx,& - ' # funeval:',nfun,' iret ',iretcode - if(iretcode.ne.10)go to 812 - enddo - if(iretcode.eq.10) then - if(iprint.gt.1) & - write(iout,*)' ... not converged again - giving up' - go to 812 - endif - endif -812 continue -! print *,'Sending results' - call MPI_SEND(iitt, 1, MPI_INTEGER, 0, tag,& - CG_COMM, ierr) - call MPI_SEND(ener, 1, MPI_DOUBLE_PRECISION, 0, tag,& - CG_COMM, ierr) - call MPI_SEND(energyx, 1, MPI_DOUBLE_PRECISION, 0, tag,& - CG_COMM, ierr) - call MPI_SEND(x(1), nvar, MPI_DOUBLE_PRECISION, 0, tag,& - CG_COMM, ierr) - call MPI_SEND(x1(1), nvar, MPI_DOUBLE_PRECISION, 0, tag,& - CG_COMM, ierr) - call MPI_SEND(i_grnum_d, 1, MPI_INTEGER, 0, tag,& - CG_COMM, ierr) - call MPI_SEND(nfun, 1, MPI_INTEGER, 0, tag,& - CG_COMM, ierr) - call MPI_SEND(i_hesnum_d, 1, MPI_INTEGER, 0, tag,& - CG_COMM, ierr) -! print *,'End sending' - go to 1001 - endif - - return - end subroutine execute_slave -#endif -!----------------------------------------------------------------------------- - subroutine heat(over) - -! implicit real*8 (a-h,o-z) -! include 'DIMENSIONS' -! include 'COMMON.MCM' -! include 'COMMON.IOUNITS' - logical :: over -! Check if there`s a need to increase temperature. - if (ntrial.gt.maxtrial) then - if (NstepH.gt.0) then - if (dabs(Tcur-TMax).lt.1.0D-7) then - if (print_mc.gt.0) & - write (iout,'(/80(1h*)/a,f8.3,a/80(1h*))') & - 'Upper limit of temperature reached. Terminating.' - over=.true. - Tcur=Tmin - else - Tcur=Tcur*TstepH - if (Tcur.gt.Tmax) Tcur=Tmax - betbol=1.0D0/(Rbol*Tcur) - if (print_mc.gt.0) & - write (iout,'(/80(1h*)/a,f8.3,a,f10.5/80(1h*))') & - 'System heated up to ',Tcur,' K; BetBol:',betbol - ntherm=ntherm+1 - ntrial=0 - over=.false. - endif - else - if (print_mc.gt.0) & - write (iout,'(a)') & - 'Maximum number of trials in a single MCM iteration exceeded.' - over=.true. - Tcur=Tmin - endif - else - over=.false. - endif - return - end subroutine heat -!----------------------------------------------------------------------------- - subroutine cool - -! implicit real*8 (a-h,o-z) -! include 'DIMENSIONS' -! include 'COMMON.MCM' -! include 'COMMON.IOUNITS' - if (nstepC.gt.0 .and. dabs(Tcur-Tmin).gt.1.0D-7) then - Tcur=Tcur/TstepC - if (Tcur.lt.Tmin) Tcur=Tmin - betbol=1.0D0/(Rbol*Tcur) - if (print_mc.gt.0) & - write (iout,'(/80(1h*)/a,f8.3,a,f10.5/80(1h*))') & - 'System cooled down up to ',Tcur,' K; BetBol:',betbol - endif - return - end subroutine cool -!----------------------------------------------------------------------------- - subroutine perturb(error,lprint,MoveType,nbond,max_phi) - - use geometry - use energy, only:nnt,nct,itype - use md_calc, only:bond_move -! implicit real*8 (a-h,o-z) -! include 'DIMENSIONS' - integer,parameter :: MMaxSideMove=100 -! include 'COMMON.MCM' -! include 'COMMON.CHAIN' -! include 'COMMON.INTERACT' -! include 'COMMON.VAR' -! include 'COMMON.GEO' -! include 'COMMON.NAMES' -! include 'COMMON.IOUNITS' -!rc include 'COMMON.DEFORM1' - logical :: error,lprint,fail - integer :: MoveType,nbond,end_select,ind_side(MMaxSideMove) - real(kind=8) :: max_phi - real(kind=8) :: psi!,gen_psi -!el external iran_num -!el integer iran_num - integer :: ifour - -!!! local variables - el - integer :: itrial,iiwin,iwindow,isctry,i,icount,j,nstart,& - nside_move,inds,indx,ii,iti - real(kind=8) :: bond_prob,theta_new - - data ifour /4/ - error=.false. - lprint=.false. -! Perturb the conformation according to a randomly selected move. - call SelectMove(MoveType) -! write (iout,*) 'MoveType=',MoveType - itrial=0 - goto (100,200,300,400,500) MoveType -!------------------------------------------------------------------------------ -! Backbone N-bond move. -! Select the number of bonds (length of the segment to perturb). - 100 continue - if (itrial.gt.1000) then - write (iout,'(a)') 'Too many attempts at multiple-bond move.' - error=.true. - return - endif - bond_prob=ran_number(0.0D0,sumpro_bond(nbm)) -! print *,'sumpro_bond(nbm)=',sumpro_bond(nbm), -! & ' Bond_prob=',Bond_Prob - do i=1,nbm-1 -! print *,i,Bond_Prob,sumpro_bond(i),sumpro_bond(i+1) - if (bond_prob.ge.sumpro_bond(i) .and. & - bond_prob.le.sumpro_bond(i+1)) then - nbond=i+1 - goto 10 - endif - enddo - write (iout,'(2a)') 'In PERTURB: Error - number of bonds',& - ' to move out of range.' - error=.true. - return - 10 continue - if (nwindow.gt.0) then -! Select the first residue to perturb - iwindow=iran_num(1,nwindow) - print *,'iwindow=',iwindow - iiwin=1 - do while (winlen(iwindow).lt.nbond) - iwindow=iran_num(1,nwindow) - iiwin=iiwin+1 - if (iiwin.gt.1000) then - write (iout,'(a)') 'Cannot select moveable residues.' - error=.true. - return - endif - enddo - nstart=iran_num(winstart(iwindow),winend(iwindow)) - else - nstart = iran_num(koniecl+2,nres-nbond-koniecl) -!d print *,'nres=',nres,' nbond=',nbond,' koniecl=',koniecl, -!d & ' nstart=',nstart - endif - psi = gen_psi() - if (psi.eq.0.0) then - error=.true. - return - endif - if (print_mc.gt.1) write (iout,'(a,i4,a,i4,a,f8.3)') & - 'PERTURB: nbond=',nbond,' nstart=',nstart,' psi=',psi*rad2deg -!d print *,'nstart=',nstart - call bond_move(nbond,nstart,psi,.false.,error) - if (error) then - write (iout,'(2a)') & - 'Could not define reference system in bond_move, ',& - 'choosing ahother segment.' - itrial=itrial+1 - goto 100 - endif - nbond_move(nbond)=nbond_move(nbond)+1 - moves(1)=moves(1)+1 - nmove=nmove+1 - return -!------------------------------------------------------------------------------ -! Backbone endmove. Perturb a SINGLE angle of a residue close to the end of -! the chain. - 200 continue - lprint=.true. -! end_select=iran_num(1,2*koniecl) -! if (end_select.gt.koniecl) then -! end_select=nphi-(end_select-koniecl) -! else -! end_select=koniecl+3 -! endif -! if (nwindow.gt.0) then -! iwin=iran_num(1,nwindow) -! i1=max0(4,winstart(iwin)) -! i2=min0(winend(imin)+2,nres) -! end_select=iran_num(i1,i2) -! else -! iselect = iran_num(1,nmov_var) -! jj = 0 -! do i=1,nphi -! if (isearch_tab(i).eq.1) jj = jj+1 -! if (jj.eq.iselect) then -! end_select=i+3 -! exit -! endif -! enddo -! endif - end_select = iran_num(4,nres) - psi=max_phi*gen_psi() - if (psi.eq.0.0D0) then - error=.true. - return - endif - phi(end_select)=pinorm(phi(end_select)+psi) - if (print_mc.gt.1) write (iout,'(a,i4,a,f8.3,a,f8.3)') & - 'End angle',end_select,' moved by ',psi*rad2deg,' new angle:',& - phi(end_select)*rad2deg -! if (end_select.gt.3) -! & theta(end_select-1)=gen_theta(itype(end_select-2), -! & phi(end_select-1),phi(end_select)) -! if (end_select.lt.nres) -! & theta(end_select)=gen_theta(itype(end_select-1), -! & phi(end_select),phi(end_select+1)) -!d print *,'nres=',nres,' end_select=',end_select -!d print *,'theta',end_select-1,theta(end_select-1) -!d print *,'theta',end_select,theta(end_select) - moves(2)=moves(2)+1 - nmove=nmove+1 - lprint=.false. - return -!------------------------------------------------------------------------------ -! Side chain move. -! Select the number of SCs to perturb. - 300 isctry=0 - 301 nside_move=iran_num(1,MaxSideMove) -! print *,'nside_move=',nside_move,' MaxSideMove',MaxSideMove -! Select the indices. - do i=1,nside_move - icount=0 - 111 inds=iran_num(nnt,nct) - icount=icount+1 - if (icount.gt.1000) then - write (iout,'(a)')'Error - cannot select side chains to move.' - error=.true. - return - endif - if (itype(inds).eq.10) goto 111 - do j=1,i-1 - if (inds.eq.ind_side(j)) goto 111 - enddo - do j=1,i-1 - if (inds.lt.ind_side(j)) then - indx=j - goto 112 - endif - enddo - indx=i - 112 do j=i,indx+1,-1 - ind_side(j)=ind_side(j-1) - enddo - 113 ind_side(indx)=inds - enddo -! Carry out perturbation. - do i=1,nside_move - ii=ind_side(i) - iti=itype(ii) - call gen_side(iti,theta(ii+1),alph(ii),omeg(ii),fail) - if (fail) then - isctry=isctry+1 - if (isctry.gt.1000) then - write (iout,'(a)') 'Too many errors in SC generation.' - error=.true. - return - endif - goto 301 - endif - if (print_mc.gt.1) write (iout,'(2a,i4,a,2f8.3)') & - 'Side chain ',restyp(iti),ii,' moved to ',& - alph(ii)*rad2deg,omeg(ii)*rad2deg - enddo - moves(3)=moves(3)+1 - nmove=nmove+1 - return -!------------------------------------------------------------------------------ -! THETA move - 400 end_select=iran_num(3,nres) - theta_new=gen_theta(itype(end_select),phi(end_select),& - phi(end_select+1)) - if (print_mc.gt.1) write (iout,'(a,i3,a,f8.3,a,f8.3)') & - 'Theta ',end_select,' moved from',theta(end_select)*rad2deg,& - ' to ',theta_new*rad2deg - theta(end_select)=theta_new - moves(4)=moves(4)+1 - nmove=nmove+1 - return -!------------------------------------------------------------------------------ -! Error returned from SelectMove. - 500 error=.true. - return - end subroutine perturb -!----------------------------------------------------------------------------- - subroutine SelectMove(MoveType) - -! implicit real*8 (a-h,o-z) -! include 'DIMENSIONS' -! include 'COMMON.MCM' -! include 'COMMON.IOUNITS' - -!!! local variables - el - integer :: i,MoveType - real(kind=8) :: what_move - - what_move=ran_number(0.0D0,sumpro_type(MaxMoveType)) - do i=1,MaxMoveType - if (what_move.ge.sumpro_type(i-1).and. & - what_move.lt.sumpro_type(i)) then - MoveType=i - return - endif - enddo - write (iout,'(a)') & - 'Fatal error in SelectMoveType: cannot select move.' - MoveType=MaxMoveType+1 - return - end subroutine SelectMove -!----------------------------------------------------------------------------- - real(kind=8) function gen_psi() - - use geometry_data, only: angmin,pi -!el implicit none - integer :: i - real(kind=8) :: x !,ran_number -! include 'COMMON.IOUNITS' -! include 'COMMON.GEO' - x=0.0D0 - do i=1,100 - x=ran_number(-pi,pi) - if (dabs(x).gt.angmin) then - gen_psi=x - return - endif - enddo - write (iout,'(a)')'From Gen_Psi: Cannot generate angle increment.' - gen_psi=0.0D0 - return - end function gen_psi -!----------------------------------------------------------------------------- - subroutine metropolis(n,xcur,xold,ecur,eold,accepted,similar,enelower) - -! implicit real*8 (a-h,o-z) -! include 'DIMENSIONS' -! include 'COMMON.MCM' -! include 'COMMON.IOUNITS' -! include 'COMMON.VAR' -! include 'COMMON.GEO' -!rc include 'COMMON.DEFORM' - integer :: n - real(kind=8) :: ecur,eold,xx,bol !,ran_number - real(kind=8),dimension(n) :: xcur,xold - real(kind=8) :: ecut1 ,ecut2 ,tola - logical :: accepted,similar,not_done,enelower - logical :: lprn - -!!! local variables - el - real(kind=8) :: xxh,difene,reldife - - data ecut1 /-1.0D-5/,ecut2 /5.0D-3/,tola/5.0D0/ -! ecut1=-5*enedif -! ecut2=50*enedif -! tola=5.0d0 -! Set lprn=.true. for debugging. - lprn=.false. - if (lprn) & -!el write(iout,*)'enedif',enedif,' ecut1',ecut1,' ecut2',ecut2 - write(iout,*)' ecut1',ecut1,' ecut2',ecut2 - similar=.false. - enelower=.false. - accepted=.false. -! Check if the conformation is similar. - difene=ecur-eold - reldife=difene/dmax1(dabs(eold),dabs(ecur),1.0D0) - if (lprn) then - write (iout,*) 'Metropolis' - write(iout,*)'ecur,eold,difene,reldife',ecur,eold,difene,reldife - endif -! If energy went down remarkably, we accept the new conformation -! unconditionally. -!jp if (reldife.lt.ecut1) then - if (difene.lt.ecut1) then - accepted=.true. - EneLower=.true. - if (lprn) write (iout,'(a)') & - 'Conformation accepted, because energy has lowered remarkably.' -! elseif (reldife.lt.ecut2 .and. dif_ang(nphi,xcur,xold).lt.tola) -!jp elseif (reldife.lt.ecut2) - elseif (difene.lt.ecut2) & - then -! Reject the conf. if energy has changed insignificantly and there is not -! much change in conformation. - if (lprn) & - write (iout,'(2a)') 'Conformation rejected, because it is',& - ' similar to the preceding one.' - accepted=.false. - similar=.true. - else -! Else carry out Metropolis test. - EneLower=.false. - xx=ran_number(0.0D0,1.0D0) - xxh=betbol*difene - if (lprn) & - write (iout,*) 'betbol=',betbol,' difene=',difene,' xxh=',xxh - if (xxh.gt.50.0D0) then - bol=0.0D0 - else - bol=exp(-xxh) - endif - if (lprn) write (iout,*) 'bol=',bol,' xx=',xx - if (bol.gt.xx) then - accepted=.true. - if (lprn) write (iout,'(a)') & - 'Conformation accepted, because it passed Metropolis test.' - else - accepted=.false. - if (lprn) write (iout,'(a)') & - 'Conformation rejected, because it did not pass Metropolis test.' - endif - endif -#ifdef AIX - call flush_(iout) -#else - call flush(iout) -#endif - return - end subroutine metropolis -!----------------------------------------------------------------------------- - integer function conf_comp(x,ene) - - use geometry_data, only: nphi -! implicit real*8 (a-h,o-z) -! include 'DIMENSIONS' -! include 'COMMON.MCM' -! include 'COMMON.VAR' -! include 'COMMON.IOUNITS' -! include 'COMMON.GEO' - real(kind=8) :: etol, angtol - real(kind=8),dimension(6*nres) :: x !(maxvar) (maxvar=6*maxres) - real(kind=8) :: difa !dif_ang, - -!!! local variables - el - integer :: ii,i - real(kind=8) :: ene - - data etol /0.1D0/, angtol /20.0D0/ - do ii=nsave,1,-1 -! write (iout,*) 'ii=',ii,'ene=',ene,esave(ii),dabs(ene-esave(ii)) - if (dabs(ene-esave(ii)).lt.etol) then - difa=dif_ang(nphi,x,varsave(1,ii)) -! do i=1,nphi -! write(iout,'(i3,3f8.3)')i,rad2deg*x(i), -! & rad2deg*varsave(i,ii) -! enddo -! write(iout,*) 'ii=',ii,' difa=',difa,' angtol=',angtol - if (difa.le.angtol) then - if (print_mc.gt.0) then - write (iout,'(a,i5,2(a,1pe15.4))') & - 'Current conformation matches #',ii,& - ' in the store array ene=',ene,' esave=',esave(ii) -! write (*,'(a,i5,a)') 'Current conformation matches #',ii, -! & ' in the store array.' - endif ! print_mc.gt.0 - if (print_mc.gt.1) then - do i=1,nphi - write(iout,'(i3,3f8.3)')i,rad2deg*x(i),& - rad2deg*varsave(i,ii) - enddo - endif ! print_mc.gt.1 - nrepm=nrepm+1 - conf_comp=ii - return - endif - endif - enddo - conf_comp=0 - return - end function conf_comp -!----------------------------------------------------------------------------- - real(kind=8) function dif_ang(n,x,y) - - use geometry_data, only: dwapi -!el implicit none - integer :: i,n - real(kind=8),dimension(n) :: x,y - real(kind=8) :: w,wa,dif,difa -!el real(kind=8) :: pinorm -! include 'COMMON.GEO' - wa=0.0D0 - difa=0.0D0 - do i=1,n - dif=dabs(pinorm(y(i)-x(i))) - if (dabs(dif-dwapi).lt.dif) dif=dabs(dif-dwapi) - w=1.0D0-(2.0D0*(i-1)/(n-1)-1.0D0)**2+1.0D0/n - wa=wa+w - difa=difa+dif*dif*w - enddo - dif_ang=rad2deg*dsqrt(difa/wa) - return - end function dif_ang -!----------------------------------------------------------------------------- - subroutine add2cache(n1,n2,ncache,nvar,SourceID,CachSrc,ecur,xcur,ecache,xcache) - -! implicit none -! include 'COMMON.GEO' -! include 'COMMON.IOUNITS' - integer :: n1,n2,ncache,nvar,SourceID,CachSrc(n2) - integer :: i,ii,j - real(kind=8) :: ecur,xcur(nvar),ecache(n2),xcache(n1,n2) -!d write (iout,*) 'Enter ADD2CACHE ncache=',ncache ,' ecur',ecur -!d write (iout,'(10f8.3)') (rad2deg*xcur(i),i=1,nvar) -!d write (iout,*) 'Old CACHE array:' -!d do i=1,ncache -!d write (iout,*) 'i=',i,' ecache=',ecache(i),' CachSrc',CachSrc(i) -!d write (iout,'(10f8.3)') (rad2deg*xcache(j,i),j=1,nvar) -!d enddo - - i=ncache - do while (i.gt.0 .and. ecur.lt.ecache(i)) - i=i-1 - enddo - i=i+1 -!d write (iout,*) 'i=',i,' ncache=',ncache - if (ncache.eq.n2) then - write (iout,*) 'Cache dimension exceeded',ncache,n2 - write (iout,*) 'Highest-energy conformation will be removed.' - ncache=ncache-1 - endif - do ii=ncache,i,-1 - ecache(ii+1)=ecache(ii) - CachSrc(ii+1)=CachSrc(ii) - do j=1,nvar - xcache(j,ii+1)=xcache(j,ii) - enddo - enddo - ecache(i)=ecur - CachSrc(i)=SourceID - do j=1,nvar - xcache(j,i)=xcur(j) - enddo - ncache=ncache+1 -!d write (iout,*) 'New CACHE array:' -!d do i=1,ncache -!d write (iout,*) 'i=',i,' ecache=',ecache(i),' CachSrc',CachSrc(i) -!d write (iout,'(10f8.3)') (rad2deg*xcache(j,i),j=1,nvar) -!d enddo - return - end subroutine add2cache -!----------------------------------------------------------------------------- - subroutine rm_from_cache(i,n1,n2,ncache,nvar,CachSrc,ecache,xcache) - -! implicit none -! include 'COMMON.GEO' -! include 'COMMON.IOUNITS' - integer :: n1,n2,ncache,nvar,CachSrc(n2) - integer :: i,ii,j - real(kind=8) :: ecache(n2),xcache(n1,n2) - -!d write (iout,*) 'Enter RM_FROM_CACHE' -!d write (iout,*) 'Old CACHE array:' -!d do ii=1,ncache -!d write (iout,*)'i=',ii,' ecache=',ecache(ii),' CachSrc',CachSrc(ii) -!d write (iout,'(10f8.3)') (rad2deg*xcache(j,ii),j=1,nvar) -!d enddo - - do ii=i+1,ncache - ecache(ii-1)=ecache(ii) - CachSrc(ii-1)=CachSrc(ii) - do j=1,nvar - xcache(j,ii-1)=xcache(j,ii) - enddo - enddo - ncache=ncache-1 -!d write (iout,*) 'New CACHE array:' -!d do i=1,ncache -!d write (iout,*) 'i=',i,' ecache=',ecache(i),' CachSrc',CachSrc(i) -!d write (iout,'(10f8.3)') (rad2deg*xcache(j,i),j=1,nvar) -!d enddo - return - end subroutine rm_from_cache -!----------------------------------------------------------------------------- -! mcm.F io_mcm -!----------------------------------------------------------------------------- - subroutine statprint(it,nfun,iretcode,etot,elowest) - - use control_data, only: MaxMoveType,minim - use control, only: tcpu - use mcm_data -! implicit real*8 (a-h,o-z) -! include 'DIMENSIONS' -! include 'COMMON.IOUNITS' -! include 'COMMON.CONTROL' -! include 'COMMON.MCM' -!el local variables - integer :: it,nfun,iretcode,i - real(kind=8) :: etot,elowest,fr_mov_i - - if (minim) then - write (iout,& - '(80(1h*)/a,i5,a,1pe14.5,a,1pe14.5/a,i3,a,i10,a,i5,a,i5)') & - 'Finished iteration #',it,' energy is',etot,& - ' lowest energy:',elowest,& - 'SUMSL return code:',iretcode,& - ' # of energy evaluations:',neneval,& - '# of temperature jumps:',ntherm,& - ' # of minima repetitions:',nrepm - else - write (iout,'(80(1h*)/a,i8,a,1pe14.5,a,1pe14.5)') & - 'Finished iteration #',it,' energy is',etot,& - ' lowest energy:',elowest - endif - write (iout,'(/4a)') & - 'Kind of move ',' total',' accepted',& - ' fraction' - write (iout,'(58(1h-))') - do i=-1,MaxMoveType - if (moves(i).eq.0) then - fr_mov_i=0.0d0 - else - fr_mov_i=dfloat(moves_acc(i))/dfloat(moves(i)) - endif - write(iout,'(a,2i15,f10.5)')MovTypID(i),moves(i),moves_acc(i),& - fr_mov_i - enddo - write (iout,'(a,2i15,f10.5)') 'total ',nmove,nacc_tot,& - dfloat(nacc_tot)/dfloat(nmove) - write (iout,'(58(1h-))') - write (iout,'(a,1pe12.4)') 'Elapsed time:',tcpu() - return - end subroutine statprint -!----------------------------------------------------------------------------- - subroutine zapis(varia,etot) - - use geometry_data, only: nres,rad2deg,nvar - use mcm_data - use MPI_data -! implicit real*8 (a-h,o-z) -! include 'DIMENSIONS' -#ifdef MPI - use MPI_data !include 'COMMON.INFO' - include 'mpif.h' -#endif -! include 'COMMON.GEO' -! include 'COMMON.VAR' -! include 'COMMON.MCM' -! include 'COMMON.IOUNITS' - integer,dimension(nsave) :: itemp - real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres) - logical :: lprint -!el local variables - integer :: j,i,maxvar - real(kind=8) :: etot - -!el allocate(esave(nsave)) !(maxsave) - - maxvar=6*nres - lprint=.false. - if (lprint) then - write (iout,'(a,i5,a,i5)') 'Enter ZAPIS NSave=',Nsave,& - ' MaxSave=',MaxSave - write (iout,'(a)') 'Current energy and conformation:' - write (iout,'(1pe14.5)') etot - write (iout,'(10f8.3)') (rad2deg*varia(i),i=1,nvar) - endif -! Shift the contents of the esave and varsave arrays if filled up. - call add2cache(6*nres,nsave,nsave,nvar,MyID,itemp,& - etot,varia,esave,varsave) - if (lprint) then - write (iout,'(a)') 'Energies and the VarSave array.' - do i=1,nsave - write (iout,'(i5,1pe14.5)') i,esave(i) - write (iout,'(10f8.3)') (rad2deg*varsave(j,i),j=1,nvar) - enddo - endif - return - end subroutine zapis -!----------------------------------------------------------------------------- -!----------------------------------------------------------------------------- - subroutine alloc_MCM_arrays - - use energy_data, only: max_ene - use MPI_data -! common.mce -! common /mce/ - allocate(entropy(-max_ene-4:max_ene)) !(-max_ene-4:max_ene) - allocate(nhist(-max_ene:max_ene)) !(-max_ene:max_ene) - allocate(nminima(maxsave)) !(maxsave) -! common /pool/ - allocate(xpool(6*nres,max_pool)) !(maxvar,max_pool)(maxvar=6*maxres) - allocate(epool(max_pool)) !(max_pool) -! commom.mcm -! common /mcm/ - if(.not.allocated(nsave_part)) allocate(nsave_part(nctasks)) !(max_cg_procs) -! common /move/ -! in io: mcmread -! real(kind=8),dimension(:),allocatable :: sumpro_type !(0:MaxMoveType) - allocate(sumpro_bond(0:nres)) !(0:maxres) - allocate(nbond_move(nres),nbond_acc(nres)) !(maxres) - allocate(moves(-1:MaxMoveType+1),moves_acc(-1:MaxMoveType+1)) !(-1:MaxMoveType+1) -! common /accept_stats/ -! allocate(nacc_part !(0:MaxProcs) !el nie uzywane??? -! common /windows/ in io: mcmread -! allocate(winstart,winend,winlen !(maxres) -! common /moveID/ -!el allocate(MovTypID(-1:MaxMoveType+1)) !(-1:MaxMoveType+1) -! common.var -! common /oldgeo/ - allocate(varsave(nres*6,maxsave)) !(maxvar,maxsave)(maxvar=6*maxres) - allocate(esave(maxsave)) !(maxsave) - allocate(Origin(maxsave)) !(maxsave) - - return - end subroutine alloc_MCM_arrays -!----------------------------------------------------------------------------- -!----------------------------------------------------------------------------- - end module mcm_md