X-Git-Url: http://mmka.chem.univ.gda.pl/gitweb/?a=blobdiff_plain;f=source%2Funres%2Fsrc_MD-NEWSC%2Fentmcm.F;fp=source%2Funres%2Fsrc_MD-NEWSC%2Fentmcm.F;h=3c2dc5a2ab4398b37baaff73226dc0f4e9a28af6;hb=7308760ff07636ef6b1ee28d8c3a67a23c14b34b;hp=0000000000000000000000000000000000000000;hpb=9a54ab407f6d0d9d564d52763b3e2136450b9ffc;p=unres.git

diff --git a/source/unres/src_MD-NEWSC/entmcm.F b/source/unres/src_MD-NEWSC/entmcm.F
new file mode 100644
index 0000000..3c2dc5a
--- /dev/null
+++ b/source/unres/src_MD-NEWSC/entmcm.F
@@ -0,0 +1,684 @@
+      subroutine entmcm
+C Does modified entropic sampling in the space of minima.
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.IOUNITS'
+#ifdef MPL
+      include 'COMMON.INFO'
+#endif
+      include 'COMMON.GEO'
+      include 'COMMON.CHAIN'
+      include 'COMMON.MCM'
+      include 'COMMON.MCE'
+      include 'COMMON.CONTACTS'
+      include 'COMMON.CONTROL'
+      include 'COMMON.VAR'
+      include 'COMMON.INTERACT'
+      include 'COMMON.THREAD'
+      include 'COMMON.NAMES'
+      logical accepted,not_done,over,ovrtim,error,lprint
+      integer MoveType,nbond
+      integer conf_comp
+      double precision RandOrPert
+      double precision varia(maxvar),elowest,ehighest,eold
+      double precision przes(3),obr(3,3)
+      double precision varold(maxvar)
+      logical non_conv
+      double precision energia(0:n_ene),energia_ave(0:n_ene)
+C
+cd    write (iout,*) 'print_mc=',print_mc
+      WhatsUp=0
+      maxtrial_iter=50
+c---------------------------------------------------------------------------
+C Initialize counters.
+c---------------------------------------------------------------------------
+C Total number of generated confs.
+      ngen=0
+C Total number of moves. In general this won't be equal to the number of
+C attempted moves, because we may want to reject some "bad" confs just by
+C overlap check.
+      nmove=0
+C Total number of shift (nbond_move(1)), spike, crankshaft, three-bond,...
+C motions.
+      do i=1,nres
+        nbond_move(i)=0
+      enddo
+C Initialize total and accepted number of moves of various kind.
+      do i=0,MaxMoveType
+        moves(i)=0
+        moves_acc(i)=0
+      enddo
+C Total number of energy evaluations.
+      neneval=0
+      nfun=0
+      indminn=-max_ene
+      indmaxx=max_ene
+      delte=0.5D0
+      facee=1.0D0/(maxacc*delte)
+      conste=dlog(facee)
+C Read entropy from previous simulations. 
+      if (ent_read) then
+        read (ientin,*) indminn,indmaxx,emin,emax 
+        print *,'indminn=',indminn,' indmaxx=',indmaxx,' emin=',emin,
+     &          ' emax=',emax
+        do i=-max_ene,max_ene
+          entropy(i)=(emin+i*delte)*betbol
+        enddo
+        read (ientin,*) (ijunk,ejunk,entropy(i),i=indminn,indmaxx)
+        indmin=indminn
+        indmax=indmaxx
+        write (iout,*) 'indminn=',indminn,' indmaxx=',indmaxx,
+     &                 ' emin=',emin,' emax=',emax
+        write (iout,'(/a)') 'Initial entropy'
+        do i=indminn,indmaxx
+          write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
+        enddo
+      endif ! ent_read
+C Read the pool of conformations
+      call read_pool
+C----------------------------------------------------------------------------
+C Entropy-sampling simulations with continually updated entropy
+C Loop thru simulations
+C----------------------------------------------------------------------------
+      DO ISWEEP=1,NSWEEP
+C----------------------------------------------------------------------------
+C Take a conformation from the pool
+C----------------------------------------------------------------------------
+      if (npool.gt.0) then
+        ii=iran_num(1,npool)
+        do i=1,nvar
+          varia(i)=xpool(i,ii)
+        enddo
+        write (iout,*) 'Took conformation',ii,' from the pool energy=',
+     &               epool(ii)
+        call var_to_geom(nvar,varia)
+C Print internal coordinates of the initial conformation
+        call intout
+      else
+        call gen_rand_conf(1,*20)
+      endif
+C----------------------------------------------------------------------------
+C Compute and print initial energies.
+C----------------------------------------------------------------------------
+      nsave=0
+#ifdef MPL
+      if (MyID.eq.MasterID) then
+        do i=1,nctasks
+          nsave_part(i)=0
+        enddo
+      endif
+#endif
+      Kwita=0
+      WhatsUp=0
+      write (iout,'(/80(1h*)/a,i2/80(1h*)/)') 'MCE iteration #',isweep
+      write (iout,'(/80(1h*)/a)') 'Initial energies:'
+      call chainbuild
+      call etotal(energia(0))
+      etot = energia(0)
+      call enerprint(energia(0))
+C Minimize the energy of the first conformation.
+      if (minim) then
+        call geom_to_var(nvar,varia)
+        call minimize(etot,varia,iretcode,nfun)
+        call etotal(energia(0))
+        etot = energia(0)
+        write (iout,'(/80(1h*)/a/80(1h*))') 
+     &    'Results of the first energy minimization:'
+        call enerprint(energia(0))
+      endif
+      if (refstr) then
+        call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup),nsup,przes,
+     &             obr,non_conv)
+        rms=dsqrt(rms)
+        call contact(.false.,ncont,icont,co)
+        frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
+        write (iout,'(a,f8.3,a,f8.3,a,f8.3)') 
+     &    'RMS deviation from the reference structure:',rms,
+     &    ' % of native contacts:',frac*100,' contact order:',co
+        write (istat,'(i5,11(1pe14.5))') 0,
+     &   (energia(print_order(i)),i=1,nprint_ene),etot,rms,frac,co
+      else
+        write (istat,'(i5,9(1pe14.5))') 0,
+     &   (energia(print_order(i)),i=1,nprint_ene),etot
+      endif
+      close(istat)
+      neneval=neneval+nfun+1
+      if (.not. ent_read) then
+C Initialize the entropy array
+        do i=-max_ene,max_ene
+         emin=etot
+C Uncomment the line below for actual entropic sampling (start with uniform
+C energy distribution).
+c        entropy(i)=0.0D0
+C Uncomment the line below for multicanonical sampling (start with Boltzmann
+C distribution).
+         entropy(i)=(emin+i*delte)*betbol 
+        enddo
+        emax=10000000.0D0
+        emin=etot
+        write (iout,'(/a)') 'Initial entropy'
+        do i=indminn,indmaxx
+          write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
+        enddo
+      endif ! ent_read
+#ifdef MPL
+      call recv_stop_sig(Kwita)
+      if (whatsup.eq.1) then
+        call send_stop_sig(-2)
+        not_done=.false.
+      else if (whatsup.le.-2) then
+        not_done=.false.
+      else if (whatsup.eq.2) then
+        not_done=.false.
+      else 
+        not_done=.true.
+      endif
+#else
+      not_done = (iretcode.ne.11)
+#endif 
+      write (iout,'(/80(1h*)/20x,a/80(1h*))')
+     &    'Enter Monte Carlo procedure.'
+      close(igeom)
+      call briefout(0,etot)
+      do i=1,nvar
+        varold(i)=varia(i)
+      enddo
+      eold=etot
+      indeold=(eold-emin)/delte
+      deix=eold-(emin+indeold*delte)
+      dent=entropy(indeold+1)-entropy(indeold)
+cd    write (iout,*) 'indeold=',indeold,' deix=',deix,' dent=',dent
+cd    write (*,*) 'Processor',MyID,' indeold=',indeold,' deix=',deix,
+cd   & ' dent=',dent
+      sold=entropy(indeold)+(dent/delte)*deix
+      elowest=etot
+      write (iout,*) 'eold=',eold,' sold=',sold,' elowest=',etot
+      write (*,*) 'Processor',MyID,' eold=',eold,' sold=',sold,
+     & ' elowest=',etot
+      if (minim) call zapis(varia,etot)
+      nminima(1)=1.0D0
+C NACC is the counter for the accepted conformations of a given processor
+      nacc=0
+C NACC_TOT counts the total number of accepted conformations
+      nacc_tot=0
+#ifdef MPL
+      if (MyID.eq.MasterID) then
+        call receive_MCM_info
+      else
+        call send_MCM_info(2)
+      endif
+#endif
+      do iene=indminn,indmaxx
+        nhist(iene)=0.0D0
+      enddo
+      do i=2,maxsave
+        nminima(i)=0.0D0
+      enddo
+C Main loop.
+c----------------------------------------------------------------------------
+      elowest=1.0D10
+      ehighest=-1.0D10
+      it=0
+      do while (not_done)
+        it=it+1
+        if (print_mc.gt.0) write (iout,'(80(1h*)/20x,a,i7)')
+     &                             'Beginning iteration #',it
+C Initialize local counter.
+        ntrial=0 ! # of generated non-overlapping confs.
+        noverlap=0 ! # of overlapping confs.
+        accepted=.false.
+        do while (.not. accepted .and. WhatsUp.eq.0 .and. Kwita.eq.0)
+          ntrial=ntrial+1
+C Retrieve the angles of previously accepted conformation
+          do j=1,nvar
+            varia(j)=varold(j)
+          enddo
+cd        write (iout,'(a)') 'Old variables:'
+cd        write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
+          call var_to_geom(nvar,varia)
+C Rebuild the chain.
+          call chainbuild
+          MoveType=0
+          nbond=0
+          lprint=.true.
+C Decide whether to generate a random conformation or perturb the old one
+          RandOrPert=ran_number(0.0D0,1.0D0)
+          if (RandOrPert.gt.RanFract) then
+            if (print_mc.gt.0) 
+     &        write (iout,'(a)') 'Perturbation-generated conformation.'
+            call perturb(error,lprint,MoveType,nbond,1.0D0)
+            if (error) goto 20
+            if (MoveType.lt.1 .or. MoveType.gt.MaxMoveType) then
+              write (iout,'(/a,i7,a/)') 'Error - unknown MoveType=',
+     &           MoveType,' returned from PERTURB.'
+              goto 20
+            endif
+            call chainbuild
+          else
+            MoveType=0
+            moves(0)=moves(0)+1
+            nstart_grow=iran_num(3,nres)
+            if (print_mc.gt.0) 
+     &        write (iout,'(2a,i3)') 'Random-generated conformation',
+     &        ' - chain regrown from residue',nstart_grow
+            call gen_rand_conf(nstart_grow,*30)
+          endif
+          call geom_to_var(nvar,varia)
+cd        write (iout,'(a)') 'New variables:'
+cd        write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
+          ngen=ngen+1
+          if (print_mc.gt.0) write (iout,'(a,i5,a,i10,a,i10)') 
+     &   'Processor',MyId,' trial move',ntrial,' total generated:',ngen
+          if (print_mc.gt.0) write (*,'(a,i5,a,i10,a,i10)') 
+     &   'Processor',MyId,' trial move',ntrial,' total generated:',ngen
+          call etotal(energia(0))
+          etot = energia(0)
+c         call enerprint(energia(0))
+c         write (iout,'(2(a,1pe14.5))') 'Etot=',Etot,' Elowest=',Elowest
+          if (etot-elowest.gt.overlap_cut) then
+            write (iout,'(a,i5,a,1pe14.5)')  'Iteration',it,
+     &      ' Overlap detected in the current conf.; energy is',etot
+            neneval=neneval+1 
+            accepted=.false.
+            noverlap=noverlap+1
+            if (noverlap.gt.maxoverlap) then
+              write (iout,'(a)') 'Too many overlapping confs.'
+              goto 20
+            endif
+          else
+            if (minim) then
+              call minimize(etot,varia,iretcode,nfun)
+cd            write (iout,'(a)') 'Variables after minimization:'
+cd            write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
+              call etotal(energia(0))
+              etot = energia(0)
+              neneval=neneval+nfun+1
+            endif
+            if (print_mc.gt.2) then
+              write (iout,'(a)') 'Total energies of trial conf:'
+              call enerprint(energia(0))
+            else if (print_mc.eq.1) then
+               write (iout,'(a,i6,a,1pe16.6)') 
+     &         'Trial conformation:',ngen,' energy:',etot
+            endif 
+C--------------------------------------------------------------------------
+C... Acceptance test
+C--------------------------------------------------------------------------
+            accepted=.false.
+            if (WhatsUp.eq.0) 
+     &        call accepting(etot,eold,scur,sold,varia,varold,
+     &                       accepted)
+            if (accepted) then
+              nacc=nacc+1
+              nacc_tot=nacc_tot+1
+              if (elowest.gt.etot) elowest=etot
+              if (ehighest.lt.etot) ehighest=etot
+              moves_acc(MoveType)=moves_acc(MoveType)+1
+              if (MoveType.eq.1) then
+                nbond_acc(nbond)=nbond_acc(nbond)+1
+              endif
+C Check against conformation repetitions.
+              irep=conf_comp(varia,etot)
+#if defined(AIX) || defined(PGI)
+              open (istat,file=statname,position='append')
+#else
+              open (istat,file=statname,access='append')
+#endif
+              if (refstr) then
+                call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup),nsup,
+     &                      przes,obr,non_conv)
+                rms=dsqrt(rms)
+                call contact(.false.,ncont,icont,co)
+                frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
+                if (print_mc.gt.0) 
+     &          write (iout,'(a,f8.3,a,f8.3,a,f8.3)') 
+     &          'RMS deviation from the reference structure:',rms,
+     &          ' % of native contacts:',frac*100,' contact order:',co
+                if (print_stat)
+     &          write (istat,'(i5,11(1pe14.5))') it,
+     &           (energia(print_order(i)),i=1,nprint_ene),etot,
+     &           rms,frac,co
+              elseif (print_stat) then
+                write (istat,'(i5,10(1pe14.5))') it,
+     &             (energia(print_order(i)),i=1,nprint_ene),etot
+              endif  
+              close(istat)
+              if (print_mc.gt.1) 
+     &          call statprint(nacc,nfun,iretcode,etot,elowest)
+C Print internal coordinates.
+              if (print_int) call briefout(nacc,etot)
+#ifdef MPL
+              if (MyID.ne.MasterID) then
+                call recv_stop_sig(Kwita)
+cd              print *,'Processor:',MyID,' STOP=',Kwita
+                if (irep.eq.0) then
+                  call send_MCM_info(2)
+                else
+                  call send_MCM_info(1)
+                endif
+              endif
+#endif
+C Store the accepted conf. and its energy.
+              eold=etot
+              sold=scur
+              do i=1,nvar
+                varold(i)=varia(i)
+              enddo
+              if (irep.eq.0) then
+                irep=nsave+1
+cd              write (iout,*) 'Accepted conformation:'
+cd              write (iout,*) (rad2deg*varia(i),i=1,nphi)
+                if (minim) call zapis(varia,etot)
+                do i=1,n_ene
+                  ener(i,nsave)=energia(i)
+                enddo
+                ener(n_ene+1,nsave)=etot
+                ener(n_ene+2,nsave)=frac
+              endif
+              nminima(irep)=nminima(irep)+1.0D0
+c             print *,'irep=',irep,' nminima=',nminima(irep)
+#ifdef MPL
+              if (Kwita.eq.0) call recv_stop_sig(kwita)
+#endif
+            endif ! accepted
+          endif ! overlap
+#ifdef MPL
+          if (MyID.eq.MasterID) then
+            call receive_MCM_info
+            if (nacc_tot.ge.maxacc) accepted=.true.
+          endif
+#endif
+          if (ntrial.gt.maxtrial_iter .and. npool.gt.0) then
+C Take a conformation from the pool
+            ii=iran_num(1,npool)
+            do i=1,nvar
+              varia(i)=xpool(i,ii)
+            enddo
+            write (iout,*) 'Iteration',it,' max. # of trials exceeded.'
+            write (iout,*) 
+     &     'Take conformation',ii,' from the pool energy=',epool(ii)
+            if (print_mc.gt.2)
+     &      write (iout,'(10f8.3)') (rad2deg*varia(i),i=1,nvar)
+            ntrial=0
+         endif ! (ntrial.gt.maxtrial_iter .and. npool.gt.0)
+   30    continue
+        enddo ! accepted
+#ifdef MPL
+        if (MyID.eq.MasterID) then
+          call receive_MCM_info
+        endif
+        if (Kwita.eq.0) call recv_stop_sig(kwita)
+#endif
+        if (ovrtim()) WhatsUp=-1
+cd      write (iout,*) 'WhatsUp=',WhatsUp,' Kwita=',Kwita
+        not_done = (nacc_tot.lt.maxacc) .and. (WhatsUp.eq.0) 
+     &         .and. (Kwita.eq.0)
+cd      write (iout,*) 'not_done=',not_done
+#ifdef MPL
+        if (Kwita.lt.0) then
+          print *,'Processor',MyID,
+     &    ' has received STOP signal =',Kwita,' in EntSamp.'
+cd        print *,'not_done=',not_done
+          if (Kwita.lt.-1) WhatsUp=Kwita
+        else if (nacc_tot.ge.maxacc) then
+          print *,'Processor',MyID,' calls send_stop_sig,',
+     &     ' because a sufficient # of confs. have been collected.'
+cd        print *,'not_done=',not_done
+          call send_stop_sig(-1)
+        else if (WhatsUp.eq.-1) then
+          print *,'Processor',MyID,
+     &               ' calls send_stop_sig because of timeout.'
+cd        print *,'not_done=',not_done
+          call send_stop_sig(-2)
+        endif
+#endif
+      enddo ! not_done
+
+C-----------------------------------------------------------------
+C... Construct energy histogram & update entropy
+C-----------------------------------------------------------------
+      go to 21
+   20 WhatsUp=-3
+#ifdef MPL
+      write (iout,*) 'Processor',MyID,
+     &       ' is broadcasting ERROR-STOP signal.'
+      write (*,*) 'Processor',MyID,
+     &       ' is broadcasting ERROR-STOP signal.'
+      call send_stop_sig(-3)
+#endif
+   21 continue
+#ifdef MPL
+      if (MyID.eq.MasterID) then
+c       call receive_MCM_results
+        call receive_energies
+#endif
+      do i=1,nsave
+        if (esave(i).lt.elowest) elowest=esave(i)
+        if (esave(i).gt.ehighest) ehighest=esave(i)
+      enddo
+      write (iout,'(a,i10)') '# of accepted confs:',nacc_tot
+      write (iout,'(a,f10.5,a,f10.5)') 'Lowest energy:',elowest,
+     & ' Highest energy',ehighest
+      if (isweep.eq.1 .and. .not.ent_read) then
+        emin=elowest
+        emax=ehighest
+        write (iout,*) 'EMAX=',emax
+        indminn=0
+        indmaxx=(ehighest-emin)/delte
+        indmin=indminn
+        indmax=indmaxx
+        do i=-max_ene,max_ene
+          entropy(i)=(emin+i*delte)*betbol
+        enddo
+        ent_read=.true.
+      else
+        indmin=(elowest-emin)/delte
+        indmax=(ehighest-emin)/delte
+        if (indmin.lt.indminn) indminn=indmin
+        if (indmax.gt.indmaxx) indmaxx=indmax
+      endif
+      write(iout,*)'indminn=',indminn,' indmaxx=',indmaxx
+C Construct energy histogram
+      do i=1,nsave
+        inde=(esave(i)-emin)/delte
+        nhist(inde)=nhist(inde)+nminima(i)
+      enddo
+C Update entropy (density of states)
+      do i=indmin,indmax
+        if (nhist(i).gt.0) then
+          entropy(i)=entropy(i)+dlog(nhist(i)+0.0D0)
+        endif
+      enddo
+Cd    do i=indmaxx+1
+Cd      entropy(i)=1.0D+10
+Cd    enddo
+      write (iout,'(/80(1h*)/a,i2/80(1h*)/)') 
+     &      'End of macroiteration',isweep
+      write (iout,'(a,f10.5,a,f10.5)') 'Elowest=',elowest,
+     &      ' Ehighest=',ehighest
+      write (iout,'(a)') 'Frequecies of minima'
+      do i=1,nsave
+        write (iout,'(i5,f5.0,f10.5)') i,nminima(i),esave(i)
+      enddo
+      write (iout,'(/a)') 'Energy histogram'
+      do i=indminn,indmaxx
+        write (iout,'(i5,2f10.5)') i,emin+i*delte,nhist(i)
+      enddo
+      write (iout,'(/a)') 'Entropy'
+      do i=indminn,indmaxx
+        write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
+      enddo
+C-----------------------------------------------------------------
+C... End of energy histogram construction
+C-----------------------------------------------------------------
+#ifdef MPL
+        entropy(-max_ene-4)=dfloat(indminn)
+        entropy(-max_ene-3)=dfloat(indmaxx)
+        entropy(-max_ene-2)=emin
+        entropy(-max_ene-1)=emax
+        call send_MCM_update
+cd      print *,entname,ientout
+        open (ientout,file=entname,status='unknown')
+        write (ientout,'(2i5,2e25.17)') indminn,indmaxx,emin,emax
+        do i=indminn,indmaxx
+          write (ientout,'(i5,f10.5,f20.15)') i,emin+i*delte,entropy(i)
+        enddo
+        close(ientout)
+      else
+        write (iout,'(a)') 'Frequecies of minima'
+        do i=1,nsave
+          write (iout,'(i5,f5.0,f10.5)') i,nminima(i),esave(i)
+        enddo
+c       call send_MCM_results
+        call send_energies
+        call receive_MCM_update
+        indminn=entropy(-max_ene-4)
+        indmaxx=entropy(-max_ene-3)
+        emin=entropy(-max_ene-2)
+        emax=entropy(-max_ene-1)
+        write (iout,*) 'Received from master:'
+        write (iout,*) 'indminn=',indminn,' indmaxx=',indmaxx,
+     &                 ' emin=',emin,' emax=',emax
+        write (iout,'(/a)') 'Entropy'
+        do i=indminn,indmaxx
+          write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
+        enddo
+      endif
+      if (WhatsUp.lt.-1) return
+#else
+      if (ovrtim() .or. WhatsUp.lt.0) return
+#endif
+
+      write (iout,'(/80(1h*)/20x,a)') 'Summary run statistics:'
+      call statprint(nacc,nfun,iretcode,etot,elowest)
+      write (iout,'(a)') 
+     & 'Statistics of multiple-bond motions. Total motions:' 
+      write (iout,'(16i5)') (nbond_move(i),i=1,Nbm)
+      write (iout,'(a)') 'Accepted motions:'
+      write (iout,'(16i5)') (nbond_acc(i),i=1,Nbm)
+      write (iout,'(a,i10)') 'Number of chain regrowths:',nregrow
+      write (iout,'(a,i10)') 'Accepted chain regrowths:',nregrow_acc
+
+C---------------------------------------------------------------------------
+      ENDDO ! ISWEEP
+C---------------------------------------------------------------------------
+
+      runtime=tcpu()
+
+      if (isweep.eq.nsweep .and. it.ge.maxacc)
+     &write (iout,'(/80(1h*)/20x,a/80(1h*)/)') 'All iterations done.'
+      return
+      end
+c------------------------------------------------------------------------------
+      subroutine accepting(ecur,eold,scur,sold,x,xold,accepted)
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.MCM'
+      include 'COMMON.MCE'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.VAR'
+#ifdef MPL
+      include 'COMMON.INFO'
+#endif
+      include 'COMMON.GEO'
+      double precision ecur,eold,xx,ran_number,bol
+      double precision x(maxvar),xold(maxvar)
+      double precision tole /1.0D-1/, tola /5.0D0/
+      logical accepted
+C Check if the conformation is similar.
+cd    write (iout,*) 'Enter ACCEPTING'
+cd    write (iout,*) 'Old PHI angles:'
+cd    write (iout,*) (rad2deg*xold(i),i=1,nphi)
+cd    write (iout,*) 'Current angles'
+cd    write (iout,*) (rad2deg*x(i),i=1,nphi)
+cd    ddif=dif_ang(nphi,x,xold)
+cd    write (iout,*) 'Angle norm:',ddif
+cd    write (iout,*) 'ecur=',ecur,' emax=',emax
+      if (ecur.gt.emax) then
+        accepted=.false.
+        if (print_mc.gt.0)
+     & write (iout,'(a)') 'Conformation rejected as too high in energy'
+        return
+      else if (dabs(ecur-eold).lt.tole .and. 
+     &       dif_ang(nphi,x,xold).lt.tola) then
+        accepted=.false.
+        if (print_mc.gt.0)
+     & write (iout,'(a)') 'Conformation rejected as too similar'
+        return
+      endif
+C Else evaluate the entropy of the conf and compare it with that of the previous
+C one.
+      indecur=(ecur-emin)/delte
+      if (iabs(indecur).gt.max_ene) then
+        write (iout,'(a,2i5)') 
+     &   'Accepting: Index out of range:',indecur
+        scur=1000.0D0 
+      else if (indecur.eq.indmaxx) then
+        scur=entropy(indecur)
+        if (print_mc.gt.0) write (iout,*)'Energy boundary reached',
+     &            indmaxx,indecur,entropy(indecur)
+      else
+        deix=ecur-(emin+indecur*delte)
+        dent=entropy(indecur+1)-entropy(indecur)
+        scur=entropy(indecur)+(dent/delte)*deix
+      endif
+cd    print *,'Processor',MyID,' ecur=',ecur,' indecur=',indecur,
+cd   & ' scur=',scur,' eold=',eold,' sold=',sold
+cd    print *,'deix=',deix,' dent=',dent,' delte=',delte
+      if (print_mc.gt.1) then
+        write(iout,*)'ecur=',ecur,' indecur=',indecur,' scur=',scur
+        write(iout,*)'eold=',eold,' sold=',sold
+      endif
+      if (scur.le.sold) then
+        accepted=.true.
+      else
+C Else carry out acceptance test
+        xx=ran_number(0.0D0,1.0D0) 
+        xxh=scur-sold
+        if (xxh.gt.50.0D0) then
+          bol=0.0D0
+        else
+          bol=exp(-xxh)
+        endif
+        if (bol.gt.xx) then
+          accepted=.true. 
+          if (print_mc.gt.0) write (iout,'(a)') 
+     &    'Conformation accepted.'
+        else
+          accepted=.false.
+          if (print_mc.gt.0) write (iout,'(a)') 
+     & 'Conformation rejected.'
+        endif
+      endif
+      return
+      end 
+c-----------------------------------------------------------------------------
+      subroutine read_pool
+      implicit real*8 (a-h,o-z)
+      include 'DIMENSIONS'
+      include 'COMMON.IOUNITS'
+      include 'COMMON.GEO'
+      include 'COMMON.MCM'
+      include 'COMMON.MCE'
+      include 'COMMON.VAR'
+      double precision varia(maxvar)
+      print '(a)','Call READ_POOL'
+      do npool=1,max_pool
+        print *,'i=',i
+        read (intin,'(i5,f10.5)',end=10,err=10) iconf,epool(npool)
+        if (epool(npool).eq.0.0D0) goto 10
+        call read_angles(intin,*10)
+        call geom_to_var(nvar,xpool(1,npool))
+      enddo
+      goto 11
+   10 npool=npool-1
+   11 write (iout,'(a,i5)') 'Number of pool conformations:',npool
+      if (print_mc.gt.2) then
+      do i=1,npool
+        write (iout,'(a,i5,a,1pe14.5)') 'Pool conformation',i,' energy',
+     &    epool(i)
+        write (iout,'(10f8.3)') (rad2deg*xpool(j,i),j=1,nvar)
+      enddo
+      endif ! (print_mc.gt.2)
+      return
+      end