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