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