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