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