source/unres/MCM_MD.F90

   1       module mcm_md
   2 !-----------------------------------------------------------------------------
   3       use io_units
   4       use names
   5       use math
   6       use geometry_data, only: nres,nvar,rad2deg
   7       use random, only: iran_num,ran_number
   8       use MD_data
   9       use MCM_data
  10       use geometry
  11       use energy
  12
  13       implicit none
  14 !-----------------------------------------------------------------------------
  15 ! Max. number of move types in MCM
  16 !      integer,parameter :: maxmovetype=4
  17 !-----------------------------------------------------------------------------
  18 ! Max. number of conformations in Master's cache array
  19       integer,parameter :: max_cache=10
  20 !-----------------------------------------------------------------------------
  21 ! Max. number of stored confs. in MC/MCM simulation
  22 !      integer,parameter :: maxsave=20
  23 !-----------------------------------------------------------------------------
  24 ! Number of threads in deformation
  25       integer,parameter :: max_thread=4, max_thread2=2*max_thread
  26 !-----------------------------------------------------------------------------
  27 ! Number of structures to compare at t=0
  28       integer,parameter :: max_threadss=8,max_threadss2=2*max_threadss
  29 !-----------------------------------------------------------------------------
  30 ! Max. number of conformations in the pool
  31       integer,parameter :: max_pool=10
  32 !-----------------------------------------------------------------------------
  33 !-----------------------------------------------------------------------------
  34 ! commom.cache
  35 !      common /cache/
  36       integer :: ncache
  37 !      integer,dimension(max_cache) :: CachSrc  nie używane
  38 !      integer,dimension(max_cache) :: isent,iused
  39 !      logical :: cache_update
  40 !      real(kind=8),dimension(max_cache) :: ecache
  41 !      real(kind=8),dimension(:,:),allocatable :: xcache !(maxvar,max_cache)
  42 !-----------------------------------------------------------------------------
  43 ! common.mce
  44 !      common /mce/
  45 !      real(kind=8) :: emin,emax
  46       real(kind=8),dimension(:),allocatable :: entropy !(-max_ene-4:max_ene)
  47       real(kind=8),dimension(:),allocatable :: nhist !(-max_ene:max_ene)
  48       real(kind=8),dimension(:),allocatable :: nminima !(maxsave)
  49 !      logical :: ent_read
  50       logical :: multican
  51       integer :: indminn,indmaxx
  52 !      common /pool/
  53       integer :: npool
  54 !      real(kind=8) :: pool_fraction
  55       real(kind=8),dimension(:,:),allocatable :: xpool !(maxvar,max_pool)
  56       real(kind=8),dimension(:),allocatable :: epool !(max_pool)
  57 !      common /mce_counters/
  58 !------------------------------------------------------------------------------
  59 !... Following COMMON block contains variables controlling motion.
  60 !------------------------------------------------------------------------------
  61 !      common /move/
  62 !      real(kind=8),dimension(0:MaxMoveType) :: sumpro_type !(0:MaxMoveType)
  63       real(kind=8),dimension(:),allocatable :: sumpro_bond !(0:maxres)
  64       integer :: koniecl,Nbm,MaxSideMove!,nmove
  65       integer,dimension(:),allocatable :: nbond_move,nbond_acc !(maxres)
  66 !      integer,dimension(-1:MaxMoveType+1) :: moves,moves_acc !(-1:MaxMoveType+1)
  67 !      common /accept_stats/
  68 !      integer :: nacc_tot
  69       integer,dimension(:),allocatable :: nacc_part !(0:MaxProcs) !el nie uzywane???
  70 !      common /windows/
  71 !      integer :: nwindow
  72 !      integer,dimension(:),allocatable :: winstart,winend,winlen !(maxres)
  73 !      common /moveID/
  74 !      character(len=16),dimension(-1:MaxMoveType+1) :: MovTypID !(-1:MaxMoveType+1)
  75 !------------------------------------------------------------------------------
  76 !... koniecl - the number of bonds to be considered "end bonds" subjected to
  77 !...          end moves;
  78 !... Nbm - The maximum length of N-bond segment to be moved;
  79 !... MaxSideMove - maximum number of side chains subjected to local moves
  80 !...               simultaneously;
  81 !... nmove - the current number of attempted moves;
  82 !... nbond_move(*) array that stores the total numbers of 2-bond,3-bond,...
  83 !...            moves;
  84 !... nendmove - number of endmoves;
  85 !... nbackmove - number of backbone moves;
  86 !... nsidemove - number of local side chain moves;
  87 !... sumpro_type(*) - array that stores the lower and upper boundary of the
  88 !...                  random-number range that determines the type of move
  89 !...                  (N-bond, backbone or side chain);
  90 !... sumpro_bond(*) - array that stores the probabilities to perform bond
  91 !...                  moves of consecutive segment length.
  92 !... winstart(*) - the starting position of the perturbation window;
  93 !... winend(*) -  the end position of the perturbation window;
  94 !... winlen(*) - length of the perturbation window;
  95 !... nwindow - the number of perturbation windows (0 - entire chain).
  96 !-----------------------------------------------------------------------------
  97 !
  98 !
  99 !-----------------------------------------------------------------------------
 100       contains
 101 !-----------------------------------------------------------------------------
 102 ! compare_s1.F
 103 !-----------------------------------------------------------------------------
 104       subroutine compare_s1(n_thr,num_thread_save,energyx,x,icomp,enetbss,&
 105                              coordss,rms_d,modif,iprint)
 106
 107 ! This subroutine compares the new conformation, whose variables are in X
 108 ! with the previously accumulated conformations whose energies and variables
 109 ! are stored in ENETBSS and COORDSS, respectively. The meaning of other
 110 ! variables is as follows:
 111 !
 112 ! N_THR - on input the previous # of accumulated confs, on output the current
 113 !         # of accumulated confs.
 114 ! N_REPEAT - an array that indicates how many times the structure has already
 115 !         been used to start the reversed-reversing procedure. Addition of
 116 !         a new structure replacement of a structure with a similar, but
 117 !         lower-energy structure resets the respective entry in N_REPEAT to zero
 118 ! I9   -  output unit
 119 ! ENERGYX,X - the energy and variables of the new conformations.
 120 ! ICOMP - comparison result:
 121 !         0 - the new structure is similar to one of the previous ones and does
 122 !             not have a remarkably lower energy and is therefore rejected;
 123 !         1 - the new structure is different and is added to the set, because
 124 !             there is still room in the COORDSS and ENETBSS arrays;
 125 !         2 - the new structure is different, but higher in energy than any
 126 !             previous one and is therefore rejected
 127 !         3 - there is no more room in the COORDSS and ENETBSS arrays, but
 128 !             the new structure is lower in energy than at least the highest-
 129 !             energy previous structure and therefore replaces it.
 130 !         9 - the new structure is similar to a number of previous structures,
 131 !             but has a remarkably lower energy than any of them; therefore
 132 !             replaces all these structures;
 133 ! MODIF - a logical variable that shows whether to include the new structure
 134 !         in the set of accumulated structures
 135
 136 !      implicit real*8 (a-h,o-z)
 137       use geometry_data
 138       use regularize_, only:fitsq
 139 !      include 'DIMENSIONS'
 140 !      include 'COMMON.GEO'
 141 !      include 'COMMON.VAR'
 142 !rc      include 'COMMON.DEFORM'
 143 !      include 'COMMON.IOUNITS'
 144 !el#ifdef UNRES
 145 !el      use geometry_data      !include 'COMMON.CHAIN'
 146 !el#endif
 147
 148       real(kind=8),dimension(6*nres) :: x,x1    !(maxvar) (maxvar=6*maxres)
 149       real(kind=8) :: przes(3),obrot(3,3)
 150       integer :: list(max_thread)
 151       logical :: non_conv,modif
 152       real(kind=8) :: enetbss(max_threadss)
 153       real(kind=8) :: coordss(6*nres,max_threadss)
 154
 155 !!! local variables - el
 156       integer :: n_thr,num_thread_save,icomp,minimize_s_flag,iprint
 157       real(kind=8) :: energyx,energyy,rms_d
 158       integer :: nlist,k,kk,j,i,iresult
 159       real(kind=8) :: enex_jp,roznica
 160
 161       nlist=0
 162 #ifdef UNRES
 163       call var_to_geom(nvar,x)
 164       write(iout,*) 'Warning calling chainbuild'
 165       call chainbuild
 166       do k=1,2*nres
 167        do kk=1,3
 168          cref(kk,k,1)=c(kk,k)
 169        enddo
 170       enddo
 171 #endif
 172 !      write(iout,*)'*ene=',energyx
 173       j=0
 174       enex_jp=-1.0d+99
 175       do i=1,n_thr
 176        do k=1,nvar
 177         x1(k)=coordss(k,i)
 178        enddo
 179        if (iprint.gt.3) then
 180        write (iout,*) 'Compare_ss, i=',i
 181        write (iout,*) 'New structure Energy:',energyx
 182        write (iout,'(10f8.3)') (rad2deg*x(k),k=1,nvar)
 183        write (iout,*) 'Template structure Energy:',enetbss(i)
 184        write (iout,'(10f8.3)') (rad2deg*x1(k),k=1,nvar)
 185        endif
 186
 187 #ifdef UNRES
 188        call var_to_geom(nvar,x1)
 189       write(iout,*) 'Warning calling chainbuild'
 190        call chainbuild
 191 !d     write(iout,*)'C and CREF'
 192 !d     write(iout,'(i5,3f10.5,5x,3f10.5)')(k,(c(j,k),j=1,3),
 193 !d   &           (cref(j,k),j=1,3),k=1,nres)
 194        call fitsq(roznica,c(1,1),cref(1,1,1),nres,przes,obrot,non_conv)
 195        if (non_conv) then
 196          print *,'Problems in FITSQ!!!'
 197          print *,'X'
 198          print '(10f8.3)',(x(k),k=1,nvar)
 199          print *,'X1'
 200          print '(10f8.3)',(x1(k),k=1,nvar)
 201          print *,'C and CREF'
 202          print '(i5,3f10.5,5x,3f10.5)',(k,(c(j,k),j=1,3),&
 203                  (cref(j,k,1),j=1,3),k=1,nres)
 204        endif
 205        roznica=dsqrt(dabs(roznica))
 206        iresult = 1
 207        if (roznica.lt.rms_d) iresult = 0
 208 #else
 209        energyy=enetbss(i)
 210 !el       call cmprs(x,x1,roznica,energyx,energyy,iresult)
 211 #endif
 212        if (iprint.gt.1) write(iout,'(i5,f10.6,$)') i,roznica
 213 !       print '(i5,f8.3)',i,roznica
 214        if(iresult.eq.0) then
 215         nlist = nlist + 1
 216         list(nlist)=i
 217         if (iprint.gt.1) write(iout,*)
 218         if(energyx.ge.enetbss(i)) then
 219          if (iprint.gt.1) &
 220             write(iout,*)'s*>> structure rejected - same as nr ',i, &
 221         ' RMS',roznica
 222          minimize_s_flag=0
 223          icomp=0
 224          go to 1106
 225         endif
 226        endif
 227        if(energyx.lt.enetbss(i).and.enex_jp.lt.enetbss(i))then
 228         j=i
 229         enex_jp=enetbss(i)
 230        endif
 231       enddo
 232       if (iprint.gt.1) write(iout,*)
 233       if(nlist.gt.0) then
 234        if (modif) then
 235          if (iprint.gt.1) &
 236           write(iout,'(a,i3,$)')'s*>> structure accepted1 - repl nr ',&
 237          list(1)
 238        else
 239          if (iprint.gt.1) &
 240           write(iout,'(a,i3)') &
 241           's*>> structure accepted1 - would repl nr ',list(1)
 242        endif
 243        icomp=9
 244        if (.not. modif) goto 1106
 245        j=list(1)
 246        enetbss(j)=energyx
 247        do i=1,nvar
 248         coordss(i,j)=x(i)
 249        enddo
 250        do j=2,nlist
 251         if (iprint.gt.1) write(iout,'(i3,$)')list(j)
 252         do kk=list(j)+1,nlist
 253          enetbss(kk-1)=enetbss(kk)
 254          do i=1,nvar
 255           coordss(i,kk-1)=coordss(i,kk)
 256          enddo
 257        enddo
 258        enddo
 259        if (iprint.gt.1) write(iout,*)
 260        go to 1106
 261       endif
 262       if(n_thr.lt.num_thread_save) then
 263        icomp=1
 264        if (modif) then
 265          if (iprint.gt.1) &
 266           write(iout,*)'s*>> structure accepted - add with nr ',n_thr+1
 267        else
 268          if (iprint.gt.1) &
 269           write(iout,*)'s*>> structure accepted - would add with nr ',&
 270             n_thr+1
 271          goto 1106
 272        endif
 273        n_thr=n_thr+1
 274        enetbss(n_thr)=energyx
 275        do i=1,nvar
 276         coordss(i,n_thr)=x(i)
 277        enddo
 278       else
 279        if(j.eq.0) then
 280         if (iprint.gt.1) &
 281          write(iout,*)'s*>> structure rejected - too high energy'
 282         icomp=2
 283         go to 1106
 284        end if
 285        icomp=3
 286        if (modif) then
 287          if (iprint.gt.1) &
 288            write(iout,*)'s*>> structure accepted - repl nr ',j
 289        else
 290          if (iprint.gt.1) &
 291            write(iout,*)'s*>> structure accepted - would repl nr ',j
 292          goto 1106
 293        endif
 294        enetbss(j)=energyx
 295        do i=1,nvar
 296         coordss(i,j)=x(i)
 297        enddo
 298       end if
 299
 300 1106  continue
 301       return
 302       end subroutine compare_s1
 303 !-----------------------------------------------------------------------------
 304 ! entmcm.F
 305 !-----------------------------------------------------------------------------
 306       subroutine entmcm
 307
 308       use energy_data
 309       use geometry_data
 310       use MPI_data, only:WhatsUp,MyID
 311       use compare_data, only: ener
 312       use control_data, only: minim,refstr
 313       use io_base
 314       use regularize_, only:fitsq
 315       use control, only: tcpu,ovrtim
 316       use compare
 317       use minimm, only:minimize
 318 ! Does modified entropic sampling in the space of minima.
 319 !      implicit real*8 (a-h,o-z)
 320 !      include 'DIMENSIONS'
 321 !      include 'COMMON.IOUNITS'
 322 #ifdef MPL
 323       use MPI_data      !include 'COMMON.INFO'
 324 #endif
 325 !      include 'COMMON.GEO'
 326 !      include 'COMMON.CHAIN'
 327 !      include 'COMMON.MCM'
 328 !      include 'COMMON.MCE'
 329 !      include 'COMMON.CONTACTS'
 330 !      include 'COMMON.CONTROL'
 331 !      include 'COMMON.VAR'
 332 !      include 'COMMON.INTERACT'
 333 !      include 'COMMON.THREAD'
 334 !      include 'COMMON.NAMES'
 335       logical :: accepted,not_done,over,error,lprint    !,ovrtim
 336       integer :: MoveType,nbond
 337 !      integer :: conf_comp
 338       real(kind=8) :: RandOrPert
 339       real(kind=8),dimension(6*nres) :: varia   !(maxvar) (maxvar=6*maxres)
 340       real(kind=8) :: elowest,ehighest,eold
 341       real(kind=8) :: przes(3),obr(3,3)
 342       real(kind=8),dimension(6*nres) :: varold  !(maxvar) (maxvar=6*maxres)
 343       logical :: non_conv
 344       real(kind=8),dimension(0:n_ene) :: energia,energia_ave
 345
 346 !!! local variables -el
 347       integer :: i,ii,kkk,it,j,nacc,ijunk,indmin,indmax,&
 348             ISWEEP,Kwita,iretcode,indeold,iene,noverlap,&
 349             irep,nstart_grow,inde
 350 #ifndef LBFGS
 351       integer :: nfun
 352 #endif
 353       real(kind=8) :: facee,conste,ejunk,etot,rms,co,frac,&
 354        deix,dent,sold,scur,runtime
 355 !
 356
 357 !      if(.not.allocated(varsave)) allocate(varsave(maxvar,maxsave)) !(maxvar,maxsave)
 358 !d    write (iout,*) 'print_mc=',print_mc
 359       WhatsUp=0
 360       maxtrial_iter=50
 361 !---------------------------------------------------------------------------
 362 ! Initialize counters.
 363 !---------------------------------------------------------------------------
 364 ! Total number of generated confs.
 365       ngen=0
 366 ! Total number of moves. In general this won't be equal to the number of
 367 ! attempted moves, because we may want to reject some "bad" confs just by
 368 ! overlap check.
 369       nmove=0
 370 ! Total number of shift (nbond_move(1)), spike, crankshaft, three-bond,...
 371 ! motions.
 372 !el      allocate(nbond_move(nres)) !(maxres)
 373
 374       do i=1,nres
 375         nbond_move(i)=0
 376       enddo
 377 ! Initialize total and accepted number of moves of various kind.
 378       do i=0,MaxMoveType
 379         moves(i)=0
 380         moves_acc(i)=0
 381       enddo
 382 ! Total number of energy evaluations.
 383       neneval=0
 384       nfun=0
 385       indminn=-max_ene
 386       indmaxx=max_ene
 387       delte=0.5D0
 388       facee=1.0D0/(maxacc*delte)
 389       conste=dlog(facee)
 390 ! Read entropy from previous simulations.
 391       if (ent_read) then
 392         read (ientin,*) indminn,indmaxx,emin,emax
 393         print *,'indminn=',indminn,' indmaxx=',indmaxx,' emin=',emin,&
 394                 ' emax=',emax
 395         do i=-max_ene,max_ene
 396           entropy(i)=(emin+i*delte)*betbol
 397         enddo
 398         read (ientin,*) (ijunk,ejunk,entropy(i),i=indminn,indmaxx)
 399         indmin=indminn
 400         indmax=indmaxx
 401         write (iout,*) 'indminn=',indminn,' indmaxx=',indmaxx,&
 402                        ' emin=',emin,' emax=',emax
 403         write (iout,'(/a)') 'Initial entropy'
 404         do i=indminn,indmaxx
 405           write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
 406         enddo
 407       endif ! ent_read
 408 ! Read the pool of conformations
 409       call read_pool
 410 !----------------------------------------------------------------------------
 411 ! Entropy-sampling simulations with continually updated entropy
 412 ! Loop thru simulations
 413 !----------------------------------------------------------------------------
 414       DO ISWEEP=1,NSWEEP
 415 !----------------------------------------------------------------------------
 416 ! Take a conformation from the pool
 417 !----------------------------------------------------------------------------
 418       if (npool.gt.0) then
 419         ii=iran_num(1,npool)
 420         do i=1,nvar
 421           varia(i)=xpool(i,ii)
 422         enddo
 423         write (iout,*) 'Took conformation',ii,' from the pool energy=',&
 424                      epool(ii)
 425         call var_to_geom(nvar,varia)
 426 ! Print internal coordinates of the initial conformation
 427         call intout
 428       else
 429         call gen_rand_conf(1,*20)
 430       endif
 431 !----------------------------------------------------------------------------
 432 ! Compute and print initial energies.
 433 !----------------------------------------------------------------------------
 434       nsave=0
 435 #ifdef MPL
 436       allocate(nsave_part(nctasks))
 437       if (MyID.eq.MasterID) then
 438         do i=1,nctasks
 439           nsave_part(i)=0
 440         enddo
 441       endif
 442 #endif
 443       Kwita=0
 444       WhatsUp=0
 445       write (iout,'(/80(1h*)/a,i2/80(1h*)/)') 'MCE iteration #',isweep
 446       write (iout,'(/80(1h*)/a)') 'Initial energies:'
 447       write(iout,*) 'Warning calling chainbuild'
 448       call chainbuild
 449       call etotal(energia)
 450       etot = energia(0)
 451       call enerprint(energia)
 452 ! Minimize the energy of the first conformation.
 453       if (minim) then
 454         call geom_to_var(nvar,varia)
 455         call minimize(etot,varia,iretcode,nfun)
 456         call etotal(energia)
 457         etot = energia(0)
 458         write (iout,'(/80(1h*)/a/80(1h*))') &
 459           'Results of the first energy minimization:'
 460         call enerprint(energia)
 461       endif
 462       if (refstr) then
 463        kkk=1
 464         call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,kkk),&
 465          nsup,przes,&
 466                    obr,non_conv)
 467         rms=dsqrt(rms)
 468         call contact(.false.,ncont,icont,co)
 469         frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
 470         write (iout,'(a,f8.3,a,f8.3,a,f8.3)') &
 471           'RMS deviation from the reference structure:',rms,&
 472           ' % of native contacts:',frac*100,' contact order:',co
 473         write (istat,'(i5,11(1pe14.5))') 0,&
 474          (energia(print_order(i)),i=1,nprint_ene),etot,rms,frac,co
 475       else
 476         write (istat,'(i5,9(1pe14.5))') 0,&
 477          (energia(print_order(i)),i=1,nprint_ene),etot
 478       endif
 479       close(istat)
 480       neneval=neneval+nfun+1
 481       if (.not. ent_read) then
 482 ! Initialize the entropy array
 483         do i=-max_ene,max_ene
 484          emin=etot
 485 ! Uncomment the line below for actual entropic sampling (start with uniform
 486 ! energy distribution).
 487 !        entropy(i)=0.0D0
 488 ! Uncomment the line below for multicanonical sampling (start with Boltzmann
 489 ! distribution).
 490          entropy(i)=(emin+i*delte)*betbol
 491         enddo
 492         emax=10000000.0D0
 493         emin=etot
 494         write (iout,'(/a)') 'Initial entropy'
 495         do i=indminn,indmaxx
 496           write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
 497         enddo
 498       endif ! ent_read
 499 #ifdef MPL
 500       call recv_stop_sig(Kwita)
 501       if (whatsup.eq.1) then
 502         call send_stop_sig(-2)
 503         not_done=.false.
 504       else if (whatsup.le.-2) then
 505         not_done=.false.
 506       else if (whatsup.eq.2) then
 507         not_done=.false.
 508       else
 509         not_done=.true.
 510       endif
 511 #else
 512       not_done = (iretcode.ne.11)
 513 #endif
 514       write (iout,'(/80(1h*)/20x,a/80(1h*))') &
 515           'Enter Monte Carlo procedure.'
 516       close(igeom)
 517       call briefout(0,etot)
 518       do i=1,nvar
 519         varold(i)=varia(i)
 520       enddo
 521       eold=etot
 522       indeold=(eold-emin)/delte
 523       deix=eold-(emin+indeold*delte)
 524       dent=entropy(indeold+1)-entropy(indeold)
 525 !d    write (iout,*) 'indeold=',indeold,' deix=',deix,' dent=',dent
 526 !d    write (*,*) 'Processor',MyID,' indeold=',indeold,' deix=',deix,
 527 !d   & ' dent=',dent
 528       sold=entropy(indeold)+(dent/delte)*deix
 529       elowest=etot
 530       write (iout,*) 'eold=',eold,' sold=',sold,' elowest=',etot
 531       write (*,*) 'Processor',MyID,' eold=',eold,' sold=',sold,&
 532        ' elowest=',etot
 533       if (minim) call zapis(varia,etot)
 534       nminima(1)=1.0D0
 535 ! NACC is the counter for the accepted conformations of a given processor
 536       nacc=0
 537 ! NACC_TOT counts the total number of accepted conformations
 538       nacc_tot=0
 539 #ifdef MPL
 540       if (MyID.eq.MasterID) then
 541         call receive_MCM_info
 542       else
 543         call send_MCM_info(2)
 544       endif
 545 #endif
 546       do iene=indminn,indmaxx
 547         nhist(iene)=0.0D0
 548       enddo
 549       do i=2,maxsave
 550         nminima(i)=0.0D0
 551       enddo
 552 ! Main loop.
 553 !----------------------------------------------------------------------------
 554       elowest=1.0D10
 555       ehighest=-1.0D10
 556       it=0
 557       do while (not_done)
 558         it=it+1
 559         if (print_mc.gt.0) write (iout,'(80(1h*)/20x,a,i7)') &
 560                                    'Beginning iteration #',it
 561 ! Initialize local counter.
 562         ntrial=0 ! # of generated non-overlapping confs.
 563         noverlap=0 ! # of overlapping confs.
 564         accepted=.false.
 565         do while (.not. accepted .and. WhatsUp.eq.0 .and. Kwita.eq.0)
 566           ntrial=ntrial+1
 567 ! Retrieve the angles of previously accepted conformation
 568           do j=1,nvar
 569             varia(j)=varold(j)
 570           enddo
 571 !d        write (iout,'(a)') 'Old variables:'
 572 !d        write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
 573           call var_to_geom(nvar,varia)
 574 ! Rebuild the chain.
 575           write(iout,*) 'Warning calling chainbuild'
 576           call chainbuild
 577           MoveType=0
 578           nbond=0
 579           lprint=.true.
 580 ! Decide whether to generate a random conformation or perturb the old one
 581           RandOrPert=ran_number(0.0D0,1.0D0)
 582           if (RandOrPert.gt.RanFract) then
 583             if (print_mc.gt.0) &
 584               write (iout,'(a)') 'Perturbation-generated conformation.'
 585             call perturb(error,lprint,MoveType,nbond,1.0D0)
 586             if (error) goto 20
 587             if (MoveType.lt.1 .or. MoveType.gt.MaxMoveType) then
 588               write (iout,'(/a,i7,a/)') 'Error - unknown MoveType=',&
 589                  MoveType,' returned from PERTURB.'
 590               goto 20
 591             endif
 592             write(iout,*) 'Warning calling chainbuild'
 593             call chainbuild
 594           else
 595             MoveType=0
 596             moves(0)=moves(0)+1
 597             nstart_grow=iran_num(3,nres)
 598             if (print_mc.gt.0) &
 599               write (iout,'(2a,i3)') 'Random-generated conformation',&
 600               ' - chain regrown from residue',nstart_grow
 601             call gen_rand_conf(nstart_grow,*30)
 602           endif
 603           call geom_to_var(nvar,varia)
 604 !d        write (iout,'(a)') 'New variables:'
 605 !d        write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
 606           ngen=ngen+1
 607           if (print_mc.gt.0) write (iout,'(a,i5,a,i10,a,i10)') &
 608          'Processor',MyId,' trial move',ntrial,' total generated:',ngen
 609           if (print_mc.gt.0) write (*,'(a,i5,a,i10,a,i10)') &
 610          'Processor',MyId,' trial move',ntrial,' total generated:',ngen
 611           call etotal(energia)
 612           etot = energia(0)
 613 !         call enerprint(energia(0))
 614 !         write (iout,'(2(a,1pe14.5))') 'Etot=',Etot,' Elowest=',Elowest
 615           if (etot-elowest.gt.overlap_cut) then
 616             write (iout,'(a,i5,a,1pe14.5)')  'Iteration',it,&
 617             ' Overlap detected in the current conf.; energy is',etot
 618             neneval=neneval+1
 619             accepted=.false.
 620             noverlap=noverlap+1
 621             if (noverlap.gt.maxoverlap) then
 622               write (iout,'(a)') 'Too many overlapping confs.'
 623               goto 20
 624             endif
 625           else
 626             if (minim) then
 627               call minimize(etot,varia,iretcode,nfun)
 628 !d            write (iout,'(a)') 'Variables after minimization:'
 629 !d            write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
 630               call etotal(energia)
 631               etot = energia(0)
 632               neneval=neneval+nfun+1
 633             endif
 634             if (print_mc.gt.2) then
 635               write (iout,'(a)') 'Total energies of trial conf:'
 636               call enerprint(energia)
 637             else if (print_mc.eq.1) then
 638                write (iout,'(a,i6,a,1pe16.6)') &
 639                'Trial conformation:',ngen,' energy:',etot
 640             endif
 641 !--------------------------------------------------------------------------
 642 !... Acceptance test
 643 !--------------------------------------------------------------------------
 644             accepted=.false.
 645             if (WhatsUp.eq.0) &
 646               call accepting(etot,eold,scur,sold,varia,varold,&
 647                              accepted)
 648             if (accepted) then
 649               nacc=nacc+1
 650               nacc_tot=nacc_tot+1
 651               if (elowest.gt.etot) elowest=etot
 652               if (ehighest.lt.etot) ehighest=etot
 653               moves_acc(MoveType)=moves_acc(MoveType)+1
 654               if (MoveType.eq.1) then
 655                 nbond_acc(nbond)=nbond_acc(nbond)+1
 656               endif
 657 ! Check against conformation repetitions.
 658               irep=conf_comp(varia,etot)
 659 #if defined(AIX) || defined(PGI)
 660               open (istat,file=statname,position='append')
 661 #else
 662               open (istat,file=statname,access='append')
 663 #endif
 664               if (refstr) then
 665                 kkk=1
 666                 call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,kkk),&
 667                            nsup,&
 668                             przes,obr,non_conv)
 669                 rms=dsqrt(rms)
 670                 call contact(.false.,ncont,icont,co)
 671                 frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
 672                 if (print_mc.gt.0) &
 673                 write (iout,'(a,f8.3,a,f8.3,a,f8.3)') &
 674                 'RMS deviation from the reference structure:',rms,&
 675                 ' % of native contacts:',frac*100,' contact order:',co
 676                 if (print_stat) &
 677                 write (istat,'(i5,11(1pe14.5))') it,&
 678                  (energia(print_order(i)),i=1,nprint_ene),etot,&
 679                  rms,frac,co
 680               elseif (print_stat) then
 681                 write (istat,'(i5,10(1pe14.5))') it,&
 682                    (energia(print_order(i)),i=1,nprint_ene),etot
 683               endif
 684               close(istat)
 685               if (print_mc.gt.1) &
 686                 call statprint(nacc,nfun,iretcode,etot,elowest)
 687 ! Print internal coordinates.
 688               if (print_int) call briefout(nacc,etot)
 689 #ifdef MPL
 690               if (MyID.ne.MasterID) then
 691                 call recv_stop_sig(Kwita)
 692 !d              print *,'Processor:',MyID,' STOP=',Kwita
 693                 if (irep.eq.0) then
 694                   call send_MCM_info(2)
 695                 else
 696                   call send_MCM_info(1)
 697                 endif
 698               endif
 699 #endif
 700 ! Store the accepted conf. and its energy.
 701               eold=etot
 702               sold=scur
 703               do i=1,nvar
 704                 varold(i)=varia(i)
 705               enddo
 706               if (irep.eq.0) then
 707                 irep=nsave+1
 708 !d              write (iout,*) 'Accepted conformation:'
 709 !d              write (iout,*) (rad2deg*varia(i),i=1,nphi)
 710                 if (minim) call zapis(varia,etot)
 711                 do i=1,n_ene
 712                   ener(i,nsave)=energia(i)
 713                 enddo
 714                 ener(n_ene+1,nsave)=etot
 715                 ener(n_ene+2,nsave)=frac
 716               endif
 717               nminima(irep)=nminima(irep)+1.0D0
 718 !             print *,'irep=',irep,' nminima=',nminima(irep)
 719 #ifdef MPL
 720               if (Kwita.eq.0) call recv_stop_sig(kwita)
 721 #endif
 722             endif ! accepted
 723           endif ! overlap
 724 #ifdef MPL
 725           if (MyID.eq.MasterID) then
 726             call receive_MCM_info
 727             if (nacc_tot.ge.maxacc) accepted=.true.
 728           endif
 729 #endif
 730           if (ntrial.gt.maxtrial_iter .and. npool.gt.0) then
 731 ! Take a conformation from the pool
 732             ii=iran_num(1,npool)
 733             do i=1,nvar
 734               varia(i)=xpool(i,ii)
 735             enddo
 736             write (iout,*) 'Iteration',it,' max. # of trials exceeded.'
 737             write (iout,*) &
 738            'Take conformation',ii,' from the pool energy=',epool(ii)
 739             if (print_mc.gt.2) &
 740             write (iout,'(10f8.3)') (rad2deg*varia(i),i=1,nvar)
 741             ntrial=0
 742          endif ! (ntrial.gt.maxtrial_iter .and. npool.gt.0)
 743    30    continue
 744         enddo ! accepted
 745 #ifdef MPL
 746         if (MyID.eq.MasterID) then
 747           call receive_MCM_info
 748         endif
 749         if (Kwita.eq.0) call recv_stop_sig(kwita)
 750 #endif
 751         if (ovrtim()) WhatsUp=-1
 752 !d      write (iout,*) 'WhatsUp=',WhatsUp,' Kwita=',Kwita
 753         not_done = (nacc_tot.lt.maxacc) .and. (WhatsUp.eq.0) &
 754                .and. (Kwita.eq.0)
 755 !d      write (iout,*) 'not_done=',not_done
 756 #ifdef MPL
 757         if (Kwita.lt.0) then
 758           print *,'Processor',MyID,&
 759           ' has received STOP signal =',Kwita,' in EntSamp.'
 760 !d        print *,'not_done=',not_done
 761           if (Kwita.lt.-1) WhatsUp=Kwita
 762         else if (nacc_tot.ge.maxacc) then
 763           print *,'Processor',MyID,' calls send_stop_sig,',&
 764            ' because a sufficient # of confs. have been collected.'
 765 !d        print *,'not_done=',not_done
 766           call send_stop_sig(-1)
 767         else if (WhatsUp.eq.-1) then
 768           print *,'Processor',MyID,&
 769                      ' calls send_stop_sig because of timeout.'
 770 !d        print *,'not_done=',not_done
 771           call send_stop_sig(-2)
 772         endif
 773 #endif
 774       enddo ! not_done
 775
 776 !-----------------------------------------------------------------
 777 !... Construct energy histogram & update entropy
 778 !-----------------------------------------------------------------
 779       go to 21
 780    20 WhatsUp=-3
 781 #ifdef MPL
 782       write (iout,*) 'Processor',MyID,&
 783              ' is broadcasting ERROR-STOP signal.'
 784       write (*,*) 'Processor',MyID,&
 785              ' is broadcasting ERROR-STOP signal.'
 786       call send_stop_sig(-3)
 787 #endif
 788    21 continue
 789 #ifdef MPL
 790       if (MyID.eq.MasterID) then
 791 !       call receive_MCM_results
 792         call receive_energies
 793 #endif
 794       do i=1,nsave
 795         if (esave(i).lt.elowest) elowest=esave(i)
 796         if (esave(i).gt.ehighest) ehighest=esave(i)
 797       enddo
 798       write (iout,'(a,i10)') '# of accepted confs:',nacc_tot
 799       write (iout,'(a,f10.5,a,f10.5)') 'Lowest energy:',elowest,&
 800        ' Highest energy',ehighest
 801       if (isweep.eq.1 .and. .not.ent_read) then
 802         emin=elowest
 803         emax=ehighest
 804         write (iout,*) 'EMAX=',emax
 805         indminn=0
 806         indmaxx=(ehighest-emin)/delte
 807         indmin=indminn
 808         indmax=indmaxx
 809         do i=-max_ene,max_ene
 810           entropy(i)=(emin+i*delte)*betbol
 811         enddo
 812         ent_read=.true.
 813       else
 814         indmin=(elowest-emin)/delte
 815         indmax=(ehighest-emin)/delte
 816         if (indmin.lt.indminn) indminn=indmin
 817         if (indmax.gt.indmaxx) indmaxx=indmax
 818       endif
 819       write(iout,*)'indminn=',indminn,' indmaxx=',indmaxx
 820 ! Construct energy histogram
 821       do i=1,nsave
 822         inde=(esave(i)-emin)/delte
 823         nhist(inde)=nhist(inde)+nminima(i)
 824       enddo
 825 ! Update entropy (density of states)
 826       do i=indmin,indmax
 827         if (nhist(i).gt.0) then
 828           entropy(i)=entropy(i)+dlog(nhist(i)+0.0D0)
 829         endif
 830       enddo
 831 !d    do i=indmaxx+1
 832 !d      entropy(i)=1.0D+10
 833 !d    enddo
 834       write (iout,'(/80(1h*)/a,i2/80(1h*)/)') &
 835             'End of macroiteration',isweep
 836       write (iout,'(a,f10.5,a,f10.5)') 'Elowest=',elowest,&
 837             ' Ehighest=',ehighest
 838       write (iout,'(a)') 'Frequecies of minima'
 839       do i=1,nsave
 840         write (iout,'(i5,f5.0,f10.5)') i,nminima(i),esave(i)
 841       enddo
 842       write (iout,'(/a)') 'Energy histogram'
 843       do i=indminn,indmaxx
 844         write (iout,'(i5,2f10.5)') i,emin+i*delte,nhist(i)
 845       enddo
 846       write (iout,'(/a)') 'Entropy'
 847       do i=indminn,indmaxx
 848         write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
 849       enddo
 850 !-----------------------------------------------------------------
 851 !... End of energy histogram construction
 852 !-----------------------------------------------------------------
 853 #ifdef MPL
 854         entropy(-max_ene-4)=dfloat(indminn)
 855         entropy(-max_ene-3)=dfloat(indmaxx)
 856         entropy(-max_ene-2)=emin
 857         entropy(-max_ene-1)=emax
 858         call send_MCM_update
 859 !d      print *,entname,ientout
 860         open (ientout,file=entname,status='unknown')
 861         write (ientout,'(2i5,2e25.17)') indminn,indmaxx,emin,emax
 862         do i=indminn,indmaxx
 863           write (ientout,'(i5,f10.5,f20.15)') i,emin+i*delte,entropy(i)
 864         enddo
 865         close(ientout)
 866       else
 867         write (iout,'(a)') 'Frequecies of minima'
 868         do i=1,nsave
 869           write (iout,'(i5,f5.0,f10.5)') i,nminima(i),esave(i)
 870         enddo
 871 !       call send_MCM_results
 872         call send_energies
 873         call receive_MCM_update
 874         indminn=entropy(-max_ene-4)
 875         indmaxx=entropy(-max_ene-3)
 876         emin=entropy(-max_ene-2)
 877         emax=entropy(-max_ene-1)
 878         write (iout,*) 'Received from master:'
 879         write (iout,*) 'indminn=',indminn,' indmaxx=',indmaxx,&
 880                        ' emin=',emin,' emax=',emax
 881         write (iout,'(/a)') 'Entropy'
 882         do i=indminn,indmaxx
 883           write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
 884         enddo
 885       endif
 886       if (WhatsUp.lt.-1) return
 887 #else
 888       if (ovrtim() .or. WhatsUp.lt.0) return
 889 #endif
 890
 891       write (iout,'(/80(1h*)/20x,a)') 'Summary run statistics:'
 892       call statprint(nacc,nfun,iretcode,etot,elowest)
 893       write (iout,'(a)') &
 894        'Statistics of multiple-bond motions. Total motions:'
 895       write (iout,'(16i5)') (nbond_move(i),i=1,Nbm)
 896       write (iout,'(a)') 'Accepted motions:'
 897       write (iout,'(16i5)') (nbond_acc(i),i=1,Nbm)
 898 !el      write (iout,'(a,i10)') 'Number of chain regrowths:',nregrow
 899 !el      write (iout,'(a,i10)') 'Accepted chain regrowths:',nregrow_acc
 900
 901 !---------------------------------------------------------------------------
 902       ENDDO ! ISWEEP
 903 !---------------------------------------------------------------------------
 904
 905       runtime=tcpu()
 906
 907       if (isweep.eq.nsweep .and. it.ge.maxacc) &
 908        write (iout,'(/80(1h*)/20x,a/80(1h*)/)') 'All iterations done.'
 909       return
 910       end subroutine entmcm
 911 !-----------------------------------------------------------------------------
 912       subroutine accepting(ecur,eold,scur,sold,x,xold,accepted)
 913
 914       use geometry_data, only: nphi
 915       use energy_data, only: max_ene
 916 !      implicit real*8 (a-h,o-z)
 917 !      include 'DIMENSIONS'
 918 !      include 'COMMON.MCM'
 919 !      include 'COMMON.MCE'
 920 !      include 'COMMON.IOUNITS'
 921 !      include 'COMMON.VAR'
 922 #ifdef MPL
 923       use MPI_data      !include 'COMMON.INFO'
 924 #endif
 925 !      include 'COMMON.GEO'
 926       real(kind=8) :: ecur,eold,xx,bol  !,ran_number
 927       real(kind=8),dimension(6*nres) :: x,xold  !(maxvar) (maxvar=6*maxres)
 928       real(kind=8) :: tole=1.0D-1, tola=5.0D0
 929       logical :: accepted
 930
 931 !!! local variables - el
 932       integer :: indecur
 933       real(kind=8) :: scur,sold,xxh,deix,dent
 934
 935 ! Check if the conformation is similar.
 936 !d    write (iout,*) 'Enter ACCEPTING'
 937 !d    write (iout,*) 'Old PHI angles:'
 938 !d    write (iout,*) (rad2deg*xold(i),i=1,nphi)
 939 !d    write (iout,*) 'Current angles'
 940 !d    write (iout,*) (rad2deg*x(i),i=1,nphi)
 941 !d    ddif=dif_ang(nphi,x,xold)
 942 !d    write (iout,*) 'Angle norm:',ddif
 943 !d    write (iout,*) 'ecur=',ecur,' emax=',emax
 944       if (ecur.gt.emax) then
 945         accepted=.false.
 946         if (print_mc.gt.0) &
 947        write (iout,'(a)') 'Conformation rejected as too high in energy'
 948         return
 949       else if (dabs(ecur-eold).lt.tole .and. &
 950              dif_ang(nphi,x,xold).lt.tola) then
 951         accepted=.false.
 952         if (print_mc.gt.0) &
 953        write (iout,'(a)') 'Conformation rejected as too similar'
 954         return
 955       endif
 956 ! Else evaluate the entropy of the conf and compare it with that of the previous
 957 ! one.
 958       indecur=(ecur-emin)/delte
 959       if (iabs(indecur).gt.max_ene) then
 960         write (iout,'(a,2i5)') &
 961          'Accepting: Index out of range:',indecur
 962         scur=1000.0D0
 963       else if (indecur.eq.indmaxx) then
 964         scur=entropy(indecur)
 965         if (print_mc.gt.0) write (iout,*)'Energy boundary reached',&
 966                   indmaxx,indecur,entropy(indecur)
 967       else
 968         deix=ecur-(emin+indecur*delte)
 969         dent=entropy(indecur+1)-entropy(indecur)
 970         scur=entropy(indecur)+(dent/delte)*deix
 971       endif
 972 !d    print *,'Processor',MyID,' ecur=',ecur,' indecur=',indecur,
 973 !d   & ' scur=',scur,' eold=',eold,' sold=',sold
 974 !d    print *,'deix=',deix,' dent=',dent,' delte=',delte
 975       if (print_mc.gt.1) then
 976         write(iout,*)'ecur=',ecur,' indecur=',indecur,' scur=',scur
 977         write(iout,*)'eold=',eold,' sold=',sold
 978       endif
 979       if (scur.le.sold) then
 980         accepted=.true.
 981       else
 982 ! Else carry out acceptance test
 983         xx=ran_number(0.0D0,1.0D0)
 984         xxh=scur-sold
 985         if (xxh.gt.50.0D0) then
 986           bol=0.0D0
 987         else
 988           bol=exp(-xxh)
 989         endif
 990         if (bol.gt.xx) then
 991           accepted=.true.
 992           if (print_mc.gt.0) write (iout,'(a)') &
 993           'Conformation accepted.'
 994         else
 995           accepted=.false.
 996           if (print_mc.gt.0) write (iout,'(a)') &
 997        'Conformation rejected.'
 998         endif
 999       endif
1000       return
1001       end subroutine accepting
1002 !-----------------------------------------------------------------------------
1003       subroutine read_pool
1004
1005       use io_base, only:read_angles
1006 !      implicit real*8 (a-h,o-z)
1007 !      include 'DIMENSIONS'
1008 !      include 'COMMON.IOUNITS'
1009 !      include 'COMMON.GEO'
1010 !      include 'COMMON.MCM'
1011 !      include 'COMMON.MCE'
1012 !      include 'COMMON.VAR'
1013       real(kind=8),dimension(6*nres) :: varia   !(maxvar) (maxvar=6*maxres)
1014
1015 !!! local variables - el
1016       integer :: j,i,iconf
1017
1018       print '(a)','Call READ_POOL'
1019       do npool=1,max_pool
1020         print *,'i=',i
1021         read (intin,'(i5,f10.5)',end=10,err=10) iconf,epool(npool)
1022         if (epool(npool).eq.0.0D0) goto 10
1023         call read_angles(intin,*10)
1024         call geom_to_var(nvar,xpool(1,npool))
1025       enddo
1026       goto 11
1027    10 npool=npool-1
1028    11 write (iout,'(a,i5)') 'Number of pool conformations:',npool
1029       if (print_mc.gt.2) then
1030       do i=1,npool
1031         write (iout,'(a,i5,a,1pe14.5)') 'Pool conformation',i,' energy',&
1032           epool(i)
1033         write (iout,'(10f8.3)') (rad2deg*xpool(j,i),j=1,nvar)
1034       enddo
1035       endif ! (print_mc.gt.2)
1036       return
1037       end subroutine read_pool
1038 !-----------------------------------------------------------------------------
1039 ! mc.F
1040 !-----------------------------------------------------------------------------
1041       subroutine monte_carlo
1042
1043       use energy_data
1044       use geometry_data
1045       use MPI_data, only:ifinish,nctasks,WhatsUp,MyID,NProcs
1046       use control_data, only:refstr !,MaxProcs
1047       use io_base
1048       use control, only:tcpu,ovrtim
1049       use regularize_, only:fitsq
1050       use compare
1051 !      use control
1052 ! Does Boltzmann and entropic sampling without energy minimization
1053 !      implicit real*8 (a-h,o-z)
1054 !      include 'DIMENSIONS'
1055 !      include 'COMMON.IOUNITS'
1056 #ifdef MPL
1057       use MPI_data      !include 'COMMON.INFO'
1058 #endif
1059 !      include 'COMMON.GEO'
1060 !      include 'COMMON.CHAIN'
1061 !      include 'COMMON.MCM'
1062 !      include 'COMMON.MCE'
1063 !      include 'COMMON.CONTACTS'
1064 !      include 'COMMON.CONTROL'
1065 !      include 'COMMON.VAR'
1066 !      include 'COMMON.INTERACT'
1067 !      include 'COMMON.THREAD'
1068 !      include 'COMMON.NAMES'
1069       logical :: accepted,not_done,over,error,lprint    !,ovrtim
1070       integer :: MoveType,nbond,nbins
1071 !      integer :: conf_comp
1072       real(kind=8) :: RandOrPert
1073       real(kind=8),dimension(6*nres) :: varia   !(maxvar) (maxvar=6*maxres)
1074       real(kind=8) :: elowest,elowest1,ehighest,ehighest1,eold
1075       real(kind=8) :: przes(3),obr(3,3)
1076       real(kind=8),dimension(6*nres) :: varold  !(maxvar) (maxvar=6*maxres)
1077       logical :: non_conv
1078       integer,dimension(-1:MaxMoveType+1,0:NProcs-1) :: moves1,moves_acc1       !(-1:MaxMoveType+1,0:MaxProcs-1)
1079 #ifdef MPL
1080       real(kind=8) :: etot_temp,etot_all(0:NProcs) !(0:MaxProcs)
1081       external d_vadd,d_vmin,d_vmax
1082       real(kind=8),dimension(-max_ene:max_ene) :: entropy1,nhist1
1083       integer,dimension(nres*(NProcs+1)) :: nbond_move1,nbond_acc1 !(nres*(MaxProcs+1))
1084       integer,dimension(2) :: itemp
1085 #endif
1086       real(kind=8),dimension(6*nres) :: var_lowest      !(maxvar) (maxvar=6*maxres)
1087       real(kind=8),dimension(0:n_ene) :: energia,energia_ave
1088 !
1089 !!! local variables - el
1090       integer :: i,j,it,ii,iproc,nacc,ISWEEP,indmax,indmin,ijunk,&
1091             Kwita,indeold,imdmax,inde,iretcode,nstart_grow,noverlap
1092 #ifndef LBFGS
1093       integer :: nfun
1094 #endif
1095       real(kind=8) :: facee,conste,ejunk,etot,sold,frac,runtime,&
1096                  frac_ave,rms_ave,etot_ave,scur,from_pool,co,rms
1097
1098       write(iout,'(a,i8,2x,a,f10.5)') &
1099        'pool_read_freq=',pool_read_freq,' pool_fraction=',pool_fraction
1100       open (istat,file=statname)
1101       WhatsUp=0
1102       indminn=-max_ene
1103       indmaxx=max_ene
1104       facee=1.0D0/(maxacc*delte)
1105 ! Number of bins in energy histogram
1106       nbins=e_up/delte-1
1107       write (iout,*) 'NBINS=',nbins
1108       conste=dlog(facee)
1109 ! Read entropy from previous simulations.
1110       if (ent_read) then
1111         read (ientin,*) indminn,indmaxx,emin,emax
1112         print *,'indminn=',indminn,' indmaxx=',indmaxx,' emin=',emin,&
1113                 ' emax=',emax
1114         do i=-max_ene,max_ene
1115           entropy(i)=0.0D0
1116         enddo
1117         read (ientin,*) (ijunk,ejunk,entropy(i),i=indminn,indmaxx)
1118         indmin=indminn
1119         indmax=indmaxx
1120         write (iout,*) 'indminn=',indminn,' indmaxx=',indmaxx,&
1121                        ' emin=',emin,' emax=',emax
1122         write (iout,'(/a)') 'Initial entropy'
1123         do i=indminn,indmaxx
1124           write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
1125         enddo
1126       endif ! ent_read
1127 ! Read the pool of conformations
1128       call read_pool
1129       elowest=1.0D+10
1130       ehighest=-1.0D+10
1131 !----------------------------------------------------------------------------
1132 ! Entropy-sampling simulations with continually updated entropy;
1133 ! set NSWEEP=1 for Boltzmann sampling.
1134 ! Loop thru simulations
1135 !----------------------------------------------------------------------------
1136       allocate(ifinish(nctasks))
1137       DO ISWEEP=1,NSWEEP
1138 !
1139 ! Initialize the IFINISH array.
1140 !
1141 #ifdef MPL
1142         do i=1,nctasks
1143           ifinish(i)=0
1144         enddo
1145 #endif
1146 !---------------------------------------------------------------------------
1147 ! Initialize counters.
1148 !---------------------------------------------------------------------------
1149 ! Total number of generated confs.
1150         ngen=0
1151 ! Total number of moves. In general this won't be equal to the number of
1152 ! attempted moves, because we may want to reject some "bad" confs just by
1153 ! overlap check.
1154         nmove=0
1155 ! Total number of shift (nbond_move(1)), spike, crankshaft, three-bond,...
1156 ! motions.
1157 !el      allocate(nbond_move(nres)) !(maxres)
1158 !el      allocate(nbond_acc(nres)) !(maxres)
1159
1160         do i=1,nres
1161           nbond_move(i)=0
1162           nbond_acc(i)=0
1163         enddo
1164 ! Initialize total and accepted number of moves of various kind.
1165         do i=-1,MaxMoveType
1166           moves(i)=0
1167           moves_acc(i)=0
1168         enddo
1169 ! Total number of energy evaluations.
1170         neneval=0
1171         nfun=0
1172 !----------------------------------------------------------------------------
1173 ! Take a conformation from the pool
1174 !----------------------------------------------------------------------------
1175       rewind(istat)
1176       write (iout,*) 'emin=',emin,' emax=',emax
1177       if (npool.gt.0) then
1178         ii=iran_num(1,npool)
1179         do i=1,nvar
1180           varia(i)=xpool(i,ii)
1181         enddo
1182         write (iout,*) 'Took conformation',ii,' from the pool energy=',&
1183                      epool(ii)
1184         call var_to_geom(nvar,varia)
1185 ! Print internal coordinates of the initial conformation
1186         call intout
1187       else if (isweep.gt.1) then
1188         if (eold.lt.emax) then
1189         do i=1,nvar
1190           varia(i)=varold(i)
1191         enddo
1192         else
1193         do i=1,nvar
1194           varia(i)=var_lowest(i)
1195         enddo
1196         endif
1197         call var_to_geom(nvar,varia)
1198       endif
1199 !----------------------------------------------------------------------------
1200 ! Compute and print initial energies.
1201 !----------------------------------------------------------------------------
1202       nsave=0
1203       Kwita=0
1204       WhatsUp=0
1205       write (iout,'(/80(1h*)/a,i2/80(1h*)/)') 'MCE iteration #',isweep
1206       write (iout,'(/80(1h*)/a)') 'Initial energies:'
1207       write(iout,*) 'Warning calling chainbuild'
1208       call chainbuild
1209       call geom_to_var(nvar,varia)
1210       call etotal(energia)
1211       etot = energia(0)
1212       call enerprint(energia)
1213       if (refstr) then
1214         call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,1),nsup,przes,&
1215                    obr,non_conv)
1216         rms=dsqrt(rms)
1217         call contact(.false.,ncont,icont,co)
1218         frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
1219         write (iout,'(a,f8.3,a,f8.3,a,f8.3)') &
1220           'RMS deviation from the reference structure:',rms,&
1221           ' % of native contacts:',frac*100,' contact order',co
1222         write (istat,'(i10,16(1pe14.5))') 0,&
1223          (energia(print_order(i)),i=1,nprint_ene),&
1224          etot,rms,frac,co
1225       else
1226         write (istat,'(i10,14(1pe14.5))') 0,&
1227          (energia(print_order(i)),i=1,nprint_ene),etot
1228       endif
1229 !     close(istat)
1230       neneval=neneval+1
1231       if (.not. ent_read) then
1232 ! Initialize the entropy array
1233 #ifdef MPL
1234 ! Collect total energies from other processors.
1235         etot_temp=etot
1236         etot_all(0)=etot
1237         call mp_gather(etot_temp,etot_all,8,MasterID,cgGroupID)
1238         if (MyID.eq.MasterID) then
1239 ! Get the lowest and the highest energy.
1240           print *,'MASTER: etot_temp: ',(etot_all(i),i=0,nprocs-1),&
1241            ' emin=',emin,' emax=',emax
1242           emin=1.0D10
1243           emax=-1.0D10
1244           do i=0,nprocs
1245             if (emin.gt.etot_all(i)) emin=etot_all(i)
1246             if (emax.lt.etot_all(i)) emax=etot_all(i)
1247           enddo
1248           emax=emin+e_up
1249         endif ! MyID.eq.MasterID
1250         etot_all(1)=emin
1251         etot_all(2)=emax
1252         print *,'Processor',MyID,' calls MP_BCAST to send/recv etot_all'
1253         call mp_bcast(etot_all(1),16,MasterID,cgGroupID)
1254         print *,'Processor',MyID,' MP_BCAST to send/recv etot_all ended'
1255         if (MyID.ne.MasterID) then
1256           print *,'Processor:',MyID,etot_all(1),etot_all(2),&
1257                 etot_all(1),etot_all(2)
1258           emin=etot_all(1)
1259           emax=etot_all(2)
1260         endif ! MyID.ne.MasterID
1261         write (iout,*) 'After MP_GATHER etot_temp=',&
1262                        etot_temp,' emin=',emin
1263 #else
1264         emin=etot
1265         emax=emin+e_up
1266         indminn=0
1267         indmin=0
1268 #endif
1269         IF (MULTICAN) THEN
1270 ! Multicanonical sampling - start from Boltzmann distribution
1271           do i=-max_ene,max_ene
1272             entropy(i)=(emin+i*delte)*betbol
1273           enddo
1274         ELSE
1275 ! Entropic sampling - start from uniform distribution of the density of states
1276           do i=-max_ene,max_ene
1277             entropy(i)=0.0D0
1278           enddo
1279         ENDIF ! MULTICAN
1280         write (iout,'(/a)') 'Initial entropy'
1281         do i=indminn,indmaxx
1282           write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
1283         enddo
1284         if (isweep.eq.1) then
1285           emax=emin+e_up
1286           indminn=0
1287           indmin=0
1288           indmaxx=indminn+nbins
1289           indmax=indmaxx
1290         endif ! isweep.eq.1
1291       endif ! .not. ent_read
1292 #ifdef MPL
1293       call recv_stop_sig(Kwita)
1294       if (whatsup.eq.1) then
1295         call send_stop_sig(-2)
1296         not_done=.false.
1297       else if (whatsup.le.-2) then
1298         not_done=.false.
1299       else if (whatsup.eq.2) then
1300         not_done=.false.
1301       else
1302         not_done=.true.
1303       endif
1304 #else
1305       not_done=.true.
1306 #endif
1307       write (iout,'(/80(1h*)/20x,a/80(1h*))') &
1308           'Enter Monte Carlo procedure.'
1309       close(igeom)
1310       call briefout(0,etot)
1311       do i=1,nvar
1312         varold(i)=varia(i)
1313       enddo
1314       eold=etot
1315       call entropia(eold,sold,indeold)
1316 ! NACC is the counter for the accepted conformations of a given processor
1317       nacc=0
1318 ! NACC_TOT counts the total number of accepted conformations
1319       nacc_tot=0
1320 ! Main loop.
1321 !----------------------------------------------------------------------------
1322 ! Zero out average energies
1323       do i=0,n_ene
1324         energia_ave(i)=0.0d0
1325       enddo
1326 ! Initialize energy histogram
1327       do i=-max_ene,max_ene
1328         nhist(i)=0.0D0
1329       enddo   ! i
1330 ! Zero out iteration counter.
1331       it=0
1332       do j=1,nvar
1333         varold(j)=varia(j)
1334        enddo
1335 ! Begin MC iteration loop.
1336       do while (not_done)
1337         it=it+1
1338 ! Initialize local counter.
1339         ntrial=0 ! # of generated non-overlapping confs.
1340         noverlap=0 ! # of overlapping confs.
1341         accepted=.false.
1342         do while (.not. accepted .and. WhatsUp.eq.0 .and. Kwita.eq.0)
1343           ntrial=ntrial+1
1344 ! Retrieve the angles of previously accepted conformation
1345           do j=1,nvar
1346             varia(j)=varold(j)
1347           enddo
1348           call var_to_geom(nvar,varia)
1349 ! Rebuild the chain.
1350           write(iout,*) 'Warning calling chainbuild'
1351           call chainbuild
1352           MoveType=0
1353           nbond=0
1354           lprint=.true.
1355 ! Decide whether to take a conformation from the pool or generate/perturb one
1356 ! randomly
1357           from_pool=ran_number(0.0D0,1.0D0)
1358           if (npool.gt.0 .and. from_pool.lt.pool_fraction) then
1359 ! Throw a dice to choose the conformation from the pool
1360             ii=iran_num(1,npool)
1361             do i=1,nvar
1362               varia(i)=xpool(i,ii)
1363             enddo
1364             call var_to_geom(nvar,varia)
1365             write(iout,*) 'Warning calling chainbuild'
1366             call chainbuild
1367 !d          call intout
1368 !d          write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
1369             if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) &
1370               write (iout,'(a,i3,a,f10.5)') &
1371            'Try conformation',ii,' from the pool energy=',epool(ii)
1372             MoveType=-1
1373             moves(-1)=moves(-1)+1
1374           else
1375 ! Decide whether to generate a random conformation or perturb the old one
1376           RandOrPert=ran_number(0.0D0,1.0D0)
1377           if (RandOrPert.gt.RanFract) then
1378             if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) &
1379               write (iout,'(a)') 'Perturbation-generated conformation.'
1380             call perturb(error,lprint,MoveType,nbond,0.1D0)
1381             if (error) goto 20
1382             if (MoveType.lt.1 .or. MoveType.gt.MaxMoveType) then
1383               write (iout,'(/a,i7,a/)') 'Error - unknown MoveType=',&
1384                  MoveType,' returned from PERTURB.'
1385               goto 20
1386             endif
1387             call chainbuild
1388           else
1389             MoveType=0
1390             moves(0)=moves(0)+1
1391             nstart_grow=iran_num(3,nres)
1392             if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) &
1393               write (iout,'(2a,i3)') 'Random-generated conformation',&
1394               ' - chain regrown from residue',nstart_grow
1395             call gen_rand_conf(nstart_grow,*30)
1396           endif
1397           call geom_to_var(nvar,varia)
1398           endif ! pool
1399 !d        write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
1400           ngen=ngen+1
1401           if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) &
1402             write (iout,'(a,i5,a,i10,a,i10)') &
1403          'Processor',MyId,' trial move',ntrial,' total generated:',ngen
1404           if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) &
1405             write (*,'(a,i5,a,i10,a,i10)') &
1406          'Processor',MyId,' trial move',ntrial,' total generated:',ngen
1407           call etotal(energia)
1408           etot = energia(0)
1409           neneval=neneval+1
1410 !d        call enerprint(energia(0))
1411 !d        write(iout,*)'it=',it,' etot=',etot
1412           if (etot-elowest.gt.overlap_cut) then
1413             if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) &
1414               write (iout,'(a,i5,a,1pe14.5)')  'Iteration',it,&
1415              ' Overlap detected in the current conf.; energy is',etot
1416             accepted=.false.
1417             noverlap=noverlap+1
1418             if (noverlap.gt.maxoverlap) then
1419               write (iout,'(a)') 'Too many overlapping confs.'
1420               goto 20
1421             endif
1422           else
1423 !--------------------------------------------------------------------------
1424 !... Acceptance test
1425 !--------------------------------------------------------------------------
1426             accepted=.false.
1427             if (WhatsUp.eq.0) &
1428             call accept_mc(it,etot,eold,scur,sold,varia,varold,accepted)
1429             if (accepted) then
1430               nacc=nacc+1
1431               nacc_tot=nacc_tot+1
1432               if (elowest.gt.etot) then
1433                 elowest=etot
1434                 do i=1,nvar
1435                   var_lowest(i)=varia(i)
1436                 enddo
1437               endif
1438               if (ehighest.lt.etot) ehighest=etot
1439               moves_acc(MoveType)=moves_acc(MoveType)+1
1440               if (MoveType.eq.1) then
1441                 nbond_acc(nbond)=nbond_acc(nbond)+1
1442               endif
1443 ! Compare with reference structure.
1444               if (refstr) then
1445                 call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,1),&
1446                             nsup,przes,obr,non_conv)
1447                 rms=dsqrt(rms)
1448                 call contact(.false.,ncont,icont,co)
1449                 frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
1450               endif ! refstr
1451 !
1452 ! Periodically save average energies and confs.
1453 !
1454               do i=0,n_ene
1455                 energia_ave(i)=energia_ave(i)+energia(i)
1456               enddo
1457               moves(MaxMoveType+1)=nmove
1458               moves_acc(MaxMoveType+1)=nacc
1459               IF ((it/save_frequency)*save_frequency.eq.it) THEN
1460                 do i=0,n_ene
1461                   energia_ave(i)=energia_ave(i)/save_frequency
1462                 enddo
1463                 etot_ave=energia_ave(0)
1464 !#ifdef AIX
1465 !                open (istat,file=statname,position='append')
1466 !#else
1467 !                open (istat,file=statname,access='append')
1468 !endif
1469                 if (print_mc.gt.0) &
1470                   write (iout,'(80(1h*)/20x,a,i20)') &
1471                                    'Iteration #',it
1472                 if (refstr .and. print_mc.gt.0)  then
1473                   write (iout,'(a,f8.3,a,f8.3,a,f8.3)') &
1474                   'RMS deviation from the reference structure:',rms,&
1475                   ' % of native contacts:',frac*100,' contact order:',co
1476                 endif
1477                 if (print_stat) then
1478                   if (refstr) then
1479                     write (istat,'(i10,10(1pe14.5))') it,&
1480                     (energia_ave(print_order(i)),i=1,nprint_ene),&
1481                       etot_ave,rms_ave,frac_ave
1482                   else
1483                     write (istat,'(i10,10(1pe14.5))') it,&
1484                     (energia_ave(print_order(i)),i=1,nprint_ene),&
1485                      etot_ave
1486                   endif
1487                 endif
1488 !               close(istat)
1489                 if (print_mc.gt.0) &
1490                   call statprint(nacc,nfun,iretcode,etot,elowest)
1491 ! Print internal coordinates.
1492                 if (print_int) call briefout(nacc,etot)
1493                 do i=0,n_ene
1494                   energia_ave(i)=0.0d0
1495                 enddo
1496               ENDIF ! ( (it/save_frequency)*save_frequency.eq.it)
1497 ! Update histogram
1498               inde=icialosc((etot-emin)/delte)
1499               nhist(inde)=nhist(inde)+1.0D0
1500 #ifdef MPL
1501               if ( (it/message_frequency)*message_frequency.eq.it &
1502                                     .and. (MyID.ne.MasterID) ) then
1503                 call recv_stop_sig(Kwita)
1504                 call send_MCM_info(message_frequency)
1505               endif
1506 #endif
1507 ! Store the accepted conf. and its energy.
1508               eold=etot
1509               sold=scur
1510               do i=1,nvar
1511                 varold(i)=varia(i)
1512               enddo
1513 #ifdef MPL
1514               if (Kwita.eq.0) call recv_stop_sig(kwita)
1515 #endif
1516             endif ! accepted
1517           endif ! overlap
1518 #ifdef MPL
1519           if (MyID.eq.MasterID .and. &
1520               (it/message_frequency)*message_frequency.eq.it) then
1521             call receive_MC_info
1522             if (nacc_tot.ge.maxacc) accepted=.true.
1523           endif
1524 #endif
1525 !         if ((ntrial.gt.maxtrial_iter
1526 !    &       .or. (it/pool_read_freq)*pool_read_freq.eq.it)
1527 !    &       .and. npool.gt.0) then
1528 ! Take a conformation from the pool
1529 !           ii=iran_num(1,npool)
1530 !           do i=1,nvar
1531 !             varold(i)=xpool(i,ii)
1532 !           enddo
1533 !           if (ntrial.gt.maxtrial_iter)
1534 !    &      write (iout,*) 'Iteration',it,' max. # of trials exceeded.'
1535 !           write (iout,*)
1536 !    &     'Take conformation',ii,' from the pool energy=',epool(ii)
1537 !           if (print_mc.gt.2)
1538 !    &      write (iout,'(10f8.3)') (rad2deg*varold(i),i=1,nvar)
1539 !           ntrial=0
1540 !           eold=epool(ii)
1541 !           call entropia(eold,sold,indeold)
1542 !           accepted=.true.
1543 !        endif ! (ntrial.gt.maxtrial_iter .and. npool.gt.0)
1544    30    continue
1545         enddo ! accepted
1546 #ifdef MPL
1547         if (MyID.eq.MasterID .and. &
1548             (it/message_frequency)*message_frequency.eq.it) then
1549           call receive_MC_info
1550         endif
1551         if (Kwita.eq.0) call recv_stop_sig(kwita)
1552 #endif
1553         if (ovrtim()) WhatsUp=-1
1554 !d      write (iout,*) 'WhatsUp=',WhatsUp,' Kwita=',Kwita
1555         not_done = (nacc_tot.lt.maxacc) .and. (WhatsUp.eq.0) &
1556                .and. (Kwita.eq.0)
1557 !d      write (iout,*) 'not_done=',not_done
1558 #ifdef MPL
1559         if (Kwita.lt.0) then
1560           print *,'Processor',MyID,&
1561           ' has received STOP signal =',Kwita,' in EntSamp.'
1562 !d        print *,'not_done=',not_done
1563           if (Kwita.lt.-1) WhatsUp=Kwita
1564           if (MyID.ne.MasterID) call send_MCM_info(-1)
1565         else if (nacc_tot.ge.maxacc) then
1566           print *,'Processor',MyID,' calls send_stop_sig,',&
1567            ' because a sufficient # of confs. have been collected.'
1568 !d        print *,'not_done=',not_done
1569           call send_stop_sig(-1)
1570           if (MyID.ne.MasterID) call send_MCM_info(-1)
1571         else if (WhatsUp.eq.-1) then
1572           print *,'Processor',MyID,&
1573                      ' calls send_stop_sig because of timeout.'
1574 !d        print *,'not_done=',not_done
1575           call send_stop_sig(-2)
1576           if (MyID.ne.MasterID) call send_MCM_info(-1)
1577         endif
1578 #endif
1579       enddo ! not_done
1580
1581 !-----------------------------------------------------------------
1582 !... Construct energy histogram & update entropy
1583 !-----------------------------------------------------------------
1584       go to 21
1585    20 WhatsUp=-3
1586 #ifdef MPL
1587       write (iout,*) 'Processor',MyID,&
1588              ' is broadcasting ERROR-STOP signal.'
1589       write (*,*) 'Processor',MyID,&
1590              ' is broadcasting ERROR-STOP signal.'
1591       call send_stop_sig(-3)
1592       if (MyID.ne.MasterID) call send_MCM_info(-1)
1593 #endif
1594    21 continue
1595       write (iout,'(/a)') 'Energy histogram'
1596       do i=-100,100
1597         write (iout,'(i5,2f20.5)') i,emin+i*delte,nhist(i)
1598       enddo
1599 #ifdef MPL
1600 ! Wait until every processor has sent complete MC info.
1601       if (MyID.eq.MasterID) then
1602         not_done=.true.
1603         do while (not_done)
1604 !         write (*,*) 'The IFINISH array:'
1605 !         write (*,*) (ifinish(i),i=1,nctasks)
1606           not_done=.false.
1607           do i=2,nctasks
1608             not_done=not_done.or.(ifinish(i).ge.0)
1609           enddo
1610           if (not_done) call receive_MC_info
1611         enddo
1612       endif
1613 ! Make collective histogram from the work of all processors.
1614       msglen=(2*max_ene+1)*8
1615       print *,&
1616        'Processor',MyID,' calls MP_REDUCE to send/receive histograms',&
1617        ' msglen=',msglen
1618       call mp_reduce(nhist,nhist1,msglen,MasterID,d_vadd,&
1619                      cgGroupID)
1620       print *,'Processor',MyID,' MP_REDUCE accomplished for histogr.'
1621       do i=-max_ene,max_ene
1622         nhist(i)=nhist1(i)
1623       enddo
1624 ! Collect min. and max. energy
1625       print *, &
1626       'Processor',MyID,' calls MP_REDUCE to send/receive energy borders'
1627       call mp_reduce(elowest,elowest1,8,MasterID,d_vmin,cgGroupID)
1628       call mp_reduce(ehighest,ehighest1,8,MasterID,d_vmax,cgGroupID)
1629       print *,'Processor',MyID,' MP_REDUCE accomplished for energies.'
1630       IF (MyID.eq.MasterID) THEN
1631         elowest=elowest1
1632         ehighest=ehighest1
1633 #endif
1634         write (iout,'(a,i10)') '# of accepted confs:',nacc_tot
1635         write (iout,'(a,f10.5,a,f10.5)') 'Lowest energy:',elowest,&
1636        ' Highest energy',ehighest
1637         indmin=icialosc((elowest-emin)/delte)
1638         imdmax=icialosc((ehighest-emin)/delte)
1639         if (indmin.lt.indminn) then
1640           emax=emin+indmin*delte+e_up
1641           indmaxx=indmin+nbins
1642           indminn=indmin
1643         endif
1644         if (.not.ent_read) ent_read=.true.
1645         write(iout,*)'indminn=',indminn,' indmaxx=',indmaxx
1646 ! Update entropy (density of states)
1647         do i=indmin,indmax
1648           if (nhist(i).gt.0) then
1649             entropy(i)=entropy(i)+dlog(nhist(i)+0.0D0)
1650           endif
1651         enddo
1652         write (iout,'(/80(1h*)/a,i2/80(1h*)/)') &
1653               'End of macroiteration',isweep
1654         write (iout,'(a,f10.5,a,f10.5)') 'Elowest=',elowest,&
1655             ' Ehighest=',ehighest
1656         write (iout,'(/a)') 'Energy histogram'
1657         do i=indminn,indmaxx
1658           write (iout,'(i5,2f20.5)') i,emin+i*delte,nhist(i)
1659         enddo
1660         write (iout,'(/a)') 'Entropy'
1661         do i=indminn,indmaxx
1662           write (iout,'(i5,2f20.5)') i,emin+i*delte,entropy(i)
1663         enddo
1664 !-----------------------------------------------------------------
1665 !... End of energy histogram construction
1666 !-----------------------------------------------------------------
1667 #ifdef MPL
1668       ELSE
1669         if (.not. ent_read) ent_read=.true.
1670       ENDIF ! MyID .eq. MaterID
1671       if (MyID.eq.MasterID) then
1672         itemp(1)=indminn
1673         itemp(2)=indmaxx
1674       endif
1675       print *,'before mp_bcast processor',MyID,' indminn=',indminn,&
1676        ' indmaxx=',indmaxx,' itemp=',itemp(1),itemp(2)
1677       call mp_bcast(itemp(1),8,MasterID,cgGroupID)
1678       call mp_bcast(emax,8,MasterID,cgGroupID)
1679       print *,'after mp_bcast processor',MyID,' indminn=',indminn,&
1680        ' indmaxx=',indmaxx,' itemp=',itemp(1),itemp(2)
1681       if (MyID .ne. MasterID) then
1682         indminn=itemp(1)
1683         indmaxx=itemp(2)
1684       endif
1685       msglen=(indmaxx-indminn+1)*8
1686       print *,'processor',MyID,' calling mp_bcast msglen=',msglen,&
1687        ' indminn=',indminn,' indmaxx=',indmaxx,' isweep=',isweep
1688       call mp_bcast(entropy(indminn),msglen,MasterID,cgGroupID)
1689       IF(MyID.eq.MasterID .and. .not. ovrtim() .and. WhatsUp.ge.0)THEN
1690         open (ientout,file=entname,status='unknown')
1691         write (ientout,'(2i5,2e25.17)') indminn,indmaxx,emin,emax
1692         do i=indminn,indmaxx
1693           write (ientout,'(i5,f10.5,f20.15)') i,emin+i*delte,entropy(i)
1694         enddo
1695         close(ientout)
1696       ELSE
1697         write (iout,*) 'Received from master:'
1698         write (iout,*) 'indminn=',indminn,' indmaxx=',indmaxx,&
1699                        ' emin=',emin,' emax=',emax
1700         write (iout,'(/a)') 'Entropy'
1701         do i=indminn,indmaxx
1702            write (iout,'(i5,2f10.5)') i,emin+i*delte,entropy(i)
1703         enddo
1704       ENDIF ! MyID.eq.MasterID
1705       print *,'Processor',MyID,' calls MP_GATHER'
1706       call mp_gather(nbond_move,nbond_move1,4*Nbm,MasterID,&
1707                      cgGroupID)
1708       call mp_gather(nbond_acc,nbond_acc1,4*Nbm,MasterID,&
1709                      cgGroupID)
1710       print *,'Processor',MyID,' MP_GATHER call accomplished'
1711       if (MyID.eq.MasterID) then
1712
1713         write (iout,'(/80(1h*)/20x,a)') 'Summary run statistics:'
1714         call statprint(nacc_tot,nfun,iretcode,etot,elowest)
1715         write (iout,'(a)') &
1716          'Statistics of multiple-bond motions. Total motions:'
1717         write (iout,'(8i10)') (nbond_move(i),i=1,Nbm)
1718         write (iout,'(a)') 'Accepted motions:'
1719         write (iout,'(8i10)') (nbond_acc(i),i=1,Nbm)
1720
1721         write (iout,'(a)') &
1722        'Statistics of multi-bond moves of respective processors:'
1723         do iproc=1,Nprocs-1
1724           do i=1,Nbm
1725             ind=iproc*nbm+i
1726             nbond_move(i)=nbond_move(i)+nbond_move1(ind)
1727             nbond_acc(i)=nbond_acc(i)+nbond_acc1(ind)
1728           enddo
1729         enddo
1730         do iproc=0,NProcs-1
1731           write (iout,*) 'Processor',iproc,' nbond_move:', &
1732               (nbond_move1(iproc*nbm+i),i=1,Nbm),&
1733               ' nbond_acc:',(nbond_acc1(iproc*nbm+i),i=1,Nbm)
1734         enddo
1735       endif
1736       call mp_gather(moves,moves1,4*(MaxMoveType+3),MasterID,&
1737                      cgGroupID)
1738       call mp_gather(moves_acc,moves_acc1,4*(MaxMoveType+3),&
1739                      MasterID,cgGroupID)
1740       if (MyID.eq.MasterID) then
1741         do iproc=1,Nprocs-1
1742           do i=-1,MaxMoveType+1
1743             moves(i)=moves(i)+moves1(i,iproc)
1744             moves_acc(i)=moves_acc(i)+moves_acc1(i,iproc)
1745           enddo
1746         enddo
1747         nmove=0
1748         do i=0,MaxMoveType+1
1749           nmove=nmove+moves(i)
1750         enddo
1751         do iproc=0,NProcs-1
1752           write (iout,*) 'Processor',iproc,' moves',&
1753            (moves1(i,iproc),i=0,MaxMoveType+1),&
1754           ' moves_acc:',(moves_acc1(i,iproc),i=0,MaxMoveType+1)
1755         enddo
1756       endif
1757 #else
1758       open (ientout,file=entname,status='unknown')
1759       write (ientout,'(2i5,2e25.17)') indminn,indmaxx,emin,emax
1760       do i=indminn,indmaxx
1761         write (ientout,'(i5,f10.5,f20.15)') i,emin+i*delte,entropy(i)
1762       enddo
1763       close(ientout)
1764 #endif
1765       write (iout,'(/80(1h*)/20x,a)') 'Summary run statistics:'
1766       call statprint(nacc_tot,nfun,iretcode,etot,elowest)
1767       write (iout,'(a)') &
1768        'Statistics of multiple-bond motions. Total motions:'
1769       write (iout,'(8i10)') (nbond_move(i),i=1,Nbm)
1770       write (iout,'(a)') 'Accepted motions:'
1771       write (iout,'(8i10)') (nbond_acc(i),i=1,Nbm)
1772       if (ovrtim() .or. WhatsUp.lt.0) return
1773
1774 !---------------------------------------------------------------------------
1775       ENDDO ! ISWEEP
1776 !---------------------------------------------------------------------------
1777
1778       runtime=tcpu()
1779
1780       if (isweep.eq.nsweep .and. it.ge.maxacc) &
1781        write (iout,'(/80(1h*)/20x,a/80(1h*)/)') 'All iterations done.'
1782       return
1783       end subroutine monte_carlo
1784 !-----------------------------------------------------------------------------
1785       subroutine accept_mc(it,ecur,eold,scur,sold,x,xold,accepted)
1786
1787 !      implicit real*8 (a-h,o-z)
1788 !      include 'DIMENSIONS'
1789 !      include 'COMMON.MCM'
1790 !      include 'COMMON.MCE'
1791 !      include 'COMMON.IOUNITS'
1792 !      include 'COMMON.VAR'
1793 #ifdef MPL
1794       use MPI_data      !include 'COMMON.INFO'
1795 #endif
1796 !      include 'COMMON.GEO'
1797       real(kind=8) :: ecur,eold,xx,bol
1798       real(kind=8),dimension(6*nres) :: x,xold  !(maxvar) (maxvar=6*maxres)
1799       logical :: accepted
1800
1801 !el local variables
1802       integer :: it,indecur
1803       real(kind=8) :: scur,sold,xxh
1804 ! Check if the conformation is similar.
1805 !d    write (iout,*) 'Enter ACCEPTING'
1806 !d    write (iout,*) 'Old PHI angles:'
1807 !d    write (iout,*) (rad2deg*xold(i),i=1,nphi)
1808 !d    write (iout,*) 'Current angles'
1809 !d    write (iout,*) (rad2deg*x(i),i=1,nphi)
1810 !d    ddif=dif_ang(nphi,x,xold)
1811 !d    write (iout,*) 'Angle norm:',ddif
1812 !d    write (iout,*) 'ecur=',ecur,' emax=',emax
1813       if (ecur.gt.emax) then
1814         accepted=.false.
1815         if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) &
1816        write (iout,'(a)') 'Conformation rejected as too high in energy'
1817         return
1818       endif
1819 ! Else evaluate the entropy of the conf and compare it with that of the previous
1820 ! one.
1821       call entropia(ecur,scur,indecur)
1822 !d    print *,'Processor',MyID,' ecur=',ecur,' indecur=',indecur,
1823 !d   & ' scur=',scur,' eold=',eold,' sold=',sold
1824 !d    print *,'deix=',deix,' dent=',dent,' delte=',delte
1825       if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) then
1826         write(iout,*)'it=',it,'ecur=',ecur,' indecur=',indecur,&
1827          ' scur=',scur
1828         write(iout,*)'eold=',eold,' sold=',sold
1829       endif
1830       if (scur.le.sold) then
1831         accepted=.true.
1832       else
1833 ! Else carry out acceptance test
1834         xx=ran_number(0.0D0,1.0D0)
1835         xxh=scur-sold
1836         if (xxh.gt.50.0D0) then
1837           bol=0.0D0
1838         else
1839           bol=exp(-xxh)
1840         endif
1841         if (bol.gt.xx) then
1842           accepted=.true.
1843           if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) &
1844              write (iout,'(a)') 'Conformation accepted.'
1845         else
1846           accepted=.false.
1847           if (print_mc.gt.0 .and. (it/print_freq)*print_freq.eq.it) &
1848              write (iout,'(a)') 'Conformation rejected.'
1849         endif
1850       endif
1851       return
1852       end subroutine accept_mc
1853 !-----------------------------------------------------------------------------
1854       integer function icialosc(x)
1855
1856       real(kind=8) :: x
1857       if (x.lt.0.0D0) then
1858         icialosc=dint(x)-1
1859       else
1860         icialosc=dint(x)
1861       endif
1862       return
1863       end function icialosc
1864 !-----------------------------------------------------------------------------
1865       subroutine entropia(ecur,scur,indecur)
1866
1867       use energy_data, only: max_ene
1868 !      implicit real*8 (a-h,o-z)
1869 !      include 'DIMENSIONS'
1870 !      include 'COMMON.MCM'
1871 !      include 'COMMON.MCE'
1872 !      include 'COMMON.IOUNITS'
1873       real(kind=8) :: ecur,scur,deix,dent
1874       integer :: indecur,it   !???el
1875
1876       indecur=icialosc((ecur-emin)/delte)
1877       if (iabs(indecur).gt.max_ene) then
1878         if ((it/print_freq)*it.eq.it) write (iout,'(a,2i5)') &
1879          'Accepting: Index out of range:',indecur
1880         scur=1000.0D0
1881       else if (indecur.ge.indmaxx) then
1882         scur=entropy(indecur)
1883         if (print_mc.gt.0 .and. (it/print_freq)*it.eq.it) &
1884           write (iout,*)'Energy boundary reached',&
1885                   indmaxx,indecur,entropy(indecur)
1886       else
1887         deix=ecur-(emin+indecur*delte)
1888         dent=entropy(indecur+1)-entropy(indecur)
1889         scur=entropy(indecur)+(dent/delte)*deix
1890       endif
1891       return
1892       end subroutine entropia
1893 !-----------------------------------------------------------------------------
1894 ! mcm.F
1895 !-----------------------------------------------------------------------------
1896       subroutine mcm_setup
1897
1898       use energy_data
1899 !      implicit real*8 (a-h,o-z)
1900 !      include 'DIMENSIONS'
1901 !      include 'COMMON.IOUNITS'
1902 !      include 'COMMON.MCM'
1903 !      include 'COMMON.CONTROL'
1904 !      include 'COMMON.INTERACT'
1905 !      include 'COMMON.NAMES'
1906 !      include 'COMMON.CHAIN'
1907 !      include 'COMMON.VAR'
1908 !
1909 !!! local variables - el
1910       integer :: i,i1,i2,it1,it2,ngly,mmm,maxwinlen
1911
1912 ! Set up variables used in MC/MCM.
1913 !
1914 !      allocate(sumpro_bond(0:nres)) !(0:maxres)
1915
1916       write (iout,'(80(1h*)/20x,a/80(1h*))') 'MCM control parameters:'
1917       write (iout,'(5(a,i7))') 'Maxacc:',maxacc,' MaxTrial:',MaxTrial,&
1918        ' MaxRepm:',MaxRepm,' MaxGen:',MaxGen,' MaxOverlap:',MaxOverlap
1919       write (iout,'(4(a,f8.1)/2(a,i3))') &
1920        'Tmin:',Tmin,' Tmax:',Tmax,' TstepH:',TstepH,&
1921        ' TstepC:',TstepC,'NstepH:',NstepH,' NstepC:',NstepC
1922       if (nwindow.gt.0) then
1923         write (iout,'(a)') 'Perturbation windows:'
1924         do i=1,nwindow
1925           i1=winstart(i)
1926           i2=winend(i)
1927           it1=itype(i1,1)
1928           it2=itype(i2,1)
1929           write (iout,'(a,i3,a,i3,a,i3)') restyp(it1,1),i1,restyp(it2,1),i2,&
1930                               ' length',winlen(i)
1931         enddo
1932       endif
1933 ! Rbolt=8.3143D-3*2.388459D-01 kcal/(mol*K)
1934       RBol=1.9858D-3
1935 ! Number of "end bonds".
1936       koniecl=0
1937 !     koniecl=nphi
1938       print *,'koniecl=',koniecl
1939       write (iout,'(a)') 'Probabilities of move types:'
1940       write (*,'(a)') 'Probabilities of move types:'
1941       do i=1,MaxMoveType
1942         write (iout,'(a,f10.5)') MovTypID(i),&
1943           sumpro_type(i)-sumpro_type(i-1)
1944         write (*,'(a,f10.5)') MovTypID(i),&
1945           sumpro_type(i)-sumpro_type(i-1)
1946       enddo
1947       write (iout,*)
1948 ! Maximum length of N-bond segment to be moved
1949 !     nbm=nres-1-(2*koniecl-1)
1950       if (nwindow.gt.0) then
1951         maxwinlen=winlen(1)
1952         do i=2,nwindow
1953           if (winlen(i).gt.maxwinlen) maxwinlen=winlen(i)
1954         enddo
1955         nbm=min0(maxwinlen,6)
1956         write (iout,'(a,i3,a,i3)') 'Nbm=',Nbm,' Maxwinlen=',Maxwinlen
1957       else
1958         nbm=min0(6,nres-2)
1959       endif
1960       sumpro_bond(0)=0.0D0
1961       sumpro_bond(1)=0.0D0
1962       do i=2,nbm
1963         sumpro_bond(i)=sumpro_bond(i-1)+1.0D0/dfloat(i)
1964       enddo
1965       write (iout,'(a)') 'The SumPro_Bond array:'
1966       write (iout,'(8f10.5)') (sumpro_bond(i),i=1,nbm)
1967       write (*,'(a)') 'The SumPro_Bond array:'
1968       write (*,'(8f10.5)') (sumpro_bond(i),i=1,nbm)
1969 ! Maximum number of side chains moved simultaneously
1970 !     print *,'nnt=',nnt,' nct=',nct
1971       ngly=0
1972       do i=nnt,nct
1973         if (itype(i,1).eq.10) ngly=ngly+1
1974       enddo
1975       mmm=nct-nnt-ngly+1
1976       if (mmm.gt.0) then
1977         MaxSideMove=min0((nct-nnt+1)/2,mmm)
1978       endif
1979 !     print *,'MaxSideMove=',MaxSideMove
1980 ! Max. number of generated confs (not used at present).
1981       maxgen=10000
1982 ! Set initial temperature
1983       Tcur=Tmin
1984       betbol=1.0D0/(Rbol*Tcur)
1985       write (iout,'(a,f8.1,a,f10.5)') 'Initial temperature:',Tcur,&
1986           ' BetBol:',betbol
1987       write (iout,*) 'RanFract=',ranfract
1988       return
1989       end subroutine mcm_setup
1990 !-----------------------------------------------------------------------------
1991 #ifndef MPI
1992       subroutine do_mcm(i_orig)
1993
1994       use geometry_data
1995       use energy_data
1996       use MPI_data, only:Whatsup
1997       use control_data, only:refstr,minim,iprint
1998       use io_base
1999       use control, only:tcpu,ovrtim
2000       use regularize_, only:fitsq
2001       use compare
2002       use minimm, only:minimize
2003 ! Monte-Carlo-with-Minimization calculations - serial code.
2004 !      implicit real*8 (a-h,o-z)
2005 !      include 'DIMENSIONS'
2006 !      include 'COMMON.IOUNITS'
2007 !      include 'COMMON.GEO'
2008 !      include 'COMMON.CHAIN'
2009 !      include 'COMMON.MCM'
2010 !      include 'COMMON.CONTACTS'
2011 !      include 'COMMON.CONTROL'
2012 !      include 'COMMON.VAR'
2013 !      include 'COMMON.INTERACT'
2014 !      include 'COMMON.CACHE'
2015 !rc      include 'COMMON.DEFORM'
2016 !rc      include 'COMMON.DEFORM1'
2017 !      include 'COMMON.NAMES'
2018       logical :: accepted,over,error,lprint,not_done,my_conf,&
2019               enelower,non_conv !,ovrtim
2020       integer :: MoveType,nbond !,conf_comp
2021       integer,dimension(max_cache) :: ifeed
2022       real(kind=8),dimension(6*nres) :: varia,varold    !(maxvar) (maxvar=6*maxres)
2023       real(kind=8) :: elowest,eold
2024       real(kind=8) :: przes(3),obr(3,3)
2025       real(kind=8) :: energia(0:n_ene)
2026       real(kind=8) :: coord1(6*nres,max_thread2),enetb1(max_threadss) !el
2027 !!! local variables - el
2028       integer :: i,nf,nacc,it,nout,j,i_orig,nfun,Kwita,iretcode,&
2029             noverlap,nstart_grow,irepet,n_thr,ii
2030       real(kind=8) :: etot,frac,rms,co,RandOrPert,&
2031             rms_deform,runtime
2032 !---------------------------------------------------------------------------
2033 ! Initialize counters.
2034 !---------------------------------------------------------------------------
2035 ! Total number of generated confs.
2036       ngen=0
2037 ! Total number of moves. In general this won't be equal to the number of
2038 ! attempted moves, because we may want to reject some "bad" confs just by
2039 ! overlap check.
2040       nmove=0
2041 ! Total number of temperature jumps.
2042       ntherm=0
2043 ! Total number of shift (nbond_move(1)), spike, crankshaft, three-bond,...
2044 ! motions.
2045 !      if(.not.allocated(varsave)) allocate(varsave(maxvar,maxsave)) !(maxvar,maxsave)
2046 !      allocate(nbond_move(nres)) !(maxres)
2047
2048       ncache=0
2049       do i=1,nres
2050         nbond_move(i)=0
2051       enddo
2052 ! Initialize total and accepted number of moves of various kind.
2053       do i=0,MaxMoveType
2054         moves(i)=0
2055         moves_acc(i)=0
2056       enddo
2057 ! Total number of energy evaluations.
2058       neneval=0
2059       nfun=0
2060       nsave=0
2061
2062       write (iout,*) 'RanFract=',RanFract
2063
2064       WhatsUp=0
2065       Kwita=0
2066
2067 !----------------------------------------------------------------------------
2068 ! Compute and print initial energies.
2069 !----------------------------------------------------------------------------
2070       call intout
2071       write (iout,'(/80(1h*)/a)') 'Initial energies:'
2072       call chainbuild
2073       nf=0
2074
2075       call etotal(energia)
2076       etot = energia(0)
2077 ! Minimize the energy of the first conformation.
2078       if (minim) then
2079         call geom_to_var(nvar,varia)
2080 !       write (iout,*) 'The VARIA array'
2081 !       write (iout,'(8f10.4)') (rad2deg*varia(i),i=1,nvar)
2082         call minimize(etot,varia,iretcode,nfun)
2083         call var_to_geom(nvar,varia)
2084         call chainbuild
2085         write (iout,*) 'etot from MINIMIZE:',etot
2086 !       write (iout,*) 'Tha VARIA array'
2087 !       write (iout,'(8f10.4)') (rad2deg*varia(i),i=1,nvar)
2088
2089         call etotal(energia)
2090         etot=energia(0)
2091         call enerprint(energia)
2092       endif
2093       if (refstr) then
2094         call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,1),nsup,przes,& !el cref(1,nstart_sup)
2095                    obr,non_conv)
2096         rms=dsqrt(rms)
2097         call contact(.false.,ncont,icont,co)
2098         frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
2099         write (iout,'(a,f8.3,a,f8.3,a,f8.3)') &
2100           'RMS deviation from the reference structure:',rms,&
2101           ' % of native contacts:',frac*100,' contact order:',co
2102         if (print_stat) &
2103         write (istat,'(i5,17(1pe14.5))') 0,&
2104          (energia(print_order(i)),i=1,nprint_ene),&
2105          etot,rms,frac,co
2106       else
2107         if (print_stat) write (istat,'(i5,16(1pe14.5))') 0,&
2108          (energia(print_order(i)),i=1,nprint_ene),etot
2109       endif
2110       if (print_stat) close(istat)
2111       neneval=neneval+nfun+1
2112       write (iout,'(/80(1h*)/20x,a/80(1h*))') &
2113           'Enter Monte Carlo procedure.'
2114       if (print_int) then
2115         close(igeom)
2116         call briefout(0,etot)
2117       endif
2118       eold=etot
2119       do i=1,nvar
2120         varold(i)=varia(i)
2121       enddo
2122       elowest=etot
2123       call zapis(varia,etot)
2124       nacc=0         ! total # of accepted confs of the current processor.
2125       nacc_tot=0     ! total # of accepted confs of all processors.
2126
2127       not_done = (iretcode.ne.11)
2128
2129 !----------------------------------------------------------------------------
2130 ! Main loop.
2131 !----------------------------------------------------------------------------
2132       it=0
2133       nout=0
2134       do while (not_done)
2135         it=it+1
2136         write (iout,'(80(1h*)/20x,a,i7)') &
2137                                    'Beginning iteration #',it
2138 ! Initialize local counter.
2139         ntrial=0 ! # of generated non-overlapping confs.
2140         accepted=.false.
2141         do while (.not. accepted)
2142
2143 ! Retrieve the angles of previously accepted conformation
2144           noverlap=0 ! # of overlapping confs.
2145           do j=1,nvar
2146             varia(j)=varold(j)
2147           enddo
2148           call var_to_geom(nvar,varia)
2149 ! Rebuild the chain.
2150           call chainbuild
2151 ! Heat up the system, if necessary.
2152           call heat(over)
2153 ! If temperature cannot be further increased, stop.
2154           if (over) goto 20
2155           MoveType=0
2156           nbond=0
2157           lprint=.true.
2158 !d        write (iout,'(a)') 'Old variables:'
2159 !d        write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
2160 ! Decide whether to generate a random conformation or perturb the old one
2161           RandOrPert=ran_number(0.0D0,1.0D0)
2162           if (RandOrPert.gt.RanFract) then
2163             if (print_mc.gt.0) &
2164               write (iout,'(a)') 'Perturbation-generated conformation.'
2165             call perturb(error,lprint,MoveType,nbond,1.0D0)
2166             if (error) goto 20
2167             if (MoveType.lt.1 .or. MoveType.gt.MaxMoveType) then
2168               write (iout,'(/a,i7,a/)') 'Error - unknown MoveType=',&
2169                  MoveType,' returned from PERTURB.'
2170               goto 20
2171             endif
2172             call chainbuild
2173           else
2174             MoveType=0
2175             moves(0)=moves(0)+1
2176             nstart_grow=iran_num(3,nres)
2177             if (print_mc.gt.0) &
2178               write (iout,'(2a,i3)') 'Random-generated conformation',&
2179               ' - chain regrown from residue',nstart_grow
2180             call gen_rand_conf(nstart_grow,*30)
2181           endif
2182           call geom_to_var(nvar,varia)
2183 !d        write (iout,'(a)') 'New variables:'
2184 !d        write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
2185           ngen=ngen+1
2186
2187           call etotal(energia)
2188           etot=energia(0)
2189 !         call enerprint(energia(0))
2190 !         write (iout,'(2(a,1pe14.5))') 'Etot=',Etot,' Elowest=',Elowest
2191           if (etot-elowest.gt.overlap_cut) then
2192             if(iprint.gt.1.or.etot.lt.1d20) &
2193              write (iout,'(a,1pe14.5)') &
2194             'Overlap detected in the current conf.; energy is',etot
2195             neneval=neneval+1
2196             accepted=.false.
2197             noverlap=noverlap+1
2198             if (noverlap.gt.maxoverlap) then
2199               write (iout,'(a)') 'Too many overlapping confs.'
2200               goto 20
2201             endif
2202           else
2203             if (minim) then
2204               call minimize(etot,varia,iretcode,nfun)
2205 !d            write (iout,*) 'etot from MINIMIZE:',etot
2206 !d            write (iout,'(a)') 'Variables after minimization:'
2207 !d            write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
2208
2209               call etotal(energia)
2210               etot = energia(0)
2211               neneval=neneval+nfun+2
2212             endif
2213 !           call enerprint(energia(0))
2214             write (iout,'(a,i6,a,1pe16.6)') 'Conformation:',ngen,&
2215             ' energy:',etot
2216 !--------------------------------------------------------------------------
2217 !... Do Metropolis test
2218 !--------------------------------------------------------------------------
2219             accepted=.false.
2220             my_conf=.false.
2221
2222             if (WhatsUp.eq.0 .and. Kwita.eq.0) then
2223               call metropolis(nvar,varia,varold,etot,eold,accepted,&
2224                             my_conf,EneLower)
2225             endif
2226             write (iout,*) 'My_Conf=',My_Conf,' EneLower=',EneLower
2227             if (accepted) then
2228
2229               nacc=nacc+1
2230               nacc_tot=nacc_tot+1
2231               if (elowest.gt.etot) elowest=etot
2232               moves_acc(MoveType)=moves_acc(MoveType)+1
2233               if (MoveType.eq.1) then
2234                 nbond_acc(nbond)=nbond_acc(nbond)+1
2235               endif
2236 ! Check against conformation repetitions.
2237               irepet=conf_comp(varia,etot)
2238               if (print_stat) then
2239 #if defined(AIX) || defined(PGI)
2240               open (istat,file=statname,position='append')
2241 #else
2242                open (istat,file=statname,access='append')
2243 #endif
2244               endif
2245               call statprint(nacc,nfun,iretcode,etot,elowest)
2246               if (refstr) then
2247                 call var_to_geom(nvar,varia)
2248                 call chainbuild
2249                 call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,1),&  !el cref(1,nstart_sup)
2250                           nsup,przes,obr,non_conv)
2251                 rms=dsqrt(rms)
2252                 call contact(.false.,ncont,icont,co)
2253                 frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
2254                 write (iout,'(a,f8.3,a,f8.3)') &
2255                 'RMS deviation from the reference structure:',rms,&
2256                 ' % of native contacts:',frac*100,' contact order',co
2257               endif ! refstr
2258               if (My_Conf) then
2259                 nout=nout+1
2260                 write (iout,*) 'Writing new conformation',nout
2261                 if (refstr) then
2262                   write (istat,'(i5,16(1pe14.5))') nout,&
2263                    (energia(print_order(i)),i=1,nprint_ene),&
2264                    etot,rms,frac
2265                 else
2266                   if (print_stat) &
2267                    write (istat,'(i5,17(1pe14.5))') nout,&
2268                     (energia(print_order(i)),i=1,nprint_ene),etot
2269                 endif ! refstr
2270                 if (print_stat) close(istat)
2271 ! Print internal coordinates.
2272                 if (print_int) call briefout(nout,etot)
2273 ! Accumulate the newly accepted conf in the coord1 array, if it is different
2274 ! from all confs that are already there.
2275                 call compare_s1(n_thr,max_thread2,etot,varia,ii,&
2276                  enetb1,coord1,rms_deform,.true.,iprint)
2277                 write (iout,*) 'After compare_ss: n_thr=',n_thr
2278                 if (ii.eq.1 .or. ii.eq.3) then
2279                   write (iout,'(8f10.4)') &
2280                       (rad2deg*coord1(i,n_thr),i=1,nvar)
2281                 endif
2282               else
2283                 write (iout,*) 'Conformation from cache, not written.'
2284               endif ! My_Conf
2285
2286               if (nrepm.gt.maxrepm) then
2287                 write (iout,'(a)') 'Too many conformation repetitions.'
2288                 goto 20
2289               endif
2290 ! Store the accepted conf. and its energy.
2291               eold=etot
2292               do i=1,nvar
2293                 varold(i)=varia(i)
2294               enddo
2295               if (irepet.eq.0) call zapis(varia,etot)
2296 ! Lower the temperature, if necessary.
2297               call cool
2298
2299             else
2300
2301               ntrial=ntrial+1
2302             endif ! accepted
2303           endif ! overlap
2304
2305    30     continue
2306         enddo ! accepted
2307 ! Check for time limit.
2308         if (ovrtim()) WhatsUp=-1
2309         not_done = (nacc_tot.lt.maxacc) .and. (WhatsUp.eq.0) &
2310              .and. (Kwita.eq.0)
2311
2312       enddo ! not_done
2313       goto 21
2314    20 WhatsUp=-3
2315
2316    21 continue
2317       runtime=tcpu()
2318       write (iout,'(/80(1h*)/20x,a)') 'Summary run statistics:'
2319       call statprint(nacc,nfun,iretcode,etot,elowest)
2320       write (iout,'(a)') &
2321        'Statistics of multiple-bond motions. Total motions:'
2322       write (iout,'(16i5)') (nbond_move(i),i=1,Nbm)
2323       write (iout,'(a)') 'Accepted motions:'
2324       write (iout,'(16i5)') (nbond_acc(i),i=1,Nbm)
2325       if (it.ge.maxacc) &
2326       write (iout,'(/80(1h*)/20x,a/80(1h*)/)') 'All iterations done.'
2327         !(maxvar) (maxvar=6*maxres)
2328       return
2329       end subroutine do_mcm
2330 #endif
2331 !-----------------------------------------------------------------------------
2332 #ifdef MPI
2333       subroutine do_mcm(i_orig)
2334
2335 ! Monte-Carlo-with-Minimization calculations - parallel code.
2336       use MPI_data
2337       use control_data, only:refstr!,tag
2338       use io_base, only:intout,briefout
2339       use control, only:ovrtim,tcpu
2340       use compare, only:contact,contact_fract
2341       use minimm, only:minimize
2342       use regularize_, only:fitsq
2343 !      use contact_, only:contact
2344 !      use minim
2345 !      implicit real*8 (a-h,o-z)
2346 !      include 'DIMENSIONS'
2347       include 'mpif.h'
2348 !      include 'COMMON.IOUNITS'
2349 !      include 'COMMON.GEO'
2350 !      include 'COMMON.CHAIN'
2351 !      include 'COMMON.MCM'
2352 !      include 'COMMON.CONTACTS'
2353 !      include 'COMMON.CONTROL'
2354 !      include 'COMMON.VAR'
2355 !      include 'COMMON.INTERACT'
2356 !      include 'COMMON.INFO'
2357 !      include 'COMMON.CACHE'
2358 !rc      include 'COMMON.DEFORM'
2359 !rc      include 'COMMON.DEFORM1'
2360 !rc      include 'COMMON.DEFORM2'
2361 !      include 'COMMON.MINIM'
2362 !      include 'COMMON.NAMES'
2363       logical :: accepted,over,error,lprint,not_done,similar,&
2364               enelower,non_conv,flag,finish     !,ovrtim
2365       integer :: MoveType,nbond !,conf_comp
2366       real(kind=8),dimension(6*nres) :: x1,varold1,varold,varia !(maxvar) (maxvar=6*maxres)
2367       real(kind=8) :: elowest,eold
2368       real(kind=8) :: przes(3),obr(3,3)
2369       integer :: iparentx(max_threadss2)
2370       integer :: iparentx1(max_threadss2)
2371       integer :: imtasks(150),imtasks_n
2372       real(kind=8) :: energia(0:n_ene)
2373
2374 !el local variables
2375       integer :: nfun,nodenum,i_orig,i,nf,nacc,it,nout,j,kkk,is,&
2376             Kwita,iretcode,noverlap,nstart_grow,ierr,iitt,&
2377             ii_grnum_d,ii_ennum_d,ii_hesnum_d,i_grnum_d,i_ennum_d,&
2378             i_hesnum_d,i_minimiz,irepet
2379       real(kind=8) :: etot,frac,eneglobal,RandOrPert,eold1,co,&
2380             runtime,rms
2381
2382 !      if(.not.allocated(varsave)) allocate(varsave(maxvar,maxsave)) !(maxvar,maxsave)
2383       print *,'Master entered DO_MCM'
2384       nodenum = nprocs
2385
2386       finish=.false.
2387       imtasks_n=0
2388       do i=1,nodenum-1
2389        imtasks(i)=0
2390       enddo
2391 !---------------------------------------------------------------------------
2392 ! Initialize counters.
2393 !---------------------------------------------------------------------------
2394 ! Total number of generated confs.
2395       ngen=0
2396 ! Total number of moves. In general this won`t be equal to the number of
2397 ! attempted moves, because we may want to reject some "bad" confs just by
2398 ! overlap check.
2399       nmove=0
2400 ! Total number of temperature jumps.
2401       ntherm=0
2402 ! Total number of shift (nbond_move(1)), spike, crankshaft, three-bond,...
2403 ! motions.
2404       allocate(nbond_move(nres)) !(maxres)
2405
2406       ncache=0
2407       do i=1,nres
2408         nbond_move(i)=0
2409       enddo
2410 ! Initialize total and accepted number of moves of various kind.
2411       do i=0,MaxMoveType
2412         moves(i)=0
2413         moves_acc(i)=0
2414       enddo
2415 ! Total number of energy evaluations.
2416       neneval=0
2417       nfun=0
2418       nsave=0
2419 !      write (iout,*) 'RanFract=',RanFract
2420       WhatsUp=0
2421       Kwita=0
2422 !----------------------------------------------------------------------------
2423 ! Compute and print initial energies.
2424 !----------------------------------------------------------------------------
2425       call intout
2426       write (iout,'(/80(1h*)/a)') 'Initial energies:'
2427       call chainbuild
2428       nf=0
2429       call etotal(energia)
2430       etot = energia(0)
2431       call enerprint(energia)
2432 ! Request energy computation from slave processors.
2433       call geom_to_var(nvar,varia)
2434 !     write (iout,*) 'The VARIA array'
2435 !     write (iout,'(8f10.4)') (rad2deg*varia(i),i=1,nvar)
2436       call minimize(etot,varia,iretcode,nfun)
2437       call var_to_geom(nvar,varia)
2438       call chainbuild
2439       write (iout,*) 'etot from MINIMIZE:',etot
2440 !     write (iout,*) 'Tha VARIA array'
2441 !     write (iout,'(8f10.4)') (rad2deg*varia(i),i=1,nvar)
2442       neneval=0
2443       eneglobal=1.0d99
2444       if (print_mc .gt. 0) write (iout,'(/80(1h*)/20x,a/80(1h*))') &
2445           'Enter Monte Carlo procedure.'
2446       if (print_mc .gt. 0) write (iout,'(i5,1pe14.5)' ) i_orig,etot
2447       eold=etot
2448       do i=1,nvar
2449         varold(i)=varia(i)
2450       enddo
2451       elowest=etot
2452       call zapis(varia,etot)
2453 ! diagnostics
2454       call var_to_geom(nvar,varia)
2455       call chainbuild
2456       call etotal(energia)
2457       if (print_mc.gt.0) write (iout,*) 'Initial energy:',etot
2458 ! end diagnostics
2459       nacc=0         ! total # of accepted confs of the current processor.
2460       nacc_tot=0     ! total # of accepted confs of all processors.
2461       not_done=.true.
2462 !----------------------------------------------------------------------------
2463 ! Main loop.
2464 !----------------------------------------------------------------------------
2465       it=0
2466       nout=0
2467       LOOP1:do while (not_done)
2468         it=it+1
2469         if (print_mc.gt.0) write (iout,'(80(1h*)/20x,a,i7)') &
2470                                    'Beginning iteration #',it
2471 ! Initialize local counter.
2472         ntrial=0 ! # of generated non-overlapping confs.
2473         noverlap=0 ! # of overlapping confs.
2474         accepted=.false.
2475         LOOP2:do while (.not. accepted)
2476
2477          LOOP3:do while (imtasks_n.lt.nodenum-1.and..not.finish)
2478           do i=1,nodenum-1
2479            if(imtasks(i).eq.0) then
2480             is=i
2481             exit
2482            endif
2483           enddo
2484 ! Retrieve the angles of previously accepted conformation
2485           do j=1,nvar
2486             varia(j)=varold(j)
2487           enddo
2488           call var_to_geom(nvar,varia)
2489 ! Rebuild the chain.
2490           call chainbuild
2491 ! Heat up the system, if necessary.
2492           call heat(over)
2493 ! If temperature cannot be further increased, stop.
2494           if (over) then
2495            finish=.true.
2496           endif
2497           MoveType=0
2498           nbond=0
2499 !          write (iout,'(a)') 'Old variables:'
2500 !          write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
2501 ! Decide whether to generate a random conformation or perturb the old one
2502           RandOrPert=ran_number(0.0D0,1.0D0)
2503           if (RandOrPert.gt.RanFract) then
2504            if (print_mc.gt.0) &
2505              write (iout,'(a)') 'Perturbation-generated conformation.'
2506            call perturb(error,lprint,MoveType,nbond,1.0D0)
2507 !           print *,'after perturb',error,finish
2508            if (error) finish = .true.
2509            if (MoveType.lt.1 .or. MoveType.gt.MaxMoveType) then
2510             write (iout,'(/a,i7,a/)') 'Error - unknown MoveType=',&
2511                MoveType,' returned from PERTURB.'
2512             finish=.true.
2513             write (*,'(/a,i7,a/)') 'Error - unknown MoveType=',&
2514                MoveType,' returned from PERTURB.'
2515            endif
2516            call chainbuild
2517           else
2518            MoveType=0
2519            moves(0)=moves(0)+1
2520            nstart_grow=iran_num(3,nres)
2521            if (print_mc.gt.0) &
2522             write (iout,'(2a,i3)') 'Random-generated conformation',&
2523             ' - chain regrown from residue',nstart_grow
2524            call gen_rand_conf(nstart_grow,*30)
2525           endif
2526           call geom_to_var(nvar,varia)
2527           ngen=ngen+1
2528 !          print *,'finish=',finish
2529           if (etot-elowest.gt.overlap_cut) then
2530            if (print_mc.gt.1) write (iout,'(a,1pe14.5)') &
2531           'Overlap detected in the current conf.; energy is',etot
2532            if(iprint.gt.1.or.etot.lt.1d19) print *,&
2533            'Overlap detected in the current conf.; energy is',etot
2534            neneval=neneval+1
2535            accepted=.false.
2536            noverlap=noverlap+1
2537            if (noverlap.gt.maxoverlap) then
2538             write (iout,*) 'Too many overlapping confs.',&
2539             ' etot, elowest, overlap_cut', etot, elowest, overlap_cut
2540             finish=.true.
2541            endif
2542           else if (.not. finish) then
2543 ! Distribute tasks to processors
2544 !           print *,'Master sending order'
2545            call MPI_SEND(12, 1, MPI_INTEGER, is, tag,&
2546                    CG_COMM, ierr)
2547 !           write (iout,*) '12: tag=',tag
2548 !           print *,'Master sent order to processor',is
2549            call MPI_SEND(it, 1, MPI_INTEGER, is, tag,&
2550                    CG_COMM, ierr)
2551 !           write (iout,*) 'it: tag=',tag
2552            call MPI_SEND(eold, 1, MPI_DOUBLE_PRECISION, is, tag,&
2553                    CG_COMM, ierr)
2554 !           write (iout,*) 'eold: tag=',tag
2555            call MPI_SEND(varia(1), nvar, MPI_DOUBLE_PRECISION, &
2556                    is, tag,&
2557                    CG_COMM, ierr)
2558 !           write (iout,*) 'varia: tag=',tag
2559            call MPI_SEND(varold(1), nvar, MPI_DOUBLE_PRECISION, &
2560                    is, tag,&
2561                    CG_COMM, ierr)
2562 !           write (iout,*) 'varold: tag=',tag
2563 #ifdef AIX
2564            call flush_(iout)
2565 #else
2566            call flush(iout)
2567 #endif
2568            imtasks(is)=1
2569            imtasks_n=imtasks_n+1
2570 ! End distribution
2571           endif ! overlap
2572          enddo LOOP3
2573
2574          flag = .false.
2575          LOOP_RECV:do while(.not.flag)
2576           do is=1, nodenum-1
2577            call MPI_IPROBE(is,tag,CG_COMM,flag,status,ierr)
2578            if(flag) then
2579             call MPI_RECV(iitt, 1, MPI_INTEGER, is, tag,&
2580                     CG_COMM, status, ierr)
2581             call MPI_RECV(eold1, 1, MPI_DOUBLE_PRECISION, is, tag,&
2582                     CG_COMM, status, ierr)
2583             call MPI_RECV(etot, 1, MPI_DOUBLE_PRECISION, is, tag,&
2584                     CG_COMM, status, ierr)
2585             call MPI_RECV(varia(1), nvar, MPI_DOUBLE_PRECISION,is,tag,&
2586                     CG_COMM, status, ierr)
2587             call MPI_RECV(varold1(1), nvar, MPI_DOUBLE_PRECISION, is, &
2588                     tag, CG_COMM, status, ierr)
2589             call MPI_RECV(ii_grnum_d, 1, MPI_INTEGER, is, tag,&
2590                     CG_COMM, status, ierr)
2591             call MPI_RECV(ii_ennum_d, 1, MPI_INTEGER, is, tag,&
2592                     CG_COMM, status, ierr)
2593             call MPI_RECV(ii_hesnum_d, 1, MPI_INTEGER, is, tag,&
2594                     CG_COMM, status, ierr)
2595             i_grnum_d=i_grnum_d+ii_grnum_d
2596             i_ennum_d=i_ennum_d+ii_ennum_d
2597             neneval = neneval+ii_ennum_d
2598             i_hesnum_d=i_hesnum_d+ii_hesnum_d
2599             i_minimiz=i_minimiz+1
2600             imtasks(is)=0
2601             imtasks_n=imtasks_n-1
2602             exit
2603            endif
2604           enddo
2605          enddo LOOP_RECV
2606
2607          if(print_mc.gt.0) write (iout,'(a,i6,a,i6,a,i6,a,1pe16.6)') &
2608             'From Worker #',is,' iitt',iitt,&
2609            ' Conformation:',ngen,' energy:',etot
2610 !--------------------------------------------------------------------------
2611 !... Do Metropolis test
2612 !--------------------------------------------------------------------------
2613          call metropolis(nvar,varia,varold1,etot,eold1,accepted,&
2614                             similar,EneLower)
2615          if(iitt.ne.it.and..not.similar) then
2616           call metropolis(nvar,varia,varold,etot,eold,accepted,&
2617                             similar,EneLower)
2618           accepted=enelower
2619          endif
2620          if(etot.lt.eneglobal)eneglobal=etot
2621 !         if(mod(it,100).eq.0)
2622          write(iout,*)'CHUJOJEB ',neneval,eneglobal
2623          if (accepted) then
2624 ! Write the accepted conformation.
2625            nout=nout+1
2626            if (refstr) then
2627              call var_to_geom(nvar,varia)
2628              call chainbuild
2629              kkk=1
2630              call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,kkk),&
2631                           nsup,przes,obr,non_conv)
2632              rms=dsqrt(rms)
2633              call contact(.false.,ncont,icont,co)
2634              frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
2635              write (iout,'(a,f8.3,a,f8.3,a,f8.3)') &
2636                'RMS deviation from the reference structure:',rms,&
2637                ' % of native contacts:',frac*100,' contact order:',co
2638            endif ! refstr
2639            if (print_mc.gt.0) &
2640             write (iout,*) 'Writing new conformation',nout
2641            if (print_stat) then
2642              call var_to_geom(nvar,varia)
2643 #if defined(AIX) || defined(PGI)
2644              open (istat,file=statname,position='append')
2645 #else
2646              open (istat,file=statname,access='append')
2647 #endif
2648              if (refstr) then
2649                write (istat,'(i5,16(1pe14.5))') nout,&
2650                 (energia(print_order(i)),i=1,nprint_ene),&
2651                 etot,rms,frac
2652              else
2653                write (istat,'(i5,16(1pe14.5))') nout,&
2654                 (energia(print_order(i)),i=1,nprint_ene),etot
2655              endif ! refstr
2656              close(istat)
2657            endif ! print_stat
2658 ! Print internal coordinates.
2659            if (print_int) call briefout(nout,etot)
2660            nacc=nacc+1
2661            nacc_tot=nacc_tot+1
2662            if (elowest.gt.etot) elowest=etot
2663            moves_acc(MoveType)=moves_acc(MoveType)+1
2664            if (MoveType.eq.1) then
2665              nbond_acc(nbond)=nbond_acc(nbond)+1
2666            endif
2667 ! Check against conformation repetitions.
2668           irepet=conf_comp(varia,etot)
2669           if (nrepm.gt.maxrepm) then
2670            if (print_mc.gt.0) &
2671             write (iout,'(a)') 'Too many conformation repetitions.'
2672             finish=.true.
2673            endif
2674 ! Store the accepted conf. and its energy.
2675            eold=etot
2676            do i=1,nvar
2677             varold(i)=varia(i)
2678            enddo
2679            if (irepet.eq.0) call zapis(varia,etot)
2680 ! Lower the temperature, if necessary.
2681            call cool
2682           else
2683            ntrial=ntrial+1
2684          endif ! accepted
2685    30    continue
2686          if(finish.and.imtasks_n.eq.0)exit LOOP2
2687         enddo LOOP2 ! accepted
2688 ! Check for time limit.
2689         not_done = (it.lt.max_mcm_it) .and. (nacc_tot.lt.maxacc)
2690         if(.not.not_done .or. finish) then
2691          if(imtasks_n.gt.0) then
2692           not_done=.true.
2693          else
2694           not_done=.false.
2695          endif
2696          finish=.true.
2697         endif
2698       enddo LOOP1 ! not_done
2699       runtime=tcpu()
2700       if (print_mc.gt.0) then
2701         write (iout,'(/80(1h*)/20x,a)') 'Summary run statistics:'
2702         call statprint(nacc,nfun,iretcode,etot,elowest)
2703         write (iout,'(a)') &
2704        'Statistics of multiple-bond motions. Total motions:'
2705         write (iout,'(16i5)') (nbond_move(i),i=1,Nbm)
2706         write (iout,'(a)') 'Accepted motions:'
2707         write (iout,'(16i5)') (nbond_acc(i),i=1,Nbm)
2708         if (it.ge.maxacc) &
2709       write (iout,'(/80(1h*)/20x,a/80(1h*)/)') 'All iterations done.'
2710       endif
2711 #ifdef AIX
2712       call flush_(iout)
2713 #else
2714       call flush(iout)
2715 #endif
2716       do is=1,nodenum-1
2717         call MPI_SEND(999, 1, MPI_INTEGER, is, tag,&
2718                    CG_COMM, ierr)
2719       enddo
2720       return
2721       end subroutine do_mcm
2722 !-----------------------------------------------------------------------------
2723       subroutine execute_slave(nodeinfo,iprint)
2724
2725       use MPI_data
2726       use minimm, only:minimize
2727 !      use minim
2728 !      implicit real*8 (a-h,o-z)
2729 !      include 'DIMENSIONS'
2730       include 'mpif.h'
2731 !      include 'COMMON.TIME1'
2732 !      include 'COMMON.IOUNITS'
2733 !rc      include 'COMMON.DEFORM'
2734 !rc      include 'COMMON.DEFORM1'
2735 !rc      include 'COMMON.DEFORM2'
2736 !      include 'COMMON.LOCAL'
2737 !      include 'COMMON.VAR'
2738 !      include 'COMMON.INFO'
2739 !      include 'COMMON.MINIM'
2740       character(len=10) :: nodeinfo
2741       real(kind=8),dimension(6*nres) :: x,x1    !(maxvar) (maxvar=6*maxres)
2742       integer :: nfun,iprint,i_switch,ierr,i_grnum_d,i_ennum_d,&
2743         i_hesnum_d,iitt,iretcode,iminrep
2744       real(kind=8) :: ener,energyx
2745
2746       nodeinfo='chujwdupe'
2747 !      print *,'Processor:',MyID,' Entering execute_slave'
2748       tag=0
2749 !      call MPI_SEND(nodeinfo, 10, MPI_CHARACTER, 0, tag,
2750 !     &              CG_COMM, ierr)
2751
2752 1001  call MPI_RECV(i_switch, 1, MPI_INTEGER, 0, tag,&
2753                     CG_COMM, status, ierr)
2754 !      write(iout,*)'12: tag=',tag
2755       if(iprint.ge.2)print *, MyID,' recv order ',i_switch
2756       if (i_switch.eq.12) then
2757        i_grnum_d=0
2758        i_ennum_d=0
2759        i_hesnum_d=0
2760        call MPI_RECV(iitt, 1, MPI_INTEGER, 0, tag,&
2761                      CG_COMM, status, ierr)
2762 !      write(iout,*)'12: tag=',tag
2763        call MPI_RECV(ener, 1, MPI_DOUBLE_PRECISION, 0, tag,&
2764                      CG_COMM, status, ierr)
2765 !      write(iout,*)'ener: tag=',tag
2766        call MPI_RECV(x(1), nvar, MPI_DOUBLE_PRECISION, 0, tag,&
2767                      CG_COMM, status, ierr)
2768 !      write(iout,*)'x: tag=',tag
2769        call MPI_RECV(x1(1), nvar, MPI_DOUBLE_PRECISION, 0, tag,&
2770                      CG_COMM, status, ierr)
2771 !      write(iout,*)'x1: tag=',tag
2772 #ifdef AIX
2773        call flush_(iout)
2774 #else
2775        call flush(iout)
2776 #endif
2777 !       print *,'calling minimize'
2778        call minimize(energyx,x,iretcode,nfun)
2779        if(iprint.gt.0) &
2780         write(iout,100)'minimized energy = ',energyx,&
2781           ' # funeval:',nfun,' iret ',iretcode
2782         write(*,100)'minimized energy = ',energyx,&
2783           ' # funeval:',nfun,' iret ',iretcode
2784  100   format(a20,f10.5,a12,i5,a6,i2)
2785        if(iretcode.eq.10) then
2786          do iminrep=2,3
2787           if(iprint.gt.1) &
2788           write(iout,*)' ... not converged - trying again ',iminrep
2789           call minimize(energyx,x,iretcode,nfun)
2790           if(iprint.gt.1) &
2791           write(iout,*)'minimized energy = ',energyx,&
2792            ' # funeval:',nfun,' iret ',iretcode
2793           if(iretcode.ne.10)go to 812
2794          enddo
2795          if(iretcode.eq.10) then
2796           if(iprint.gt.1) &
2797           write(iout,*)' ... not converged again - giving up'
2798           go to 812
2799          endif
2800        endif
2801 812    continue
2802 !       print *,'Sending results'
2803        call MPI_SEND(iitt, 1, MPI_INTEGER, 0, tag,&
2804                     CG_COMM, ierr)
2805        call MPI_SEND(ener, 1, MPI_DOUBLE_PRECISION, 0, tag,&
2806                     CG_COMM, ierr)
2807        call MPI_SEND(energyx, 1, MPI_DOUBLE_PRECISION, 0, tag,&
2808                     CG_COMM, ierr)
2809        call MPI_SEND(x(1), nvar, MPI_DOUBLE_PRECISION, 0, tag,&
2810                     CG_COMM, ierr)
2811        call MPI_SEND(x1(1), nvar, MPI_DOUBLE_PRECISION, 0, tag,&
2812                     CG_COMM, ierr)
2813        call MPI_SEND(i_grnum_d, 1, MPI_INTEGER, 0, tag,&
2814                     CG_COMM, ierr)
2815        call MPI_SEND(nfun, 1, MPI_INTEGER, 0, tag,&
2816                     CG_COMM, ierr)
2817        call MPI_SEND(i_hesnum_d, 1, MPI_INTEGER, 0, tag,&
2818                     CG_COMM, ierr)
2819 !       print *,'End sending'
2820        go to 1001
2821       endif
2822
2823       return
2824       end subroutine execute_slave
2825 #endif
2826 !-----------------------------------------------------------------------------
2827       subroutine heat(over)
2828
2829 !      implicit real*8 (a-h,o-z)
2830 !      include 'DIMENSIONS'
2831 !      include 'COMMON.MCM'
2832 !      include 'COMMON.IOUNITS'
2833       logical :: over
2834 ! Check if there`s a need to increase temperature.
2835       if (ntrial.gt.maxtrial) then
2836         if (NstepH.gt.0) then
2837           if (dabs(Tcur-TMax).lt.1.0D-7) then
2838             if (print_mc.gt.0) &
2839             write (iout,'(/80(1h*)/a,f8.3,a/80(1h*))') &
2840             'Upper limit of temperature reached. Terminating.'
2841             over=.true.
2842             Tcur=Tmin
2843           else
2844             Tcur=Tcur*TstepH
2845             if (Tcur.gt.Tmax) Tcur=Tmax
2846             betbol=1.0D0/(Rbol*Tcur)
2847             if (print_mc.gt.0) &
2848             write (iout,'(/80(1h*)/a,f8.3,a,f10.5/80(1h*))') &
2849             'System heated up to ',Tcur,' K; BetBol:',betbol
2850             ntherm=ntherm+1
2851             ntrial=0
2852             over=.false.
2853           endif
2854         else
2855          if (print_mc.gt.0) &
2856          write (iout,'(a)') &
2857        'Maximum number of trials in a single MCM iteration exceeded.'
2858          over=.true.
2859          Tcur=Tmin
2860         endif
2861       else
2862         over=.false.
2863       endif
2864       return
2865       end subroutine heat
2866 !-----------------------------------------------------------------------------
2867       subroutine cool
2868
2869 !      implicit real*8 (a-h,o-z)
2870 !      include 'DIMENSIONS'
2871 !      include 'COMMON.MCM'
2872 !      include 'COMMON.IOUNITS'
2873       if (nstepC.gt.0 .and. dabs(Tcur-Tmin).gt.1.0D-7) then
2874         Tcur=Tcur/TstepC
2875         if (Tcur.lt.Tmin) Tcur=Tmin
2876         betbol=1.0D0/(Rbol*Tcur)
2877         if (print_mc.gt.0) &
2878         write (iout,'(/80(1h*)/a,f8.3,a,f10.5/80(1h*))') &
2879         'System cooled down up to ',Tcur,' K; BetBol:',betbol
2880       endif
2881       return
2882       end subroutine cool
2883 !-----------------------------------------------------------------------------
2884       subroutine perturb(error,lprint,MoveType,nbond,max_phi)
2885
2886       use geometry
2887       use energy, only:nnt,nct,itype
2888       use md_calc, only:bond_move
2889 !      implicit real*8 (a-h,o-z)
2890 !      include 'DIMENSIONS'
2891       integer,parameter :: MMaxSideMove=100
2892 !      include 'COMMON.MCM'
2893 !      include 'COMMON.CHAIN'
2894 !      include 'COMMON.INTERACT'
2895 !      include 'COMMON.VAR'
2896 !      include 'COMMON.GEO'
2897 !      include 'COMMON.NAMES'
2898 !      include 'COMMON.IOUNITS'
2899 !rc      include 'COMMON.DEFORM1'
2900       logical :: error,lprint,fail
2901       integer :: MoveType,nbond,end_select,ind_side(MMaxSideMove)
2902       real(kind=8) :: max_phi
2903       real(kind=8) :: psi!,gen_psi
2904 !el      external iran_num
2905 !el      integer iran_num
2906       integer :: ifour
2907
2908 !!! local variables - el
2909       integer :: itrial,iiwin,iwindow,isctry,i,icount,j,nstart,&
2910             nside_move,inds,indx,ii,iti
2911       real(kind=8) :: bond_prob,theta_new
2912
2913       data ifour /4/
2914       error=.false.
2915       lprint=.false.
2916 ! Perturb the conformation according to a randomly selected move.
2917       call SelectMove(MoveType)
2918 !      write (iout,*) 'MoveType=',MoveType
2919       itrial=0
2920       goto (100,200,300,400,500) MoveType
2921 !------------------------------------------------------------------------------
2922 ! Backbone N-bond move.
2923 ! Select the number of bonds (length of the segment to perturb).
2924   100 continue
2925       if (itrial.gt.1000) then
2926         write (iout,'(a)') 'Too many attempts at multiple-bond move.'
2927         error=.true.
2928         return
2929       endif
2930       bond_prob=ran_number(0.0D0,sumpro_bond(nbm))
2931 !     print *,'sumpro_bond(nbm)=',sumpro_bond(nbm),
2932 !    & ' Bond_prob=',Bond_Prob
2933       do i=1,nbm-1
2934 !       print *,i,Bond_Prob,sumpro_bond(i),sumpro_bond(i+1)
2935         if (bond_prob.ge.sumpro_bond(i) .and. &
2936                      bond_prob.le.sumpro_bond(i+1)) then
2937           nbond=i+1
2938           goto 10
2939         endif
2940       enddo
2941       write (iout,'(2a)') 'In PERTURB: Error - number of bonds',&
2942                           ' to move out of range.'
2943       error=.true.
2944       return
2945    10 continue
2946       if (nwindow.gt.0) then
2947 ! Select the first residue to perturb
2948         iwindow=iran_num(1,nwindow)
2949         print *,'iwindow=',iwindow
2950         iiwin=1
2951         do while (winlen(iwindow).lt.nbond)
2952           iwindow=iran_num(1,nwindow)
2953           iiwin=iiwin+1
2954           if (iiwin.gt.1000) then
2955              write (iout,'(a)') 'Cannot select moveable residues.'
2956              error=.true.
2957              return
2958           endif
2959         enddo
2960         nstart=iran_num(winstart(iwindow),winend(iwindow))
2961       else
2962         nstart = iran_num(koniecl+2,nres-nbond-koniecl)
2963 !d      print *,'nres=',nres,' nbond=',nbond,' koniecl=',koniecl,
2964 !d   &        ' nstart=',nstart
2965       endif
2966       psi = gen_psi()
2967       if (psi.eq.0.0) then
2968         error=.true.
2969         return
2970       endif
2971       if (print_mc.gt.1) write (iout,'(a,i4,a,i4,a,f8.3)') &
2972        'PERTURB: nbond=',nbond,' nstart=',nstart,' psi=',psi*rad2deg
2973 !d    print *,'nstart=',nstart
2974       call bond_move(nbond,nstart,psi,.false.,error)
2975       if (error) then
2976         write (iout,'(2a)') &
2977        'Could not define reference system in bond_move, ',&
2978        'choosing ahother segment.'
2979         itrial=itrial+1
2980         goto 100
2981       endif
2982       nbond_move(nbond)=nbond_move(nbond)+1
2983       moves(1)=moves(1)+1
2984       nmove=nmove+1
2985       return
2986 !------------------------------------------------------------------------------
2987 ! Backbone endmove. Perturb a SINGLE angle of a residue close to the end of
2988 ! the chain.
2989   200 continue
2990       lprint=.true.
2991 !     end_select=iran_num(1,2*koniecl)
2992 !     if (end_select.gt.koniecl) then
2993 !       end_select=nphi-(end_select-koniecl)
2994 !     else
2995 !       end_select=koniecl+3
2996 !     endif
2997 !     if (nwindow.gt.0) then
2998 !       iwin=iran_num(1,nwindow)
2999 !       i1=max0(4,winstart(iwin))
3000 !       i2=min0(winend(imin)+2,nres)
3001 !       end_select=iran_num(i1,i2)
3002 !     else
3003 !      iselect = iran_num(1,nmov_var)
3004 !      jj = 0
3005 !      do i=1,nphi
3006 !        if (isearch_tab(i).eq.1) jj = jj+1
3007 !        if (jj.eq.iselect) then
3008 !          end_select=i+3
3009 !          exit
3010 !        endif
3011 !      enddo
3012 !     endif
3013       end_select = iran_num(4,nres)
3014       psi=max_phi*gen_psi()
3015       if (psi.eq.0.0D0) then
3016         error=.true.
3017         return
3018       endif
3019       phi(end_select)=pinorm(phi(end_select)+psi)
3020       if (print_mc.gt.1) write (iout,'(a,i4,a,f8.3,a,f8.3)') &
3021        'End angle',end_select,' moved by ',psi*rad2deg,' new angle:',&
3022        phi(end_select)*rad2deg
3023 !     if (end_select.gt.3)
3024 !    &   theta(end_select-1)=gen_theta(itype(end_select-2),
3025 !    &                              phi(end_select-1),phi(end_select))
3026 !     if (end_select.lt.nres)
3027 !    &    theta(end_select)=gen_theta(itype(end_select-1),
3028 !    &                              phi(end_select),phi(end_select+1))
3029 !d    print *,'nres=',nres,' end_select=',end_select
3030 !d    print *,'theta',end_select-1,theta(end_select-1)
3031 !d    print *,'theta',end_select,theta(end_select)
3032       moves(2)=moves(2)+1
3033       nmove=nmove+1
3034       lprint=.false.
3035       return
3036 !------------------------------------------------------------------------------
3037 ! Side chain move.
3038 ! Select the number of SCs to perturb.
3039   300 isctry=0
3040   301 nside_move=iran_num(1,MaxSideMove)
3041 !     print *,'nside_move=',nside_move,' MaxSideMove',MaxSideMove
3042 ! Select the indices.
3043       do i=1,nside_move
3044         icount=0
3045   111   inds=iran_num(nnt,nct)
3046         icount=icount+1
3047         if (icount.gt.1000) then
3048           write (iout,'(a)')'Error - cannot select side chains to move.'
3049           error=.true.
3050           return
3051         endif
3052         if (itype(inds,1).eq.10) goto 111
3053         do j=1,i-1
3054           if (inds.eq.ind_side(j)) goto 111
3055         enddo
3056         do j=1,i-1
3057           if (inds.lt.ind_side(j)) then
3058             indx=j
3059             goto 112
3060           endif
3061         enddo
3062         indx=i
3063   112   do j=i,indx+1,-1
3064           ind_side(j)=ind_side(j-1)
3065         enddo
3066   113   ind_side(indx)=inds
3067       enddo
3068 ! Carry out perturbation.
3069       do i=1,nside_move
3070         ii=ind_side(i)
3071         iti=itype(ii,1)
3072         call gen_side(iti,theta(ii+1),alph(ii),omeg(ii),fail,1)
3073         if (fail) then
3074           isctry=isctry+1
3075           if (isctry.gt.1000) then
3076             write (iout,'(a)') 'Too many errors in SC generation.'
3077             error=.true.
3078             return
3079           endif
3080           goto 301
3081         endif
3082         if (print_mc.gt.1) write (iout,'(2a,i4,a,2f8.3)') &
3083          'Side chain ',restyp(iti,1),ii,' moved to ',&
3084          alph(ii)*rad2deg,omeg(ii)*rad2deg
3085       enddo
3086       moves(3)=moves(3)+1
3087       nmove=nmove+1
3088       return
3089 !------------------------------------------------------------------------------
3090 ! THETA move
3091   400 end_select=iran_num(3,nres)
3092       theta_new=gen_theta(itype(end_select,1),phi(end_select),&
3093                           phi(end_select+1),1)
3094       if (print_mc.gt.1) write (iout,'(a,i3,a,f8.3,a,f8.3)') &
3095        'Theta ',end_select,' moved from',theta(end_select)*rad2deg,&
3096        ' to ',theta_new*rad2deg
3097       theta(end_select)=theta_new
3098       moves(4)=moves(4)+1
3099       nmove=nmove+1
3100       return
3101 !------------------------------------------------------------------------------
3102 ! Error returned from SelectMove.
3103   500 error=.true.
3104       return
3105       end subroutine perturb
3106 !-----------------------------------------------------------------------------
3107       subroutine SelectMove(MoveType)
3108
3109 !      implicit real*8 (a-h,o-z)
3110 !      include 'DIMENSIONS'
3111 !      include 'COMMON.MCM'
3112 !      include 'COMMON.IOUNITS'
3113
3114 !!! local variables - el
3115       integer :: i,MoveType
3116       real(kind=8) :: what_move
3117
3118       what_move=ran_number(0.0D0,sumpro_type(MaxMoveType))
3119       do i=1,MaxMoveType
3120         if (what_move.ge.sumpro_type(i-1).and. &
3121                   what_move.lt.sumpro_type(i)) then
3122           MoveType=i
3123           return
3124         endif
3125       enddo
3126       write (iout,'(a)') &
3127        'Fatal error in SelectMoveType: cannot select move.'
3128       MoveType=MaxMoveType+1
3129       return
3130       end subroutine SelectMove
3131 !-----------------------------------------------------------------------------
3132       real(kind=8) function gen_psi()
3133
3134       use geometry_data, only: angmin,pi
3135 !el      implicit none
3136       integer :: i
3137       real(kind=8) :: x !,ran_number
3138 !      include 'COMMON.IOUNITS'
3139 !      include 'COMMON.GEO'
3140       x=0.0D0
3141       do i=1,100
3142         x=ran_number(-pi,pi)
3143         if (dabs(x).gt.angmin) then
3144           gen_psi=x
3145           return
3146         endif
3147       enddo
3148       write (iout,'(a)')'From Gen_Psi: Cannot generate angle increment.'
3149       gen_psi=0.0D0
3150       return
3151       end function gen_psi
3152 !-----------------------------------------------------------------------------
3153       subroutine metropolis(n,xcur,xold,ecur,eold,accepted,similar,enelower)
3154
3155 !      implicit real*8 (a-h,o-z)
3156 !      include 'DIMENSIONS'
3157 !      include 'COMMON.MCM'
3158 !      include 'COMMON.IOUNITS'
3159 !      include 'COMMON.VAR'
3160 !      include 'COMMON.GEO'
3161 !rc      include 'COMMON.DEFORM'
3162       integer :: n
3163       real(kind=8) :: ecur,eold,xx,bol  !,ran_number
3164       real(kind=8),dimension(n) :: xcur,xold
3165       real(kind=8) :: ecut1 ,ecut2 ,tola
3166       logical :: accepted,similar,not_done,enelower
3167       logical :: lprn
3168
3169 !!! local variables - el
3170       real(kind=8) :: xxh,difene,reldife
3171
3172       data ecut1 /-1.0D-5/,ecut2 /5.0D-3/,tola/5.0D0/
3173 !      ecut1=-5*enedif
3174 !      ecut2=50*enedif
3175 !      tola=5.0d0
3176 ! Set lprn=.true. for debugging.
3177       lprn=.false.
3178       if (lprn) &
3179 !el      write(iout,*)'enedif',enedif,' ecut1',ecut1,' ecut2',ecut2
3180       write(iout,*)' ecut1',ecut1,' ecut2',ecut2
3181       similar=.false.
3182       enelower=.false.
3183       accepted=.false.
3184 ! Check if the conformation is similar.
3185       difene=ecur-eold
3186       reldife=difene/dmax1(dabs(eold),dabs(ecur),1.0D0)
3187       if (lprn) then
3188         write (iout,*) 'Metropolis'
3189         write(iout,*)'ecur,eold,difene,reldife',ecur,eold,difene,reldife
3190       endif
3191 ! If energy went down remarkably, we accept the new conformation
3192 ! unconditionally.
3193 !jp      if (reldife.lt.ecut1) then
3194       if (difene.lt.ecut1) then
3195         accepted=.true.
3196         EneLower=.true.
3197         if (lprn) write (iout,'(a)') &
3198          'Conformation accepted, because energy has lowered remarkably.'
3199 !      elseif (reldife.lt.ecut2 .and. dif_ang(nphi,xcur,xold).lt.tola)
3200 !jp      elseif (reldife.lt.ecut2)
3201       elseif (difene.lt.ecut2) &
3202        then
3203 ! Reject the conf. if energy has changed insignificantly and there is not
3204 ! much change in conformation.
3205         if (lprn) &
3206          write (iout,'(2a)') 'Conformation rejected, because it is',&
3207             ' similar to the preceding one.'
3208         accepted=.false.
3209         similar=.true.
3210       else
3211 ! Else carry out Metropolis test.
3212         EneLower=.false.
3213         xx=ran_number(0.0D0,1.0D0)
3214         xxh=betbol*difene
3215         if (lprn) &
3216           write (iout,*) 'betbol=',betbol,' difene=',difene,' xxh=',xxh
3217         if (xxh.gt.50.0D0) then
3218           bol=0.0D0
3219         else
3220           bol=exp(-xxh)
3221         endif
3222         if (lprn) write (iout,*) 'bol=',bol,' xx=',xx
3223         if (bol.gt.xx) then
3224           accepted=.true.
3225           if (lprn) write (iout,'(a)') &
3226           'Conformation accepted, because it passed Metropolis test.'
3227         else
3228           accepted=.false.
3229           if (lprn) write (iout,'(a)') &
3230        'Conformation rejected, because it did not pass Metropolis test.'
3231         endif
3232       endif
3233 #ifdef AIX
3234       call flush_(iout)
3235 #else
3236       call flush(iout)
3237 #endif
3238       return
3239       end subroutine metropolis
3240 !-----------------------------------------------------------------------------
3241       integer function conf_comp(x,ene)
3242
3243       use geometry_data, only: nphi
3244 !      implicit real*8 (a-h,o-z)
3245 !      include 'DIMENSIONS'
3246 !      include 'COMMON.MCM'
3247 !      include 'COMMON.VAR'
3248 !      include 'COMMON.IOUNITS'
3249 !      include 'COMMON.GEO'
3250       real(kind=8) :: etol, angtol
3251       real(kind=8),dimension(6*nres) :: x       !(maxvar) (maxvar=6*maxres)
3252       real(kind=8) :: difa      !dif_ang,
3253
3254 !!! local variables - el
3255       integer :: ii,i
3256       real(kind=8) :: ene
3257
3258       data etol /0.1D0/, angtol /20.0D0/
3259       do ii=nsave,1,-1
3260 !       write (iout,*) 'ii=',ii,'ene=',ene,esave(ii),dabs(ene-esave(ii))
3261         if (dabs(ene-esave(ii)).lt.etol) then
3262           difa=dif_ang(nphi,x,varsave(1,ii))
3263 !         do i=1,nphi
3264 !           write(iout,'(i3,3f8.3)')i,rad2deg*x(i),
3265 !    &          rad2deg*varsave(i,ii)
3266 !         enddo
3267 !         write(iout,*) 'ii=',ii,' difa=',difa,' angtol=',angtol
3268           if (difa.le.angtol) then
3269             if (print_mc.gt.0) then
3270             write (iout,'(a,i5,2(a,1pe15.4))') &
3271             'Current conformation matches #',ii,&
3272             ' in the store array ene=',ene,' esave=',esave(ii)
3273 !           write (*,'(a,i5,a)') 'Current conformation matches #',ii,
3274 !    &      ' in the store array.'
3275             endif ! print_mc.gt.0
3276             if (print_mc.gt.1) then
3277             do i=1,nphi
3278               write(iout,'(i3,3f8.3)')i,rad2deg*x(i),&
3279                   rad2deg*varsave(i,ii)
3280             enddo
3281             endif ! print_mc.gt.1
3282             nrepm=nrepm+1
3283             conf_comp=ii
3284             return
3285           endif
3286         endif
3287       enddo
3288       conf_comp=0
3289       return
3290       end function conf_comp
3291 !-----------------------------------------------------------------------------
3292       real(kind=8) function dif_ang(n,x,y)
3293
3294       use geometry_data, only: dwapi
3295 !el      implicit none
3296       integer :: i,n
3297       real(kind=8),dimension(n) :: x,y
3298       real(kind=8) :: w,wa,dif,difa
3299 !el      real(kind=8) :: pinorm
3300 !      include 'COMMON.GEO'
3301       wa=0.0D0
3302       difa=0.0D0
3303       do i=1,n
3304         dif=dabs(pinorm(y(i)-x(i)))
3305         if (dabs(dif-dwapi).lt.dif) dif=dabs(dif-dwapi)
3306         w=1.0D0-(2.0D0*(i-1)/(n-1)-1.0D0)**2+1.0D0/n
3307         wa=wa+w
3308         difa=difa+dif*dif*w
3309       enddo
3310       dif_ang=rad2deg*dsqrt(difa/wa)
3311       return
3312       end function dif_ang
3313 !-----------------------------------------------------------------------------
3314       subroutine add2cache(n1,n2,ncache,nvar,SourceID,CachSrc,ecur,xcur,ecache,xcache)
3315
3316 !      implicit none
3317 !      include 'COMMON.GEO'
3318 !      include 'COMMON.IOUNITS'
3319       integer :: n1,n2,ncache,nvar,SourceID,CachSrc(n2)
3320       integer :: i,ii,j
3321       real(kind=8) :: ecur,xcur(nvar),ecache(n2),xcache(n1,n2)
3322 !d    write (iout,*) 'Enter ADD2CACHE ncache=',ncache ,' ecur',ecur
3323 !d    write (iout,'(10f8.3)') (rad2deg*xcur(i),i=1,nvar)
3324 !d    write (iout,*) 'Old CACHE array:'
3325 !d    do i=1,ncache
3326 !d      write (iout,*) 'i=',i,' ecache=',ecache(i),' CachSrc',CachSrc(i)
3327 !d      write (iout,'(10f8.3)') (rad2deg*xcache(j,i),j=1,nvar)
3328 !d    enddo
3329
3330       i=ncache
3331       do while (i.gt.0 .and. ecur.lt.ecache(i))
3332         i=i-1
3333       enddo
3334       i=i+1
3335 !d    write (iout,*) 'i=',i,' ncache=',ncache
3336       if (ncache.eq.n2) then
3337         write (iout,*) 'Cache dimension exceeded',ncache,n2
3338         write (iout,*) 'Highest-energy conformation will be removed.'
3339         ncache=ncache-1
3340       endif
3341       do ii=ncache,i,-1
3342         ecache(ii+1)=ecache(ii)
3343         CachSrc(ii+1)=CachSrc(ii)
3344         do j=1,nvar
3345           xcache(j,ii+1)=xcache(j,ii)
3346         enddo
3347       enddo
3348       ecache(i)=ecur
3349       CachSrc(i)=SourceID
3350       do j=1,nvar
3351         xcache(j,i)=xcur(j)
3352       enddo
3353       ncache=ncache+1
3354 !d    write (iout,*) 'New CACHE array:'
3355 !d    do i=1,ncache
3356 !d      write (iout,*) 'i=',i,' ecache=',ecache(i),' CachSrc',CachSrc(i)
3357 !d      write (iout,'(10f8.3)') (rad2deg*xcache(j,i),j=1,nvar)
3358 !d    enddo
3359       return
3360       end subroutine add2cache
3361 !-----------------------------------------------------------------------------
3362       subroutine rm_from_cache(i,n1,n2,ncache,nvar,CachSrc,ecache,xcache)
3363
3364 !      implicit none
3365 !      include 'COMMON.GEO'
3366 !      include 'COMMON.IOUNITS'
3367       integer :: n1,n2,ncache,nvar,CachSrc(n2)
3368       integer :: i,ii,j
3369       real(kind=8) :: ecache(n2),xcache(n1,n2)
3370
3371 !d    write (iout,*) 'Enter RM_FROM_CACHE'
3372 !d    write (iout,*) 'Old CACHE array:'
3373 !d    do ii=1,ncache
3374 !d    write (iout,*)'i=',ii,' ecache=',ecache(ii),' CachSrc',CachSrc(ii)
3375 !d      write (iout,'(10f8.3)') (rad2deg*xcache(j,ii),j=1,nvar)
3376 !d    enddo
3377
3378       do ii=i+1,ncache
3379         ecache(ii-1)=ecache(ii)
3380         CachSrc(ii-1)=CachSrc(ii)
3381         do j=1,nvar
3382           xcache(j,ii-1)=xcache(j,ii)
3383         enddo
3384       enddo
3385       ncache=ncache-1
3386 !d    write (iout,*) 'New CACHE array:'
3387 !d    do i=1,ncache
3388 !d      write (iout,*) 'i=',i,' ecache=',ecache(i),' CachSrc',CachSrc(i)
3389 !d      write (iout,'(10f8.3)') (rad2deg*xcache(j,i),j=1,nvar)
3390 !d    enddo
3391       return
3392       end subroutine rm_from_cache
3393 !-----------------------------------------------------------------------------
3394 ! mcm.F         io_mcm
3395 !-----------------------------------------------------------------------------
3396       subroutine statprint(it,nfun,iretcode,etot,elowest)
3397
3398       use control_data, only: MaxMoveType,minim
3399       use control, only: tcpu
3400       use mcm_data
3401 !      implicit real*8 (a-h,o-z)
3402 !      include 'DIMENSIONS'
3403 !      include 'COMMON.IOUNITS'
3404 !      include 'COMMON.CONTROL'
3405 !      include 'COMMON.MCM'
3406 !el local variables
3407       integer :: it,nfun,iretcode,i
3408       real(kind=8) :: etot,elowest,fr_mov_i
3409
3410       if (minim) then
3411         write (iout,&
3412         '(80(1h*)/a,i5,a,1pe14.5,a,1pe14.5/a,i3,a,i10,a,i5,a,i5)') &
3413         'Finished iteration #',it,' energy is',etot,&
3414         ' lowest energy:',elowest,&
3415         'SUMSL return code:',iretcode,&
3416         ' # of energy evaluations:',neneval,&
3417         '# of temperature jumps:',ntherm,&
3418         ' # of minima repetitions:',nrepm
3419       else
3420         write (iout,'(80(1h*)/a,i8,a,1pe14.5,a,1pe14.5)') &
3421         'Finished iteration #',it,' energy is',etot,&
3422         ' lowest energy:',elowest
3423       endif
3424       write (iout,'(/4a)') &
3425        'Kind of move   ','           total','       accepted',&
3426        '  fraction'
3427       write (iout,'(58(1h-))')
3428       do i=-1,MaxMoveType
3429         if (moves(i).eq.0) then
3430           fr_mov_i=0.0d0
3431         else
3432           fr_mov_i=dfloat(moves_acc(i))/dfloat(moves(i))
3433         endif
3434         write(iout,'(a,2i15,f10.5)')MovTypID(i),moves(i),moves_acc(i),&
3435                fr_mov_i
3436       enddo
3437       write (iout,'(a,2i15,f10.5)') 'total           ',nmove,nacc_tot,&
3438                dfloat(nacc_tot)/dfloat(nmove)
3439       write (iout,'(58(1h-))')
3440       write (iout,'(a,1pe12.4)') 'Elapsed time:',tcpu()
3441       return
3442       end subroutine statprint
3443 !-----------------------------------------------------------------------------
3444       subroutine zapis(varia,etot)
3445
3446       use geometry_data, only: nres,rad2deg,nvar
3447       use mcm_data
3448       use MPI_data
3449 !      implicit real*8 (a-h,o-z)
3450 !      include 'DIMENSIONS'
3451 #ifdef MPI
3452       use MPI_data      !include 'COMMON.INFO'
3453       include 'mpif.h'
3454 #endif
3455 !      include 'COMMON.GEO'
3456 !      include 'COMMON.VAR'
3457 !      include 'COMMON.MCM'
3458 !      include 'COMMON.IOUNITS'
3459       integer,dimension(nsave) :: itemp
3460       real(kind=8),dimension(6*nres) :: varia   !(maxvar)       (maxvar=6*maxres)
3461       logical :: lprint
3462 !el local variables
3463       integer :: j,i,maxvar
3464       real(kind=8) :: etot
3465
3466 !el      allocate(esave(nsave)) !(maxsave)
3467
3468       maxvar=6*nres
3469       lprint=.false.
3470       if (lprint) then
3471       write (iout,'(a,i5,a,i5)') 'Enter ZAPIS NSave=',Nsave,&
3472         ' MaxSave=',MaxSave
3473       write (iout,'(a)') 'Current energy and conformation:'
3474       write (iout,'(1pe14.5)') etot
3475       write (iout,'(10f8.3)') (rad2deg*varia(i),i=1,nvar)
3476       endif
3477 ! Shift the contents of the esave and varsave arrays if filled up.
3478       call add2cache(6*nres,nsave,nsave,nvar,MyID,itemp,&
3479                      etot,varia,esave,varsave)
3480       if (lprint) then
3481       write (iout,'(a)') 'Energies and the VarSave array.'
3482       do i=1,nsave
3483         write (iout,'(i5,1pe14.5)') i,esave(i)
3484         write (iout,'(10f8.3)') (rad2deg*varsave(j,i),j=1,nvar)
3485       enddo
3486       endif
3487       return
3488       end subroutine zapis
3489 !-----------------------------------------------------------------------------
3490 !-----------------------------------------------------------------------------
3491       subroutine alloc_MCM_arrays
3492
3493       use energy_data, only: max_ene
3494       use MPI_data
3495 ! common.mce
3496 !      common /mce/
3497       allocate(entropy(-max_ene-4:max_ene)) !(-max_ene-4:max_ene)
3498       allocate(nhist(-max_ene:max_ene)) !(-max_ene:max_ene)
3499       allocate(nminima(maxsave)) !(maxsave)
3500 !      common /pool/
3501       allocate(xpool(6*nres,max_pool)) !(maxvar,max_pool)(maxvar=6*maxres)
3502       allocate(epool(max_pool)) !(max_pool)
3503 ! commom.mcm
3504 !      common /mcm/
3505       if(.not.allocated(nsave_part)) allocate(nsave_part(nctasks)) !(max_cg_procs)
3506 !      common /move/
3507 ! in io: mcmread
3508 !      real(kind=8),dimension(:),allocatable :: sumpro_type !(0:MaxMoveType)
3509       allocate(sumpro_bond(0:nres)) !(0:maxres)
3510       allocate(nbond_move(nres),nbond_acc(nres)) !(maxres)
3511       allocate(moves(-1:MaxMoveType+1),moves_acc(-1:MaxMoveType+1)) !(-1:MaxMoveType+1)
3512 !      common /accept_stats/
3513 !      allocate(nacc_part !(0:MaxProcs) !el nie uzywane???
3514 !      common /windows/ in io: mcmread
3515 !      allocate(winstart,winend,winlen !(maxres)
3516 !      common /moveID/
3517 !el      allocate(MovTypID(-1:MaxMoveType+1)) !(-1:MaxMoveType+1)
3518 ! common.var
3519 !      common /oldgeo/
3520       allocate(varsave(nres*6,maxsave))  !(maxvar,maxsave)(maxvar=6*maxres)
3521       allocate(esave(maxsave)) !(maxsave)
3522       allocate(Origin(maxsave)) !(maxsave)
3523
3524       return
3525       end subroutine alloc_MCM_arrays
3526 !-----------------------------------------------------------------------------
3527 !-----------------------------------------------------------------------------
3528       end module mcm_md