2 !-----------------------------------------------------------------------------
6 use geometry_data, only: nres,nvar,rad2deg
7 use random, only: iran_num,ran_number
14 !-----------------------------------------------------------------------------
15 ! Max. number of move types in MCM
16 ! integer,parameter :: maxmovetype=4
17 !-----------------------------------------------------------------------------
18 ! Max. number of conformations in Master's cache array
19 integer,parameter :: max_cache=10
20 !-----------------------------------------------------------------------------
21 ! Max. number of stored confs. in MC/MCM simulation
22 ! integer,parameter :: maxsave=20
23 !-----------------------------------------------------------------------------
24 ! Number of threads in deformation
25 integer,parameter :: max_thread=4, max_thread2=2*max_thread
26 !-----------------------------------------------------------------------------
27 ! Number of structures to compare at t=0
28 integer,parameter :: max_threadss=8,max_threadss2=2*max_threadss
29 !-----------------------------------------------------------------------------
30 ! Max. number of conformations in the pool
31 integer,parameter :: max_pool=10
32 !-----------------------------------------------------------------------------
33 !-----------------------------------------------------------------------------
37 ! integer,dimension(max_cache) :: CachSrc nie używane
38 ! integer,dimension(max_cache) :: isent,iused
39 ! logical :: cache_update
40 ! real(kind=8),dimension(max_cache) :: ecache
41 ! real(kind=8),dimension(:,:),allocatable :: xcache !(maxvar,max_cache)
42 !-----------------------------------------------------------------------------
45 ! real(kind=8) :: emin,emax
46 real(kind=8),dimension(:),allocatable :: entropy !(-max_ene-4:max_ene)
47 real(kind=8),dimension(:),allocatable :: nhist !(-max_ene:max_ene)
48 real(kind=8),dimension(:),allocatable :: nminima !(maxsave)
51 integer :: indminn,indmaxx
54 ! real(kind=8) :: pool_fraction
55 real(kind=8),dimension(:,:),allocatable :: xpool !(maxvar,max_pool)
56 real(kind=8),dimension(:),allocatable :: epool !(max_pool)
57 ! common /mce_counters/
58 !------------------------------------------------------------------------------
59 !... Following COMMON block contains variables controlling motion.
60 !------------------------------------------------------------------------------
62 ! real(kind=8),dimension(0:MaxMoveType) :: sumpro_type !(0:MaxMoveType)
63 real(kind=8),dimension(:),allocatable :: sumpro_bond !(0:maxres)
64 integer :: koniecl,Nbm,MaxSideMove!,nmove
65 integer,dimension(:),allocatable :: nbond_move,nbond_acc !(maxres)
66 ! integer,dimension(-1:MaxMoveType+1) :: moves,moves_acc !(-1:MaxMoveType+1)
67 ! common /accept_stats/
69 integer,dimension(:),allocatable :: nacc_part !(0:MaxProcs) !el nie uzywane???
72 ! integer,dimension(:),allocatable :: winstart,winend,winlen !(maxres)
74 ! character(len=16),dimension(-1:MaxMoveType+1) :: MovTypID !(-1:MaxMoveType+1)
75 !------------------------------------------------------------------------------
76 !... koniecl - the number of bonds to be considered "end bonds" subjected to
78 !... Nbm - The maximum length of N-bond segment to be moved;
79 !... MaxSideMove - maximum number of side chains subjected to local moves
81 !... nmove - the current number of attempted moves;
82 !... nbond_move(*) array that stores the total numbers of 2-bond,3-bond,...
84 !... nendmove - number of endmoves;
85 !... nbackmove - number of backbone moves;
86 !... nsidemove - number of local side chain moves;
87 !... sumpro_type(*) - array that stores the lower and upper boundary of the
88 !... random-number range that determines the type of move
89 !... (N-bond, backbone or side chain);
90 !... sumpro_bond(*) - array that stores the probabilities to perform bond
91 !... moves of consecutive segment length.
92 !... winstart(*) - the starting position of the perturbation window;
93 !... winend(*) - the end position of the perturbation window;
94 !... winlen(*) - length of the perturbation window;
95 !... nwindow - the number of perturbation windows (0 - entire chain).
96 !-----------------------------------------------------------------------------
99 !-----------------------------------------------------------------------------
101 !-----------------------------------------------------------------------------
103 !-----------------------------------------------------------------------------
104 subroutine compare_s1(n_thr,num_thread_save,energyx,x,icomp,enetbss,&
105 coordss,rms_d,modif,iprint)
107 ! This subroutine compares the new conformation, whose variables are in X
108 ! with the previously accumulated conformations whose energies and variables
109 ! are stored in ENETBSS and COORDSS, respectively. The meaning of other
110 ! variables is as follows:
112 ! N_THR - on input the previous # of accumulated confs, on output the current
113 ! # of accumulated confs.
114 ! N_REPEAT - an array that indicates how many times the structure has already
115 ! been used to start the reversed-reversing procedure. Addition of
116 ! a new structure replacement of a structure with a similar, but
117 ! lower-energy structure resets the respective entry in N_REPEAT to zero
119 ! ENERGYX,X - the energy and variables of the new conformations.
120 ! ICOMP - comparison result:
121 ! 0 - the new structure is similar to one of the previous ones and does
122 ! not have a remarkably lower energy and is therefore rejected;
123 ! 1 - the new structure is different and is added to the set, because
124 ! there is still room in the COORDSS and ENETBSS arrays;
125 ! 2 - the new structure is different, but higher in energy than any
126 ! previous one and is therefore rejected
127 ! 3 - there is no more room in the COORDSS and ENETBSS arrays, but
128 ! the new structure is lower in energy than at least the highest-
129 ! energy previous structure and therefore replaces it.
130 ! 9 - the new structure is similar to a number of previous structures,
131 ! but has a remarkably lower energy than any of them; therefore
132 ! replaces all these structures;
133 ! MODIF - a logical variable that shows whether to include the new structure
134 ! in the set of accumulated structures
136 ! implicit real*8 (a-h,o-z)
138 use regularize_, only:fitsq
139 ! include 'DIMENSIONS'
140 ! include 'COMMON.GEO'
141 ! include 'COMMON.VAR'
142 !rc include 'COMMON.DEFORM'
143 ! include 'COMMON.IOUNITS'
145 !el use geometry_data !include 'COMMON.CHAIN'
148 real(kind=8),dimension(6*nres) :: x,x1 !(maxvar) (maxvar=6*maxres)
149 real(kind=8) :: przes(3),obrot(3,3)
150 integer :: list(max_thread)
151 logical :: non_conv,modif
152 real(kind=8) :: enetbss(max_threadss)
153 real(kind=8) :: coordss(6*nres,max_threadss)
155 !!! local variables - el
156 integer :: n_thr,num_thread_save,icomp,minimize_s_flag,iprint
157 real(kind=8) :: energyx,energyy,rms_d
158 integer :: nlist,k,kk,j,i,iresult
159 real(kind=8) :: enex_jp,roznica
163 call var_to_geom(nvar,x)
164 write(iout,*) 'Warning calling chainbuild'
172 ! write(iout,*)'*ene=',energyx
179 if (iprint.gt.3) then
180 write (iout,*) 'Compare_ss, i=',i
181 write (iout,*) 'New structure Energy:',energyx
182 write (iout,'(10f8.3)') (rad2deg*x(k),k=1,nvar)
183 write (iout,*) 'Template structure Energy:',enetbss(i)
184 write (iout,'(10f8.3)') (rad2deg*x1(k),k=1,nvar)
188 call var_to_geom(nvar,x1)
189 write(iout,*) 'Warning calling chainbuild'
191 !d write(iout,*)'C and CREF'
192 !d write(iout,'(i5,3f10.5,5x,3f10.5)')(k,(c(j,k),j=1,3),
193 !d & (cref(j,k),j=1,3),k=1,nres)
194 call fitsq(roznica,c(1,1),cref(1,1,1),nres,przes,obrot,non_conv)
196 print *,'Problems in FITSQ!!!'
198 print '(10f8.3)',(x(k),k=1,nvar)
200 print '(10f8.3)',(x1(k),k=1,nvar)
202 print '(i5,3f10.5,5x,3f10.5)',(k,(c(j,k),j=1,3),&
203 (cref(j,k,1),j=1,3),k=1,nres)
205 roznica=dsqrt(dabs(roznica))
207 if (roznica.lt.rms_d) iresult = 0
210 !el call cmprs(x,x1,roznica,energyx,energyy,iresult)
212 if (iprint.gt.1) write(iout,'(i5,f10.6,$)') i,roznica
213 ! print '(i5,f8.3)',i,roznica
214 if(iresult.eq.0) then
217 if (iprint.gt.1) write(iout,*)
218 if(energyx.ge.enetbss(i)) then
220 write(iout,*)'s*>> structure rejected - same as nr ',i, &
227 if(energyx.lt.enetbss(i).and.enex_jp.lt.enetbss(i))then
232 if (iprint.gt.1) write(iout,*)
236 write(iout,'(a,i3,$)')'s*>> structure accepted1 - repl nr ',&
240 write(iout,'(a,i3)') &
241 's*>> structure accepted1 - would repl nr ',list(1)
244 if (.not. modif) goto 1106
251 if (iprint.gt.1) write(iout,'(i3,$)')list(j)
252 do kk=list(j)+1,nlist
253 enetbss(kk-1)=enetbss(kk)
255 coordss(i,kk-1)=coordss(i,kk)
259 if (iprint.gt.1) write(iout,*)
262 if(n_thr.lt.num_thread_save) then
266 write(iout,*)'s*>> structure accepted - add with nr ',n_thr+1
269 write(iout,*)'s*>> structure accepted - would add with nr ',&
274 enetbss(n_thr)=energyx
276 coordss(i,n_thr)=x(i)
281 write(iout,*)'s*>> structure rejected - too high energy'
288 write(iout,*)'s*>> structure accepted - repl nr ',j
291 write(iout,*)'s*>> structure accepted - would repl nr ',j
302 end subroutine compare_s1
303 !-----------------------------------------------------------------------------
305 !-----------------------------------------------------------------------------
310 use MPI_data, only:WhatsUp,MyID
311 use compare_data, only: ener
312 use control_data, only: minim,refstr
314 use regularize_, only:fitsq
315 use control, only: tcpu,ovrtim
317 use minimm, only:minimize
318 ! Does modified entropic sampling in the space of minima.
319 ! implicit real*8 (a-h,o-z)
320 ! include 'DIMENSIONS'
321 ! include 'COMMON.IOUNITS'
323 use MPI_data !include 'COMMON.INFO'
325 ! include 'COMMON.GEO'
326 ! include 'COMMON.CHAIN'
327 ! include 'COMMON.MCM'
328 ! include 'COMMON.MCE'
329 ! include 'COMMON.CONTACTS'
330 ! include 'COMMON.CONTROL'
331 ! include 'COMMON.VAR'
332 ! include 'COMMON.INTERACT'
333 ! include 'COMMON.THREAD'
334 ! include 'COMMON.NAMES'
335 logical :: accepted,not_done,over,error,lprint !,ovrtim
336 integer :: MoveType,nbond
337 ! integer :: conf_comp
338 real(kind=8) :: RandOrPert
339 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
340 real(kind=8) :: elowest,ehighest,eold
341 real(kind=8) :: przes(3),obr(3,3)
342 real(kind=8),dimension(6*nres) :: varold !(maxvar) (maxvar=6*maxres)
344 real(kind=8),dimension(0:n_ene) :: energia,energia_ave
346 !!! local variables -el
347 integer :: i,ii,kkk,it,j,nacc,nfun,ijunk,indmin,indmax,&
348 ISWEEP,Kwita,iretcode,indeold,iene,noverlap,&
349 irep,nstart_grow,inde
350 real(kind=8) :: facee,conste,ejunk,etot,rms,co,frac,&
351 deix,dent,sold,scur,runtime
354 ! if(.not.allocated(varsave)) allocate(varsave(maxvar,maxsave)) !(maxvar,maxsave)
355 !d write (iout,*) 'print_mc=',print_mc
358 !---------------------------------------------------------------------------
359 ! Initialize counters.
360 !---------------------------------------------------------------------------
361 ! Total number of generated confs.
363 ! Total number of moves. In general this won't be equal to the number of
364 ! attempted moves, because we may want to reject some "bad" confs just by
367 ! Total number of shift (nbond_move(1)), spike, crankshaft, three-bond,...
369 !el allocate(nbond_move(nres)) !(maxres)
374 ! Initialize total and accepted number of moves of various kind.
379 ! Total number of energy evaluations.
385 facee=1.0D0/(maxacc*delte)
387 ! Read entropy from previous simulations.
389 read (ientin,*) indminn,indmaxx,emin,emax
390 print *,'indminn=',indminn,' indmaxx=',indmaxx,' emin=',emin,&
392 do i=-max_ene,max_ene
393 entropy(i)=(emin+i*delte)*betbol
395 read (ientin,*) (ijunk,ejunk,entropy(i),i=indminn,indmaxx)
398 write (iout,*) 'indminn=',indminn,' indmaxx=',indmaxx,&
399 ' emin=',emin,' emax=',emax
400 write (iout,'(/a)') 'Initial entropy'
402 write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
405 ! Read the pool of conformations
407 !----------------------------------------------------------------------------
408 ! Entropy-sampling simulations with continually updated entropy
409 ! Loop thru simulations
410 !----------------------------------------------------------------------------
412 !----------------------------------------------------------------------------
413 ! Take a conformation from the pool
414 !----------------------------------------------------------------------------
420 write (iout,*) 'Took conformation',ii,' from the pool energy=',&
422 call var_to_geom(nvar,varia)
423 ! Print internal coordinates of the initial conformation
426 call gen_rand_conf(1,*20)
428 !----------------------------------------------------------------------------
429 ! Compute and print initial energies.
430 !----------------------------------------------------------------------------
433 allocate(nsave_part(nctasks))
434 if (MyID.eq.MasterID) then
442 write (iout,'(/80(1h*)/a,i2/80(1h*)/)') 'MCE iteration #',isweep
443 write (iout,'(/80(1h*)/a)') 'Initial energies:'
444 write(iout,*) 'Warning calling chainbuild'
448 call enerprint(energia)
449 ! Minimize the energy of the first conformation.
451 call geom_to_var(nvar,varia)
452 call minimize(etot,varia,iretcode,nfun)
455 write (iout,'(/80(1h*)/a/80(1h*))') &
456 'Results of the first energy minimization:'
457 call enerprint(energia)
461 call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,kkk),&
465 call contact(.false.,ncont,icont,co)
466 frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
467 write (iout,'(a,f8.3,a,f8.3,a,f8.3)') &
468 'RMS deviation from the reference structure:',rms,&
469 ' % of native contacts:',frac*100,' contact order:',co
470 write (istat,'(i5,11(1pe14.5))') 0,&
471 (energia(print_order(i)),i=1,nprint_ene),etot,rms,frac,co
473 write (istat,'(i5,9(1pe14.5))') 0,&
474 (energia(print_order(i)),i=1,nprint_ene),etot
477 neneval=neneval+nfun+1
478 if (.not. ent_read) then
479 ! Initialize the entropy array
480 do i=-max_ene,max_ene
482 ! Uncomment the line below for actual entropic sampling (start with uniform
483 ! energy distribution).
485 ! Uncomment the line below for multicanonical sampling (start with Boltzmann
487 entropy(i)=(emin+i*delte)*betbol
491 write (iout,'(/a)') 'Initial entropy'
493 write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
497 call recv_stop_sig(Kwita)
498 if (whatsup.eq.1) then
499 call send_stop_sig(-2)
501 else if (whatsup.le.-2) then
503 else if (whatsup.eq.2) then
509 not_done = (iretcode.ne.11)
511 write (iout,'(/80(1h*)/20x,a/80(1h*))') &
512 'Enter Monte Carlo procedure.'
514 call briefout(0,etot)
519 indeold=(eold-emin)/delte
520 deix=eold-(emin+indeold*delte)
521 dent=entropy(indeold+1)-entropy(indeold)
522 !d write (iout,*) 'indeold=',indeold,' deix=',deix,' dent=',dent
523 !d write (*,*) 'Processor',MyID,' indeold=',indeold,' deix=',deix,
525 sold=entropy(indeold)+(dent/delte)*deix
527 write (iout,*) 'eold=',eold,' sold=',sold,' elowest=',etot
528 write (*,*) 'Processor',MyID,' eold=',eold,' sold=',sold,&
530 if (minim) call zapis(varia,etot)
532 ! NACC is the counter for the accepted conformations of a given processor
534 ! NACC_TOT counts the total number of accepted conformations
537 if (MyID.eq.MasterID) then
538 call receive_MCM_info
540 call send_MCM_info(2)
543 do iene=indminn,indmaxx
550 !----------------------------------------------------------------------------
556 if (print_mc.gt.0) write (iout,'(80(1h*)/20x,a,i7)') &
557 'Beginning iteration #',it
558 ! Initialize local counter.
559 ntrial=0 ! # of generated non-overlapping confs.
560 noverlap=0 ! # of overlapping confs.
562 do while (.not. accepted .and. WhatsUp.eq.0 .and. Kwita.eq.0)
564 ! Retrieve the angles of previously accepted conformation
568 !d write (iout,'(a)') 'Old variables:'
569 !d write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
570 call var_to_geom(nvar,varia)
572 write(iout,*) 'Warning calling chainbuild'
577 ! Decide whether to generate a random conformation or perturb the old one
578 RandOrPert=ran_number(0.0D0,1.0D0)
579 if (RandOrPert.gt.RanFract) then
581 write (iout,'(a)') 'Perturbation-generated conformation.'
582 call perturb(error,lprint,MoveType,nbond,1.0D0)
584 if (MoveType.lt.1 .or. MoveType.gt.MaxMoveType) then
585 write (iout,'(/a,i7,a/)') 'Error - unknown MoveType=',&
586 MoveType,' returned from PERTURB.'
589 write(iout,*) 'Warning calling chainbuild'
594 nstart_grow=iran_num(3,nres)
596 write (iout,'(2a,i3)') 'Random-generated conformation',&
597 ' - chain regrown from residue',nstart_grow
598 call gen_rand_conf(nstart_grow,*30)
600 call geom_to_var(nvar,varia)
601 !d write (iout,'(a)') 'New variables:'
602 !d write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
604 if (print_mc.gt.0) write (iout,'(a,i5,a,i10,a,i10)') &
605 'Processor',MyId,' trial move',ntrial,' total generated:',ngen
606 if (print_mc.gt.0) write (*,'(a,i5,a,i10,a,i10)') &
607 'Processor',MyId,' trial move',ntrial,' total generated:',ngen
610 ! call enerprint(energia(0))
611 ! write (iout,'(2(a,1pe14.5))') 'Etot=',Etot,' Elowest=',Elowest
612 if (etot-elowest.gt.overlap_cut) then
613 write (iout,'(a,i5,a,1pe14.5)') 'Iteration',it,&
614 ' Overlap detected in the current conf.; energy is',etot
618 if (noverlap.gt.maxoverlap) then
619 write (iout,'(a)') 'Too many overlapping confs.'
624 call minimize(etot,varia,iretcode,nfun)
625 !d write (iout,'(a)') 'Variables after minimization:'
626 !d write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
629 neneval=neneval+nfun+1
631 if (print_mc.gt.2) then
632 write (iout,'(a)') 'Total energies of trial conf:'
633 call enerprint(energia)
634 else if (print_mc.eq.1) then
635 write (iout,'(a,i6,a,1pe16.6)') &
636 'Trial conformation:',ngen,' energy:',etot
638 !--------------------------------------------------------------------------
640 !--------------------------------------------------------------------------
643 call accepting(etot,eold,scur,sold,varia,varold,&
648 if (elowest.gt.etot) elowest=etot
649 if (ehighest.lt.etot) ehighest=etot
650 moves_acc(MoveType)=moves_acc(MoveType)+1
651 if (MoveType.eq.1) then
652 nbond_acc(nbond)=nbond_acc(nbond)+1
654 ! Check against conformation repetitions.
655 irep=conf_comp(varia,etot)
656 #if defined(AIX) || defined(PGI)
657 open (istat,file=statname,position='append')
659 open (istat,file=statname,access='append')
663 call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,kkk),&
667 call contact(.false.,ncont,icont,co)
668 frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
670 write (iout,'(a,f8.3,a,f8.3,a,f8.3)') &
671 'RMS deviation from the reference structure:',rms,&
672 ' % of native contacts:',frac*100,' contact order:',co
674 write (istat,'(i5,11(1pe14.5))') it,&
675 (energia(print_order(i)),i=1,nprint_ene),etot,&
677 elseif (print_stat) then
678 write (istat,'(i5,10(1pe14.5))') it,&
679 (energia(print_order(i)),i=1,nprint_ene),etot
683 call statprint(nacc,nfun,iretcode,etot,elowest)
684 ! Print internal coordinates.
685 if (print_int) call briefout(nacc,etot)
687 if (MyID.ne.MasterID) then
688 call recv_stop_sig(Kwita)
689 !d print *,'Processor:',MyID,' STOP=',Kwita
691 call send_MCM_info(2)
693 call send_MCM_info(1)
697 ! Store the accepted conf. and its energy.
705 !d write (iout,*) 'Accepted conformation:'
706 !d write (iout,*) (rad2deg*varia(i),i=1,nphi)
707 if (minim) call zapis(varia,etot)
709 ener(i,nsave)=energia(i)
711 ener(n_ene+1,nsave)=etot
712 ener(n_ene+2,nsave)=frac
714 nminima(irep)=nminima(irep)+1.0D0
715 ! print *,'irep=',irep,' nminima=',nminima(irep)
717 if (Kwita.eq.0) call recv_stop_sig(kwita)
722 if (MyID.eq.MasterID) then
723 call receive_MCM_info
724 if (nacc_tot.ge.maxacc) accepted=.true.
727 if (ntrial.gt.maxtrial_iter .and. npool.gt.0) then
728 ! Take a conformation from the pool
733 write (iout,*) 'Iteration',it,' max. # of trials exceeded.'
735 'Take conformation',ii,' from the pool energy=',epool(ii)
737 write (iout,'(10f8.3)') (rad2deg*varia(i),i=1,nvar)
739 endif ! (ntrial.gt.maxtrial_iter .and. npool.gt.0)
743 if (MyID.eq.MasterID) then
744 call receive_MCM_info
746 if (Kwita.eq.0) call recv_stop_sig(kwita)
748 if (ovrtim()) WhatsUp=-1
749 !d write (iout,*) 'WhatsUp=',WhatsUp,' Kwita=',Kwita
750 not_done = (nacc_tot.lt.maxacc) .and. (WhatsUp.eq.0) &
752 !d write (iout,*) 'not_done=',not_done
755 print *,'Processor',MyID,&
756 ' has received STOP signal =',Kwita,' in EntSamp.'
757 !d print *,'not_done=',not_done
758 if (Kwita.lt.-1) WhatsUp=Kwita
759 else if (nacc_tot.ge.maxacc) then
760 print *,'Processor',MyID,' calls send_stop_sig,',&
761 ' because a sufficient # of confs. have been collected.'
762 !d print *,'not_done=',not_done
763 call send_stop_sig(-1)
764 else if (WhatsUp.eq.-1) then
765 print *,'Processor',MyID,&
766 ' calls send_stop_sig because of timeout.'
767 !d print *,'not_done=',not_done
768 call send_stop_sig(-2)
773 !-----------------------------------------------------------------
774 !... Construct energy histogram & update entropy
775 !-----------------------------------------------------------------
779 write (iout,*) 'Processor',MyID,&
780 ' is broadcasting ERROR-STOP signal.'
781 write (*,*) 'Processor',MyID,&
782 ' is broadcasting ERROR-STOP signal.'
783 call send_stop_sig(-3)
787 if (MyID.eq.MasterID) then
788 ! call receive_MCM_results
789 call receive_energies
792 if (esave(i).lt.elowest) elowest=esave(i)
793 if (esave(i).gt.ehighest) ehighest=esave(i)
795 write (iout,'(a,i10)') '# of accepted confs:',nacc_tot
796 write (iout,'(a,f10.5,a,f10.5)') 'Lowest energy:',elowest,&
797 ' Highest energy',ehighest
798 if (isweep.eq.1 .and. .not.ent_read) then
801 write (iout,*) 'EMAX=',emax
803 indmaxx=(ehighest-emin)/delte
806 do i=-max_ene,max_ene
807 entropy(i)=(emin+i*delte)*betbol
811 indmin=(elowest-emin)/delte
812 indmax=(ehighest-emin)/delte
813 if (indmin.lt.indminn) indminn=indmin
814 if (indmax.gt.indmaxx) indmaxx=indmax
816 write(iout,*)'indminn=',indminn,' indmaxx=',indmaxx
817 ! Construct energy histogram
819 inde=(esave(i)-emin)/delte
820 nhist(inde)=nhist(inde)+nminima(i)
822 ! Update entropy (density of states)
824 if (nhist(i).gt.0) then
825 entropy(i)=entropy(i)+dlog(nhist(i)+0.0D0)
829 !d entropy(i)=1.0D+10
831 write (iout,'(/80(1h*)/a,i2/80(1h*)/)') &
832 'End of macroiteration',isweep
833 write (iout,'(a,f10.5,a,f10.5)') 'Elowest=',elowest,&
834 ' Ehighest=',ehighest
835 write (iout,'(a)') 'Frequecies of minima'
837 write (iout,'(i5,f5.0,f10.5)') i,nminima(i),esave(i)
839 write (iout,'(/a)') 'Energy histogram'
841 write (iout,'(i5,2f10.5)') i,emin+i*delte,nhist(i)
843 write (iout,'(/a)') 'Entropy'
845 write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
847 !-----------------------------------------------------------------
848 !... End of energy histogram construction
849 !-----------------------------------------------------------------
851 entropy(-max_ene-4)=dfloat(indminn)
852 entropy(-max_ene-3)=dfloat(indmaxx)
853 entropy(-max_ene-2)=emin
854 entropy(-max_ene-1)=emax
856 !d print *,entname,ientout
857 open (ientout,file=entname,status='unknown')
858 write (ientout,'(2i5,2e25.17)') indminn,indmaxx,emin,emax
860 write (ientout,'(i5,f10.5,f20.15)') i,emin+i*delte,entropy(i)
864 write (iout,'(a)') 'Frequecies of minima'
866 write (iout,'(i5,f5.0,f10.5)') i,nminima(i),esave(i)
868 ! call send_MCM_results
870 call receive_MCM_update
871 indminn=entropy(-max_ene-4)
872 indmaxx=entropy(-max_ene-3)
873 emin=entropy(-max_ene-2)
874 emax=entropy(-max_ene-1)
875 write (iout,*) 'Received from master:'
876 write (iout,*) 'indminn=',indminn,' indmaxx=',indmaxx,&
877 ' emin=',emin,' emax=',emax
878 write (iout,'(/a)') 'Entropy'
880 write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
883 if (WhatsUp.lt.-1) return
885 if (ovrtim() .or. WhatsUp.lt.0) return
888 write (iout,'(/80(1h*)/20x,a)') 'Summary run statistics:'
889 call statprint(nacc,nfun,iretcode,etot,elowest)
891 'Statistics of multiple-bond motions. Total motions:'
892 write (iout,'(16i5)') (nbond_move(i),i=1,Nbm)
893 write (iout,'(a)') 'Accepted motions:'
894 write (iout,'(16i5)') (nbond_acc(i),i=1,Nbm)
895 !el write (iout,'(a,i10)') 'Number of chain regrowths:',nregrow
896 !el write (iout,'(a,i10)') 'Accepted chain regrowths:',nregrow_acc
898 !---------------------------------------------------------------------------
900 !---------------------------------------------------------------------------
904 if (isweep.eq.nsweep .and. it.ge.maxacc) &
905 write (iout,'(/80(1h*)/20x,a/80(1h*)/)') 'All iterations done.'
907 end subroutine entmcm
908 !-----------------------------------------------------------------------------
909 subroutine accepting(ecur,eold,scur,sold,x,xold,accepted)
911 use geometry_data, only: nphi
912 use energy_data, only: max_ene
913 ! implicit real*8 (a-h,o-z)
914 ! include 'DIMENSIONS'
915 ! include 'COMMON.MCM'
916 ! include 'COMMON.MCE'
917 ! include 'COMMON.IOUNITS'
918 ! include 'COMMON.VAR'
920 use MPI_data !include 'COMMON.INFO'
922 ! include 'COMMON.GEO'
923 real(kind=8) :: ecur,eold,xx,bol !,ran_number
924 real(kind=8),dimension(6*nres) :: x,xold !(maxvar) (maxvar=6*maxres)
925 real(kind=8) :: tole=1.0D-1, tola=5.0D0
928 !!! local variables - el
930 real(kind=8) :: scur,sold,xxh,deix,dent
932 ! Check if the conformation is similar.
933 !d write (iout,*) 'Enter ACCEPTING'
934 !d write (iout,*) 'Old PHI angles:'
935 !d write (iout,*) (rad2deg*xold(i),i=1,nphi)
936 !d write (iout,*) 'Current angles'
937 !d write (iout,*) (rad2deg*x(i),i=1,nphi)
938 !d ddif=dif_ang(nphi,x,xold)
939 !d write (iout,*) 'Angle norm:',ddif
940 !d write (iout,*) 'ecur=',ecur,' emax=',emax
941 if (ecur.gt.emax) then
944 write (iout,'(a)') 'Conformation rejected as too high in energy'
946 else if (dabs(ecur-eold).lt.tole .and. &
947 dif_ang(nphi,x,xold).lt.tola) then
950 write (iout,'(a)') 'Conformation rejected as too similar'
953 ! Else evaluate the entropy of the conf and compare it with that of the previous
955 indecur=(ecur-emin)/delte
956 if (iabs(indecur).gt.max_ene) then
957 write (iout,'(a,2i5)') &
958 'Accepting: Index out of range:',indecur
960 else if (indecur.eq.indmaxx) then
961 scur=entropy(indecur)
962 if (print_mc.gt.0) write (iout,*)'Energy boundary reached',&
963 indmaxx,indecur,entropy(indecur)
965 deix=ecur-(emin+indecur*delte)
966 dent=entropy(indecur+1)-entropy(indecur)
967 scur=entropy(indecur)+(dent/delte)*deix
969 !d print *,'Processor',MyID,' ecur=',ecur,' indecur=',indecur,
970 !d & ' scur=',scur,' eold=',eold,' sold=',sold
971 !d print *,'deix=',deix,' dent=',dent,' delte=',delte
972 if (print_mc.gt.1) then
973 write(iout,*)'ecur=',ecur,' indecur=',indecur,' scur=',scur
974 write(iout,*)'eold=',eold,' sold=',sold
976 if (scur.le.sold) then
979 ! Else carry out acceptance test
980 xx=ran_number(0.0D0,1.0D0)
982 if (xxh.gt.50.0D0) then
989 if (print_mc.gt.0) write (iout,'(a)') &
990 'Conformation accepted.'
993 if (print_mc.gt.0) write (iout,'(a)') &
994 'Conformation rejected.'
998 end subroutine accepting
999 !-----------------------------------------------------------------------------
1000 subroutine read_pool
1002 use io_base, only:read_angles
1003 ! implicit real*8 (a-h,o-z)
1004 ! include 'DIMENSIONS'
1005 ! include 'COMMON.IOUNITS'
1006 ! include 'COMMON.GEO'
1007 ! include 'COMMON.MCM'
1008 ! include 'COMMON.MCE'
1009 ! include 'COMMON.VAR'
1010 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
1012 !!! local variables - el
1013 integer :: j,i,iconf
1015 print '(a)','Call READ_POOL'
1018 read (intin,'(i5,f10.5)',end=10,err=10) iconf,epool(npool)
1019 if (epool(npool).eq.0.0D0) goto 10
1020 call read_angles(intin,*10)
1021 call geom_to_var(nvar,xpool(1,npool))
1025 11 write (iout,'(a,i5)') 'Number of pool conformations:',npool
1026 if (print_mc.gt.2) then
1028 write (iout,'(a,i5,a,1pe14.5)') 'Pool conformation',i,' energy',&
1030 write (iout,'(10f8.3)') (rad2deg*xpool(j,i),j=1,nvar)
1032 endif ! (print_mc.gt.2)
1034 end subroutine read_pool
1035 !-----------------------------------------------------------------------------
1037 !-----------------------------------------------------------------------------
1038 subroutine monte_carlo
1042 use MPI_data, only:ifinish,nctasks,WhatsUp,MyID,NProcs
1043 use control_data, only:refstr !,MaxProcs
1045 use control, only:tcpu,ovrtim
1046 use regularize_, only:fitsq
1049 ! Does Boltzmann and entropic sampling without energy minimization
1050 ! implicit real*8 (a-h,o-z)
1051 ! include 'DIMENSIONS'
1052 ! include 'COMMON.IOUNITS'
1054 use MPI_data !include 'COMMON.INFO'
1056 ! include 'COMMON.GEO'
1057 ! include 'COMMON.CHAIN'
1058 ! include 'COMMON.MCM'
1059 ! include 'COMMON.MCE'
1060 ! include 'COMMON.CONTACTS'
1061 ! include 'COMMON.CONTROL'
1062 ! include 'COMMON.VAR'
1063 ! include 'COMMON.INTERACT'
1064 ! include 'COMMON.THREAD'
1065 ! include 'COMMON.NAMES'
1066 logical :: accepted,not_done,over,error,lprint !,ovrtim
1067 integer :: MoveType,nbond,nbins
1068 ! integer :: conf_comp
1069 real(kind=8) :: RandOrPert
1070 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
1071 real(kind=8) :: elowest,elowest1,ehighest,ehighest1,eold
1072 real(kind=8) :: przes(3),obr(3,3)
1073 real(kind=8),dimension(6*nres) :: varold !(maxvar) (maxvar=6*maxres)
1075 integer,dimension(-1:MaxMoveType+1,0:NProcs-1) :: moves1,moves_acc1 !(-1:MaxMoveType+1,0:MaxProcs-1)
1077 real(kind=8) :: etot_temp,etot_all(0:NProcs) !(0:MaxProcs)
1078 external d_vadd,d_vmin,d_vmax
1079 real(kind=8),dimension(-max_ene:max_ene) :: entropy1,nhist1
1080 integer,dimension(nres*(NProcs+1)) :: nbond_move1,nbond_acc1 !(nres*(MaxProcs+1))
1081 integer,dimension(2) :: itemp
1083 real(kind=8),dimension(6*nres) :: var_lowest !(maxvar) (maxvar=6*maxres)
1084 real(kind=8),dimension(0:n_ene) :: energia,energia_ave
1086 !!! local variables - el
1087 integer :: i,j,it,ii,iproc,nacc,ISWEEP,nfun,indmax,indmin,ijunk,&
1088 Kwita,indeold,imdmax,inde,iretcode,nstart_grow,noverlap
1089 real(kind=8) :: facee,conste,ejunk,etot,sold,frac,runtime,&
1090 frac_ave,rms_ave,etot_ave,scur,from_pool,co,rms
1092 write(iout,'(a,i8,2x,a,f10.5)') &
1093 'pool_read_freq=',pool_read_freq,' pool_fraction=',pool_fraction
1094 open (istat,file=statname)
1098 facee=1.0D0/(maxacc*delte)
1099 ! Number of bins in energy histogram
1101 write (iout,*) 'NBINS=',nbins
1103 ! Read entropy from previous simulations.
1105 read (ientin,*) indminn,indmaxx,emin,emax
1106 print *,'indminn=',indminn,' indmaxx=',indmaxx,' emin=',emin,&
1108 do i=-max_ene,max_ene
1111 read (ientin,*) (ijunk,ejunk,entropy(i),i=indminn,indmaxx)
1114 write (iout,*) 'indminn=',indminn,' indmaxx=',indmaxx,&
1115 ' emin=',emin,' emax=',emax
1116 write (iout,'(/a)') 'Initial entropy'
1117 do i=indminn,indmaxx
1118 write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
1121 ! Read the pool of conformations
1125 !----------------------------------------------------------------------------
1126 ! Entropy-sampling simulations with continually updated entropy;
1127 ! set NSWEEP=1 for Boltzmann sampling.
1128 ! Loop thru simulations
1129 !----------------------------------------------------------------------------
1130 allocate(ifinish(nctasks))
1133 ! Initialize the IFINISH array.
1140 !---------------------------------------------------------------------------
1141 ! Initialize counters.
1142 !---------------------------------------------------------------------------
1143 ! Total number of generated confs.
1145 ! Total number of moves. In general this won't be equal to the number of
1146 ! attempted moves, because we may want to reject some "bad" confs just by
1149 ! Total number of shift (nbond_move(1)), spike, crankshaft, three-bond,...
1151 !el allocate(nbond_move(nres)) !(maxres)
1152 !el allocate(nbond_acc(nres)) !(maxres)
1158 ! Initialize total and accepted number of moves of various kind.
1163 ! Total number of energy evaluations.
1166 !----------------------------------------------------------------------------
1167 ! Take a conformation from the pool
1168 !----------------------------------------------------------------------------
1170 write (iout,*) 'emin=',emin,' emax=',emax
1171 if (npool.gt.0) then
1172 ii=iran_num(1,npool)
1174 varia(i)=xpool(i,ii)
1176 write (iout,*) 'Took conformation',ii,' from the pool energy=',&
1178 call var_to_geom(nvar,varia)
1179 ! Print internal coordinates of the initial conformation
1181 else if (isweep.gt.1) then
1182 if (eold.lt.emax) then
1188 varia(i)=var_lowest(i)
1191 call var_to_geom(nvar,varia)
1193 !----------------------------------------------------------------------------
1194 ! Compute and print initial energies.
1195 !----------------------------------------------------------------------------
1199 write (iout,'(/80(1h*)/a,i2/80(1h*)/)') 'MCE iteration #',isweep
1200 write (iout,'(/80(1h*)/a)') 'Initial energies:'
1201 write(iout,*) 'Warning calling chainbuild'
1203 call geom_to_var(nvar,varia)
1204 call etotal(energia)
1206 call enerprint(energia)
1208 call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,1),nsup,przes,&
1211 call contact(.false.,ncont,icont,co)
1212 frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
1213 write (iout,'(a,f8.3,a,f8.3,a,f8.3)') &
1214 'RMS deviation from the reference structure:',rms,&
1215 ' % of native contacts:',frac*100,' contact order',co
1216 write (istat,'(i10,16(1pe14.5))') 0,&
1217 (energia(print_order(i)),i=1,nprint_ene),&
1220 write (istat,'(i10,14(1pe14.5))') 0,&
1221 (energia(print_order(i)),i=1,nprint_ene),etot
1225 if (.not. ent_read) then
1226 ! Initialize the entropy array
1228 ! Collect total energies from other processors.
1231 call mp_gather(etot_temp,etot_all,8,MasterID,cgGroupID)
1232 if (MyID.eq.MasterID) then
1233 ! Get the lowest and the highest energy.
1234 print *,'MASTER: etot_temp: ',(etot_all(i),i=0,nprocs-1),&
1235 ' emin=',emin,' emax=',emax
1239 if (emin.gt.etot_all(i)) emin=etot_all(i)
1240 if (emax.lt.etot_all(i)) emax=etot_all(i)
1243 endif ! MyID.eq.MasterID
1246 print *,'Processor',MyID,' calls MP_BCAST to send/recv etot_all'
1247 call mp_bcast(etot_all(1),16,MasterID,cgGroupID)
1248 print *,'Processor',MyID,' MP_BCAST to send/recv etot_all ended'
1249 if (MyID.ne.MasterID) then
1250 print *,'Processor:',MyID,etot_all(1),etot_all(2),&
1251 etot_all(1),etot_all(2)
1254 endif ! MyID.ne.MasterID
1255 write (iout,*) 'After MP_GATHER etot_temp=',&
1256 etot_temp,' emin=',emin
1264 ! Multicanonical sampling - start from Boltzmann distribution
1265 do i=-max_ene,max_ene
1266 entropy(i)=(emin+i*delte)*betbol
1269 ! Entropic sampling - start from uniform distribution of the density of states
1270 do i=-max_ene,max_ene
1274 write (iout,'(/a)') 'Initial entropy'
1275 do i=indminn,indmaxx
1276 write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
1278 if (isweep.eq.1) then
1282 indmaxx=indminn+nbins
1285 endif ! .not. ent_read
1287 call recv_stop_sig(Kwita)
1288 if (whatsup.eq.1) then
1289 call send_stop_sig(-2)
1291 else if (whatsup.le.-2) then
1293 else if (whatsup.eq.2) then
1301 write (iout,'(/80(1h*)/20x,a/80(1h*))') &
1302 'Enter Monte Carlo procedure.'
1304 call briefout(0,etot)
1309 call entropia(eold,sold,indeold)
1310 ! NACC is the counter for the accepted conformations of a given processor
1312 ! NACC_TOT counts the total number of accepted conformations
1315 !----------------------------------------------------------------------------
1316 ! Zero out average energies
1318 energia_ave(i)=0.0d0
1320 ! Initialize energy histogram
1321 do i=-max_ene,max_ene
1324 ! Zero out iteration counter.
1329 ! Begin MC iteration loop.
1332 ! Initialize local counter.
1333 ntrial=0 ! # of generated non-overlapping confs.
1334 noverlap=0 ! # of overlapping confs.
1336 do while (.not. accepted .and. WhatsUp.eq.0 .and. Kwita.eq.0)
1338 ! Retrieve the angles of previously accepted conformation
1342 call var_to_geom(nvar,varia)
1343 ! Rebuild the chain.
1344 write(iout,*) 'Warning calling chainbuild'
1349 ! Decide whether to take a conformation from the pool or generate/perturb one
1351 from_pool=ran_number(0.0D0,1.0D0)
1352 if (npool.gt.0 .and. from_pool.lt.pool_fraction) then
1353 ! Throw a dice to choose the conformation from the pool
1354 ii=iran_num(1,npool)
1356 varia(i)=xpool(i,ii)
1358 call var_to_geom(nvar,varia)
1359 write(iout,*) 'Warning calling chainbuild'
1362 !d write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
1363 if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) &
1364 write (iout,'(a,i3,a,f10.5)') &
1365 'Try conformation',ii,' from the pool energy=',epool(ii)
1367 moves(-1)=moves(-1)+1
1369 ! Decide whether to generate a random conformation or perturb the old one
1370 RandOrPert=ran_number(0.0D0,1.0D0)
1371 if (RandOrPert.gt.RanFract) then
1372 if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) &
1373 write (iout,'(a)') 'Perturbation-generated conformation.'
1374 call perturb(error,lprint,MoveType,nbond,0.1D0)
1376 if (MoveType.lt.1 .or. MoveType.gt.MaxMoveType) then
1377 write (iout,'(/a,i7,a/)') 'Error - unknown MoveType=',&
1378 MoveType,' returned from PERTURB.'
1385 nstart_grow=iran_num(3,nres)
1386 if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) &
1387 write (iout,'(2a,i3)') 'Random-generated conformation',&
1388 ' - chain regrown from residue',nstart_grow
1389 call gen_rand_conf(nstart_grow,*30)
1391 call geom_to_var(nvar,varia)
1393 !d write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
1395 if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) &
1396 write (iout,'(a,i5,a,i10,a,i10)') &
1397 'Processor',MyId,' trial move',ntrial,' total generated:',ngen
1398 if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) &
1399 write (*,'(a,i5,a,i10,a,i10)') &
1400 'Processor',MyId,' trial move',ntrial,' total generated:',ngen
1401 call etotal(energia)
1404 !d call enerprint(energia(0))
1405 !d write(iout,*)'it=',it,' etot=',etot
1406 if (etot-elowest.gt.overlap_cut) then
1407 if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) &
1408 write (iout,'(a,i5,a,1pe14.5)') 'Iteration',it,&
1409 ' Overlap detected in the current conf.; energy is',etot
1412 if (noverlap.gt.maxoverlap) then
1413 write (iout,'(a)') 'Too many overlapping confs.'
1417 !--------------------------------------------------------------------------
1418 !... Acceptance test
1419 !--------------------------------------------------------------------------
1422 call accept_mc(it,etot,eold,scur,sold,varia,varold,accepted)
1426 if (elowest.gt.etot) then
1429 var_lowest(i)=varia(i)
1432 if (ehighest.lt.etot) ehighest=etot
1433 moves_acc(MoveType)=moves_acc(MoveType)+1
1434 if (MoveType.eq.1) then
1435 nbond_acc(nbond)=nbond_acc(nbond)+1
1437 ! Compare with reference structure.
1439 call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,1),&
1440 nsup,przes,obr,non_conv)
1442 call contact(.false.,ncont,icont,co)
1443 frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
1446 ! Periodically save average energies and confs.
1449 energia_ave(i)=energia_ave(i)+energia(i)
1451 moves(MaxMoveType+1)=nmove
1452 moves_acc(MaxMoveType+1)=nacc
1453 IF ((it/save_frequency)*save_frequency.eq.it) THEN
1455 energia_ave(i)=energia_ave(i)/save_frequency
1457 etot_ave=energia_ave(0)
1459 ! open (istat,file=statname,position='append')
1461 ! open (istat,file=statname,access='append')
1463 if (print_mc.gt.0) &
1464 write (iout,'(80(1h*)/20x,a,i20)') &
1466 if (refstr .and. print_mc.gt.0) then
1467 write (iout,'(a,f8.3,a,f8.3,a,f8.3)') &
1468 'RMS deviation from the reference structure:',rms,&
1469 ' % of native contacts:',frac*100,' contact order:',co
1471 if (print_stat) then
1473 write (istat,'(i10,10(1pe14.5))') it,&
1474 (energia_ave(print_order(i)),i=1,nprint_ene),&
1475 etot_ave,rms_ave,frac_ave
1477 write (istat,'(i10,10(1pe14.5))') it,&
1478 (energia_ave(print_order(i)),i=1,nprint_ene),&
1483 if (print_mc.gt.0) &
1484 call statprint(nacc,nfun,iretcode,etot,elowest)
1485 ! Print internal coordinates.
1486 if (print_int) call briefout(nacc,etot)
1488 energia_ave(i)=0.0d0
1490 ENDIF ! ( (it/save_frequency)*save_frequency.eq.it)
1492 inde=icialosc((etot-emin)/delte)
1493 nhist(inde)=nhist(inde)+1.0D0
1495 if ( (it/message_frequency)*message_frequency.eq.it &
1496 .and. (MyID.ne.MasterID) ) then
1497 call recv_stop_sig(Kwita)
1498 call send_MCM_info(message_frequency)
1501 ! Store the accepted conf. and its energy.
1508 if (Kwita.eq.0) call recv_stop_sig(kwita)
1513 if (MyID.eq.MasterID .and. &
1514 (it/message_frequency)*message_frequency.eq.it) then
1515 call receive_MC_info
1516 if (nacc_tot.ge.maxacc) accepted=.true.
1519 ! if ((ntrial.gt.maxtrial_iter
1520 ! & .or. (it/pool_read_freq)*pool_read_freq.eq.it)
1521 ! & .and. npool.gt.0) then
1522 ! Take a conformation from the pool
1523 ! ii=iran_num(1,npool)
1525 ! varold(i)=xpool(i,ii)
1527 ! if (ntrial.gt.maxtrial_iter)
1528 ! & write (iout,*) 'Iteration',it,' max. # of trials exceeded.'
1530 ! & 'Take conformation',ii,' from the pool energy=',epool(ii)
1531 ! if (print_mc.gt.2)
1532 ! & write (iout,'(10f8.3)') (rad2deg*varold(i),i=1,nvar)
1535 ! call entropia(eold,sold,indeold)
1537 ! endif ! (ntrial.gt.maxtrial_iter .and. npool.gt.0)
1541 if (MyID.eq.MasterID .and. &
1542 (it/message_frequency)*message_frequency.eq.it) then
1543 call receive_MC_info
1545 if (Kwita.eq.0) call recv_stop_sig(kwita)
1547 if (ovrtim()) WhatsUp=-1
1548 !d write (iout,*) 'WhatsUp=',WhatsUp,' Kwita=',Kwita
1549 not_done = (nacc_tot.lt.maxacc) .and. (WhatsUp.eq.0) &
1551 !d write (iout,*) 'not_done=',not_done
1553 if (Kwita.lt.0) then
1554 print *,'Processor',MyID,&
1555 ' has received STOP signal =',Kwita,' in EntSamp.'
1556 !d print *,'not_done=',not_done
1557 if (Kwita.lt.-1) WhatsUp=Kwita
1558 if (MyID.ne.MasterID) call send_MCM_info(-1)
1559 else if (nacc_tot.ge.maxacc) then
1560 print *,'Processor',MyID,' calls send_stop_sig,',&
1561 ' because a sufficient # of confs. have been collected.'
1562 !d print *,'not_done=',not_done
1563 call send_stop_sig(-1)
1564 if (MyID.ne.MasterID) call send_MCM_info(-1)
1565 else if (WhatsUp.eq.-1) then
1566 print *,'Processor',MyID,&
1567 ' calls send_stop_sig because of timeout.'
1568 !d print *,'not_done=',not_done
1569 call send_stop_sig(-2)
1570 if (MyID.ne.MasterID) call send_MCM_info(-1)
1575 !-----------------------------------------------------------------
1576 !... Construct energy histogram & update entropy
1577 !-----------------------------------------------------------------
1581 write (iout,*) 'Processor',MyID,&
1582 ' is broadcasting ERROR-STOP signal.'
1583 write (*,*) 'Processor',MyID,&
1584 ' is broadcasting ERROR-STOP signal.'
1585 call send_stop_sig(-3)
1586 if (MyID.ne.MasterID) call send_MCM_info(-1)
1589 write (iout,'(/a)') 'Energy histogram'
1591 write (iout,'(i5,2f20.5)') i,emin+i*delte,nhist(i)
1594 ! Wait until every processor has sent complete MC info.
1595 if (MyID.eq.MasterID) then
1598 ! write (*,*) 'The IFINISH array:'
1599 ! write (*,*) (ifinish(i),i=1,nctasks)
1602 not_done=not_done.or.(ifinish(i).ge.0)
1604 if (not_done) call receive_MC_info
1607 ! Make collective histogram from the work of all processors.
1608 msglen=(2*max_ene+1)*8
1610 'Processor',MyID,' calls MP_REDUCE to send/receive histograms',&
1612 call mp_reduce(nhist,nhist1,msglen,MasterID,d_vadd,&
1614 print *,'Processor',MyID,' MP_REDUCE accomplished for histogr.'
1615 do i=-max_ene,max_ene
1618 ! Collect min. and max. energy
1620 'Processor',MyID,' calls MP_REDUCE to send/receive energy borders'
1621 call mp_reduce(elowest,elowest1,8,MasterID,d_vmin,cgGroupID)
1622 call mp_reduce(ehighest,ehighest1,8,MasterID,d_vmax,cgGroupID)
1623 print *,'Processor',MyID,' MP_REDUCE accomplished for energies.'
1624 IF (MyID.eq.MasterID) THEN
1628 write (iout,'(a,i10)') '# of accepted confs:',nacc_tot
1629 write (iout,'(a,f10.5,a,f10.5)') 'Lowest energy:',elowest,&
1630 ' Highest energy',ehighest
1631 indmin=icialosc((elowest-emin)/delte)
1632 imdmax=icialosc((ehighest-emin)/delte)
1633 if (indmin.lt.indminn) then
1634 emax=emin+indmin*delte+e_up
1635 indmaxx=indmin+nbins
1638 if (.not.ent_read) ent_read=.true.
1639 write(iout,*)'indminn=',indminn,' indmaxx=',indmaxx
1640 ! Update entropy (density of states)
1642 if (nhist(i).gt.0) then
1643 entropy(i)=entropy(i)+dlog(nhist(i)+0.0D0)
1646 write (iout,'(/80(1h*)/a,i2/80(1h*)/)') &
1647 'End of macroiteration',isweep
1648 write (iout,'(a,f10.5,a,f10.5)') 'Elowest=',elowest,&
1649 ' Ehighest=',ehighest
1650 write (iout,'(/a)') 'Energy histogram'
1651 do i=indminn,indmaxx
1652 write (iout,'(i5,2f20.5)') i,emin+i*delte,nhist(i)
1654 write (iout,'(/a)') 'Entropy'
1655 do i=indminn,indmaxx
1656 write (iout,'(i5,2f20.5)') i,emin+i*delte,entropy(i)
1658 !-----------------------------------------------------------------
1659 !... End of energy histogram construction
1660 !-----------------------------------------------------------------
1663 if (.not. ent_read) ent_read=.true.
1664 ENDIF ! MyID .eq. MaterID
1665 if (MyID.eq.MasterID) then
1669 print *,'before mp_bcast processor',MyID,' indminn=',indminn,&
1670 ' indmaxx=',indmaxx,' itemp=',itemp(1),itemp(2)
1671 call mp_bcast(itemp(1),8,MasterID,cgGroupID)
1672 call mp_bcast(emax,8,MasterID,cgGroupID)
1673 print *,'after mp_bcast processor',MyID,' indminn=',indminn,&
1674 ' indmaxx=',indmaxx,' itemp=',itemp(1),itemp(2)
1675 if (MyID .ne. MasterID) then
1679 msglen=(indmaxx-indminn+1)*8
1680 print *,'processor',MyID,' calling mp_bcast msglen=',msglen,&
1681 ' indminn=',indminn,' indmaxx=',indmaxx,' isweep=',isweep
1682 call mp_bcast(entropy(indminn),msglen,MasterID,cgGroupID)
1683 IF(MyID.eq.MasterID .and. .not. ovrtim() .and. WhatsUp.ge.0)THEN
1684 open (ientout,file=entname,status='unknown')
1685 write (ientout,'(2i5,2e25.17)') indminn,indmaxx,emin,emax
1686 do i=indminn,indmaxx
1687 write (ientout,'(i5,f10.5,f20.15)') i,emin+i*delte,entropy(i)
1691 write (iout,*) 'Received from master:'
1692 write (iout,*) 'indminn=',indminn,' indmaxx=',indmaxx,&
1693 ' emin=',emin,' emax=',emax
1694 write (iout,'(/a)') 'Entropy'
1695 do i=indminn,indmaxx
1696 write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
1698 ENDIF ! MyID.eq.MasterID
1699 print *,'Processor',MyID,' calls MP_GATHER'
1700 call mp_gather(nbond_move,nbond_move1,4*Nbm,MasterID,&
1702 call mp_gather(nbond_acc,nbond_acc1,4*Nbm,MasterID,&
1704 print *,'Processor',MyID,' MP_GATHER call accomplished'
1705 if (MyID.eq.MasterID) then
1707 write (iout,'(/80(1h*)/20x,a)') 'Summary run statistics:'
1708 call statprint(nacc_tot,nfun,iretcode,etot,elowest)
1709 write (iout,'(a)') &
1710 'Statistics of multiple-bond motions. Total motions:'
1711 write (iout,'(8i10)') (nbond_move(i),i=1,Nbm)
1712 write (iout,'(a)') 'Accepted motions:'
1713 write (iout,'(8i10)') (nbond_acc(i),i=1,Nbm)
1715 write (iout,'(a)') &
1716 'Statistics of multi-bond moves of respective processors:'
1720 nbond_move(i)=nbond_move(i)+nbond_move1(ind)
1721 nbond_acc(i)=nbond_acc(i)+nbond_acc1(ind)
1725 write (iout,*) 'Processor',iproc,' nbond_move:', &
1726 (nbond_move1(iproc*nbm+i),i=1,Nbm),&
1727 ' nbond_acc:',(nbond_acc1(iproc*nbm+i),i=1,Nbm)
1730 call mp_gather(moves,moves1,4*(MaxMoveType+3),MasterID,&
1732 call mp_gather(moves_acc,moves_acc1,4*(MaxMoveType+3),&
1734 if (MyID.eq.MasterID) then
1736 do i=-1,MaxMoveType+1
1737 moves(i)=moves(i)+moves1(i,iproc)
1738 moves_acc(i)=moves_acc(i)+moves_acc1(i,iproc)
1742 do i=0,MaxMoveType+1
1743 nmove=nmove+moves(i)
1746 write (iout,*) 'Processor',iproc,' moves',&
1747 (moves1(i,iproc),i=0,MaxMoveType+1),&
1748 ' moves_acc:',(moves_acc1(i,iproc),i=0,MaxMoveType+1)
1752 open (ientout,file=entname,status='unknown')
1753 write (ientout,'(2i5,2e25.17)') indminn,indmaxx,emin,emax
1754 do i=indminn,indmaxx
1755 write (ientout,'(i5,f10.5,f20.15)') i,emin+i*delte,entropy(i)
1759 write (iout,'(/80(1h*)/20x,a)') 'Summary run statistics:'
1760 call statprint(nacc_tot,nfun,iretcode,etot,elowest)
1761 write (iout,'(a)') &
1762 'Statistics of multiple-bond motions. Total motions:'
1763 write (iout,'(8i10)') (nbond_move(i),i=1,Nbm)
1764 write (iout,'(a)') 'Accepted motions:'
1765 write (iout,'(8i10)') (nbond_acc(i),i=1,Nbm)
1766 if (ovrtim() .or. WhatsUp.lt.0) return
1768 !---------------------------------------------------------------------------
1770 !---------------------------------------------------------------------------
1774 if (isweep.eq.nsweep .and. it.ge.maxacc) &
1775 write (iout,'(/80(1h*)/20x,a/80(1h*)/)') 'All iterations done.'
1777 end subroutine monte_carlo
1778 !-----------------------------------------------------------------------------
1779 subroutine accept_mc(it,ecur,eold,scur,sold,x,xold,accepted)
1781 ! implicit real*8 (a-h,o-z)
1782 ! include 'DIMENSIONS'
1783 ! include 'COMMON.MCM'
1784 ! include 'COMMON.MCE'
1785 ! include 'COMMON.IOUNITS'
1786 ! include 'COMMON.VAR'
1788 use MPI_data !include 'COMMON.INFO'
1790 ! include 'COMMON.GEO'
1791 real(kind=8) :: ecur,eold,xx,bol
1792 real(kind=8),dimension(6*nres) :: x,xold !(maxvar) (maxvar=6*maxres)
1796 integer :: it,indecur
1797 real(kind=8) :: scur,sold,xxh
1798 ! Check if the conformation is similar.
1799 !d write (iout,*) 'Enter ACCEPTING'
1800 !d write (iout,*) 'Old PHI angles:'
1801 !d write (iout,*) (rad2deg*xold(i),i=1,nphi)
1802 !d write (iout,*) 'Current angles'
1803 !d write (iout,*) (rad2deg*x(i),i=1,nphi)
1804 !d ddif=dif_ang(nphi,x,xold)
1805 !d write (iout,*) 'Angle norm:',ddif
1806 !d write (iout,*) 'ecur=',ecur,' emax=',emax
1807 if (ecur.gt.emax) then
1809 if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) &
1810 write (iout,'(a)') 'Conformation rejected as too high in energy'
1813 ! Else evaluate the entropy of the conf and compare it with that of the previous
1815 call entropia(ecur,scur,indecur)
1816 !d print *,'Processor',MyID,' ecur=',ecur,' indecur=',indecur,
1817 !d & ' scur=',scur,' eold=',eold,' sold=',sold
1818 !d print *,'deix=',deix,' dent=',dent,' delte=',delte
1819 if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) then
1820 write(iout,*)'it=',it,'ecur=',ecur,' indecur=',indecur,&
1822 write(iout,*)'eold=',eold,' sold=',sold
1824 if (scur.le.sold) then
1827 ! Else carry out acceptance test
1828 xx=ran_number(0.0D0,1.0D0)
1830 if (xxh.gt.50.0D0) then
1837 if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) &
1838 write (iout,'(a)') 'Conformation accepted.'
1841 if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) &
1842 write (iout,'(a)') 'Conformation rejected.'
1846 end subroutine accept_mc
1847 !-----------------------------------------------------------------------------
1848 integer function icialosc(x)
1851 if (x.lt.0.0D0) then
1857 end function icialosc
1858 !-----------------------------------------------------------------------------
1859 subroutine entropia(ecur,scur,indecur)
1861 use energy_data, only: max_ene
1862 ! implicit real*8 (a-h,o-z)
1863 ! include 'DIMENSIONS'
1864 ! include 'COMMON.MCM'
1865 ! include 'COMMON.MCE'
1866 ! include 'COMMON.IOUNITS'
1867 real(kind=8) :: ecur,scur,deix,dent
1868 integer :: indecur,it !???el
1870 indecur=icialosc((ecur-emin)/delte)
1871 if (iabs(indecur).gt.max_ene) then
1872 if ((it/print_freq)*it.eq.it) write (iout,'(a,2i5)') &
1873 'Accepting: Index out of range:',indecur
1875 else if (indecur.ge.indmaxx) then
1876 scur=entropy(indecur)
1877 if (print_mc.gt.0 .and. (it/print_freq)*it.eq.it) &
1878 write (iout,*)'Energy boundary reached',&
1879 indmaxx,indecur,entropy(indecur)
1881 deix=ecur-(emin+indecur*delte)
1882 dent=entropy(indecur+1)-entropy(indecur)
1883 scur=entropy(indecur)+(dent/delte)*deix
1886 end subroutine entropia
1887 !-----------------------------------------------------------------------------
1889 !-----------------------------------------------------------------------------
1890 subroutine mcm_setup
1893 ! implicit real*8 (a-h,o-z)
1894 ! include 'DIMENSIONS'
1895 ! include 'COMMON.IOUNITS'
1896 ! include 'COMMON.MCM'
1897 ! include 'COMMON.CONTROL'
1898 ! include 'COMMON.INTERACT'
1899 ! include 'COMMON.NAMES'
1900 ! include 'COMMON.CHAIN'
1901 ! include 'COMMON.VAR'
1903 !!! local variables - el
1904 integer :: i,i1,i2,it1,it2,ngly,mmm,maxwinlen
1906 ! Set up variables used in MC/MCM.
1908 ! allocate(sumpro_bond(0:nres)) !(0:maxres)
1910 write (iout,'(80(1h*)/20x,a/80(1h*))') 'MCM control parameters:'
1911 write (iout,'(5(a,i7))') 'Maxacc:',maxacc,' MaxTrial:',MaxTrial,&
1912 ' MaxRepm:',MaxRepm,' MaxGen:',MaxGen,' MaxOverlap:',MaxOverlap
1913 write (iout,'(4(a,f8.1)/2(a,i3))') &
1914 'Tmin:',Tmin,' Tmax:',Tmax,' TstepH:',TstepH,&
1915 ' TstepC:',TstepC,'NstepH:',NstepH,' NstepC:',NstepC
1916 if (nwindow.gt.0) then
1917 write (iout,'(a)') 'Perturbation windows:'
1923 write (iout,'(a,i3,a,i3,a,i3)') restyp(it1,1),i1,restyp(it2,1),i2,&
1927 ! Rbolt=8.3143D-3*2.388459D-01 kcal/(mol*K)
1929 ! Number of "end bonds".
1932 print *,'koniecl=',koniecl
1933 write (iout,'(a)') 'Probabilities of move types:'
1934 write (*,'(a)') 'Probabilities of move types:'
1936 write (iout,'(a,f10.5)') MovTypID(i),&
1937 sumpro_type(i)-sumpro_type(i-1)
1938 write (*,'(a,f10.5)') MovTypID(i),&
1939 sumpro_type(i)-sumpro_type(i-1)
1942 ! Maximum length of N-bond segment to be moved
1943 ! nbm=nres-1-(2*koniecl-1)
1944 if (nwindow.gt.0) then
1947 if (winlen(i).gt.maxwinlen) maxwinlen=winlen(i)
1949 nbm=min0(maxwinlen,6)
1950 write (iout,'(a,i3,a,i3)') 'Nbm=',Nbm,' Maxwinlen=',Maxwinlen
1954 sumpro_bond(0)=0.0D0
1955 sumpro_bond(1)=0.0D0
1957 sumpro_bond(i)=sumpro_bond(i-1)+1.0D0/dfloat(i)
1959 write (iout,'(a)') 'The SumPro_Bond array:'
1960 write (iout,'(8f10.5)') (sumpro_bond(i),i=1,nbm)
1961 write (*,'(a)') 'The SumPro_Bond array:'
1962 write (*,'(8f10.5)') (sumpro_bond(i),i=1,nbm)
1963 ! Maximum number of side chains moved simultaneously
1964 ! print *,'nnt=',nnt,' nct=',nct
1967 if (itype(i,1).eq.10) ngly=ngly+1
1971 MaxSideMove=min0((nct-nnt+1)/2,mmm)
1973 ! print *,'MaxSideMove=',MaxSideMove
1974 ! Max. number of generated confs (not used at present).
1976 ! Set initial temperature
1978 betbol=1.0D0/(Rbol*Tcur)
1979 write (iout,'(a,f8.1,a,f10.5)') 'Initial temperature:',Tcur,&
1981 write (iout,*) 'RanFract=',ranfract
1983 end subroutine mcm_setup
1984 !-----------------------------------------------------------------------------
1986 subroutine do_mcm(i_orig)
1990 use MPI_data, only:Whatsup
1991 use control_data, only:refstr,minim,iprint
1993 use control, only:tcpu,ovrtim
1994 use regularize_, only:fitsq
1996 use minimm, only:minimize
1997 ! Monte-Carlo-with-Minimization calculations - serial code.
1998 ! implicit real*8 (a-h,o-z)
1999 ! include 'DIMENSIONS'
2000 ! include 'COMMON.IOUNITS'
2001 ! include 'COMMON.GEO'
2002 ! include 'COMMON.CHAIN'
2003 ! include 'COMMON.MCM'
2004 ! include 'COMMON.CONTACTS'
2005 ! include 'COMMON.CONTROL'
2006 ! include 'COMMON.VAR'
2007 ! include 'COMMON.INTERACT'
2008 ! include 'COMMON.CACHE'
2009 !rc include 'COMMON.DEFORM'
2010 !rc include 'COMMON.DEFORM1'
2011 ! include 'COMMON.NAMES'
2012 logical :: accepted,over,error,lprint,not_done,my_conf,&
2013 enelower,non_conv !,ovrtim
2014 integer :: MoveType,nbond !,conf_comp
2015 integer,dimension(max_cache) :: ifeed
2016 real(kind=8),dimension(6*nres) :: varia,varold !(maxvar) (maxvar=6*maxres)
2017 real(kind=8) :: elowest,eold
2018 real(kind=8) :: przes(3),obr(3,3)
2019 real(kind=8) :: energia(0:n_ene)
2020 real(kind=8) :: coord1(6*nres,max_thread2),enetb1(max_threadss) !el
2021 !!! local variables - el
2022 integer :: i,nf,nacc,it,nout,j,i_orig,nfun,Kwita,iretcode,&
2023 noverlap,nstart_grow,irepet,n_thr,ii
2024 real(kind=8) :: etot,frac,rms,co,RandOrPert,&
2026 !---------------------------------------------------------------------------
2027 ! Initialize counters.
2028 !---------------------------------------------------------------------------
2029 ! Total number of generated confs.
2031 ! Total number of moves. In general this won't be equal to the number of
2032 ! attempted moves, because we may want to reject some "bad" confs just by
2035 ! Total number of temperature jumps.
2037 ! Total number of shift (nbond_move(1)), spike, crankshaft, three-bond,...
2039 ! if(.not.allocated(varsave)) allocate(varsave(maxvar,maxsave)) !(maxvar,maxsave)
2040 ! allocate(nbond_move(nres)) !(maxres)
2046 ! Initialize total and accepted number of moves of various kind.
2051 ! Total number of energy evaluations.
2056 write (iout,*) 'RanFract=',RanFract
2061 !----------------------------------------------------------------------------
2062 ! Compute and print initial energies.
2063 !----------------------------------------------------------------------------
2065 write (iout,'(/80(1h*)/a)') 'Initial energies:'
2069 call etotal(energia)
2071 ! Minimize the energy of the first conformation.
2073 call geom_to_var(nvar,varia)
2074 ! write (iout,*) 'The VARIA array'
2075 ! write (iout,'(8f10.4)') (rad2deg*varia(i),i=1,nvar)
2076 call minimize(etot,varia,iretcode,nfun)
2077 call var_to_geom(nvar,varia)
2079 write (iout,*) 'etot from MINIMIZE:',etot
2080 ! write (iout,*) 'Tha VARIA array'
2081 ! write (iout,'(8f10.4)') (rad2deg*varia(i),i=1,nvar)
2083 call etotal(energia)
2085 call enerprint(energia)
2088 call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,1),nsup,przes,& !el cref(1,nstart_sup)
2091 call contact(.false.,ncont,icont,co)
2092 frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
2093 write (iout,'(a,f8.3,a,f8.3,a,f8.3)') &
2094 'RMS deviation from the reference structure:',rms,&
2095 ' % of native contacts:',frac*100,' contact order:',co
2097 write (istat,'(i5,17(1pe14.5))') 0,&
2098 (energia(print_order(i)),i=1,nprint_ene),&
2101 if (print_stat) write (istat,'(i5,16(1pe14.5))') 0,&
2102 (energia(print_order(i)),i=1,nprint_ene),etot
2104 if (print_stat) close(istat)
2105 neneval=neneval+nfun+1
2106 write (iout,'(/80(1h*)/20x,a/80(1h*))') &
2107 'Enter Monte Carlo procedure.'
2110 call briefout(0,etot)
2117 call zapis(varia,etot)
2118 nacc=0 ! total # of accepted confs of the current processor.
2119 nacc_tot=0 ! total # of accepted confs of all processors.
2121 not_done = (iretcode.ne.11)
2123 !----------------------------------------------------------------------------
2125 !----------------------------------------------------------------------------
2130 write (iout,'(80(1h*)/20x,a,i7)') &
2131 'Beginning iteration #',it
2132 ! Initialize local counter.
2133 ntrial=0 ! # of generated non-overlapping confs.
2135 do while (.not. accepted)
2137 ! Retrieve the angles of previously accepted conformation
2138 noverlap=0 ! # of overlapping confs.
2142 call var_to_geom(nvar,varia)
2143 ! Rebuild the chain.
2145 ! Heat up the system, if necessary.
2147 ! If temperature cannot be further increased, stop.
2152 !d write (iout,'(a)') 'Old variables:'
2153 !d write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
2154 ! Decide whether to generate a random conformation or perturb the old one
2155 RandOrPert=ran_number(0.0D0,1.0D0)
2156 if (RandOrPert.gt.RanFract) then
2157 if (print_mc.gt.0) &
2158 write (iout,'(a)') 'Perturbation-generated conformation.'
2159 call perturb(error,lprint,MoveType,nbond,1.0D0)
2161 if (MoveType.lt.1 .or. MoveType.gt.MaxMoveType) then
2162 write (iout,'(/a,i7,a/)') 'Error - unknown MoveType=',&
2163 MoveType,' returned from PERTURB.'
2170 nstart_grow=iran_num(3,nres)
2171 if (print_mc.gt.0) &
2172 write (iout,'(2a,i3)') 'Random-generated conformation',&
2173 ' - chain regrown from residue',nstart_grow
2174 call gen_rand_conf(nstart_grow,*30)
2176 call geom_to_var(nvar,varia)
2177 !d write (iout,'(a)') 'New variables:'
2178 !d write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
2181 call etotal(energia)
2183 ! call enerprint(energia(0))
2184 ! write (iout,'(2(a,1pe14.5))') 'Etot=',Etot,' Elowest=',Elowest
2185 if (etot-elowest.gt.overlap_cut) then
2186 if(iprint.gt.1.or.etot.lt.1d20) &
2187 write (iout,'(a,1pe14.5)') &
2188 'Overlap detected in the current conf.; energy is',etot
2192 if (noverlap.gt.maxoverlap) then
2193 write (iout,'(a)') 'Too many overlapping confs.'
2198 call minimize(etot,varia,iretcode,nfun)
2199 !d write (iout,*) 'etot from MINIMIZE:',etot
2200 !d write (iout,'(a)') 'Variables after minimization:'
2201 !d write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
2203 call etotal(energia)
2205 neneval=neneval+nfun+2
2207 ! call enerprint(energia(0))
2208 write (iout,'(a,i6,a,1pe16.6)') 'Conformation:',ngen,&
2210 !--------------------------------------------------------------------------
2211 !... Do Metropolis test
2212 !--------------------------------------------------------------------------
2216 if (WhatsUp.eq.0 .and. Kwita.eq.0) then
2217 call metropolis(nvar,varia,varold,etot,eold,accepted,&
2220 write (iout,*) 'My_Conf=',My_Conf,' EneLower=',EneLower
2225 if (elowest.gt.etot) elowest=etot
2226 moves_acc(MoveType)=moves_acc(MoveType)+1
2227 if (MoveType.eq.1) then
2228 nbond_acc(nbond)=nbond_acc(nbond)+1
2230 ! Check against conformation repetitions.
2231 irepet=conf_comp(varia,etot)
2232 if (print_stat) then
2233 #if defined(AIX) || defined(PGI)
2234 open (istat,file=statname,position='append')
2236 open (istat,file=statname,access='append')
2239 call statprint(nacc,nfun,iretcode,etot,elowest)
2241 call var_to_geom(nvar,varia)
2243 call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,1),& !el cref(1,nstart_sup)
2244 nsup,przes,obr,non_conv)
2246 call contact(.false.,ncont,icont,co)
2247 frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
2248 write (iout,'(a,f8.3,a,f8.3)') &
2249 'RMS deviation from the reference structure:',rms,&
2250 ' % of native contacts:',frac*100,' contact order',co
2254 write (iout,*) 'Writing new conformation',nout
2256 write (istat,'(i5,16(1pe14.5))') nout,&
2257 (energia(print_order(i)),i=1,nprint_ene),&
2261 write (istat,'(i5,17(1pe14.5))') nout,&
2262 (energia(print_order(i)),i=1,nprint_ene),etot
2264 if (print_stat) close(istat)
2265 ! Print internal coordinates.
2266 if (print_int) call briefout(nout,etot)
2267 ! Accumulate the newly accepted conf in the coord1 array, if it is different
2268 ! from all confs that are already there.
2269 call compare_s1(n_thr,max_thread2,etot,varia,ii,&
2270 enetb1,coord1,rms_deform,.true.,iprint)
2271 write (iout,*) 'After compare_ss: n_thr=',n_thr
2272 if (ii.eq.1 .or. ii.eq.3) then
2273 write (iout,'(8f10.4)') &
2274 (rad2deg*coord1(i,n_thr),i=1,nvar)
2277 write (iout,*) 'Conformation from cache, not written.'
2280 if (nrepm.gt.maxrepm) then
2281 write (iout,'(a)') 'Too many conformation repetitions.'
2284 ! Store the accepted conf. and its energy.
2289 if (irepet.eq.0) call zapis(varia,etot)
2290 ! Lower the temperature, if necessary.
2301 ! Check for time limit.
2302 if (ovrtim()) WhatsUp=-1
2303 not_done = (nacc_tot.lt.maxacc) .and. (WhatsUp.eq.0) &
2312 write (iout,'(/80(1h*)/20x,a)') 'Summary run statistics:'
2313 call statprint(nacc,nfun,iretcode,etot,elowest)
2314 write (iout,'(a)') &
2315 'Statistics of multiple-bond motions. Total motions:'
2316 write (iout,'(16i5)') (nbond_move(i),i=1,Nbm)
2317 write (iout,'(a)') 'Accepted motions:'
2318 write (iout,'(16i5)') (nbond_acc(i),i=1,Nbm)
2320 write (iout,'(/80(1h*)/20x,a/80(1h*)/)') 'All iterations done.'
2321 !(maxvar) (maxvar=6*maxres)
2323 end subroutine do_mcm
2325 !-----------------------------------------------------------------------------
2327 subroutine do_mcm(i_orig)
2329 ! Monte-Carlo-with-Minimization calculations - parallel code.
2331 use control_data, only:refstr!,tag
2332 use io_base, only:intout,briefout
2333 use control, only:ovrtim,tcpu
2334 use compare, only:contact,contact_fract
2335 use minimm, only:minimize
2336 use regularize_, only:fitsq
2337 ! use contact_, only:contact
2339 ! implicit real*8 (a-h,o-z)
2340 ! include 'DIMENSIONS'
2342 ! include 'COMMON.IOUNITS'
2343 ! include 'COMMON.GEO'
2344 ! include 'COMMON.CHAIN'
2345 ! include 'COMMON.MCM'
2346 ! include 'COMMON.CONTACTS'
2347 ! include 'COMMON.CONTROL'
2348 ! include 'COMMON.VAR'
2349 ! include 'COMMON.INTERACT'
2350 ! include 'COMMON.INFO'
2351 ! include 'COMMON.CACHE'
2352 !rc include 'COMMON.DEFORM'
2353 !rc include 'COMMON.DEFORM1'
2354 !rc include 'COMMON.DEFORM2'
2355 ! include 'COMMON.MINIM'
2356 ! include 'COMMON.NAMES'
2357 logical :: accepted,over,error,lprint,not_done,similar,&
2358 enelower,non_conv,flag,finish !,ovrtim
2359 integer :: MoveType,nbond !,conf_comp
2360 real(kind=8),dimension(6*nres) :: x1,varold1,varold,varia !(maxvar) (maxvar=6*maxres)
2361 real(kind=8) :: elowest,eold
2362 real(kind=8) :: przes(3),obr(3,3)
2363 integer :: iparentx(max_threadss2)
2364 integer :: iparentx1(max_threadss2)
2365 integer :: imtasks(150),imtasks_n
2366 real(kind=8) :: energia(0:n_ene)
2369 integer :: nfun,nodenum,i_orig,i,nf,nacc,it,nout,j,kkk,is,&
2370 Kwita,iretcode,noverlap,nstart_grow,ierr,iitt,&
2371 ii_grnum_d,ii_ennum_d,ii_hesnum_d,i_grnum_d,i_ennum_d,&
2372 i_hesnum_d,i_minimiz,irepet
2373 real(kind=8) :: etot,frac,eneglobal,RandOrPert,eold1,co,&
2376 ! if(.not.allocated(varsave)) allocate(varsave(maxvar,maxsave)) !(maxvar,maxsave)
2377 print *,'Master entered DO_MCM'
2385 !---------------------------------------------------------------------------
2386 ! Initialize counters.
2387 !---------------------------------------------------------------------------
2388 ! Total number of generated confs.
2390 ! Total number of moves. In general this won`t be equal to the number of
2391 ! attempted moves, because we may want to reject some "bad" confs just by
2394 ! Total number of temperature jumps.
2396 ! Total number of shift (nbond_move(1)), spike, crankshaft, three-bond,...
2398 allocate(nbond_move(nres)) !(maxres)
2404 ! Initialize total and accepted number of moves of various kind.
2409 ! Total number of energy evaluations.
2413 ! write (iout,*) 'RanFract=',RanFract
2416 !----------------------------------------------------------------------------
2417 ! Compute and print initial energies.
2418 !----------------------------------------------------------------------------
2420 write (iout,'(/80(1h*)/a)') 'Initial energies:'
2423 call etotal(energia)
2425 call enerprint(energia)
2426 ! Request energy computation from slave processors.
2427 call geom_to_var(nvar,varia)
2428 ! write (iout,*) 'The VARIA array'
2429 ! write (iout,'(8f10.4)') (rad2deg*varia(i),i=1,nvar)
2430 call minimize(etot,varia,iretcode,nfun)
2431 call var_to_geom(nvar,varia)
2433 write (iout,*) 'etot from MINIMIZE:',etot
2434 ! write (iout,*) 'Tha VARIA array'
2435 ! write (iout,'(8f10.4)') (rad2deg*varia(i),i=1,nvar)
2438 if (print_mc .gt. 0) write (iout,'(/80(1h*)/20x,a/80(1h*))') &
2439 'Enter Monte Carlo procedure.'
2440 if (print_mc .gt. 0) write (iout,'(i5,1pe14.5)' ) i_orig,etot
2446 call zapis(varia,etot)
2448 call var_to_geom(nvar,varia)
2450 call etotal(energia)
2451 if (print_mc.gt.0) write (iout,*) 'Initial energy:',etot
2453 nacc=0 ! total # of accepted confs of the current processor.
2454 nacc_tot=0 ! total # of accepted confs of all processors.
2456 !----------------------------------------------------------------------------
2458 !----------------------------------------------------------------------------
2461 LOOP1:do while (not_done)
2463 if (print_mc.gt.0) write (iout,'(80(1h*)/20x,a,i7)') &
2464 'Beginning iteration #',it
2465 ! Initialize local counter.
2466 ntrial=0 ! # of generated non-overlapping confs.
2467 noverlap=0 ! # of overlapping confs.
2469 LOOP2:do while (.not. accepted)
2471 LOOP3:do while (imtasks_n.lt.nodenum-1.and..not.finish)
2473 if(imtasks(i).eq.0) then
2478 ! Retrieve the angles of previously accepted conformation
2482 call var_to_geom(nvar,varia)
2483 ! Rebuild the chain.
2485 ! Heat up the system, if necessary.
2487 ! If temperature cannot be further increased, stop.
2493 ! write (iout,'(a)') 'Old variables:'
2494 ! write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
2495 ! Decide whether to generate a random conformation or perturb the old one
2496 RandOrPert=ran_number(0.0D0,1.0D0)
2497 if (RandOrPert.gt.RanFract) then
2498 if (print_mc.gt.0) &
2499 write (iout,'(a)') 'Perturbation-generated conformation.'
2500 call perturb(error,lprint,MoveType,nbond,1.0D0)
2501 ! print *,'after perturb',error,finish
2502 if (error) finish = .true.
2503 if (MoveType.lt.1 .or. MoveType.gt.MaxMoveType) then
2504 write (iout,'(/a,i7,a/)') 'Error - unknown MoveType=',&
2505 MoveType,' returned from PERTURB.'
2507 write (*,'(/a,i7,a/)') 'Error - unknown MoveType=',&
2508 MoveType,' returned from PERTURB.'
2514 nstart_grow=iran_num(3,nres)
2515 if (print_mc.gt.0) &
2516 write (iout,'(2a,i3)') 'Random-generated conformation',&
2517 ' - chain regrown from residue',nstart_grow
2518 call gen_rand_conf(nstart_grow,*30)
2520 call geom_to_var(nvar,varia)
2522 ! print *,'finish=',finish
2523 if (etot-elowest.gt.overlap_cut) then
2524 if (print_mc.gt.1) write (iout,'(a,1pe14.5)') &
2525 'Overlap detected in the current conf.; energy is',etot
2526 if(iprint.gt.1.or.etot.lt.1d19) print *,&
2527 'Overlap detected in the current conf.; energy is',etot
2531 if (noverlap.gt.maxoverlap) then
2532 write (iout,*) 'Too many overlapping confs.',&
2533 ' etot, elowest, overlap_cut', etot, elowest, overlap_cut
2536 else if (.not. finish) then
2537 ! Distribute tasks to processors
2538 ! print *,'Master sending order'
2539 call MPI_SEND(12, 1, MPI_INTEGER, is, tag,&
2541 ! write (iout,*) '12: tag=',tag
2542 ! print *,'Master sent order to processor',is
2543 call MPI_SEND(it, 1, MPI_INTEGER, is, tag,&
2545 ! write (iout,*) 'it: tag=',tag
2546 call MPI_SEND(eold, 1, MPI_DOUBLE_PRECISION, is, tag,&
2548 ! write (iout,*) 'eold: tag=',tag
2549 call MPI_SEND(varia(1), nvar, MPI_DOUBLE_PRECISION, &
2552 ! write (iout,*) 'varia: tag=',tag
2553 call MPI_SEND(varold(1), nvar, MPI_DOUBLE_PRECISION, &
2556 ! write (iout,*) 'varold: tag=',tag
2563 imtasks_n=imtasks_n+1
2569 LOOP_RECV:do while(.not.flag)
2571 call MPI_IPROBE(is,tag,CG_COMM,flag,status,ierr)
2573 call MPI_RECV(iitt, 1, MPI_INTEGER, is, tag,&
2574 CG_COMM, status, ierr)
2575 call MPI_RECV(eold1, 1, MPI_DOUBLE_PRECISION, is, tag,&
2576 CG_COMM, status, ierr)
2577 call MPI_RECV(etot, 1, MPI_DOUBLE_PRECISION, is, tag,&
2578 CG_COMM, status, ierr)
2579 call MPI_RECV(varia(1), nvar, MPI_DOUBLE_PRECISION,is,tag,&
2580 CG_COMM, status, ierr)
2581 call MPI_RECV(varold1(1), nvar, MPI_DOUBLE_PRECISION, is, &
2582 tag, CG_COMM, status, ierr)
2583 call MPI_RECV(ii_grnum_d, 1, MPI_INTEGER, is, tag,&
2584 CG_COMM, status, ierr)
2585 call MPI_RECV(ii_ennum_d, 1, MPI_INTEGER, is, tag,&
2586 CG_COMM, status, ierr)
2587 call MPI_RECV(ii_hesnum_d, 1, MPI_INTEGER, is, tag,&
2588 CG_COMM, status, ierr)
2589 i_grnum_d=i_grnum_d+ii_grnum_d
2590 i_ennum_d=i_ennum_d+ii_ennum_d
2591 neneval = neneval+ii_ennum_d
2592 i_hesnum_d=i_hesnum_d+ii_hesnum_d
2593 i_minimiz=i_minimiz+1
2595 imtasks_n=imtasks_n-1
2601 if(print_mc.gt.0) write (iout,'(a,i6,a,i6,a,i6,a,1pe16.6)') &
2602 'From Worker #',is,' iitt',iitt,&
2603 ' Conformation:',ngen,' energy:',etot
2604 !--------------------------------------------------------------------------
2605 !... Do Metropolis test
2606 !--------------------------------------------------------------------------
2607 call metropolis(nvar,varia,varold1,etot,eold1,accepted,&
2609 if(iitt.ne.it.and..not.similar) then
2610 call metropolis(nvar,varia,varold,etot,eold,accepted,&
2614 if(etot.lt.eneglobal)eneglobal=etot
2615 ! if(mod(it,100).eq.0)
2616 write(iout,*)'CHUJOJEB ',neneval,eneglobal
2618 ! Write the accepted conformation.
2621 call var_to_geom(nvar,varia)
2624 call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,kkk),&
2625 nsup,przes,obr,non_conv)
2627 call contact(.false.,ncont,icont,co)
2628 frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
2629 write (iout,'(a,f8.3,a,f8.3,a,f8.3)') &
2630 'RMS deviation from the reference structure:',rms,&
2631 ' % of native contacts:',frac*100,' contact order:',co
2633 if (print_mc.gt.0) &
2634 write (iout,*) 'Writing new conformation',nout
2635 if (print_stat) then
2636 call var_to_geom(nvar,varia)
2637 #if defined(AIX) || defined(PGI)
2638 open (istat,file=statname,position='append')
2640 open (istat,file=statname,access='append')
2643 write (istat,'(i5,16(1pe14.5))') nout,&
2644 (energia(print_order(i)),i=1,nprint_ene),&
2647 write (istat,'(i5,16(1pe14.5))') nout,&
2648 (energia(print_order(i)),i=1,nprint_ene),etot
2652 ! Print internal coordinates.
2653 if (print_int) call briefout(nout,etot)
2656 if (elowest.gt.etot) elowest=etot
2657 moves_acc(MoveType)=moves_acc(MoveType)+1
2658 if (MoveType.eq.1) then
2659 nbond_acc(nbond)=nbond_acc(nbond)+1
2661 ! Check against conformation repetitions.
2662 irepet=conf_comp(varia,etot)
2663 if (nrepm.gt.maxrepm) then
2664 if (print_mc.gt.0) &
2665 write (iout,'(a)') 'Too many conformation repetitions.'
2668 ! Store the accepted conf. and its energy.
2673 if (irepet.eq.0) call zapis(varia,etot)
2674 ! Lower the temperature, if necessary.
2680 if(finish.and.imtasks_n.eq.0)exit LOOP2
2681 enddo LOOP2 ! accepted
2682 ! Check for time limit.
2683 not_done = (it.lt.max_mcm_it) .and. (nacc_tot.lt.maxacc)
2684 if(.not.not_done .or. finish) then
2685 if(imtasks_n.gt.0) then
2692 enddo LOOP1 ! not_done
2694 if (print_mc.gt.0) then
2695 write (iout,'(/80(1h*)/20x,a)') 'Summary run statistics:'
2696 call statprint(nacc,nfun,iretcode,etot,elowest)
2697 write (iout,'(a)') &
2698 'Statistics of multiple-bond motions. Total motions:'
2699 write (iout,'(16i5)') (nbond_move(i),i=1,Nbm)
2700 write (iout,'(a)') 'Accepted motions:'
2701 write (iout,'(16i5)') (nbond_acc(i),i=1,Nbm)
2703 write (iout,'(/80(1h*)/20x,a/80(1h*)/)') 'All iterations done.'
2711 call MPI_SEND(999, 1, MPI_INTEGER, is, tag,&
2715 end subroutine do_mcm
2716 !-----------------------------------------------------------------------------
2717 subroutine execute_slave(nodeinfo,iprint)
2720 use minimm, only:minimize
2722 ! implicit real*8 (a-h,o-z)
2723 ! include 'DIMENSIONS'
2725 ! include 'COMMON.TIME1'
2726 ! include 'COMMON.IOUNITS'
2727 !rc include 'COMMON.DEFORM'
2728 !rc include 'COMMON.DEFORM1'
2729 !rc include 'COMMON.DEFORM2'
2730 ! include 'COMMON.LOCAL'
2731 ! include 'COMMON.VAR'
2732 ! include 'COMMON.INFO'
2733 ! include 'COMMON.MINIM'
2734 character(len=10) :: nodeinfo
2735 real(kind=8),dimension(6*nres) :: x,x1 !(maxvar) (maxvar=6*maxres)
2736 integer :: nfun,iprint,i_switch,ierr,i_grnum_d,i_ennum_d,&
2737 i_hesnum_d,iitt,iretcode,iminrep
2738 real(kind=8) :: ener,energyx
2740 nodeinfo='chujwdupe'
2741 ! print *,'Processor:',MyID,' Entering execute_slave'
2743 ! call MPI_SEND(nodeinfo, 10, MPI_CHARACTER, 0, tag,
2746 1001 call MPI_RECV(i_switch, 1, MPI_INTEGER, 0, tag,&
2747 CG_COMM, status, ierr)
2748 ! write(iout,*)'12: tag=',tag
2749 if(iprint.ge.2)print *, MyID,' recv order ',i_switch
2750 if (i_switch.eq.12) then
2754 call MPI_RECV(iitt, 1, MPI_INTEGER, 0, tag,&
2755 CG_COMM, status, ierr)
2756 ! write(iout,*)'12: tag=',tag
2757 call MPI_RECV(ener, 1, MPI_DOUBLE_PRECISION, 0, tag,&
2758 CG_COMM, status, ierr)
2759 ! write(iout,*)'ener: tag=',tag
2760 call MPI_RECV(x(1), nvar, MPI_DOUBLE_PRECISION, 0, tag,&
2761 CG_COMM, status, ierr)
2762 ! write(iout,*)'x: tag=',tag
2763 call MPI_RECV(x1(1), nvar, MPI_DOUBLE_PRECISION, 0, tag,&
2764 CG_COMM, status, ierr)
2765 ! write(iout,*)'x1: tag=',tag
2771 ! print *,'calling minimize'
2772 call minimize(energyx,x,iretcode,nfun)
2774 write(iout,100)'minimized energy = ',energyx,&
2775 ' # funeval:',nfun,' iret ',iretcode
2776 write(*,100)'minimized energy = ',energyx,&
2777 ' # funeval:',nfun,' iret ',iretcode
2778 100 format(a20,f10.5,a12,i5,a6,i2)
2779 if(iretcode.eq.10) then
2782 write(iout,*)' ... not converged - trying again ',iminrep
2783 call minimize(energyx,x,iretcode,nfun)
2785 write(iout,*)'minimized energy = ',energyx,&
2786 ' # funeval:',nfun,' iret ',iretcode
2787 if(iretcode.ne.10)go to 812
2789 if(iretcode.eq.10) then
2791 write(iout,*)' ... not converged again - giving up'
2796 ! print *,'Sending results'
2797 call MPI_SEND(iitt, 1, MPI_INTEGER, 0, tag,&
2799 call MPI_SEND(ener, 1, MPI_DOUBLE_PRECISION, 0, tag,&
2801 call MPI_SEND(energyx, 1, MPI_DOUBLE_PRECISION, 0, tag,&
2803 call MPI_SEND(x(1), nvar, MPI_DOUBLE_PRECISION, 0, tag,&
2805 call MPI_SEND(x1(1), nvar, MPI_DOUBLE_PRECISION, 0, tag,&
2807 call MPI_SEND(i_grnum_d, 1, MPI_INTEGER, 0, tag,&
2809 call MPI_SEND(nfun, 1, MPI_INTEGER, 0, tag,&
2811 call MPI_SEND(i_hesnum_d, 1, MPI_INTEGER, 0, tag,&
2813 ! print *,'End sending'
2818 end subroutine execute_slave
2820 !-----------------------------------------------------------------------------
2821 subroutine heat(over)
2823 ! implicit real*8 (a-h,o-z)
2824 ! include 'DIMENSIONS'
2825 ! include 'COMMON.MCM'
2826 ! include 'COMMON.IOUNITS'
2828 ! Check if there`s a need to increase temperature.
2829 if (ntrial.gt.maxtrial) then
2830 if (NstepH.gt.0) then
2831 if (dabs(Tcur-TMax).lt.1.0D-7) then
2832 if (print_mc.gt.0) &
2833 write (iout,'(/80(1h*)/a,f8.3,a/80(1h*))') &
2834 'Upper limit of temperature reached. Terminating.'
2839 if (Tcur.gt.Tmax) Tcur=Tmax
2840 betbol=1.0D0/(Rbol*Tcur)
2841 if (print_mc.gt.0) &
2842 write (iout,'(/80(1h*)/a,f8.3,a,f10.5/80(1h*))') &
2843 'System heated up to ',Tcur,' K; BetBol:',betbol
2849 if (print_mc.gt.0) &
2850 write (iout,'(a)') &
2851 'Maximum number of trials in a single MCM iteration exceeded.'
2860 !-----------------------------------------------------------------------------
2863 ! implicit real*8 (a-h,o-z)
2864 ! include 'DIMENSIONS'
2865 ! include 'COMMON.MCM'
2866 ! include 'COMMON.IOUNITS'
2867 if (nstepC.gt.0 .and. dabs(Tcur-Tmin).gt.1.0D-7) then
2869 if (Tcur.lt.Tmin) Tcur=Tmin
2870 betbol=1.0D0/(Rbol*Tcur)
2871 if (print_mc.gt.0) &
2872 write (iout,'(/80(1h*)/a,f8.3,a,f10.5/80(1h*))') &
2873 'System cooled down up to ',Tcur,' K; BetBol:',betbol
2877 !-----------------------------------------------------------------------------
2878 subroutine perturb(error,lprint,MoveType,nbond,max_phi)
2881 use energy, only:nnt,nct,itype
2882 use md_calc, only:bond_move
2883 ! implicit real*8 (a-h,o-z)
2884 ! include 'DIMENSIONS'
2885 integer,parameter :: MMaxSideMove=100
2886 ! include 'COMMON.MCM'
2887 ! include 'COMMON.CHAIN'
2888 ! include 'COMMON.INTERACT'
2889 ! include 'COMMON.VAR'
2890 ! include 'COMMON.GEO'
2891 ! include 'COMMON.NAMES'
2892 ! include 'COMMON.IOUNITS'
2893 !rc include 'COMMON.DEFORM1'
2894 logical :: error,lprint,fail
2895 integer :: MoveType,nbond,end_select,ind_side(MMaxSideMove)
2896 real(kind=8) :: max_phi
2897 real(kind=8) :: psi!,gen_psi
2898 !el external iran_num
2899 !el integer iran_num
2902 !!! local variables - el
2903 integer :: itrial,iiwin,iwindow,isctry,i,icount,j,nstart,&
2904 nside_move,inds,indx,ii,iti
2905 real(kind=8) :: bond_prob,theta_new
2910 ! Perturb the conformation according to a randomly selected move.
2911 call SelectMove(MoveType)
2912 ! write (iout,*) 'MoveType=',MoveType
2914 goto (100,200,300,400,500) MoveType
2915 !------------------------------------------------------------------------------
2916 ! Backbone N-bond move.
2917 ! Select the number of bonds (length of the segment to perturb).
2919 if (itrial.gt.1000) then
2920 write (iout,'(a)') 'Too many attempts at multiple-bond move.'
2924 bond_prob=ran_number(0.0D0,sumpro_bond(nbm))
2925 ! print *,'sumpro_bond(nbm)=',sumpro_bond(nbm),
2926 ! & ' Bond_prob=',Bond_Prob
2928 ! print *,i,Bond_Prob,sumpro_bond(i),sumpro_bond(i+1)
2929 if (bond_prob.ge.sumpro_bond(i) .and. &
2930 bond_prob.le.sumpro_bond(i+1)) then
2935 write (iout,'(2a)') 'In PERTURB: Error - number of bonds',&
2936 ' to move out of range.'
2940 if (nwindow.gt.0) then
2941 ! Select the first residue to perturb
2942 iwindow=iran_num(1,nwindow)
2943 print *,'iwindow=',iwindow
2945 do while (winlen(iwindow).lt.nbond)
2946 iwindow=iran_num(1,nwindow)
2948 if (iiwin.gt.1000) then
2949 write (iout,'(a)') 'Cannot select moveable residues.'
2954 nstart=iran_num(winstart(iwindow),winend(iwindow))
2956 nstart = iran_num(koniecl+2,nres-nbond-koniecl)
2957 !d print *,'nres=',nres,' nbond=',nbond,' koniecl=',koniecl,
2958 !d & ' nstart=',nstart
2961 if (psi.eq.0.0) then
2965 if (print_mc.gt.1) write (iout,'(a,i4,a,i4,a,f8.3)') &
2966 'PERTURB: nbond=',nbond,' nstart=',nstart,' psi=',psi*rad2deg
2967 !d print *,'nstart=',nstart
2968 call bond_move(nbond,nstart,psi,.false.,error)
2970 write (iout,'(2a)') &
2971 'Could not define reference system in bond_move, ',&
2972 'choosing ahother segment.'
2976 nbond_move(nbond)=nbond_move(nbond)+1
2980 !------------------------------------------------------------------------------
2981 ! Backbone endmove. Perturb a SINGLE angle of a residue close to the end of
2985 ! end_select=iran_num(1,2*koniecl)
2986 ! if (end_select.gt.koniecl) then
2987 ! end_select=nphi-(end_select-koniecl)
2989 ! end_select=koniecl+3
2991 ! if (nwindow.gt.0) then
2992 ! iwin=iran_num(1,nwindow)
2993 ! i1=max0(4,winstart(iwin))
2994 ! i2=min0(winend(imin)+2,nres)
2995 ! end_select=iran_num(i1,i2)
2997 ! iselect = iran_num(1,nmov_var)
3000 ! if (isearch_tab(i).eq.1) jj = jj+1
3001 ! if (jj.eq.iselect) then
3007 end_select = iran_num(4,nres)
3008 psi=max_phi*gen_psi()
3009 if (psi.eq.0.0D0) then
3013 phi(end_select)=pinorm(phi(end_select)+psi)
3014 if (print_mc.gt.1) write (iout,'(a,i4,a,f8.3,a,f8.3)') &
3015 'End angle',end_select,' moved by ',psi*rad2deg,' new angle:',&
3016 phi(end_select)*rad2deg
3017 ! if (end_select.gt.3)
3018 ! & theta(end_select-1)=gen_theta(itype(end_select-2),
3019 ! & phi(end_select-1),phi(end_select))
3020 ! if (end_select.lt.nres)
3021 ! & theta(end_select)=gen_theta(itype(end_select-1),
3022 ! & phi(end_select),phi(end_select+1))
3023 !d print *,'nres=',nres,' end_select=',end_select
3024 !d print *,'theta',end_select-1,theta(end_select-1)
3025 !d print *,'theta',end_select,theta(end_select)
3030 !------------------------------------------------------------------------------
3032 ! Select the number of SCs to perturb.
3034 301 nside_move=iran_num(1,MaxSideMove)
3035 ! print *,'nside_move=',nside_move,' MaxSideMove',MaxSideMove
3036 ! Select the indices.
3039 111 inds=iran_num(nnt,nct)
3041 if (icount.gt.1000) then
3042 write (iout,'(a)')'Error - cannot select side chains to move.'
3046 if (itype(inds,1).eq.10) goto 111
3048 if (inds.eq.ind_side(j)) goto 111
3051 if (inds.lt.ind_side(j)) then
3057 112 do j=i,indx+1,-1
3058 ind_side(j)=ind_side(j-1)
3060 113 ind_side(indx)=inds
3062 ! Carry out perturbation.
3066 call gen_side(iti,theta(ii+1),alph(ii),omeg(ii),fail)
3069 if (isctry.gt.1000) then
3070 write (iout,'(a)') 'Too many errors in SC generation.'
3076 if (print_mc.gt.1) write (iout,'(2a,i4,a,2f8.3)') &
3077 'Side chain ',restyp(iti,1),ii,' moved to ',&
3078 alph(ii)*rad2deg,omeg(ii)*rad2deg
3083 !------------------------------------------------------------------------------
3085 400 end_select=iran_num(3,nres)
3086 theta_new=gen_theta(itype(end_select,1),phi(end_select),&
3088 if (print_mc.gt.1) write (iout,'(a,i3,a,f8.3,a,f8.3)') &
3089 'Theta ',end_select,' moved from',theta(end_select)*rad2deg,&
3090 ' to ',theta_new*rad2deg
3091 theta(end_select)=theta_new
3095 !------------------------------------------------------------------------------
3096 ! Error returned from SelectMove.
3099 end subroutine perturb
3100 !-----------------------------------------------------------------------------
3101 subroutine SelectMove(MoveType)
3103 ! implicit real*8 (a-h,o-z)
3104 ! include 'DIMENSIONS'
3105 ! include 'COMMON.MCM'
3106 ! include 'COMMON.IOUNITS'
3108 !!! local variables - el
3109 integer :: i,MoveType
3110 real(kind=8) :: what_move
3112 what_move=ran_number(0.0D0,sumpro_type(MaxMoveType))
3114 if (what_move.ge.sumpro_type(i-1).and. &
3115 what_move.lt.sumpro_type(i)) then
3120 write (iout,'(a)') &
3121 'Fatal error in SelectMoveType: cannot select move.'
3122 MoveType=MaxMoveType+1
3124 end subroutine SelectMove
3125 !-----------------------------------------------------------------------------
3126 real(kind=8) function gen_psi()
3128 use geometry_data, only: angmin,pi
3131 real(kind=8) :: x !,ran_number
3132 ! include 'COMMON.IOUNITS'
3133 ! include 'COMMON.GEO'
3136 x=ran_number(-pi,pi)
3137 if (dabs(x).gt.angmin) then
3142 write (iout,'(a)')'From Gen_Psi: Cannot generate angle increment.'
3145 end function gen_psi
3146 !-----------------------------------------------------------------------------
3147 subroutine metropolis(n,xcur,xold,ecur,eold,accepted,similar,enelower)
3149 ! implicit real*8 (a-h,o-z)
3150 ! include 'DIMENSIONS'
3151 ! include 'COMMON.MCM'
3152 ! include 'COMMON.IOUNITS'
3153 ! include 'COMMON.VAR'
3154 ! include 'COMMON.GEO'
3155 !rc include 'COMMON.DEFORM'
3157 real(kind=8) :: ecur,eold,xx,bol !,ran_number
3158 real(kind=8),dimension(n) :: xcur,xold
3159 real(kind=8) :: ecut1 ,ecut2 ,tola
3160 logical :: accepted,similar,not_done,enelower
3163 !!! local variables - el
3164 real(kind=8) :: xxh,difene,reldife
3166 data ecut1 /-1.0D-5/,ecut2 /5.0D-3/,tola/5.0D0/
3170 ! Set lprn=.true. for debugging.
3173 !el write(iout,*)'enedif',enedif,' ecut1',ecut1,' ecut2',ecut2
3174 write(iout,*)' ecut1',ecut1,' ecut2',ecut2
3178 ! Check if the conformation is similar.
3180 reldife=difene/dmax1(dabs(eold),dabs(ecur),1.0D0)
3182 write (iout,*) 'Metropolis'
3183 write(iout,*)'ecur,eold,difene,reldife',ecur,eold,difene,reldife
3185 ! If energy went down remarkably, we accept the new conformation
3187 !jp if (reldife.lt.ecut1) then
3188 if (difene.lt.ecut1) then
3191 if (lprn) write (iout,'(a)') &
3192 'Conformation accepted, because energy has lowered remarkably.'
3193 ! elseif (reldife.lt.ecut2 .and. dif_ang(nphi,xcur,xold).lt.tola)
3194 !jp elseif (reldife.lt.ecut2)
3195 elseif (difene.lt.ecut2) &
3197 ! Reject the conf. if energy has changed insignificantly and there is not
3198 ! much change in conformation.
3200 write (iout,'(2a)') 'Conformation rejected, because it is',&
3201 ' similar to the preceding one.'
3205 ! Else carry out Metropolis test.
3207 xx=ran_number(0.0D0,1.0D0)
3210 write (iout,*) 'betbol=',betbol,' difene=',difene,' xxh=',xxh
3211 if (xxh.gt.50.0D0) then
3216 if (lprn) write (iout,*) 'bol=',bol,' xx=',xx
3219 if (lprn) write (iout,'(a)') &
3220 'Conformation accepted, because it passed Metropolis test.'
3223 if (lprn) write (iout,'(a)') &
3224 'Conformation rejected, because it did not pass Metropolis test.'
3233 end subroutine metropolis
3234 !-----------------------------------------------------------------------------
3235 integer function conf_comp(x,ene)
3237 use geometry_data, only: nphi
3238 ! implicit real*8 (a-h,o-z)
3239 ! include 'DIMENSIONS'
3240 ! include 'COMMON.MCM'
3241 ! include 'COMMON.VAR'
3242 ! include 'COMMON.IOUNITS'
3243 ! include 'COMMON.GEO'
3244 real(kind=8) :: etol, angtol
3245 real(kind=8),dimension(6*nres) :: x !(maxvar) (maxvar=6*maxres)
3246 real(kind=8) :: difa !dif_ang,
3248 !!! local variables - el
3252 data etol /0.1D0/, angtol /20.0D0/
3254 ! write (iout,*) 'ii=',ii,'ene=',ene,esave(ii),dabs(ene-esave(ii))
3255 if (dabs(ene-esave(ii)).lt.etol) then
3256 difa=dif_ang(nphi,x,varsave(1,ii))
3258 ! write(iout,'(i3,3f8.3)')i,rad2deg*x(i),
3259 ! & rad2deg*varsave(i,ii)
3261 ! write(iout,*) 'ii=',ii,' difa=',difa,' angtol=',angtol
3262 if (difa.le.angtol) then
3263 if (print_mc.gt.0) then
3264 write (iout,'(a,i5,2(a,1pe15.4))') &
3265 'Current conformation matches #',ii,&
3266 ' in the store array ene=',ene,' esave=',esave(ii)
3267 ! write (*,'(a,i5,a)') 'Current conformation matches #',ii,
3268 ! & ' in the store array.'
3269 endif ! print_mc.gt.0
3270 if (print_mc.gt.1) then
3272 write(iout,'(i3,3f8.3)')i,rad2deg*x(i),&
3273 rad2deg*varsave(i,ii)
3275 endif ! print_mc.gt.1
3284 end function conf_comp
3285 !-----------------------------------------------------------------------------
3286 real(kind=8) function dif_ang(n,x,y)
3288 use geometry_data, only: dwapi
3291 real(kind=8),dimension(n) :: x,y
3292 real(kind=8) :: w,wa,dif,difa
3293 !el real(kind=8) :: pinorm
3294 ! include 'COMMON.GEO'
3298 dif=dabs(pinorm(y(i)-x(i)))
3299 if (dabs(dif-dwapi).lt.dif) dif=dabs(dif-dwapi)
3300 w=1.0D0-(2.0D0*(i-1)/(n-1)-1.0D0)**2+1.0D0/n
3304 dif_ang=rad2deg*dsqrt(difa/wa)
3306 end function dif_ang
3307 !-----------------------------------------------------------------------------
3308 subroutine add2cache(n1,n2,ncache,nvar,SourceID,CachSrc,ecur,xcur,ecache,xcache)
3311 ! include 'COMMON.GEO'
3312 ! include 'COMMON.IOUNITS'
3313 integer :: n1,n2,ncache,nvar,SourceID,CachSrc(n2)
3315 real(kind=8) :: ecur,xcur(nvar),ecache(n2),xcache(n1,n2)
3316 !d write (iout,*) 'Enter ADD2CACHE ncache=',ncache ,' ecur',ecur
3317 !d write (iout,'(10f8.3)') (rad2deg*xcur(i),i=1,nvar)
3318 !d write (iout,*) 'Old CACHE array:'
3320 !d write (iout,*) 'i=',i,' ecache=',ecache(i),' CachSrc',CachSrc(i)
3321 !d write (iout,'(10f8.3)') (rad2deg*xcache(j,i),j=1,nvar)
3325 do while (i.gt.0 .and. ecur.lt.ecache(i))
3329 !d write (iout,*) 'i=',i,' ncache=',ncache
3330 if (ncache.eq.n2) then
3331 write (iout,*) 'Cache dimension exceeded',ncache,n2
3332 write (iout,*) 'Highest-energy conformation will be removed.'
3336 ecache(ii+1)=ecache(ii)
3337 CachSrc(ii+1)=CachSrc(ii)
3339 xcache(j,ii+1)=xcache(j,ii)
3348 !d write (iout,*) 'New CACHE array:'
3350 !d write (iout,*) 'i=',i,' ecache=',ecache(i),' CachSrc',CachSrc(i)
3351 !d write (iout,'(10f8.3)') (rad2deg*xcache(j,i),j=1,nvar)
3354 end subroutine add2cache
3355 !-----------------------------------------------------------------------------
3356 subroutine rm_from_cache(i,n1,n2,ncache,nvar,CachSrc,ecache,xcache)
3359 ! include 'COMMON.GEO'
3360 ! include 'COMMON.IOUNITS'
3361 integer :: n1,n2,ncache,nvar,CachSrc(n2)
3363 real(kind=8) :: ecache(n2),xcache(n1,n2)
3365 !d write (iout,*) 'Enter RM_FROM_CACHE'
3366 !d write (iout,*) 'Old CACHE array:'
3368 !d write (iout,*)'i=',ii,' ecache=',ecache(ii),' CachSrc',CachSrc(ii)
3369 !d write (iout,'(10f8.3)') (rad2deg*xcache(j,ii),j=1,nvar)
3373 ecache(ii-1)=ecache(ii)
3374 CachSrc(ii-1)=CachSrc(ii)
3376 xcache(j,ii-1)=xcache(j,ii)
3380 !d write (iout,*) 'New CACHE array:'
3382 !d write (iout,*) 'i=',i,' ecache=',ecache(i),' CachSrc',CachSrc(i)
3383 !d write (iout,'(10f8.3)') (rad2deg*xcache(j,i),j=1,nvar)
3386 end subroutine rm_from_cache
3387 !-----------------------------------------------------------------------------
3389 !-----------------------------------------------------------------------------
3390 subroutine statprint(it,nfun,iretcode,etot,elowest)
3392 use control_data, only: MaxMoveType,minim
3393 use control, only: tcpu
3395 ! implicit real*8 (a-h,o-z)
3396 ! include 'DIMENSIONS'
3397 ! include 'COMMON.IOUNITS'
3398 ! include 'COMMON.CONTROL'
3399 ! include 'COMMON.MCM'
3401 integer :: it,nfun,iretcode,i
3402 real(kind=8) :: etot,elowest,fr_mov_i
3406 '(80(1h*)/a,i5,a,1pe14.5,a,1pe14.5/a,i3,a,i10,a,i5,a,i5)') &
3407 'Finished iteration #',it,' energy is',etot,&
3408 ' lowest energy:',elowest,&
3409 'SUMSL return code:',iretcode,&
3410 ' # of energy evaluations:',neneval,&
3411 '# of temperature jumps:',ntherm,&
3412 ' # of minima repetitions:',nrepm
3414 write (iout,'(80(1h*)/a,i8,a,1pe14.5,a,1pe14.5)') &
3415 'Finished iteration #',it,' energy is',etot,&
3416 ' lowest energy:',elowest
3418 write (iout,'(/4a)') &
3419 'Kind of move ',' total',' accepted',&
3421 write (iout,'(58(1h-))')
3423 if (moves(i).eq.0) then
3426 fr_mov_i=dfloat(moves_acc(i))/dfloat(moves(i))
3428 write(iout,'(a,2i15,f10.5)')MovTypID(i),moves(i),moves_acc(i),&
3431 write (iout,'(a,2i15,f10.5)') 'total ',nmove,nacc_tot,&
3432 dfloat(nacc_tot)/dfloat(nmove)
3433 write (iout,'(58(1h-))')
3434 write (iout,'(a,1pe12.4)') 'Elapsed time:',tcpu()
3436 end subroutine statprint
3437 !-----------------------------------------------------------------------------
3438 subroutine zapis(varia,etot)
3440 use geometry_data, only: nres,rad2deg,nvar
3443 ! implicit real*8 (a-h,o-z)
3444 ! include 'DIMENSIONS'
3446 use MPI_data !include 'COMMON.INFO'
3449 ! include 'COMMON.GEO'
3450 ! include 'COMMON.VAR'
3451 ! include 'COMMON.MCM'
3452 ! include 'COMMON.IOUNITS'
3453 integer,dimension(nsave) :: itemp
3454 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
3457 integer :: j,i,maxvar
3458 real(kind=8) :: etot
3460 !el allocate(esave(nsave)) !(maxsave)
3465 write (iout,'(a,i5,a,i5)') 'Enter ZAPIS NSave=',Nsave,&
3467 write (iout,'(a)') 'Current energy and conformation:'
3468 write (iout,'(1pe14.5)') etot
3469 write (iout,'(10f8.3)') (rad2deg*varia(i),i=1,nvar)
3471 ! Shift the contents of the esave and varsave arrays if filled up.
3472 call add2cache(6*nres,nsave,nsave,nvar,MyID,itemp,&
3473 etot,varia,esave,varsave)
3475 write (iout,'(a)') 'Energies and the VarSave array.'
3477 write (iout,'(i5,1pe14.5)') i,esave(i)
3478 write (iout,'(10f8.3)') (rad2deg*varsave(j,i),j=1,nvar)
3482 end subroutine zapis
3483 !-----------------------------------------------------------------------------
3484 !-----------------------------------------------------------------------------
3485 subroutine alloc_MCM_arrays
3487 use energy_data, only: max_ene
3491 allocate(entropy(-max_ene-4:max_ene)) !(-max_ene-4:max_ene)
3492 allocate(nhist(-max_ene:max_ene)) !(-max_ene:max_ene)
3493 allocate(nminima(maxsave)) !(maxsave)
3495 allocate(xpool(6*nres,max_pool)) !(maxvar,max_pool)(maxvar=6*maxres)
3496 allocate(epool(max_pool)) !(max_pool)
3499 if(.not.allocated(nsave_part)) allocate(nsave_part(nctasks)) !(max_cg_procs)
3502 ! real(kind=8),dimension(:),allocatable :: sumpro_type !(0:MaxMoveType)
3503 allocate(sumpro_bond(0:nres)) !(0:maxres)
3504 allocate(nbond_move(nres),nbond_acc(nres)) !(maxres)
3505 allocate(moves(-1:MaxMoveType+1),moves_acc(-1:MaxMoveType+1)) !(-1:MaxMoveType+1)
3506 ! common /accept_stats/
3507 ! allocate(nacc_part !(0:MaxProcs) !el nie uzywane???
3508 ! common /windows/ in io: mcmread
3509 ! allocate(winstart,winend,winlen !(maxres)
3511 !el allocate(MovTypID(-1:MaxMoveType+1)) !(-1:MaxMoveType+1)
3514 allocate(varsave(nres*6,maxsave)) !(maxvar,maxsave)(maxvar=6*maxres)
3515 allocate(esave(maxsave)) !(maxsave)
3516 allocate(Origin(maxsave)) !(maxsave)
3519 end subroutine alloc_MCM_arrays
3520 !-----------------------------------------------------------------------------
3521 !-----------------------------------------------------------------------------