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