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