source/unres/MCM_MD.F90

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