2 !-----------------------------------------------------------------------------
6 use geometry_data, only: nres,nvar,rad2deg
7 use random, only: iran_num,ran_number
14 !-----------------------------------------------------------------------------
15 ! Max. number of conformations in Master's cache array
16 integer,parameter :: max_cache=10
17 !-----------------------------------------------------------------------------
18 ! Number of threads in deformation
19 integer,parameter :: max_thread=4, max_thread2=2*max_thread
20 !-----------------------------------------------------------------------------
21 ! Number of structures to compare at t=0
22 integer,parameter :: max_threadss=8,max_threadss2=2*max_threadss
23 !-----------------------------------------------------------------------------
24 ! Max. number of conformations in the pool
25 integer,parameter :: max_pool=10
26 !-----------------------------------------------------------------------------
27 !-----------------------------------------------------------------------------
31 !-----------------------------------------------------------------------------
34 ! real(kind=8) :: emin,emax
35 real(kind=8),dimension(:),allocatable :: entropy !(-max_ene-4:max_ene)
36 real(kind=8),dimension(:),allocatable :: nhist !(-max_ene:max_ene)
37 real(kind=8),dimension(:),allocatable :: nminima !(maxsave)
40 integer :: indminn,indmaxx
43 ! real(kind=8) :: pool_fraction
44 real(kind=8),dimension(:,:),allocatable :: xpool !(maxvar,max_pool)
45 real(kind=8),dimension(:),allocatable :: epool !(max_pool)
46 ! common /mce_counters/
47 !------------------------------------------------------------------------------
48 !... Following COMMON block contains variables controlling motion.
49 !------------------------------------------------------------------------------
51 ! real(kind=8),dimension(0:MaxMoveType) :: sumpro_type !(0:MaxMoveType)
52 real(kind=8),dimension(:),allocatable :: sumpro_bond !(0:maxres)
53 integer :: koniecl,Nbm,MaxSideMove!,nmove
54 integer,dimension(:),allocatable :: nbond_move,nbond_acc !(maxres)
55 ! integer,dimension(-1:MaxMoveType+1) :: moves,moves_acc !(-1:MaxMoveType+1)
56 ! common /accept_stats/
58 integer,dimension(:),allocatable :: nacc_part !(0:MaxProcs) !el nie uzywane???
61 !------------------------------------------------------------------------------
62 !... koniecl - the number of bonds to be considered "end bonds" subjected to
64 !... Nbm - The maximum length of N-bond segment to be moved;
65 !... MaxSideMove - maximum number of side chains subjected to local moves
67 !... nmove - the current number of attempted moves;
68 !... nbond_move(*) array that stores the total numbers of 2-bond,3-bond,...
70 !... nendmove - number of endmoves;
71 !... nbackmove - number of backbone moves;
72 !... nsidemove - number of local side chain moves;
73 !... sumpro_type(*) - array that stores the lower and upper boundary of the
74 !... random-number range that determines the type of move
75 !... (N-bond, backbone or side chain);
76 !... sumpro_bond(*) - array that stores the probabilities to perform bond
77 !... moves of consecutive segment length.
78 !... winstart(*) - the starting position of the perturbation window;
79 !... winend(*) - the end position of the perturbation window;
80 !... winlen(*) - length of the perturbation window;
81 !... nwindow - the number of perturbation windows (0 - entire chain).
82 !-----------------------------------------------------------------------------
85 !-----------------------------------------------------------------------------
87 !-----------------------------------------------------------------------------
89 !-----------------------------------------------------------------------------
90 subroutine compare_s1(n_thr,num_thread_save,energyx,x,icomp,enetbss,&
91 coordss,rms_d,modif,iprint)
93 ! This subroutine compares the new conformation, whose variables are in X
94 ! with the previously accumulated conformations whose energies and variables
95 ! are stored in ENETBSS and COORDSS, respectively. The meaning of other
96 ! variables is as follows:
98 ! N_THR - on input the previous # of accumulated confs, on output the current
99 ! # of accumulated confs.
100 ! N_REPEAT - an array that indicates how many times the structure has already
101 ! been used to start the reversed-reversing procedure. Addition of
102 ! a new structure replacement of a structure with a similar, but
103 ! lower-energy structure resets the respective entry in N_REPEAT to zero
105 ! ENERGYX,X - the energy and variables of the new conformations.
106 ! ICOMP - comparison result:
107 ! 0 - the new structure is similar to one of the previous ones and does
108 ! not have a remarkably lower energy and is therefore rejected;
109 ! 1 - the new structure is different and is added to the set, because
110 ! there is still room in the COORDSS and ENETBSS arrays;
111 ! 2 - the new structure is different, but higher in energy than any
112 ! previous one and is therefore rejected
113 ! 3 - there is no more room in the COORDSS and ENETBSS arrays, but
114 ! the new structure is lower in energy than at least the highest-
115 ! energy previous structure and therefore replaces it.
116 ! 9 - the new structure is similar to a number of previous structures,
117 ! but has a remarkably lower energy than any of them; therefore
118 ! replaces all these structures;
119 ! MODIF - a logical variable that shows whether to include the new structure
120 ! in the set of accumulated structures
122 ! implicit real*8 (a-h,o-z)
124 use regularize_, only:fitsq
125 ! include 'DIMENSIONS'
126 ! include 'COMMON.GEO'
127 ! include 'COMMON.VAR'
128 !rc include 'COMMON.DEFORM'
129 ! include 'COMMON.IOUNITS'
131 !el use geometry_data !include 'COMMON.CHAIN'
134 real(kind=8),dimension(6*nres) :: x,x1 !(maxvar) (maxvar=6*maxres)
135 real(kind=8) :: przes(3),obrot(3,3)
136 integer :: list(max_thread)
137 logical :: non_conv,modif
138 real(kind=8) :: enetbss(max_threadss)
139 real(kind=8) :: coordss(6*nres,max_threadss)
141 !!! local variables - el
142 integer :: n_thr,num_thread_save,icomp,minimize_s_flag,iprint
143 real(kind=8) :: energyx,energyy,rms_d
144 integer :: nlist,k,kk,j,i,iresult
145 real(kind=8) :: enex_jp,roznica
149 call var_to_geom(nvar,x)
157 ! write(iout,*)'*ene=',energyx
164 if (iprint.gt.3) then
165 write (iout,*) 'Compare_ss, i=',i
166 write (iout,*) 'New structure Energy:',energyx
167 write (iout,'(10f8.3)') (rad2deg*x(k),k=1,nvar)
168 write (iout,*) 'Template structure Energy:',enetbss(i)
169 write (iout,'(10f8.3)') (rad2deg*x1(k),k=1,nvar)
173 call var_to_geom(nvar,x1)
175 !d write(iout,*)'C and CREF'
176 !d write(iout,'(i5,3f10.5,5x,3f10.5)')(k,(c(j,k),j=1,3),
177 !d & (cref(j,k),j=1,3),k=1,nres)
178 call fitsq(roznica,c(1,1),cref(1,1,1),nres,przes,obrot,non_conv)
180 print *,'Problems in FITSQ!!!'
182 print '(10f8.3)',(x(k),k=1,nvar)
184 print '(10f8.3)',(x1(k),k=1,nvar)
186 print '(i5,3f10.5,5x,3f10.5)',(k,(c(j,k),j=1,3),&
187 (cref(j,k,1),j=1,3),k=1,nres)
189 roznica=dsqrt(dabs(roznica))
191 if (roznica.lt.rms_d) iresult = 0
194 !el call cmprs(x,x1,roznica,energyx,energyy,iresult)
196 if (iprint.gt.1) write(iout,'(i5,f10.6,$)') i,roznica
197 ! print '(i5,f8.3)',i,roznica
198 if(iresult.eq.0) then
201 if (iprint.gt.1) write(iout,*)
202 if(energyx.ge.enetbss(i)) then
204 write(iout,*)'s*>> structure rejected - same as nr ',i, &
211 if(energyx.lt.enetbss(i).and.enex_jp.lt.enetbss(i))then
216 if (iprint.gt.1) write(iout,*)
220 write(iout,'(a,i3,$)')'s*>> structure accepted1 - repl nr ',&
224 write(iout,'(a,i3)') &
225 's*>> structure accepted1 - would repl nr ',list(1)
228 if (.not. modif) goto 1106
235 if (iprint.gt.1) write(iout,'(i3,$)')list(j)
236 do kk=list(j)+1,nlist
237 enetbss(kk-1)=enetbss(kk)
239 coordss(i,kk-1)=coordss(i,kk)
243 if (iprint.gt.1) write(iout,*)
246 if(n_thr.lt.num_thread_save) then
250 write(iout,*)'s*>> structure accepted - add with nr ',n_thr+1
253 write(iout,*)'s*>> structure accepted - would add with nr ',&
258 enetbss(n_thr)=energyx
260 coordss(i,n_thr)=x(i)
265 write(iout,*)'s*>> structure rejected - too high energy'
272 write(iout,*)'s*>> structure accepted - repl nr ',j
275 write(iout,*)'s*>> structure accepted - would repl nr ',j
286 end subroutine compare_s1
287 !-----------------------------------------------------------------------------
289 !-----------------------------------------------------------------------------
294 use MPI_data, only:WhatsUp,MyID
295 use compare_data, only: ener
296 use control_data, only: minim,refstr
298 use regularize_, only:fitsq
299 use control, only: tcpu,ovrtim
301 use minimm, only:minimize
302 ! Does modified entropic sampling in the space of minima.
303 ! implicit real*8 (a-h,o-z)
304 ! include 'DIMENSIONS'
305 ! include 'COMMON.IOUNITS'
307 use MPI_data !include 'COMMON.INFO'
309 ! include 'COMMON.GEO'
310 ! include 'COMMON.CHAIN'
311 ! include 'COMMON.MCM'
312 ! include 'COMMON.MCE'
313 ! include 'COMMON.CONTACTS'
314 ! include 'COMMON.CONTROL'
315 ! include 'COMMON.VAR'
316 ! include 'COMMON.INTERACT'
317 ! include 'COMMON.THREAD'
318 ! include 'COMMON.NAMES'
319 logical :: accepted,not_done,over,error,lprint !,ovrtim
320 integer :: MoveType,nbond
321 ! integer :: conf_comp
322 real(kind=8) :: RandOrPert
323 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
324 real(kind=8) :: elowest,ehighest,eold
325 real(kind=8) :: przes(3),obr(3,3)
326 real(kind=8),dimension(6*nres) :: varold !(maxvar) (maxvar=6*maxres)
328 real(kind=8),dimension(0:n_ene) :: energia,energia_ave
330 !!! local variables -el
331 integer :: i,ii,kkk,it,j,nacc,nfun,ijunk,indmin,indmax,&
332 ISWEEP,Kwita,iretcode,indeold,iene,noverlap,&
333 irep,nstart_grow,inde
334 real(kind=8) :: facee,conste,ejunk,etot,rms,co,frac,&
335 deix,dent,sold,scur,runtime
338 ! if(.not.allocated(varsave)) allocate(varsave(maxvar,maxsave)) !(maxvar,maxsave)
339 !d write (iout,*) 'print_mc=',print_mc
342 !---------------------------------------------------------------------------
343 ! Initialize counters.
344 !---------------------------------------------------------------------------
345 ! Total number of generated confs.
347 ! Total number of moves. In general this won't be equal to the number of
348 ! attempted moves, because we may want to reject some "bad" confs just by
351 ! Total number of shift (nbond_move(1)), spike, crankshaft, three-bond,...
353 !el allocate(nbond_move(nres)) !(maxres)
358 ! Initialize total and accepted number of moves of various kind.
363 ! Total number of energy evaluations.
369 facee=1.0D0/(maxacc*delte)
371 ! Read entropy from previous simulations.
373 read (ientin,*) indminn,indmaxx,emin,emax
374 print *,'indminn=',indminn,' indmaxx=',indmaxx,' emin=',emin,&
376 do i=-max_ene,max_ene
377 entropy(i)=(emin+i*delte)*betbol
379 read (ientin,*) (ijunk,ejunk,entropy(i),i=indminn,indmaxx)
382 write (iout,*) 'indminn=',indminn,' indmaxx=',indmaxx,&
383 ' emin=',emin,' emax=',emax
384 write (iout,'(/a)') 'Initial entropy'
386 write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
389 ! Read the pool of conformations
391 !----------------------------------------------------------------------------
392 ! Entropy-sampling simulations with continually updated entropy
393 ! Loop thru simulations
394 !----------------------------------------------------------------------------
396 !----------------------------------------------------------------------------
397 ! Take a conformation from the pool
398 !----------------------------------------------------------------------------
404 write (iout,*) 'Took conformation',ii,' from the pool energy=',&
406 call var_to_geom(nvar,varia)
407 ! Print internal coordinates of the initial conformation
410 call gen_rand_conf(1,*20)
412 !----------------------------------------------------------------------------
413 ! Compute and print initial energies.
414 !----------------------------------------------------------------------------
417 allocate(nsave_part(nctasks))
418 if (MyID.eq.MasterID) then
426 write (iout,'(/80(1h*)/a,i2/80(1h*)/)') 'MCE iteration #',isweep
427 write (iout,'(/80(1h*)/a)') 'Initial energies:'
431 call enerprint(energia)
432 ! Minimize the energy of the first conformation.
434 call geom_to_var(nvar,varia)
435 call minimize(etot,varia,iretcode,nfun)
438 write (iout,'(/80(1h*)/a/80(1h*))') &
439 'Results of the first energy minimization:'
440 call enerprint(energia)
444 call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,kkk),&
448 call contact(.false.,ncont,icont,co)
449 frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
450 write (iout,'(a,f8.3,a,f8.3,a,f8.3)') &
451 'RMS deviation from the reference structure:',rms,&
452 ' % of native contacts:',frac*100,' contact order:',co
453 write (istat,'(i5,11(1pe14.5))') 0,&
454 (energia(print_order(i)),i=1,nprint_ene),etot,rms,frac,co
456 write (istat,'(i5,9(1pe14.5))') 0,&
457 (energia(print_order(i)),i=1,nprint_ene),etot
460 neneval=neneval+nfun+1
461 if (.not. ent_read) then
462 ! Initialize the entropy array
463 do i=-max_ene,max_ene
465 ! Uncomment the line below for actual entropic sampling (start with uniform
466 ! energy distribution).
468 ! Uncomment the line below for multicanonical sampling (start with Boltzmann
470 entropy(i)=(emin+i*delte)*betbol
474 write (iout,'(/a)') 'Initial entropy'
476 write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
480 call recv_stop_sig(Kwita)
481 if (whatsup.eq.1) then
482 call send_stop_sig(-2)
484 else if (whatsup.le.-2) then
486 else if (whatsup.eq.2) then
492 not_done = (iretcode.ne.11)
494 write (iout,'(/80(1h*)/20x,a/80(1h*))') &
495 'Enter Monte Carlo procedure.'
497 call briefout(0,etot)
502 indeold=(eold-emin)/delte
503 deix=eold-(emin+indeold*delte)
504 dent=entropy(indeold+1)-entropy(indeold)
505 !d write (iout,*) 'indeold=',indeold,' deix=',deix,' dent=',dent
506 !d write (*,*) 'Processor',MyID,' indeold=',indeold,' deix=',deix,
508 sold=entropy(indeold)+(dent/delte)*deix
510 write (iout,*) 'eold=',eold,' sold=',sold,' elowest=',etot
511 write (*,*) 'Processor',MyID,' eold=',eold,' sold=',sold,&
513 if (minim) call zapis(varia,etot)
515 ! NACC is the counter for the accepted conformations of a given processor
517 ! NACC_TOT counts the total number of accepted conformations
520 if (MyID.eq.MasterID) then
521 call receive_MCM_info
523 call send_MCM_info(2)
526 do iene=indminn,indmaxx
533 !----------------------------------------------------------------------------
539 if (print_mc.gt.0) write (iout,'(80(1h*)/20x,a,i7)') &
540 'Beginning iteration #',it
541 ! Initialize local counter.
542 ntrial=0 ! # of generated non-overlapping confs.
543 noverlap=0 ! # of overlapping confs.
545 do while (.not. accepted .and. WhatsUp.eq.0 .and. Kwita.eq.0)
547 ! Retrieve the angles of previously accepted conformation
551 !d write (iout,'(a)') 'Old variables:'
552 !d write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
553 call var_to_geom(nvar,varia)
559 ! Decide whether to generate a random conformation or perturb the old one
560 RandOrPert=ran_number(0.0D0,1.0D0)
561 if (RandOrPert.gt.RanFract) then
563 write (iout,'(a)') 'Perturbation-generated conformation.'
564 call perturb(error,lprint,MoveType,nbond,1.0D0)
566 if (MoveType.lt.1 .or. MoveType.gt.MaxMoveType) then
567 write (iout,'(/a,i7,a/)') 'Error - unknown MoveType=',&
568 MoveType,' returned from PERTURB.'
575 nstart_grow=iran_num(3,nres)
577 write (iout,'(2a,i3)') 'Random-generated conformation',&
578 ' - chain regrown from residue',nstart_grow
579 call gen_rand_conf(nstart_grow,*30)
581 call geom_to_var(nvar,varia)
582 !d write (iout,'(a)') 'New variables:'
583 !d write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
585 if (print_mc.gt.0) write (iout,'(a,i5,a,i10,a,i10)') &
586 'Processor',MyId,' trial move',ntrial,' total generated:',ngen
587 if (print_mc.gt.0) write (*,'(a,i5,a,i10,a,i10)') &
588 'Processor',MyId,' trial move',ntrial,' total generated:',ngen
591 ! call enerprint(energia(0))
592 ! write (iout,'(2(a,1pe14.5))') 'Etot=',Etot,' Elowest=',Elowest
593 if (etot-elowest.gt.overlap_cut) then
594 write (iout,'(a,i5,a,1pe14.5)') 'Iteration',it,&
595 ' Overlap detected in the current conf.; energy is',etot
599 if (noverlap.gt.maxoverlap) then
600 write (iout,'(a)') 'Too many overlapping confs.'
605 call minimize(etot,varia,iretcode,nfun)
606 !d write (iout,'(a)') 'Variables after minimization:'
607 !d write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
610 neneval=neneval+nfun+1
612 if (print_mc.gt.2) then
613 write (iout,'(a)') 'Total energies of trial conf:'
614 call enerprint(energia)
615 else if (print_mc.eq.1) then
616 write (iout,'(a,i6,a,1pe16.6)') &
617 'Trial conformation:',ngen,' energy:',etot
619 !--------------------------------------------------------------------------
621 !--------------------------------------------------------------------------
624 call accepting(etot,eold,scur,sold,varia,varold,&
629 if (elowest.gt.etot) elowest=etot
630 if (ehighest.lt.etot) ehighest=etot
631 moves_acc(MoveType)=moves_acc(MoveType)+1
632 if (MoveType.eq.1) then
633 nbond_acc(nbond)=nbond_acc(nbond)+1
635 ! Check against conformation repetitions.
636 irep=conf_comp(varia,etot)
637 #if defined(AIX) || defined(PGI)
638 open (istat,file=statname,position='append')
640 open (istat,file=statname,access='append')
644 call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,kkk),&
648 call contact(.false.,ncont,icont,co)
649 frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
651 write (iout,'(a,f8.3,a,f8.3,a,f8.3)') &
652 'RMS deviation from the reference structure:',rms,&
653 ' % of native contacts:',frac*100,' contact order:',co
655 write (istat,'(i5,11(1pe14.5))') it,&
656 (energia(print_order(i)),i=1,nprint_ene),etot,&
658 elseif (print_stat) then
659 write (istat,'(i5,10(1pe14.5))') it,&
660 (energia(print_order(i)),i=1,nprint_ene),etot
664 call statprint(nacc,nfun,iretcode,etot,elowest)
665 ! Print internal coordinates.
666 if (print_int) call briefout(nacc,etot)
668 if (MyID.ne.MasterID) then
669 call recv_stop_sig(Kwita)
670 !d print *,'Processor:',MyID,' STOP=',Kwita
672 call send_MCM_info(2)
674 call send_MCM_info(1)
678 ! Store the accepted conf. and its energy.
686 !d write (iout,*) 'Accepted conformation:'
687 !d write (iout,*) (rad2deg*varia(i),i=1,nphi)
688 if (minim) call zapis(varia,etot)
690 ener(i,nsave)=energia(i)
692 ener(n_ene+1,nsave)=etot
693 ener(n_ene+2,nsave)=frac
695 nminima(irep)=nminima(irep)+1.0D0
696 ! print *,'irep=',irep,' nminima=',nminima(irep)
698 if (Kwita.eq.0) call recv_stop_sig(kwita)
703 if (MyID.eq.MasterID) then
704 call receive_MCM_info
705 if (nacc_tot.ge.maxacc) accepted=.true.
708 if (ntrial.gt.maxtrial_iter .and. npool.gt.0) then
709 ! Take a conformation from the pool
714 write (iout,*) 'Iteration',it,' max. # of trials exceeded.'
716 'Take conformation',ii,' from the pool energy=',epool(ii)
718 write (iout,'(10f8.3)') (rad2deg*varia(i),i=1,nvar)
720 endif ! (ntrial.gt.maxtrial_iter .and. npool.gt.0)
724 if (MyID.eq.MasterID) then
725 call receive_MCM_info
727 if (Kwita.eq.0) call recv_stop_sig(kwita)
729 if (ovrtim()) WhatsUp=-1
730 !d write (iout,*) 'WhatsUp=',WhatsUp,' Kwita=',Kwita
731 not_done = (nacc_tot.lt.maxacc) .and. (WhatsUp.eq.0) &
733 !d write (iout,*) 'not_done=',not_done
736 print *,'Processor',MyID,&
737 ' has received STOP signal =',Kwita,' in EntSamp.'
738 !d print *,'not_done=',not_done
739 if (Kwita.lt.-1) WhatsUp=Kwita
740 else if (nacc_tot.ge.maxacc) then
741 print *,'Processor',MyID,' calls send_stop_sig,',&
742 ' because a sufficient # of confs. have been collected.'
743 !d print *,'not_done=',not_done
744 call send_stop_sig(-1)
745 else if (WhatsUp.eq.-1) then
746 print *,'Processor',MyID,&
747 ' calls send_stop_sig because of timeout.'
748 !d print *,'not_done=',not_done
749 call send_stop_sig(-2)
754 !-----------------------------------------------------------------
755 !... Construct energy histogram & update entropy
756 !-----------------------------------------------------------------
760 write (iout,*) 'Processor',MyID,&
761 ' is broadcasting ERROR-STOP signal.'
762 write (*,*) 'Processor',MyID,&
763 ' is broadcasting ERROR-STOP signal.'
764 call send_stop_sig(-3)
768 if (MyID.eq.MasterID) then
769 ! call receive_MCM_results
770 call receive_energies
773 if (esave(i).lt.elowest) elowest=esave(i)
774 if (esave(i).gt.ehighest) ehighest=esave(i)
776 write (iout,'(a,i10)') '# of accepted confs:',nacc_tot
777 write (iout,'(a,f10.5,a,f10.5)') 'Lowest energy:',elowest,&
778 ' Highest energy',ehighest
779 if (isweep.eq.1 .and. .not.ent_read) then
782 write (iout,*) 'EMAX=',emax
784 indmaxx=(ehighest-emin)/delte
787 do i=-max_ene,max_ene
788 entropy(i)=(emin+i*delte)*betbol
792 indmin=(elowest-emin)/delte
793 indmax=(ehighest-emin)/delte
794 if (indmin.lt.indminn) indminn=indmin
795 if (indmax.gt.indmaxx) indmaxx=indmax
797 write(iout,*)'indminn=',indminn,' indmaxx=',indmaxx
798 ! Construct energy histogram
800 inde=(esave(i)-emin)/delte
801 nhist(inde)=nhist(inde)+nminima(i)
803 ! Update entropy (density of states)
805 if (nhist(i).gt.0) then
806 entropy(i)=entropy(i)+dlog(nhist(i)+0.0D0)
810 !d entropy(i)=1.0D+10
812 write (iout,'(/80(1h*)/a,i2/80(1h*)/)') &
813 'End of macroiteration',isweep
814 write (iout,'(a,f10.5,a,f10.5)') 'Elowest=',elowest,&
815 ' Ehighest=',ehighest
816 write (iout,'(a)') 'Frequecies of minima'
818 write (iout,'(i5,f5.0,f10.5)') i,nminima(i),esave(i)
820 write (iout,'(/a)') 'Energy histogram'
822 write (iout,'(i5,2f10.5)') i,emin+i*delte,nhist(i)
824 write (iout,'(/a)') 'Entropy'
826 write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
828 !-----------------------------------------------------------------
829 !... End of energy histogram construction
830 !-----------------------------------------------------------------
832 entropy(-max_ene-4)=dfloat(indminn)
833 entropy(-max_ene-3)=dfloat(indmaxx)
834 entropy(-max_ene-2)=emin
835 entropy(-max_ene-1)=emax
837 !d print *,entname,ientout
838 open (ientout,file=entname,status='unknown')
839 write (ientout,'(2i5,2e25.17)') indminn,indmaxx,emin,emax
841 write (ientout,'(i5,f10.5,f20.15)') i,emin+i*delte,entropy(i)
845 write (iout,'(a)') 'Frequecies of minima'
847 write (iout,'(i5,f5.0,f10.5)') i,nminima(i),esave(i)
849 ! call send_MCM_results
851 call receive_MCM_update
852 indminn=entropy(-max_ene-4)
853 indmaxx=entropy(-max_ene-3)
854 emin=entropy(-max_ene-2)
855 emax=entropy(-max_ene-1)
856 write (iout,*) 'Received from master:'
857 write (iout,*) 'indminn=',indminn,' indmaxx=',indmaxx,&
858 ' emin=',emin,' emax=',emax
859 write (iout,'(/a)') 'Entropy'
861 write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
864 if (WhatsUp.lt.-1) return
866 if (ovrtim() .or. WhatsUp.lt.0) return
869 write (iout,'(/80(1h*)/20x,a)') 'Summary run statistics:'
870 call statprint(nacc,nfun,iretcode,etot,elowest)
872 'Statistics of multiple-bond motions. Total motions:'
873 write (iout,'(16i5)') (nbond_move(i),i=1,Nbm)
874 write (iout,'(a)') 'Accepted motions:'
875 write (iout,'(16i5)') (nbond_acc(i),i=1,Nbm)
876 !el write (iout,'(a,i10)') 'Number of chain regrowths:',nregrow
877 !el write (iout,'(a,i10)') 'Accepted chain regrowths:',nregrow_acc
879 !---------------------------------------------------------------------------
881 !---------------------------------------------------------------------------
885 if (isweep.eq.nsweep .and. it.ge.maxacc) &
886 write (iout,'(/80(1h*)/20x,a/80(1h*)/)') 'All iterations done.'
888 end subroutine entmcm
889 !-----------------------------------------------------------------------------
890 subroutine accepting(ecur,eold,scur,sold,x,xold,accepted)
892 use geometry_data, only: nphi
893 use energy_data, only: max_ene
894 ! implicit real*8 (a-h,o-z)
895 ! include 'DIMENSIONS'
896 ! include 'COMMON.MCM'
897 ! include 'COMMON.MCE'
898 ! include 'COMMON.IOUNITS'
899 ! include 'COMMON.VAR'
901 use MPI_data !include 'COMMON.INFO'
903 ! include 'COMMON.GEO'
904 real(kind=8) :: ecur,eold,xx,bol !,ran_number
905 real(kind=8),dimension(6*nres) :: x,xold !(maxvar) (maxvar=6*maxres)
906 real(kind=8) :: tole=1.0D-1, tola=5.0D0
909 !!! local variables - el
911 real(kind=8) :: scur,sold,xxh,deix,dent
913 ! Check if the conformation is similar.
914 !d write (iout,*) 'Enter ACCEPTING'
915 !d write (iout,*) 'Old PHI angles:'
916 !d write (iout,*) (rad2deg*xold(i),i=1,nphi)
917 !d write (iout,*) 'Current angles'
918 !d write (iout,*) (rad2deg*x(i),i=1,nphi)
919 !d ddif=dif_ang(nphi,x,xold)
920 !d write (iout,*) 'Angle norm:',ddif
921 !d write (iout,*) 'ecur=',ecur,' emax=',emax
922 if (ecur.gt.emax) then
925 write (iout,'(a)') 'Conformation rejected as too high in energy'
927 else if (dabs(ecur-eold).lt.tole .and. &
928 dif_ang(nphi,x,xold).lt.tola) then
931 write (iout,'(a)') 'Conformation rejected as too similar'
934 ! Else evaluate the entropy of the conf and compare it with that of the previous
936 indecur=(ecur-emin)/delte
937 if (iabs(indecur).gt.max_ene) then
938 write (iout,'(a,2i5)') &
939 'Accepting: Index out of range:',indecur
941 else if (indecur.eq.indmaxx) then
942 scur=entropy(indecur)
943 if (print_mc.gt.0) write (iout,*)'Energy boundary reached',&
944 indmaxx,indecur,entropy(indecur)
946 deix=ecur-(emin+indecur*delte)
947 dent=entropy(indecur+1)-entropy(indecur)
948 scur=entropy(indecur)+(dent/delte)*deix
950 !d print *,'Processor',MyID,' ecur=',ecur,' indecur=',indecur,
951 !d & ' scur=',scur,' eold=',eold,' sold=',sold
952 !d print *,'deix=',deix,' dent=',dent,' delte=',delte
953 if (print_mc.gt.1) then
954 write(iout,*)'ecur=',ecur,' indecur=',indecur,' scur=',scur
955 write(iout,*)'eold=',eold,' sold=',sold
957 if (scur.le.sold) then
960 ! Else carry out acceptance test
961 xx=ran_number(0.0D0,1.0D0)
963 if (xxh.gt.50.0D0) then
970 if (print_mc.gt.0) write (iout,'(a)') &
971 'Conformation accepted.'
974 if (print_mc.gt.0) write (iout,'(a)') &
975 'Conformation rejected.'
979 end subroutine accepting
980 !-----------------------------------------------------------------------------
983 use io_base, only:read_angles
984 ! implicit real*8 (a-h,o-z)
985 ! include 'DIMENSIONS'
986 ! include 'COMMON.IOUNITS'
987 ! include 'COMMON.GEO'
988 ! include 'COMMON.MCM'
989 ! include 'COMMON.MCE'
990 ! include 'COMMON.VAR'
991 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
993 !!! local variables - el
996 print '(a)','Call READ_POOL'
999 read (intin,'(i5,f10.5)',end=10,err=10) iconf,epool(npool)
1000 if (epool(npool).eq.0.0D0) goto 10
1001 call read_angles(intin,*10)
1002 call geom_to_var(nvar,xpool(1,npool))
1006 11 write (iout,'(a,i5)') 'Number of pool conformations:',npool
1007 if (print_mc.gt.2) then
1009 write (iout,'(a,i5,a,1pe14.5)') 'Pool conformation',i,' energy',&
1011 write (iout,'(10f8.3)') (rad2deg*xpool(j,i),j=1,nvar)
1013 endif ! (print_mc.gt.2)
1015 end subroutine read_pool
1016 !-----------------------------------------------------------------------------
1018 !-----------------------------------------------------------------------------
1019 subroutine monte_carlo
1023 use MPI_data, only:ifinish,nctasks,WhatsUp,MyID
1024 use control_data, only:refstr,MaxProcs
1026 use control, only:tcpu,ovrtim
1027 use regularize_, only:fitsq
1029 ! Does Boltzmann and entropic sampling without energy minimization
1030 ! implicit real*8 (a-h,o-z)
1031 ! include 'DIMENSIONS'
1032 ! include 'COMMON.IOUNITS'
1034 use MPI_data !include 'COMMON.INFO'
1036 ! include 'COMMON.GEO'
1037 ! include 'COMMON.CHAIN'
1038 ! include 'COMMON.MCM'
1039 ! include 'COMMON.MCE'
1040 ! include 'COMMON.CONTACTS'
1041 ! include 'COMMON.CONTROL'
1042 ! include 'COMMON.VAR'
1043 ! include 'COMMON.INTERACT'
1044 ! include 'COMMON.THREAD'
1045 ! include 'COMMON.NAMES'
1046 logical :: accepted,not_done,over,error,lprint !,ovrtim
1047 integer :: MoveType,nbond,nbins
1048 ! integer :: conf_comp
1049 real(kind=8) :: RandOrPert
1050 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
1051 real(kind=8) :: elowest,elowest1,ehighest,ehighest1,eold
1052 real(kind=8) :: przes(3),obr(3,3)
1053 real(kind=8),dimension(6*nres) :: varold !(maxvar) (maxvar=6*maxres)
1055 integer,dimension(-1:MaxMoveType+1,0:MaxProcs-1) :: moves1,moves_acc1 !(-1:MaxMoveType+1,0:MaxProcs-1)
1057 real(kind=8) :: etot_temp,etot_all(0:MaxProcs)
1058 external d_vadd,d_vmin,d_vmax
1059 real(kind=8),dimension(-max_ene:max_ene) :: entropy1,nhist1
1060 integer,dimension(nres*(MaxProcs+1)) :: nbond_move1,nbond_acc1
1061 integer,dimension(2) :: itemp
1063 real(kind=8),dimension(6*nres) :: var_lowest !(maxvar) (maxvar=6*maxres)
1064 real(kind=8),dimension(0:n_ene) :: energia,energia_ave
1066 !!! local variables - el
1067 integer :: i,j,it,ii,iproc,nacc,ISWEEP,nfun,indmax,indmin,ijunk,&
1068 Kwita,indeold,imdmax,inde,iretcode,nstart_grow,noverlap
1069 real(kind=8) :: facee,conste,ejunk,etot,sold,frac,runtime,&
1070 frac_ave,rms_ave,etot_ave,scur,from_pool,co,rms
1072 write(iout,'(a,i8,2x,a,f10.5)') &
1073 'pool_read_freq=',pool_read_freq,' pool_fraction=',pool_fraction
1074 open (istat,file=statname)
1078 facee=1.0D0/(maxacc*delte)
1079 ! Number of bins in energy histogram
1081 write (iout,*) 'NBINS=',nbins
1083 ! Read entropy from previous simulations.
1085 read (ientin,*) indminn,indmaxx,emin,emax
1086 print *,'indminn=',indminn,' indmaxx=',indmaxx,' emin=',emin,&
1088 do i=-max_ene,max_ene
1091 read (ientin,*) (ijunk,ejunk,entropy(i),i=indminn,indmaxx)
1094 write (iout,*) 'indminn=',indminn,' indmaxx=',indmaxx,&
1095 ' emin=',emin,' emax=',emax
1096 write (iout,'(/a)') 'Initial entropy'
1097 do i=indminn,indmaxx
1098 write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
1101 ! Read the pool of conformations
1105 !----------------------------------------------------------------------------
1106 ! Entropy-sampling simulations with continually updated entropy;
1107 ! set NSWEEP=1 for Boltzmann sampling.
1108 ! Loop thru simulations
1109 !----------------------------------------------------------------------------
1110 allocate(ifinish(nctasks))
1113 ! Initialize the IFINISH array.
1120 !---------------------------------------------------------------------------
1121 ! Initialize counters.
1122 !---------------------------------------------------------------------------
1123 ! Total number of generated confs.
1125 ! Total number of moves. In general this won't be equal to the number of
1126 ! attempted moves, because we may want to reject some "bad" confs just by
1129 ! Total number of shift (nbond_move(1)), spike, crankshaft, three-bond,...
1131 !el allocate(nbond_move(nres)) !(maxres)
1132 !el allocate(nbond_acc(nres)) !(maxres)
1138 ! Initialize total and accepted number of moves of various kind.
1143 ! Total number of energy evaluations.
1146 !----------------------------------------------------------------------------
1147 ! Take a conformation from the pool
1148 !----------------------------------------------------------------------------
1150 write (iout,*) 'emin=',emin,' emax=',emax
1151 if (npool.gt.0) then
1152 ii=iran_num(1,npool)
1154 varia(i)=xpool(i,ii)
1156 write (iout,*) 'Took conformation',ii,' from the pool energy=',&
1158 call var_to_geom(nvar,varia)
1159 ! Print internal coordinates of the initial conformation
1161 else if (isweep.gt.1) then
1162 if (eold.lt.emax) then
1168 varia(i)=var_lowest(i)
1171 call var_to_geom(nvar,varia)
1173 !----------------------------------------------------------------------------
1174 ! Compute and print initial energies.
1175 !----------------------------------------------------------------------------
1179 write (iout,'(/80(1h*)/a,i2/80(1h*)/)') 'MCE iteration #',isweep
1180 write (iout,'(/80(1h*)/a)') 'Initial energies:'
1182 call geom_to_var(nvar,varia)
1183 call etotal(energia)
1185 call enerprint(energia)
1187 call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,1),nsup,przes,&
1190 call contact(.false.,ncont,icont,co)
1191 frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
1192 write (iout,'(a,f8.3,a,f8.3,a,f8.3)') &
1193 'RMS deviation from the reference structure:',rms,&
1194 ' % of native contacts:',frac*100,' contact order',co
1195 write (istat,'(i10,16(1pe14.5))') 0,&
1196 (energia(print_order(i)),i=1,nprint_ene),&
1199 write (istat,'(i10,14(1pe14.5))') 0,&
1200 (energia(print_order(i)),i=1,nprint_ene),etot
1204 if (.not. ent_read) then
1205 ! Initialize the entropy array
1207 ! Collect total energies from other processors.
1210 call mp_gather(etot_temp,etot_all,8,MasterID,cgGroupID)
1211 if (MyID.eq.MasterID) then
1212 ! Get the lowest and the highest energy.
1213 print *,'MASTER: etot_temp: ',(etot_all(i),i=0,nprocs-1),&
1214 ' emin=',emin,' emax=',emax
1218 if (emin.gt.etot_all(i)) emin=etot_all(i)
1219 if (emax.lt.etot_all(i)) emax=etot_all(i)
1222 endif ! MyID.eq.MasterID
1225 print *,'Processor',MyID,' calls MP_BCAST to send/recv etot_all'
1226 call mp_bcast(etot_all(1),16,MasterID,cgGroupID)
1227 print *,'Processor',MyID,' MP_BCAST to send/recv etot_all ended'
1228 if (MyID.ne.MasterID) then
1229 print *,'Processor:',MyID,etot_all(1),etot_all(2),&
1230 etot_all(1),etot_all(2)
1233 endif ! MyID.ne.MasterID
1234 write (iout,*) 'After MP_GATHER etot_temp=',&
1235 etot_temp,' emin=',emin
1243 ! Multicanonical sampling - start from Boltzmann distribution
1244 do i=-max_ene,max_ene
1245 entropy(i)=(emin+i*delte)*betbol
1248 ! Entropic sampling - start from uniform distribution of the density of states
1249 do i=-max_ene,max_ene
1253 write (iout,'(/a)') 'Initial entropy'
1254 do i=indminn,indmaxx
1255 write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
1257 if (isweep.eq.1) then
1261 indmaxx=indminn+nbins
1264 endif ! .not. ent_read
1266 call recv_stop_sig(Kwita)
1267 if (whatsup.eq.1) then
1268 call send_stop_sig(-2)
1270 else if (whatsup.le.-2) then
1272 else if (whatsup.eq.2) then
1280 write (iout,'(/80(1h*)/20x,a/80(1h*))') &
1281 'Enter Monte Carlo procedure.'
1283 call briefout(0,etot)
1288 call entropia(eold,sold,indeold)
1289 ! NACC is the counter for the accepted conformations of a given processor
1291 ! NACC_TOT counts the total number of accepted conformations
1294 !----------------------------------------------------------------------------
1295 ! Zero out average energies
1297 energia_ave(i)=0.0d0
1299 ! Initialize energy histogram
1300 do i=-max_ene,max_ene
1303 ! Zero out iteration counter.
1308 ! Begin MC iteration loop.
1311 ! Initialize local counter.
1312 ntrial=0 ! # of generated non-overlapping confs.
1313 noverlap=0 ! # of overlapping confs.
1315 do while (.not. accepted .and. WhatsUp.eq.0 .and. Kwita.eq.0)
1317 ! Retrieve the angles of previously accepted conformation
1321 call var_to_geom(nvar,varia)
1322 ! Rebuild the chain.
1327 ! Decide whether to take a conformation from the pool or generate/perturb one
1329 from_pool=ran_number(0.0D0,1.0D0)
1330 if (npool.gt.0 .and. from_pool.lt.pool_fraction) then
1331 ! Throw a dice to choose the conformation from the pool
1332 ii=iran_num(1,npool)
1334 varia(i)=xpool(i,ii)
1336 call var_to_geom(nvar,varia)
1339 !d write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
1340 if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) &
1341 write (iout,'(a,i3,a,f10.5)') &
1342 'Try conformation',ii,' from the pool energy=',epool(ii)
1344 moves(-1)=moves(-1)+1
1346 ! Decide whether to generate a random conformation or perturb the old one
1347 RandOrPert=ran_number(0.0D0,1.0D0)
1348 if (RandOrPert.gt.RanFract) then
1349 if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) &
1350 write (iout,'(a)') 'Perturbation-generated conformation.'
1351 call perturb(error,lprint,MoveType,nbond,0.1D0)
1353 if (MoveType.lt.1 .or. MoveType.gt.MaxMoveType) then
1354 write (iout,'(/a,i7,a/)') 'Error - unknown MoveType=',&
1355 MoveType,' returned from PERTURB.'
1362 nstart_grow=iran_num(3,nres)
1363 if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) &
1364 write (iout,'(2a,i3)') 'Random-generated conformation',&
1365 ' - chain regrown from residue',nstart_grow
1366 call gen_rand_conf(nstart_grow,*30)
1368 call geom_to_var(nvar,varia)
1370 !d write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
1372 if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) &
1373 write (iout,'(a,i5,a,i10,a,i10)') &
1374 'Processor',MyId,' trial move',ntrial,' total generated:',ngen
1375 if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) &
1376 write (*,'(a,i5,a,i10,a,i10)') &
1377 'Processor',MyId,' trial move',ntrial,' total generated:',ngen
1378 call etotal(energia)
1381 !d call enerprint(energia(0))
1382 !d write(iout,*)'it=',it,' etot=',etot
1383 if (etot-elowest.gt.overlap_cut) then
1384 if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) &
1385 write (iout,'(a,i5,a,1pe14.5)') 'Iteration',it,&
1386 ' Overlap detected in the current conf.; energy is',etot
1389 if (noverlap.gt.maxoverlap) then
1390 write (iout,'(a)') 'Too many overlapping confs.'
1394 !--------------------------------------------------------------------------
1395 !... Acceptance test
1396 !--------------------------------------------------------------------------
1399 call accept_mc(it,etot,eold,scur,sold,varia,varold,accepted)
1403 if (elowest.gt.etot) then
1406 var_lowest(i)=varia(i)
1409 if (ehighest.lt.etot) ehighest=etot
1410 moves_acc(MoveType)=moves_acc(MoveType)+1
1411 if (MoveType.eq.1) then
1412 nbond_acc(nbond)=nbond_acc(nbond)+1
1414 ! Compare with reference structure.
1416 call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,1),&
1417 nsup,przes,obr,non_conv)
1419 call contact(.false.,ncont,icont,co)
1420 frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
1423 ! Periodically save average energies and confs.
1426 energia_ave(i)=energia_ave(i)+energia(i)
1428 moves(MaxMoveType+1)=nmove
1429 moves_acc(MaxMoveType+1)=nacc
1430 IF ((it/save_frequency)*save_frequency.eq.it) THEN
1432 energia_ave(i)=energia_ave(i)/save_frequency
1434 etot_ave=energia_ave(0)
1436 ! open (istat,file=statname,position='append')
1438 ! open (istat,file=statname,access='append')
1440 if (print_mc.gt.0) &
1441 write (iout,'(80(1h*)/20x,a,i20)') &
1443 if (refstr .and. print_mc.gt.0) then
1444 write (iout,'(a,f8.3,a,f8.3,a,f8.3)') &
1445 'RMS deviation from the reference structure:',rms,&
1446 ' % of native contacts:',frac*100,' contact order:',co
1448 if (print_stat) then
1450 write (istat,'(i10,10(1pe14.5))') it,&
1451 (energia_ave(print_order(i)),i=1,nprint_ene),&
1452 etot_ave,rms_ave,frac_ave
1454 write (istat,'(i10,10(1pe14.5))') it,&
1455 (energia_ave(print_order(i)),i=1,nprint_ene),&
1460 if (print_mc.gt.0) &
1461 call statprint(nacc,nfun,iretcode,etot,elowest)
1462 ! Print internal coordinates.
1463 if (print_int) call briefout(nacc,etot)
1465 energia_ave(i)=0.0d0
1467 ENDIF ! ( (it/save_frequency)*save_frequency.eq.it)
1469 inde=icialosc((etot-emin)/delte)
1470 nhist(inde)=nhist(inde)+1.0D0
1472 if ( (it/message_frequency)*message_frequency.eq.it &
1473 .and. (MyID.ne.MasterID) ) then
1474 call recv_stop_sig(Kwita)
1475 call send_MCM_info(message_frequency)
1478 ! Store the accepted conf. and its energy.
1485 if (Kwita.eq.0) call recv_stop_sig(kwita)
1490 if (MyID.eq.MasterID .and. &
1491 (it/message_frequency)*message_frequency.eq.it) then
1492 call receive_MC_info
1493 if (nacc_tot.ge.maxacc) accepted=.true.
1496 ! if ((ntrial.gt.maxtrial_iter
1497 ! & .or. (it/pool_read_freq)*pool_read_freq.eq.it)
1498 ! & .and. npool.gt.0) then
1499 ! Take a conformation from the pool
1500 ! ii=iran_num(1,npool)
1502 ! varold(i)=xpool(i,ii)
1504 ! if (ntrial.gt.maxtrial_iter)
1505 ! & write (iout,*) 'Iteration',it,' max. # of trials exceeded.'
1507 ! & 'Take conformation',ii,' from the pool energy=',epool(ii)
1508 ! if (print_mc.gt.2)
1509 ! & write (iout,'(10f8.3)') (rad2deg*varold(i),i=1,nvar)
1512 ! call entropia(eold,sold,indeold)
1514 ! endif ! (ntrial.gt.maxtrial_iter .and. npool.gt.0)
1518 if (MyID.eq.MasterID .and. &
1519 (it/message_frequency)*message_frequency.eq.it) then
1520 call receive_MC_info
1522 if (Kwita.eq.0) call recv_stop_sig(kwita)
1524 if (ovrtim()) WhatsUp=-1
1525 !d write (iout,*) 'WhatsUp=',WhatsUp,' Kwita=',Kwita
1526 not_done = (nacc_tot.lt.maxacc) .and. (WhatsUp.eq.0) &
1528 !d write (iout,*) 'not_done=',not_done
1530 if (Kwita.lt.0) then
1531 print *,'Processor',MyID,&
1532 ' has received STOP signal =',Kwita,' in EntSamp.'
1533 !d print *,'not_done=',not_done
1534 if (Kwita.lt.-1) WhatsUp=Kwita
1535 if (MyID.ne.MasterID) call send_MCM_info(-1)
1536 else if (nacc_tot.ge.maxacc) then
1537 print *,'Processor',MyID,' calls send_stop_sig,',&
1538 ' because a sufficient # of confs. have been collected.'
1539 !d print *,'not_done=',not_done
1540 call send_stop_sig(-1)
1541 if (MyID.ne.MasterID) call send_MCM_info(-1)
1542 else if (WhatsUp.eq.-1) then
1543 print *,'Processor',MyID,&
1544 ' calls send_stop_sig because of timeout.'
1545 !d print *,'not_done=',not_done
1546 call send_stop_sig(-2)
1547 if (MyID.ne.MasterID) call send_MCM_info(-1)
1552 !-----------------------------------------------------------------
1553 !... Construct energy histogram & update entropy
1554 !-----------------------------------------------------------------
1558 write (iout,*) 'Processor',MyID,&
1559 ' is broadcasting ERROR-STOP signal.'
1560 write (*,*) 'Processor',MyID,&
1561 ' is broadcasting ERROR-STOP signal.'
1562 call send_stop_sig(-3)
1563 if (MyID.ne.MasterID) call send_MCM_info(-1)
1566 write (iout,'(/a)') 'Energy histogram'
1568 write (iout,'(i5,2f20.5)') i,emin+i*delte,nhist(i)
1571 ! Wait until every processor has sent complete MC info.
1572 if (MyID.eq.MasterID) then
1575 ! write (*,*) 'The IFINISH array:'
1576 ! write (*,*) (ifinish(i),i=1,nctasks)
1579 not_done=not_done.or.(ifinish(i).ge.0)
1581 if (not_done) call receive_MC_info
1584 ! Make collective histogram from the work of all processors.
1585 msglen=(2*max_ene+1)*8
1587 'Processor',MyID,' calls MP_REDUCE to send/receive histograms',&
1589 call mp_reduce(nhist,nhist1,msglen,MasterID,d_vadd,&
1591 print *,'Processor',MyID,' MP_REDUCE accomplished for histogr.'
1592 do i=-max_ene,max_ene
1595 ! Collect min. and max. energy
1597 'Processor',MyID,' calls MP_REDUCE to send/receive energy borders'
1598 call mp_reduce(elowest,elowest1,8,MasterID,d_vmin,cgGroupID)
1599 call mp_reduce(ehighest,ehighest1,8,MasterID,d_vmax,cgGroupID)
1600 print *,'Processor',MyID,' MP_REDUCE accomplished for energies.'
1601 IF (MyID.eq.MasterID) THEN
1605 write (iout,'(a,i10)') '# of accepted confs:',nacc_tot
1606 write (iout,'(a,f10.5,a,f10.5)') 'Lowest energy:',elowest,&
1607 ' Highest energy',ehighest
1608 indmin=icialosc((elowest-emin)/delte)
1609 imdmax=icialosc((ehighest-emin)/delte)
1610 if (indmin.lt.indminn) then
1611 emax=emin+indmin*delte+e_up
1612 indmaxx=indmin+nbins
1615 if (.not.ent_read) ent_read=.true.
1616 write(iout,*)'indminn=',indminn,' indmaxx=',indmaxx
1617 ! Update entropy (density of states)
1619 if (nhist(i).gt.0) then
1620 entropy(i)=entropy(i)+dlog(nhist(i)+0.0D0)
1623 write (iout,'(/80(1h*)/a,i2/80(1h*)/)') &
1624 'End of macroiteration',isweep
1625 write (iout,'(a,f10.5,a,f10.5)') 'Elowest=',elowest,&
1626 ' Ehighest=',ehighest
1627 write (iout,'(/a)') 'Energy histogram'
1628 do i=indminn,indmaxx
1629 write (iout,'(i5,2f20.5)') i,emin+i*delte,nhist(i)
1631 write (iout,'(/a)') 'Entropy'
1632 do i=indminn,indmaxx
1633 write (iout,'(i5,2f20.5)') i,emin+i*delte,entropy(i)
1635 !-----------------------------------------------------------------
1636 !... End of energy histogram construction
1637 !-----------------------------------------------------------------
1640 if (.not. ent_read) ent_read=.true.
1641 ENDIF ! MyID .eq. MaterID
1642 if (MyID.eq.MasterID) then
1646 print *,'before mp_bcast processor',MyID,' indminn=',indminn,&
1647 ' indmaxx=',indmaxx,' itemp=',itemp(1),itemp(2)
1648 call mp_bcast(itemp(1),8,MasterID,cgGroupID)
1649 call mp_bcast(emax,8,MasterID,cgGroupID)
1650 print *,'after mp_bcast processor',MyID,' indminn=',indminn,&
1651 ' indmaxx=',indmaxx,' itemp=',itemp(1),itemp(2)
1652 if (MyID .ne. MasterID) then
1656 msglen=(indmaxx-indminn+1)*8
1657 print *,'processor',MyID,' calling mp_bcast msglen=',msglen,&
1658 ' indminn=',indminn,' indmaxx=',indmaxx,' isweep=',isweep
1659 call mp_bcast(entropy(indminn),msglen,MasterID,cgGroupID)
1660 IF(MyID.eq.MasterID .and. .not. ovrtim() .and. WhatsUp.ge.0)THEN
1661 open (ientout,file=entname,status='unknown')
1662 write (ientout,'(2i5,2e25.17)') indminn,indmaxx,emin,emax
1663 do i=indminn,indmaxx
1664 write (ientout,'(i5,f10.5,f20.15)') i,emin+i*delte,entropy(i)
1668 write (iout,*) 'Received from master:'
1669 write (iout,*) 'indminn=',indminn,' indmaxx=',indmaxx,&
1670 ' emin=',emin,' emax=',emax
1671 write (iout,'(/a)') 'Entropy'
1672 do i=indminn,indmaxx
1673 write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
1675 ENDIF ! MyID.eq.MasterID
1676 print *,'Processor',MyID,' calls MP_GATHER'
1677 call mp_gather(nbond_move,nbond_move1,4*Nbm,MasterID,&
1679 call mp_gather(nbond_acc,nbond_acc1,4*Nbm,MasterID,&
1681 print *,'Processor',MyID,' MP_GATHER call accomplished'
1682 if (MyID.eq.MasterID) then
1684 write (iout,'(/80(1h*)/20x,a)') 'Summary run statistics:'
1685 call statprint(nacc_tot,nfun,iretcode,etot,elowest)
1686 write (iout,'(a)') &
1687 'Statistics of multiple-bond motions. Total motions:'
1688 write (iout,'(8i10)') (nbond_move(i),i=1,Nbm)
1689 write (iout,'(a)') 'Accepted motions:'
1690 write (iout,'(8i10)') (nbond_acc(i),i=1,Nbm)
1692 write (iout,'(a)') &
1693 'Statistics of multi-bond moves of respective processors:'
1697 nbond_move(i)=nbond_move(i)+nbond_move1(ind)
1698 nbond_acc(i)=nbond_acc(i)+nbond_acc1(ind)
1702 write (iout,*) 'Processor',iproc,' nbond_move:', &
1703 (nbond_move1(iproc*nbm+i),i=1,Nbm),&
1704 ' nbond_acc:',(nbond_acc1(iproc*nbm+i),i=1,Nbm)
1707 call mp_gather(moves,moves1,4*(MaxMoveType+3),MasterID,&
1709 call mp_gather(moves_acc,moves_acc1,4*(MaxMoveType+3),&
1711 if (MyID.eq.MasterID) then
1713 do i=-1,MaxMoveType+1
1714 moves(i)=moves(i)+moves1(i,iproc)
1715 moves_acc(i)=moves_acc(i)+moves_acc1(i,iproc)
1719 do i=0,MaxMoveType+1
1720 nmove=nmove+moves(i)
1723 write (iout,*) 'Processor',iproc,' moves',&
1724 (moves1(i,iproc),i=0,MaxMoveType+1),&
1725 ' moves_acc:',(moves_acc1(i,iproc),i=0,MaxMoveType+1)
1729 open (ientout,file=entname,status='unknown')
1730 write (ientout,'(2i5,2e25.17)') indminn,indmaxx,emin,emax
1731 do i=indminn,indmaxx
1732 write (ientout,'(i5,f10.5,f20.15)') i,emin+i*delte,entropy(i)
1736 write (iout,'(/80(1h*)/20x,a)') 'Summary run statistics:'
1737 call statprint(nacc_tot,nfun,iretcode,etot,elowest)
1738 write (iout,'(a)') &
1739 'Statistics of multiple-bond motions. Total motions:'
1740 write (iout,'(8i10)') (nbond_move(i),i=1,Nbm)
1741 write (iout,'(a)') 'Accepted motions:'
1742 write (iout,'(8i10)') (nbond_acc(i),i=1,Nbm)
1743 if (ovrtim() .or. WhatsUp.lt.0) return
1745 !---------------------------------------------------------------------------
1747 !---------------------------------------------------------------------------
1751 if (isweep.eq.nsweep .and. it.ge.maxacc) &
1752 write (iout,'(/80(1h*)/20x,a/80(1h*)/)') 'All iterations done.'
1754 end subroutine monte_carlo
1755 !-----------------------------------------------------------------------------
1756 subroutine accept_mc(it,ecur,eold,scur,sold,x,xold,accepted)
1758 ! implicit real*8 (a-h,o-z)
1759 ! include 'DIMENSIONS'
1760 ! include 'COMMON.MCM'
1761 ! include 'COMMON.MCE'
1762 ! include 'COMMON.IOUNITS'
1763 ! include 'COMMON.VAR'
1765 use MPI_data !include 'COMMON.INFO'
1767 ! include 'COMMON.GEO'
1768 real(kind=8) :: ecur,eold,xx,bol
1769 real(kind=8),dimension(6*nres) :: x,xold !(maxvar) (maxvar=6*maxres)
1773 integer :: it,indecur
1774 real(kind=8) :: scur,sold,xxh
1775 ! Check if the conformation is similar.
1776 !d write (iout,*) 'Enter ACCEPTING'
1777 !d write (iout,*) 'Old PHI angles:'
1778 !d write (iout,*) (rad2deg*xold(i),i=1,nphi)
1779 !d write (iout,*) 'Current angles'
1780 !d write (iout,*) (rad2deg*x(i),i=1,nphi)
1781 !d ddif=dif_ang(nphi,x,xold)
1782 !d write (iout,*) 'Angle norm:',ddif
1783 !d write (iout,*) 'ecur=',ecur,' emax=',emax
1784 if (ecur.gt.emax) then
1786 if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) &
1787 write (iout,'(a)') 'Conformation rejected as too high in energy'
1790 ! Else evaluate the entropy of the conf and compare it with that of the previous
1792 call entropia(ecur,scur,indecur)
1793 !d print *,'Processor',MyID,' ecur=',ecur,' indecur=',indecur,
1794 !d & ' scur=',scur,' eold=',eold,' sold=',sold
1795 !d print *,'deix=',deix,' dent=',dent,' delte=',delte
1796 if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) then
1797 write(iout,*)'it=',it,'ecur=',ecur,' indecur=',indecur,&
1799 write(iout,*)'eold=',eold,' sold=',sold
1801 if (scur.le.sold) then
1804 ! Else carry out acceptance test
1805 xx=ran_number(0.0D0,1.0D0)
1807 if (xxh.gt.50.0D0) then
1814 if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) &
1815 write (iout,'(a)') 'Conformation accepted.'
1818 if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) &
1819 write (iout,'(a)') 'Conformation rejected.'
1823 end subroutine accept_mc
1824 !-----------------------------------------------------------------------------
1825 integer function icialosc(x)
1828 if (x.lt.0.0D0) then
1834 end function icialosc
1835 !-----------------------------------------------------------------------------
1836 subroutine entropia(ecur,scur,indecur)
1838 use energy_data, only: max_ene
1839 ! implicit real*8 (a-h,o-z)
1840 ! include 'DIMENSIONS'
1841 ! include 'COMMON.MCM'
1842 ! include 'COMMON.MCE'
1843 ! include 'COMMON.IOUNITS'
1844 real(kind=8) :: ecur,scur,deix,dent
1845 integer :: indecur,it !???el
1847 indecur=icialosc((ecur-emin)/delte)
1848 if (iabs(indecur).gt.max_ene) then
1849 if ((it/print_freq)*it.eq.it) write (iout,'(a,2i5)') &
1850 'Accepting: Index out of range:',indecur
1852 else if (indecur.ge.indmaxx) then
1853 scur=entropy(indecur)
1854 if (print_mc.gt.0 .and. (it/print_freq)*it.eq.it) &
1855 write (iout,*)'Energy boundary reached',&
1856 indmaxx,indecur,entropy(indecur)
1858 deix=ecur-(emin+indecur*delte)
1859 dent=entropy(indecur+1)-entropy(indecur)
1860 scur=entropy(indecur)+(dent/delte)*deix
1863 end subroutine entropia
1864 !-----------------------------------------------------------------------------
1866 !-----------------------------------------------------------------------------
1867 subroutine mcm_setup
1870 ! implicit real*8 (a-h,o-z)
1871 ! include 'DIMENSIONS'
1872 ! include 'COMMON.IOUNITS'
1873 ! include 'COMMON.MCM'
1874 ! include 'COMMON.CONTROL'
1875 ! include 'COMMON.INTERACT'
1876 ! include 'COMMON.NAMES'
1877 ! include 'COMMON.CHAIN'
1878 ! include 'COMMON.VAR'
1880 !!! local variables - el
1881 integer :: i,i1,i2,it1,it2,ngly,mmm,maxwinlen
1883 ! Set up variables used in MC/MCM.
1885 ! allocate(sumpro_bond(0:nres)) !(0:maxres)
1887 write (iout,'(80(1h*)/20x,a/80(1h*))') 'MCM control parameters:'
1888 write (iout,'(5(a,i7))') 'Maxacc:',maxacc,' MaxTrial:',MaxTrial,&
1889 ' MaxRepm:',MaxRepm,' MaxGen:',MaxGen,' MaxOverlap:',MaxOverlap
1890 write (iout,'(4(a,f8.1)/2(a,i3))') &
1891 'Tmin:',Tmin,' Tmax:',Tmax,' TstepH:',TstepH,&
1892 ' TstepC:',TstepC,'NstepH:',NstepH,' NstepC:',NstepC
1893 if (nwindow.gt.0) then
1894 write (iout,'(a)') 'Perturbation windows:'
1900 write (iout,'(a,i3,a,i3,a,i3)') restyp(it1),i1,restyp(it2),i2,&
1904 ! Rbolt=8.3143D-3*2.388459D-01 kcal/(mol*K)
1906 ! Number of "end bonds".
1909 print *,'koniecl=',koniecl
1910 write (iout,'(a)') 'Probabilities of move types:'
1911 write (*,'(a)') 'Probabilities of move types:'
1913 write (iout,'(a,f10.5)') MovTypID(i),&
1914 sumpro_type(i)-sumpro_type(i-1)
1915 write (*,'(a,f10.5)') MovTypID(i),&
1916 sumpro_type(i)-sumpro_type(i-1)
1919 ! Maximum length of N-bond segment to be moved
1920 ! nbm=nres-1-(2*koniecl-1)
1921 if (nwindow.gt.0) then
1924 if (winlen(i).gt.maxwinlen) maxwinlen=winlen(i)
1926 nbm=min0(maxwinlen,6)
1927 write (iout,'(a,i3,a,i3)') 'Nbm=',Nbm,' Maxwinlen=',Maxwinlen
1931 sumpro_bond(0)=0.0D0
1932 sumpro_bond(1)=0.0D0
1934 sumpro_bond(i)=sumpro_bond(i-1)+1.0D0/dfloat(i)
1936 write (iout,'(a)') 'The SumPro_Bond array:'
1937 write (iout,'(8f10.5)') (sumpro_bond(i),i=1,nbm)
1938 write (*,'(a)') 'The SumPro_Bond array:'
1939 write (*,'(8f10.5)') (sumpro_bond(i),i=1,nbm)
1940 ! Maximum number of side chains moved simultaneously
1941 ! print *,'nnt=',nnt,' nct=',nct
1944 if (itype(i).eq.10) ngly=ngly+1
1948 MaxSideMove=min0((nct-nnt+1)/2,mmm)
1950 ! print *,'MaxSideMove=',MaxSideMove
1951 ! Max. number of generated confs (not used at present).
1953 ! Set initial temperature
1955 betbol=1.0D0/(Rbol*Tcur)
1956 write (iout,'(a,f8.1,a,f10.5)') 'Initial temperature:',Tcur,&
1958 write (iout,*) 'RanFract=',ranfract
1960 end subroutine mcm_setup
1961 !-----------------------------------------------------------------------------
1963 subroutine do_mcm(i_orig)
1967 use MPI_data, only:Whatsup
1968 use control_data, only:refstr,minim,iprint
1970 use control, only:tcpu,ovrtim
1971 use regularize_, only:fitsq
1973 use minimm, only:minimize
1974 ! Monte-Carlo-with-Minimization calculations - serial code.
1975 ! implicit real*8 (a-h,o-z)
1976 ! include 'DIMENSIONS'
1977 ! include 'COMMON.IOUNITS'
1978 ! include 'COMMON.GEO'
1979 ! include 'COMMON.CHAIN'
1980 ! include 'COMMON.MCM'
1981 ! include 'COMMON.CONTACTS'
1982 ! include 'COMMON.CONTROL'
1983 ! include 'COMMON.VAR'
1984 ! include 'COMMON.INTERACT'
1985 ! include 'COMMON.CACHE'
1986 !rc include 'COMMON.DEFORM'
1987 !rc include 'COMMON.DEFORM1'
1988 ! include 'COMMON.NAMES'
1989 logical :: accepted,over,error,lprint,not_done,my_conf,&
1990 enelower,non_conv !,ovrtim
1991 integer :: MoveType,nbond !,conf_comp
1992 integer,dimension(max_cache) :: ifeed
1993 real(kind=8),dimension(6*nres) :: varia,varold !(maxvar) (maxvar=6*maxres)
1994 real(kind=8) :: elowest,eold
1995 real(kind=8) :: przes(3),obr(3,3)
1996 real(kind=8) :: energia(0:n_ene)
1997 real(kind=8) :: coord1(6*nres,max_thread2),enetb1(max_threadss) !el
1998 !!! local variables - el
1999 integer :: i,nf,nacc,it,nout,j,i_orig,nfun,Kwita,iretcode,&
2000 noverlap,nstart_grow,irepet,n_thr,ii
2001 real(kind=8) :: etot,frac,rms,co,RandOrPert,&
2003 !---------------------------------------------------------------------------
2004 ! Initialize counters.
2005 !---------------------------------------------------------------------------
2006 ! Total number of generated confs.
2008 ! Total number of moves. In general this won't be equal to the number of
2009 ! attempted moves, because we may want to reject some "bad" confs just by
2012 ! Total number of temperature jumps.
2014 ! Total number of shift (nbond_move(1)), spike, crankshaft, three-bond,...
2016 ! if(.not.allocated(varsave)) allocate(varsave(maxvar,maxsave)) !(maxvar,maxsave)
2017 ! allocate(nbond_move(nres)) !(maxres)
2023 ! Initialize total and accepted number of moves of various kind.
2028 ! Total number of energy evaluations.
2033 write (iout,*) 'RanFract=',RanFract
2038 !----------------------------------------------------------------------------
2039 ! Compute and print initial energies.
2040 !----------------------------------------------------------------------------
2042 write (iout,'(/80(1h*)/a)') 'Initial energies:'
2046 call etotal(energia)
2048 ! Minimize the energy of the first conformation.
2050 call geom_to_var(nvar,varia)
2051 ! write (iout,*) 'The VARIA array'
2052 ! write (iout,'(8f10.4)') (rad2deg*varia(i),i=1,nvar)
2053 call minimize(etot,varia,iretcode,nfun)
2054 call var_to_geom(nvar,varia)
2056 write (iout,*) 'etot from MINIMIZE:',etot
2057 ! write (iout,*) 'Tha VARIA array'
2058 ! write (iout,'(8f10.4)') (rad2deg*varia(i),i=1,nvar)
2060 call etotal(energia)
2062 call enerprint(energia)
2065 call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,1),nsup,przes,& !el cref(1,nstart_sup)
2068 call contact(.false.,ncont,icont,co)
2069 frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
2070 write (iout,'(a,f8.3,a,f8.3,a,f8.3)') &
2071 'RMS deviation from the reference structure:',rms,&
2072 ' % of native contacts:',frac*100,' contact order:',co
2074 write (istat,'(i5,17(1pe14.5))') 0,&
2075 (energia(print_order(i)),i=1,nprint_ene),&
2078 if (print_stat) write (istat,'(i5,16(1pe14.5))') 0,&
2079 (energia(print_order(i)),i=1,nprint_ene),etot
2081 if (print_stat) close(istat)
2082 neneval=neneval+nfun+1
2083 write (iout,'(/80(1h*)/20x,a/80(1h*))') &
2084 'Enter Monte Carlo procedure.'
2087 call briefout(0,etot)
2094 call zapis(varia,etot)
2095 nacc=0 ! total # of accepted confs of the current processor.
2096 nacc_tot=0 ! total # of accepted confs of all processors.
2098 not_done = (iretcode.ne.11)
2100 !----------------------------------------------------------------------------
2102 !----------------------------------------------------------------------------
2107 write (iout,'(80(1h*)/20x,a,i7)') &
2108 'Beginning iteration #',it
2109 ! Initialize local counter.
2110 ntrial=0 ! # of generated non-overlapping confs.
2112 do while (.not. accepted)
2114 ! Retrieve the angles of previously accepted conformation
2115 noverlap=0 ! # of overlapping confs.
2119 call var_to_geom(nvar,varia)
2120 ! Rebuild the chain.
2122 ! Heat up the system, if necessary.
2124 ! If temperature cannot be further increased, stop.
2129 !d write (iout,'(a)') 'Old variables:'
2130 !d write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
2131 ! Decide whether to generate a random conformation or perturb the old one
2132 RandOrPert=ran_number(0.0D0,1.0D0)
2133 if (RandOrPert.gt.RanFract) then
2134 if (print_mc.gt.0) &
2135 write (iout,'(a)') 'Perturbation-generated conformation.'
2136 call perturb(error,lprint,MoveType,nbond,1.0D0)
2138 if (MoveType.lt.1 .or. MoveType.gt.MaxMoveType) then
2139 write (iout,'(/a,i7,a/)') 'Error - unknown MoveType=',&
2140 MoveType,' returned from PERTURB.'
2147 nstart_grow=iran_num(3,nres)
2148 if (print_mc.gt.0) &
2149 write (iout,'(2a,i3)') 'Random-generated conformation',&
2150 ' - chain regrown from residue',nstart_grow
2151 call gen_rand_conf(nstart_grow,*30)
2153 call geom_to_var(nvar,varia)
2154 !d write (iout,'(a)') 'New variables:'
2155 !d write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
2158 call etotal(energia)
2160 ! call enerprint(energia(0))
2161 ! write (iout,'(2(a,1pe14.5))') 'Etot=',Etot,' Elowest=',Elowest
2162 if (etot-elowest.gt.overlap_cut) then
2163 if(iprint.gt.1.or.etot.lt.1d20) &
2164 write (iout,'(a,1pe14.5)') &
2165 'Overlap detected in the current conf.; energy is',etot
2169 if (noverlap.gt.maxoverlap) then
2170 write (iout,'(a)') 'Too many overlapping confs.'
2175 call minimize(etot,varia,iretcode,nfun)
2176 !d write (iout,*) 'etot from MINIMIZE:',etot
2177 !d write (iout,'(a)') 'Variables after minimization:'
2178 !d write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
2180 call etotal(energia)
2182 neneval=neneval+nfun+2
2184 ! call enerprint(energia(0))
2185 write (iout,'(a,i6,a,1pe16.6)') 'Conformation:',ngen,&
2187 !--------------------------------------------------------------------------
2188 !... Do Metropolis test
2189 !--------------------------------------------------------------------------
2193 if (WhatsUp.eq.0 .and. Kwita.eq.0) then
2194 call metropolis(nvar,varia,varold,etot,eold,accepted,&
2197 write (iout,*) 'My_Conf=',My_Conf,' EneLower=',EneLower
2202 if (elowest.gt.etot) elowest=etot
2203 moves_acc(MoveType)=moves_acc(MoveType)+1
2204 if (MoveType.eq.1) then
2205 nbond_acc(nbond)=nbond_acc(nbond)+1
2207 ! Check against conformation repetitions.
2208 irepet=conf_comp(varia,etot)
2209 if (print_stat) then
2210 #if defined(AIX) || defined(PGI)
2211 open (istat,file=statname,position='append')
2213 open (istat,file=statname,access='append')
2216 call statprint(nacc,nfun,iretcode,etot,elowest)
2218 call var_to_geom(nvar,varia)
2220 call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,1),& !el cref(1,nstart_sup)
2221 nsup,przes,obr,non_conv)
2223 call contact(.false.,ncont,icont,co)
2224 frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
2225 write (iout,'(a,f8.3,a,f8.3)') &
2226 'RMS deviation from the reference structure:',rms,&
2227 ' % of native contacts:',frac*100,' contact order',co
2231 write (iout,*) 'Writing new conformation',nout
2233 write (istat,'(i5,16(1pe14.5))') nout,&
2234 (energia(print_order(i)),i=1,nprint_ene),&
2238 write (istat,'(i5,17(1pe14.5))') nout,&
2239 (energia(print_order(i)),i=1,nprint_ene),etot
2241 if (print_stat) close(istat)
2242 ! Print internal coordinates.
2243 if (print_int) call briefout(nout,etot)
2244 ! Accumulate the newly accepted conf in the coord1 array, if it is different
2245 ! from all confs that are already there.
2246 call compare_s1(n_thr,max_thread2,etot,varia,ii,&
2247 enetb1,coord1,rms_deform,.true.,iprint)
2248 write (iout,*) 'After compare_ss: n_thr=',n_thr
2249 if (ii.eq.1 .or. ii.eq.3) then
2250 write (iout,'(8f10.4)') &
2251 (rad2deg*coord1(i,n_thr),i=1,nvar)
2254 write (iout,*) 'Conformation from cache, not written.'
2257 if (nrepm.gt.maxrepm) then
2258 write (iout,'(a)') 'Too many conformation repetitions.'
2261 ! Store the accepted conf. and its energy.
2266 if (irepet.eq.0) call zapis(varia,etot)
2267 ! Lower the temperature, if necessary.
2278 ! Check for time limit.
2279 if (ovrtim()) WhatsUp=-1
2280 not_done = (nacc_tot.lt.maxacc) .and. (WhatsUp.eq.0) &
2289 write (iout,'(/80(1h*)/20x,a)') 'Summary run statistics:'
2290 call statprint(nacc,nfun,iretcode,etot,elowest)
2291 write (iout,'(a)') &
2292 'Statistics of multiple-bond motions. Total motions:'
2293 write (iout,'(16i5)') (nbond_move(i),i=1,Nbm)
2294 write (iout,'(a)') 'Accepted motions:'
2295 write (iout,'(16i5)') (nbond_acc(i),i=1,Nbm)
2297 write (iout,'(/80(1h*)/20x,a/80(1h*)/)') 'All iterations done.'
2298 !(maxvar) (maxvar=6*maxres)
2300 end subroutine do_mcm
2302 !-----------------------------------------------------------------------------
2304 subroutine do_mcm(i_orig)
2306 ! Monte-Carlo-with-Minimization calculations - parallel code.
2308 use control_data, only:refstr!,tag
2309 use io_base, only:intout,briefout
2310 use control, only:ovrtim,tcpu
2311 use compare, only:contact,contact_fract
2312 use minimm, only:minimize
2313 use regularize_, only:fitsq
2314 ! implicit real*8 (a-h,o-z)
2315 ! include 'DIMENSIONS'
2317 ! include 'COMMON.IOUNITS'
2318 ! include 'COMMON.GEO'
2319 ! include 'COMMON.CHAIN'
2320 ! include 'COMMON.MCM'
2321 ! include 'COMMON.CONTACTS'
2322 ! include 'COMMON.CONTROL'
2323 ! include 'COMMON.VAR'
2324 ! include 'COMMON.INTERACT'
2325 ! include 'COMMON.INFO'
2326 ! include 'COMMON.CACHE'
2327 !rc include 'COMMON.DEFORM'
2328 !rc include 'COMMON.DEFORM1'
2329 !rc include 'COMMON.DEFORM2'
2330 ! include 'COMMON.MINIM'
2331 ! include 'COMMON.NAMES'
2332 logical :: accepted,over,error,lprint,not_done,similar,&
2333 enelower,non_conv,flag,finish !,ovrtim
2334 integer :: MoveType,nbond !,conf_comp
2335 real(kind=8),dimension(6*nres) :: x1,varold1,varold,varia !(maxvar) (maxvar=6*maxres)
2336 real(kind=8) :: elowest,eold
2337 real(kind=8) :: przes(3),obr(3,3)
2338 integer :: iparentx(max_threadss2)
2339 integer :: iparentx1(max_threadss2)
2340 integer :: imtasks(150),imtasks_n
2341 real(kind=8) :: energia(0:n_ene)
2344 integer :: nfun,nodenum,i_orig,i,nf,nacc,it,nout,j,kkk,is,&
2345 Kwita,iretcode,noverlap,nstart_grow,ierr,iitt,&
2346 ii_grnum_d,ii_ennum_d,ii_hesnum_d,i_grnum_d,i_ennum_d,&
2347 i_hesnum_d,i_minimiz,irepet
2348 real(kind=8) :: etot,frac,eneglobal,RandOrPert,eold1,co,&
2351 ! if(.not.allocated(varsave)) allocate(varsave(maxvar,maxsave)) !(maxvar,maxsave)
2352 print *,'Master entered DO_MCM'
2360 !---------------------------------------------------------------------------
2361 ! Initialize counters.
2362 !---------------------------------------------------------------------------
2363 ! Total number of generated confs.
2365 ! Total number of moves. In general this won`t be equal to the number of
2366 ! attempted moves, because we may want to reject some "bad" confs just by
2369 ! Total number of temperature jumps.
2371 ! Total number of shift (nbond_move(1)), spike, crankshaft, three-bond,...
2373 allocate(nbond_move(nres)) !(maxres)
2379 ! Initialize total and accepted number of moves of various kind.
2384 ! Total number of energy evaluations.
2388 ! write (iout,*) 'RanFract=',RanFract
2391 !----------------------------------------------------------------------------
2392 ! Compute and print initial energies.
2393 !----------------------------------------------------------------------------
2395 write (iout,'(/80(1h*)/a)') 'Initial energies:'
2398 call etotal(energia)
2400 call enerprint(energia)
2401 ! Request energy computation from slave processors.
2402 call geom_to_var(nvar,varia)
2403 ! write (iout,*) 'The VARIA array'
2404 ! write (iout,'(8f10.4)') (rad2deg*varia(i),i=1,nvar)
2405 call minimize(etot,varia,iretcode,nfun)
2406 call var_to_geom(nvar,varia)
2408 write (iout,*) 'etot from MINIMIZE:',etot
2409 ! write (iout,*) 'Tha VARIA array'
2410 ! write (iout,'(8f10.4)') (rad2deg*varia(i),i=1,nvar)
2413 if (print_mc .gt. 0) write (iout,'(/80(1h*)/20x,a/80(1h*))') &
2414 'Enter Monte Carlo procedure.'
2415 if (print_mc .gt. 0) write (iout,'(i5,1pe14.5)' ) i_orig,etot
2421 call zapis(varia,etot)
2423 call var_to_geom(nvar,varia)
2425 call etotal(energia)
2426 if (print_mc.gt.0) write (iout,*) 'Initial energy:',etot
2428 nacc=0 ! total # of accepted confs of the current processor.
2429 nacc_tot=0 ! total # of accepted confs of all processors.
2431 !----------------------------------------------------------------------------
2433 !----------------------------------------------------------------------------
2436 LOOP1:do while (not_done)
2438 if (print_mc.gt.0) write (iout,'(80(1h*)/20x,a,i7)') &
2439 'Beginning iteration #',it
2440 ! Initialize local counter.
2441 ntrial=0 ! # of generated non-overlapping confs.
2442 noverlap=0 ! # of overlapping confs.
2444 LOOP2:do while (.not. accepted)
2446 LOOP3:do while (imtasks_n.lt.nodenum-1.and..not.finish)
2448 if(imtasks(i).eq.0) then
2453 ! Retrieve the angles of previously accepted conformation
2457 call var_to_geom(nvar,varia)
2458 ! Rebuild the chain.
2460 ! Heat up the system, if necessary.
2462 ! If temperature cannot be further increased, stop.
2468 ! write (iout,'(a)') 'Old variables:'
2469 ! write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
2470 ! Decide whether to generate a random conformation or perturb the old one
2471 RandOrPert=ran_number(0.0D0,1.0D0)
2472 if (RandOrPert.gt.RanFract) then
2473 if (print_mc.gt.0) &
2474 write (iout,'(a)') 'Perturbation-generated conformation.'
2475 call perturb(error,lprint,MoveType,nbond,1.0D0)
2476 ! print *,'after perturb',error,finish
2477 if (error) finish = .true.
2478 if (MoveType.lt.1 .or. MoveType.gt.MaxMoveType) then
2479 write (iout,'(/a,i7,a/)') 'Error - unknown MoveType=',&
2480 MoveType,' returned from PERTURB.'
2482 write (*,'(/a,i7,a/)') 'Error - unknown MoveType=',&
2483 MoveType,' returned from PERTURB.'
2489 nstart_grow=iran_num(3,nres)
2490 if (print_mc.gt.0) &
2491 write (iout,'(2a,i3)') 'Random-generated conformation',&
2492 ' - chain regrown from residue',nstart_grow
2493 call gen_rand_conf(nstart_grow,*30)
2495 call geom_to_var(nvar,varia)
2497 ! print *,'finish=',finish
2498 if (etot-elowest.gt.overlap_cut) then
2499 if (print_mc.gt.1) write (iout,'(a,1pe14.5)') &
2500 'Overlap detected in the current conf.; energy is',etot
2501 if(iprint.gt.1.or.etot.lt.1d19) print *,&
2502 'Overlap detected in the current conf.; energy is',etot
2506 if (noverlap.gt.maxoverlap) then
2507 write (iout,*) 'Too many overlapping confs.',&
2508 ' etot, elowest, overlap_cut', etot, elowest, overlap_cut
2511 else if (.not. finish) then
2512 ! Distribute tasks to processors
2513 ! print *,'Master sending order'
2514 call MPI_SEND(12, 1, MPI_INTEGER, is, tag,&
2516 ! write (iout,*) '12: tag=',tag
2517 ! print *,'Master sent order to processor',is
2518 call MPI_SEND(it, 1, MPI_INTEGER, is, tag,&
2520 ! write (iout,*) 'it: tag=',tag
2521 call MPI_SEND(eold, 1, MPI_DOUBLE_PRECISION, is, tag,&
2523 ! write (iout,*) 'eold: tag=',tag
2524 call MPI_SEND(varia(1), nvar, MPI_DOUBLE_PRECISION, &
2527 ! write (iout,*) 'varia: tag=',tag
2528 call MPI_SEND(varold(1), nvar, MPI_DOUBLE_PRECISION, &
2531 ! write (iout,*) 'varold: tag=',tag
2538 imtasks_n=imtasks_n+1
2544 LOOP_RECV:do while(.not.flag)
2546 call MPI_IPROBE(is,tag,CG_COMM,flag,status,ierr)
2548 call MPI_RECV(iitt, 1, MPI_INTEGER, is, tag,&
2549 CG_COMM, status, ierr)
2550 call MPI_RECV(eold1, 1, MPI_DOUBLE_PRECISION, is, tag,&
2551 CG_COMM, status, ierr)
2552 call MPI_RECV(etot, 1, MPI_DOUBLE_PRECISION, is, tag,&
2553 CG_COMM, status, ierr)
2554 call MPI_RECV(varia(1), nvar, MPI_DOUBLE_PRECISION,is,tag,&
2555 CG_COMM, status, ierr)
2556 call MPI_RECV(varold1(1), nvar, MPI_DOUBLE_PRECISION, is, &
2557 tag, CG_COMM, status, ierr)
2558 call MPI_RECV(ii_grnum_d, 1, MPI_INTEGER, is, tag,&
2559 CG_COMM, status, ierr)
2560 call MPI_RECV(ii_ennum_d, 1, MPI_INTEGER, is, tag,&
2561 CG_COMM, status, ierr)
2562 call MPI_RECV(ii_hesnum_d, 1, MPI_INTEGER, is, tag,&
2563 CG_COMM, status, ierr)
2564 i_grnum_d=i_grnum_d+ii_grnum_d
2565 i_ennum_d=i_ennum_d+ii_ennum_d
2566 neneval = neneval+ii_ennum_d
2567 i_hesnum_d=i_hesnum_d+ii_hesnum_d
2568 i_minimiz=i_minimiz+1
2570 imtasks_n=imtasks_n-1
2576 if(print_mc.gt.0) write (iout,'(a,i6,a,i6,a,i6,a,1pe16.6)') &
2577 'From Worker #',is,' iitt',iitt,&
2578 ' Conformation:',ngen,' energy:',etot
2579 !--------------------------------------------------------------------------
2580 !... Do Metropolis test
2581 !--------------------------------------------------------------------------
2582 call metropolis(nvar,varia,varold1,etot,eold1,accepted,&
2584 if(iitt.ne.it.and..not.similar) then
2585 call metropolis(nvar,varia,varold,etot,eold,accepted,&
2589 if(etot.lt.eneglobal)eneglobal=etot
2590 ! if(mod(it,100).eq.0)
2591 write(iout,*)'CHUJOJEB ',neneval,eneglobal
2593 ! Write the accepted conformation.
2596 call var_to_geom(nvar,varia)
2599 call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,kkk),&
2600 nsup,przes,obr,non_conv)
2602 call contact(.false.,ncont,icont,co)
2603 frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
2604 write (iout,'(a,f8.3,a,f8.3,a,f8.3)') &
2605 'RMS deviation from the reference structure:',rms,&
2606 ' % of native contacts:',frac*100,' contact order:',co
2608 if (print_mc.gt.0) &
2609 write (iout,*) 'Writing new conformation',nout
2610 if (print_stat) then
2611 call var_to_geom(nvar,varia)
2612 #if defined(AIX) || defined(PGI)
2613 open (istat,file=statname,position='append')
2615 open (istat,file=statname,access='append')
2618 write (istat,'(i5,16(1pe14.5))') nout,&
2619 (energia(print_order(i)),i=1,nprint_ene),&
2622 write (istat,'(i5,16(1pe14.5))') nout,&
2623 (energia(print_order(i)),i=1,nprint_ene),etot
2627 ! Print internal coordinates.
2628 if (print_int) call briefout(nout,etot)
2631 if (elowest.gt.etot) elowest=etot
2632 moves_acc(MoveType)=moves_acc(MoveType)+1
2633 if (MoveType.eq.1) then
2634 nbond_acc(nbond)=nbond_acc(nbond)+1
2636 ! Check against conformation repetitions.
2637 irepet=conf_comp(varia,etot)
2638 if (nrepm.gt.maxrepm) then
2639 if (print_mc.gt.0) &
2640 write (iout,'(a)') 'Too many conformation repetitions.'
2643 ! Store the accepted conf. and its energy.
2648 if (irepet.eq.0) call zapis(varia,etot)
2649 ! Lower the temperature, if necessary.
2655 if(finish.and.imtasks_n.eq.0)exit LOOP2
2656 enddo LOOP2 ! accepted
2657 ! Check for time limit.
2658 not_done = (it.lt.max_mcm_it) .and. (nacc_tot.lt.maxacc)
2659 if(.not.not_done .or. finish) then
2660 if(imtasks_n.gt.0) then
2667 enddo LOOP1 ! not_done
2669 if (print_mc.gt.0) then
2670 write (iout,'(/80(1h*)/20x,a)') 'Summary run statistics:'
2671 call statprint(nacc,nfun,iretcode,etot,elowest)
2672 write (iout,'(a)') &
2673 'Statistics of multiple-bond motions. Total motions:'
2674 write (iout,'(16i5)') (nbond_move(i),i=1,Nbm)
2675 write (iout,'(a)') 'Accepted motions:'
2676 write (iout,'(16i5)') (nbond_acc(i),i=1,Nbm)
2678 write (iout,'(/80(1h*)/20x,a/80(1h*)/)') 'All iterations done.'
2686 call MPI_SEND(999, 1, MPI_INTEGER, is, tag,&
2690 end subroutine do_mcm
2691 !-----------------------------------------------------------------------------
2692 subroutine execute_slave(nodeinfo,iprint)
2695 use minimm, only:minimize
2696 ! implicit real*8 (a-h,o-z)
2697 ! include 'DIMENSIONS'
2699 ! include 'COMMON.TIME1'
2700 ! include 'COMMON.IOUNITS'
2701 !rc include 'COMMON.DEFORM'
2702 !rc include 'COMMON.DEFORM1'
2703 !rc include 'COMMON.DEFORM2'
2704 ! include 'COMMON.LOCAL'
2705 ! include 'COMMON.VAR'
2706 ! include 'COMMON.INFO'
2707 ! include 'COMMON.MINIM'
2708 character(len=10) :: nodeinfo
2709 real(kind=8),dimension(6*nres) :: x,x1 !(maxvar) (maxvar=6*maxres)
2710 integer :: nfun,iprint,i_switch,ierr,i_grnum_d,i_ennum_d,&
2711 i_hesnum_d,iitt,iretcode,iminrep
2712 real(kind=8) :: ener,energyx
2714 nodeinfo='chujwdupe'
2715 ! print *,'Processor:',MyID,' Entering execute_slave'
2717 ! call MPI_SEND(nodeinfo, 10, MPI_CHARACTER, 0, tag,
2720 1001 call MPI_RECV(i_switch, 1, MPI_INTEGER, 0, tag,&
2721 CG_COMM, status, ierr)
2722 ! write(iout,*)'12: tag=',tag
2723 if(iprint.ge.2)print *, MyID,' recv order ',i_switch
2724 if (i_switch.eq.12) then
2728 call MPI_RECV(iitt, 1, MPI_INTEGER, 0, tag,&
2729 CG_COMM, status, ierr)
2730 ! write(iout,*)'12: tag=',tag
2731 call MPI_RECV(ener, 1, MPI_DOUBLE_PRECISION, 0, tag,&
2732 CG_COMM, status, ierr)
2733 ! write(iout,*)'ener: tag=',tag
2734 call MPI_RECV(x(1), nvar, MPI_DOUBLE_PRECISION, 0, tag,&
2735 CG_COMM, status, ierr)
2736 ! write(iout,*)'x: tag=',tag
2737 call MPI_RECV(x1(1), nvar, MPI_DOUBLE_PRECISION, 0, tag,&
2738 CG_COMM, status, ierr)
2739 ! write(iout,*)'x1: tag=',tag
2745 ! print *,'calling minimize'
2746 call minimize(energyx,x,iretcode,nfun)
2748 write(iout,100)'minimized energy = ',energyx,&
2749 ' # funeval:',nfun,' iret ',iretcode
2750 write(*,100)'minimized energy = ',energyx,&
2751 ' # funeval:',nfun,' iret ',iretcode
2752 100 format(a20,f10.5,a12,i5,a6,i2)
2753 if(iretcode.eq.10) then
2756 write(iout,*)' ... not converged - trying again ',iminrep
2757 call minimize(energyx,x,iretcode,nfun)
2759 write(iout,*)'minimized energy = ',energyx,&
2760 ' # funeval:',nfun,' iret ',iretcode
2761 if(iretcode.ne.10)go to 812
2763 if(iretcode.eq.10) then
2765 write(iout,*)' ... not converged again - giving up'
2770 ! print *,'Sending results'
2771 call MPI_SEND(iitt, 1, MPI_INTEGER, 0, tag,&
2773 call MPI_SEND(ener, 1, MPI_DOUBLE_PRECISION, 0, tag,&
2775 call MPI_SEND(energyx, 1, MPI_DOUBLE_PRECISION, 0, tag,&
2777 call MPI_SEND(x(1), nvar, MPI_DOUBLE_PRECISION, 0, tag,&
2779 call MPI_SEND(x1(1), nvar, MPI_DOUBLE_PRECISION, 0, tag,&
2781 call MPI_SEND(i_grnum_d, 1, MPI_INTEGER, 0, tag,&
2783 call MPI_SEND(nfun, 1, MPI_INTEGER, 0, tag,&
2785 call MPI_SEND(i_hesnum_d, 1, MPI_INTEGER, 0, tag,&
2787 ! print *,'End sending'
2792 end subroutine execute_slave
2794 !-----------------------------------------------------------------------------
2795 subroutine heat(over)
2797 ! implicit real*8 (a-h,o-z)
2798 ! include 'DIMENSIONS'
2799 ! include 'COMMON.MCM'
2800 ! include 'COMMON.IOUNITS'
2802 ! Check if there`s a need to increase temperature.
2803 if (ntrial.gt.maxtrial) then
2804 if (NstepH.gt.0) then
2805 if (dabs(Tcur-TMax).lt.1.0D-7) then
2806 if (print_mc.gt.0) &
2807 write (iout,'(/80(1h*)/a,f8.3,a/80(1h*))') &
2808 'Upper limit of temperature reached. Terminating.'
2813 if (Tcur.gt.Tmax) Tcur=Tmax
2814 betbol=1.0D0/(Rbol*Tcur)
2815 if (print_mc.gt.0) &
2816 write (iout,'(/80(1h*)/a,f8.3,a,f10.5/80(1h*))') &
2817 'System heated up to ',Tcur,' K; BetBol:',betbol
2823 if (print_mc.gt.0) &
2824 write (iout,'(a)') &
2825 'Maximum number of trials in a single MCM iteration exceeded.'
2834 !-----------------------------------------------------------------------------
2837 ! implicit real*8 (a-h,o-z)
2838 ! include 'DIMENSIONS'
2839 ! include 'COMMON.MCM'
2840 ! include 'COMMON.IOUNITS'
2841 if (nstepC.gt.0 .and. dabs(Tcur-Tmin).gt.1.0D-7) then
2843 if (Tcur.lt.Tmin) Tcur=Tmin
2844 betbol=1.0D0/(Rbol*Tcur)
2845 if (print_mc.gt.0) &
2846 write (iout,'(/80(1h*)/a,f8.3,a,f10.5/80(1h*))') &
2847 'System cooled down up to ',Tcur,' K; BetBol:',betbol
2851 !-----------------------------------------------------------------------------
2852 subroutine perturb(error,lprint,MoveType,nbond,max_phi)
2855 use energy, only:nnt,nct,itype
2856 use md_calc, only:bond_move
2857 ! implicit real*8 (a-h,o-z)
2858 ! include 'DIMENSIONS'
2859 integer,parameter :: MMaxSideMove=100
2860 ! include 'COMMON.MCM'
2861 ! include 'COMMON.CHAIN'
2862 ! include 'COMMON.INTERACT'
2863 ! include 'COMMON.VAR'
2864 ! include 'COMMON.GEO'
2865 ! include 'COMMON.NAMES'
2866 ! include 'COMMON.IOUNITS'
2867 !rc include 'COMMON.DEFORM1'
2868 logical :: error,lprint,fail
2869 integer :: MoveType,nbond,end_select,ind_side(MMaxSideMove)
2870 real(kind=8) :: max_phi
2871 real(kind=8) :: psi!,gen_psi
2872 !el external iran_num
2873 !el integer iran_num
2876 !!! local variables - el
2877 integer :: itrial,iiwin,iwindow,isctry,i,icount,j,nstart,&
2878 nside_move,inds,indx,ii,iti
2879 real(kind=8) :: bond_prob,theta_new
2884 ! Perturb the conformation according to a randomly selected move.
2885 call SelectMove(MoveType)
2886 ! write (iout,*) 'MoveType=',MoveType
2888 goto (100,200,300,400,500) MoveType
2889 !------------------------------------------------------------------------------
2890 ! Backbone N-bond move.
2891 ! Select the number of bonds (length of the segment to perturb).
2893 if (itrial.gt.1000) then
2894 write (iout,'(a)') 'Too many attempts at multiple-bond move.'
2898 bond_prob=ran_number(0.0D0,sumpro_bond(nbm))
2899 ! print *,'sumpro_bond(nbm)=',sumpro_bond(nbm),
2900 ! & ' Bond_prob=',Bond_Prob
2902 ! print *,i,Bond_Prob,sumpro_bond(i),sumpro_bond(i+1)
2903 if (bond_prob.ge.sumpro_bond(i) .and. &
2904 bond_prob.le.sumpro_bond(i+1)) then
2909 write (iout,'(2a)') 'In PERTURB: Error - number of bonds',&
2910 ' to move out of range.'
2914 if (nwindow.gt.0) then
2915 ! Select the first residue to perturb
2916 iwindow=iran_num(1,nwindow)
2917 print *,'iwindow=',iwindow
2919 do while (winlen(iwindow).lt.nbond)
2920 iwindow=iran_num(1,nwindow)
2922 if (iiwin.gt.1000) then
2923 write (iout,'(a)') 'Cannot select moveable residues.'
2928 nstart=iran_num(winstart(iwindow),winend(iwindow))
2930 nstart = iran_num(koniecl+2,nres-nbond-koniecl)
2931 !d print *,'nres=',nres,' nbond=',nbond,' koniecl=',koniecl,
2932 !d & ' nstart=',nstart
2935 if (psi.eq.0.0) then
2939 if (print_mc.gt.1) write (iout,'(a,i4,a,i4,a,f8.3)') &
2940 'PERTURB: nbond=',nbond,' nstart=',nstart,' psi=',psi*rad2deg
2941 !d print *,'nstart=',nstart
2942 call bond_move(nbond,nstart,psi,.false.,error)
2944 write (iout,'(2a)') &
2945 'Could not define reference system in bond_move, ',&
2946 'choosing ahother segment.'
2950 nbond_move(nbond)=nbond_move(nbond)+1
2954 !------------------------------------------------------------------------------
2955 ! Backbone endmove. Perturb a SINGLE angle of a residue close to the end of
2959 ! end_select=iran_num(1,2*koniecl)
2960 ! if (end_select.gt.koniecl) then
2961 ! end_select=nphi-(end_select-koniecl)
2963 ! end_select=koniecl+3
2965 ! if (nwindow.gt.0) then
2966 ! iwin=iran_num(1,nwindow)
2967 ! i1=max0(4,winstart(iwin))
2968 ! i2=min0(winend(imin)+2,nres)
2969 ! end_select=iran_num(i1,i2)
2971 ! iselect = iran_num(1,nmov_var)
2974 ! if (isearch_tab(i).eq.1) jj = jj+1
2975 ! if (jj.eq.iselect) then
2981 end_select = iran_num(4,nres)
2982 psi=max_phi*gen_psi()
2983 if (psi.eq.0.0D0) then
2987 phi(end_select)=pinorm(phi(end_select)+psi)
2988 if (print_mc.gt.1) write (iout,'(a,i4,a,f8.3,a,f8.3)') &
2989 'End angle',end_select,' moved by ',psi*rad2deg,' new angle:',&
2990 phi(end_select)*rad2deg
2991 ! if (end_select.gt.3)
2992 ! & theta(end_select-1)=gen_theta(itype(end_select-2),
2993 ! & phi(end_select-1),phi(end_select))
2994 ! if (end_select.lt.nres)
2995 ! & theta(end_select)=gen_theta(itype(end_select-1),
2996 ! & phi(end_select),phi(end_select+1))
2997 !d print *,'nres=',nres,' end_select=',end_select
2998 !d print *,'theta',end_select-1,theta(end_select-1)
2999 !d print *,'theta',end_select,theta(end_select)
3004 !------------------------------------------------------------------------------
3006 ! Select the number of SCs to perturb.
3008 301 nside_move=iran_num(1,MaxSideMove)
3009 ! print *,'nside_move=',nside_move,' MaxSideMove',MaxSideMove
3010 ! Select the indices.
3013 111 inds=iran_num(nnt,nct)
3015 if (icount.gt.1000) then
3016 write (iout,'(a)')'Error - cannot select side chains to move.'
3020 if (itype(inds).eq.10) goto 111
3022 if (inds.eq.ind_side(j)) goto 111
3025 if (inds.lt.ind_side(j)) then
3031 112 do j=i,indx+1,-1
3032 ind_side(j)=ind_side(j-1)
3034 113 ind_side(indx)=inds
3036 ! Carry out perturbation.
3040 call gen_side(iti,theta(ii+1),alph(ii),omeg(ii),fail)
3043 if (isctry.gt.1000) then
3044 write (iout,'(a)') 'Too many errors in SC generation.'
3050 if (print_mc.gt.1) write (iout,'(2a,i4,a,2f8.3)') &
3051 'Side chain ',restyp(iti),ii,' moved to ',&
3052 alph(ii)*rad2deg,omeg(ii)*rad2deg
3057 !------------------------------------------------------------------------------
3059 400 end_select=iran_num(3,nres)
3060 theta_new=gen_theta(itype(end_select),phi(end_select),&
3062 if (print_mc.gt.1) write (iout,'(a,i3,a,f8.3,a,f8.3)') &
3063 'Theta ',end_select,' moved from',theta(end_select)*rad2deg,&
3064 ' to ',theta_new*rad2deg
3065 theta(end_select)=theta_new
3069 !------------------------------------------------------------------------------
3070 ! Error returned from SelectMove.
3073 end subroutine perturb
3074 !-----------------------------------------------------------------------------
3075 subroutine SelectMove(MoveType)
3077 ! implicit real*8 (a-h,o-z)
3078 ! include 'DIMENSIONS'
3079 ! include 'COMMON.MCM'
3080 ! include 'COMMON.IOUNITS'
3082 !!! local variables - el
3083 integer :: i,MoveType
3084 real(kind=8) :: what_move
3086 what_move=ran_number(0.0D0,sumpro_type(MaxMoveType))
3088 if (what_move.ge.sumpro_type(i-1).and. &
3089 what_move.lt.sumpro_type(i)) then
3094 write (iout,'(a)') &
3095 'Fatal error in SelectMoveType: cannot select move.'
3096 MoveType=MaxMoveType+1
3098 end subroutine SelectMove
3099 !-----------------------------------------------------------------------------
3100 real(kind=8) function gen_psi()
3102 use geometry_data, only: angmin,pi
3105 real(kind=8) :: x !,ran_number
3106 ! include 'COMMON.IOUNITS'
3107 ! include 'COMMON.GEO'
3110 x=ran_number(-pi,pi)
3111 if (dabs(x).gt.angmin) then
3116 write (iout,'(a)')'From Gen_Psi: Cannot generate angle increment.'
3119 end function gen_psi
3120 !-----------------------------------------------------------------------------
3121 subroutine metropolis(n,xcur,xold,ecur,eold,accepted,similar,enelower)
3123 ! implicit real*8 (a-h,o-z)
3124 ! include 'DIMENSIONS'
3125 ! include 'COMMON.MCM'
3126 ! include 'COMMON.IOUNITS'
3127 ! include 'COMMON.VAR'
3128 ! include 'COMMON.GEO'
3129 !rc include 'COMMON.DEFORM'
3131 real(kind=8) :: ecur,eold,xx,bol !,ran_number
3132 real(kind=8),dimension(n) :: xcur,xold
3133 real(kind=8) :: ecut1 ,ecut2 ,tola
3134 logical :: accepted,similar,not_done,enelower
3137 !!! local variables - el
3138 real(kind=8) :: xxh,difene,reldife
3140 data ecut1 /-1.0D-5/,ecut2 /5.0D-3/,tola/5.0D0/
3144 ! Set lprn=.true. for debugging.
3147 !el write(iout,*)'enedif',enedif,' ecut1',ecut1,' ecut2',ecut2
3148 write(iout,*)' ecut1',ecut1,' ecut2',ecut2
3152 ! Check if the conformation is similar.
3154 reldife=difene/dmax1(dabs(eold),dabs(ecur),1.0D0)
3156 write (iout,*) 'Metropolis'
3157 write(iout,*)'ecur,eold,difene,reldife',ecur,eold,difene,reldife
3159 ! If energy went down remarkably, we accept the new conformation
3161 !jp if (reldife.lt.ecut1) then
3162 if (difene.lt.ecut1) then
3165 if (lprn) write (iout,'(a)') &
3166 'Conformation accepted, because energy has lowered remarkably.'
3167 ! elseif (reldife.lt.ecut2 .and. dif_ang(nphi,xcur,xold).lt.tola)
3168 !jp elseif (reldife.lt.ecut2)
3169 elseif (difene.lt.ecut2) &
3171 ! Reject the conf. if energy has changed insignificantly and there is not
3172 ! much change in conformation.
3174 write (iout,'(2a)') 'Conformation rejected, because it is',&
3175 ' similar to the preceding one.'
3179 ! Else carry out Metropolis test.
3181 xx=ran_number(0.0D0,1.0D0)
3184 write (iout,*) 'betbol=',betbol,' difene=',difene,' xxh=',xxh
3185 if (xxh.gt.50.0D0) then
3190 if (lprn) write (iout,*) 'bol=',bol,' xx=',xx
3193 if (lprn) write (iout,'(a)') &
3194 'Conformation accepted, because it passed Metropolis test.'
3197 if (lprn) write (iout,'(a)') &
3198 'Conformation rejected, because it did not pass Metropolis test.'
3207 end subroutine metropolis
3208 !-----------------------------------------------------------------------------
3209 integer function conf_comp(x,ene)
3211 use geometry_data, only: nphi
3212 ! implicit real*8 (a-h,o-z)
3213 ! include 'DIMENSIONS'
3214 ! include 'COMMON.MCM'
3215 ! include 'COMMON.VAR'
3216 ! include 'COMMON.IOUNITS'
3217 ! include 'COMMON.GEO'
3218 real(kind=8) :: etol, angtol
3219 real(kind=8),dimension(6*nres) :: x !(maxvar) (maxvar=6*maxres)
3220 real(kind=8) :: difa !dif_ang,
3222 !!! local variables - el
3226 data etol /0.1D0/, angtol /20.0D0/
3228 ! write (iout,*) 'ii=',ii,'ene=',ene,esave(ii),dabs(ene-esave(ii))
3229 if (dabs(ene-esave(ii)).lt.etol) then
3230 difa=dif_ang(nphi,x,varsave(1,ii))
3232 ! write(iout,'(i3,3f8.3)')i,rad2deg*x(i),
3233 ! & rad2deg*varsave(i,ii)
3235 ! write(iout,*) 'ii=',ii,' difa=',difa,' angtol=',angtol
3236 if (difa.le.angtol) then
3237 if (print_mc.gt.0) then
3238 write (iout,'(a,i5,2(a,1pe15.4))') &
3239 'Current conformation matches #',ii,&
3240 ' in the store array ene=',ene,' esave=',esave(ii)
3241 ! write (*,'(a,i5,a)') 'Current conformation matches #',ii,
3242 ! & ' in the store array.'
3243 endif ! print_mc.gt.0
3244 if (print_mc.gt.1) then
3246 write(iout,'(i3,3f8.3)')i,rad2deg*x(i),&
3247 rad2deg*varsave(i,ii)
3249 endif ! print_mc.gt.1
3258 end function conf_comp
3259 !-----------------------------------------------------------------------------
3260 real(kind=8) function dif_ang(n,x,y)
3262 use geometry_data, only: dwapi
3265 real(kind=8),dimension(n) :: x,y
3266 real(kind=8) :: w,wa,dif,difa
3267 !el real(kind=8) :: pinorm
3268 ! include 'COMMON.GEO'
3272 dif=dabs(pinorm(y(i)-x(i)))
3273 if (dabs(dif-dwapi).lt.dif) dif=dabs(dif-dwapi)
3274 w=1.0D0-(2.0D0*(i-1)/(n-1)-1.0D0)**2+1.0D0/n
3278 dif_ang=rad2deg*dsqrt(difa/wa)
3280 end function dif_ang
3281 !-----------------------------------------------------------------------------
3282 subroutine add2cache(n1,n2,ncache,nvar,SourceID,CachSrc,ecur,xcur,ecache,xcache)
3285 ! include 'COMMON.GEO'
3286 ! include 'COMMON.IOUNITS'
3287 integer :: n1,n2,ncache,nvar,SourceID,CachSrc(n2)
3289 real(kind=8) :: ecur,xcur(nvar),ecache(n2),xcache(n1,n2)
3290 !d write (iout,*) 'Enter ADD2CACHE ncache=',ncache ,' ecur',ecur
3291 !d write (iout,'(10f8.3)') (rad2deg*xcur(i),i=1,nvar)
3292 !d write (iout,*) 'Old CACHE array:'
3294 !d write (iout,*) 'i=',i,' ecache=',ecache(i),' CachSrc',CachSrc(i)
3295 !d write (iout,'(10f8.3)') (rad2deg*xcache(j,i),j=1,nvar)
3299 do while (i.gt.0 .and. ecur.lt.ecache(i))
3303 !d write (iout,*) 'i=',i,' ncache=',ncache
3304 if (ncache.eq.n2) then
3305 write (iout,*) 'Cache dimension exceeded',ncache,n2
3306 write (iout,*) 'Highest-energy conformation will be removed.'
3310 ecache(ii+1)=ecache(ii)
3311 CachSrc(ii+1)=CachSrc(ii)
3313 xcache(j,ii+1)=xcache(j,ii)
3322 !d write (iout,*) 'New CACHE array:'
3324 !d write (iout,*) 'i=',i,' ecache=',ecache(i),' CachSrc',CachSrc(i)
3325 !d write (iout,'(10f8.3)') (rad2deg*xcache(j,i),j=1,nvar)
3328 end subroutine add2cache
3329 !-----------------------------------------------------------------------------
3330 subroutine rm_from_cache(i,n1,n2,ncache,nvar,CachSrc,ecache,xcache)
3333 ! include 'COMMON.GEO'
3334 ! include 'COMMON.IOUNITS'
3335 integer :: n1,n2,ncache,nvar,CachSrc(n2)
3337 real(kind=8) :: ecache(n2),xcache(n1,n2)
3339 !d write (iout,*) 'Enter RM_FROM_CACHE'
3340 !d write (iout,*) 'Old CACHE array:'
3342 !d write (iout,*)'i=',ii,' ecache=',ecache(ii),' CachSrc',CachSrc(ii)
3343 !d write (iout,'(10f8.3)') (rad2deg*xcache(j,ii),j=1,nvar)
3347 ecache(ii-1)=ecache(ii)
3348 CachSrc(ii-1)=CachSrc(ii)
3350 xcache(j,ii-1)=xcache(j,ii)
3354 !d write (iout,*) 'New CACHE array:'
3356 !d write (iout,*) 'i=',i,' ecache=',ecache(i),' CachSrc',CachSrc(i)
3357 !d write (iout,'(10f8.3)') (rad2deg*xcache(j,i),j=1,nvar)
3360 end subroutine rm_from_cache
3361 !-----------------------------------------------------------------------------
3363 !-----------------------------------------------------------------------------
3364 subroutine statprint(it,nfun,iretcode,etot,elowest)
3366 use control_data, only: MaxMoveType,minim
3367 use control, only: tcpu
3369 ! implicit real*8 (a-h,o-z)
3370 ! include 'DIMENSIONS'
3371 ! include 'COMMON.IOUNITS'
3372 ! include 'COMMON.CONTROL'
3373 ! include 'COMMON.MCM'
3375 integer :: it,nfun,iretcode,i
3376 real(kind=8) :: etot,elowest,fr_mov_i
3380 '(80(1h*)/a,i5,a,1pe14.5,a,1pe14.5/a,i3,a,i10,a,i5,a,i5)') &
3381 'Finished iteration #',it,' energy is',etot,&
3382 ' lowest energy:',elowest,&
3383 'SUMSL return code:',iretcode,&
3384 ' # of energy evaluations:',neneval,&
3385 '# of temperature jumps:',ntherm,&
3386 ' # of minima repetitions:',nrepm
3388 write (iout,'(80(1h*)/a,i8,a,1pe14.5,a,1pe14.5)') &
3389 'Finished iteration #',it,' energy is',etot,&
3390 ' lowest energy:',elowest
3392 write (iout,'(/4a)') &
3393 'Kind of move ',' total',' accepted',&
3395 write (iout,'(58(1h-))')
3397 if (moves(i).eq.0) then
3400 fr_mov_i=dfloat(moves_acc(i))/dfloat(moves(i))
3402 write(iout,'(a,2i15,f10.5)')MovTypID(i),moves(i),moves_acc(i),&
3405 write (iout,'(a,2i15,f10.5)') 'total ',nmove,nacc_tot,&
3406 dfloat(nacc_tot)/dfloat(nmove)
3407 write (iout,'(58(1h-))')
3408 write (iout,'(a,1pe12.4)') 'Elapsed time:',tcpu()
3410 end subroutine statprint
3411 !-----------------------------------------------------------------------------
3412 subroutine zapis(varia,etot)
3414 use geometry_data, only: nres,rad2deg,nvar
3417 ! implicit real*8 (a-h,o-z)
3418 ! include 'DIMENSIONS'
3420 use MPI_data !include 'COMMON.INFO'
3423 ! include 'COMMON.GEO'
3424 ! include 'COMMON.VAR'
3425 ! include 'COMMON.MCM'
3426 ! include 'COMMON.IOUNITS'
3427 integer,dimension(nsave) :: itemp
3428 real(kind=8),dimension(6*nres) :: varia !(maxvar) (maxvar=6*maxres)
3431 integer :: j,i,maxvar
3432 real(kind=8) :: etot
3434 !el allocate(esave(nsave)) !(maxsave)
3439 write (iout,'(a,i5,a,i5)') 'Enter ZAPIS NSave=',Nsave,&
3441 write (iout,'(a)') 'Current energy and conformation:'
3442 write (iout,'(1pe14.5)') etot
3443 write (iout,'(10f8.3)') (rad2deg*varia(i),i=1,nvar)
3445 ! Shift the contents of the esave and varsave arrays if filled up.
3446 call add2cache(6*nres,nsave,nsave,nvar,MyID,itemp,&
3447 etot,varia,esave,varsave)
3449 write (iout,'(a)') 'Energies and the VarSave array.'
3451 write (iout,'(i5,1pe14.5)') i,esave(i)
3452 write (iout,'(10f8.3)') (rad2deg*varsave(j,i),j=1,nvar)
3456 end subroutine zapis
3457 !-----------------------------------------------------------------------------
3458 !-----------------------------------------------------------------------------
3459 subroutine alloc_MCM_arrays
3461 use energy_data, only: max_ene
3465 allocate(entropy(-max_ene-4:max_ene)) !(-max_ene-4:max_ene)
3466 allocate(nhist(-max_ene:max_ene)) !(-max_ene:max_ene)
3467 allocate(nminima(maxsave)) !(maxsave)
3469 allocate(xpool(6*nres,max_pool)) !(maxvar,max_pool)(maxvar=6*maxres)
3470 allocate(epool(max_pool)) !(max_pool)
3473 if(.not.allocated(nsave_part)) allocate(nsave_part(nctasks)) !(max_cg_procs)
3476 ! real(kind=8),dimension(:),allocatable :: sumpro_type !(0:MaxMoveType)
3477 allocate(sumpro_bond(0:nres)) !(0:maxres)
3478 allocate(nbond_move(nres),nbond_acc(nres)) !(maxres)
3479 allocate(moves(-1:MaxMoveType+1),moves_acc(-1:MaxMoveType+1)) !(-1:MaxMoveType+1)
3480 ! common /accept_stats/
3481 ! allocate(nacc_part !(0:MaxProcs) !el nie uzywane???
3482 ! common /windows/ in io: mcmread
3483 ! allocate(winstart,winend,winlen !(maxres)
3485 !el allocate(MovTypID(-1:MaxMoveType+1)) !(-1:MaxMoveType+1)
3488 allocate(varsave(nres*6,maxsave)) !(maxvar,maxsave)(maxvar=6*maxres)
3489 allocate(esave(maxsave)) !(maxsave)
3490 allocate(Origin(maxsave)) !(maxsave)
3493 end subroutine alloc_MCM_arrays
3494 !-----------------------------------------------------------------------------
3495 !-----------------------------------------------------------------------------