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