source/unres/src_MD-M/mcm.F

   1       subroutine mcm_setup
   2       implicit real*8 (a-h,o-z)
   3       include 'DIMENSIONS'
   4       include 'COMMON.IOUNITS'
   5       include 'COMMON.MCM'
   6       include 'COMMON.CONTROL'
   7       include 'COMMON.INTERACT'
   8       include 'COMMON.NAMES'
   9       include 'COMMON.CHAIN'
  10       include 'COMMON.VAR'
  11 C
  12 C Set up variables used in MC/MCM.
  13 C
  14       write (iout,'(80(1h*)/20x,a/80(1h*))') 'MCM control parameters:'
  15       write (iout,'(5(a,i7))') 'Maxacc:',maxacc,' MaxTrial:',MaxTrial,
  16      & ' MaxRepm:',MaxRepm,' MaxGen:',MaxGen,' MaxOverlap:',MaxOverlap
  17       write (iout,'(4(a,f8.1)/2(a,i3))')
  18      & 'Tmin:',Tmin,' Tmax:',Tmax,' TstepH:',TstepH,
  19      & ' TstepC:',TstepC,'NstepH:',NstepH,' NstepC:',NstepC
  20       if (nwindow.gt.0) then
  21         write (iout,'(a)') 'Perturbation windows:'
  22         do i=1,nwindow
  23           i1=winstart(i)
  24           i2=winend(i)
  25           it1=itype(i1)
  26           it2=itype(i2)
  27           write (iout,'(a,i3,a,i3,a,i3)') restyp(it1),i1,restyp(it2),i2,
  28      &                        ' length',winlen(i)
  29         enddo
  30       endif
  31 C Rbolt=8.3143D-3*2.388459D-01 kcal/(mol*K)
  32       RBol=1.9858D-3
  33 C Number of "end bonds".
  34       koniecl=0
  35 c     koniecl=nphi
  36       print *,'koniecl=',koniecl
  37       write (iout,'(a)') 'Probabilities of move types:'
  38       write (*,'(a)') 'Probabilities of move types:'
  39       do i=1,MaxMoveType
  40         write (iout,'(a,f10.5)') MovTypID(i),
  41      &    sumpro_type(i)-sumpro_type(i-1)
  42         write (*,'(a,f10.5)') MovTypID(i),
  43      &    sumpro_type(i)-sumpro_type(i-1)
  44       enddo
  45       write (iout,*)
  46 C Maximum length of N-bond segment to be moved
  47 c     nbm=nres-1-(2*koniecl-1)
  48       if (nwindow.gt.0) then
  49         maxwinlen=winlen(1)
  50         do i=2,nwindow
  51           if (winlen(i).gt.maxwinlen) maxwinlen=winlen(i)
  52         enddo
  53         nbm=min0(maxwinlen,6)
  54         write (iout,'(a,i3,a,i3)') 'Nbm=',Nbm,' Maxwinlen=',Maxwinlen
  55       else
  56         nbm=min0(6,nres-2)
  57       endif
  58       sumpro_bond(0)=0.0D0
  59       sumpro_bond(1)=0.0D0
  60       do i=2,nbm
  61         sumpro_bond(i)=sumpro_bond(i-1)+1.0D0/dfloat(i)
  62       enddo
  63       write (iout,'(a)') 'The SumPro_Bond array:'
  64       write (iout,'(8f10.5)') (sumpro_bond(i),i=1,nbm)
  65       write (*,'(a)') 'The SumPro_Bond array:'
  66       write (*,'(8f10.5)') (sumpro_bond(i),i=1,nbm)
  67 C Maximum number of side chains moved simultaneously
  68 c     print *,'nnt=',nnt,' nct=',nct
  69       ngly=0
  70       do i=nnt,nct
  71         if (itype(i).eq.10) ngly=ngly+1
  72       enddo
  73       mmm=nct-nnt-ngly+1
  74       if (mmm.gt.0) then
  75         MaxSideMove=min0((nct-nnt+1)/2,mmm)
  76       endif
  77 c     print *,'MaxSideMove=',MaxSideMove
  78 C Max. number of generated confs (not used at present).
  79       maxgen=10000
  80 C Set initial temperature
  81       Tcur=Tmin
  82       betbol=1.0D0/(Rbol*Tcur)
  83       write (iout,'(a,f8.1,a,f10.5)') 'Initial temperature:',Tcur,
  84      &    ' BetBol:',betbol
  85       write (iout,*) 'RanFract=',ranfract
  86       return
  87       end
  88 c------------------------------------------------------------------------------
  89 #ifndef MPI
  90       subroutine do_mcm(i_orig)
  91 C Monte-Carlo-with-Minimization calculations - serial code.
  92       implicit real*8 (a-h,o-z)
  93       include 'DIMENSIONS'
  94       include 'COMMON.IOUNITS'
  95       include 'COMMON.GEO'
  96       include 'COMMON.CHAIN'
  97       include 'COMMON.MCM'
  98       include 'COMMON.CONTACTS'
  99       include 'COMMON.CONTROL'
 100       include 'COMMON.VAR'
 101       include 'COMMON.INTERACT'
 102       include 'COMMON.CACHE'
 103 crc      include 'COMMON.DEFORM'
 104 crc      include 'COMMON.DEFORM1'
 105       include 'COMMON.NAMES'
 106       logical accepted,over,ovrtim,error,lprint,not_done,my_conf,
 107      &        enelower,non_conv
 108       integer MoveType,nbond,conf_comp
 109       integer ifeed(max_cache)
 110       double precision varia(maxvar),varold(maxvar),elowest,eold,
 111      & przes(3),obr(3,3)
 112       double precision energia(0:n_ene)
 113
 114 C---------------------------------------------------------------------------
 115 C Initialize counters.
 116 C---------------------------------------------------------------------------
 117 C Total number of generated confs.
 118       ngen=0
 119 C Total number of moves. In general this won't be equal to the number of
 120 C attempted moves, because we may want to reject some "bad" confs just by
 121 C overlap check.
 122       nmove=0
 123 C Total number of temperature jumps.
 124       ntherm=0
 125 C Total number of shift (nbond_move(1)), spike, crankshaft, three-bond,...
 126 C motions.
 127       ncache=0
 128       do i=1,nres
 129         nbond_move(i)=0
 130       enddo
 131 C Initialize total and accepted number of moves of various kind.
 132       do i=0,MaxMoveType
 133         moves(i)=0
 134         moves_acc(i)=0
 135       enddo
 136 C Total number of energy evaluations.
 137       neneval=0
 138       nfun=0
 139       nsave=0
 140
 141       write (iout,*) 'RanFract=',RanFract
 142
 143       WhatsUp=0
 144       Kwita=0
 145
 146 c----------------------------------------------------------------------------
 147 C Compute and print initial energies.
 148 c----------------------------------------------------------------------------
 149       call intout
 150       write (iout,'(/80(1h*)/a)') 'Initial energies:'
 151       call chainbuild
 152       nf=0
 153
 154       call etotal(energia(0))
 155       etot = energia(0)
 156 C Minimize the energy of the first conformation.
 157       if (minim) then
 158         call geom_to_var(nvar,varia)
 159 !       write (iout,*) 'The VARIA array'
 160 !       write (iout,'(8f10.4)') (rad2deg*varia(i),i=1,nvar)
 161         call minimize(etot,varia,iretcode,nfun)
 162         call var_to_geom(nvar,varia)
 163         call chainbuild
 164         write (iout,*) 'etot from MINIMIZE:',etot
 165 !       write (iout,*) 'Tha VARIA array'
 166 !       write (iout,'(8f10.4)') (rad2deg*varia(i),i=1,nvar)
 167
 168         call etotal(energia(0))
 169         etot=energia(0)
 170         call enerprint(energia(0))
 171       endif
 172       if (refstr) then
 173         call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup),nsup,przes,
 174      &             obr,non_conv)
 175         rms=dsqrt(rms)
 176         call contact(.false.,ncont,icont,co)
 177         frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
 178         write (iout,'(a,f8.3,a,f8.3,a,f8.3)')
 179      &    'RMS deviation from the reference structure:',rms,
 180      &    ' % of native contacts:',frac*100,' contact order:',co
 181         if (print_stat)
 182      &  write (istat,'(i5,17(1pe14.5))') 0,
 183      &   (energia(print_order(i)),i=1,nprint_ene),
 184      &   etot,rms,frac,co
 185       else
 186         if (print_stat) write (istat,'(i5,16(1pe14.5))') 0,
 187      &   (energia(print_order(i)),i=1,nprint_ene),etot
 188       endif
 189       if (print_stat) close(istat)
 190       neneval=neneval+nfun+1
 191       write (iout,'(/80(1h*)/20x,a/80(1h*))')
 192      &    'Enter Monte Carlo procedure.'
 193       if (print_int) then
 194         close(igeom)
 195         call briefout(0,etot)
 196       endif
 197       eold=etot
 198       do i=1,nvar
 199         varold(i)=varia(i)
 200       enddo
 201       elowest=etot
 202       call zapis(varia,etot)
 203       nacc=0         ! total # of accepted confs of the current processor.
 204       nacc_tot=0     ! total # of accepted confs of all processors.
 205
 206       not_done = (iretcode.ne.11)
 207
 208 C----------------------------------------------------------------------------
 209 C Main loop.
 210 c----------------------------------------------------------------------------
 211       it=0
 212       nout=0
 213       do while (not_done)
 214         it=it+1
 215         write (iout,'(80(1h*)/20x,a,i7)')
 216      &                             'Beginning iteration #',it
 217 C Initialize local counter.
 218         ntrial=0 ! # of generated non-overlapping confs.
 219         accepted=.false.
 220         do while (.not. accepted)
 221
 222 C Retrieve the angles of previously accepted conformation
 223           noverlap=0 ! # of overlapping confs.
 224           do j=1,nvar
 225             varia(j)=varold(j)
 226           enddo
 227           call var_to_geom(nvar,varia)
 228 C Rebuild the chain.
 229           call chainbuild
 230 C Heat up the system, if necessary.
 231           call heat(over)
 232 C If temperature cannot be further increased, stop.
 233           if (over) goto 20
 234           MoveType=0
 235           nbond=0
 236           lprint=.true.
 237 cd        write (iout,'(a)') 'Old variables:'
 238 cd        write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
 239 C Decide whether to generate a random conformation or perturb the old one
 240           RandOrPert=ran_number(0.0D0,1.0D0)
 241           if (RandOrPert.gt.RanFract) then
 242             if (print_mc.gt.0)
 243      &        write (iout,'(a)') 'Perturbation-generated conformation.'
 244             call perturb(error,lprint,MoveType,nbond,1.0D0)
 245             if (error) goto 20
 246             if (MoveType.lt.1 .or. MoveType.gt.MaxMoveType) then
 247               write (iout,'(/a,i7,a/)') 'Error - unknown MoveType=',
 248      &           MoveType,' returned from PERTURB.'
 249               goto 20
 250             endif
 251             call chainbuild
 252           else
 253             MoveType=0
 254             moves(0)=moves(0)+1
 255             nstart_grow=iran_num(3,nres)
 256             if (print_mc.gt.0)
 257      &        write (iout,'(2a,i3)') 'Random-generated conformation',
 258      &        ' - chain regrown from residue',nstart_grow
 259             call gen_rand_conf(nstart_grow,*30)
 260           endif
 261           call geom_to_var(nvar,varia)
 262 cd        write (iout,'(a)') 'New variables:'
 263 cd        write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
 264           ngen=ngen+1
 265
 266           call etotal(energia(0))
 267           etot=energia(0)
 268 c         call enerprint(energia(0))
 269 c         write (iout,'(2(a,1pe14.5))') 'Etot=',Etot,' Elowest=',Elowest
 270           if (etot-elowest.gt.overlap_cut) then
 271             if(iprint.gt.1.or.etot.lt.1d20)
 272      &       write (iout,'(a,1pe14.5)')
 273      &      'Overlap detected in the current conf.; energy is',etot
 274             neneval=neneval+1
 275             accepted=.false.
 276             noverlap=noverlap+1
 277             if (noverlap.gt.maxoverlap) then
 278               write (iout,'(a)') 'Too many overlapping confs.'
 279               goto 20
 280             endif
 281           else
 282             if (minim) then
 283               call minimize(etot,varia,iretcode,nfun)
 284 cd            write (iout,*) 'etot from MINIMIZE:',etot
 285 cd            write (iout,'(a)') 'Variables after minimization:'
 286 cd            write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
 287
 288               call etotal(energia(0))
 289               etot = energia(0)
 290               neneval=neneval+nfun+2
 291             endif
 292 c           call enerprint(energia(0))
 293             write (iout,'(a,i6,a,1pe16.6)') 'Conformation:',ngen,
 294      &      ' energy:',etot
 295 C--------------------------------------------------------------------------
 296 C... Do Metropolis test
 297 C--------------------------------------------------------------------------
 298             accepted=.false.
 299             my_conf=.false.
 300
 301             if (WhatsUp.eq.0 .and. Kwita.eq.0) then
 302               call metropolis(nvar,varia,varold,etot,eold,accepted,
 303      &                      my_conf,EneLower,it)
 304             endif
 305             write (iout,*) 'My_Conf=',My_Conf,' EneLower=',EneLower
 306             if (accepted) then
 307
 308               nacc=nacc+1
 309               nacc_tot=nacc_tot+1
 310               if (elowest.gt.etot) elowest=etot
 311               moves_acc(MoveType)=moves_acc(MoveType)+1
 312               if (MoveType.eq.1) then
 313                 nbond_acc(nbond)=nbond_acc(nbond)+1
 314               endif
 315 C Check against conformation repetitions.
 316               irepet=conf_comp(varia,etot)
 317               if (print_stat) then
 318 #if defined(AIX) || defined(PGI)
 319               open (istat,file=statname,position='append')
 320 #else
 321                open (istat,file=statname,access='append')
 322 #endif
 323               endif
 324               call statprint(nacc,nfun,iretcode,etot,elowest)
 325               if (refstr) then
 326                 call var_to_geom(nvar,varia)
 327                 call chainbuild
 328                 call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup),
 329      &                    nsup,przes,obr,non_conv)
 330                 rms=dsqrt(rms)
 331                 call contact(.false.,ncont,icont,co)
 332                 frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
 333                 write (iout,'(a,f8.3,a,f8.3)')
 334      &          'RMS deviation from the reference structure:',rms,
 335      &          ' % of native contacts:',frac*100,' contact order',co
 336               endif ! refstr
 337               if (My_Conf) then
 338                 nout=nout+1
 339                 write (iout,*) 'Writing new conformation',nout
 340                 if (refstr) then
 341                   write (istat,'(i5,16(1pe14.5))') nout,
 342      &             (energia(print_order(i)),i=1,nprint_ene),
 343      &             etot,rms,frac
 344                 else
 345                   if (print_stat)
 346      &             write (istat,'(i5,17(1pe14.5))') nout,
 347      &              (energia(print_order(i)),i=1,nprint_ene),etot
 348                 endif ! refstr
 349                 if (print_stat) close(istat)
 350 C Print internal coordinates.
 351                 if (print_int) call briefout(nout,etot)
 352 C Accumulate the newly accepted conf in the coord1 array, if it is different
 353 C from all confs that are already there.
 354                 call compare_s1(n_thr,max_thread2,etot,varia,ii,
 355      &           enetb1,coord1,rms_deform,.true.,iprint)
 356                 write (iout,*) 'After compare_ss: n_thr=',n_thr
 357                 if (ii.eq.1 .or. ii.eq.3) then
 358                   write (iout,'(8f10.4)')
 359      &                (rad2deg*coord1(i,n_thr),i=1,nvar)
 360                 endif
 361               else
 362                 write (iout,*) 'Conformation from cache, not written.'
 363               endif ! My_Conf
 364
 365               if (nrepm.gt.maxrepm) then
 366                 write (iout,'(a)') 'Too many conformation repetitions.'
 367                 goto 20
 368               endif
 369 C Store the accepted conf. and its energy.
 370               eold=etot
 371               do i=1,nvar
 372                 varold(i)=varia(i)
 373               enddo
 374               if (irepet.eq.0) call zapis(varia,etot)
 375 C Lower the temperature, if necessary.
 376               call cool
 377
 378             else
 379
 380               ntrial=ntrial+1
 381             endif ! accepted
 382           endif ! overlap
 383
 384    30     continue
 385         enddo ! accepted
 386 C Check for time limit.
 387         if (ovrtim()) WhatsUp=-1
 388         not_done = (nacc_tot.lt.maxacc) .and. (WhatsUp.eq.0)
 389      &       .and. (Kwita.eq.0)
 390
 391       enddo ! not_done
 392       goto 21
 393    20 WhatsUp=-3
 394
 395    21 continue
 396       runtime=tcpu()
 397       write (iout,'(/80(1h*)/20x,a)') 'Summary run statistics:'
 398       call statprint(nacc,nfun,iretcode,etot,elowest)
 399       write (iout,'(a)')
 400      & 'Statistics of multiple-bond motions. Total motions:'
 401       write (iout,'(16i5)') (nbond_move(i),i=1,Nbm)
 402       write (iout,'(a)') 'Accepted motions:'
 403       write (iout,'(16i5)') (nbond_acc(i),i=1,Nbm)
 404       if (it.ge.maxacc)
 405      &write (iout,'(/80(1h*)/20x,a/80(1h*)/)') 'All iterations done.'
 406
 407       return
 408       end
 409 #endif
 410 #ifdef MPI
 411 c------------------------------------------------------------------------------
 412       subroutine do_mcm(i_orig)
 413 C Monte-Carlo-with-Minimization calculations - parallel code.
 414       implicit real*8 (a-h,o-z)
 415       include 'DIMENSIONS'
 416       include 'mpif.h'
 417       include 'COMMON.IOUNITS'
 418       include 'COMMON.GEO'
 419       include 'COMMON.CHAIN'
 420       include 'COMMON.MCM'
 421       include 'COMMON.CONTACTS'
 422       include 'COMMON.CONTROL'
 423       include 'COMMON.VAR'
 424       include 'COMMON.INTERACT'
 425       include 'COMMON.INFO'
 426       include 'COMMON.CACHE'
 427 crc      include 'COMMON.DEFORM'
 428 crc      include 'COMMON.DEFORM1'
 429 crc      include 'COMMON.DEFORM2'
 430       include 'COMMON.MINIM'
 431       include 'COMMON.NAMES'
 432       logical accepted,over,ovrtim,error,lprint,not_done,similar,
 433      &        enelower,non_conv,flag,finish
 434       integer MoveType,nbond,conf_comp
 435       double precision varia(maxvar),varold(maxvar),elowest,eold,
 436      & x1(maxvar), varold1(maxvar), przes(3),obr(3,3)
 437       integer iparentx(max_threadss2)
 438       integer iparentx1(max_threadss2)
 439       integer imtasks(150),imtasks_n
 440       double precision energia(0:n_ene)
 441
 442       print *,'Master entered DO_MCM'
 443       nodenum = nprocs
 444
 445       finish=.false.
 446       imtasks_n=0
 447       do i=1,nodenum-1
 448        imtasks(i)=0
 449       enddo
 450 C---------------------------------------------------------------------------
 451 C Initialize counters.
 452 C---------------------------------------------------------------------------
 453 C Total number of generated confs.
 454       ngen=0
 455 C Total number of moves. In general this won`t be equal to the number of
 456 C attempted moves, because we may want to reject some "bad" confs just by
 457 C overlap check.
 458       nmove=0
 459 C Total number of temperature jumps.
 460       ntherm=0
 461 C Total number of shift (nbond_move(1)), spike, crankshaft, three-bond,...
 462 C motions.
 463       ncache=0
 464       do i=1,nres
 465         nbond_move(i)=0
 466       enddo
 467 C Initialize total and accepted number of moves of various kind.
 468       do i=0,MaxMoveType
 469         moves(i)=0
 470         moves_acc(i)=0
 471       enddo
 472 C Total number of energy evaluations.
 473       neneval=0
 474       nfun=0
 475       nsave=0
 476 c      write (iout,*) 'RanFract=',RanFract
 477       WhatsUp=0
 478       Kwita=0
 479 c----------------------------------------------------------------------------
 480 C Compute and print initial energies.
 481 c----------------------------------------------------------------------------
 482       call intout
 483       write (iout,'(/80(1h*)/a)') 'Initial energies:'
 484       call chainbuild
 485       nf=0
 486       call etotal(energia(0))
 487       etot = energia(0)
 488       call enerprint(energia(0))
 489 C Request energy computation from slave processors.
 490       call geom_to_var(nvar,varia)
 491 !     write (iout,*) 'The VARIA array'
 492 !     write (iout,'(8f10.4)') (rad2deg*varia(i),i=1,nvar)
 493       call minimize(etot,varia,iretcode,nfun)
 494       call var_to_geom(nvar,varia)
 495       call chainbuild
 496       write (iout,*) 'etot from MINIMIZE:',etot
 497 !     write (iout,*) 'Tha VARIA array'
 498 !     write (iout,'(8f10.4)') (rad2deg*varia(i),i=1,nvar)
 499       neneval=0
 500       eneglobal=1.0d99
 501       if (print_mc .gt. 0) write (iout,'(/80(1h*)/20x,a/80(1h*))')
 502      &    'Enter Monte Carlo procedure.'
 503       if (print_mc .gt. 0) write (iout,'(i5,1pe14.5)' ) i_orig,etot
 504       eold=etot
 505       do i=1,nvar
 506         varold(i)=varia(i)
 507       enddo
 508       elowest=etot
 509       call zapis(varia,etot)
 510 c diagnostics
 511       call var_to_geom(nvar,varia)
 512       call chainbuild
 513       call etotal(energia(0))
 514       if (print_mc.gt.0) write (iout,*) 'Initial energy:',etot
 515 c end diagnostics
 516       nacc=0         ! total # of accepted confs of the current processor.
 517       nacc_tot=0     ! total # of accepted confs of all processors.
 518       not_done=.true.
 519 C----------------------------------------------------------------------------
 520 C Main loop.
 521 c----------------------------------------------------------------------------
 522       it=0
 523       nout=0
 524       LOOP1:do while (not_done)
 525         it=it+1
 526         if (print_mc.gt.0) write (iout,'(80(1h*)/20x,a,i7)')
 527      &                             'Beginning iteration #',it
 528 C Initialize local counter.
 529         ntrial=0 ! # of generated non-overlapping confs.
 530         noverlap=0 ! # of overlapping confs.
 531         accepted=.false.
 532         LOOP2:do while (.not. accepted)
 533
 534          LOOP3:do while (imtasks_n.lt.nodenum-1.and..not.finish)
 535           do i=1,nodenum-1
 536            if(imtasks(i).eq.0) then
 537             is=i
 538             exit
 539            endif
 540           enddo
 541 C Retrieve the angles of previously accepted conformation
 542           do j=1,nvar
 543             varia(j)=varold(j)
 544           enddo
 545           call var_to_geom(nvar,varia)
 546 C Rebuild the chain.
 547           call chainbuild
 548 C Heat up the system, if necessary.
 549           call heat(over)
 550 C If temperature cannot be further increased, stop.
 551           if (over) then
 552            finish=.true.
 553           endif
 554           MoveType=0
 555           nbond=0
 556 c          write (iout,'(a)') 'Old variables:'
 557 c          write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
 558 C Decide whether to generate a random conformation or perturb the old one
 559           RandOrPert=ran_number(0.0D0,1.0D0)
 560           if (RandOrPert.gt.RanFract) then
 561            if (print_mc.gt.0)
 562      &       write (iout,'(a)') 'Perturbation-generated conformation.'
 563            call perturb(error,lprint,MoveType,nbond,1.0D0)
 564 c           print *,'after perturb',error,finish
 565            if (error) finish = .true.
 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             finish=.true.
 570             write (*,'(/a,i7,a/)') 'Error - unknown MoveType=',
 571      &         MoveType,' returned from PERTURB.'
 572            endif
 573            call chainbuild
 574           else
 575            MoveType=0
 576            moves(0)=moves(0)+1
 577            nstart_grow=iran_num(3,nres)
 578            if (print_mc.gt.0)
 579      &      write (iout,'(2a,i3)') 'Random-generated conformation',
 580      &      ' - chain regrown from residue',nstart_grow
 581            call gen_rand_conf(nstart_grow,*30)
 582           endif
 583           call geom_to_var(nvar,varia)
 584           ngen=ngen+1
 585 c          print *,'finish=',finish
 586           if (etot-elowest.gt.overlap_cut) then
 587            if (print_mc.gt.1) write (iout,'(a,1pe14.5)')
 588      &    'Overlap detected in the current conf.; energy is',etot
 589            if(iprint.gt.1.or.etot.lt.1d19) print *,
 590      &     'Overlap detected in the current conf.; energy is',etot
 591            neneval=neneval+1
 592            accepted=.false.
 593            noverlap=noverlap+1
 594            if (noverlap.gt.maxoverlap) then
 595             write (iout,*) 'Too many overlapping confs.',
 596      &      ' etot, elowest, overlap_cut', etot, elowest, overlap_cut
 597             finish=.true.
 598            endif
 599           else if (.not. finish) then
 600 C Distribute tasks to processors
 601 c           print *,'Master sending order'
 602            call MPI_SEND(12, 1, MPI_INTEGER, is, tag,
 603      &             CG_COMM, ierr)
 604 c           write (iout,*) '12: tag=',tag
 605 c           print *,'Master sent order to processor',is
 606            call MPI_SEND(it, 1, MPI_INTEGER, is, tag,
 607      &             CG_COMM, ierr)
 608 c           write (iout,*) 'it: tag=',tag
 609            call MPI_SEND(eold, 1, MPI_DOUBLE_PRECISION, is, tag,
 610      &             CG_COMM, ierr)
 611 c           write (iout,*) 'eold: tag=',tag
 612            call MPI_SEND(varia(1), nvar, MPI_DOUBLE_PRECISION,
 613      &             is, tag,
 614      &             CG_COMM, ierr)
 615 c           write (iout,*) 'varia: tag=',tag
 616            call MPI_SEND(varold(1), nvar, MPI_DOUBLE_PRECISION,
 617      &             is, tag,
 618      &             CG_COMM, ierr)
 619 c           write (iout,*) 'varold: tag=',tag
 620 #ifdef AIX
 621            call flush_(iout)
 622 #else
 623            call flush(iout)
 624 #endif
 625            imtasks(is)=1
 626            imtasks_n=imtasks_n+1
 627 C End distribution
 628           endif ! overlap
 629          enddo LOOP3
 630
 631          flag = .false.
 632          LOOP_RECV:do while(.not.flag)
 633           do is=1, nodenum-1
 634            call MPI_IPROBE(is,tag,CG_COMM,flag,status,ierr)
 635            if(flag) then
 636             call MPI_RECV(iitt, 1, MPI_INTEGER, is, tag,
 637      &              CG_COMM, status, ierr)
 638             call MPI_RECV(eold1, 1, MPI_DOUBLE_PRECISION, is, tag,
 639      &              CG_COMM, status, ierr)
 640             call MPI_RECV(etot, 1, MPI_DOUBLE_PRECISION, is, tag,
 641      &              CG_COMM, status, ierr)
 642             call MPI_RECV(varia(1), nvar, MPI_DOUBLE_PRECISION,is,tag,
 643      &              CG_COMM, status, ierr)
 644             call MPI_RECV(varold1(1), nvar, MPI_DOUBLE_PRECISION, is,
 645      &              tag, CG_COMM, status, ierr)
 646             call MPI_RECV(ii_grnum_d, 1, MPI_INTEGER, is, tag,
 647      &              CG_COMM, status, ierr)
 648             call MPI_RECV(ii_ennum_d, 1, MPI_INTEGER, is, tag,
 649      &              CG_COMM, status, ierr)
 650             call MPI_RECV(ii_hesnum_d, 1, MPI_INTEGER, is, tag,
 651      &              CG_COMM, status, ierr)
 652             i_grnum_d=i_grnum_d+ii_grnum_d
 653             i_ennum_d=i_ennum_d+ii_ennum_d
 654             neneval = neneval+ii_ennum_d
 655             i_hesnum_d=i_hesnum_d+ii_hesnum_d
 656             i_minimiz=i_minimiz+1
 657             imtasks(is)=0
 658             imtasks_n=imtasks_n-1
 659             exit
 660            endif
 661           enddo
 662          enddo LOOP_RECV
 663
 664          if(print_mc.gt.0) write (iout,'(a,i6,a,i6,a,i6,a,1pe16.6)')
 665      &      'From Worker #',is,' iitt',iitt,
 666      &     ' Conformation:',ngen,' energy:',etot
 667 C--------------------------------------------------------------------------
 668 C... Do Metropolis test
 669 C--------------------------------------------------------------------------
 670          call metropolis(nvar,varia,varold1,etot,eold1,accepted,
 671      &                      similar,EneLower)
 672          if(iitt.ne.it.and..not.similar) then
 673           call metropolis(nvar,varia,varold,etot,eold,accepted,
 674      &                      similar,EneLower)
 675           accepted=enelower
 676          endif
 677          if(etot.lt.eneglobal)eneglobal=etot
 678 c         if(mod(it,100).eq.0)
 679          write(iout,*)'CHUJOJEB ',neneval,eneglobal
 680          if (accepted) then
 681 C Write the accepted conformation.
 682            nout=nout+1
 683            if (refstr) then
 684              call var_to_geom(nvar,varia)
 685              call chainbuild
 686              kkk=1
 687              call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,kkk),
 688      &                    nsup,przes,obr,non_conv)
 689              rms=dsqrt(rms)
 690              call contact(.false.,ncont,icont,co)
 691              frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
 692              write (iout,'(a,f8.3,a,f8.3,a,f8.3)')
 693      &         'RMS deviation from the reference structure:',rms,
 694      &         ' % of native contacts:',frac*100,' contact order:',co
 695            endif ! refstr
 696            if (print_mc.gt.0)
 697      &      write (iout,*) 'Writing new conformation',nout
 698            if (print_stat) then
 699              call var_to_geom(nvar,varia)
 700 #if defined(AIX) || defined(PGI)
 701              open (istat,file=statname,position='append')
 702 #else
 703              open (istat,file=statname,access='append')
 704 #endif
 705              if (refstr) then
 706                write (istat,'(i5,16(1pe14.5))') nout,
 707      &          (energia(print_order(i)),i=1,nprint_ene),
 708      &          etot,rms,frac
 709              else
 710                write (istat,'(i5,16(1pe14.5))') nout,
 711      &          (energia(print_order(i)),i=1,nprint_ene),etot
 712              endif ! refstr
 713              close(istat)
 714            endif ! print_stat
 715 C Print internal coordinates.
 716            if (print_int) call briefout(nout,etot)
 717            nacc=nacc+1
 718            nacc_tot=nacc_tot+1
 719            if (elowest.gt.etot) elowest=etot
 720            moves_acc(MoveType)=moves_acc(MoveType)+1
 721            if (MoveType.eq.1) then
 722              nbond_acc(nbond)=nbond_acc(nbond)+1
 723            endif
 724 C Check against conformation repetitions.
 725           irepet=conf_comp(varia,etot)
 726           if (nrepm.gt.maxrepm) then
 727            if (print_mc.gt.0)
 728      &      write (iout,'(a)') 'Too many conformation repetitions.'
 729             finish=.true.
 730            endif
 731 C Store the accepted conf. and its energy.
 732            eold=etot
 733            do i=1,nvar
 734             varold(i)=varia(i)
 735            enddo
 736            if (irepet.eq.0) call zapis(varia,etot)
 737 C Lower the temperature, if necessary.
 738            call cool
 739           else
 740            ntrial=ntrial+1
 741          endif ! accepted
 742    30    continue
 743          if(finish.and.imtasks_n.eq.0)exit LOOP2
 744         enddo LOOP2 ! accepted
 745 C Check for time limit.
 746         not_done = (it.lt.max_mcm_it) .and. (nacc_tot.lt.maxacc)
 747         if(.not.not_done .or. finish) then
 748          if(imtasks_n.gt.0) then
 749           not_done=.true.
 750          else
 751           not_done=.false.
 752          endif
 753          finish=.true.
 754         endif
 755       enddo LOOP1 ! not_done
 756       runtime=tcpu()
 757       if (print_mc.gt.0) then
 758         write (iout,'(/80(1h*)/20x,a)') 'Summary run statistics:'
 759         call statprint(nacc,nfun,iretcode,etot,elowest)
 760         write (iout,'(a)')
 761      & 'Statistics of multiple-bond motions. Total motions:'
 762         write (iout,'(16i5)') (nbond_move(i),i=1,Nbm)
 763         write (iout,'(a)') 'Accepted motions:'
 764         write (iout,'(16i5)') (nbond_acc(i),i=1,Nbm)
 765         if (it.ge.maxacc)
 766      &write (iout,'(/80(1h*)/20x,a/80(1h*)/)') 'All iterations done.'
 767       endif
 768 #ifdef AIX
 769       call flush_(iout)
 770 #else
 771       call flush(iout)
 772 #endif
 773       do is=1,nodenum-1
 774         call MPI_SEND(999, 1, MPI_INTEGER, is, tag,
 775      &             CG_COMM, ierr)
 776       enddo
 777       return
 778       end
 779 c------------------------------------------------------------------------------
 780       subroutine execute_slave(nodeinfo,iprint)
 781       implicit real*8 (a-h,o-z)
 782       include 'DIMENSIONS'
 783       include 'mpif.h'
 784       include 'COMMON.TIME1'
 785       include 'COMMON.IOUNITS'
 786 crc      include 'COMMON.DEFORM'
 787 crc      include 'COMMON.DEFORM1'
 788 crc      include 'COMMON.DEFORM2'
 789       include 'COMMON.LOCAL'
 790       include 'COMMON.VAR'
 791       include 'COMMON.INFO'
 792       include 'COMMON.MINIM'
 793       character*10 nodeinfo
 794       double precision x(maxvar),x1(maxvar)
 795       nodeinfo='chujwdupe'
 796 c      print *,'Processor:',MyID,' Entering execute_slave'
 797       tag=0
 798 c      call MPI_SEND(nodeinfo, 10, MPI_CHARACTER, 0, tag,
 799 c     &              CG_COMM, ierr)
 800
 801 1001  call MPI_RECV(i_switch, 1, MPI_INTEGER, 0, tag,
 802      &              CG_COMM, status, ierr)
 803 c      write(iout,*)'12: tag=',tag
 804       if(iprint.ge.2)print *, MyID,' recv order ',i_switch
 805       if (i_switch.eq.12) then
 806        i_grnum_d=0
 807        i_ennum_d=0
 808        i_hesnum_d=0
 809        call MPI_RECV(iitt, 1, MPI_INTEGER, 0, tag,
 810      &               CG_COMM, status, ierr)
 811 c      write(iout,*)'12: tag=',tag
 812        call MPI_RECV(ener, 1, MPI_DOUBLE_PRECISION, 0, tag,
 813      &               CG_COMM, status, ierr)
 814 c      write(iout,*)'ener: tag=',tag
 815        call MPI_RECV(x(1), nvar, MPI_DOUBLE_PRECISION, 0, tag,
 816      &               CG_COMM, status, ierr)
 817 c      write(iout,*)'x: tag=',tag
 818        call MPI_RECV(x1(1), nvar, MPI_DOUBLE_PRECISION, 0, tag,
 819      &               CG_COMM, status, ierr)
 820 c      write(iout,*)'x1: tag=',tag
 821 #ifdef AIX
 822        call flush_(iout)
 823 #else
 824        call flush(iout)
 825 #endif
 826 c       print *,'calling minimize'
 827        call minimize(energyx,x,iretcode,nfun)
 828        if(iprint.gt.0)
 829      &  write(iout,100)'minimized energy = ',energyx,
 830      &    ' # funeval:',nfun,' iret ',iretcode
 831         write(*,100)'minimized energy = ',energyx,
 832      &    ' # funeval:',nfun,' iret ',iretcode
 833  100   format(a20,f10.5,a12,i5,a6,i2)
 834        if(iretcode.eq.10) then
 835          do iminrep=2,3
 836           if(iprint.gt.1)
 837      &    write(iout,*)' ... not converged - trying again ',iminrep
 838           call minimize(energyx,x,iretcode,nfun)
 839           if(iprint.gt.1)
 840      &    write(iout,*)'minimized energy = ',energyx,
 841      &     ' # funeval:',nfun,' iret ',iretcode
 842           if(iretcode.ne.10)go to 812
 843          enddo
 844          if(iretcode.eq.10) then
 845           if(iprint.gt.1)
 846      &    write(iout,*)' ... not converged again - giving up'
 847           go to 812
 848          endif
 849        endif
 850 812    continue
 851 c       print *,'Sending results'
 852        call MPI_SEND(iitt, 1, MPI_INTEGER, 0, tag,
 853      &              CG_COMM, ierr)
 854        call MPI_SEND(ener, 1, MPI_DOUBLE_PRECISION, 0, tag,
 855      &              CG_COMM, ierr)
 856        call MPI_SEND(energyx, 1, MPI_DOUBLE_PRECISION, 0, tag,
 857      &              CG_COMM, ierr)
 858        call MPI_SEND(x(1), nvar, MPI_DOUBLE_PRECISION, 0, tag,
 859      &              CG_COMM, ierr)
 860        call MPI_SEND(x1(1), nvar, MPI_DOUBLE_PRECISION, 0, tag,
 861      &              CG_COMM, ierr)
 862        call MPI_SEND(i_grnum_d, 1, MPI_INTEGER, 0, tag,
 863      &              CG_COMM, ierr)
 864        call MPI_SEND(nfun, 1, MPI_INTEGER, 0, tag,
 865      &              CG_COMM, ierr)
 866        call MPI_SEND(i_hesnum_d, 1, MPI_INTEGER, 0, tag,
 867      &              CG_COMM, ierr)
 868 c       print *,'End sending'
 869        go to 1001
 870       endif
 871
 872       return
 873       end
 874 #endif
 875 c------------------------------------------------------------------------------
 876       subroutine statprint(it,nfun,iretcode,etot,elowest)
 877       implicit real*8 (a-h,o-z)
 878       include 'DIMENSIONS'
 879       include 'COMMON.IOUNITS'
 880       include 'COMMON.CONTROL'
 881       include 'COMMON.MCM'
 882       if (minim) then
 883         write (iout,
 884      &  '(80(1h*)/a,i5,a,1pe14.5,a,1pe14.5/a,i3,a,i10,a,i5,a,i5)')
 885      &  'Finished iteration #',it,' energy is',etot,
 886      &  ' lowest energy:',elowest,
 887      &  'SUMSL return code:',iretcode,
 888      &  ' # of energy evaluations:',neneval,
 889      &  '# of temperature jumps:',ntherm,
 890      &  ' # of minima repetitions:',nrepm
 891       else
 892         write (iout,'(80(1h*)/a,i8,a,1pe14.5,a,1pe14.5)')
 893      &  'Finished iteration #',it,' energy is',etot,
 894      &  ' lowest energy:',elowest
 895       endif
 896       write (iout,'(/4a)')
 897      & 'Kind of move   ','           total','       accepted',
 898      & '  fraction'
 899       write (iout,'(58(1h-))')
 900       do i=-1,MaxMoveType
 901         if (moves(i).eq.0) then
 902           fr_mov_i=0.0d0
 903         else
 904           fr_mov_i=dfloat(moves_acc(i))/dfloat(moves(i))
 905         endif
 906         write(iout,'(a,2i15,f10.5)')MovTypID(i),moves(i),moves_acc(i),
 907      &         fr_mov_i
 908       enddo
 909       write (iout,'(a,2i15,f10.5)') 'total           ',nmove,nacc_tot,
 910      &         dfloat(nacc_tot)/dfloat(nmove)
 911       write (iout,'(58(1h-))')
 912       write (iout,'(a,1pe12.4)') 'Elapsed time:',tcpu()
 913       return
 914       end
 915 c------------------------------------------------------------------------------
 916       subroutine heat(over)
 917       implicit real*8 (a-h,o-z)
 918       include 'DIMENSIONS'
 919       include 'COMMON.MCM'
 920       include 'COMMON.IOUNITS'
 921       logical over
 922 C Check if there`s a need to increase temperature.
 923       if (ntrial.gt.maxtrial) then
 924         if (NstepH.gt.0) then
 925           if (dabs(Tcur-TMax).lt.1.0D-7) then
 926             if (print_mc.gt.0)
 927      &      write (iout,'(/80(1h*)/a,f8.3,a/80(1h*))')
 928      &      'Upper limit of temperature reached. Terminating.'
 929             over=.true.
 930             Tcur=Tmin
 931           else
 932             Tcur=Tcur*TstepH
 933             if (Tcur.gt.Tmax) Tcur=Tmax
 934             betbol=1.0D0/(Rbol*Tcur)
 935             if (print_mc.gt.0)
 936      &      write (iout,'(/80(1h*)/a,f8.3,a,f10.5/80(1h*))')
 937      &      'System heated up to ',Tcur,' K; BetBol:',betbol
 938             ntherm=ntherm+1
 939             ntrial=0
 940             over=.false.
 941           endif
 942         else
 943          if (print_mc.gt.0)
 944      &   write (iout,'(a)')
 945      & 'Maximum number of trials in a single MCM iteration exceeded.'
 946          over=.true.
 947          Tcur=Tmin
 948         endif
 949       else
 950         over=.false.
 951       endif
 952       return
 953       end
 954 c------------------------------------------------------------------------------
 955       subroutine cool
 956       implicit real*8 (a-h,o-z)
 957       include 'DIMENSIONS'
 958       include 'COMMON.MCM'
 959       include 'COMMON.IOUNITS'
 960       if (nstepC.gt.0 .and. dabs(Tcur-Tmin).gt.1.0D-7) then
 961         Tcur=Tcur/TstepC
 962         if (Tcur.lt.Tmin) Tcur=Tmin
 963         betbol=1.0D0/(Rbol*Tcur)
 964         if (print_mc.gt.0)
 965      &  write (iout,'(/80(1h*)/a,f8.3,a,f10.5/80(1h*))')
 966      &  'System cooled down up to ',Tcur,' K; BetBol:',betbol
 967       endif
 968       return
 969       end
 970 C---------------------------------------------------------------------------
 971       subroutine zapis(varia,etot)
 972       implicit real*8 (a-h,o-z)
 973       include 'DIMENSIONS'
 974 #ifdef MP
 975       include 'mpif.h'
 976       include 'COMMON.INFO'
 977 #endif
 978       include 'COMMON.GEO'
 979       include 'COMMON.VAR'
 980       include 'COMMON.MCM'
 981       include 'COMMON.IOUNITS'
 982       integer itemp(maxsave)
 983       double precision varia(maxvar)
 984       logical lprint
 985       lprint=.false.
 986       if (lprint) then
 987       write (iout,'(a,i5,a,i5)') 'Enter ZAPIS NSave=',Nsave,
 988      &  ' MaxSave=',MaxSave
 989       write (iout,'(a)') 'Current energy and conformation:'
 990       write (iout,'(1pe14.5)') etot
 991       write (iout,'(10f8.3)') (rad2deg*varia(i),i=1,nvar)
 992       endif
 993 C Shift the contents of the esave and varsave arrays if filled up.
 994       call add2cache(maxvar,maxsave,nsave,nvar,MyID,itemp,
 995      &               etot,varia,esave,varsave)
 996       if (lprint) then
 997       write (iout,'(a)') 'Energies and the VarSave array.'
 998       do i=1,nsave
 999         write (iout,'(i5,1pe14.5)') i,esave(i)
1000         write (iout,'(10f8.3)') (rad2deg*varsave(j,i),j=1,nvar)
1001       enddo
1002       endif
1003       return
1004       end
1005 C---------------------------------------------------------------------------
1006       subroutine perturb(error,lprint,MoveType,nbond,max_phi)
1007       implicit real*8 (a-h,o-z)
1008       include 'DIMENSIONS'
1009       parameter (MMaxSideMove=100)
1010       include 'COMMON.MCM'
1011       include 'COMMON.CHAIN'
1012       include 'COMMON.INTERACT'
1013       include 'COMMON.VAR'
1014       include 'COMMON.GEO'
1015       include 'COMMON.NAMES'
1016       include 'COMMON.IOUNITS'
1017 crc      include 'COMMON.DEFORM1'
1018       logical error,lprint,fail
1019       integer MoveType,nbond,end_select,ind_side(MMaxSideMove)
1020       double precision max_phi
1021       double precision psi,gen_psi
1022       external iran_num
1023       integer iran_num
1024       integer ifour
1025       data ifour /4/
1026       error=.false.
1027       lprint=.false.
1028 C Perturb the conformation according to a randomly selected move.
1029       call SelectMove(MoveType)
1030 c      write (iout,*) 'MoveType=',MoveType
1031       itrial=0
1032       goto (100,200,300,400,500) MoveType
1033 C------------------------------------------------------------------------------
1034 C Backbone N-bond move.
1035 C Select the number of bonds (length of the segment to perturb).
1036   100 continue
1037       if (itrial.gt.1000) then
1038         write (iout,'(a)') 'Too many attempts at multiple-bond move.'
1039         error=.true.
1040         return
1041       endif
1042       bond_prob=ran_number(0.0D0,sumpro_bond(nbm))
1043 c     print *,'sumpro_bond(nbm)=',sumpro_bond(nbm),
1044 c    & ' Bond_prob=',Bond_Prob
1045       do i=1,nbm-1
1046 c       print *,i,Bond_Prob,sumpro_bond(i),sumpro_bond(i+1)
1047         if (bond_prob.ge.sumpro_bond(i) .and.
1048      &               bond_prob.le.sumpro_bond(i+1)) then
1049           nbond=i+1
1050           goto 10
1051         endif
1052       enddo
1053       write (iout,'(2a)') 'In PERTURB: Error - number of bonds',
1054      &                    ' to move out of range.'
1055       error=.true.
1056       return
1057    10 continue
1058       if (nwindow.gt.0) then
1059 C Select the first residue to perturb
1060         iwindow=iran_num(1,nwindow)
1061         print *,'iwindow=',iwindow
1062         iiwin=1
1063         do while (winlen(iwindow).lt.nbond)
1064           iwindow=iran_num(1,nwindow)
1065           iiwin=iiwin+1
1066           if (iiwin.gt.1000) then
1067              write (iout,'(a)') 'Cannot select moveable residues.'
1068              error=.true.
1069              return
1070           endif
1071         enddo
1072         nstart=iran_num(winstart(iwindow),winend(iwindow))
1073       else
1074         nstart = iran_num(koniecl+2,nres-nbond-koniecl)
1075 cd      print *,'nres=',nres,' nbond=',nbond,' koniecl=',koniecl,
1076 cd   &        ' nstart=',nstart
1077       endif
1078       psi = gen_psi()
1079       if (psi.eq.0.0) then
1080         error=.true.
1081         return
1082       endif
1083       if (print_mc.gt.1) write (iout,'(a,i4,a,i4,a,f8.3)')
1084      & 'PERTURB: nbond=',nbond,' nstart=',nstart,' psi=',psi*rad2deg
1085 cd    print *,'nstart=',nstart
1086       call bond_move(nbond,nstart,psi,.false.,error)
1087       if (error) then
1088         write (iout,'(2a)')
1089      & 'Could not define reference system in bond_move, ',
1090      & 'choosing ahother segment.'
1091         itrial=itrial+1
1092         goto 100
1093       endif
1094       nbond_move(nbond)=nbond_move(nbond)+1
1095       moves(1)=moves(1)+1
1096       nmove=nmove+1
1097       return
1098 C------------------------------------------------------------------------------
1099 C Backbone endmove. Perturb a SINGLE angle of a residue close to the end of
1100 C the chain.
1101   200 continue
1102       lprint=.true.
1103 c     end_select=iran_num(1,2*koniecl)
1104 c     if (end_select.gt.koniecl) then
1105 c       end_select=nphi-(end_select-koniecl)
1106 c     else
1107 c       end_select=koniecl+3
1108 c     endif
1109 c     if (nwindow.gt.0) then
1110 c       iwin=iran_num(1,nwindow)
1111 c       i1=max0(4,winstart(iwin))
1112 c       i2=min0(winend(imin)+2,nres)
1113 c       end_select=iran_num(i1,i2)
1114 c     else
1115 c      iselect = iran_num(1,nmov_var)
1116 c      jj = 0
1117 c      do i=1,nphi
1118 c        if (isearch_tab(i).eq.1) jj = jj+1
1119 c        if (jj.eq.iselect) then
1120 c          end_select=i+3
1121 c          exit
1122 c        endif
1123 c      enddo
1124 c     endif
1125       end_select = iran_num(4,nres)
1126       psi=max_phi*gen_psi()
1127       if (psi.eq.0.0D0) then
1128         error=.true.
1129         return
1130       endif
1131       phi(end_select)=pinorm(phi(end_select)+psi)
1132       if (print_mc.gt.1) write (iout,'(a,i4,a,f8.3,a,f8.3)')
1133      & 'End angle',end_select,' moved by ',psi*rad2deg,' new angle:',
1134      & phi(end_select)*rad2deg
1135 c     if (end_select.gt.3)
1136 c    &   theta(end_select-1)=gen_theta(itype(end_select-2),
1137 c    &                              phi(end_select-1),phi(end_select))
1138 c     if (end_select.lt.nres)
1139 c    &    theta(end_select)=gen_theta(itype(end_select-1),
1140 c    &                              phi(end_select),phi(end_select+1))
1141 cd    print *,'nres=',nres,' end_select=',end_select
1142 cd    print *,'theta',end_select-1,theta(end_select-1)
1143 cd    print *,'theta',end_select,theta(end_select)
1144       moves(2)=moves(2)+1
1145       nmove=nmove+1
1146       lprint=.false.
1147       return
1148 C------------------------------------------------------------------------------
1149 C Side chain move.
1150 C Select the number of SCs to perturb.
1151   300 isctry=0
1152   301 nside_move=iran_num(1,MaxSideMove)
1153 c     print *,'nside_move=',nside_move,' MaxSideMove',MaxSideMove
1154 C Select the indices.
1155       do i=1,nside_move
1156         icount=0
1157   111   inds=iran_num(nnt,nct)
1158         icount=icount+1
1159         if (icount.gt.1000) then
1160           write (iout,'(a)')'Error - cannot select side chains to move.'
1161           error=.true.
1162           return
1163         endif
1164         if (itype(inds).eq.10) goto 111
1165         do j=1,i-1
1166           if (inds.eq.ind_side(j)) goto 111
1167         enddo
1168         do j=1,i-1
1169           if (inds.lt.ind_side(j)) then
1170             indx=j
1171             goto 112
1172           endif
1173         enddo
1174         indx=i
1175   112   do j=i,indx+1,-1
1176           ind_side(j)=ind_side(j-1)
1177         enddo
1178   113   ind_side(indx)=inds
1179       enddo
1180 C Carry out perturbation.
1181       do i=1,nside_move
1182         ii=ind_side(i)
1183         iti=itype(ii)
1184         call gen_side(iti,theta(ii+1),alph(ii),omeg(ii),fail)
1185         if (fail) then
1186           isctry=isctry+1
1187           if (isctry.gt.1000) then
1188             write (iout,'(a)') 'Too many errors in SC generation.'
1189             error=.true.
1190             return
1191           endif
1192           goto 301
1193         endif
1194         if (print_mc.gt.1) write (iout,'(2a,i4,a,2f8.3)')
1195      &   'Side chain ',restyp(iti),ii,' moved to ',
1196      &   alph(ii)*rad2deg,omeg(ii)*rad2deg
1197       enddo
1198       moves(3)=moves(3)+1
1199       nmove=nmove+1
1200       return
1201 C------------------------------------------------------------------------------
1202 C THETA move
1203   400 end_select=iran_num(3,nres)
1204       theta_new=gen_theta(itype(end_select),phi(end_select),
1205      &                    phi(end_select+1))
1206       if (print_mc.gt.1) write (iout,'(a,i3,a,f8.3,a,f8.3)')
1207      & 'Theta ',end_select,' moved from',theta(end_select)*rad2deg,
1208      & ' to ',theta_new*rad2deg
1209       theta(end_select)=theta_new
1210       moves(4)=moves(4)+1
1211       nmove=nmove+1
1212       return
1213 C------------------------------------------------------------------------------
1214 C Error returned from SelectMove.
1215   500 error=.true.
1216       return
1217       end
1218 C------------------------------------------------------------------------------
1219       subroutine SelectMove(MoveType)
1220       implicit real*8 (a-h,o-z)
1221       include 'DIMENSIONS'
1222       include 'COMMON.MCM'
1223       include 'COMMON.IOUNITS'
1224       what_move=ran_number(0.0D0,sumpro_type(MaxMoveType))
1225       do i=1,MaxMoveType
1226         if (what_move.ge.sumpro_type(i-1).and.
1227      &            what_move.lt.sumpro_type(i)) then
1228           MoveType=i
1229           return
1230         endif
1231       enddo
1232       write (iout,'(a)')
1233      & 'Fatal error in SelectMoveType: cannot select move.'
1234       MoveType=MaxMoveType+1
1235       return
1236       end
1237 c----------------------------------------------------------------------------
1238       double precision function gen_psi()
1239       implicit none
1240       integer i
1241       double precision x,ran_number
1242       include 'COMMON.IOUNITS'
1243       include 'COMMON.GEO'
1244       x=0.0D0
1245       do i=1,100
1246         x=ran_number(-pi,pi)
1247         if (dabs(x).gt.angmin) then
1248           gen_psi=x
1249           return
1250         endif
1251       enddo
1252       write (iout,'(a)')'From Gen_Psi: Cannot generate angle increment.'
1253       gen_psi=0.0D0
1254       return
1255       end
1256 c----------------------------------------------------------------------------
1257       subroutine metropolis(n,xcur,xold,ecur,eold,accepted,similar,
1258      &                      enelower)
1259       implicit real*8 (a-h,o-z)
1260       include 'DIMENSIONS'
1261       include 'COMMON.MCM'
1262       include 'COMMON.IOUNITS'
1263       include 'COMMON.VAR'
1264       include 'COMMON.GEO'
1265 crc      include 'COMMON.DEFORM'
1266       double precision ecur,eold,xx,ran_number,bol
1267       double precision xcur(n),xold(n)
1268       double precision ecut1 ,ecut2 ,tola
1269       logical accepted,similar,not_done,enelower
1270       logical lprn
1271       data ecut1 /-1.0D-5/,ecut2 /5.0D-3/,tola/5.0D0/
1272 !      ecut1=-5*enedif
1273 !      ecut2=50*enedif
1274 !      tola=5.0d0
1275 C Set lprn=.true. for debugging.
1276       lprn=.false.
1277       if (lprn)
1278      &write(iout,*)'enedif',enedif,' ecut1',ecut1,' ecut2',ecut2
1279       similar=.false.
1280       enelower=.false.
1281       accepted=.false.
1282 C Check if the conformation is similar.
1283       difene=ecur-eold
1284       reldife=difene/dmax1(dabs(eold),dabs(ecur),1.0D0)
1285       if (lprn) then
1286         write (iout,*) 'Metropolis'
1287         write(iout,*)'ecur,eold,difene,reldife',ecur,eold,difene,reldife
1288       endif
1289 C If energy went down remarkably, we accept the new conformation
1290 C unconditionally.
1291 cjp      if (reldife.lt.ecut1) then
1292       if (difene.lt.ecut1) then
1293         accepted=.true.
1294         EneLower=.true.
1295         if (lprn) write (iout,'(a)')
1296      &   'Conformation accepted, because energy has lowered remarkably.'
1297 !      elseif (reldife.lt.ecut2 .and. dif_ang(nphi,xcur,xold).lt.tola)
1298 cjp      elseif (reldife.lt.ecut2)
1299       elseif (difene.lt.ecut2)
1300      & then
1301 C Reject the conf. if energy has changed insignificantly and there is not
1302 C much change in conformation.
1303         if (lprn)
1304      &   write (iout,'(2a)') 'Conformation rejected, because it is',
1305      &      ' similar to the preceding one.'
1306         accepted=.false.
1307         similar=.true.
1308       else
1309 C Else carry out Metropolis test.
1310         EneLower=.false.
1311         xx=ran_number(0.0D0,1.0D0)
1312         xxh=betbol*difene
1313         if (lprn)
1314      &    write (iout,*) 'betbol=',betbol,' difene=',difene,' xxh=',xxh
1315         if (xxh.gt.50.0D0) then
1316           bol=0.0D0
1317         else
1318           bol=exp(-xxh)
1319         endif
1320         if (lprn) write (iout,*) 'bol=',bol,' xx=',xx
1321         if (bol.gt.xx) then
1322           accepted=.true.
1323           if (lprn) write (iout,'(a)')
1324      &    'Conformation accepted, because it passed Metropolis test.'
1325         else
1326           accepted=.false.
1327           if (lprn) write (iout,'(a)')
1328      & 'Conformation rejected, because it did not pass Metropolis test.'
1329         endif
1330       endif
1331 #ifdef AIX
1332       call flush_(iout)
1333 #else
1334       call flush(iout)
1335 #endif
1336       return
1337       end
1338 c------------------------------------------------------------------------------
1339       integer function conf_comp(x,ene)
1340       implicit real*8 (a-h,o-z)
1341       include 'DIMENSIONS'
1342       include 'COMMON.MCM'
1343       include 'COMMON.VAR'
1344       include 'COMMON.IOUNITS'
1345       include 'COMMON.GEO'
1346       double precision etol , angtol
1347       double precision x(maxvar)
1348       double precision dif_ang,difa
1349       data etol /0.1D0/, angtol /20.0D0/
1350       do ii=nsave,1,-1
1351 c       write (iout,*) 'ii=',ii,'ene=',ene,esave(ii),dabs(ene-esave(ii))
1352         if (dabs(ene-esave(ii)).lt.etol) then
1353           difa=dif_ang(nphi,x,varsave(1,ii))
1354 c         do i=1,nphi
1355 c           write(iout,'(i3,3f8.3)')i,rad2deg*x(i),
1356 c    &          rad2deg*varsave(i,ii)
1357 c         enddo
1358 c         write(iout,*) 'ii=',ii,' difa=',difa,' angtol=',angtol
1359           if (difa.le.angtol) then
1360             if (print_mc.gt.0) then
1361             write (iout,'(a,i5,2(a,1pe15.4))')
1362      &      'Current conformation matches #',ii,
1363      &      ' in the store array ene=',ene,' esave=',esave(ii)
1364 c           write (*,'(a,i5,a)') 'Current conformation matches #',ii,
1365 c    &      ' in the store array.'
1366             endif ! print_mc.gt.0
1367             if (print_mc.gt.1) then
1368             do i=1,nphi
1369               write(iout,'(i3,3f8.3)')i,rad2deg*x(i),
1370      &            rad2deg*varsave(i,ii)
1371             enddo
1372             endif ! print_mc.gt.1
1373             nrepm=nrepm+1
1374             conf_comp=ii
1375             return
1376           endif
1377         endif
1378       enddo
1379       conf_comp=0
1380       return
1381       end
1382 C----------------------------------------------------------------------------
1383       double precision function dif_ang(n,x,y)
1384       implicit none
1385       integer i,n
1386       double precision x(n),y(n)
1387       double precision w,wa,dif,difa
1388       double precision pinorm
1389       include 'COMMON.GEO'
1390       wa=0.0D0
1391       difa=0.0D0
1392       do i=1,n
1393         dif=dabs(pinorm(y(i)-x(i)))
1394         if (dabs(dif-dwapi).lt.dif) dif=dabs(dif-dwapi)
1395         w=1.0D0-(2.0D0*(i-1)/(n-1)-1.0D0)**2+1.0D0/n
1396         wa=wa+w
1397         difa=difa+dif*dif*w
1398       enddo
1399       dif_ang=rad2deg*dsqrt(difa/wa)
1400       return
1401       end
1402 c--------------------------------------------------------------------------
1403       subroutine add2cache(n1,n2,ncache,nvar,SourceID,CachSrc,
1404      &                     ecur,xcur,ecache,xcache)
1405       implicit none
1406       include 'COMMON.GEO'
1407       include 'COMMON.IOUNITS'
1408       integer n1,n2,ncache,nvar,SourceID,CachSrc(n2)
1409       integer i,ii,j
1410       double precision ecur,xcur(nvar),ecache(n2),xcache(n1,n2)
1411 cd    write (iout,*) 'Enter ADD2CACHE ncache=',ncache ,' ecur',ecur
1412 cd    write (iout,'(10f8.3)') (rad2deg*xcur(i),i=1,nvar)
1413 cd    write (iout,*) 'Old CACHE array:'
1414 cd    do i=1,ncache
1415 cd      write (iout,*) 'i=',i,' ecache=',ecache(i),' CachSrc',CachSrc(i)
1416 cd      write (iout,'(10f8.3)') (rad2deg*xcache(j,i),j=1,nvar)
1417 cd    enddo
1418
1419       i=ncache
1420       do while (i.gt.0 .and. ecur.lt.ecache(i))
1421         i=i-1
1422       enddo
1423       i=i+1
1424 cd    write (iout,*) 'i=',i,' ncache=',ncache
1425       if (ncache.eq.n2) then
1426         write (iout,*) 'Cache dimension exceeded',ncache,n2
1427         write (iout,*) 'Highest-energy conformation will be removed.'
1428         ncache=ncache-1
1429       endif
1430       do ii=ncache,i,-1
1431         ecache(ii+1)=ecache(ii)
1432         CachSrc(ii+1)=CachSrc(ii)
1433         do j=1,nvar
1434           xcache(j,ii+1)=xcache(j,ii)
1435         enddo
1436       enddo
1437       ecache(i)=ecur
1438       CachSrc(i)=SourceID
1439       do j=1,nvar
1440         xcache(j,i)=xcur(j)
1441       enddo
1442       ncache=ncache+1
1443 cd    write (iout,*) 'New CACHE array:'
1444 cd    do i=1,ncache
1445 cd      write (iout,*) 'i=',i,' ecache=',ecache(i),' CachSrc',CachSrc(i)
1446 cd      write (iout,'(10f8.3)') (rad2deg*xcache(j,i),j=1,nvar)
1447 cd    enddo
1448       return
1449       end
1450 c--------------------------------------------------------------------------
1451       subroutine rm_from_cache(i,n1,n2,ncache,nvar,CachSrc,ecache,
1452      &                         xcache)
1453       implicit none
1454       include 'COMMON.GEO'
1455       include 'COMMON.IOUNITS'
1456       integer n1,n2,ncache,nvar,CachSrc(n2)
1457       integer i,ii,j
1458       double precision ecache(n2),xcache(n1,n2)
1459
1460 cd    write (iout,*) 'Enter RM_FROM_CACHE'
1461 cd    write (iout,*) 'Old CACHE array:'
1462 cd    do ii=1,ncache
1463 cd    write (iout,*)'i=',ii,' ecache=',ecache(ii),' CachSrc',CachSrc(ii)
1464 cd      write (iout,'(10f8.3)') (rad2deg*xcache(j,ii),j=1,nvar)
1465 cd    enddo
1466
1467       do ii=i+1,ncache
1468         ecache(ii-1)=ecache(ii)
1469         CachSrc(ii-1)=CachSrc(ii)
1470         do j=1,nvar
1471           xcache(j,ii-1)=xcache(j,ii)
1472         enddo
1473       enddo
1474       ncache=ncache-1
1475 cd    write (iout,*) 'New CACHE array:'
1476 cd    do i=1,ncache
1477 cd      write (iout,*) 'i=',i,' ecache=',ecache(i),' CachSrc',CachSrc(i)
1478 cd      write (iout,'(10f8.3)') (rad2deg*xcache(j,i),j=1,nvar)
1479 cd    enddo
1480       return
1481       end