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