wham and cluster_wham Adam's new constr_dist multichain
[unres.git] / source / unres / src_MD-M / entmcm.F
1       subroutine entmcm
2 C Does modified entropic sampling in the space of minima.
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
21       integer conf_comp
22       double precision RandOrPert
23       double precision varia(maxvar),elowest,ehighest,eold
24       double precision przes(3),obr(3,3)
25       double precision varold(maxvar)
26       logical non_conv
27       double precision energia(0:n_ene),energia_ave(0:n_ene)
28 C
29 cd    write (iout,*) 'print_mc=',print_mc
30       WhatsUp=0
31       maxtrial_iter=50
32 c---------------------------------------------------------------------------
33 C Initialize counters.
34 c---------------------------------------------------------------------------
35 C Total number of generated confs.
36       ngen=0
37 C Total number of moves. In general this won't be equal to the number of
38 C attempted moves, because we may want to reject some "bad" confs just by
39 C overlap check.
40       nmove=0
41 C Total number of shift (nbond_move(1)), spike, crankshaft, three-bond,...
42 C motions.
43       do i=1,nres
44         nbond_move(i)=0
45       enddo
46 C Initialize total and accepted number of moves of various kind.
47       do i=0,MaxMoveType
48         moves(i)=0
49         moves_acc(i)=0
50       enddo
51 C Total number of energy evaluations.
52       neneval=0
53       nfun=0
54       indminn=-max_ene
55       indmaxx=max_ene
56       delte=0.5D0
57       facee=1.0D0/(maxacc*delte)
58       conste=dlog(facee)
59 C Read entropy from previous simulations. 
60       if (ent_read) then
61         read (ientin,*) indminn,indmaxx,emin,emax 
62         print *,'indminn=',indminn,' indmaxx=',indmaxx,' emin=',emin,
63      &          ' emax=',emax
64         do i=-max_ene,max_ene
65           entropy(i)=(emin+i*delte)*betbol
66         enddo
67         read (ientin,*) (ijunk,ejunk,entropy(i),i=indminn,indmaxx)
68         indmin=indminn
69         indmax=indmaxx
70         write (iout,*) 'indminn=',indminn,' indmaxx=',indmaxx,
71      &                 ' emin=',emin,' emax=',emax
72         write (iout,'(/a)') 'Initial entropy'
73         do i=indminn,indmaxx
74           write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
75         enddo
76       endif ! ent_read
77 C Read the pool of conformations
78       call read_pool
79 C----------------------------------------------------------------------------
80 C Entropy-sampling simulations with continually updated entropy
81 C Loop thru simulations
82 C----------------------------------------------------------------------------
83       DO ISWEEP=1,NSWEEP
84 C----------------------------------------------------------------------------
85 C Take a conformation from the pool
86 C----------------------------------------------------------------------------
87       if (npool.gt.0) then
88         ii=iran_num(1,npool)
89         do i=1,nvar
90           varia(i)=xpool(i,ii)
91         enddo
92         write (iout,*) 'Took conformation',ii,' from the pool energy=',
93      &               epool(ii)
94         call var_to_geom(nvar,varia)
95 C Print internal coordinates of the initial conformation
96         call intout
97       else
98         call gen_rand_conf(1,*20)
99       endif
100 C----------------------------------------------------------------------------
101 C Compute and print initial energies.
102 C----------------------------------------------------------------------------
103       nsave=0
104 #ifdef MPL
105       if (MyID.eq.MasterID) then
106         do i=1,nctasks
107           nsave_part(i)=0
108         enddo
109       endif
110 #endif
111       Kwita=0
112       WhatsUp=0
113       write (iout,'(/80(1h*)/a,i2/80(1h*)/)') 'MCE iteration #',isweep
114       write (iout,'(/80(1h*)/a)') 'Initial energies:'
115       call chainbuild
116       call etotal(energia(0))
117       etot = energia(0)
118       call enerprint(energia(0))
119 C Minimize the energy of the first conformation.
120       if (minim) then
121         call geom_to_var(nvar,varia)
122         call minimize(etot,varia,iretcode,nfun)
123         call etotal(energia(0))
124         etot = energia(0)
125         write (iout,'(/80(1h*)/a/80(1h*))') 
126      &    'Results of the first energy minimization:'
127         call enerprint(energia(0))
128       endif
129       if (refstr) then
130        kkk=1
131         call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,kkk),
132      &nsup,przes,
133      &             obr,non_conv)
134         rms=dsqrt(rms)
135         call contact(.false.,ncont,icont,co)
136         frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
137         write (iout,'(a,f8.3,a,f8.3,a,f8.3)') 
138      &    'RMS deviation from the reference structure:',rms,
139      &    ' % of native contacts:',frac*100,' contact order:',co
140         write (istat,'(i5,11(1pe14.5))') 0,
141      &   (energia(print_order(i)),i=1,nprint_ene),etot,rms,frac,co
142       else
143         write (istat,'(i5,9(1pe14.5))') 0,
144      &   (energia(print_order(i)),i=1,nprint_ene),etot
145       endif
146       close(istat)
147       neneval=neneval+nfun+1
148       if (.not. ent_read) then
149 C Initialize the entropy array
150         do i=-max_ene,max_ene
151          emin=etot
152 C Uncomment the line below for actual entropic sampling (start with uniform
153 C energy distribution).
154 c        entropy(i)=0.0D0
155 C Uncomment the line below for multicanonical sampling (start with Boltzmann
156 C distribution).
157          entropy(i)=(emin+i*delte)*betbol 
158         enddo
159         emax=10000000.0D0
160         emin=etot
161         write (iout,'(/a)') 'Initial entropy'
162         do i=indminn,indmaxx
163           write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
164         enddo
165       endif ! ent_read
166 #ifdef MPL
167       call recv_stop_sig(Kwita)
168       if (whatsup.eq.1) then
169         call send_stop_sig(-2)
170         not_done=.false.
171       else if (whatsup.le.-2) then
172         not_done=.false.
173       else if (whatsup.eq.2) then
174         not_done=.false.
175       else 
176         not_done=.true.
177       endif
178 #else
179       not_done = (iretcode.ne.11)
180 #endif 
181       write (iout,'(/80(1h*)/20x,a/80(1h*))')
182      &    'Enter Monte Carlo procedure.'
183       close(igeom)
184       call briefout(0,etot)
185       do i=1,nvar
186         varold(i)=varia(i)
187       enddo
188       eold=etot
189       indeold=(eold-emin)/delte
190       deix=eold-(emin+indeold*delte)
191       dent=entropy(indeold+1)-entropy(indeold)
192 cd    write (iout,*) 'indeold=',indeold,' deix=',deix,' dent=',dent
193 cd    write (*,*) 'Processor',MyID,' indeold=',indeold,' deix=',deix,
194 cd   & ' dent=',dent
195       sold=entropy(indeold)+(dent/delte)*deix
196       elowest=etot
197       write (iout,*) 'eold=',eold,' sold=',sold,' elowest=',etot
198       write (*,*) 'Processor',MyID,' eold=',eold,' sold=',sold,
199      & ' elowest=',etot
200       if (minim) call zapis(varia,etot)
201       nminima(1)=1.0D0
202 C NACC is the counter for the accepted conformations of a given processor
203       nacc=0
204 C NACC_TOT counts the total number of accepted conformations
205       nacc_tot=0
206 #ifdef MPL
207       if (MyID.eq.MasterID) then
208         call receive_MCM_info
209       else
210         call send_MCM_info(2)
211       endif
212 #endif
213       do iene=indminn,indmaxx
214         nhist(iene)=0.0D0
215       enddo
216       do i=2,maxsave
217         nminima(i)=0.0D0
218       enddo
219 C Main loop.
220 c----------------------------------------------------------------------------
221       elowest=1.0D10
222       ehighest=-1.0D10
223       it=0
224       do while (not_done)
225         it=it+1
226         if (print_mc.gt.0) write (iout,'(80(1h*)/20x,a,i7)')
227      &                             'Beginning iteration #',it
228 C Initialize local counter.
229         ntrial=0 ! # of generated non-overlapping confs.
230         noverlap=0 ! # of overlapping confs.
231         accepted=.false.
232         do while (.not. accepted .and. WhatsUp.eq.0 .and. Kwita.eq.0)
233           ntrial=ntrial+1
234 C Retrieve the angles of previously accepted conformation
235           do j=1,nvar
236             varia(j)=varold(j)
237           enddo
238 cd        write (iout,'(a)') 'Old variables:'
239 cd        write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
240           call var_to_geom(nvar,varia)
241 C Rebuild the chain.
242           call chainbuild
243           MoveType=0
244           nbond=0
245           lprint=.true.
246 C Decide whether to generate a random conformation or perturb the old one
247           RandOrPert=ran_number(0.0D0,1.0D0)
248           if (RandOrPert.gt.RanFract) then
249             if (print_mc.gt.0) 
250      &        write (iout,'(a)') 'Perturbation-generated conformation.'
251             call perturb(error,lprint,MoveType,nbond,1.0D0)
252             if (error) goto 20
253             if (MoveType.lt.1 .or. MoveType.gt.MaxMoveType) then
254               write (iout,'(/a,i7,a/)') 'Error - unknown MoveType=',
255      &           MoveType,' returned from PERTURB.'
256               goto 20
257             endif
258             call chainbuild
259           else
260             MoveType=0
261             moves(0)=moves(0)+1
262             nstart_grow=iran_num(3,nres)
263             if (print_mc.gt.0) 
264      &        write (iout,'(2a,i3)') 'Random-generated conformation',
265      &        ' - chain regrown from residue',nstart_grow
266             call gen_rand_conf(nstart_grow,*30)
267           endif
268           call geom_to_var(nvar,varia)
269 cd        write (iout,'(a)') 'New variables:'
270 cd        write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
271           ngen=ngen+1
272           if (print_mc.gt.0) write (iout,'(a,i5,a,i10,a,i10)') 
273      &   'Processor',MyId,' trial move',ntrial,' total generated:',ngen
274           if (print_mc.gt.0) write (*,'(a,i5,a,i10,a,i10)') 
275      &   'Processor',MyId,' trial move',ntrial,' total generated:',ngen
276           call etotal(energia(0))
277           etot = energia(0)
278 c         call enerprint(energia(0))
279 c         write (iout,'(2(a,1pe14.5))') 'Etot=',Etot,' Elowest=',Elowest
280           if (etot-elowest.gt.overlap_cut) then
281             write (iout,'(a,i5,a,1pe14.5)')  'Iteration',it,
282      &      ' Overlap detected in the current conf.; energy is',etot
283             neneval=neneval+1 
284             accepted=.false.
285             noverlap=noverlap+1
286             if (noverlap.gt.maxoverlap) then
287               write (iout,'(a)') 'Too many overlapping confs.'
288               goto 20
289             endif
290           else
291             if (minim) then
292               call minimize(etot,varia,iretcode,nfun)
293 cd            write (iout,'(a)') 'Variables after minimization:'
294 cd            write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
295               call etotal(energia(0))
296               etot = energia(0)
297               neneval=neneval+nfun+1
298             endif
299             if (print_mc.gt.2) then
300               write (iout,'(a)') 'Total energies of trial conf:'
301               call enerprint(energia(0))
302             else if (print_mc.eq.1) then
303                write (iout,'(a,i6,a,1pe16.6)') 
304      &         'Trial conformation:',ngen,' energy:',etot
305             endif 
306 C--------------------------------------------------------------------------
307 C... Acceptance test
308 C--------------------------------------------------------------------------
309             accepted=.false.
310             if (WhatsUp.eq.0) 
311      &        call accepting(etot,eold,scur,sold,varia,varold,
312      &                       accepted)
313             if (accepted) then
314               nacc=nacc+1
315               nacc_tot=nacc_tot+1
316               if (elowest.gt.etot) elowest=etot
317               if (ehighest.lt.etot) ehighest=etot
318               moves_acc(MoveType)=moves_acc(MoveType)+1
319               if (MoveType.eq.1) then
320                 nbond_acc(nbond)=nbond_acc(nbond)+1
321               endif
322 C Check against conformation repetitions.
323               irep=conf_comp(varia,etot)
324 #if defined(AIX) || defined(PGI)
325               open (istat,file=statname,position='append')
326 #else
327               open (istat,file=statname,access='append')
328 #endif
329               if (refstr) then
330                 kkk=1
331                 call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,kkk),
332      &                     nsup,
333      &                      przes,obr,non_conv)
334                 rms=dsqrt(rms)
335                 call contact(.false.,ncont,icont,co)
336                 frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
337                 if (print_mc.gt.0) 
338      &          write (iout,'(a,f8.3,a,f8.3,a,f8.3)') 
339      &          'RMS deviation from the reference structure:',rms,
340      &          ' % of native contacts:',frac*100,' contact order:',co
341                 if (print_stat)
342      &          write (istat,'(i5,11(1pe14.5))') it,
343      &           (energia(print_order(i)),i=1,nprint_ene),etot,
344      &           rms,frac,co
345               elseif (print_stat) then
346                 write (istat,'(i5,10(1pe14.5))') it,
347      &             (energia(print_order(i)),i=1,nprint_ene),etot
348               endif  
349               close(istat)
350               if (print_mc.gt.1) 
351      &          call statprint(nacc,nfun,iretcode,etot,elowest)
352 C Print internal coordinates.
353               if (print_int) call briefout(nacc,etot)
354 #ifdef MPL
355               if (MyID.ne.MasterID) then
356                 call recv_stop_sig(Kwita)
357 cd              print *,'Processor:',MyID,' STOP=',Kwita
358                 if (irep.eq.0) then
359                   call send_MCM_info(2)
360                 else
361                   call send_MCM_info(1)
362                 endif
363               endif
364 #endif
365 C Store the accepted conf. and its energy.
366               eold=etot
367               sold=scur
368               do i=1,nvar
369                 varold(i)=varia(i)
370               enddo
371               if (irep.eq.0) then
372                 irep=nsave+1
373 cd              write (iout,*) 'Accepted conformation:'
374 cd              write (iout,*) (rad2deg*varia(i),i=1,nphi)
375                 if (minim) call zapis(varia,etot)
376                 do i=1,n_ene
377                   ener(i,nsave)=energia(i)
378                 enddo
379                 ener(n_ene+1,nsave)=etot
380                 ener(n_ene+2,nsave)=frac
381               endif
382               nminima(irep)=nminima(irep)+1.0D0
383 c             print *,'irep=',irep,' nminima=',nminima(irep)
384 #ifdef MPL
385               if (Kwita.eq.0) call recv_stop_sig(kwita)
386 #endif
387             endif ! accepted
388           endif ! overlap
389 #ifdef MPL
390           if (MyID.eq.MasterID) then
391             call receive_MCM_info
392             if (nacc_tot.ge.maxacc) accepted=.true.
393           endif
394 #endif
395           if (ntrial.gt.maxtrial_iter .and. npool.gt.0) then
396 C Take a conformation from the pool
397             ii=iran_num(1,npool)
398             do i=1,nvar
399               varia(i)=xpool(i,ii)
400             enddo
401             write (iout,*) 'Iteration',it,' max. # of trials exceeded.'
402             write (iout,*) 
403      &     'Take conformation',ii,' from the pool energy=',epool(ii)
404             if (print_mc.gt.2)
405      &      write (iout,'(10f8.3)') (rad2deg*varia(i),i=1,nvar)
406             ntrial=0
407          endif ! (ntrial.gt.maxtrial_iter .and. npool.gt.0)
408    30    continue
409         enddo ! accepted
410 #ifdef MPL
411         if (MyID.eq.MasterID) then
412           call receive_MCM_info
413         endif
414         if (Kwita.eq.0) call recv_stop_sig(kwita)
415 #endif
416         if (ovrtim()) WhatsUp=-1
417 cd      write (iout,*) 'WhatsUp=',WhatsUp,' Kwita=',Kwita
418         not_done = (nacc_tot.lt.maxacc) .and. (WhatsUp.eq.0) 
419      &         .and. (Kwita.eq.0)
420 cd      write (iout,*) 'not_done=',not_done
421 #ifdef MPL
422         if (Kwita.lt.0) then
423           print *,'Processor',MyID,
424      &    ' has received STOP signal =',Kwita,' in EntSamp.'
425 cd        print *,'not_done=',not_done
426           if (Kwita.lt.-1) WhatsUp=Kwita
427         else if (nacc_tot.ge.maxacc) then
428           print *,'Processor',MyID,' calls send_stop_sig,',
429      &     ' because a sufficient # of confs. have been collected.'
430 cd        print *,'not_done=',not_done
431           call send_stop_sig(-1)
432         else if (WhatsUp.eq.-1) then
433           print *,'Processor',MyID,
434      &               ' calls send_stop_sig because of timeout.'
435 cd        print *,'not_done=',not_done
436           call send_stop_sig(-2)
437         endif
438 #endif
439       enddo ! not_done
440
441 C-----------------------------------------------------------------
442 C... Construct energy histogram & update entropy
443 C-----------------------------------------------------------------
444       go to 21
445    20 WhatsUp=-3
446 #ifdef MPL
447       write (iout,*) 'Processor',MyID,
448      &       ' is broadcasting ERROR-STOP signal.'
449       write (*,*) 'Processor',MyID,
450      &       ' is broadcasting ERROR-STOP signal.'
451       call send_stop_sig(-3)
452 #endif
453    21 continue
454 #ifdef MPL
455       if (MyID.eq.MasterID) then
456 c       call receive_MCM_results
457         call receive_energies
458 #endif
459       do i=1,nsave
460         if (esave(i).lt.elowest) elowest=esave(i)
461         if (esave(i).gt.ehighest) ehighest=esave(i)
462       enddo
463       write (iout,'(a,i10)') '# of accepted confs:',nacc_tot
464       write (iout,'(a,f10.5,a,f10.5)') 'Lowest energy:',elowest,
465      & ' Highest energy',ehighest
466       if (isweep.eq.1 .and. .not.ent_read) then
467         emin=elowest
468         emax=ehighest
469         write (iout,*) 'EMAX=',emax
470         indminn=0
471         indmaxx=(ehighest-emin)/delte
472         indmin=indminn
473         indmax=indmaxx
474         do i=-max_ene,max_ene
475           entropy(i)=(emin+i*delte)*betbol
476         enddo
477         ent_read=.true.
478       else
479         indmin=(elowest-emin)/delte
480         indmax=(ehighest-emin)/delte
481         if (indmin.lt.indminn) indminn=indmin
482         if (indmax.gt.indmaxx) indmaxx=indmax
483       endif
484       write(iout,*)'indminn=',indminn,' indmaxx=',indmaxx
485 C Construct energy histogram
486       do i=1,nsave
487         inde=(esave(i)-emin)/delte
488         nhist(inde)=nhist(inde)+nminima(i)
489       enddo
490 C Update entropy (density of states)
491       do i=indmin,indmax
492         if (nhist(i).gt.0) then
493           entropy(i)=entropy(i)+dlog(nhist(i)+0.0D0)
494         endif
495       enddo
496 Cd    do i=indmaxx+1
497 Cd      entropy(i)=1.0D+10
498 Cd    enddo
499       write (iout,'(/80(1h*)/a,i2/80(1h*)/)') 
500      &      'End of macroiteration',isweep
501       write (iout,'(a,f10.5,a,f10.5)') 'Elowest=',elowest,
502      &      ' Ehighest=',ehighest
503       write (iout,'(a)') 'Frequecies of minima'
504       do i=1,nsave
505         write (iout,'(i5,f5.0,f10.5)') i,nminima(i),esave(i)
506       enddo
507       write (iout,'(/a)') 'Energy histogram'
508       do i=indminn,indmaxx
509         write (iout,'(i5,2f10.5)') i,emin+i*delte,nhist(i)
510       enddo
511       write (iout,'(/a)') 'Entropy'
512       do i=indminn,indmaxx
513         write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
514       enddo
515 C-----------------------------------------------------------------
516 C... End of energy histogram construction
517 C-----------------------------------------------------------------
518 #ifdef MPL
519         entropy(-max_ene-4)=dfloat(indminn)
520         entropy(-max_ene-3)=dfloat(indmaxx)
521         entropy(-max_ene-2)=emin
522         entropy(-max_ene-1)=emax
523         call send_MCM_update
524 cd      print *,entname,ientout
525         open (ientout,file=entname,status='unknown')
526         write (ientout,'(2i5,2e25.17)') indminn,indmaxx,emin,emax
527         do i=indminn,indmaxx
528           write (ientout,'(i5,f10.5,f20.15)') i,emin+i*delte,entropy(i)
529         enddo
530         close(ientout)
531       else
532         write (iout,'(a)') 'Frequecies of minima'
533         do i=1,nsave
534           write (iout,'(i5,f5.0,f10.5)') i,nminima(i),esave(i)
535         enddo
536 c       call send_MCM_results
537         call send_energies
538         call receive_MCM_update
539         indminn=entropy(-max_ene-4)
540         indmaxx=entropy(-max_ene-3)
541         emin=entropy(-max_ene-2)
542         emax=entropy(-max_ene-1)
543         write (iout,*) 'Received from master:'
544         write (iout,*) 'indminn=',indminn,' indmaxx=',indmaxx,
545      &                 ' emin=',emin,' emax=',emax
546         write (iout,'(/a)') 'Entropy'
547         do i=indminn,indmaxx
548           write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
549         enddo
550       endif
551       if (WhatsUp.lt.-1) return
552 #else
553       if (ovrtim() .or. WhatsUp.lt.0) return
554 #endif
555
556       write (iout,'(/80(1h*)/20x,a)') 'Summary run statistics:'
557       call statprint(nacc,nfun,iretcode,etot,elowest)
558       write (iout,'(a)') 
559      & 'Statistics of multiple-bond motions. Total motions:' 
560       write (iout,'(16i5)') (nbond_move(i),i=1,Nbm)
561       write (iout,'(a)') 'Accepted motions:'
562       write (iout,'(16i5)') (nbond_acc(i),i=1,Nbm)
563       write (iout,'(a,i10)') 'Number of chain regrowths:',nregrow
564       write (iout,'(a,i10)') 'Accepted chain regrowths:',nregrow_acc
565
566 C---------------------------------------------------------------------------
567       ENDDO ! ISWEEP
568 C---------------------------------------------------------------------------
569
570       runtime=tcpu()
571
572       if (isweep.eq.nsweep .and. it.ge.maxacc)
573      &write (iout,'(/80(1h*)/20x,a/80(1h*)/)') 'All iterations done.'
574       return
575       end
576 c------------------------------------------------------------------------------
577       subroutine accepting(ecur,eold,scur,sold,x,xold,accepted)
578       implicit real*8 (a-h,o-z)
579       include 'DIMENSIONS'
580       include 'COMMON.MCM'
581       include 'COMMON.MCE'
582       include 'COMMON.IOUNITS'
583       include 'COMMON.VAR'
584 #ifdef MPL
585       include 'COMMON.INFO'
586 #endif
587       include 'COMMON.GEO'
588       double precision ecur,eold,xx,ran_number,bol
589       double precision x(maxvar),xold(maxvar)
590       double precision tole /1.0D-1/, tola /5.0D0/
591       logical accepted
592 C Check if the conformation is similar.
593 cd    write (iout,*) 'Enter ACCEPTING'
594 cd    write (iout,*) 'Old PHI angles:'
595 cd    write (iout,*) (rad2deg*xold(i),i=1,nphi)
596 cd    write (iout,*) 'Current angles'
597 cd    write (iout,*) (rad2deg*x(i),i=1,nphi)
598 cd    ddif=dif_ang(nphi,x,xold)
599 cd    write (iout,*) 'Angle norm:',ddif
600 cd    write (iout,*) 'ecur=',ecur,' emax=',emax
601       if (ecur.gt.emax) then
602         accepted=.false.
603         if (print_mc.gt.0)
604      & write (iout,'(a)') 'Conformation rejected as too high in energy'
605         return
606       else if (dabs(ecur-eold).lt.tole .and. 
607      &       dif_ang(nphi,x,xold).lt.tola) then
608         accepted=.false.
609         if (print_mc.gt.0)
610      & write (iout,'(a)') 'Conformation rejected as too similar'
611         return
612       endif
613 C Else evaluate the entropy of the conf and compare it with that of the previous
614 C one.
615       indecur=(ecur-emin)/delte
616       if (iabs(indecur).gt.max_ene) then
617         write (iout,'(a,2i5)') 
618      &   'Accepting: Index out of range:',indecur
619         scur=1000.0D0 
620       else if (indecur.eq.indmaxx) then
621         scur=entropy(indecur)
622         if (print_mc.gt.0) write (iout,*)'Energy boundary reached',
623      &            indmaxx,indecur,entropy(indecur)
624       else
625         deix=ecur-(emin+indecur*delte)
626         dent=entropy(indecur+1)-entropy(indecur)
627         scur=entropy(indecur)+(dent/delte)*deix
628       endif
629 cd    print *,'Processor',MyID,' ecur=',ecur,' indecur=',indecur,
630 cd   & ' scur=',scur,' eold=',eold,' sold=',sold
631 cd    print *,'deix=',deix,' dent=',dent,' delte=',delte
632       if (print_mc.gt.1) then
633         write(iout,*)'ecur=',ecur,' indecur=',indecur,' scur=',scur
634         write(iout,*)'eold=',eold,' sold=',sold
635       endif
636       if (scur.le.sold) then
637         accepted=.true.
638       else
639 C Else carry out acceptance test
640         xx=ran_number(0.0D0,1.0D0) 
641         xxh=scur-sold
642         if (xxh.gt.50.0D0) then
643           bol=0.0D0
644         else
645           bol=exp(-xxh)
646         endif
647         if (bol.gt.xx) then
648           accepted=.true. 
649           if (print_mc.gt.0) write (iout,'(a)') 
650      &    'Conformation accepted.'
651         else
652           accepted=.false.
653           if (print_mc.gt.0) write (iout,'(a)') 
654      & 'Conformation rejected.'
655         endif
656       endif
657       return
658       end 
659 c-----------------------------------------------------------------------------
660       subroutine read_pool
661       implicit real*8 (a-h,o-z)
662       include 'DIMENSIONS'
663       include 'COMMON.IOUNITS'
664       include 'COMMON.GEO'
665       include 'COMMON.MCM'
666       include 'COMMON.MCE'
667       include 'COMMON.VAR'
668       double precision varia(maxvar)
669       print '(a)','Call READ_POOL'
670       do npool=1,max_pool
671         print *,'i=',i
672         read (intin,'(i5,f10.5)',end=10,err=10) iconf,epool(npool)
673         if (epool(npool).eq.0.0D0) goto 10
674         call read_angles(intin,*10)
675         call geom_to_var(nvar,xpool(1,npool))
676       enddo
677       goto 11
678    10 npool=npool-1
679    11 write (iout,'(a,i5)') 'Number of pool conformations:',npool
680       if (print_mc.gt.2) then
681       do i=1,npool
682         write (iout,'(a,i5,a,1pe14.5)') 'Pool conformation',i,' energy',
683      &    epool(i)
684         write (iout,'(10f8.3)') (rad2deg*xpool(j,i),j=1,nvar)
685       enddo
686       endif ! (print_mc.gt.2)
687       return
688       end