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)
171 ! write(iout,*)'*ene=',energyx
178 if (iprint.gt.3) then
179 write (iout,*) 'Compare_ss, i=',i
180 write (iout,*) 'New structure Energy:',energyx
181 write (iout,'(10f8.3)') (rad2deg*x(k),k=1,nvar)
182 write (iout,*) 'Template structure Energy:',enetbss(i)
183 write (iout,'(10f8.3)') (rad2deg*x1(k),k=1,nvar)
187 call var_to_geom(nvar,x1)
189 !d write(iout,*)'C and CREF'
190 !d write(iout,'(i5,3f10.5,5x,3f10.5)')(k,(c(j,k),j=1,3),
191 !d & (cref(j,k),j=1,3),k=1,nres)
192 call fitsq(roznica,c(1,1),cref(1,1,1),nres,przes,obrot,non_conv)
194 print *,'Problems in FITSQ!!!'
196 print '(10f8.3)',(x(k),k=1,nvar)
198 print '(10f8.3)',(x1(k),k=1,nvar)
200 print '(i5,3f10.5,5x,3f10.5)',(k,(c(j,k),j=1,3),&
201 (cref(j,k,1),j=1,3),k=1,nres)
203 roznica=dsqrt(dabs(roznica))
205 if (roznica.lt.rms_d) iresult = 0
208 !el call cmprs(x,x1,roznica,energyx,energyy,iresult)
210 if (iprint.gt.1) write(iout,'(i5,f10.6,$)') i,roznica
211 ! print '(i5,f8.3)',i,roznica
212 if(iresult.eq.0) then
215 if (iprint.gt.1) write(iout,*)
216 if(energyx.ge.enetbss(i)) then
218 write(iout,*)'s*>> structure rejected - same as nr ',i, &
225 if(energyx.lt.enetbss(i).and.enex_jp.lt.enetbss(i))then
230 if (iprint.gt.1) write(iout,*)
234 write(iout,'(a,i3,$)')'s*>> structure accepted1 - repl nr ',&
238 write(iout,'(a,i3)') &
239 's*>> structure accepted1 - would repl nr ',list(1)
242 if (.not. modif) goto 1106
249 if (iprint.gt.1) write(iout,'(i3,$)')list(j)
250 do kk=list(j)+1,nlist
251 enetbss(kk-1)=enetbss(kk)
253 coordss(i,kk-1)=coordss(i,kk)
257 if (iprint.gt.1) write(iout,*)
260 if(n_thr.lt.num_thread_save) then
264 write(iout,*)'s*>> structure accepted - add with nr ',n_thr+1
267 write(iout,*)'s*>> structure accepted - would add with nr ',&
272 enetbss(n_thr)=energyx
274 coordss(i,n_thr)=x(i)
279 write(iout,*)'s*>> structure rejected - too high energy'
286 write(iout,*)'s*>> structure accepted - repl nr ',j
289 write(iout,*)'s*>> structure accepted - would repl nr ',j
300 end subroutine compare_s1
301 !-----------------------------------------------------------------------------
303 !-----------------------------------------------------------------------------
308 use MPI_data, only:WhatsUp,MyID
309 use compare_data, only: ener
310 use control_data, only: minim,refstr
312 use regularize_, only:fitsq
313 use control, only: tcpu,ovrtim
315 use minimm, only:minimize
316 ! Does modified entropic sampling in the space of minima.
317 ! implicit real*8 (a-h,o-z)
318 ! include 'DIMENSIONS'
319 ! include 'COMMON.IOUNITS'
321 use MPI_data !include 'COMMON.INFO'
323 ! include 'COMMON.GEO'
324 ! include 'COMMON.CHAIN'
325 ! include 'COMMON.MCM'
326 ! include 'COMMON.MCE'
327 ! include 'COMMON.CONTACTS'
328 ! include 'COMMON.CONTROL'
329 ! include 'COMMON.VAR'
330 ! include 'COMMON.INTERACT'
331 ! include 'COMMON.THREAD'
332 ! include 'COMMON.NAMES'
333 logical :: accepted,not_done,over,error,lprint !,ovrtim
334 integer :: MoveType,nbond
335 ! integer :: conf_comp
336 real(kind=8) :: RandOrPert
337 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
338 real(kind=8) :: elowest,ehighest,eold
339 real(kind=8) :: przes(3),obr(3,3)
340 real(kind=8),dimension(6*nres) :: varold !(maxvar) (maxvar=6*maxres)
342 real(kind=8),dimension(0:n_ene) :: energia,energia_ave
344 !!! local variables -el
345 integer :: i,ii,kkk,it,j,nacc,nfun,ijunk,indmin,indmax,&
346 ISWEEP,Kwita,iretcode,indeold,iene,noverlap,&
347 irep,nstart_grow,inde
348 real(kind=8) :: facee,conste,ejunk,etot,rms,co,frac,&
349 deix,dent,sold,scur,runtime
352 ! if(.not.allocated(varsave)) allocate(varsave(maxvar,maxsave)) !(maxvar,maxsave)
353 !d write (iout,*) 'print_mc=',print_mc
356 !---------------------------------------------------------------------------
357 ! Initialize counters.
358 !---------------------------------------------------------------------------
359 ! Total number of generated confs.
361 ! Total number of moves. In general this won't be equal to the number of
362 ! attempted moves, because we may want to reject some "bad" confs just by
365 ! Total number of shift (nbond_move(1)), spike, crankshaft, three-bond,...
367 !el allocate(nbond_move(nres)) !(maxres)
372 ! Initialize total and accepted number of moves of various kind.
377 ! Total number of energy evaluations.
383 facee=1.0D0/(maxacc*delte)
385 ! Read entropy from previous simulations.
387 read (ientin,*) indminn,indmaxx,emin,emax
388 print *,'indminn=',indminn,' indmaxx=',indmaxx,' emin=',emin,&
390 do i=-max_ene,max_ene
391 entropy(i)=(emin+i*delte)*betbol
393 read (ientin,*) (ijunk,ejunk,entropy(i),i=indminn,indmaxx)
396 write (iout,*) 'indminn=',indminn,' indmaxx=',indmaxx,&
397 ' emin=',emin,' emax=',emax
398 write (iout,'(/a)') 'Initial entropy'
400 write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
403 ! Read the pool of conformations
405 !----------------------------------------------------------------------------
406 ! Entropy-sampling simulations with continually updated entropy
407 ! Loop thru simulations
408 !----------------------------------------------------------------------------
410 !----------------------------------------------------------------------------
411 ! Take a conformation from the pool
412 !----------------------------------------------------------------------------
418 write (iout,*) 'Took conformation',ii,' from the pool energy=',&
420 call var_to_geom(nvar,varia)
421 ! Print internal coordinates of the initial conformation
424 call gen_rand_conf(1,*20)
426 !----------------------------------------------------------------------------
427 ! Compute and print initial energies.
428 !----------------------------------------------------------------------------
431 allocate(nsave_part(nctasks))
432 if (MyID.eq.MasterID) then
440 write (iout,'(/80(1h*)/a,i2/80(1h*)/)') 'MCE iteration #',isweep
441 write (iout,'(/80(1h*)/a)') 'Initial energies:'
445 call enerprint(energia)
446 ! Minimize the energy of the first conformation.
448 call geom_to_var(nvar,varia)
449 call minimize(etot,varia,iretcode,nfun)
452 write (iout,'(/80(1h*)/a/80(1h*))') &
453 'Results of the first energy minimization:'
454 call enerprint(energia)
458 call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,kkk),&
462 call contact(.false.,ncont,icont,co)
463 frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
464 write (iout,'(a,f8.3,a,f8.3,a,f8.3)') &
465 'RMS deviation from the reference structure:',rms,&
466 ' % of native contacts:',frac*100,' contact order:',co
467 write (istat,'(i5,11(1pe14.5))') 0,&
468 (energia(print_order(i)),i=1,nprint_ene),etot,rms,frac,co
470 write (istat,'(i5,9(1pe14.5))') 0,&
471 (energia(print_order(i)),i=1,nprint_ene),etot
474 neneval=neneval+nfun+1
475 if (.not. ent_read) then
476 ! Initialize the entropy array
477 do i=-max_ene,max_ene
479 ! Uncomment the line below for actual entropic sampling (start with uniform
480 ! energy distribution).
482 ! Uncomment the line below for multicanonical sampling (start with Boltzmann
484 entropy(i)=(emin+i*delte)*betbol
488 write (iout,'(/a)') 'Initial entropy'
490 write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
494 call recv_stop_sig(Kwita)
495 if (whatsup.eq.1) then
496 call send_stop_sig(-2)
498 else if (whatsup.le.-2) then
500 else if (whatsup.eq.2) then
506 not_done = (iretcode.ne.11)
508 write (iout,'(/80(1h*)/20x,a/80(1h*))') &
509 'Enter Monte Carlo procedure.'
511 call briefout(0,etot)
516 indeold=(eold-emin)/delte
517 deix=eold-(emin+indeold*delte)
518 dent=entropy(indeold+1)-entropy(indeold)
519 !d write (iout,*) 'indeold=',indeold,' deix=',deix,' dent=',dent
520 !d write (*,*) 'Processor',MyID,' indeold=',indeold,' deix=',deix,
522 sold=entropy(indeold)+(dent/delte)*deix
524 write (iout,*) 'eold=',eold,' sold=',sold,' elowest=',etot
525 write (*,*) 'Processor',MyID,' eold=',eold,' sold=',sold,&
527 if (minim) call zapis(varia,etot)
529 ! NACC is the counter for the accepted conformations of a given processor
531 ! NACC_TOT counts the total number of accepted conformations
534 if (MyID.eq.MasterID) then
535 call receive_MCM_info
537 call send_MCM_info(2)
540 do iene=indminn,indmaxx
547 !----------------------------------------------------------------------------
553 if (print_mc.gt.0) write (iout,'(80(1h*)/20x,a,i7)') &
554 'Beginning iteration #',it
555 ! Initialize local counter.
556 ntrial=0 ! # of generated non-overlapping confs.
557 noverlap=0 ! # of overlapping confs.
559 do while (.not. accepted .and. WhatsUp.eq.0 .and. Kwita.eq.0)
561 ! Retrieve the angles of previously accepted conformation
565 !d write (iout,'(a)') 'Old variables:'
566 !d write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
567 call var_to_geom(nvar,varia)
573 ! Decide whether to generate a random conformation or perturb the old one
574 RandOrPert=ran_number(0.0D0,1.0D0)
575 if (RandOrPert.gt.RanFract) then
577 write (iout,'(a)') 'Perturbation-generated conformation.'
578 call perturb(error,lprint,MoveType,nbond,1.0D0)
580 if (MoveType.lt.1 .or. MoveType.gt.MaxMoveType) then
581 write (iout,'(/a,i7,a/)') 'Error - unknown MoveType=',&
582 MoveType,' returned from PERTURB.'
589 nstart_grow=iran_num(3,nres)
591 write (iout,'(2a,i3)') 'Random-generated conformation',&
592 ' - chain regrown from residue',nstart_grow
593 call gen_rand_conf(nstart_grow,*30)
595 call geom_to_var(nvar,varia)
596 !d write (iout,'(a)') 'New variables:'
597 !d write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
599 if (print_mc.gt.0) write (iout,'(a,i5,a,i10,a,i10)') &
600 'Processor',MyId,' trial move',ntrial,' total generated:',ngen
601 if (print_mc.gt.0) write (*,'(a,i5,a,i10,a,i10)') &
602 'Processor',MyId,' trial move',ntrial,' total generated:',ngen
605 ! call enerprint(energia(0))
606 ! write (iout,'(2(a,1pe14.5))') 'Etot=',Etot,' Elowest=',Elowest
607 if (etot-elowest.gt.overlap_cut) then
608 write (iout,'(a,i5,a,1pe14.5)') 'Iteration',it,&
609 ' Overlap detected in the current conf.; energy is',etot
613 if (noverlap.gt.maxoverlap) then
614 write (iout,'(a)') 'Too many overlapping confs.'
619 call minimize(etot,varia,iretcode,nfun)
620 !d write (iout,'(a)') 'Variables after minimization:'
621 !d write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
624 neneval=neneval+nfun+1
626 if (print_mc.gt.2) then
627 write (iout,'(a)') 'Total energies of trial conf:'
628 call enerprint(energia)
629 else if (print_mc.eq.1) then
630 write (iout,'(a,i6,a,1pe16.6)') &
631 'Trial conformation:',ngen,' energy:',etot
633 !--------------------------------------------------------------------------
635 !--------------------------------------------------------------------------
638 call accepting(etot,eold,scur,sold,varia,varold,&
643 if (elowest.gt.etot) elowest=etot
644 if (ehighest.lt.etot) ehighest=etot
645 moves_acc(MoveType)=moves_acc(MoveType)+1
646 if (MoveType.eq.1) then
647 nbond_acc(nbond)=nbond_acc(nbond)+1
649 ! Check against conformation repetitions.
650 irep=conf_comp(varia,etot)
651 #if defined(AIX) || defined(PGI)
652 open (istat,file=statname,position='append')
654 open (istat,file=statname,access='append')
658 call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,kkk),&
662 call contact(.false.,ncont,icont,co)
663 frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
665 write (iout,'(a,f8.3,a,f8.3,a,f8.3)') &
666 'RMS deviation from the reference structure:',rms,&
667 ' % of native contacts:',frac*100,' contact order:',co
669 write (istat,'(i5,11(1pe14.5))') it,&
670 (energia(print_order(i)),i=1,nprint_ene),etot,&
672 elseif (print_stat) then
673 write (istat,'(i5,10(1pe14.5))') it,&
674 (energia(print_order(i)),i=1,nprint_ene),etot
678 call statprint(nacc,nfun,iretcode,etot,elowest)
679 ! Print internal coordinates.
680 if (print_int) call briefout(nacc,etot)
682 if (MyID.ne.MasterID) then
683 call recv_stop_sig(Kwita)
684 !d print *,'Processor:',MyID,' STOP=',Kwita
686 call send_MCM_info(2)
688 call send_MCM_info(1)
692 ! Store the accepted conf. and its energy.
700 !d write (iout,*) 'Accepted conformation:'
701 !d write (iout,*) (rad2deg*varia(i),i=1,nphi)
702 if (minim) call zapis(varia,etot)
704 ener(i,nsave)=energia(i)
706 ener(n_ene+1,nsave)=etot
707 ener(n_ene+2,nsave)=frac
709 nminima(irep)=nminima(irep)+1.0D0
710 ! print *,'irep=',irep,' nminima=',nminima(irep)
712 if (Kwita.eq.0) call recv_stop_sig(kwita)
717 if (MyID.eq.MasterID) then
718 call receive_MCM_info
719 if (nacc_tot.ge.maxacc) accepted=.true.
722 if (ntrial.gt.maxtrial_iter .and. npool.gt.0) then
723 ! Take a conformation from the pool
728 write (iout,*) 'Iteration',it,' max. # of trials exceeded.'
730 'Take conformation',ii,' from the pool energy=',epool(ii)
732 write (iout,'(10f8.3)') (rad2deg*varia(i),i=1,nvar)
734 endif ! (ntrial.gt.maxtrial_iter .and. npool.gt.0)
738 if (MyID.eq.MasterID) then
739 call receive_MCM_info
741 if (Kwita.eq.0) call recv_stop_sig(kwita)
743 if (ovrtim()) WhatsUp=-1
744 !d write (iout,*) 'WhatsUp=',WhatsUp,' Kwita=',Kwita
745 not_done = (nacc_tot.lt.maxacc) .and. (WhatsUp.eq.0) &
747 !d write (iout,*) 'not_done=',not_done
750 print *,'Processor',MyID,&
751 ' has received STOP signal =',Kwita,' in EntSamp.'
752 !d print *,'not_done=',not_done
753 if (Kwita.lt.-1) WhatsUp=Kwita
754 else if (nacc_tot.ge.maxacc) then
755 print *,'Processor',MyID,' calls send_stop_sig,',&
756 ' because a sufficient # of confs. have been collected.'
757 !d print *,'not_done=',not_done
758 call send_stop_sig(-1)
759 else if (WhatsUp.eq.-1) then
760 print *,'Processor',MyID,&
761 ' calls send_stop_sig because of timeout.'
762 !d print *,'not_done=',not_done
763 call send_stop_sig(-2)
768 !-----------------------------------------------------------------
769 !... Construct energy histogram & update entropy
770 !-----------------------------------------------------------------
774 write (iout,*) 'Processor',MyID,&
775 ' is broadcasting ERROR-STOP signal.'
776 write (*,*) 'Processor',MyID,&
777 ' is broadcasting ERROR-STOP signal.'
778 call send_stop_sig(-3)
782 if (MyID.eq.MasterID) then
783 ! call receive_MCM_results
784 call receive_energies
787 if (esave(i).lt.elowest) elowest=esave(i)
788 if (esave(i).gt.ehighest) ehighest=esave(i)
790 write (iout,'(a,i10)') '# of accepted confs:',nacc_tot
791 write (iout,'(a,f10.5,a,f10.5)') 'Lowest energy:',elowest,&
792 ' Highest energy',ehighest
793 if (isweep.eq.1 .and. .not.ent_read) then
796 write (iout,*) 'EMAX=',emax
798 indmaxx=(ehighest-emin)/delte
801 do i=-max_ene,max_ene
802 entropy(i)=(emin+i*delte)*betbol
806 indmin=(elowest-emin)/delte
807 indmax=(ehighest-emin)/delte
808 if (indmin.lt.indminn) indminn=indmin
809 if (indmax.gt.indmaxx) indmaxx=indmax
811 write(iout,*)'indminn=',indminn,' indmaxx=',indmaxx
812 ! Construct energy histogram
814 inde=(esave(i)-emin)/delte
815 nhist(inde)=nhist(inde)+nminima(i)
817 ! Update entropy (density of states)
819 if (nhist(i).gt.0) then
820 entropy(i)=entropy(i)+dlog(nhist(i)+0.0D0)
824 !d entropy(i)=1.0D+10
826 write (iout,'(/80(1h*)/a,i2/80(1h*)/)') &
827 'End of macroiteration',isweep
828 write (iout,'(a,f10.5,a,f10.5)') 'Elowest=',elowest,&
829 ' Ehighest=',ehighest
830 write (iout,'(a)') 'Frequecies of minima'
832 write (iout,'(i5,f5.0,f10.5)') i,nminima(i),esave(i)
834 write (iout,'(/a)') 'Energy histogram'
836 write (iout,'(i5,2f10.5)') i,emin+i*delte,nhist(i)
838 write (iout,'(/a)') 'Entropy'
840 write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
842 !-----------------------------------------------------------------
843 !... End of energy histogram construction
844 !-----------------------------------------------------------------
846 entropy(-max_ene-4)=dfloat(indminn)
847 entropy(-max_ene-3)=dfloat(indmaxx)
848 entropy(-max_ene-2)=emin
849 entropy(-max_ene-1)=emax
851 !d print *,entname,ientout
852 open (ientout,file=entname,status='unknown')
853 write (ientout,'(2i5,2e25.17)') indminn,indmaxx,emin,emax
855 write (ientout,'(i5,f10.5,f20.15)') i,emin+i*delte,entropy(i)
859 write (iout,'(a)') 'Frequecies of minima'
861 write (iout,'(i5,f5.0,f10.5)') i,nminima(i),esave(i)
863 ! call send_MCM_results
865 call receive_MCM_update
866 indminn=entropy(-max_ene-4)
867 indmaxx=entropy(-max_ene-3)
868 emin=entropy(-max_ene-2)
869 emax=entropy(-max_ene-1)
870 write (iout,*) 'Received from master:'
871 write (iout,*) 'indminn=',indminn,' indmaxx=',indmaxx,&
872 ' emin=',emin,' emax=',emax
873 write (iout,'(/a)') 'Entropy'
875 write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
878 if (WhatsUp.lt.-1) return
880 if (ovrtim() .or. WhatsUp.lt.0) return
883 write (iout,'(/80(1h*)/20x,a)') 'Summary run statistics:'
884 call statprint(nacc,nfun,iretcode,etot,elowest)
886 'Statistics of multiple-bond motions. Total motions:'
887 write (iout,'(16i5)') (nbond_move(i),i=1,Nbm)
888 write (iout,'(a)') 'Accepted motions:'
889 write (iout,'(16i5)') (nbond_acc(i),i=1,Nbm)
890 !el write (iout,'(a,i10)') 'Number of chain regrowths:',nregrow
891 !el write (iout,'(a,i10)') 'Accepted chain regrowths:',nregrow_acc
893 !---------------------------------------------------------------------------
895 !---------------------------------------------------------------------------
899 if (isweep.eq.nsweep .and. it.ge.maxacc) &
900 write (iout,'(/80(1h*)/20x,a/80(1h*)/)') 'All iterations done.'
902 end subroutine entmcm
903 !-----------------------------------------------------------------------------
904 subroutine accepting(ecur,eold,scur,sold,x,xold,accepted)
906 use geometry_data, only: nphi
907 use energy_data, only: max_ene
908 ! implicit real*8 (a-h,o-z)
909 ! include 'DIMENSIONS'
910 ! include 'COMMON.MCM'
911 ! include 'COMMON.MCE'
912 ! include 'COMMON.IOUNITS'
913 ! include 'COMMON.VAR'
915 use MPI_data !include 'COMMON.INFO'
917 ! include 'COMMON.GEO'
918 real(kind=8) :: ecur,eold,xx,bol !,ran_number
919 real(kind=8),dimension(6*nres) :: x,xold !(maxvar) (maxvar=6*maxres)
920 real(kind=8) :: tole=1.0D-1, tola=5.0D0
923 !!! local variables - el
925 real(kind=8) :: scur,sold,xxh,deix,dent
927 ! Check if the conformation is similar.
928 !d write (iout,*) 'Enter ACCEPTING'
929 !d write (iout,*) 'Old PHI angles:'
930 !d write (iout,*) (rad2deg*xold(i),i=1,nphi)
931 !d write (iout,*) 'Current angles'
932 !d write (iout,*) (rad2deg*x(i),i=1,nphi)
933 !d ddif=dif_ang(nphi,x,xold)
934 !d write (iout,*) 'Angle norm:',ddif
935 !d write (iout,*) 'ecur=',ecur,' emax=',emax
936 if (ecur.gt.emax) then
939 write (iout,'(a)') 'Conformation rejected as too high in energy'
941 else if (dabs(ecur-eold).lt.tole .and. &
942 dif_ang(nphi,x,xold).lt.tola) then
945 write (iout,'(a)') 'Conformation rejected as too similar'
948 ! Else evaluate the entropy of the conf and compare it with that of the previous
950 indecur=(ecur-emin)/delte
951 if (iabs(indecur).gt.max_ene) then
952 write (iout,'(a,2i5)') &
953 'Accepting: Index out of range:',indecur
955 else if (indecur.eq.indmaxx) then
956 scur=entropy(indecur)
957 if (print_mc.gt.0) write (iout,*)'Energy boundary reached',&
958 indmaxx,indecur,entropy(indecur)
960 deix=ecur-(emin+indecur*delte)
961 dent=entropy(indecur+1)-entropy(indecur)
962 scur=entropy(indecur)+(dent/delte)*deix
964 !d print *,'Processor',MyID,' ecur=',ecur,' indecur=',indecur,
965 !d & ' scur=',scur,' eold=',eold,' sold=',sold
966 !d print *,'deix=',deix,' dent=',dent,' delte=',delte
967 if (print_mc.gt.1) then
968 write(iout,*)'ecur=',ecur,' indecur=',indecur,' scur=',scur
969 write(iout,*)'eold=',eold,' sold=',sold
971 if (scur.le.sold) then
974 ! Else carry out acceptance test
975 xx=ran_number(0.0D0,1.0D0)
977 if (xxh.gt.50.0D0) then
984 if (print_mc.gt.0) write (iout,'(a)') &
985 'Conformation accepted.'
988 if (print_mc.gt.0) write (iout,'(a)') &
989 'Conformation rejected.'
993 end subroutine accepting
994 !-----------------------------------------------------------------------------
997 use io_base, only:read_angles
998 ! implicit real*8 (a-h,o-z)
999 ! include 'DIMENSIONS'
1000 ! include 'COMMON.IOUNITS'
1001 ! include 'COMMON.GEO'
1002 ! include 'COMMON.MCM'
1003 ! include 'COMMON.MCE'
1004 ! include 'COMMON.VAR'
1005 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
1007 !!! local variables - el
1008 integer :: j,i,iconf
1010 print '(a)','Call READ_POOL'
1013 read (intin,'(i5,f10.5)',end=10,err=10) iconf,epool(npool)
1014 if (epool(npool).eq.0.0D0) goto 10
1015 call read_angles(intin,*10)
1016 call geom_to_var(nvar,xpool(1,npool))
1020 11 write (iout,'(a,i5)') 'Number of pool conformations:',npool
1021 if (print_mc.gt.2) then
1023 write (iout,'(a,i5,a,1pe14.5)') 'Pool conformation',i,' energy',&
1025 write (iout,'(10f8.3)') (rad2deg*xpool(j,i),j=1,nvar)
1027 endif ! (print_mc.gt.2)
1029 end subroutine read_pool
1030 !-----------------------------------------------------------------------------
1032 !-----------------------------------------------------------------------------
1033 subroutine monte_carlo
1037 use MPI_data, only:ifinish,nctasks,WhatsUp,MyID
1038 use control_data, only:refstr,MaxProcs
1040 use control, only:tcpu,ovrtim
1041 use regularize_, only:fitsq
1044 ! Does Boltzmann and entropic sampling without energy minimization
1045 ! implicit real*8 (a-h,o-z)
1046 ! include 'DIMENSIONS'
1047 ! include 'COMMON.IOUNITS'
1049 use MPI_data !include 'COMMON.INFO'
1051 ! include 'COMMON.GEO'
1052 ! include 'COMMON.CHAIN'
1053 ! include 'COMMON.MCM'
1054 ! include 'COMMON.MCE'
1055 ! include 'COMMON.CONTACTS'
1056 ! include 'COMMON.CONTROL'
1057 ! include 'COMMON.VAR'
1058 ! include 'COMMON.INTERACT'
1059 ! include 'COMMON.THREAD'
1060 ! include 'COMMON.NAMES'
1061 logical :: accepted,not_done,over,error,lprint !,ovrtim
1062 integer :: MoveType,nbond,nbins
1063 ! integer :: conf_comp
1064 real(kind=8) :: RandOrPert
1065 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
1066 real(kind=8) :: elowest,elowest1,ehighest,ehighest1,eold
1067 real(kind=8) :: przes(3),obr(3,3)
1068 real(kind=8),dimension(6*nres) :: varold !(maxvar) (maxvar=6*maxres)
1070 integer,dimension(-1:MaxMoveType+1,0:MaxProcs-1) :: moves1,moves_acc1 !(-1:MaxMoveType+1,0:MaxProcs-1)
1072 real(kind=8) :: etot_temp,etot_all(0:MaxProcs)
1073 external d_vadd,d_vmin,d_vmax
1074 real(kind=8),dimension(-max_ene:max_ene) :: entropy1,nhist1
1075 integer,dimension(nres*(MaxProcs+1)) :: nbond_move1,nbond_acc1
1076 integer,dimension(2) :: itemp
1078 real(kind=8),dimension(6*nres) :: var_lowest !(maxvar) (maxvar=6*maxres)
1079 real(kind=8),dimension(0:n_ene) :: energia,energia_ave
1081 !!! local variables - el
1082 integer :: i,j,it,ii,iproc,nacc,ISWEEP,nfun,indmax,indmin,ijunk,&
1083 Kwita,indeold,imdmax,inde,iretcode,nstart_grow,noverlap
1084 real(kind=8) :: facee,conste,ejunk,etot,sold,frac,runtime,&
1085 frac_ave,rms_ave,etot_ave,scur,from_pool,co,rms
1087 write(iout,'(a,i8,2x,a,f10.5)') &
1088 'pool_read_freq=',pool_read_freq,' pool_fraction=',pool_fraction
1089 open (istat,file=statname)
1093 facee=1.0D0/(maxacc*delte)
1094 ! Number of bins in energy histogram
1096 write (iout,*) 'NBINS=',nbins
1098 ! Read entropy from previous simulations.
1100 read (ientin,*) indminn,indmaxx,emin,emax
1101 print *,'indminn=',indminn,' indmaxx=',indmaxx,' emin=',emin,&
1103 do i=-max_ene,max_ene
1106 read (ientin,*) (ijunk,ejunk,entropy(i),i=indminn,indmaxx)
1109 write (iout,*) 'indminn=',indminn,' indmaxx=',indmaxx,&
1110 ' emin=',emin,' emax=',emax
1111 write (iout,'(/a)') 'Initial entropy'
1112 do i=indminn,indmaxx
1113 write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
1116 ! Read the pool of conformations
1120 !----------------------------------------------------------------------------
1121 ! Entropy-sampling simulations with continually updated entropy;
1122 ! set NSWEEP=1 for Boltzmann sampling.
1123 ! Loop thru simulations
1124 !----------------------------------------------------------------------------
1125 allocate(ifinish(nctasks))
1128 ! Initialize the IFINISH array.
1135 !---------------------------------------------------------------------------
1136 ! Initialize counters.
1137 !---------------------------------------------------------------------------
1138 ! Total number of generated confs.
1140 ! Total number of moves. In general this won't be equal to the number of
1141 ! attempted moves, because we may want to reject some "bad" confs just by
1144 ! Total number of shift (nbond_move(1)), spike, crankshaft, three-bond,...
1146 !el allocate(nbond_move(nres)) !(maxres)
1147 !el allocate(nbond_acc(nres)) !(maxres)
1153 ! Initialize total and accepted number of moves of various kind.
1158 ! Total number of energy evaluations.
1161 !----------------------------------------------------------------------------
1162 ! Take a conformation from the pool
1163 !----------------------------------------------------------------------------
1165 write (iout,*) 'emin=',emin,' emax=',emax
1166 if (npool.gt.0) then
1167 ii=iran_num(1,npool)
1169 varia(i)=xpool(i,ii)
1171 write (iout,*) 'Took conformation',ii,' from the pool energy=',&
1173 call var_to_geom(nvar,varia)
1174 ! Print internal coordinates of the initial conformation
1176 else if (isweep.gt.1) then
1177 if (eold.lt.emax) then
1183 varia(i)=var_lowest(i)
1186 call var_to_geom(nvar,varia)
1188 !----------------------------------------------------------------------------
1189 ! Compute and print initial energies.
1190 !----------------------------------------------------------------------------
1194 write (iout,'(/80(1h*)/a,i2/80(1h*)/)') 'MCE iteration #',isweep
1195 write (iout,'(/80(1h*)/a)') 'Initial energies:'
1197 call geom_to_var(nvar,varia)
1198 call etotal(energia)
1200 call enerprint(energia)
1202 call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,1),nsup,przes,&
1205 call contact(.false.,ncont,icont,co)
1206 frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
1207 write (iout,'(a,f8.3,a,f8.3,a,f8.3)') &
1208 'RMS deviation from the reference structure:',rms,&
1209 ' % of native contacts:',frac*100,' contact order',co
1210 write (istat,'(i10,16(1pe14.5))') 0,&
1211 (energia(print_order(i)),i=1,nprint_ene),&
1214 write (istat,'(i10,14(1pe14.5))') 0,&
1215 (energia(print_order(i)),i=1,nprint_ene),etot
1219 if (.not. ent_read) then
1220 ! Initialize the entropy array
1222 ! Collect total energies from other processors.
1225 call mp_gather(etot_temp,etot_all,8,MasterID,cgGroupID)
1226 if (MyID.eq.MasterID) then
1227 ! Get the lowest and the highest energy.
1228 print *,'MASTER: etot_temp: ',(etot_all(i),i=0,nprocs-1),&
1229 ' emin=',emin,' emax=',emax
1233 if (emin.gt.etot_all(i)) emin=etot_all(i)
1234 if (emax.lt.etot_all(i)) emax=etot_all(i)
1237 endif ! MyID.eq.MasterID
1240 print *,'Processor',MyID,' calls MP_BCAST to send/recv etot_all'
1241 call mp_bcast(etot_all(1),16,MasterID,cgGroupID)
1242 print *,'Processor',MyID,' MP_BCAST to send/recv etot_all ended'
1243 if (MyID.ne.MasterID) then
1244 print *,'Processor:',MyID,etot_all(1),etot_all(2),&
1245 etot_all(1),etot_all(2)
1248 endif ! MyID.ne.MasterID
1249 write (iout,*) 'After MP_GATHER etot_temp=',&
1250 etot_temp,' emin=',emin
1258 ! Multicanonical sampling - start from Boltzmann distribution
1259 do i=-max_ene,max_ene
1260 entropy(i)=(emin+i*delte)*betbol
1263 ! Entropic sampling - start from uniform distribution of the density of states
1264 do i=-max_ene,max_ene
1268 write (iout,'(/a)') 'Initial entropy'
1269 do i=indminn,indmaxx
1270 write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
1272 if (isweep.eq.1) then
1276 indmaxx=indminn+nbins
1279 endif ! .not. ent_read
1281 call recv_stop_sig(Kwita)
1282 if (whatsup.eq.1) then
1283 call send_stop_sig(-2)
1285 else if (whatsup.le.-2) then
1287 else if (whatsup.eq.2) then
1295 write (iout,'(/80(1h*)/20x,a/80(1h*))') &
1296 'Enter Monte Carlo procedure.'
1298 call briefout(0,etot)
1303 call entropia(eold,sold,indeold)
1304 ! NACC is the counter for the accepted conformations of a given processor
1306 ! NACC_TOT counts the total number of accepted conformations
1309 !----------------------------------------------------------------------------
1310 ! Zero out average energies
1312 energia_ave(i)=0.0d0
1314 ! Initialize energy histogram
1315 do i=-max_ene,max_ene
1318 ! Zero out iteration counter.
1323 ! Begin MC iteration loop.
1326 ! Initialize local counter.
1327 ntrial=0 ! # of generated non-overlapping confs.
1328 noverlap=0 ! # of overlapping confs.
1330 do while (.not. accepted .and. WhatsUp.eq.0 .and. Kwita.eq.0)
1332 ! Retrieve the angles of previously accepted conformation
1336 call var_to_geom(nvar,varia)
1337 ! Rebuild the chain.
1342 ! Decide whether to take a conformation from the pool or generate/perturb one
1344 from_pool=ran_number(0.0D0,1.0D0)
1345 if (npool.gt.0 .and. from_pool.lt.pool_fraction) then
1346 ! Throw a dice to choose the conformation from the pool
1347 ii=iran_num(1,npool)
1349 varia(i)=xpool(i,ii)
1351 call var_to_geom(nvar,varia)
1354 !d write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
1355 if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) &
1356 write (iout,'(a,i3,a,f10.5)') &
1357 'Try conformation',ii,' from the pool energy=',epool(ii)
1359 moves(-1)=moves(-1)+1
1361 ! Decide whether to generate a random conformation or perturb the old one
1362 RandOrPert=ran_number(0.0D0,1.0D0)
1363 if (RandOrPert.gt.RanFract) then
1364 if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) &
1365 write (iout,'(a)') 'Perturbation-generated conformation.'
1366 call perturb(error,lprint,MoveType,nbond,0.1D0)
1368 if (MoveType.lt.1 .or. MoveType.gt.MaxMoveType) then
1369 write (iout,'(/a,i7,a/)') 'Error - unknown MoveType=',&
1370 MoveType,' returned from PERTURB.'
1377 nstart_grow=iran_num(3,nres)
1378 if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) &
1379 write (iout,'(2a,i3)') 'Random-generated conformation',&
1380 ' - chain regrown from residue',nstart_grow
1381 call gen_rand_conf(nstart_grow,*30)
1383 call geom_to_var(nvar,varia)
1385 !d write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
1387 if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) &
1388 write (iout,'(a,i5,a,i10,a,i10)') &
1389 'Processor',MyId,' trial move',ntrial,' total generated:',ngen
1390 if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) &
1391 write (*,'(a,i5,a,i10,a,i10)') &
1392 'Processor',MyId,' trial move',ntrial,' total generated:',ngen
1393 call etotal(energia)
1396 !d call enerprint(energia(0))
1397 !d write(iout,*)'it=',it,' etot=',etot
1398 if (etot-elowest.gt.overlap_cut) then
1399 if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) &
1400 write (iout,'(a,i5,a,1pe14.5)') 'Iteration',it,&
1401 ' Overlap detected in the current conf.; energy is',etot
1404 if (noverlap.gt.maxoverlap) then
1405 write (iout,'(a)') 'Too many overlapping confs.'
1409 !--------------------------------------------------------------------------
1410 !... Acceptance test
1411 !--------------------------------------------------------------------------
1414 call accept_mc(it,etot,eold,scur,sold,varia,varold,accepted)
1418 if (elowest.gt.etot) then
1421 var_lowest(i)=varia(i)
1424 if (ehighest.lt.etot) ehighest=etot
1425 moves_acc(MoveType)=moves_acc(MoveType)+1
1426 if (MoveType.eq.1) then
1427 nbond_acc(nbond)=nbond_acc(nbond)+1
1429 ! Compare with reference structure.
1431 call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,1),&
1432 nsup,przes,obr,non_conv)
1434 call contact(.false.,ncont,icont,co)
1435 frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
1438 ! Periodically save average energies and confs.
1441 energia_ave(i)=energia_ave(i)+energia(i)
1443 moves(MaxMoveType+1)=nmove
1444 moves_acc(MaxMoveType+1)=nacc
1445 IF ((it/save_frequency)*save_frequency.eq.it) THEN
1447 energia_ave(i)=energia_ave(i)/save_frequency
1449 etot_ave=energia_ave(0)
1451 ! open (istat,file=statname,position='append')
1453 ! open (istat,file=statname,access='append')
1455 if (print_mc.gt.0) &
1456 write (iout,'(80(1h*)/20x,a,i20)') &
1458 if (refstr .and. print_mc.gt.0) then
1459 write (iout,'(a,f8.3,a,f8.3,a,f8.3)') &
1460 'RMS deviation from the reference structure:',rms,&
1461 ' % of native contacts:',frac*100,' contact order:',co
1463 if (print_stat) then
1465 write (istat,'(i10,10(1pe14.5))') it,&
1466 (energia_ave(print_order(i)),i=1,nprint_ene),&
1467 etot_ave,rms_ave,frac_ave
1469 write (istat,'(i10,10(1pe14.5))') it,&
1470 (energia_ave(print_order(i)),i=1,nprint_ene),&
1475 if (print_mc.gt.0) &
1476 call statprint(nacc,nfun,iretcode,etot,elowest)
1477 ! Print internal coordinates.
1478 if (print_int) call briefout(nacc,etot)
1480 energia_ave(i)=0.0d0
1482 ENDIF ! ( (it/save_frequency)*save_frequency.eq.it)
1484 inde=icialosc((etot-emin)/delte)
1485 nhist(inde)=nhist(inde)+1.0D0
1487 if ( (it/message_frequency)*message_frequency.eq.it &
1488 .and. (MyID.ne.MasterID) ) then
1489 call recv_stop_sig(Kwita)
1490 call send_MCM_info(message_frequency)
1493 ! Store the accepted conf. and its energy.
1500 if (Kwita.eq.0) call recv_stop_sig(kwita)
1505 if (MyID.eq.MasterID .and. &
1506 (it/message_frequency)*message_frequency.eq.it) then
1507 call receive_MC_info
1508 if (nacc_tot.ge.maxacc) accepted=.true.
1511 ! if ((ntrial.gt.maxtrial_iter
1512 ! & .or. (it/pool_read_freq)*pool_read_freq.eq.it)
1513 ! & .and. npool.gt.0) then
1514 ! Take a conformation from the pool
1515 ! ii=iran_num(1,npool)
1517 ! varold(i)=xpool(i,ii)
1519 ! if (ntrial.gt.maxtrial_iter)
1520 ! & write (iout,*) 'Iteration',it,' max. # of trials exceeded.'
1522 ! & 'Take conformation',ii,' from the pool energy=',epool(ii)
1523 ! if (print_mc.gt.2)
1524 ! & write (iout,'(10f8.3)') (rad2deg*varold(i),i=1,nvar)
1527 ! call entropia(eold,sold,indeold)
1529 ! endif ! (ntrial.gt.maxtrial_iter .and. npool.gt.0)
1533 if (MyID.eq.MasterID .and. &
1534 (it/message_frequency)*message_frequency.eq.it) then
1535 call receive_MC_info
1537 if (Kwita.eq.0) call recv_stop_sig(kwita)
1539 if (ovrtim()) WhatsUp=-1
1540 !d write (iout,*) 'WhatsUp=',WhatsUp,' Kwita=',Kwita
1541 not_done = (nacc_tot.lt.maxacc) .and. (WhatsUp.eq.0) &
1543 !d write (iout,*) 'not_done=',not_done
1545 if (Kwita.lt.0) then
1546 print *,'Processor',MyID,&
1547 ' has received STOP signal =',Kwita,' in EntSamp.'
1548 !d print *,'not_done=',not_done
1549 if (Kwita.lt.-1) WhatsUp=Kwita
1550 if (MyID.ne.MasterID) call send_MCM_info(-1)
1551 else if (nacc_tot.ge.maxacc) then
1552 print *,'Processor',MyID,' calls send_stop_sig,',&
1553 ' because a sufficient # of confs. have been collected.'
1554 !d print *,'not_done=',not_done
1555 call send_stop_sig(-1)
1556 if (MyID.ne.MasterID) call send_MCM_info(-1)
1557 else if (WhatsUp.eq.-1) then
1558 print *,'Processor',MyID,&
1559 ' calls send_stop_sig because of timeout.'
1560 !d print *,'not_done=',not_done
1561 call send_stop_sig(-2)
1562 if (MyID.ne.MasterID) call send_MCM_info(-1)
1567 !-----------------------------------------------------------------
1568 !... Construct energy histogram & update entropy
1569 !-----------------------------------------------------------------
1573 write (iout,*) 'Processor',MyID,&
1574 ' is broadcasting ERROR-STOP signal.'
1575 write (*,*) 'Processor',MyID,&
1576 ' is broadcasting ERROR-STOP signal.'
1577 call send_stop_sig(-3)
1578 if (MyID.ne.MasterID) call send_MCM_info(-1)
1581 write (iout,'(/a)') 'Energy histogram'
1583 write (iout,'(i5,2f20.5)') i,emin+i*delte,nhist(i)
1586 ! Wait until every processor has sent complete MC info.
1587 if (MyID.eq.MasterID) then
1590 ! write (*,*) 'The IFINISH array:'
1591 ! write (*,*) (ifinish(i),i=1,nctasks)
1594 not_done=not_done.or.(ifinish(i).ge.0)
1596 if (not_done) call receive_MC_info
1599 ! Make collective histogram from the work of all processors.
1600 msglen=(2*max_ene+1)*8
1602 'Processor',MyID,' calls MP_REDUCE to send/receive histograms',&
1604 call mp_reduce(nhist,nhist1,msglen,MasterID,d_vadd,&
1606 print *,'Processor',MyID,' MP_REDUCE accomplished for histogr.'
1607 do i=-max_ene,max_ene
1610 ! Collect min. and max. energy
1612 'Processor',MyID,' calls MP_REDUCE to send/receive energy borders'
1613 call mp_reduce(elowest,elowest1,8,MasterID,d_vmin,cgGroupID)
1614 call mp_reduce(ehighest,ehighest1,8,MasterID,d_vmax,cgGroupID)
1615 print *,'Processor',MyID,' MP_REDUCE accomplished for energies.'
1616 IF (MyID.eq.MasterID) THEN
1620 write (iout,'(a,i10)') '# of accepted confs:',nacc_tot
1621 write (iout,'(a,f10.5,a,f10.5)') 'Lowest energy:',elowest,&
1622 ' Highest energy',ehighest
1623 indmin=icialosc((elowest-emin)/delte)
1624 imdmax=icialosc((ehighest-emin)/delte)
1625 if (indmin.lt.indminn) then
1626 emax=emin+indmin*delte+e_up
1627 indmaxx=indmin+nbins
1630 if (.not.ent_read) ent_read=.true.
1631 write(iout,*)'indminn=',indminn,' indmaxx=',indmaxx
1632 ! Update entropy (density of states)
1634 if (nhist(i).gt.0) then
1635 entropy(i)=entropy(i)+dlog(nhist(i)+0.0D0)
1638 write (iout,'(/80(1h*)/a,i2/80(1h*)/)') &
1639 'End of macroiteration',isweep
1640 write (iout,'(a,f10.5,a,f10.5)') 'Elowest=',elowest,&
1641 ' Ehighest=',ehighest
1642 write (iout,'(/a)') 'Energy histogram'
1643 do i=indminn,indmaxx
1644 write (iout,'(i5,2f20.5)') i,emin+i*delte,nhist(i)
1646 write (iout,'(/a)') 'Entropy'
1647 do i=indminn,indmaxx
1648 write (iout,'(i5,2f20.5)') i,emin+i*delte,entropy(i)
1650 !-----------------------------------------------------------------
1651 !... End of energy histogram construction
1652 !-----------------------------------------------------------------
1655 if (.not. ent_read) ent_read=.true.
1656 ENDIF ! MyID .eq. MaterID
1657 if (MyID.eq.MasterID) then
1661 print *,'before mp_bcast processor',MyID,' indminn=',indminn,&
1662 ' indmaxx=',indmaxx,' itemp=',itemp(1),itemp(2)
1663 call mp_bcast(itemp(1),8,MasterID,cgGroupID)
1664 call mp_bcast(emax,8,MasterID,cgGroupID)
1665 print *,'after mp_bcast processor',MyID,' indminn=',indminn,&
1666 ' indmaxx=',indmaxx,' itemp=',itemp(1),itemp(2)
1667 if (MyID .ne. MasterID) then
1671 msglen=(indmaxx-indminn+1)*8
1672 print *,'processor',MyID,' calling mp_bcast msglen=',msglen,&
1673 ' indminn=',indminn,' indmaxx=',indmaxx,' isweep=',isweep
1674 call mp_bcast(entropy(indminn),msglen,MasterID,cgGroupID)
1675 IF(MyID.eq.MasterID .and. .not. ovrtim() .and. WhatsUp.ge.0)THEN
1676 open (ientout,file=entname,status='unknown')
1677 write (ientout,'(2i5,2e25.17)') indminn,indmaxx,emin,emax
1678 do i=indminn,indmaxx
1679 write (ientout,'(i5,f10.5,f20.15)') i,emin+i*delte,entropy(i)
1683 write (iout,*) 'Received from master:'
1684 write (iout,*) 'indminn=',indminn,' indmaxx=',indmaxx,&
1685 ' emin=',emin,' emax=',emax
1686 write (iout,'(/a)') 'Entropy'
1687 do i=indminn,indmaxx
1688 write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
1690 ENDIF ! MyID.eq.MasterID
1691 print *,'Processor',MyID,' calls MP_GATHER'
1692 call mp_gather(nbond_move,nbond_move1,4*Nbm,MasterID,&
1694 call mp_gather(nbond_acc,nbond_acc1,4*Nbm,MasterID,&
1696 print *,'Processor',MyID,' MP_GATHER call accomplished'
1697 if (MyID.eq.MasterID) then
1699 write (iout,'(/80(1h*)/20x,a)') 'Summary run statistics:'
1700 call statprint(nacc_tot,nfun,iretcode,etot,elowest)
1701 write (iout,'(a)') &
1702 'Statistics of multiple-bond motions. Total motions:'
1703 write (iout,'(8i10)') (nbond_move(i),i=1,Nbm)
1704 write (iout,'(a)') 'Accepted motions:'
1705 write (iout,'(8i10)') (nbond_acc(i),i=1,Nbm)
1707 write (iout,'(a)') &
1708 'Statistics of multi-bond moves of respective processors:'
1712 nbond_move(i)=nbond_move(i)+nbond_move1(ind)
1713 nbond_acc(i)=nbond_acc(i)+nbond_acc1(ind)
1717 write (iout,*) 'Processor',iproc,' nbond_move:', &
1718 (nbond_move1(iproc*nbm+i),i=1,Nbm),&
1719 ' nbond_acc:',(nbond_acc1(iproc*nbm+i),i=1,Nbm)
1722 call mp_gather(moves,moves1,4*(MaxMoveType+3),MasterID,&
1724 call mp_gather(moves_acc,moves_acc1,4*(MaxMoveType+3),&
1726 if (MyID.eq.MasterID) then
1728 do i=-1,MaxMoveType+1
1729 moves(i)=moves(i)+moves1(i,iproc)
1730 moves_acc(i)=moves_acc(i)+moves_acc1(i,iproc)
1734 do i=0,MaxMoveType+1
1735 nmove=nmove+moves(i)
1738 write (iout,*) 'Processor',iproc,' moves',&
1739 (moves1(i,iproc),i=0,MaxMoveType+1),&
1740 ' moves_acc:',(moves_acc1(i,iproc),i=0,MaxMoveType+1)
1744 open (ientout,file=entname,status='unknown')
1745 write (ientout,'(2i5,2e25.17)') indminn,indmaxx,emin,emax
1746 do i=indminn,indmaxx
1747 write (ientout,'(i5,f10.5,f20.15)') i,emin+i*delte,entropy(i)
1751 write (iout,'(/80(1h*)/20x,a)') 'Summary run statistics:'
1752 call statprint(nacc_tot,nfun,iretcode,etot,elowest)
1753 write (iout,'(a)') &
1754 'Statistics of multiple-bond motions. Total motions:'
1755 write (iout,'(8i10)') (nbond_move(i),i=1,Nbm)
1756 write (iout,'(a)') 'Accepted motions:'
1757 write (iout,'(8i10)') (nbond_acc(i),i=1,Nbm)
1758 if (ovrtim() .or. WhatsUp.lt.0) return
1760 !---------------------------------------------------------------------------
1762 !---------------------------------------------------------------------------
1766 if (isweep.eq.nsweep .and. it.ge.maxacc) &
1767 write (iout,'(/80(1h*)/20x,a/80(1h*)/)') 'All iterations done.'
1769 end subroutine monte_carlo
1770 !-----------------------------------------------------------------------------
1771 subroutine accept_mc(it,ecur,eold,scur,sold,x,xold,accepted)
1773 ! implicit real*8 (a-h,o-z)
1774 ! include 'DIMENSIONS'
1775 ! include 'COMMON.MCM'
1776 ! include 'COMMON.MCE'
1777 ! include 'COMMON.IOUNITS'
1778 ! include 'COMMON.VAR'
1780 use MPI_data !include 'COMMON.INFO'
1782 ! include 'COMMON.GEO'
1783 real(kind=8) :: ecur,eold,xx,bol
1784 real(kind=8),dimension(6*nres) :: x,xold !(maxvar) (maxvar=6*maxres)
1788 integer :: it,indecur
1789 real(kind=8) :: scur,sold,xxh
1790 ! Check if the conformation is similar.
1791 !d write (iout,*) 'Enter ACCEPTING'
1792 !d write (iout,*) 'Old PHI angles:'
1793 !d write (iout,*) (rad2deg*xold(i),i=1,nphi)
1794 !d write (iout,*) 'Current angles'
1795 !d write (iout,*) (rad2deg*x(i),i=1,nphi)
1796 !d ddif=dif_ang(nphi,x,xold)
1797 !d write (iout,*) 'Angle norm:',ddif
1798 !d write (iout,*) 'ecur=',ecur,' emax=',emax
1799 if (ecur.gt.emax) then
1801 if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) &
1802 write (iout,'(a)') 'Conformation rejected as too high in energy'
1805 ! Else evaluate the entropy of the conf and compare it with that of the previous
1807 call entropia(ecur,scur,indecur)
1808 !d print *,'Processor',MyID,' ecur=',ecur,' indecur=',indecur,
1809 !d & ' scur=',scur,' eold=',eold,' sold=',sold
1810 !d print *,'deix=',deix,' dent=',dent,' delte=',delte
1811 if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) then
1812 write(iout,*)'it=',it,'ecur=',ecur,' indecur=',indecur,&
1814 write(iout,*)'eold=',eold,' sold=',sold
1816 if (scur.le.sold) then
1819 ! Else carry out acceptance test
1820 xx=ran_number(0.0D0,1.0D0)
1822 if (xxh.gt.50.0D0) then
1829 if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) &
1830 write (iout,'(a)') 'Conformation accepted.'
1833 if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) &
1834 write (iout,'(a)') 'Conformation rejected.'
1838 end subroutine accept_mc
1839 !-----------------------------------------------------------------------------
1840 integer function icialosc(x)
1843 if (x.lt.0.0D0) then
1849 end function icialosc
1850 !-----------------------------------------------------------------------------
1851 subroutine entropia(ecur,scur,indecur)
1853 use energy_data, only: max_ene
1854 ! implicit real*8 (a-h,o-z)
1855 ! include 'DIMENSIONS'
1856 ! include 'COMMON.MCM'
1857 ! include 'COMMON.MCE'
1858 ! include 'COMMON.IOUNITS'
1859 real(kind=8) :: ecur,scur,deix,dent
1860 integer :: indecur,it !???el
1862 indecur=icialosc((ecur-emin)/delte)
1863 if (iabs(indecur).gt.max_ene) then
1864 if ((it/print_freq)*it.eq.it) write (iout,'(a,2i5)') &
1865 'Accepting: Index out of range:',indecur
1867 else if (indecur.ge.indmaxx) then
1868 scur=entropy(indecur)
1869 if (print_mc.gt.0 .and. (it/print_freq)*it.eq.it) &
1870 write (iout,*)'Energy boundary reached',&
1871 indmaxx,indecur,entropy(indecur)
1873 deix=ecur-(emin+indecur*delte)
1874 dent=entropy(indecur+1)-entropy(indecur)
1875 scur=entropy(indecur)+(dent/delte)*deix
1878 end subroutine entropia
1879 !-----------------------------------------------------------------------------
1881 !-----------------------------------------------------------------------------
1882 subroutine mcm_setup
1885 ! implicit real*8 (a-h,o-z)
1886 ! include 'DIMENSIONS'
1887 ! include 'COMMON.IOUNITS'
1888 ! include 'COMMON.MCM'
1889 ! include 'COMMON.CONTROL'
1890 ! include 'COMMON.INTERACT'
1891 ! include 'COMMON.NAMES'
1892 ! include 'COMMON.CHAIN'
1893 ! include 'COMMON.VAR'
1895 !!! local variables - el
1896 integer :: i,i1,i2,it1,it2,ngly,mmm,maxwinlen
1898 ! Set up variables used in MC/MCM.
1900 ! allocate(sumpro_bond(0:nres)) !(0:maxres)
1902 write (iout,'(80(1h*)/20x,a/80(1h*))') 'MCM control parameters:'
1903 write (iout,'(5(a,i7))') 'Maxacc:',maxacc,' MaxTrial:',MaxTrial,&
1904 ' MaxRepm:',MaxRepm,' MaxGen:',MaxGen,' MaxOverlap:',MaxOverlap
1905 write (iout,'(4(a,f8.1)/2(a,i3))') &
1906 'Tmin:',Tmin,' Tmax:',Tmax,' TstepH:',TstepH,&
1907 ' TstepC:',TstepC,'NstepH:',NstepH,' NstepC:',NstepC
1908 if (nwindow.gt.0) then
1909 write (iout,'(a)') 'Perturbation windows:'
1915 write (iout,'(a,i3,a,i3,a,i3)') restyp(it1),i1,restyp(it2),i2,&
1919 ! Rbolt=8.3143D-3*2.388459D-01 kcal/(mol*K)
1921 ! Number of "end bonds".
1924 print *,'koniecl=',koniecl
1925 write (iout,'(a)') 'Probabilities of move types:'
1926 write (*,'(a)') 'Probabilities of move types:'
1928 write (iout,'(a,f10.5)') MovTypID(i),&
1929 sumpro_type(i)-sumpro_type(i-1)
1930 write (*,'(a,f10.5)') MovTypID(i),&
1931 sumpro_type(i)-sumpro_type(i-1)
1934 ! Maximum length of N-bond segment to be moved
1935 ! nbm=nres-1-(2*koniecl-1)
1936 if (nwindow.gt.0) then
1939 if (winlen(i).gt.maxwinlen) maxwinlen=winlen(i)
1941 nbm=min0(maxwinlen,6)
1942 write (iout,'(a,i3,a,i3)') 'Nbm=',Nbm,' Maxwinlen=',Maxwinlen
1946 sumpro_bond(0)=0.0D0
1947 sumpro_bond(1)=0.0D0
1949 sumpro_bond(i)=sumpro_bond(i-1)+1.0D0/dfloat(i)
1951 write (iout,'(a)') 'The SumPro_Bond array:'
1952 write (iout,'(8f10.5)') (sumpro_bond(i),i=1,nbm)
1953 write (*,'(a)') 'The SumPro_Bond array:'
1954 write (*,'(8f10.5)') (sumpro_bond(i),i=1,nbm)
1955 ! Maximum number of side chains moved simultaneously
1956 ! print *,'nnt=',nnt,' nct=',nct
1959 if (itype(i).eq.10) ngly=ngly+1
1963 MaxSideMove=min0((nct-nnt+1)/2,mmm)
1965 ! print *,'MaxSideMove=',MaxSideMove
1966 ! Max. number of generated confs (not used at present).
1968 ! Set initial temperature
1970 betbol=1.0D0/(Rbol*Tcur)
1971 write (iout,'(a,f8.1,a,f10.5)') 'Initial temperature:',Tcur,&
1973 write (iout,*) 'RanFract=',ranfract
1975 end subroutine mcm_setup
1976 !-----------------------------------------------------------------------------
1978 subroutine do_mcm(i_orig)
1982 use MPI_data, only:Whatsup
1983 use control_data, only:refstr,minim,iprint
1985 use control, only:tcpu,ovrtim
1986 use regularize_, only:fitsq
1988 use minimm, only:minimize
1989 ! Monte-Carlo-with-Minimization calculations - serial code.
1990 ! implicit real*8 (a-h,o-z)
1991 ! include 'DIMENSIONS'
1992 ! include 'COMMON.IOUNITS'
1993 ! include 'COMMON.GEO'
1994 ! include 'COMMON.CHAIN'
1995 ! include 'COMMON.MCM'
1996 ! include 'COMMON.CONTACTS'
1997 ! include 'COMMON.CONTROL'
1998 ! include 'COMMON.VAR'
1999 ! include 'COMMON.INTERACT'
2000 ! include 'COMMON.CACHE'
2001 !rc include 'COMMON.DEFORM'
2002 !rc include 'COMMON.DEFORM1'
2003 ! include 'COMMON.NAMES'
2004 logical :: accepted,over,error,lprint,not_done,my_conf,&
2005 enelower,non_conv !,ovrtim
2006 integer :: MoveType,nbond !,conf_comp
2007 integer,dimension(max_cache) :: ifeed
2008 real(kind=8),dimension(6*nres) :: varia,varold !(maxvar) (maxvar=6*maxres)
2009 real(kind=8) :: elowest,eold
2010 real(kind=8) :: przes(3),obr(3,3)
2011 real(kind=8) :: energia(0:n_ene)
2012 real(kind=8) :: coord1(6*nres,max_thread2),enetb1(max_threadss) !el
2013 !!! local variables - el
2014 integer :: i,nf,nacc,it,nout,j,i_orig,nfun,Kwita,iretcode,&
2015 noverlap,nstart_grow,irepet,n_thr,ii
2016 real(kind=8) :: etot,frac,rms,co,RandOrPert,&
2018 !---------------------------------------------------------------------------
2019 ! Initialize counters.
2020 !---------------------------------------------------------------------------
2021 ! Total number of generated confs.
2023 ! Total number of moves. In general this won't be equal to the number of
2024 ! attempted moves, because we may want to reject some "bad" confs just by
2027 ! Total number of temperature jumps.
2029 ! Total number of shift (nbond_move(1)), spike, crankshaft, three-bond,...
2031 ! if(.not.allocated(varsave)) allocate(varsave(maxvar,maxsave)) !(maxvar,maxsave)
2032 ! allocate(nbond_move(nres)) !(maxres)
2038 ! Initialize total and accepted number of moves of various kind.
2043 ! Total number of energy evaluations.
2048 write (iout,*) 'RanFract=',RanFract
2053 !----------------------------------------------------------------------------
2054 ! Compute and print initial energies.
2055 !----------------------------------------------------------------------------
2057 write (iout,'(/80(1h*)/a)') 'Initial energies:'
2061 call etotal(energia)
2063 ! Minimize the energy of the first conformation.
2065 call geom_to_var(nvar,varia)
2066 ! write (iout,*) 'The VARIA array'
2067 ! write (iout,'(8f10.4)') (rad2deg*varia(i),i=1,nvar)
2068 call minimize(etot,varia,iretcode,nfun)
2069 call var_to_geom(nvar,varia)
2071 write (iout,*) 'etot from MINIMIZE:',etot
2072 ! write (iout,*) 'Tha VARIA array'
2073 ! write (iout,'(8f10.4)') (rad2deg*varia(i),i=1,nvar)
2075 call etotal(energia)
2077 call enerprint(energia)
2080 call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,1),nsup,przes,& !el cref(1,nstart_sup)
2083 call contact(.false.,ncont,icont,co)
2084 frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
2085 write (iout,'(a,f8.3,a,f8.3,a,f8.3)') &
2086 'RMS deviation from the reference structure:',rms,&
2087 ' % of native contacts:',frac*100,' contact order:',co
2089 write (istat,'(i5,17(1pe14.5))') 0,&
2090 (energia(print_order(i)),i=1,nprint_ene),&
2093 if (print_stat) write (istat,'(i5,16(1pe14.5))') 0,&
2094 (energia(print_order(i)),i=1,nprint_ene),etot
2096 if (print_stat) close(istat)
2097 neneval=neneval+nfun+1
2098 write (iout,'(/80(1h*)/20x,a/80(1h*))') &
2099 'Enter Monte Carlo procedure.'
2102 call briefout(0,etot)
2109 call zapis(varia,etot)
2110 nacc=0 ! total # of accepted confs of the current processor.
2111 nacc_tot=0 ! total # of accepted confs of all processors.
2113 not_done = (iretcode.ne.11)
2115 !----------------------------------------------------------------------------
2117 !----------------------------------------------------------------------------
2122 write (iout,'(80(1h*)/20x,a,i7)') &
2123 'Beginning iteration #',it
2124 ! Initialize local counter.
2125 ntrial=0 ! # of generated non-overlapping confs.
2127 do while (.not. accepted)
2129 ! Retrieve the angles of previously accepted conformation
2130 noverlap=0 ! # of overlapping confs.
2134 call var_to_geom(nvar,varia)
2135 ! Rebuild the chain.
2137 ! Heat up the system, if necessary.
2139 ! If temperature cannot be further increased, stop.
2144 !d write (iout,'(a)') 'Old variables:'
2145 !d write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
2146 ! Decide whether to generate a random conformation or perturb the old one
2147 RandOrPert=ran_number(0.0D0,1.0D0)
2148 if (RandOrPert.gt.RanFract) then
2149 if (print_mc.gt.0) &
2150 write (iout,'(a)') 'Perturbation-generated conformation.'
2151 call perturb(error,lprint,MoveType,nbond,1.0D0)
2153 if (MoveType.lt.1 .or. MoveType.gt.MaxMoveType) then
2154 write (iout,'(/a,i7,a/)') 'Error - unknown MoveType=',&
2155 MoveType,' returned from PERTURB.'
2162 nstart_grow=iran_num(3,nres)
2163 if (print_mc.gt.0) &
2164 write (iout,'(2a,i3)') 'Random-generated conformation',&
2165 ' - chain regrown from residue',nstart_grow
2166 call gen_rand_conf(nstart_grow,*30)
2168 call geom_to_var(nvar,varia)
2169 !d write (iout,'(a)') 'New variables:'
2170 !d write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
2173 call etotal(energia)
2175 ! call enerprint(energia(0))
2176 ! write (iout,'(2(a,1pe14.5))') 'Etot=',Etot,' Elowest=',Elowest
2177 if (etot-elowest.gt.overlap_cut) then
2178 if(iprint.gt.1.or.etot.lt.1d20) &
2179 write (iout,'(a,1pe14.5)') &
2180 'Overlap detected in the current conf.; energy is',etot
2184 if (noverlap.gt.maxoverlap) then
2185 write (iout,'(a)') 'Too many overlapping confs.'
2190 call minimize(etot,varia,iretcode,nfun)
2191 !d write (iout,*) 'etot from MINIMIZE:',etot
2192 !d write (iout,'(a)') 'Variables after minimization:'
2193 !d write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
2195 call etotal(energia)
2197 neneval=neneval+nfun+2
2199 ! call enerprint(energia(0))
2200 write (iout,'(a,i6,a,1pe16.6)') 'Conformation:',ngen,&
2202 !--------------------------------------------------------------------------
2203 !... Do Metropolis test
2204 !--------------------------------------------------------------------------
2208 if (WhatsUp.eq.0 .and. Kwita.eq.0) then
2209 call metropolis(nvar,varia,varold,etot,eold,accepted,&
2212 write (iout,*) 'My_Conf=',My_Conf,' EneLower=',EneLower
2217 if (elowest.gt.etot) elowest=etot
2218 moves_acc(MoveType)=moves_acc(MoveType)+1
2219 if (MoveType.eq.1) then
2220 nbond_acc(nbond)=nbond_acc(nbond)+1
2222 ! Check against conformation repetitions.
2223 irepet=conf_comp(varia,etot)
2224 if (print_stat) then
2225 #if defined(AIX) || defined(PGI)
2226 open (istat,file=statname,position='append')
2228 open (istat,file=statname,access='append')
2231 call statprint(nacc,nfun,iretcode,etot,elowest)
2233 call var_to_geom(nvar,varia)
2235 call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,1),& !el cref(1,nstart_sup)
2236 nsup,przes,obr,non_conv)
2238 call contact(.false.,ncont,icont,co)
2239 frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
2240 write (iout,'(a,f8.3,a,f8.3)') &
2241 'RMS deviation from the reference structure:',rms,&
2242 ' % of native contacts:',frac*100,' contact order',co
2246 write (iout,*) 'Writing new conformation',nout
2248 write (istat,'(i5,16(1pe14.5))') nout,&
2249 (energia(print_order(i)),i=1,nprint_ene),&
2253 write (istat,'(i5,17(1pe14.5))') nout,&
2254 (energia(print_order(i)),i=1,nprint_ene),etot
2256 if (print_stat) close(istat)
2257 ! Print internal coordinates.
2258 if (print_int) call briefout(nout,etot)
2259 ! Accumulate the newly accepted conf in the coord1 array, if it is different
2260 ! from all confs that are already there.
2261 call compare_s1(n_thr,max_thread2,etot,varia,ii,&
2262 enetb1,coord1,rms_deform,.true.,iprint)
2263 write (iout,*) 'After compare_ss: n_thr=',n_thr
2264 if (ii.eq.1 .or. ii.eq.3) then
2265 write (iout,'(8f10.4)') &
2266 (rad2deg*coord1(i,n_thr),i=1,nvar)
2269 write (iout,*) 'Conformation from cache, not written.'
2272 if (nrepm.gt.maxrepm) then
2273 write (iout,'(a)') 'Too many conformation repetitions.'
2276 ! Store the accepted conf. and its energy.
2281 if (irepet.eq.0) call zapis(varia,etot)
2282 ! Lower the temperature, if necessary.
2293 ! Check for time limit.
2294 if (ovrtim()) WhatsUp=-1
2295 not_done = (nacc_tot.lt.maxacc) .and. (WhatsUp.eq.0) &
2304 write (iout,'(/80(1h*)/20x,a)') 'Summary run statistics:'
2305 call statprint(nacc,nfun,iretcode,etot,elowest)
2306 write (iout,'(a)') &
2307 'Statistics of multiple-bond motions. Total motions:'
2308 write (iout,'(16i5)') (nbond_move(i),i=1,Nbm)
2309 write (iout,'(a)') 'Accepted motions:'
2310 write (iout,'(16i5)') (nbond_acc(i),i=1,Nbm)
2312 write (iout,'(/80(1h*)/20x,a/80(1h*)/)') 'All iterations done.'
2313 !(maxvar) (maxvar=6*maxres)
2315 end subroutine do_mcm
2317 !-----------------------------------------------------------------------------
2319 subroutine do_mcm(i_orig)
2321 ! Monte-Carlo-with-Minimization calculations - parallel code.
2323 use control_data, only:refstr!,tag
2324 use io_base, only:intout,briefout
2325 use control, only:ovrtim,tcpu
2326 use compare, only:contact,contact_fract
2327 use minimm, only:minimize
2328 use regularize_, only:fitsq
2329 ! use contact_, only:contact
2331 ! implicit real*8 (a-h,o-z)
2332 ! include 'DIMENSIONS'
2334 ! include 'COMMON.IOUNITS'
2335 ! include 'COMMON.GEO'
2336 ! include 'COMMON.CHAIN'
2337 ! include 'COMMON.MCM'
2338 ! include 'COMMON.CONTACTS'
2339 ! include 'COMMON.CONTROL'
2340 ! include 'COMMON.VAR'
2341 ! include 'COMMON.INTERACT'
2342 ! include 'COMMON.INFO'
2343 ! include 'COMMON.CACHE'
2344 !rc include 'COMMON.DEFORM'
2345 !rc include 'COMMON.DEFORM1'
2346 !rc include 'COMMON.DEFORM2'
2347 ! include 'COMMON.MINIM'
2348 ! include 'COMMON.NAMES'
2349 logical :: accepted,over,error,lprint,not_done,similar,&
2350 enelower,non_conv,flag,finish !,ovrtim
2351 integer :: MoveType,nbond !,conf_comp
2352 real(kind=8),dimension(6*nres) :: x1,varold1,varold,varia !(maxvar) (maxvar=6*maxres)
2353 real(kind=8) :: elowest,eold
2354 real(kind=8) :: przes(3),obr(3,3)
2355 integer :: iparentx(max_threadss2)
2356 integer :: iparentx1(max_threadss2)
2357 integer :: imtasks(150),imtasks_n
2358 real(kind=8) :: energia(0:n_ene)
2361 integer :: nfun,nodenum,i_orig,i,nf,nacc,it,nout,j,kkk,is,&
2362 Kwita,iretcode,noverlap,nstart_grow,ierr,iitt,&
2363 ii_grnum_d,ii_ennum_d,ii_hesnum_d,i_grnum_d,i_ennum_d,&
2364 i_hesnum_d,i_minimiz,irepet
2365 real(kind=8) :: etot,frac,eneglobal,RandOrPert,eold1,co,&
2368 ! if(.not.allocated(varsave)) allocate(varsave(maxvar,maxsave)) !(maxvar,maxsave)
2369 print *,'Master entered DO_MCM'
2377 !---------------------------------------------------------------------------
2378 ! Initialize counters.
2379 !---------------------------------------------------------------------------
2380 ! Total number of generated confs.
2382 ! Total number of moves. In general this won`t be equal to the number of
2383 ! attempted moves, because we may want to reject some "bad" confs just by
2386 ! Total number of temperature jumps.
2388 ! Total number of shift (nbond_move(1)), spike, crankshaft, three-bond,...
2390 allocate(nbond_move(nres)) !(maxres)
2396 ! Initialize total and accepted number of moves of various kind.
2401 ! Total number of energy evaluations.
2405 ! write (iout,*) 'RanFract=',RanFract
2408 !----------------------------------------------------------------------------
2409 ! Compute and print initial energies.
2410 !----------------------------------------------------------------------------
2412 write (iout,'(/80(1h*)/a)') 'Initial energies:'
2415 call etotal(energia)
2417 call enerprint(energia)
2418 ! Request energy computation from slave processors.
2419 call geom_to_var(nvar,varia)
2420 ! write (iout,*) 'The VARIA array'
2421 ! write (iout,'(8f10.4)') (rad2deg*varia(i),i=1,nvar)
2422 call minimize(etot,varia,iretcode,nfun)
2423 call var_to_geom(nvar,varia)
2425 write (iout,*) 'etot from MINIMIZE:',etot
2426 ! write (iout,*) 'Tha VARIA array'
2427 ! write (iout,'(8f10.4)') (rad2deg*varia(i),i=1,nvar)
2430 if (print_mc .gt. 0) write (iout,'(/80(1h*)/20x,a/80(1h*))') &
2431 'Enter Monte Carlo procedure.'
2432 if (print_mc .gt. 0) write (iout,'(i5,1pe14.5)' ) i_orig,etot
2438 call zapis(varia,etot)
2440 call var_to_geom(nvar,varia)
2442 call etotal(energia)
2443 if (print_mc.gt.0) write (iout,*) 'Initial energy:',etot
2445 nacc=0 ! total # of accepted confs of the current processor.
2446 nacc_tot=0 ! total # of accepted confs of all processors.
2448 !----------------------------------------------------------------------------
2450 !----------------------------------------------------------------------------
2453 LOOP1:do while (not_done)
2455 if (print_mc.gt.0) write (iout,'(80(1h*)/20x,a,i7)') &
2456 'Beginning iteration #',it
2457 ! Initialize local counter.
2458 ntrial=0 ! # of generated non-overlapping confs.
2459 noverlap=0 ! # of overlapping confs.
2461 LOOP2:do while (.not. accepted)
2463 LOOP3:do while (imtasks_n.lt.nodenum-1.and..not.finish)
2465 if(imtasks(i).eq.0) then
2470 ! Retrieve the angles of previously accepted conformation
2474 call var_to_geom(nvar,varia)
2475 ! Rebuild the chain.
2477 ! Heat up the system, if necessary.
2479 ! If temperature cannot be further increased, stop.
2485 ! write (iout,'(a)') 'Old variables:'
2486 ! write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
2487 ! Decide whether to generate a random conformation or perturb the old one
2488 RandOrPert=ran_number(0.0D0,1.0D0)
2489 if (RandOrPert.gt.RanFract) then
2490 if (print_mc.gt.0) &
2491 write (iout,'(a)') 'Perturbation-generated conformation.'
2492 call perturb(error,lprint,MoveType,nbond,1.0D0)
2493 ! print *,'after perturb',error,finish
2494 if (error) finish = .true.
2495 if (MoveType.lt.1 .or. MoveType.gt.MaxMoveType) then
2496 write (iout,'(/a,i7,a/)') 'Error - unknown MoveType=',&
2497 MoveType,' returned from PERTURB.'
2499 write (*,'(/a,i7,a/)') 'Error - unknown MoveType=',&
2500 MoveType,' returned from PERTURB.'
2506 nstart_grow=iran_num(3,nres)
2507 if (print_mc.gt.0) &
2508 write (iout,'(2a,i3)') 'Random-generated conformation',&
2509 ' - chain regrown from residue',nstart_grow
2510 call gen_rand_conf(nstart_grow,*30)
2512 call geom_to_var(nvar,varia)
2514 ! print *,'finish=',finish
2515 if (etot-elowest.gt.overlap_cut) then
2516 if (print_mc.gt.1) write (iout,'(a,1pe14.5)') &
2517 'Overlap detected in the current conf.; energy is',etot
2518 if(iprint.gt.1.or.etot.lt.1d19) print *,&
2519 'Overlap detected in the current conf.; energy is',etot
2523 if (noverlap.gt.maxoverlap) then
2524 write (iout,*) 'Too many overlapping confs.',&
2525 ' etot, elowest, overlap_cut', etot, elowest, overlap_cut
2528 else if (.not. finish) then
2529 ! Distribute tasks to processors
2530 ! print *,'Master sending order'
2531 call MPI_SEND(12, 1, MPI_INTEGER, is, tag,&
2533 ! write (iout,*) '12: tag=',tag
2534 ! print *,'Master sent order to processor',is
2535 call MPI_SEND(it, 1, MPI_INTEGER, is, tag,&
2537 ! write (iout,*) 'it: tag=',tag
2538 call MPI_SEND(eold, 1, MPI_DOUBLE_PRECISION, is, tag,&
2540 ! write (iout,*) 'eold: tag=',tag
2541 call MPI_SEND(varia(1), nvar, MPI_DOUBLE_PRECISION, &
2544 ! write (iout,*) 'varia: tag=',tag
2545 call MPI_SEND(varold(1), nvar, MPI_DOUBLE_PRECISION, &
2548 ! write (iout,*) 'varold: tag=',tag
2555 imtasks_n=imtasks_n+1
2561 LOOP_RECV:do while(.not.flag)
2563 call MPI_IPROBE(is,tag,CG_COMM,flag,status,ierr)
2565 call MPI_RECV(iitt, 1, MPI_INTEGER, is, tag,&
2566 CG_COMM, status, ierr)
2567 call MPI_RECV(eold1, 1, MPI_DOUBLE_PRECISION, is, tag,&
2568 CG_COMM, status, ierr)
2569 call MPI_RECV(etot, 1, MPI_DOUBLE_PRECISION, is, tag,&
2570 CG_COMM, status, ierr)
2571 call MPI_RECV(varia(1), nvar, MPI_DOUBLE_PRECISION,is,tag,&
2572 CG_COMM, status, ierr)
2573 call MPI_RECV(varold1(1), nvar, MPI_DOUBLE_PRECISION, is, &
2574 tag, CG_COMM, status, ierr)
2575 call MPI_RECV(ii_grnum_d, 1, MPI_INTEGER, is, tag,&
2576 CG_COMM, status, ierr)
2577 call MPI_RECV(ii_ennum_d, 1, MPI_INTEGER, is, tag,&
2578 CG_COMM, status, ierr)
2579 call MPI_RECV(ii_hesnum_d, 1, MPI_INTEGER, is, tag,&
2580 CG_COMM, status, ierr)
2581 i_grnum_d=i_grnum_d+ii_grnum_d
2582 i_ennum_d=i_ennum_d+ii_ennum_d
2583 neneval = neneval+ii_ennum_d
2584 i_hesnum_d=i_hesnum_d+ii_hesnum_d
2585 i_minimiz=i_minimiz+1
2587 imtasks_n=imtasks_n-1
2593 if(print_mc.gt.0) write (iout,'(a,i6,a,i6,a,i6,a,1pe16.6)') &
2594 'From Worker #',is,' iitt',iitt,&
2595 ' Conformation:',ngen,' energy:',etot
2596 !--------------------------------------------------------------------------
2597 !... Do Metropolis test
2598 !--------------------------------------------------------------------------
2599 call metropolis(nvar,varia,varold1,etot,eold1,accepted,&
2601 if(iitt.ne.it.and..not.similar) then
2602 call metropolis(nvar,varia,varold,etot,eold,accepted,&
2606 if(etot.lt.eneglobal)eneglobal=etot
2607 ! if(mod(it,100).eq.0)
2608 write(iout,*)'CHUJOJEB ',neneval,eneglobal
2610 ! Write the accepted conformation.
2613 call var_to_geom(nvar,varia)
2616 call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,kkk),&
2617 nsup,przes,obr,non_conv)
2619 call contact(.false.,ncont,icont,co)
2620 frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
2621 write (iout,'(a,f8.3,a,f8.3,a,f8.3)') &
2622 'RMS deviation from the reference structure:',rms,&
2623 ' % of native contacts:',frac*100,' contact order:',co
2625 if (print_mc.gt.0) &
2626 write (iout,*) 'Writing new conformation',nout
2627 if (print_stat) then
2628 call var_to_geom(nvar,varia)
2629 #if defined(AIX) || defined(PGI)
2630 open (istat,file=statname,position='append')
2632 open (istat,file=statname,access='append')
2635 write (istat,'(i5,16(1pe14.5))') nout,&
2636 (energia(print_order(i)),i=1,nprint_ene),&
2639 write (istat,'(i5,16(1pe14.5))') nout,&
2640 (energia(print_order(i)),i=1,nprint_ene),etot
2644 ! Print internal coordinates.
2645 if (print_int) call briefout(nout,etot)
2648 if (elowest.gt.etot) elowest=etot
2649 moves_acc(MoveType)=moves_acc(MoveType)+1
2650 if (MoveType.eq.1) then
2651 nbond_acc(nbond)=nbond_acc(nbond)+1
2653 ! Check against conformation repetitions.
2654 irepet=conf_comp(varia,etot)
2655 if (nrepm.gt.maxrepm) then
2656 if (print_mc.gt.0) &
2657 write (iout,'(a)') 'Too many conformation repetitions.'
2660 ! Store the accepted conf. and its energy.
2665 if (irepet.eq.0) call zapis(varia,etot)
2666 ! Lower the temperature, if necessary.
2672 if(finish.and.imtasks_n.eq.0)exit LOOP2
2673 enddo LOOP2 ! accepted
2674 ! Check for time limit.
2675 not_done = (it.lt.max_mcm_it) .and. (nacc_tot.lt.maxacc)
2676 if(.not.not_done .or. finish) then
2677 if(imtasks_n.gt.0) then
2684 enddo LOOP1 ! not_done
2686 if (print_mc.gt.0) then
2687 write (iout,'(/80(1h*)/20x,a)') 'Summary run statistics:'
2688 call statprint(nacc,nfun,iretcode,etot,elowest)
2689 write (iout,'(a)') &
2690 'Statistics of multiple-bond motions. Total motions:'
2691 write (iout,'(16i5)') (nbond_move(i),i=1,Nbm)
2692 write (iout,'(a)') 'Accepted motions:'
2693 write (iout,'(16i5)') (nbond_acc(i),i=1,Nbm)
2695 write (iout,'(/80(1h*)/20x,a/80(1h*)/)') 'All iterations done.'
2703 call MPI_SEND(999, 1, MPI_INTEGER, is, tag,&
2707 end subroutine do_mcm
2708 !-----------------------------------------------------------------------------
2709 subroutine execute_slave(nodeinfo,iprint)
2712 use minimm, only:minimize
2714 ! implicit real*8 (a-h,o-z)
2715 ! include 'DIMENSIONS'
2717 ! include 'COMMON.TIME1'
2718 ! include 'COMMON.IOUNITS'
2719 !rc include 'COMMON.DEFORM'
2720 !rc include 'COMMON.DEFORM1'
2721 !rc include 'COMMON.DEFORM2'
2722 ! include 'COMMON.LOCAL'
2723 ! include 'COMMON.VAR'
2724 ! include 'COMMON.INFO'
2725 ! include 'COMMON.MINIM'
2726 character(len=10) :: nodeinfo
2727 real(kind=8),dimension(6*nres) :: x,x1 !(maxvar) (maxvar=6*maxres)
2728 integer :: nfun,iprint,i_switch,ierr,i_grnum_d,i_ennum_d,&
2729 i_hesnum_d,iitt,iretcode,iminrep
2730 real(kind=8) :: ener,energyx
2732 nodeinfo='chujwdupe'
2733 ! print *,'Processor:',MyID,' Entering execute_slave'
2735 ! call MPI_SEND(nodeinfo, 10, MPI_CHARACTER, 0, tag,
2738 1001 call MPI_RECV(i_switch, 1, MPI_INTEGER, 0, tag,&
2739 CG_COMM, status, ierr)
2740 ! write(iout,*)'12: tag=',tag
2741 if(iprint.ge.2)print *, MyID,' recv order ',i_switch
2742 if (i_switch.eq.12) then
2746 call MPI_RECV(iitt, 1, MPI_INTEGER, 0, tag,&
2747 CG_COMM, status, ierr)
2748 ! write(iout,*)'12: tag=',tag
2749 call MPI_RECV(ener, 1, MPI_DOUBLE_PRECISION, 0, tag,&
2750 CG_COMM, status, ierr)
2751 ! write(iout,*)'ener: tag=',tag
2752 call MPI_RECV(x(1), nvar, MPI_DOUBLE_PRECISION, 0, tag,&
2753 CG_COMM, status, ierr)
2754 ! write(iout,*)'x: tag=',tag
2755 call MPI_RECV(x1(1), nvar, MPI_DOUBLE_PRECISION, 0, tag,&
2756 CG_COMM, status, ierr)
2757 ! write(iout,*)'x1: tag=',tag
2763 ! print *,'calling minimize'
2764 call minimize(energyx,x,iretcode,nfun)
2766 write(iout,100)'minimized energy = ',energyx,&
2767 ' # funeval:',nfun,' iret ',iretcode
2768 write(*,100)'minimized energy = ',energyx,&
2769 ' # funeval:',nfun,' iret ',iretcode
2770 100 format(a20,f10.5,a12,i5,a6,i2)
2771 if(iretcode.eq.10) then
2774 write(iout,*)' ... not converged - trying again ',iminrep
2775 call minimize(energyx,x,iretcode,nfun)
2777 write(iout,*)'minimized energy = ',energyx,&
2778 ' # funeval:',nfun,' iret ',iretcode
2779 if(iretcode.ne.10)go to 812
2781 if(iretcode.eq.10) then
2783 write(iout,*)' ... not converged again - giving up'
2788 ! print *,'Sending results'
2789 call MPI_SEND(iitt, 1, MPI_INTEGER, 0, tag,&
2791 call MPI_SEND(ener, 1, MPI_DOUBLE_PRECISION, 0, tag,&
2793 call MPI_SEND(energyx, 1, MPI_DOUBLE_PRECISION, 0, tag,&
2795 call MPI_SEND(x(1), nvar, MPI_DOUBLE_PRECISION, 0, tag,&
2797 call MPI_SEND(x1(1), nvar, MPI_DOUBLE_PRECISION, 0, tag,&
2799 call MPI_SEND(i_grnum_d, 1, MPI_INTEGER, 0, tag,&
2801 call MPI_SEND(nfun, 1, MPI_INTEGER, 0, tag,&
2803 call MPI_SEND(i_hesnum_d, 1, MPI_INTEGER, 0, tag,&
2805 ! print *,'End sending'
2810 end subroutine execute_slave
2812 !-----------------------------------------------------------------------------
2813 subroutine heat(over)
2815 ! implicit real*8 (a-h,o-z)
2816 ! include 'DIMENSIONS'
2817 ! include 'COMMON.MCM'
2818 ! include 'COMMON.IOUNITS'
2820 ! Check if there`s a need to increase temperature.
2821 if (ntrial.gt.maxtrial) then
2822 if (NstepH.gt.0) then
2823 if (dabs(Tcur-TMax).lt.1.0D-7) then
2824 if (print_mc.gt.0) &
2825 write (iout,'(/80(1h*)/a,f8.3,a/80(1h*))') &
2826 'Upper limit of temperature reached. Terminating.'
2831 if (Tcur.gt.Tmax) Tcur=Tmax
2832 betbol=1.0D0/(Rbol*Tcur)
2833 if (print_mc.gt.0) &
2834 write (iout,'(/80(1h*)/a,f8.3,a,f10.5/80(1h*))') &
2835 'System heated up to ',Tcur,' K; BetBol:',betbol
2841 if (print_mc.gt.0) &
2842 write (iout,'(a)') &
2843 'Maximum number of trials in a single MCM iteration exceeded.'
2852 !-----------------------------------------------------------------------------
2855 ! implicit real*8 (a-h,o-z)
2856 ! include 'DIMENSIONS'
2857 ! include 'COMMON.MCM'
2858 ! include 'COMMON.IOUNITS'
2859 if (nstepC.gt.0 .and. dabs(Tcur-Tmin).gt.1.0D-7) then
2861 if (Tcur.lt.Tmin) Tcur=Tmin
2862 betbol=1.0D0/(Rbol*Tcur)
2863 if (print_mc.gt.0) &
2864 write (iout,'(/80(1h*)/a,f8.3,a,f10.5/80(1h*))') &
2865 'System cooled down up to ',Tcur,' K; BetBol:',betbol
2869 !-----------------------------------------------------------------------------
2870 subroutine perturb(error,lprint,MoveType,nbond,max_phi)
2873 use energy, only:nnt,nct,itype
2874 use md_calc, only:bond_move
2875 ! implicit real*8 (a-h,o-z)
2876 ! include 'DIMENSIONS'
2877 integer,parameter :: MMaxSideMove=100
2878 ! include 'COMMON.MCM'
2879 ! include 'COMMON.CHAIN'
2880 ! include 'COMMON.INTERACT'
2881 ! include 'COMMON.VAR'
2882 ! include 'COMMON.GEO'
2883 ! include 'COMMON.NAMES'
2884 ! include 'COMMON.IOUNITS'
2885 !rc include 'COMMON.DEFORM1'
2886 logical :: error,lprint,fail
2887 integer :: MoveType,nbond,end_select,ind_side(MMaxSideMove)
2888 real(kind=8) :: max_phi
2889 real(kind=8) :: psi!,gen_psi
2890 !el external iran_num
2891 !el integer iran_num
2894 !!! local variables - el
2895 integer :: itrial,iiwin,iwindow,isctry,i,icount,j,nstart,&
2896 nside_move,inds,indx,ii,iti
2897 real(kind=8) :: bond_prob,theta_new
2902 ! Perturb the conformation according to a randomly selected move.
2903 call SelectMove(MoveType)
2904 ! write (iout,*) 'MoveType=',MoveType
2906 goto (100,200,300,400,500) MoveType
2907 !------------------------------------------------------------------------------
2908 ! Backbone N-bond move.
2909 ! Select the number of bonds (length of the segment to perturb).
2911 if (itrial.gt.1000) then
2912 write (iout,'(a)') 'Too many attempts at multiple-bond move.'
2916 bond_prob=ran_number(0.0D0,sumpro_bond(nbm))
2917 ! print *,'sumpro_bond(nbm)=',sumpro_bond(nbm),
2918 ! & ' Bond_prob=',Bond_Prob
2920 ! print *,i,Bond_Prob,sumpro_bond(i),sumpro_bond(i+1)
2921 if (bond_prob.ge.sumpro_bond(i) .and. &
2922 bond_prob.le.sumpro_bond(i+1)) then
2927 write (iout,'(2a)') 'In PERTURB: Error - number of bonds',&
2928 ' to move out of range.'
2932 if (nwindow.gt.0) then
2933 ! Select the first residue to perturb
2934 iwindow=iran_num(1,nwindow)
2935 print *,'iwindow=',iwindow
2937 do while (winlen(iwindow).lt.nbond)
2938 iwindow=iran_num(1,nwindow)
2940 if (iiwin.gt.1000) then
2941 write (iout,'(a)') 'Cannot select moveable residues.'
2946 nstart=iran_num(winstart(iwindow),winend(iwindow))
2948 nstart = iran_num(koniecl+2,nres-nbond-koniecl)
2949 !d print *,'nres=',nres,' nbond=',nbond,' koniecl=',koniecl,
2950 !d & ' nstart=',nstart
2953 if (psi.eq.0.0) then
2957 if (print_mc.gt.1) write (iout,'(a,i4,a,i4,a,f8.3)') &
2958 'PERTURB: nbond=',nbond,' nstart=',nstart,' psi=',psi*rad2deg
2959 !d print *,'nstart=',nstart
2960 call bond_move(nbond,nstart,psi,.false.,error)
2962 write (iout,'(2a)') &
2963 'Could not define reference system in bond_move, ',&
2964 'choosing ahother segment.'
2968 nbond_move(nbond)=nbond_move(nbond)+1
2972 !------------------------------------------------------------------------------
2973 ! Backbone endmove. Perturb a SINGLE angle of a residue close to the end of
2977 ! end_select=iran_num(1,2*koniecl)
2978 ! if (end_select.gt.koniecl) then
2979 ! end_select=nphi-(end_select-koniecl)
2981 ! end_select=koniecl+3
2983 ! if (nwindow.gt.0) then
2984 ! iwin=iran_num(1,nwindow)
2985 ! i1=max0(4,winstart(iwin))
2986 ! i2=min0(winend(imin)+2,nres)
2987 ! end_select=iran_num(i1,i2)
2989 ! iselect = iran_num(1,nmov_var)
2992 ! if (isearch_tab(i).eq.1) jj = jj+1
2993 ! if (jj.eq.iselect) then
2999 end_select = iran_num(4,nres)
3000 psi=max_phi*gen_psi()
3001 if (psi.eq.0.0D0) then
3005 phi(end_select)=pinorm(phi(end_select)+psi)
3006 if (print_mc.gt.1) write (iout,'(a,i4,a,f8.3,a,f8.3)') &
3007 'End angle',end_select,' moved by ',psi*rad2deg,' new angle:',&
3008 phi(end_select)*rad2deg
3009 ! if (end_select.gt.3)
3010 ! & theta(end_select-1)=gen_theta(itype(end_select-2),
3011 ! & phi(end_select-1),phi(end_select))
3012 ! if (end_select.lt.nres)
3013 ! & theta(end_select)=gen_theta(itype(end_select-1),
3014 ! & phi(end_select),phi(end_select+1))
3015 !d print *,'nres=',nres,' end_select=',end_select
3016 !d print *,'theta',end_select-1,theta(end_select-1)
3017 !d print *,'theta',end_select,theta(end_select)
3022 !------------------------------------------------------------------------------
3024 ! Select the number of SCs to perturb.
3026 301 nside_move=iran_num(1,MaxSideMove)
3027 ! print *,'nside_move=',nside_move,' MaxSideMove',MaxSideMove
3028 ! Select the indices.
3031 111 inds=iran_num(nnt,nct)
3033 if (icount.gt.1000) then
3034 write (iout,'(a)')'Error - cannot select side chains to move.'
3038 if (itype(inds).eq.10) goto 111
3040 if (inds.eq.ind_side(j)) goto 111
3043 if (inds.lt.ind_side(j)) then
3049 112 do j=i,indx+1,-1
3050 ind_side(j)=ind_side(j-1)
3052 113 ind_side(indx)=inds
3054 ! Carry out perturbation.
3058 call gen_side(iti,theta(ii+1),alph(ii),omeg(ii),fail)
3061 if (isctry.gt.1000) then
3062 write (iout,'(a)') 'Too many errors in SC generation.'
3068 if (print_mc.gt.1) write (iout,'(2a,i4,a,2f8.3)') &
3069 'Side chain ',restyp(iti),ii,' moved to ',&
3070 alph(ii)*rad2deg,omeg(ii)*rad2deg
3075 !------------------------------------------------------------------------------
3077 400 end_select=iran_num(3,nres)
3078 theta_new=gen_theta(itype(end_select),phi(end_select),&
3080 if (print_mc.gt.1) write (iout,'(a,i3,a,f8.3,a,f8.3)') &
3081 'Theta ',end_select,' moved from',theta(end_select)*rad2deg,&
3082 ' to ',theta_new*rad2deg
3083 theta(end_select)=theta_new
3087 !------------------------------------------------------------------------------
3088 ! Error returned from SelectMove.
3091 end subroutine perturb
3092 !-----------------------------------------------------------------------------
3093 subroutine SelectMove(MoveType)
3095 ! implicit real*8 (a-h,o-z)
3096 ! include 'DIMENSIONS'
3097 ! include 'COMMON.MCM'
3098 ! include 'COMMON.IOUNITS'
3100 !!! local variables - el
3101 integer :: i,MoveType
3102 real(kind=8) :: what_move
3104 what_move=ran_number(0.0D0,sumpro_type(MaxMoveType))
3106 if (what_move.ge.sumpro_type(i-1).and. &
3107 what_move.lt.sumpro_type(i)) then
3112 write (iout,'(a)') &
3113 'Fatal error in SelectMoveType: cannot select move.'
3114 MoveType=MaxMoveType+1
3116 end subroutine SelectMove
3117 !-----------------------------------------------------------------------------
3118 real(kind=8) function gen_psi()
3120 use geometry_data, only: angmin,pi
3123 real(kind=8) :: x !,ran_number
3124 ! include 'COMMON.IOUNITS'
3125 ! include 'COMMON.GEO'
3128 x=ran_number(-pi,pi)
3129 if (dabs(x).gt.angmin) then
3134 write (iout,'(a)')'From Gen_Psi: Cannot generate angle increment.'
3137 end function gen_psi
3138 !-----------------------------------------------------------------------------
3139 subroutine metropolis(n,xcur,xold,ecur,eold,accepted,similar,enelower)
3141 ! implicit real*8 (a-h,o-z)
3142 ! include 'DIMENSIONS'
3143 ! include 'COMMON.MCM'
3144 ! include 'COMMON.IOUNITS'
3145 ! include 'COMMON.VAR'
3146 ! include 'COMMON.GEO'
3147 !rc include 'COMMON.DEFORM'
3149 real(kind=8) :: ecur,eold,xx,bol !,ran_number
3150 real(kind=8),dimension(n) :: xcur,xold
3151 real(kind=8) :: ecut1 ,ecut2 ,tola
3152 logical :: accepted,similar,not_done,enelower
3155 !!! local variables - el
3156 real(kind=8) :: xxh,difene,reldife
3158 data ecut1 /-1.0D-5/,ecut2 /5.0D-3/,tola/5.0D0/
3162 ! Set lprn=.true. for debugging.
3165 !el write(iout,*)'enedif',enedif,' ecut1',ecut1,' ecut2',ecut2
3166 write(iout,*)' ecut1',ecut1,' ecut2',ecut2
3170 ! Check if the conformation is similar.
3172 reldife=difene/dmax1(dabs(eold),dabs(ecur),1.0D0)
3174 write (iout,*) 'Metropolis'
3175 write(iout,*)'ecur,eold,difene,reldife',ecur,eold,difene,reldife
3177 ! If energy went down remarkably, we accept the new conformation
3179 !jp if (reldife.lt.ecut1) then
3180 if (difene.lt.ecut1) then
3183 if (lprn) write (iout,'(a)') &
3184 'Conformation accepted, because energy has lowered remarkably.'
3185 ! elseif (reldife.lt.ecut2 .and. dif_ang(nphi,xcur,xold).lt.tola)
3186 !jp elseif (reldife.lt.ecut2)
3187 elseif (difene.lt.ecut2) &
3189 ! Reject the conf. if energy has changed insignificantly and there is not
3190 ! much change in conformation.
3192 write (iout,'(2a)') 'Conformation rejected, because it is',&
3193 ' similar to the preceding one.'
3197 ! Else carry out Metropolis test.
3199 xx=ran_number(0.0D0,1.0D0)
3202 write (iout,*) 'betbol=',betbol,' difene=',difene,' xxh=',xxh
3203 if (xxh.gt.50.0D0) then
3208 if (lprn) write (iout,*) 'bol=',bol,' xx=',xx
3211 if (lprn) write (iout,'(a)') &
3212 'Conformation accepted, because it passed Metropolis test.'
3215 if (lprn) write (iout,'(a)') &
3216 'Conformation rejected, because it did not pass Metropolis test.'
3225 end subroutine metropolis
3226 !-----------------------------------------------------------------------------
3227 integer function conf_comp(x,ene)
3229 use geometry_data, only: nphi
3230 ! implicit real*8 (a-h,o-z)
3231 ! include 'DIMENSIONS'
3232 ! include 'COMMON.MCM'
3233 ! include 'COMMON.VAR'
3234 ! include 'COMMON.IOUNITS'
3235 ! include 'COMMON.GEO'
3236 real(kind=8) :: etol, angtol
3237 real(kind=8),dimension(6*nres) :: x !(maxvar) (maxvar=6*maxres)
3238 real(kind=8) :: difa !dif_ang,
3240 !!! local variables - el
3244 data etol /0.1D0/, angtol /20.0D0/
3246 ! write (iout,*) 'ii=',ii,'ene=',ene,esave(ii),dabs(ene-esave(ii))
3247 if (dabs(ene-esave(ii)).lt.etol) then
3248 difa=dif_ang(nphi,x,varsave(1,ii))
3250 ! write(iout,'(i3,3f8.3)')i,rad2deg*x(i),
3251 ! & rad2deg*varsave(i,ii)
3253 ! write(iout,*) 'ii=',ii,' difa=',difa,' angtol=',angtol
3254 if (difa.le.angtol) then
3255 if (print_mc.gt.0) then
3256 write (iout,'(a,i5,2(a,1pe15.4))') &
3257 'Current conformation matches #',ii,&
3258 ' in the store array ene=',ene,' esave=',esave(ii)
3259 ! write (*,'(a,i5,a)') 'Current conformation matches #',ii,
3260 ! & ' in the store array.'
3261 endif ! print_mc.gt.0
3262 if (print_mc.gt.1) then
3264 write(iout,'(i3,3f8.3)')i,rad2deg*x(i),&
3265 rad2deg*varsave(i,ii)
3267 endif ! print_mc.gt.1
3276 end function conf_comp
3277 !-----------------------------------------------------------------------------
3278 real(kind=8) function dif_ang(n,x,y)
3280 use geometry_data, only: dwapi
3283 real(kind=8),dimension(n) :: x,y
3284 real(kind=8) :: w,wa,dif,difa
3285 !el real(kind=8) :: pinorm
3286 ! include 'COMMON.GEO'
3290 dif=dabs(pinorm(y(i)-x(i)))
3291 if (dabs(dif-dwapi).lt.dif) dif=dabs(dif-dwapi)
3292 w=1.0D0-(2.0D0*(i-1)/(n-1)-1.0D0)**2+1.0D0/n
3296 dif_ang=rad2deg*dsqrt(difa/wa)
3298 end function dif_ang
3299 !-----------------------------------------------------------------------------
3300 subroutine add2cache(n1,n2,ncache,nvar,SourceID,CachSrc,ecur,xcur,ecache,xcache)
3303 ! include 'COMMON.GEO'
3304 ! include 'COMMON.IOUNITS'
3305 integer :: n1,n2,ncache,nvar,SourceID,CachSrc(n2)
3307 real(kind=8) :: ecur,xcur(nvar),ecache(n2),xcache(n1,n2)
3308 !d write (iout,*) 'Enter ADD2CACHE ncache=',ncache ,' ecur',ecur
3309 !d write (iout,'(10f8.3)') (rad2deg*xcur(i),i=1,nvar)
3310 !d write (iout,*) 'Old CACHE array:'
3312 !d write (iout,*) 'i=',i,' ecache=',ecache(i),' CachSrc',CachSrc(i)
3313 !d write (iout,'(10f8.3)') (rad2deg*xcache(j,i),j=1,nvar)
3317 do while (i.gt.0 .and. ecur.lt.ecache(i))
3321 !d write (iout,*) 'i=',i,' ncache=',ncache
3322 if (ncache.eq.n2) then
3323 write (iout,*) 'Cache dimension exceeded',ncache,n2
3324 write (iout,*) 'Highest-energy conformation will be removed.'
3328 ecache(ii+1)=ecache(ii)
3329 CachSrc(ii+1)=CachSrc(ii)
3331 xcache(j,ii+1)=xcache(j,ii)
3340 !d write (iout,*) 'New CACHE array:'
3342 !d write (iout,*) 'i=',i,' ecache=',ecache(i),' CachSrc',CachSrc(i)
3343 !d write (iout,'(10f8.3)') (rad2deg*xcache(j,i),j=1,nvar)
3346 end subroutine add2cache
3347 !-----------------------------------------------------------------------------
3348 subroutine rm_from_cache(i,n1,n2,ncache,nvar,CachSrc,ecache,xcache)
3351 ! include 'COMMON.GEO'
3352 ! include 'COMMON.IOUNITS'
3353 integer :: n1,n2,ncache,nvar,CachSrc(n2)
3355 real(kind=8) :: ecache(n2),xcache(n1,n2)
3357 !d write (iout,*) 'Enter RM_FROM_CACHE'
3358 !d write (iout,*) 'Old CACHE array:'
3360 !d write (iout,*)'i=',ii,' ecache=',ecache(ii),' CachSrc',CachSrc(ii)
3361 !d write (iout,'(10f8.3)') (rad2deg*xcache(j,ii),j=1,nvar)
3365 ecache(ii-1)=ecache(ii)
3366 CachSrc(ii-1)=CachSrc(ii)
3368 xcache(j,ii-1)=xcache(j,ii)
3372 !d write (iout,*) 'New CACHE array:'
3374 !d write (iout,*) 'i=',i,' ecache=',ecache(i),' CachSrc',CachSrc(i)
3375 !d write (iout,'(10f8.3)') (rad2deg*xcache(j,i),j=1,nvar)
3378 end subroutine rm_from_cache
3379 !-----------------------------------------------------------------------------
3381 !-----------------------------------------------------------------------------
3382 subroutine statprint(it,nfun,iretcode,etot,elowest)
3384 use control_data, only: MaxMoveType,minim
3385 use control, only: tcpu
3387 ! implicit real*8 (a-h,o-z)
3388 ! include 'DIMENSIONS'
3389 ! include 'COMMON.IOUNITS'
3390 ! include 'COMMON.CONTROL'
3391 ! include 'COMMON.MCM'
3393 integer :: it,nfun,iretcode,i
3394 real(kind=8) :: etot,elowest,fr_mov_i
3398 '(80(1h*)/a,i5,a,1pe14.5,a,1pe14.5/a,i3,a,i10,a,i5,a,i5)') &
3399 'Finished iteration #',it,' energy is',etot,&
3400 ' lowest energy:',elowest,&
3401 'SUMSL return code:',iretcode,&
3402 ' # of energy evaluations:',neneval,&
3403 '# of temperature jumps:',ntherm,&
3404 ' # of minima repetitions:',nrepm
3406 write (iout,'(80(1h*)/a,i8,a,1pe14.5,a,1pe14.5)') &
3407 'Finished iteration #',it,' energy is',etot,&
3408 ' lowest energy:',elowest
3410 write (iout,'(/4a)') &
3411 'Kind of move ',' total',' accepted',&
3413 write (iout,'(58(1h-))')
3415 if (moves(i).eq.0) then
3418 fr_mov_i=dfloat(moves_acc(i))/dfloat(moves(i))
3420 write(iout,'(a,2i15,f10.5)')MovTypID(i),moves(i),moves_acc(i),&
3423 write (iout,'(a,2i15,f10.5)') 'total ',nmove,nacc_tot,&
3424 dfloat(nacc_tot)/dfloat(nmove)
3425 write (iout,'(58(1h-))')
3426 write (iout,'(a,1pe12.4)') 'Elapsed time:',tcpu()
3428 end subroutine statprint
3429 !-----------------------------------------------------------------------------
3430 subroutine zapis(varia,etot)
3432 use geometry_data, only: nres,rad2deg,nvar
3435 ! implicit real*8 (a-h,o-z)
3436 ! include 'DIMENSIONS'
3438 use MPI_data !include 'COMMON.INFO'
3441 ! include 'COMMON.GEO'
3442 ! include 'COMMON.VAR'
3443 ! include 'COMMON.MCM'
3444 ! include 'COMMON.IOUNITS'
3445 integer,dimension(nsave) :: itemp
3446 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
3449 integer :: j,i,maxvar
3450 real(kind=8) :: etot
3452 !el allocate(esave(nsave)) !(maxsave)
3457 write (iout,'(a,i5,a,i5)') 'Enter ZAPIS NSave=',Nsave,&
3459 write (iout,'(a)') 'Current energy and conformation:'
3460 write (iout,'(1pe14.5)') etot
3461 write (iout,'(10f8.3)') (rad2deg*varia(i),i=1,nvar)
3463 ! Shift the contents of the esave and varsave arrays if filled up.
3464 call add2cache(6*nres,nsave,nsave,nvar,MyID,itemp,&
3465 etot,varia,esave,varsave)
3467 write (iout,'(a)') 'Energies and the VarSave array.'
3469 write (iout,'(i5,1pe14.5)') i,esave(i)
3470 write (iout,'(10f8.3)') (rad2deg*varsave(j,i),j=1,nvar)
3474 end subroutine zapis
3475 !-----------------------------------------------------------------------------
3476 !-----------------------------------------------------------------------------
3477 subroutine alloc_MCM_arrays
3479 use energy_data, only: max_ene
3483 allocate(entropy(-max_ene-4:max_ene)) !(-max_ene-4:max_ene)
3484 allocate(nhist(-max_ene:max_ene)) !(-max_ene:max_ene)
3485 allocate(nminima(maxsave)) !(maxsave)
3487 allocate(xpool(6*nres,max_pool)) !(maxvar,max_pool)(maxvar=6*maxres)
3488 allocate(epool(max_pool)) !(max_pool)
3491 if(.not.allocated(nsave_part)) allocate(nsave_part(nctasks)) !(max_cg_procs)
3494 ! real(kind=8),dimension(:),allocatable :: sumpro_type !(0:MaxMoveType)
3495 allocate(sumpro_bond(0:nres)) !(0:maxres)
3496 allocate(nbond_move(nres),nbond_acc(nres)) !(maxres)
3497 allocate(moves(-1:MaxMoveType+1),moves_acc(-1:MaxMoveType+1)) !(-1:MaxMoveType+1)
3498 ! common /accept_stats/
3499 ! allocate(nacc_part !(0:MaxProcs) !el nie uzywane???
3500 ! common /windows/ in io: mcmread
3501 ! allocate(winstart,winend,winlen !(maxres)
3503 !el allocate(MovTypID(-1:MaxMoveType+1)) !(-1:MaxMoveType+1)
3506 allocate(varsave(nres*6,maxsave)) !(maxvar,maxsave)(maxvar=6*maxres)
3507 allocate(esave(maxsave)) !(maxsave)
3508 allocate(Origin(maxsave)) !(maxsave)
3511 end subroutine alloc_MCM_arrays
3512 !-----------------------------------------------------------------------------
3513 !-----------------------------------------------------------------------------