source/unres/src_MD-M/mc.F

   1       subroutine monte_carlo
   2 C Does Boltzmann and entropic sampling without energy minimization
   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,nbins
  21       integer conf_comp
  22       double precision RandOrPert
  23       double precision varia(maxvar),elowest,elowest1,
  24      &                 ehighest,ehighest1,eold
  25       double precision przes(3),obr(3,3)
  26       double precision varold(maxvar)
  27       logical non_conv
  28       integer moves1(-1:MaxMoveType+1,0:MaxProcs-1),
  29      &        moves_acc1(-1:MaxMoveType+1,0:MaxProcs-1)
  30 #ifdef MPL
  31       double precision etot_temp,etot_all(0:MaxProcs)
  32       external d_vadd,d_vmin,d_vmax
  33       double precision entropy1(-max_ene:max_ene),
  34      &                 nhist1(-max_ene:max_ene)
  35       integer nbond_move1(maxres*(MaxProcs+1)),
  36      &   nbond_acc1(maxres*(MaxProcs+1)),itemp(2)
  37 #endif
  38       double precision var_lowest(maxvar)
  39       double precision energia(0:n_ene),energia_ave(0:n_ene)
  40 C
  41       write(iout,'(a,i8,2x,a,f10.5)')
  42      & 'pool_read_freq=',pool_read_freq,' pool_fraction=',pool_fraction
  43       open (istat,file=statname)
  44       WhatsUp=0
  45       indminn=-max_ene
  46       indmaxx=max_ene
  47       facee=1.0D0/(maxacc*delte)
  48 C Number of bins in energy histogram
  49       nbins=e_up/delte-1
  50       write (iout,*) 'NBINS=',nbins
  51       conste=dlog(facee)
  52 C Read entropy from previous simulations.
  53       if (ent_read) then
  54         read (ientin,*) indminn,indmaxx,emin,emax
  55         print *,'indminn=',indminn,' indmaxx=',indmaxx,' emin=',emin,
  56      &          ' emax=',emax
  57         do i=-max_ene,max_ene
  58           entropy(i)=0.0D0
  59         enddo
  60         read (ientin,*) (ijunk,ejunk,entropy(i),i=indminn,indmaxx)
  61         indmin=indminn
  62         indmax=indmaxx
  63         write (iout,*) 'indminn=',indminn,' indmaxx=',indmaxx,
  64      &                 ' emin=',emin,' emax=',emax
  65         write (iout,'(/a)') 'Initial entropy'
  66         do i=indminn,indmaxx
  67           write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
  68         enddo
  69       endif ! ent_read
  70 C Read the pool of conformations
  71       call read_pool
  72       elowest=1.0D+10
  73       ehighest=-1.0D+10
  74 C----------------------------------------------------------------------------
  75 C Entropy-sampling simulations with continually updated entropy;
  76 C set NSWEEP=1 for Boltzmann sampling.
  77 C Loop thru simulations
  78 C----------------------------------------------------------------------------
  79       DO ISWEEP=1,NSWEEP
  80 C
  81 C Initialize the IFINISH array.
  82 C
  83 #ifdef MPL
  84         do i=1,nctasks
  85           ifinish(i)=0
  86         enddo
  87 #endif
  88 c---------------------------------------------------------------------------
  89 C Initialize counters.
  90 c---------------------------------------------------------------------------
  91 C Total number of generated confs.
  92         ngen=0
  93 C Total number of moves. In general this won't be equal to the number of
  94 C attempted moves, because we may want to reject some "bad" confs just by
  95 C overlap check.
  96         nmove=0
  97 C Total number of shift (nbond_move(1)), spike, crankshaft, three-bond,...
  98 C motions.
  99         do i=1,nres
 100           nbond_move(i)=0
 101           nbond_acc(i)=0
 102         enddo
 103 C Initialize total and accepted number of moves of various kind.
 104         do i=-1,MaxMoveType
 105           moves(i)=0
 106           moves_acc(i)=0
 107         enddo
 108 C Total number of energy evaluations.
 109         neneval=0
 110         nfun=0
 111 C----------------------------------------------------------------------------
 112 C Take a conformation from the pool
 113 C----------------------------------------------------------------------------
 114       rewind(istat)
 115       write (iout,*) 'emin=',emin,' emax=',emax
 116       if (npool.gt.0) then
 117         ii=iran_num(1,npool)
 118         do i=1,nvar
 119           varia(i)=xpool(i,ii)
 120         enddo
 121         write (iout,*) 'Took conformation',ii,' from the pool energy=',
 122      &               epool(ii)
 123         call var_to_geom(nvar,varia)
 124 C Print internal coordinates of the initial conformation
 125         call intout
 126       else if (isweep.gt.1) then
 127         if (eold.lt.emax) then
 128         do i=1,nvar
 129           varia(i)=varold(i)
 130         enddo
 131         else
 132         do i=1,nvar
 133           varia(i)=var_lowest(i)
 134         enddo
 135         endif
 136         call var_to_geom(nvar,varia)
 137       endif
 138 C----------------------------------------------------------------------------
 139 C Compute and print initial energies.
 140 C----------------------------------------------------------------------------
 141       nsave=0
 142       Kwita=0
 143       WhatsUp=0
 144       write (iout,'(/80(1h*)/a,i2/80(1h*)/)') 'MCE iteration #',isweep
 145       write (iout,'(/80(1h*)/a)') 'Initial energies:'
 146       call chainbuild
 147       call geom_to_var(nvar,varia)
 148       call etotal(energia(0))
 149       etot = energia(0)
 150       call enerprint(energia(0))
 151       if (refstr) then
 152         call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,1),nsup,przes,
 153      &             obr,non_conv)
 154         rms=dsqrt(rms)
 155         call contact(.false.,ncont,icont,co)
 156         frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
 157         write (iout,'(a,f8.3,a,f8.3,a,f8.3)')
 158      &    'RMS deviation from the reference structure:',rms,
 159      &    ' % of native contacts:',frac*100,' contact order',co
 160         write (istat,'(i10,16(1pe14.5))') 0,
 161      &   (energia(print_order(i)),i=1,nprint_ene),
 162      &   etot,rms,frac,co
 163       else
 164         write (istat,'(i10,14(1pe14.5))') 0,
 165      &   (energia(print_order(i)),i=1,nprint_ene),etot
 166       endif
 167 c     close(istat)
 168       neneval=neneval+1
 169       if (.not. ent_read) then
 170 C Initialize the entropy array
 171 #ifdef MPL
 172 C Collect total energies from other processors.
 173         etot_temp=etot
 174         etot_all(0)=etot
 175         call mp_gather(etot_temp,etot_all,8,MasterID,cgGroupID)
 176         if (MyID.eq.MasterID) then
 177 C Get the lowest and the highest energy.
 178           print *,'MASTER: etot_temp: ',(etot_all(i),i=0,nprocs-1),
 179      &     ' emin=',emin,' emax=',emax
 180           emin=1.0D10
 181           emax=-1.0D10
 182           do i=0,nprocs
 183             if (emin.gt.etot_all(i)) emin=etot_all(i)
 184             if (emax.lt.etot_all(i)) emax=etot_all(i)
 185           enddo
 186           emax=emin+e_up
 187         endif ! MyID.eq.MasterID
 188         etot_all(1)=emin
 189         etot_all(2)=emax
 190         print *,'Processor',MyID,' calls MP_BCAST to send/recv etot_all'
 191         call mp_bcast(etot_all(1),16,MasterID,cgGroupID)
 192         print *,'Processor',MyID,' MP_BCAST to send/recv etot_all ended'
 193         if (MyID.ne.MasterID) then
 194           print *,'Processor:',MyID,etot_all(1),etot_all(2),
 195      &          etot_all(1),etot_all(2)
 196           emin=etot_all(1)
 197           emax=etot_all(2)
 198         endif ! MyID.ne.MasterID
 199         write (iout,*) 'After MP_GATHER etot_temp=',
 200      &                 etot_temp,' emin=',emin
 201 #else
 202         emin=etot
 203         emax=emin+e_up
 204         indminn=0
 205         indmin=0
 206 #endif
 207         IF (MULTICAN) THEN
 208 C Multicanonical sampling - start from Boltzmann distribution
 209           do i=-max_ene,max_ene
 210             entropy(i)=(emin+i*delte)*betbol
 211           enddo
 212         ELSE
 213 C Entropic sampling - start from uniform distribution of the density of states
 214           do i=-max_ene,max_ene
 215             entropy(i)=0.0D0
 216           enddo
 217         ENDIF ! MULTICAN
 218         write (iout,'(/a)') 'Initial entropy'
 219         do i=indminn,indmaxx
 220           write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
 221         enddo
 222         if (isweep.eq.1) then
 223           emax=emin+e_up
 224           indminn=0
 225           indmin=0
 226           indmaxx=indminn+nbins
 227           indmax=indmaxx
 228         endif ! isweep.eq.1
 229       endif ! .not. ent_read
 230 #ifdef MPL
 231       call recv_stop_sig(Kwita)
 232       if (whatsup.eq.1) then
 233         call send_stop_sig(-2)
 234         not_done=.false.
 235       else if (whatsup.le.-2) then
 236         not_done=.false.
 237       else if (whatsup.eq.2) then
 238         not_done=.false.
 239       else
 240         not_done=.true.
 241       endif
 242 #else
 243       not_done=.true.
 244 #endif
 245       write (iout,'(/80(1h*)/20x,a/80(1h*))')
 246      &    'Enter Monte Carlo procedure.'
 247       close(igeom)
 248       call briefout(0,etot)
 249       do i=1,nvar
 250         varold(i)=varia(i)
 251       enddo
 252       eold=etot
 253       call entropia(eold,sold,indeold)
 254 C NACC is the counter for the accepted conformations of a given processor
 255       nacc=0
 256 C NACC_TOT counts the total number of accepted conformations
 257       nacc_tot=0
 258 C Main loop.
 259 c----------------------------------------------------------------------------
 260 C Zero out average energies
 261       do i=0,n_ene
 262         energia_ave(i)=0.0d0
 263       enddo
 264 C Initialize energy histogram
 265       do i=-max_ene,max_ene
 266         nhist(i)=0.0D0
 267       enddo   ! i
 268 C Zero out iteration counter.
 269       it=0
 270       do j=1,nvar
 271         varold(j)=varia(j)
 272        enddo
 273 C Begin MC iteration loop.
 274       do while (not_done)
 275         it=it+1
 276 C Initialize local counter.
 277         ntrial=0 ! # of generated non-overlapping confs.
 278         noverlap=0 ! # of overlapping confs.
 279         accepted=.false.
 280         do while (.not. accepted .and. WhatsUp.eq.0 .and. Kwita.eq.0)
 281           ntrial=ntrial+1
 282 C Retrieve the angles of previously accepted conformation
 283           do j=1,nvar
 284             varia(j)=varold(j)
 285           enddo
 286           call var_to_geom(nvar,varia)
 287 C Rebuild the chain.
 288           call chainbuild
 289           MoveType=0
 290           nbond=0
 291           lprint=.true.
 292 C Decide whether to take a conformation from the pool or generate/perturb one
 293 C randomly
 294           from_pool=ran_number(0.0D0,1.0D0)
 295           if (npool.gt.0 .and. from_pool.lt.pool_fraction) then
 296 C Throw a dice to choose the conformation from the pool
 297             ii=iran_num(1,npool)
 298             do i=1,nvar
 299               varia(i)=xpool(i,ii)
 300             enddo
 301             call var_to_geom(nvar,varia)
 302             call chainbuild
 303 cd          call intout
 304 cd          write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
 305             if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it)
 306      &        write (iout,'(a,i3,a,f10.5)')
 307      &     'Try conformation',ii,' from the pool energy=',epool(ii)
 308             MoveType=-1
 309             moves(-1)=moves(-1)+1
 310           else
 311 C Decide whether to generate a random conformation or perturb the old one
 312           RandOrPert=ran_number(0.0D0,1.0D0)
 313           if (RandOrPert.gt.RanFract) then
 314             if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it)
 315      &        write (iout,'(a)') 'Perturbation-generated conformation.'
 316             call perturb(error,lprint,MoveType,nbond,0.1D0)
 317             if (error) goto 20
 318             if (MoveType.lt.1 .or. MoveType.gt.MaxMoveType) then
 319               write (iout,'(/a,i7,a/)') 'Error - unknown MoveType=',
 320      &           MoveType,' returned from PERTURB.'
 321               goto 20
 322             endif
 323             call chainbuild
 324           else
 325             MoveType=0
 326             moves(0)=moves(0)+1
 327             nstart_grow=iran_num(3,nres)
 328             if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it)
 329      &        write (iout,'(2a,i3)') 'Random-generated conformation',
 330      &        ' - chain regrown from residue',nstart_grow
 331             nrestmp=nres
 332             call gen_rand_conf(nstart_grow,nrestmp,*30)
 333           endif
 334           call geom_to_var(nvar,varia)
 335           endif ! pool
 336 Cd        write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
 337           ngen=ngen+1
 338           if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it)
 339      &      write (iout,'(a,i5,a,i10,a,i10)')
 340      &   'Processor',MyId,' trial move',ntrial,' total generated:',ngen
 341           if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it)
 342      &      write (*,'(a,i5,a,i10,a,i10)')
 343      &   'Processor',MyId,' trial move',ntrial,' total generated:',ngen
 344           call etotal(energia(0))
 345           etot = energia(0)
 346           neneval=neneval+1
 347 cd        call enerprint(energia(0))
 348 cd        write(iout,*)'it=',it,' etot=',etot
 349           if (etot-elowest.gt.overlap_cut) then
 350             if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it)
 351      &        write (iout,'(a,i5,a,1pe14.5)')  'Iteration',it,
 352      &       ' Overlap detected in the current conf.; energy is',etot
 353             accepted=.false.
 354             noverlap=noverlap+1
 355             if (noverlap.gt.maxoverlap) then
 356               write (iout,'(a)') 'Too many overlapping confs.'
 357               goto 20
 358             endif
 359           else
 360 C--------------------------------------------------------------------------
 361 C... Acceptance test
 362 C--------------------------------------------------------------------------
 363             accepted=.false.
 364             if (WhatsUp.eq.0)
 365      &      call accept_mc(it,etot,eold,scur,sold,varia,varold,accepted)
 366             if (accepted) then
 367               nacc=nacc+1
 368               nacc_tot=nacc_tot+1
 369               if (elowest.gt.etot) then
 370                 elowest=etot
 371                 do i=1,nvar
 372                   var_lowest(i)=varia(i)
 373                 enddo
 374               endif
 375               if (ehighest.lt.etot) ehighest=etot
 376               moves_acc(MoveType)=moves_acc(MoveType)+1
 377               if (MoveType.eq.1) then
 378                 nbond_acc(nbond)=nbond_acc(nbond)+1
 379               endif
 380 C Compare with reference structure.
 381               if (refstr) then
 382                 call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,1),
 383      &                      nsup,przes,obr,non_conv)
 384                 rms=dsqrt(rms)
 385                 call contact(.false.,ncont,icont,co)
 386                 frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
 387               endif ! refstr
 388 C
 389 C Periodically save average energies and confs.
 390 C
 391               do i=0,n_ene
 392                 energia_ave(i)=energia_ave(i)+energia(i)
 393               enddo
 394               moves(MaxMoveType+1)=nmove
 395               moves_acc(MaxMoveType+1)=nacc
 396               IF ((it/save_frequency)*save_frequency.eq.it) THEN
 397                 do i=0,n_ene
 398                   energia_ave(i)=energia_ave(i)/save_frequency
 399                 enddo
 400                 etot_ave=energia_ave(0)
 401 C#ifdef AIX
 402 C                open (istat,file=statname,position='append')
 403 C#else
 404 C                open (istat,file=statname,access='append')
 405 Cendif
 406                 if (print_mc.gt.0)
 407      &            write (iout,'(80(1h*)/20x,a,i20)')
 408      &                             'Iteration #',it
 409                 if (refstr .and. print_mc.gt.0)  then
 410                   write (iout,'(a,f8.3,a,f8.3,a,f8.3)')
 411      &            'RMS deviation from the reference structure:',rms,
 412      &            ' % of native contacts:',frac*100,' contact order:',co
 413                 endif
 414                 if (print_stat) then
 415                   if (refstr) then
 416                     write (istat,'(i10,10(1pe14.5))') it,
 417      &              (energia_ave(print_order(i)),i=1,nprint_ene),
 418      &                etot_ave,rms_ave,frac_ave
 419                   else
 420                     write (istat,'(i10,10(1pe14.5))') it,
 421      &              (energia_ave(print_order(i)),i=1,nprint_ene),
 422      &               etot_ave
 423                   endif
 424                 endif
 425 c               close(istat)
 426                 if (print_mc.gt.0)
 427      &            call statprint(nacc,nfun,iretcode,etot,elowest)
 428 C Print internal coordinates.
 429                 if (print_int) call briefout(nacc,etot)
 430                 do i=0,n_ene
 431                   energia_ave(i)=0.0d0
 432                 enddo
 433               ENDIF ! ( (it/save_frequency)*save_frequency.eq.it)
 434 C Update histogram
 435               inde=icialosc((etot-emin)/delte)
 436               nhist(inde)=nhist(inde)+1.0D0
 437 #ifdef MPL
 438               if ( (it/message_frequency)*message_frequency.eq.it
 439      &                              .and. (MyID.ne.MasterID) ) then
 440                 call recv_stop_sig(Kwita)
 441                 call send_MCM_info(message_frequency)
 442               endif
 443 #endif
 444 C Store the accepted conf. and its energy.
 445               eold=etot
 446               sold=scur
 447               do i=1,nvar
 448                 varold(i)=varia(i)
 449               enddo
 450 #ifdef MPL
 451               if (Kwita.eq.0) call recv_stop_sig(kwita)
 452 #endif
 453             endif ! accepted
 454           endif ! overlap
 455 #ifdef MPL
 456           if (MyID.eq.MasterID .and.
 457      &        (it/message_frequency)*message_frequency.eq.it) then
 458             call receive_MC_info
 459             if (nacc_tot.ge.maxacc) accepted=.true.
 460           endif
 461 #endif
 462 C         if ((ntrial.gt.maxtrial_iter
 463 C    &       .or. (it/pool_read_freq)*pool_read_freq.eq.it)
 464 C    &       .and. npool.gt.0) then
 465 C Take a conformation from the pool
 466 C           ii=iran_num(1,npool)
 467 C           do i=1,nvar
 468 C             varold(i)=xpool(i,ii)
 469 C           enddo
 470 C           if (ntrial.gt.maxtrial_iter)
 471 C    &      write (iout,*) 'Iteration',it,' max. # of trials exceeded.'
 472 C           write (iout,*)
 473 C    &     'Take conformation',ii,' from the pool energy=',epool(ii)
 474 C           if (print_mc.gt.2)
 475 C    &      write (iout,'(10f8.3)') (rad2deg*varold(i),i=1,nvar)
 476 C           ntrial=0
 477 C           eold=epool(ii)
 478 C           call entropia(eold,sold,indeold)
 479 C           accepted=.true.
 480 C        endif ! (ntrial.gt.maxtrial_iter .and. npool.gt.0)
 481    30    continue
 482         enddo ! accepted
 483 #ifdef MPL
 484         if (MyID.eq.MasterID .and.
 485      &      (it/message_frequency)*message_frequency.eq.it) then
 486           call receive_MC_info
 487         endif
 488         if (Kwita.eq.0) call recv_stop_sig(kwita)
 489 #endif
 490         if (ovrtim()) WhatsUp=-1
 491 cd      write (iout,*) 'WhatsUp=',WhatsUp,' Kwita=',Kwita
 492         not_done = (nacc_tot.lt.maxacc) .and. (WhatsUp.eq.0)
 493      &         .and. (Kwita.eq.0)
 494 cd      write (iout,*) 'not_done=',not_done
 495 #ifdef MPL
 496         if (Kwita.lt.0) then
 497           print *,'Processor',MyID,
 498      &    ' has received STOP signal =',Kwita,' in EntSamp.'
 499 cd        print *,'not_done=',not_done
 500           if (Kwita.lt.-1) WhatsUp=Kwita
 501           if (MyID.ne.MasterID) call send_MCM_info(-1)
 502         else if (nacc_tot.ge.maxacc) then
 503           print *,'Processor',MyID,' calls send_stop_sig,',
 504      &     ' because a sufficient # of confs. have been collected.'
 505 cd        print *,'not_done=',not_done
 506           call send_stop_sig(-1)
 507           if (MyID.ne.MasterID) call send_MCM_info(-1)
 508         else if (WhatsUp.eq.-1) then
 509           print *,'Processor',MyID,
 510      &               ' calls send_stop_sig because of timeout.'
 511 cd        print *,'not_done=',not_done
 512           call send_stop_sig(-2)
 513           if (MyID.ne.MasterID) call send_MCM_info(-1)
 514         endif
 515 #endif
 516       enddo ! not_done
 517
 518 C-----------------------------------------------------------------
 519 C... Construct energy histogram & update entropy
 520 C-----------------------------------------------------------------
 521       go to 21
 522    20 WhatsUp=-3
 523 #ifdef MPL
 524       write (iout,*) 'Processor',MyID,
 525      &       ' is broadcasting ERROR-STOP signal.'
 526       write (*,*) 'Processor',MyID,
 527      &       ' is broadcasting ERROR-STOP signal.'
 528       call send_stop_sig(-3)
 529       if (MyID.ne.MasterID) call send_MCM_info(-1)
 530 #endif
 531    21 continue
 532       write (iout,'(/a)') 'Energy histogram'
 533       do i=-100,100
 534         write (iout,'(i5,2f20.5)') i,emin+i*delte,nhist(i)
 535       enddo
 536 #ifdef MPL
 537 C Wait until every processor has sent complete MC info.
 538       if (MyID.eq.MasterID) then
 539         not_done=.true.
 540         do while (not_done)
 541 C         write (*,*) 'The IFINISH array:'
 542 C         write (*,*) (ifinish(i),i=1,nctasks)
 543           not_done=.false.
 544           do i=2,nctasks
 545             not_done=not_done.or.(ifinish(i).ge.0)
 546           enddo
 547           if (not_done) call receive_MC_info
 548         enddo
 549       endif
 550 C Make collective histogram from the work of all processors.
 551       msglen=(2*max_ene+1)*8
 552       print *,
 553      & 'Processor',MyID,' calls MP_REDUCE to send/receive histograms',
 554      & ' msglen=',msglen
 555       call mp_reduce(nhist,nhist1,msglen,MasterID,d_vadd,
 556      &               cgGroupID)
 557       print *,'Processor',MyID,' MP_REDUCE accomplished for histogr.'
 558       do i=-max_ene,max_ene
 559         nhist(i)=nhist1(i)
 560       enddo
 561 C Collect min. and max. energy
 562       print *,
 563      &'Processor',MyID,' calls MP_REDUCE to send/receive energy borders'
 564       call mp_reduce(elowest,elowest1,8,MasterID,d_vmin,cgGroupID)
 565       call mp_reduce(ehighest,ehighest1,8,MasterID,d_vmax,cgGroupID)
 566       print *,'Processor',MyID,' MP_REDUCE accomplished for energies.'
 567       IF (MyID.eq.MasterID) THEN
 568         elowest=elowest1
 569         ehighest=ehighest1
 570 #endif
 571         write (iout,'(a,i10)') '# of accepted confs:',nacc_tot
 572         write (iout,'(a,f10.5,a,f10.5)') 'Lowest energy:',elowest,
 573      & ' Highest energy',ehighest
 574         indmin=icialosc((elowest-emin)/delte)
 575         imdmax=icialosc((ehighest-emin)/delte)
 576         if (indmin.lt.indminn) then
 577           emax=emin+indmin*delte+e_up
 578           indmaxx=indmin+nbins
 579           indminn=indmin
 580         endif
 581         if (.not.ent_read) ent_read=.true.
 582         write(iout,*)'indminn=',indminn,' indmaxx=',indmaxx
 583 C Update entropy (density of states)
 584         do i=indmin,indmax
 585           if (nhist(i).gt.0) then
 586             entropy(i)=entropy(i)+dlog(nhist(i)+0.0D0)
 587           endif
 588         enddo
 589         write (iout,'(/80(1h*)/a,i2/80(1h*)/)')
 590      &        'End of macroiteration',isweep
 591         write (iout,'(a,f10.5,a,f10.5)') 'Elowest=',elowest,
 592      &      ' Ehighest=',ehighest
 593         write (iout,'(/a)') 'Energy histogram'
 594         do i=indminn,indmaxx
 595           write (iout,'(i5,2f20.5)') i,emin+i*delte,nhist(i)
 596         enddo
 597         write (iout,'(/a)') 'Entropy'
 598         do i=indminn,indmaxx
 599           write (iout,'(i5,2f20.5)') i,emin+i*delte,entropy(i)
 600         enddo
 601 C-----------------------------------------------------------------
 602 C... End of energy histogram construction
 603 C-----------------------------------------------------------------
 604 #ifdef MPL
 605       ELSE
 606         if (.not. ent_read) ent_read=.true.
 607       ENDIF ! MyID .eq. MaterID
 608       if (MyID.eq.MasterID) then
 609         itemp(1)=indminn
 610         itemp(2)=indmaxx
 611       endif
 612       print *,'before mp_bcast processor',MyID,' indminn=',indminn,
 613      & ' indmaxx=',indmaxx,' itemp=',itemp(1),itemp(2)
 614       call mp_bcast(itemp(1),8,MasterID,cgGroupID)
 615       call mp_bcast(emax,8,MasterID,cgGroupID)
 616       print *,'after mp_bcast processor',MyID,' indminn=',indminn,
 617      & ' indmaxx=',indmaxx,' itemp=',itemp(1),itemp(2)
 618       if (MyID .ne. MasterID) then
 619         indminn=itemp(1)
 620         indmaxx=itemp(2)
 621       endif
 622       msglen=(indmaxx-indminn+1)*8
 623       print *,'processor',MyID,' calling mp_bcast msglen=',msglen,
 624      & ' indminn=',indminn,' indmaxx=',indmaxx,' isweep=',isweep
 625       call mp_bcast(entropy(indminn),msglen,MasterID,cgGroupID)
 626       IF(MyID.eq.MasterID .and. .not. ovrtim() .and. WhatsUp.ge.0)THEN
 627         open (ientout,file=entname,status='unknown')
 628         write (ientout,'(2i5,2e25.17)') indminn,indmaxx,emin,emax
 629         do i=indminn,indmaxx
 630           write (ientout,'(i5,f10.5,f20.15)') i,emin+i*delte,entropy(i)
 631         enddo
 632         close(ientout)
 633       ELSE
 634         write (iout,*) 'Received from master:'
 635         write (iout,*) 'indminn=',indminn,' indmaxx=',indmaxx,
 636      &                 ' emin=',emin,' emax=',emax
 637         write (iout,'(/a)') 'Entropy'
 638         do i=indminn,indmaxx
 639            write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
 640         enddo
 641       ENDIF ! MyID.eq.MasterID
 642       print *,'Processor',MyID,' calls MP_GATHER'
 643       call mp_gather(nbond_move,nbond_move1,4*Nbm,MasterID,
 644      &               cgGroupID)
 645       call mp_gather(nbond_acc,nbond_acc1,4*Nbm,MasterID,
 646      &               cgGroupID)
 647       print *,'Processor',MyID,' MP_GATHER call accomplished'
 648       if (MyID.eq.MasterID) then
 649
 650         write (iout,'(/80(1h*)/20x,a)') 'Summary run statistics:'
 651         call statprint(nacc_tot,nfun,iretcode,etot,elowest)
 652         write (iout,'(a)')
 653      &   'Statistics of multiple-bond motions. Total motions:'
 654         write (iout,'(8i10)') (nbond_move(i),i=1,Nbm)
 655         write (iout,'(a)') 'Accepted motions:'
 656         write (iout,'(8i10)') (nbond_acc(i),i=1,Nbm)
 657
 658         write (iout,'(a)')
 659      & 'Statistics of multi-bond moves of respective processors:'
 660         do iproc=1,Nprocs-1
 661           do i=1,Nbm
 662             ind=iproc*nbm+i
 663             nbond_move(i)=nbond_move(i)+nbond_move1(ind)
 664             nbond_acc(i)=nbond_acc(i)+nbond_acc1(ind)
 665           enddo
 666         enddo
 667         do iproc=0,NProcs-1
 668           write (iout,*) 'Processor',iproc,' nbond_move:',
 669      &        (nbond_move1(iproc*nbm+i),i=1,Nbm),
 670      &        ' nbond_acc:',(nbond_acc1(iproc*nbm+i),i=1,Nbm)
 671         enddo
 672       endif
 673       call mp_gather(moves,moves1,4*(MaxMoveType+3),MasterID,
 674      &               cgGroupID)
 675       call mp_gather(moves_acc,moves_acc1,4*(MaxMoveType+3),
 676      &               MasterID,cgGroupID)
 677       if (MyID.eq.MasterID) then
 678         do iproc=1,Nprocs-1
 679           do i=-1,MaxMoveType+1
 680             moves(i)=moves(i)+moves1(i,iproc)
 681             moves_acc(i)=moves_acc(i)+moves_acc1(i,iproc)
 682           enddo
 683         enddo
 684         nmove=0
 685         do i=0,MaxMoveType+1
 686           nmove=nmove+moves(i)
 687         enddo
 688         do iproc=0,NProcs-1
 689           write (iout,*) 'Processor',iproc,' moves',
 690      &     (moves1(i,iproc),i=0,MaxMoveType+1),
 691      &    ' moves_acc:',(moves_acc1(i,iproc),i=0,MaxMoveType+1)
 692         enddo
 693       endif
 694 #else
 695       open (ientout,file=entname,status='unknown')
 696       write (ientout,'(2i5,2e25.17)') indminn,indmaxx,emin,emax
 697       do i=indminn,indmaxx
 698         write (ientout,'(i5,f10.5,f20.15)') i,emin+i*delte,entropy(i)
 699       enddo
 700       close(ientout)
 701 #endif
 702       write (iout,'(/80(1h*)/20x,a)') 'Summary run statistics:'
 703       call statprint(nacc_tot,nfun,iretcode,etot,elowest)
 704       write (iout,'(a)')
 705      & 'Statistics of multiple-bond motions. Total motions:'
 706       write (iout,'(8i10)') (nbond_move(i),i=1,Nbm)
 707       write (iout,'(a)') 'Accepted motions:'
 708       write (iout,'(8i10)') (nbond_acc(i),i=1,Nbm)
 709       if (ovrtim() .or. WhatsUp.lt.0) return
 710
 711 C---------------------------------------------------------------------------
 712       ENDDO ! ISWEEP
 713 C---------------------------------------------------------------------------
 714
 715       runtime=tcpu()
 716
 717       if (isweep.eq.nsweep .and. it.ge.maxacc)
 718      &write (iout,'(/80(1h*)/20x,a/80(1h*)/)') 'All iterations done.'
 719       return
 720       end
 721 c------------------------------------------------------------------------------
 722       subroutine accept_mc(it,ecur,eold,scur,sold,x,xold,accepted)
 723       implicit real*8 (a-h,o-z)
 724       include 'DIMENSIONS'
 725       include 'COMMON.MCM'
 726       include 'COMMON.MCE'
 727       include 'COMMON.IOUNITS'
 728       include 'COMMON.VAR'
 729 #ifdef MPL
 730       include 'COMMON.INFO'
 731 #endif
 732       include 'COMMON.GEO'
 733       double precision ecur,eold,xx,ran_number,bol
 734       double precision x(maxvar),xold(maxvar)
 735       logical accepted
 736 C Check if the conformation is similar.
 737 cd    write (iout,*) 'Enter ACCEPTING'
 738 cd    write (iout,*) 'Old PHI angles:'
 739 cd    write (iout,*) (rad2deg*xold(i),i=1,nphi)
 740 cd    write (iout,*) 'Current angles'
 741 cd    write (iout,*) (rad2deg*x(i),i=1,nphi)
 742 cd    ddif=dif_ang(nphi,x,xold)
 743 cd    write (iout,*) 'Angle norm:',ddif
 744 cd    write (iout,*) 'ecur=',ecur,' emax=',emax
 745       if (ecur.gt.emax) then
 746         accepted=.false.
 747         if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it)
 748      & write (iout,'(a)') 'Conformation rejected as too high in energy'
 749         return
 750       endif
 751 C Else evaluate the entropy of the conf and compare it with that of the previous
 752 C one.
 753       call entropia(ecur,scur,indecur)
 754 cd    print *,'Processor',MyID,' ecur=',ecur,' indecur=',indecur,
 755 cd   & ' scur=',scur,' eold=',eold,' sold=',sold
 756 cd    print *,'deix=',deix,' dent=',dent,' delte=',delte
 757       if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) then
 758         write(iout,*)'it=',it,'ecur=',ecur,' indecur=',indecur,
 759      &   ' scur=',scur
 760         write(iout,*)'eold=',eold,' sold=',sold
 761       endif
 762       if (scur.le.sold) then
 763         accepted=.true.
 764       else
 765 C Else carry out acceptance test
 766         xx=ran_number(0.0D0,1.0D0)
 767         xxh=scur-sold
 768         if (xxh.gt.50.0D0) then
 769           bol=0.0D0
 770         else
 771           bol=exp(-xxh)
 772         endif
 773         if (bol.gt.xx) then
 774           accepted=.true.
 775           if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it)
 776      &       write (iout,'(a)') 'Conformation accepted.'
 777         else
 778           accepted=.false.
 779           if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it)
 780      &       write (iout,'(a)') 'Conformation rejected.'
 781         endif
 782       endif
 783       return
 784       end
 785 c--------------------------------------------------------------------------
 786       integer function icialosc(x)
 787       double precision x
 788       if (x.lt.0.0D0) then
 789         icialosc=dint(x)-1
 790       else
 791         icialosc=dint(x)
 792       endif
 793       return
 794       end
 795 c--------------------------------------------------------------------------
 796       subroutine entropia(ecur,scur,indecur)
 797       implicit real*8 (a-h,o-z)
 798       include 'DIMENSIONS'
 799       include 'COMMON.MCM'
 800       include 'COMMON.MCE'
 801       include 'COMMON.IOUNITS'
 802       double precision ecur,scur
 803       integer indecur
 804       indecur=icialosc((ecur-emin)/delte)
 805       if (iabs(indecur).gt.max_ene) then
 806         if ((it/print_freq)*it.eq.it) write (iout,'(a,2i5)')
 807      &   'Accepting: Index out of range:',indecur
 808         scur=1000.0D0
 809       else if (indecur.ge.indmaxx) then
 810         scur=entropy(indecur)
 811         if (print_mc.gt.0 .and. (it/print_freq)*it.eq.it)
 812      &    write (iout,*)'Energy boundary reached',
 813      &            indmaxx,indecur,entropy(indecur)
 814       else
 815         deix=ecur-(emin+indecur*delte)
 816         dent=entropy(indecur+1)-entropy(indecur)
 817         scur=entropy(indecur)+(dent/delte)*deix
 818       endif
 819       return
 820       end