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