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