Merge branch 'homology' of mmka.chem.univ.gda.pl:unres into homology
[unres.git] / source / unres / src_MD-M / mc.F
1       subroutine monte_carlo
2 C Does Boltzmann and entropic sampling without energy minimization
3       implicit real*8 (a-h,o-z)
4       include 'DIMENSIONS'
5       include 'COMMON.IOUNITS'
6 #ifdef MPL
7       include 'COMMON.INFO'
8 #endif
9       include 'COMMON.GEO'
10       include 'COMMON.CHAIN'
11       include 'COMMON.MCM'
12       include 'COMMON.MCE'
13       include 'COMMON.CONTACTS'
14       include 'COMMON.CONTROL'
15       include 'COMMON.VAR'
16       include 'COMMON.INTERACT'
17       include 'COMMON.THREAD'
18       include 'COMMON.NAMES'
19       logical accepted,not_done,over,ovrtim,error,lprint
20       integer MoveType,nbond,nbins
21       integer conf_comp
22       double precision RandOrPert
23       double precision varia(maxvar),elowest,elowest1,
24      &                 ehighest,ehighest1,eold
25       double precision przes(3),obr(3,3)
26       double precision varold(maxvar)
27       logical non_conv
28       integer moves1(-1:MaxMoveType+1,0:MaxProcs-1),
29      &        moves_acc1(-1:MaxMoveType+1,0:MaxProcs-1)
30 #ifdef MPL
31       double precision etot_temp,etot_all(0:MaxProcs)
32       external d_vadd,d_vmin,d_vmax
33       double precision entropy1(-max_ene:max_ene),
34      &                 nhist1(-max_ene:max_ene)
35       integer nbond_move1(maxres*(MaxProcs+1)),
36      &   nbond_acc1(maxres*(MaxProcs+1)),itemp(2)
37 #endif
38       double precision var_lowest(maxvar)
39       double precision energia(0:n_ene),energia_ave(0:n_ene)
40 C
41       write(iout,'(a,i8,2x,a,f10.5)')
42      & 'pool_read_freq=',pool_read_freq,' pool_fraction=',pool_fraction
43       open (istat,file=statname)
44       WhatsUp=0
45       indminn=-max_ene
46       indmaxx=max_ene
47       facee=1.0D0/(maxacc*delte)
48 C Number of bins in energy histogram
49       nbins=e_up/delte-1
50       write (iout,*) 'NBINS=',nbins
51       conste=dlog(facee)
52 C Read entropy from previous simulations. 
53       if (ent_read) then
54         read (ientin,*) indminn,indmaxx,emin,emax 
55         print *,'indminn=',indminn,' indmaxx=',indmaxx,' emin=',emin,
56      &          ' emax=',emax
57         do i=-max_ene,max_ene
58           entropy(i)=0.0D0
59         enddo
60         read (ientin,*) (ijunk,ejunk,entropy(i),i=indminn,indmaxx)
61         indmin=indminn
62         indmax=indmaxx
63         write (iout,*) 'indminn=',indminn,' indmaxx=',indmaxx,
64      &                 ' emin=',emin,' emax=',emax
65         write (iout,'(/a)') 'Initial entropy'
66         do i=indminn,indmaxx
67           write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
68         enddo
69       endif ! ent_read
70 C Read the pool of conformations
71       call read_pool
72       elowest=1.0D+10
73       ehighest=-1.0D+10
74 C----------------------------------------------------------------------------
75 C Entropy-sampling simulations with continually updated entropy;
76 C set NSWEEP=1 for Boltzmann sampling.
77 C Loop thru simulations
78 C----------------------------------------------------------------------------
79       DO ISWEEP=1,NSWEEP
80 C
81 C Initialize the IFINISH array.
82 C
83 #ifdef MPL
84         do i=1,nctasks
85           ifinish(i)=0
86         enddo
87 #endif
88 c---------------------------------------------------------------------------
89 C Initialize counters.
90 c---------------------------------------------------------------------------
91 C Total number of generated confs.
92         ngen=0
93 C Total number of moves. In general this won't be equal to the number of
94 C attempted moves, because we may want to reject some "bad" confs just by
95 C overlap check.
96         nmove=0
97 C Total number of shift (nbond_move(1)), spike, crankshaft, three-bond,...
98 C motions.
99         do i=1,nres
100           nbond_move(i)=0
101           nbond_acc(i)=0
102         enddo
103 C Initialize total and accepted number of moves of various kind.
104         do i=-1,MaxMoveType
105           moves(i)=0
106           moves_acc(i)=0
107         enddo
108 C Total number of energy evaluations.
109         neneval=0
110         nfun=0
111 C----------------------------------------------------------------------------
112 C Take a conformation from the pool
113 C----------------------------------------------------------------------------
114       rewind(istat)
115       write (iout,*) 'emin=',emin,' emax=',emax
116       if (npool.gt.0) then
117         ii=iran_num(1,npool)
118         do i=1,nvar
119           varia(i)=xpool(i,ii)
120         enddo
121         write (iout,*) 'Took conformation',ii,' from the pool energy=',
122      &               epool(ii)
123         call var_to_geom(nvar,varia)
124 C Print internal coordinates of the initial conformation
125         call intout
126       else if (isweep.gt.1) then
127         if (eold.lt.emax) then
128         do i=1,nvar
129           varia(i)=varold(i)
130         enddo
131         else
132         do i=1,nvar
133           varia(i)=var_lowest(i)
134         enddo
135         endif
136         call var_to_geom(nvar,varia)
137       endif
138 C----------------------------------------------------------------------------
139 C Compute and print initial energies.
140 C----------------------------------------------------------------------------
141       nsave=0
142       Kwita=0
143       WhatsUp=0
144       write (iout,'(/80(1h*)/a,i2/80(1h*)/)') 'MCE iteration #',isweep
145       write (iout,'(/80(1h*)/a)') 'Initial energies:'
146       call chainbuild
147       call geom_to_var(nvar,varia)
148       call etotal(energia(0))
149       etot = energia(0)
150       call enerprint(energia(0))
151       if (refstr) then
152         call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,1),nsup,przes,
153      &             obr,non_conv)
154         rms=dsqrt(rms)
155         call contact(.false.,ncont,icont,co)
156         frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
157         write (iout,'(a,f8.3,a,f8.3,a,f8.3)') 
158      &    'RMS deviation from the reference structure:',rms,
159      &    ' % of native contacts:',frac*100,' contact order',co
160         write (istat,'(i10,16(1pe14.5))') 0,
161      &   (energia(print_order(i)),i=1,nprint_ene),
162      &   etot,rms,frac,co
163       else
164         write (istat,'(i10,14(1pe14.5))') 0,
165      &   (energia(print_order(i)),i=1,nprint_ene),etot
166       endif
167 c     close(istat)
168       neneval=neneval+1
169       if (.not. ent_read) then
170 C Initialize the entropy array
171 #ifdef MPL
172 C Collect total energies from other processors.
173         etot_temp=etot
174         etot_all(0)=etot
175         call mp_gather(etot_temp,etot_all,8,MasterID,cgGroupID)
176         if (MyID.eq.MasterID) then
177 C Get the lowest and the highest energy. 
178           print *,'MASTER: etot_temp: ',(etot_all(i),i=0,nprocs-1),
179      &     ' emin=',emin,' emax=',emax
180           emin=1.0D10
181           emax=-1.0D10
182           do i=0,nprocs
183             if (emin.gt.etot_all(i)) emin=etot_all(i)
184             if (emax.lt.etot_all(i)) emax=etot_all(i)
185           enddo
186           emax=emin+e_up
187         endif ! MyID.eq.MasterID
188         etot_all(1)=emin
189         etot_all(2)=emax
190         print *,'Processor',MyID,' calls MP_BCAST to send/recv etot_all'
191         call mp_bcast(etot_all(1),16,MasterID,cgGroupID)
192         print *,'Processor',MyID,' MP_BCAST to send/recv etot_all ended'
193         if (MyID.ne.MasterID) then
194           print *,'Processor:',MyID,etot_all(1),etot_all(2),
195      &          etot_all(1),etot_all(2)
196           emin=etot_all(1)
197           emax=etot_all(2)
198         endif ! MyID.ne.MasterID
199         write (iout,*) 'After MP_GATHER etot_temp=',
200      &                 etot_temp,' emin=',emin
201 #else
202         emin=etot
203         emax=emin+e_up
204         indminn=0
205         indmin=0
206 #endif
207         IF (MULTICAN) THEN
208 C Multicanonical sampling - start from Boltzmann distribution
209           do i=-max_ene,max_ene
210             entropy(i)=(emin+i*delte)*betbol 
211           enddo
212         ELSE
213 C Entropic sampling - start from uniform distribution of the density of states
214           do i=-max_ene,max_ene
215             entropy(i)=0.0D0
216           enddo
217         ENDIF ! MULTICAN
218         write (iout,'(/a)') 'Initial entropy'
219         do i=indminn,indmaxx
220           write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
221         enddo
222         if (isweep.eq.1) then
223           emax=emin+e_up
224           indminn=0
225           indmin=0
226           indmaxx=indminn+nbins
227           indmax=indmaxx
228         endif ! isweep.eq.1
229       endif ! .not. ent_read
230 #ifdef MPL
231       call recv_stop_sig(Kwita)
232       if (whatsup.eq.1) then
233         call send_stop_sig(-2)
234         not_done=.false.
235       else if (whatsup.le.-2) then
236         not_done=.false.
237       else if (whatsup.eq.2) then
238         not_done=.false.
239       else 
240         not_done=.true.
241       endif
242 #else
243       not_done=.true.
244 #endif 
245       write (iout,'(/80(1h*)/20x,a/80(1h*))')
246      &    'Enter Monte Carlo procedure.'
247       close(igeom)
248       call briefout(0,etot)
249       do i=1,nvar
250         varold(i)=varia(i)
251       enddo
252       eold=etot
253       call entropia(eold,sold,indeold)
254 C NACC is the counter for the accepted conformations of a given processor
255       nacc=0
256 C NACC_TOT counts the total number of accepted conformations
257       nacc_tot=0
258 C Main loop.
259 c----------------------------------------------------------------------------
260 C Zero out average energies
261       do i=0,n_ene
262         energia_ave(i)=0.0d0
263       enddo
264 C Initialize energy histogram
265       do i=-max_ene,max_ene
266         nhist(i)=0.0D0
267       enddo   ! i
268 C Zero out iteration counter.
269       it=0
270       do j=1,nvar
271         varold(j)=varia(j)
272        enddo
273 C Begin MC iteration loop.
274       do while (not_done)
275         it=it+1
276 C Initialize local counter.
277         ntrial=0 ! # of generated non-overlapping confs.
278         noverlap=0 ! # of overlapping confs.
279         accepted=.false.
280         do while (.not. accepted .and. WhatsUp.eq.0 .and. Kwita.eq.0)
281           ntrial=ntrial+1
282 C Retrieve the angles of previously accepted conformation
283           do j=1,nvar
284             varia(j)=varold(j)
285           enddo
286           call var_to_geom(nvar,varia)
287 C Rebuild the chain.
288           call chainbuild
289           MoveType=0
290           nbond=0
291           lprint=.true.
292 C Decide whether to take a conformation from the pool or generate/perturb one
293 C randomly
294           from_pool=ran_number(0.0D0,1.0D0)
295           if (npool.gt.0 .and. from_pool.lt.pool_fraction) then
296 C Throw a dice to choose the conformation from the pool
297             ii=iran_num(1,npool)
298             do i=1,nvar
299               varia(i)=xpool(i,ii)
300             enddo
301             call var_to_geom(nvar,varia)
302             call chainbuild  
303 cd          call intout
304 cd          write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
305             if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it)
306      &        write (iout,'(a,i3,a,f10.5)') 
307      &     'Try conformation',ii,' from the pool energy=',epool(ii)
308             MoveType=-1
309             moves(-1)=moves(-1)+1
310           else
311 C Decide whether to generate a random conformation or perturb the old one
312           RandOrPert=ran_number(0.0D0,1.0D0)
313           if (RandOrPert.gt.RanFract) then
314             if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) 
315      &        write (iout,'(a)') 'Perturbation-generated conformation.'
316             call perturb(error,lprint,MoveType,nbond,0.1D0)
317             if (error) goto 20
318             if (MoveType.lt.1 .or. MoveType.gt.MaxMoveType) then
319               write (iout,'(/a,i7,a/)') 'Error - unknown MoveType=',
320      &           MoveType,' returned from PERTURB.'
321               goto 20
322             endif
323             call chainbuild
324           else
325             MoveType=0
326             moves(0)=moves(0)+1
327             nstart_grow=iran_num(3,nres)
328             if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) 
329      &        write (iout,'(2a,i3)') 'Random-generated conformation',
330      &        ' - chain regrown from residue',nstart_grow
331             call gen_rand_conf(nstart_grow,*30)
332           endif
333           call geom_to_var(nvar,varia)
334           endif ! pool
335 Cd        write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
336           ngen=ngen+1
337           if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) 
338      &      write (iout,'(a,i5,a,i10,a,i10)') 
339      &   'Processor',MyId,' trial move',ntrial,' total generated:',ngen
340           if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) 
341      &      write (*,'(a,i5,a,i10,a,i10)') 
342      &   'Processor',MyId,' trial move',ntrial,' total generated:',ngen
343           call etotal(energia(0))
344           etot = energia(0)
345           neneval=neneval+1 
346 cd        call enerprint(energia(0))
347 cd        write(iout,*)'it=',it,' etot=',etot
348           if (etot-elowest.gt.overlap_cut) then
349             if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) 
350      &        write (iout,'(a,i5,a,1pe14.5)')  'Iteration',it,
351      &       ' Overlap detected in the current conf.; energy is',etot
352             accepted=.false.
353             noverlap=noverlap+1
354             if (noverlap.gt.maxoverlap) then
355               write (iout,'(a)') 'Too many overlapping confs.'
356               goto 20
357             endif
358           else
359 C--------------------------------------------------------------------------
360 C... Acceptance test
361 C--------------------------------------------------------------------------
362             accepted=.false.
363             if (WhatsUp.eq.0) 
364      &      call accept_mc(it,etot,eold,scur,sold,varia,varold,accepted)
365             if (accepted) then
366               nacc=nacc+1
367               nacc_tot=nacc_tot+1
368               if (elowest.gt.etot) then
369                 elowest=etot
370                 do i=1,nvar
371                   var_lowest(i)=varia(i)
372                 enddo
373               endif
374               if (ehighest.lt.etot) ehighest=etot
375               moves_acc(MoveType)=moves_acc(MoveType)+1
376               if (MoveType.eq.1) then
377                 nbond_acc(nbond)=nbond_acc(nbond)+1
378               endif
379 C Compare with reference structure.
380               if (refstr) then
381                 call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,1),
382      &                      nsup,przes,obr,non_conv)
383                 rms=dsqrt(rms)
384                 call contact(.false.,ncont,icont,co)
385                 frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
386               endif ! refstr
387 C
388 C Periodically save average energies and confs.
389 C
390               do i=0,n_ene
391                 energia_ave(i)=energia_ave(i)+energia(i)
392               enddo
393               moves(MaxMoveType+1)=nmove
394               moves_acc(MaxMoveType+1)=nacc
395               IF ((it/save_frequency)*save_frequency.eq.it) THEN
396                 do i=0,n_ene
397                   energia_ave(i)=energia_ave(i)/save_frequency
398                 enddo
399                 etot_ave=energia_ave(0)
400 C#ifdef AIX
401 C                open (istat,file=statname,position='append')
402 C#else
403 C                open (istat,file=statname,access='append')
404 Cendif
405                 if (print_mc.gt.0) 
406      &            write (iout,'(80(1h*)/20x,a,i20)')
407      &                             'Iteration #',it
408                 if (refstr .and. print_mc.gt.0)  then
409                   write (iout,'(a,f8.3,a,f8.3,a,f8.3)') 
410      &            'RMS deviation from the reference structure:',rms,
411      &            ' % of native contacts:',frac*100,' contact order:',co
412                 endif
413                 if (print_stat) then
414                   if (refstr) then
415                     write (istat,'(i10,10(1pe14.5))') it,
416      &              (energia_ave(print_order(i)),i=1,nprint_ene),
417      &                etot_ave,rms_ave,frac_ave
418                   else
419                     write (istat,'(i10,10(1pe14.5))') it,
420      &              (energia_ave(print_order(i)),i=1,nprint_ene),
421      &               etot_ave
422                   endif
423                 endif 
424 c               close(istat)
425                 if (print_mc.gt.0) 
426      &            call statprint(nacc,nfun,iretcode,etot,elowest)
427 C Print internal coordinates.
428                 if (print_int) call briefout(nacc,etot)
429                 do i=0,n_ene
430                   energia_ave(i)=0.0d0
431                 enddo
432               ENDIF ! ( (it/save_frequency)*save_frequency.eq.it)
433 C Update histogram
434               inde=icialosc((etot-emin)/delte)
435               nhist(inde)=nhist(inde)+1.0D0
436 #ifdef MPL
437               if ( (it/message_frequency)*message_frequency.eq.it
438      &                              .and. (MyID.ne.MasterID) ) then
439                 call recv_stop_sig(Kwita)
440                 call send_MCM_info(message_frequency)
441               endif
442 #endif
443 C Store the accepted conf. and its energy.
444               eold=etot
445               sold=scur
446               do i=1,nvar
447                 varold(i)=varia(i)
448               enddo
449 #ifdef MPL
450               if (Kwita.eq.0) call recv_stop_sig(kwita)
451 #endif
452             endif ! accepted
453           endif ! overlap
454 #ifdef MPL
455           if (MyID.eq.MasterID .and. 
456      &        (it/message_frequency)*message_frequency.eq.it) then
457             call receive_MC_info
458             if (nacc_tot.ge.maxacc) accepted=.true.
459           endif
460 #endif
461 C         if ((ntrial.gt.maxtrial_iter 
462 C    &       .or. (it/pool_read_freq)*pool_read_freq.eq.it) 
463 C    &       .and. npool.gt.0) then
464 C Take a conformation from the pool
465 C           ii=iran_num(1,npool)
466 C           do i=1,nvar
467 C             varold(i)=xpool(i,ii)
468 C           enddo
469 C           if (ntrial.gt.maxtrial_iter) 
470 C    &      write (iout,*) 'Iteration',it,' max. # of trials exceeded.'
471 C           write (iout,*) 
472 C    &     'Take conformation',ii,' from the pool energy=',epool(ii)
473 C           if (print_mc.gt.2)
474 C    &      write (iout,'(10f8.3)') (rad2deg*varold(i),i=1,nvar)
475 C           ntrial=0
476 C           eold=epool(ii)
477 C           call entropia(eold,sold,indeold)
478 C           accepted=.true.
479 C        endif ! (ntrial.gt.maxtrial_iter .and. npool.gt.0)
480    30    continue
481         enddo ! accepted
482 #ifdef MPL
483         if (MyID.eq.MasterID .and.
484      &      (it/message_frequency)*message_frequency.eq.it) then
485           call receive_MC_info
486         endif
487         if (Kwita.eq.0) call recv_stop_sig(kwita)
488 #endif
489         if (ovrtim()) WhatsUp=-1
490 cd      write (iout,*) 'WhatsUp=',WhatsUp,' Kwita=',Kwita
491         not_done = (nacc_tot.lt.maxacc) .and. (WhatsUp.eq.0) 
492      &         .and. (Kwita.eq.0)
493 cd      write (iout,*) 'not_done=',not_done
494 #ifdef MPL
495         if (Kwita.lt.0) then
496           print *,'Processor',MyID,
497      &    ' has received STOP signal =',Kwita,' in EntSamp.'
498 cd        print *,'not_done=',not_done
499           if (Kwita.lt.-1) WhatsUp=Kwita
500           if (MyID.ne.MasterID) call send_MCM_info(-1)
501         else if (nacc_tot.ge.maxacc) then
502           print *,'Processor',MyID,' calls send_stop_sig,',
503      &     ' because a sufficient # of confs. have been collected.'
504 cd        print *,'not_done=',not_done
505           call send_stop_sig(-1)
506           if (MyID.ne.MasterID) call send_MCM_info(-1)
507         else if (WhatsUp.eq.-1) then
508           print *,'Processor',MyID,
509      &               ' calls send_stop_sig because of timeout.'
510 cd        print *,'not_done=',not_done
511           call send_stop_sig(-2)
512           if (MyID.ne.MasterID) call send_MCM_info(-1)
513         endif
514 #endif
515       enddo ! not_done
516
517 C-----------------------------------------------------------------
518 C... Construct energy histogram & update entropy
519 C-----------------------------------------------------------------
520       go to 21
521    20 WhatsUp=-3
522 #ifdef MPL
523       write (iout,*) 'Processor',MyID,
524      &       ' is broadcasting ERROR-STOP signal.'
525       write (*,*) 'Processor',MyID,
526      &       ' is broadcasting ERROR-STOP signal.'
527       call send_stop_sig(-3)
528       if (MyID.ne.MasterID) call send_MCM_info(-1)
529 #endif
530    21 continue
531       write (iout,'(/a)') 'Energy histogram'
532       do i=-100,100
533         write (iout,'(i5,2f20.5)') i,emin+i*delte,nhist(i)
534       enddo
535 #ifdef MPL
536 C Wait until every processor has sent complete MC info.
537       if (MyID.eq.MasterID) then
538         not_done=.true.
539         do while (not_done)
540 C         write (*,*) 'The IFINISH array:'
541 C         write (*,*) (ifinish(i),i=1,nctasks)
542           not_done=.false.
543           do i=2,nctasks
544             not_done=not_done.or.(ifinish(i).ge.0)
545           enddo
546           if (not_done) call receive_MC_info
547         enddo
548       endif
549 C Make collective histogram from the work of all processors.
550       msglen=(2*max_ene+1)*8
551       print *,
552      & 'Processor',MyID,' calls MP_REDUCE to send/receive histograms',
553      & ' msglen=',msglen
554       call mp_reduce(nhist,nhist1,msglen,MasterID,d_vadd,
555      &               cgGroupID)
556       print *,'Processor',MyID,' MP_REDUCE accomplished for histogr.'
557       do i=-max_ene,max_ene
558         nhist(i)=nhist1(i)
559       enddo
560 C Collect min. and max. energy
561       print *,
562      &'Processor',MyID,' calls MP_REDUCE to send/receive energy borders'
563       call mp_reduce(elowest,elowest1,8,MasterID,d_vmin,cgGroupID)
564       call mp_reduce(ehighest,ehighest1,8,MasterID,d_vmax,cgGroupID)
565       print *,'Processor',MyID,' MP_REDUCE accomplished for energies.'
566       IF (MyID.eq.MasterID) THEN
567         elowest=elowest1
568         ehighest=ehighest1
569 #endif
570         write (iout,'(a,i10)') '# of accepted confs:',nacc_tot
571         write (iout,'(a,f10.5,a,f10.5)') 'Lowest energy:',elowest,
572      & ' Highest energy',ehighest
573         indmin=icialosc((elowest-emin)/delte)
574         imdmax=icialosc((ehighest-emin)/delte)
575         if (indmin.lt.indminn) then 
576           emax=emin+indmin*delte+e_up
577           indmaxx=indmin+nbins
578           indminn=indmin
579         endif
580         if (.not.ent_read) ent_read=.true.
581         write(iout,*)'indminn=',indminn,' indmaxx=',indmaxx
582 C Update entropy (density of states)
583         do i=indmin,indmax
584           if (nhist(i).gt.0) then
585             entropy(i)=entropy(i)+dlog(nhist(i)+0.0D0)
586           endif
587         enddo
588         write (iout,'(/80(1h*)/a,i2/80(1h*)/)') 
589      &        'End of macroiteration',isweep
590         write (iout,'(a,f10.5,a,f10.5)') 'Elowest=',elowest,
591      &      ' Ehighest=',ehighest
592         write (iout,'(/a)') 'Energy histogram'
593         do i=indminn,indmaxx
594           write (iout,'(i5,2f20.5)') i,emin+i*delte,nhist(i)
595         enddo
596         write (iout,'(/a)') 'Entropy'
597         do i=indminn,indmaxx
598           write (iout,'(i5,2f20.5)') i,emin+i*delte,entropy(i)
599         enddo
600 C-----------------------------------------------------------------
601 C... End of energy histogram construction
602 C-----------------------------------------------------------------
603 #ifdef MPL
604       ELSE
605         if (.not. ent_read) ent_read=.true.
606       ENDIF ! MyID .eq. MaterID
607       if (MyID.eq.MasterID) then
608         itemp(1)=indminn
609         itemp(2)=indmaxx
610       endif
611       print *,'before mp_bcast processor',MyID,' indminn=',indminn,
612      & ' indmaxx=',indmaxx,' itemp=',itemp(1),itemp(2)
613       call mp_bcast(itemp(1),8,MasterID,cgGroupID)
614       call mp_bcast(emax,8,MasterID,cgGroupID)
615       print *,'after mp_bcast processor',MyID,' indminn=',indminn,
616      & ' indmaxx=',indmaxx,' itemp=',itemp(1),itemp(2)
617       if (MyID .ne. MasterID) then
618         indminn=itemp(1)
619         indmaxx=itemp(2)
620       endif
621       msglen=(indmaxx-indminn+1)*8
622       print *,'processor',MyID,' calling mp_bcast msglen=',msglen,
623      & ' indminn=',indminn,' indmaxx=',indmaxx,' isweep=',isweep
624       call mp_bcast(entropy(indminn),msglen,MasterID,cgGroupID)
625       IF(MyID.eq.MasterID .and. .not. ovrtim() .and. WhatsUp.ge.0)THEN
626         open (ientout,file=entname,status='unknown')
627         write (ientout,'(2i5,2e25.17)') indminn,indmaxx,emin,emax
628         do i=indminn,indmaxx
629           write (ientout,'(i5,f10.5,f20.15)') i,emin+i*delte,entropy(i)
630         enddo
631         close(ientout)
632       ELSE
633         write (iout,*) 'Received from master:'
634         write (iout,*) 'indminn=',indminn,' indmaxx=',indmaxx,
635      &                 ' emin=',emin,' emax=',emax
636         write (iout,'(/a)') 'Entropy'
637         do i=indminn,indmaxx
638            write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
639         enddo
640       ENDIF ! MyID.eq.MasterID
641       print *,'Processor',MyID,' calls MP_GATHER'
642       call mp_gather(nbond_move,nbond_move1,4*Nbm,MasterID,
643      &               cgGroupID)
644       call mp_gather(nbond_acc,nbond_acc1,4*Nbm,MasterID,
645      &               cgGroupID)
646       print *,'Processor',MyID,' MP_GATHER call accomplished'
647       if (MyID.eq.MasterID) then
648
649         write (iout,'(/80(1h*)/20x,a)') 'Summary run statistics:'
650         call statprint(nacc_tot,nfun,iretcode,etot,elowest)
651         write (iout,'(a)') 
652      &   'Statistics of multiple-bond motions. Total motions:' 
653         write (iout,'(8i10)') (nbond_move(i),i=1,Nbm)
654         write (iout,'(a)') 'Accepted motions:'
655         write (iout,'(8i10)') (nbond_acc(i),i=1,Nbm)
656
657         write (iout,'(a)') 
658      & 'Statistics of multi-bond moves of respective processors:'
659         do iproc=1,Nprocs-1
660           do i=1,Nbm
661             ind=iproc*nbm+i
662             nbond_move(i)=nbond_move(i)+nbond_move1(ind)
663             nbond_acc(i)=nbond_acc(i)+nbond_acc1(ind)
664           enddo
665         enddo
666         do iproc=0,NProcs-1
667           write (iout,*) 'Processor',iproc,' nbond_move:', 
668      &        (nbond_move1(iproc*nbm+i),i=1,Nbm),
669      &        ' nbond_acc:',(nbond_acc1(iproc*nbm+i),i=1,Nbm)
670         enddo
671       endif
672       call mp_gather(moves,moves1,4*(MaxMoveType+3),MasterID,
673      &               cgGroupID)
674       call mp_gather(moves_acc,moves_acc1,4*(MaxMoveType+3),
675      &               MasterID,cgGroupID)
676       if (MyID.eq.MasterID) then
677         do iproc=1,Nprocs-1 
678           do i=-1,MaxMoveType+1
679             moves(i)=moves(i)+moves1(i,iproc)
680             moves_acc(i)=moves_acc(i)+moves_acc1(i,iproc)
681           enddo
682         enddo
683         nmove=0
684         do i=0,MaxMoveType+1
685           nmove=nmove+moves(i)
686         enddo
687         do iproc=0,NProcs-1
688           write (iout,*) 'Processor',iproc,' moves',
689      &     (moves1(i,iproc),i=0,MaxMoveType+1),
690      &    ' moves_acc:',(moves_acc1(i,iproc),i=0,MaxMoveType+1)
691         enddo   
692       endif
693 #else
694       open (ientout,file=entname,status='unknown')
695       write (ientout,'(2i5,2e25.17)') indminn,indmaxx,emin,emax
696       do i=indminn,indmaxx
697         write (ientout,'(i5,f10.5,f20.15)') i,emin+i*delte,entropy(i)
698       enddo
699       close(ientout)
700 #endif
701       write (iout,'(/80(1h*)/20x,a)') 'Summary run statistics:'
702       call statprint(nacc_tot,nfun,iretcode,etot,elowest)
703       write (iout,'(a)') 
704      & 'Statistics of multiple-bond motions. Total motions:' 
705       write (iout,'(8i10)') (nbond_move(i),i=1,Nbm)
706       write (iout,'(a)') 'Accepted motions:'
707       write (iout,'(8i10)') (nbond_acc(i),i=1,Nbm)
708       if (ovrtim() .or. WhatsUp.lt.0) return
709
710 C---------------------------------------------------------------------------
711       ENDDO ! ISWEEP
712 C---------------------------------------------------------------------------
713
714       runtime=tcpu()
715
716       if (isweep.eq.nsweep .and. it.ge.maxacc)
717      &write (iout,'(/80(1h*)/20x,a/80(1h*)/)') 'All iterations done.'
718       return
719       end
720 c------------------------------------------------------------------------------
721       subroutine accept_mc(it,ecur,eold,scur,sold,x,xold,accepted)
722       implicit real*8 (a-h,o-z)
723       include 'DIMENSIONS'
724       include 'COMMON.MCM'
725       include 'COMMON.MCE'
726       include 'COMMON.IOUNITS'
727       include 'COMMON.VAR'
728 #ifdef MPL
729       include 'COMMON.INFO'
730 #endif
731       include 'COMMON.GEO'
732       double precision ecur,eold,xx,ran_number,bol
733       double precision x(maxvar),xold(maxvar)
734       logical accepted
735 C Check if the conformation is similar.
736 cd    write (iout,*) 'Enter ACCEPTING'
737 cd    write (iout,*) 'Old PHI angles:'
738 cd    write (iout,*) (rad2deg*xold(i),i=1,nphi)
739 cd    write (iout,*) 'Current angles'
740 cd    write (iout,*) (rad2deg*x(i),i=1,nphi)
741 cd    ddif=dif_ang(nphi,x,xold)
742 cd    write (iout,*) 'Angle norm:',ddif
743 cd    write (iout,*) 'ecur=',ecur,' emax=',emax
744       if (ecur.gt.emax) then
745         accepted=.false.
746         if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it)
747      & write (iout,'(a)') 'Conformation rejected as too high in energy'
748         return
749       endif
750 C Else evaluate the entropy of the conf and compare it with that of the previous
751 C one.
752       call entropia(ecur,scur,indecur)
753 cd    print *,'Processor',MyID,' ecur=',ecur,' indecur=',indecur,
754 cd   & ' scur=',scur,' eold=',eold,' sold=',sold
755 cd    print *,'deix=',deix,' dent=',dent,' delte=',delte
756       if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) then
757         write(iout,*)'it=',it,'ecur=',ecur,' indecur=',indecur,
758      &   ' scur=',scur
759         write(iout,*)'eold=',eold,' sold=',sold
760       endif
761       if (scur.le.sold) then
762         accepted=.true.
763       else
764 C Else carry out acceptance test
765         xx=ran_number(0.0D0,1.0D0) 
766         xxh=scur-sold
767         if (xxh.gt.50.0D0) then
768           bol=0.0D0
769         else
770           bol=exp(-xxh)
771         endif
772         if (bol.gt.xx) then
773           accepted=.true. 
774           if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) 
775      &       write (iout,'(a)') 'Conformation accepted.'
776         else
777           accepted=.false.
778           if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) 
779      &       write (iout,'(a)') 'Conformation rejected.'
780         endif
781       endif
782       return
783       end 
784 c--------------------------------------------------------------------------
785       integer function icialosc(x)
786       double precision x
787       if (x.lt.0.0D0) then
788         icialosc=dint(x)-1
789       else
790         icialosc=dint(x)
791       endif
792       return
793       end 
794 c--------------------------------------------------------------------------
795       subroutine entropia(ecur,scur,indecur)
796       implicit real*8 (a-h,o-z)
797       include 'DIMENSIONS'
798       include 'COMMON.MCM'
799       include 'COMMON.MCE'
800       include 'COMMON.IOUNITS'
801       double precision ecur,scur
802       integer indecur
803       indecur=icialosc((ecur-emin)/delte)
804       if (iabs(indecur).gt.max_ene) then
805         if ((it/print_freq)*it.eq.it) write (iout,'(a,2i5)') 
806      &   'Accepting: Index out of range:',indecur
807         scur=1000.0D0 
808       else if (indecur.ge.indmaxx) then
809         scur=entropy(indecur)
810         if (print_mc.gt.0 .and. (it/print_freq)*it.eq.it) 
811      &    write (iout,*)'Energy boundary reached',
812      &            indmaxx,indecur,entropy(indecur)
813       else
814         deix=ecur-(emin+indecur*delte)
815         dent=entropy(indecur+1)-entropy(indecur)
816         scur=entropy(indecur)+(dent/delte)*deix
817       endif
818       return
819       end