+++ /dev/null
- 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
projects / unres4.git / blobdiff