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