source/unres/src_MD/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),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             call gen_rand_conf(nstart_grow,*30)
 332           endif
 333           call geom_to_var(nvar,varia)
 334           endif ! pool
 335 Cd        write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
 336           ngen=ngen+1
 337           if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it)
 338      &      write (iout,'(a,i5,a,i10,a,i10)')
 339      &   'Processor',MyId,' trial move',ntrial,' total generated:',ngen
 340           if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it)
 341      &      write (*,'(a,i5,a,i10,a,i10)')
 342      &   'Processor',MyId,' trial move',ntrial,' total generated:',ngen
 343           call etotal(energia(0))
 344           etot = energia(0)
 345           neneval=neneval+1
 346 cd        call enerprint(energia(0))
 347 cd        write(iout,*)'it=',it,' etot=',etot
 348           if (etot-elowest.gt.overlap_cut) then
 349             if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it)
 350      &        write (iout,'(a,i5,a,1pe14.5)')  'Iteration',it,
 351      &       ' Overlap detected in the current conf.; energy is',etot
 352             accepted=.false.
 353             noverlap=noverlap+1
 354             if (noverlap.gt.maxoverlap) then
 355               write (iout,'(a)') 'Too many overlapping confs.'
 356               goto 20
 357             endif
 358           else
 359 C--------------------------------------------------------------------------
 360 C... Acceptance test
 361 C--------------------------------------------------------------------------
 362             accepted=.false.
 363             if (WhatsUp.eq.0)
 364      &      call accept_mc(it,etot,eold,scur,sold,varia,varold,accepted)
 365             if (accepted) then
 366               nacc=nacc+1
 367               nacc_tot=nacc_tot+1
 368               if (elowest.gt.etot) then
 369                 elowest=etot
 370                 do i=1,nvar
 371                   var_lowest(i)=varia(i)
 372                 enddo
 373               endif
 374               if (ehighest.lt.etot) ehighest=etot
 375               moves_acc(MoveType)=moves_acc(MoveType)+1
 376               if (MoveType.eq.1) then
 377                 nbond_acc(nbond)=nbond_acc(nbond)+1
 378               endif
 379 C Compare with reference structure.
 380               if (refstr) then
 381                 call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup),
 382      &                      nsup,przes,obr,non_conv)
 383                 rms=dsqrt(rms)
 384                 call contact(.false.,ncont,icont,co)
 385                 frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
 386               endif ! refstr
 387 C
 388 C Periodically save average energies and confs.
 389 C
 390               do i=0,n_ene
 391                 energia_ave(i)=energia_ave(i)+energia(i)
 392               enddo
 393               moves(MaxMoveType+1)=nmove
 394               moves_acc(MaxMoveType+1)=nacc
 395               IF ((it/save_frequency)*save_frequency.eq.it) THEN
 396                 do i=0,n_ene
 397                   energia_ave(i)=energia_ave(i)/save_frequency
 398                 enddo
 399                 etot_ave=energia_ave(0)
 400 C#ifdef AIX
 401 C                open (istat,file=statname,position='append')
 402 C#else
 403 C                open (istat,file=statname,access='append')
 404 Cendif
 405                 if (print_mc.gt.0)
 406      &            write (iout,'(80(1h*)/20x,a,i20)')
 407      &                             'Iteration #',it
 408                 if (refstr .and. print_mc.gt.0)  then
 409                   write (iout,'(a,f8.3,a,f8.3,a,f8.3)')
 410      &            'RMS deviation from the reference structure:',rms,
 411      &            ' % of native contacts:',frac*100,' contact order:',co
 412                 endif
 413                 if (print_stat) then
 414                   if (refstr) then
 415                     write (istat,'(i10,10(1pe14.5))') it,
 416      &              (energia_ave(print_order(i)),i=1,nprint_ene),
 417      &                etot_ave,rms_ave,frac_ave
 418                   else
 419                     write (istat,'(i10,10(1pe14.5))') it,
 420      &              (energia_ave(print_order(i)),i=1,nprint_ene),
 421      &               etot_ave
 422                   endif
 423                 endif
 424 c               close(istat)
 425                 if (print_mc.gt.0)
 426      &            call statprint(nacc,nfun,iretcode,etot,elowest)
 427 C Print internal coordinates.
 428                 if (print_int) call briefout(nacc,etot)
 429                 do i=0,n_ene
 430                   energia_ave(i)=0.0d0
 431                 enddo
 432               ENDIF ! ( (it/save_frequency)*save_frequency.eq.it)
 433 C Update histogram
 434               inde=icialosc((etot-emin)/delte)
 435               nhist(inde)=nhist(inde)+1.0D0
 436 #ifdef MPL
 437               if ( (it/message_frequency)*message_frequency.eq.it
 438      &                              .and. (MyID.ne.MasterID) ) then
 439                 call recv_stop_sig(Kwita)
 440                 call send_MCM_info(message_frequency)
 441               endif
 442 #endif
 443 C Store the accepted conf. and its energy.
 444               eold=etot
 445               sold=scur
 446               do i=1,nvar
 447                 varold(i)=varia(i)
 448               enddo
 449 #ifdef MPL
 450               if (Kwita.eq.0) call recv_stop_sig(kwita)
 451 #endif
 452             endif ! accepted
 453           endif ! overlap
 454 #ifdef MPL
 455           if (MyID.eq.MasterID .and.
 456      &        (it/message_frequency)*message_frequency.eq.it) then
 457             call receive_MC_info
 458             if (nacc_tot.ge.maxacc) accepted=.true.
 459           endif
 460 #endif
 461 C         if ((ntrial.gt.maxtrial_iter
 462 C    &       .or. (it/pool_read_freq)*pool_read_freq.eq.it)
 463 C    &       .and. npool.gt.0) then
 464 C Take a conformation from the pool
 465 C           ii=iran_num(1,npool)
 466 C           do i=1,nvar
 467 C             varold(i)=xpool(i,ii)
 468 C           enddo
 469 C           if (ntrial.gt.maxtrial_iter)
 470 C    &      write (iout,*) 'Iteration',it,' max. # of trials exceeded.'
 471 C           write (iout,*)
 472 C    &     'Take conformation',ii,' from the pool energy=',epool(ii)
 473 C           if (print_mc.gt.2)
 474 C    &      write (iout,'(10f8.3)') (rad2deg*varold(i),i=1,nvar)
 475 C           ntrial=0
 476 C           eold=epool(ii)
 477 C           call entropia(eold,sold,indeold)
 478 C           accepted=.true.
 479 C        endif ! (ntrial.gt.maxtrial_iter .and. npool.gt.0)
 480    30    continue
 481         enddo ! accepted
 482 #ifdef MPL
 483         if (MyID.eq.MasterID .and.
 484      &      (it/message_frequency)*message_frequency.eq.it) then
 485           call receive_MC_info
 486         endif
 487         if (Kwita.eq.0) call recv_stop_sig(kwita)
 488 #endif
 489         if (ovrtim()) WhatsUp=-1
 490 cd      write (iout,*) 'WhatsUp=',WhatsUp,' Kwita=',Kwita
 491         not_done = (nacc_tot.lt.maxacc) .and. (WhatsUp.eq.0)
 492      &         .and. (Kwita.eq.0)
 493 cd      write (iout,*) 'not_done=',not_done
 494 #ifdef MPL
 495         if (Kwita.lt.0) then
 496           print *,'Processor',MyID,
 497      &    ' has received STOP signal =',Kwita,' in EntSamp.'
 498 cd        print *,'not_done=',not_done
 499           if (Kwita.lt.-1) WhatsUp=Kwita
 500           if (MyID.ne.MasterID) call send_MCM_info(-1)
 501         else if (nacc_tot.ge.maxacc) then
 502           print *,'Processor',MyID,' calls send_stop_sig,',
 503      &     ' because a sufficient # of confs. have been collected.'
 504 cd        print *,'not_done=',not_done
 505           call send_stop_sig(-1)
 506           if (MyID.ne.MasterID) call send_MCM_info(-1)
 507         else if (WhatsUp.eq.-1) then
 508           print *,'Processor',MyID,
 509      &               ' calls send_stop_sig because of timeout.'
 510 cd        print *,'not_done=',not_done
 511           call send_stop_sig(-2)
 512           if (MyID.ne.MasterID) call send_MCM_info(-1)
 513         endif
 514 #endif
 515       enddo ! not_done
 516
 517 C-----------------------------------------------------------------
 518 C... Construct energy histogram & update entropy
 519 C-----------------------------------------------------------------
 520       go to 21
 521    20 WhatsUp=-3
 522 #ifdef MPL
 523       write (iout,*) 'Processor',MyID,
 524      &       ' is broadcasting ERROR-STOP signal.'
 525       write (*,*) 'Processor',MyID,
 526      &       ' is broadcasting ERROR-STOP signal.'
 527       call send_stop_sig(-3)
 528       if (MyID.ne.MasterID) call send_MCM_info(-1)
 529 #endif
 530    21 continue
 531       write (iout,'(/a)') 'Energy histogram'
 532       do i=-100,100
 533         write (iout,'(i5,2f20.5)') i,emin+i*delte,nhist(i)
 534       enddo
 535 #ifdef MPL
 536 C Wait until every processor has sent complete MC info.
 537       if (MyID.eq.MasterID) then
 538         not_done=.true.
 539         do while (not_done)
 540 C         write (*,*) 'The IFINISH array:'
 541 C         write (*,*) (ifinish(i),i=1,nctasks)
 542           not_done=.false.
 543           do i=2,nctasks
 544             not_done=not_done.or.(ifinish(i).ge.0)
 545           enddo
 546           if (not_done) call receive_MC_info
 547         enddo
 548       endif
 549 C Make collective histogram from the work of all processors.
 550       msglen=(2*max_ene+1)*8
 551       print *,
 552      & 'Processor',MyID,' calls MP_REDUCE to send/receive histograms',
 553      & ' msglen=',msglen
 554       call mp_reduce(nhist,nhist1,msglen,MasterID,d_vadd,
 555      &               cgGroupID)
 556       print *,'Processor',MyID,' MP_REDUCE accomplished for histogr.'
 557       do i=-max_ene,max_ene
 558         nhist(i)=nhist1(i)
 559       enddo
 560 C Collect min. and max. energy
 561       print *,
 562      &'Processor',MyID,' calls MP_REDUCE to send/receive energy borders'
 563       call mp_reduce(elowest,elowest1,8,MasterID,d_vmin,cgGroupID)
 564       call mp_reduce(ehighest,ehighest1,8,MasterID,d_vmax,cgGroupID)
 565       print *,'Processor',MyID,' MP_REDUCE accomplished for energies.'
 566       IF (MyID.eq.MasterID) THEN
 567         elowest=elowest1
 568         ehighest=ehighest1
 569 #endif
 570         write (iout,'(a,i10)') '# of accepted confs:',nacc_tot
 571         write (iout,'(a,f10.5,a,f10.5)') 'Lowest energy:',elowest,
 572      & ' Highest energy',ehighest
 573         indmin=icialosc((elowest-emin)/delte)
 574         imdmax=icialosc((ehighest-emin)/delte)
 575         if (indmin.lt.indminn) then
 576           emax=emin+indmin*delte+e_up
 577           indmaxx=indmin+nbins
 578           indminn=indmin
 579         endif
 580         if (.not.ent_read) ent_read=.true.
 581         write(iout,*)'indminn=',indminn,' indmaxx=',indmaxx
 582 C Update entropy (density of states)
 583         do i=indmin,indmax
 584           if (nhist(i).gt.0) then
 585             entropy(i)=entropy(i)+dlog(nhist(i)+0.0D0)
 586           endif
 587         enddo
 588         write (iout,'(/80(1h*)/a,i2/80(1h*)/)')
 589      &        'End of macroiteration',isweep
 590         write (iout,'(a,f10.5,a,f10.5)') 'Elowest=',elowest,
 591      &      ' Ehighest=',ehighest
 592         write (iout,'(/a)') 'Energy histogram'
 593         do i=indminn,indmaxx
 594           write (iout,'(i5,2f20.5)') i,emin+i*delte,nhist(i)
 595         enddo
 596         write (iout,'(/a)') 'Entropy'
 597         do i=indminn,indmaxx
 598           write (iout,'(i5,2f20.5)') i,emin+i*delte,entropy(i)
 599         enddo
 600 C-----------------------------------------------------------------
 601 C... End of energy histogram construction
 602 C-----------------------------------------------------------------
 603 #ifdef MPL
 604       ELSE
 605         if (.not. ent_read) ent_read=.true.
 606       ENDIF ! MyID .eq. MaterID
 607       if (MyID.eq.MasterID) then
 608         itemp(1)=indminn
 609         itemp(2)=indmaxx
 610       endif
 611       print *,'before mp_bcast processor',MyID,' indminn=',indminn,
 612      & ' indmaxx=',indmaxx,' itemp=',itemp(1),itemp(2)
 613       call mp_bcast(itemp(1),8,MasterID,cgGroupID)
 614       call mp_bcast(emax,8,MasterID,cgGroupID)
 615       print *,'after mp_bcast processor',MyID,' indminn=',indminn,
 616      & ' indmaxx=',indmaxx,' itemp=',itemp(1),itemp(2)
 617       if (MyID .ne. MasterID) then
 618         indminn=itemp(1)
 619         indmaxx=itemp(2)
 620       endif
 621       msglen=(indmaxx-indminn+1)*8
 622       print *,'processor',MyID,' calling mp_bcast msglen=',msglen,
 623      & ' indminn=',indminn,' indmaxx=',indmaxx,' isweep=',isweep
 624       call mp_bcast(entropy(indminn),msglen,MasterID,cgGroupID)
 625       IF(MyID.eq.MasterID .and. .not. ovrtim() .and. WhatsUp.ge.0)THEN
 626         open (ientout,file=entname,status='unknown')
 627         write (ientout,'(2i5,2e25.17)') indminn,indmaxx,emin,emax
 628         do i=indminn,indmaxx
 629           write (ientout,'(i5,f10.5,f20.15)') i,emin+i*delte,entropy(i)
 630         enddo
 631         close(ientout)
 632       ELSE
 633         write (iout,*) 'Received from master:'
 634         write (iout,*) 'indminn=',indminn,' indmaxx=',indmaxx,
 635      &                 ' emin=',emin,' emax=',emax
 636         write (iout,'(/a)') 'Entropy'
 637         do i=indminn,indmaxx
 638            write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
 639         enddo
 640       ENDIF ! MyID.eq.MasterID
 641       print *,'Processor',MyID,' calls MP_GATHER'
 642       call mp_gather(nbond_move,nbond_move1,4*Nbm,MasterID,
 643      &               cgGroupID)
 644       call mp_gather(nbond_acc,nbond_acc1,4*Nbm,MasterID,
 645      &               cgGroupID)
 646       print *,'Processor',MyID,' MP_GATHER call accomplished'
 647       if (MyID.eq.MasterID) then
 648
 649         write (iout,'(/80(1h*)/20x,a)') 'Summary run statistics:'
 650         call statprint(nacc_tot,nfun,iretcode,etot,elowest)
 651         write (iout,'(a)')
 652      &   'Statistics of multiple-bond motions. Total motions:'
 653         write (iout,'(8i10)') (nbond_move(i),i=1,Nbm)
 654         write (iout,'(a)') 'Accepted motions:'
 655         write (iout,'(8i10)') (nbond_acc(i),i=1,Nbm)
 656
 657         write (iout,'(a)')
 658      & 'Statistics of multi-bond moves of respective processors:'
 659         do iproc=1,Nprocs-1
 660           do i=1,Nbm
 661             ind=iproc*nbm+i
 662             nbond_move(i)=nbond_move(i)+nbond_move1(ind)
 663             nbond_acc(i)=nbond_acc(i)+nbond_acc1(ind)
 664           enddo
 665         enddo
 666         do iproc=0,NProcs-1
 667           write (iout,*) 'Processor',iproc,' nbond_move:',
 668      &        (nbond_move1(iproc*nbm+i),i=1,Nbm),
 669      &        ' nbond_acc:',(nbond_acc1(iproc*nbm+i),i=1,Nbm)
 670         enddo
 671       endif
 672       call mp_gather(moves,moves1,4*(MaxMoveType+3),MasterID,
 673      &               cgGroupID)
 674       call mp_gather(moves_acc,moves_acc1,4*(MaxMoveType+3),
 675      &               MasterID,cgGroupID)
 676       if (MyID.eq.MasterID) then
 677         do iproc=1,Nprocs-1
 678           do i=-1,MaxMoveType+1
 679             moves(i)=moves(i)+moves1(i,iproc)
 680             moves_acc(i)=moves_acc(i)+moves_acc1(i,iproc)
 681           enddo
 682         enddo
 683         nmove=0
 684         do i=0,MaxMoveType+1
 685           nmove=nmove+moves(i)
 686         enddo
 687         do iproc=0,NProcs-1
 688           write (iout,*) 'Processor',iproc,' moves',
 689      &     (moves1(i,iproc),i=0,MaxMoveType+1),
 690      &    ' moves_acc:',(moves_acc1(i,iproc),i=0,MaxMoveType+1)
 691         enddo
 692       endif
 693 #else
 694       open (ientout,file=entname,status='unknown')
 695       write (ientout,'(2i5,2e25.17)') indminn,indmaxx,emin,emax
 696       do i=indminn,indmaxx
 697         write (ientout,'(i5,f10.5,f20.15)') i,emin+i*delte,entropy(i)
 698       enddo
 699       close(ientout)
 700 #endif
 701       write (iout,'(/80(1h*)/20x,a)') 'Summary run statistics:'
 702       call statprint(nacc_tot,nfun,iretcode,etot,elowest)
 703       write (iout,'(a)')
 704      & 'Statistics of multiple-bond motions. Total motions:'
 705       write (iout,'(8i10)') (nbond_move(i),i=1,Nbm)
 706       write (iout,'(a)') 'Accepted motions:'
 707       write (iout,'(8i10)') (nbond_acc(i),i=1,Nbm)
 708       if (ovrtim() .or. WhatsUp.lt.0) return
 709
 710 C---------------------------------------------------------------------------
 711       ENDDO ! ISWEEP
 712 C---------------------------------------------------------------------------
 713
 714       runtime=tcpu()
 715
 716       if (isweep.eq.nsweep .and. it.ge.maxacc)
 717      &write (iout,'(/80(1h*)/20x,a/80(1h*)/)') 'All iterations done.'
 718       return
 719       end
 720 c------------------------------------------------------------------------------
 721       subroutine accept_mc(it,ecur,eold,scur,sold,x,xold,accepted)
 722       implicit real*8 (a-h,o-z)
 723       include 'DIMENSIONS'
 724       include 'COMMON.MCM'
 725       include 'COMMON.MCE'
 726       include 'COMMON.IOUNITS'
 727       include 'COMMON.VAR'
 728 #ifdef MPL
 729       include 'COMMON.INFO'
 730 #endif
 731       include 'COMMON.GEO'
 732       double precision ecur,eold,xx,ran_number,bol
 733       double precision x(maxvar),xold(maxvar)
 734       logical accepted
 735 C Check if the conformation is similar.
 736 cd    write (iout,*) 'Enter ACCEPTING'
 737 cd    write (iout,*) 'Old PHI angles:'
 738 cd    write (iout,*) (rad2deg*xold(i),i=1,nphi)
 739 cd    write (iout,*) 'Current angles'
 740 cd    write (iout,*) (rad2deg*x(i),i=1,nphi)
 741 cd    ddif=dif_ang(nphi,x,xold)
 742 cd    write (iout,*) 'Angle norm:',ddif
 743 cd    write (iout,*) 'ecur=',ecur,' emax=',emax
 744       if (ecur.gt.emax) then
 745         accepted=.false.
 746         if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it)
 747      & write (iout,'(a)') 'Conformation rejected as too high in energy'
 748         return
 749       endif
 750 C Else evaluate the entropy of the conf and compare it with that of the previous
 751 C one.
 752       call entropia(ecur,scur,indecur)
 753 cd    print *,'Processor',MyID,' ecur=',ecur,' indecur=',indecur,
 754 cd   & ' scur=',scur,' eold=',eold,' sold=',sold
 755 cd    print *,'deix=',deix,' dent=',dent,' delte=',delte
 756       if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) then
 757         write(iout,*)'it=',it,'ecur=',ecur,' indecur=',indecur,
 758      &   ' scur=',scur
 759         write(iout,*)'eold=',eold,' sold=',sold
 760       endif
 761       if (scur.le.sold) then
 762         accepted=.true.
 763       else
 764 C Else carry out acceptance test
 765         xx=ran_number(0.0D0,1.0D0)
 766         xxh=scur-sold
 767         if (xxh.gt.50.0D0) then
 768           bol=0.0D0
 769         else
 770           bol=exp(-xxh)
 771         endif
 772         if (bol.gt.xx) then
 773           accepted=.true.
 774           if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it)
 775      &       write (iout,'(a)') 'Conformation accepted.'
 776         else
 777           accepted=.false.
 778           if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it)
 779      &       write (iout,'(a)') 'Conformation rejected.'
 780         endif
 781       endif
 782       return
 783       end
 784 c--------------------------------------------------------------------------
 785       integer function icialosc(x)
 786       double precision x
 787       if (x.lt.0.0D0) then
 788         icialosc=dint(x)-1
 789       else
 790         icialosc=dint(x)
 791       endif
 792       return
 793       end
 794 c--------------------------------------------------------------------------
 795       subroutine entropia(ecur,scur,indecur)
 796       implicit real*8 (a-h,o-z)
 797       include 'DIMENSIONS'
 798       include 'COMMON.MCM'
 799       include 'COMMON.MCE'
 800       include 'COMMON.IOUNITS'
 801       double precision ecur,scur
 802       integer indecur
 803       indecur=icialosc((ecur-emin)/delte)
 804       if (iabs(indecur).gt.max_ene) then
 805         if ((it/print_freq)*it.eq.it) write (iout,'(a,2i5)')
 806      &   'Accepting: Index out of range:',indecur
 807         scur=1000.0D0
 808       else if (indecur.ge.indmaxx) then
 809         scur=entropy(indecur)
 810         if (print_mc.gt.0 .and. (it/print_freq)*it.eq.it)
 811      &    write (iout,*)'Energy boundary reached',
 812      &            indmaxx,indecur,entropy(indecur)
 813       else
 814         deix=ecur-(emin+indecur*delte)
 815         dent=entropy(indecur+1)-entropy(indecur)
 816         scur=entropy(indecur)+(dent/delte)*deix
 817       endif
 818       return
 819       end