source/unres/src_MD-M-SAXS.safe/mcm.F

   1       subroutine mcm_setup
   2       implicit real*8 (a-h,o-z)
   3       include 'DIMENSIONS'
   4       include 'COMMON.IOUNITS'
   5       include 'COMMON.MCM'
   6       include 'COMMON.CONTROL'
   7       include 'COMMON.INTERACT'
   8       include 'COMMON.NAMES'
   9       include 'COMMON.CHAIN'
  10       include 'COMMON.VAR'
  11 C
  12 C Set up variables used in MC/MCM.
  13 C
  14       write (iout,'(80(1h*)/20x,a/80(1h*))') 'MCM control parameters:'
  15       write (iout,'(5(a,i7))') 'Maxacc:',maxacc,' MaxTrial:',MaxTrial,
  16      & ' MaxRepm:',MaxRepm,' MaxGen:',MaxGen,' MaxOverlap:',MaxOverlap
  17       write (iout,'(4(a,f8.1)/2(a,i3))')
  18      & 'Tmin:',Tmin,' Tmax:',Tmax,' TstepH:',TstepH,
  19      & ' TstepC:',TstepC,'NstepH:',NstepH,' NstepC:',NstepC
  20       if (nwindow.gt.0) then
  21         write (iout,'(a)') 'Perturbation windows:'
  22         do i=1,nwindow
  23           i1=winstart(i)
  24           i2=winend(i)
  25           it1=itype(i1)
  26           it2=itype(i2)
  27           write (iout,'(a,i3,a,i3,a,i3)') restyp(it1),i1,restyp(it2),i2,
  28      &                        ' length',winlen(i)
  29         enddo
  30       endif
  31 C Rbolt=8.3143D-3*2.388459D-01 kcal/(mol*K)
  32       RBol=1.9858D-3
  33 C Number of "end bonds".
  34       koniecl=0
  35 c     koniecl=nphi
  36       print *,'koniecl=',koniecl
  37       write (iout,'(a)') 'Probabilities of move types:'
  38       write (*,'(a)') 'Probabilities of move types:'
  39       do i=1,MaxMoveType
  40         write (iout,'(a,f10.5)') MovTypID(i),
  41      &    sumpro_type(i)-sumpro_type(i-1)
  42         write (*,'(a,f10.5)') MovTypID(i),
  43      &    sumpro_type(i)-sumpro_type(i-1)
  44       enddo
  45       write (iout,*)
  46 C Maximum length of N-bond segment to be moved
  47 c     nbm=nres-1-(2*koniecl-1)
  48       if (nwindow.gt.0) then
  49         maxwinlen=winlen(1)
  50         do i=2,nwindow
  51           if (winlen(i).gt.maxwinlen) maxwinlen=winlen(i)
  52         enddo
  53         nbm=min0(maxwinlen,6)
  54         write (iout,'(a,i3,a,i3)') 'Nbm=',Nbm,' Maxwinlen=',Maxwinlen
  55       else
  56         nbm=min0(6,nres-2)
  57       endif
  58       sumpro_bond(0)=0.0D0
  59       sumpro_bond(1)=0.0D0
  60       do i=2,nbm
  61         sumpro_bond(i)=sumpro_bond(i-1)+1.0D0/dfloat(i)
  62       enddo
  63       write (iout,'(a)') 'The SumPro_Bond array:'
  64       write (iout,'(8f10.5)') (sumpro_bond(i),i=1,nbm)
  65       write (*,'(a)') 'The SumPro_Bond array:'
  66       write (*,'(8f10.5)') (sumpro_bond(i),i=1,nbm)
  67 C Maximum number of side chains moved simultaneously
  68 c     print *,'nnt=',nnt,' nct=',nct
  69       ngly=0
  70       do i=nnt,nct
  71         if (itype(i).eq.10) ngly=ngly+1
  72       enddo
  73       mmm=nct-nnt-ngly+1
  74       if (mmm.gt.0) then
  75         MaxSideMove=min0((nct-nnt+1)/2,mmm)
  76       endif
  77 c     print *,'MaxSideMove=',MaxSideMove
  78 C Max. number of generated confs (not used at present).
  79       maxgen=10000
  80 C Set initial temperature
  81       Tcur=Tmin
  82       betbol=1.0D0/(Rbol*Tcur)
  83       write (iout,'(a,f8.1,a,f10.5)') 'Initial temperature:',Tcur,
  84      &    ' BetBol:',betbol
  85       write (iout,*) 'RanFract=',ranfract
  86       return
  87       end
  88 c------------------------------------------------------------------------------
  89 #ifndef MPI
  90       subroutine do_mcm(i_orig)
  91 C Monte-Carlo-with-Minimization calculations - serial code.
  92       implicit real*8 (a-h,o-z)
  93       include 'DIMENSIONS'
  94       include 'COMMON.IOUNITS'
  95       include 'COMMON.GEO'
  96       include 'COMMON.CHAIN'
  97       include 'COMMON.MCM'
  98       include 'COMMON.CONTACTS'
  99       include 'COMMON.CONTROL'
 100       include 'COMMON.VAR'
 101       include 'COMMON.INTERACT'
 102       include 'COMMON.CACHE'
 103 crc      include 'COMMON.DEFORM'
 104 crc      include 'COMMON.DEFORM1'
 105       include 'COMMON.NAMES'
 106       logical accepted,over,ovrtim,error,lprint,not_done,my_conf,
 107      &        enelower,non_conv
 108       integer MoveType,nbond,conf_comp
 109       integer ifeed(max_cache)
 110       double precision varia(maxvar),varold(maxvar),elowest,eold,
 111      & przes(3),obr(3,3)
 112       double precision energia(0:n_ene)
 113       double precision coord1(maxres,3)
 114
 115
 116 C---------------------------------------------------------------------------
 117 C Initialize counters.
 118 C---------------------------------------------------------------------------
 119 C Total number of generated confs.
 120       ngen=0
 121 C Total number of moves. In general this won't be equal to the number of
 122 C attempted moves, because we may want to reject some "bad" confs just by
 123 C overlap check.
 124       nmove=0
 125 C Total number of temperature jumps.
 126       ntherm=0
 127 C Total number of shift (nbond_move(1)), spike, crankshaft, three-bond,...
 128 C motions.
 129       ncache=0
 130       do i=1,nres
 131         nbond_move(i)=0
 132       enddo
 133 C Initialize total and accepted number of moves of various kind.
 134       do i=0,MaxMoveType
 135         moves(i)=0
 136         moves_acc(i)=0
 137       enddo
 138 C Total number of energy evaluations.
 139       neneval=0
 140       nfun=0
 141       nsave=0
 142
 143       write (iout,*) 'RanFract=',RanFract
 144
 145       WhatsUp=0
 146       Kwita=0
 147
 148 c----------------------------------------------------------------------------
 149 C Compute and print initial energies.
 150 c----------------------------------------------------------------------------
 151       call intout
 152       write (iout,'(/80(1h*)/a)') 'Initial energies:'
 153       call chainbuild
 154       nf=0
 155
 156       call etotal(energia(0))
 157       etot = energia(0)
 158 C Minimize the energy of the first conformation.
 159       if (minim) then
 160         call geom_to_var(nvar,varia)
 161 !       write (iout,*) 'The VARIA array'
 162 !       write (iout,'(8f10.4)') (rad2deg*varia(i),i=1,nvar)
 163         call minimize(etot,varia,iretcode,nfun)
 164         call var_to_geom(nvar,varia)
 165         call chainbuild
 166         write (iout,*) 'etot from MINIMIZE:',etot
 167 !       write (iout,*) 'Tha VARIA array'
 168 !       write (iout,'(8f10.4)') (rad2deg*varia(i),i=1,nvar)
 169
 170         call etotal(energia(0))
 171         etot=energia(0)
 172         call enerprint(energia(0))
 173       endif
 174       if (refstr) then
 175         call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,1),nsup,przes,
 176      &             obr,non_conv)
 177         rms=dsqrt(rms)
 178         call contact(.false.,ncont,icont,co)
 179         frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
 180         write (iout,'(a,f8.3,a,f8.3,a,f8.3)')
 181      &    'RMS deviation from the reference structure:',rms,
 182      &    ' % of native contacts:',frac*100,' contact order:',co
 183         if (print_stat)
 184      &  write (istat,'(i5,17(1pe14.5))') 0,
 185      &   (energia(print_order(i)),i=1,nprint_ene),
 186      &   etot,rms,frac,co
 187       else
 188         if (print_stat) write (istat,'(i5,16(1pe14.5))') 0,
 189      &   (energia(print_order(i)),i=1,nprint_ene),etot
 190       endif
 191       if (print_stat) close(istat)
 192       neneval=neneval+nfun+1
 193       write (iout,'(/80(1h*)/20x,a/80(1h*))')
 194      &    'Enter Monte Carlo procedure.'
 195       if (print_int) then
 196         close(igeom)
 197         call briefout(0,etot)
 198       endif
 199       eold=etot
 200       do i=1,nvar
 201         varold(i)=varia(i)
 202       enddo
 203       elowest=etot
 204       call zapis(varia,etot)
 205       nacc=0         ! total # of accepted confs of the current processor.
 206       nacc_tot=0     ! total # of accepted confs of all processors.
 207
 208       not_done = (iretcode.ne.11)
 209
 210 C----------------------------------------------------------------------------
 211 C Main loop.
 212 c----------------------------------------------------------------------------
 213       it=0
 214       nout=0
 215       do while (not_done)
 216         it=it+1
 217         write (iout,'(80(1h*)/20x,a,i7)')
 218      &                             'Beginning iteration #',it
 219 C Initialize local counter.
 220         ntrial=0 ! # of generated non-overlapping confs.
 221         accepted=.false.
 222         do while (.not. accepted)
 223
 224 C Retrieve the angles of previously accepted conformation
 225           noverlap=0 ! # of overlapping confs.
 226           do j=1,nvar
 227             varia(j)=varold(j)
 228           enddo
 229           call var_to_geom(nvar,varia)
 230 C Rebuild the chain.
 231           call chainbuild
 232 C Heat up the system, if necessary.
 233           call heat(over)
 234 C If temperature cannot be further increased, stop.
 235           if (over) goto 20
 236           MoveType=0
 237           nbond=0
 238           lprint=.true.
 239 cd        write (iout,'(a)') 'Old variables:'
 240 cd        write (iout,'(10f8.1)') (rad2deg*varia(i),i=1,nvar)
 241 C Decide whether to generate a random conformation or perturb the old one
 242           RandOrPert=ran_number(0.0D0,1.0D0)
 243           if (RandOrPert.gt.RanFract) then
 244             if (print_mc.gt.0)
 245      &        write (iout,'(a)') 'Perturbation-generated conformation.'
 246             call perturb(error,lprint,MoveType,nbond,1.0D0)
 247             if (error) goto 20
 248             if (MoveType.lt.1 .or. MoveType.gt.MaxMoveType) then
 249               write (iout,'(/a,i7,a/)') 'Error - unknown MoveType=',
 250      &           MoveType,' returned from PERTURB.'
 251               goto 20
 252             endif
 253             call chainbuild
 254           else
 255             MoveType=0
 256             moves(0)=moves(0)+1
 257             nstart_grow=iran_num(3,nres)
 258             if (print_mc.gt.0)
 259      &        write (iout,'(2a,i3)') 'Random-generated conformation',
 260      &        ' - chain regrown from residue',nstart_grow
 261             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,1),
 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              kkk=1
 689              call fitsq(rms,c(1,nstart_seq),cref(1,nstart_sup,kkk),
 690      &                    nsup,przes,obr,non_conv)
 691              rms=dsqrt(rms)
 692              call contact(.false.,ncont,icont,co)
 693              frac=contact_fract(ncont,ncont_ref,icont,icont_ref)
 694              write (iout,'(a,f8.3,a,f8.3,a,f8.3)')
 695      &         'RMS deviation from the reference structure:',rms,
 696      &         ' % of native contacts:',frac*100,' contact order:',co
 697            endif ! refstr
 698            if (print_mc.gt.0)
 699      &      write (iout,*) 'Writing new conformation',nout
 700            if (print_stat) then
 701              call var_to_geom(nvar,varia)
 702 #if defined(AIX) || defined(PGI) || defined(CRAY)
 703              open (istat,file=statname,position='append')
 704 #else
 705              open (istat,file=statname,access='append')
 706 #endif
 707              if (refstr) then
 708                write (istat,'(i5,16(1pe14.5))') nout,
 709      &          (energia(print_order(i)),i=1,nprint_ene),
 710      &          etot,rms,frac
 711              else
 712                write (istat,'(i5,16(1pe14.5))') nout,
 713      &          (energia(print_order(i)),i=1,nprint_ene),etot
 714              endif ! refstr
 715              close(istat)
 716            endif ! print_stat
 717 C Print internal coordinates.
 718            if (print_int) call briefout(nout,etot)
 719            nacc=nacc+1
 720            nacc_tot=nacc_tot+1
 721            if (elowest.gt.etot) elowest=etot
 722            moves_acc(MoveType)=moves_acc(MoveType)+1
 723            if (MoveType.eq.1) then
 724              nbond_acc(nbond)=nbond_acc(nbond)+1
 725            endif
 726 C Check against conformation repetitions.
 727           irepet=conf_comp(varia,etot)
 728           if (nrepm.gt.maxrepm) then
 729            if (print_mc.gt.0)
 730      &      write (iout,'(a)') 'Too many conformation repetitions.'
 731             finish=.true.
 732            endif
 733 C Store the accepted conf. and its energy.
 734            eold=etot
 735            do i=1,nvar
 736             varold(i)=varia(i)
 737            enddo
 738            if (irepet.eq.0) call zapis(varia,etot)
 739 C Lower the temperature, if necessary.
 740            call cool
 741           else
 742            ntrial=ntrial+1
 743          endif ! accepted
 744    30    continue
 745          if(finish.and.imtasks_n.eq.0)exit LOOP2
 746         enddo LOOP2 ! accepted
 747 C Check for time limit.
 748         not_done = (it.lt.max_mcm_it) .and. (nacc_tot.lt.maxacc)
 749         if(.not.not_done .or. finish) then
 750          if(imtasks_n.gt.0) then
 751           not_done=.true.
 752          else
 753           not_done=.false.
 754          endif
 755          finish=.true.
 756         endif
 757       enddo LOOP1 ! not_done
 758       runtime=tcpu()
 759       if (print_mc.gt.0) then
 760         write (iout,'(/80(1h*)/20x,a)') 'Summary run statistics:'
 761         call statprint(nacc,nfun,iretcode,etot,elowest)
 762         write (iout,'(a)')
 763      & 'Statistics of multiple-bond motions. Total motions:'
 764         write (iout,'(16i5)') (nbond_move(i),i=1,Nbm)
 765         write (iout,'(a)') 'Accepted motions:'
 766         write (iout,'(16i5)') (nbond_acc(i),i=1,Nbm)
 767         if (it.ge.maxacc)
 768      &write (iout,'(/80(1h*)/20x,a/80(1h*)/)') 'All iterations done.'
 769       endif
 770 #ifdef AIX
 771       call flush_(iout)
 772 #else
 773       call flush(iout)
 774 #endif
 775       do is=1,nodenum-1
 776         call MPI_SEND(999, 1, MPI_INTEGER, is, tag,
 777      &             CG_COMM, ierr)
 778       enddo
 779       return
 780       end
 781 c------------------------------------------------------------------------------
 782       subroutine execute_slave(nodeinfo,iprint)
 783       implicit real*8 (a-h,o-z)
 784       include 'DIMENSIONS'
 785       include 'mpif.h'
 786       include 'COMMON.TIME1'
 787       include 'COMMON.IOUNITS'
 788 crc      include 'COMMON.DEFORM'
 789 crc      include 'COMMON.DEFORM1'
 790 crc      include 'COMMON.DEFORM2'
 791       include 'COMMON.LOCAL'
 792       include 'COMMON.VAR'
 793       include 'COMMON.INFO'
 794       include 'COMMON.MINIM'
 795       character*10 nodeinfo
 796       double precision x(maxvar),x1(maxvar)
 797       nodeinfo='chujwdupe'
 798 c      print *,'Processor:',MyID,' Entering execute_slave'
 799       tag=0
 800 c      call MPI_SEND(nodeinfo, 10, MPI_CHARACTER, 0, tag,
 801 c     &              CG_COMM, ierr)
 802
 803 1001  call MPI_RECV(i_switch, 1, MPI_INTEGER, 0, tag,
 804      &              CG_COMM, status, ierr)
 805 c      write(iout,*)'12: tag=',tag
 806       if(iprint.ge.2)print *, MyID,' recv order ',i_switch
 807       if (i_switch.eq.12) then
 808        i_grnum_d=0
 809        i_ennum_d=0
 810        i_hesnum_d=0
 811        call MPI_RECV(iitt, 1, MPI_INTEGER, 0, tag,
 812      &               CG_COMM, status, ierr)
 813 c      write(iout,*)'12: tag=',tag
 814        call MPI_RECV(ener, 1, MPI_DOUBLE_PRECISION, 0, tag,
 815      &               CG_COMM, status, ierr)
 816 c      write(iout,*)'ener: tag=',tag
 817        call MPI_RECV(x(1), nvar, MPI_DOUBLE_PRECISION, 0, tag,
 818      &               CG_COMM, status, ierr)
 819 c      write(iout,*)'x: tag=',tag
 820        call MPI_RECV(x1(1), nvar, MPI_DOUBLE_PRECISION, 0, tag,
 821      &               CG_COMM, status, ierr)
 822 c      write(iout,*)'x1: tag=',tag
 823 #ifdef AIX
 824        call flush_(iout)
 825 #else
 826        call flush(iout)
 827 #endif
 828 c       print *,'calling minimize'
 829        call minimize(energyx,x,iretcode,nfun)
 830        if(iprint.gt.0)
 831      &  write(iout,100)'minimized energy = ',energyx,
 832      &    ' # funeval:',nfun,' iret ',iretcode
 833         write(*,100)'minimized energy = ',energyx,
 834      &    ' # funeval:',nfun,' iret ',iretcode
 835  100   format(a20,f10.5,a12,i5,a6,i2)
 836        if(iretcode.eq.10) then
 837          do iminrep=2,3
 838           if(iprint.gt.1)
 839      &    write(iout,*)' ... not converged - trying again ',iminrep
 840           call minimize(energyx,x,iretcode,nfun)
 841           if(iprint.gt.1)
 842      &    write(iout,*)'minimized energy = ',energyx,
 843      &     ' # funeval:',nfun,' iret ',iretcode
 844           if(iretcode.ne.10)go to 812
 845          enddo
 846          if(iretcode.eq.10) then
 847           if(iprint.gt.1)
 848      &    write(iout,*)' ... not converged again - giving up'
 849           go to 812
 850          endif
 851        endif
 852 812    continue
 853 c       print *,'Sending results'
 854        call MPI_SEND(iitt, 1, MPI_INTEGER, 0, tag,
 855      &              CG_COMM, ierr)
 856        call MPI_SEND(ener, 1, MPI_DOUBLE_PRECISION, 0, tag,
 857      &              CG_COMM, ierr)
 858        call MPI_SEND(energyx, 1, MPI_DOUBLE_PRECISION, 0, tag,
 859      &              CG_COMM, ierr)
 860        call MPI_SEND(x(1), nvar, MPI_DOUBLE_PRECISION, 0, tag,
 861      &              CG_COMM, ierr)
 862        call MPI_SEND(x1(1), nvar, MPI_DOUBLE_PRECISION, 0, tag,
 863      &              CG_COMM, ierr)
 864        call MPI_SEND(i_grnum_d, 1, MPI_INTEGER, 0, tag,
 865      &              CG_COMM, ierr)
 866        call MPI_SEND(nfun, 1, MPI_INTEGER, 0, tag,
 867      &              CG_COMM, ierr)
 868        call MPI_SEND(i_hesnum_d, 1, MPI_INTEGER, 0, tag,
 869      &              CG_COMM, ierr)
 870 c       print *,'End sending'
 871        go to 1001
 872       endif
 873
 874       return
 875       end
 876 #endif
 877 c------------------------------------------------------------------------------
 878       subroutine statprint(it,nfun,iretcode,etot,elowest)
 879       implicit real*8 (a-h,o-z)
 880       include 'DIMENSIONS'
 881       include 'COMMON.IOUNITS'
 882       include 'COMMON.CONTROL'
 883       include 'COMMON.MCM'
 884       if (minim) then
 885         write (iout,
 886      &  '(80(1h*)/a,i5,a,1pe14.5,a,1pe14.5/a,i3,a,i10,a,i5,a,i5)')
 887      &  'Finished iteration #',it,' energy is',etot,
 888      &  ' lowest energy:',elowest,
 889      &  'SUMSL return code:',iretcode,
 890      &  ' # of energy evaluations:',neneval,
 891      &  '# of temperature jumps:',ntherm,
 892      &  ' # of minima repetitions:',nrepm
 893       else
 894         write (iout,'(80(1h*)/a,i8,a,1pe14.5,a,1pe14.5)')
 895      &  'Finished iteration #',it,' energy is',etot,
 896      &  ' lowest energy:',elowest
 897       endif
 898       write (iout,'(/4a)')
 899      & 'Kind of move   ','           total','       accepted',
 900      & '  fraction'
 901       write (iout,'(58(1h-))')
 902       do i=-1,MaxMoveType
 903         if (moves(i).eq.0) then
 904           fr_mov_i=0.0d0
 905         else
 906           fr_mov_i=dfloat(moves_acc(i))/dfloat(moves(i))
 907         endif
 908         write(iout,'(a,2i15,f10.5)')MovTypID(i),moves(i),moves_acc(i),
 909      &         fr_mov_i
 910       enddo
 911       write (iout,'(a,2i15,f10.5)') 'total           ',nmove,nacc_tot,
 912      &         dfloat(nacc_tot)/dfloat(nmove)
 913       write (iout,'(58(1h-))')
 914       write (iout,'(a,1pe12.4)') 'Elapsed time:',tcpu()
 915       return
 916       end
 917 c------------------------------------------------------------------------------
 918       subroutine heat(over)
 919       implicit real*8 (a-h,o-z)
 920       include 'DIMENSIONS'
 921       include 'COMMON.MCM'
 922       include 'COMMON.IOUNITS'
 923       logical over
 924 C Check if there`s a need to increase temperature.
 925       if (ntrial.gt.maxtrial) then
 926         if (NstepH.gt.0) then
 927           if (dabs(Tcur-TMax).lt.1.0D-7) then
 928             if (print_mc.gt.0)
 929      &      write (iout,'(/80(1h*)/a,f8.3,a/80(1h*))')
 930      &      'Upper limit of temperature reached. Terminating.'
 931             over=.true.
 932             Tcur=Tmin
 933           else
 934             Tcur=Tcur*TstepH
 935             if (Tcur.gt.Tmax) Tcur=Tmax
 936             betbol=1.0D0/(Rbol*Tcur)
 937             if (print_mc.gt.0)
 938      &      write (iout,'(/80(1h*)/a,f8.3,a,f10.5/80(1h*))')
 939      &      'System heated up to ',Tcur,' K; BetBol:',betbol
 940             ntherm=ntherm+1
 941             ntrial=0
 942             over=.false.
 943           endif
 944         else
 945          if (print_mc.gt.0)
 946      &   write (iout,'(a)')
 947      & 'Maximum number of trials in a single MCM iteration exceeded.'
 948          over=.true.
 949          Tcur=Tmin
 950         endif
 951       else
 952         over=.false.
 953       endif
 954       return
 955       end
 956 c------------------------------------------------------------------------------
 957       subroutine cool
 958       implicit real*8 (a-h,o-z)
 959       include 'DIMENSIONS'
 960       include 'COMMON.MCM'
 961       include 'COMMON.IOUNITS'
 962       if (nstepC.gt.0 .and. dabs(Tcur-Tmin).gt.1.0D-7) then
 963         Tcur=Tcur/TstepC
 964         if (Tcur.lt.Tmin) Tcur=Tmin
 965         betbol=1.0D0/(Rbol*Tcur)
 966         if (print_mc.gt.0)
 967      &  write (iout,'(/80(1h*)/a,f8.3,a,f10.5/80(1h*))')
 968      &  'System cooled down up to ',Tcur,' K; BetBol:',betbol
 969       endif
 970       return
 971       end
 972 C---------------------------------------------------------------------------
 973       subroutine zapis(varia,etot)
 974       implicit real*8 (a-h,o-z)
 975       include 'DIMENSIONS'
 976 #ifdef MP
 977       include 'mpif.h'
 978       include 'COMMON.INFO'
 979 #endif
 980       include 'COMMON.GEO'
 981       include 'COMMON.VAR'
 982       include 'COMMON.MCM'
 983       include 'COMMON.IOUNITS'
 984       integer itemp(maxsave)
 985       double precision varia(maxvar)
 986       logical lprint
 987       lprint=.false.
 988       if (lprint) then
 989       write (iout,'(a,i5,a,i5)') 'Enter ZAPIS NSave=',Nsave,
 990      &  ' MaxSave=',MaxSave
 991       write (iout,'(a)') 'Current energy and conformation:'
 992       write (iout,'(1pe14.5)') etot
 993       write (iout,'(10f8.3)') (rad2deg*varia(i),i=1,nvar)
 994       endif
 995 C Shift the contents of the esave and varsave arrays if filled up.
 996       call add2cache(maxvar,maxsave,nsave,nvar,MyID,itemp,
 997      &               etot,varia,esave,varsave)
 998       if (lprint) then
 999       write (iout,'(a)') 'Energies and the VarSave array.'
1000       do i=1,nsave
1001         write (iout,'(i5,1pe14.5)') i,esave(i)
1002         write (iout,'(10f8.3)') (rad2deg*varsave(j,i),j=1,nvar)
1003       enddo
1004       endif
1005       return
1006       end
1007 C---------------------------------------------------------------------------
1008       subroutine perturb(error,lprint,MoveType,nbond,max_phi)
1009       implicit real*8 (a-h,o-z)
1010       include 'DIMENSIONS'
1011       parameter (MMaxSideMove=100)
1012       include 'COMMON.MCM'
1013       include 'COMMON.CHAIN'
1014       include 'COMMON.INTERACT'
1015       include 'COMMON.VAR'
1016       include 'COMMON.GEO'
1017       include 'COMMON.NAMES'
1018       include 'COMMON.IOUNITS'
1019 crc      include 'COMMON.DEFORM1'
1020       logical error,lprint,fail
1021       integer MoveType,nbond,end_select,ind_side(MMaxSideMove)
1022       double precision max_phi
1023       double precision psi,gen_psi
1024       external iran_num
1025       integer iran_num
1026       integer ifour
1027       data ifour /4/
1028       error=.false.
1029       lprint=.false.
1030 C Perturb the conformation according to a randomly selected move.
1031       call SelectMove(MoveType)
1032 c      write (iout,*) 'MoveType=',MoveType
1033       itrial=0
1034       goto (100,200,300,400,500) MoveType
1035 C------------------------------------------------------------------------------
1036 C Backbone N-bond move.
1037 C Select the number of bonds (length of the segment to perturb).
1038   100 continue
1039       if (itrial.gt.1000) then
1040         write (iout,'(a)') 'Too many attempts at multiple-bond move.'
1041         error=.true.
1042         return
1043       endif
1044       bond_prob=ran_number(0.0D0,sumpro_bond(nbm))
1045 c     print *,'sumpro_bond(nbm)=',sumpro_bond(nbm),
1046 c    & ' Bond_prob=',Bond_Prob
1047       do i=1,nbm-1
1048 c       print *,i,Bond_Prob,sumpro_bond(i),sumpro_bond(i+1)
1049         if (bond_prob.ge.sumpro_bond(i) .and.
1050      &               bond_prob.le.sumpro_bond(i+1)) then
1051           nbond=i+1
1052           goto 10
1053         endif
1054       enddo
1055       write (iout,'(2a)') 'In PERTURB: Error - number of bonds',
1056      &                    ' to move out of range.'
1057       error=.true.
1058       return
1059    10 continue
1060       if (nwindow.gt.0) then
1061 C Select the first residue to perturb
1062         iwindow=iran_num(1,nwindow)
1063         print *,'iwindow=',iwindow
1064         iiwin=1
1065         do while (winlen(iwindow).lt.nbond)
1066           iwindow=iran_num(1,nwindow)
1067           iiwin=iiwin+1
1068           if (iiwin.gt.1000) then
1069              write (iout,'(a)') 'Cannot select moveable residues.'
1070              error=.true.
1071              return
1072           endif
1073         enddo
1074         nstart=iran_num(winstart(iwindow),winend(iwindow))
1075       else
1076         nstart = iran_num(koniecl+2,nres-nbond-koniecl)
1077 cd      print *,'nres=',nres,' nbond=',nbond,' koniecl=',koniecl,
1078 cd   &        ' nstart=',nstart
1079       endif
1080       psi = gen_psi()
1081       if (psi.eq.0.0) then
1082         error=.true.
1083         return
1084       endif
1085       if (print_mc.gt.1) write (iout,'(a,i4,a,i4,a,f8.3)')
1086      & 'PERTURB: nbond=',nbond,' nstart=',nstart,' psi=',psi*rad2deg
1087 cd    print *,'nstart=',nstart
1088       call bond_move(nbond,nstart,psi,.false.,error)
1089       if (error) then
1090         write (iout,'(2a)')
1091      & 'Could not define reference system in bond_move, ',
1092      & 'choosing ahother segment.'
1093         itrial=itrial+1
1094         goto 100
1095       endif
1096       nbond_move(nbond)=nbond_move(nbond)+1
1097       moves(1)=moves(1)+1
1098       nmove=nmove+1
1099       return
1100 C------------------------------------------------------------------------------
1101 C Backbone endmove. Perturb a SINGLE angle of a residue close to the end of
1102 C the chain.
1103   200 continue
1104       lprint=.true.
1105 c     end_select=iran_num(1,2*koniecl)
1106 c     if (end_select.gt.koniecl) then
1107 c       end_select=nphi-(end_select-koniecl)
1108 c     else
1109 c       end_select=koniecl+3
1110 c     endif
1111 c     if (nwindow.gt.0) then
1112 c       iwin=iran_num(1,nwindow)
1113 c       i1=max0(4,winstart(iwin))
1114 c       i2=min0(winend(imin)+2,nres)
1115 c       end_select=iran_num(i1,i2)
1116 c     else
1117 c      iselect = iran_num(1,nmov_var)
1118 c      jj = 0
1119 c      do i=1,nphi
1120 c        if (isearch_tab(i).eq.1) jj = jj+1
1121 c        if (jj.eq.iselect) then
1122 c          end_select=i+3
1123 c          exit
1124 c        endif
1125 c      enddo
1126 c     endif
1127       end_select = iran_num(4,nres)
1128       psi=max_phi*gen_psi()
1129       if (psi.eq.0.0D0) then
1130         error=.true.
1131         return
1132       endif
1133       phi(end_select)=pinorm(phi(end_select)+psi)
1134       if (print_mc.gt.1) write (iout,'(a,i4,a,f8.3,a,f8.3)')
1135      & 'End angle',end_select,' moved by ',psi*rad2deg,' new angle:',
1136      & phi(end_select)*rad2deg
1137 c     if (end_select.gt.3)
1138 c    &   theta(end_select-1)=gen_theta(itype(end_select-2),
1139 c    &                              phi(end_select-1),phi(end_select))
1140 c     if (end_select.lt.nres)
1141 c    &    theta(end_select)=gen_theta(itype(end_select-1),
1142 c    &                              phi(end_select),phi(end_select+1))
1143 cd    print *,'nres=',nres,' end_select=',end_select
1144 cd    print *,'theta',end_select-1,theta(end_select-1)
1145 cd    print *,'theta',end_select,theta(end_select)
1146       moves(2)=moves(2)+1
1147       nmove=nmove+1
1148       lprint=.false.
1149       return
1150 C------------------------------------------------------------------------------
1151 C Side chain move.
1152 C Select the number of SCs to perturb.
1153   300 isctry=0
1154   301 nside_move=iran_num(1,MaxSideMove)
1155 c     print *,'nside_move=',nside_move,' MaxSideMove',MaxSideMove
1156 C Select the indices.
1157       do i=1,nside_move
1158         icount=0
1159   111   inds=iran_num(nnt,nct)
1160         icount=icount+1
1161         if (icount.gt.1000) then
1162           write (iout,'(a)')'Error - cannot select side chains to move.'
1163           error=.true.
1164           return
1165         endif
1166         if (itype(inds).eq.10) goto 111
1167         do j=1,i-1
1168           if (inds.eq.ind_side(j)) goto 111
1169         enddo
1170         do j=1,i-1
1171           if (inds.lt.ind_side(j)) then
1172             indx=j
1173             goto 112
1174           endif
1175         enddo
1176         indx=i
1177   112   do j=i,indx+1,-1
1178           ind_side(j)=ind_side(j-1)
1179         enddo
1180   113   ind_side(indx)=inds
1181       enddo
1182 C Carry out perturbation.
1183       do i=1,nside_move
1184         ii=ind_side(i)
1185         iti=itype(ii)
1186         call gen_side(iti,theta(ii+1),alph(ii),omeg(ii),fail)
1187         if (fail) then
1188           isctry=isctry+1
1189           if (isctry.gt.1000) then
1190             write (iout,'(a)') 'Too many errors in SC generation.'
1191             error=.true.
1192             return
1193           endif
1194           goto 301
1195         endif
1196         if (print_mc.gt.1) write (iout,'(2a,i4,a,2f8.3)')
1197      &   'Side chain ',restyp(iti),ii,' moved to ',
1198      &   alph(ii)*rad2deg,omeg(ii)*rad2deg
1199       enddo
1200       moves(3)=moves(3)+1
1201       nmove=nmove+1
1202       return
1203 C------------------------------------------------------------------------------
1204 C THETA move
1205   400 end_select=iran_num(3,nres)
1206       theta_new=gen_theta(itype(end_select),phi(end_select),
1207      &                    phi(end_select+1))
1208       if (print_mc.gt.1) write (iout,'(a,i3,a,f8.3,a,f8.3)')
1209      & 'Theta ',end_select,' moved from',theta(end_select)*rad2deg,
1210      & ' to ',theta_new*rad2deg
1211       theta(end_select)=theta_new
1212       moves(4)=moves(4)+1
1213       nmove=nmove+1
1214       return
1215 C------------------------------------------------------------------------------
1216 C Error returned from SelectMove.
1217   500 error=.true.
1218       return
1219       end
1220 C------------------------------------------------------------------------------
1221       subroutine SelectMove(MoveType)
1222       implicit real*8 (a-h,o-z)
1223       include 'DIMENSIONS'
1224       include 'COMMON.MCM'
1225       include 'COMMON.IOUNITS'
1226       what_move=ran_number(0.0D0,sumpro_type(MaxMoveType))
1227       do i=1,MaxMoveType
1228         if (what_move.ge.sumpro_type(i-1).and.
1229      &            what_move.lt.sumpro_type(i)) then
1230           MoveType=i
1231           return
1232         endif
1233       enddo
1234       write (iout,'(a)')
1235      & 'Fatal error in SelectMoveType: cannot select move.'
1236       MoveType=MaxMoveType+1
1237       return
1238       end
1239 c----------------------------------------------------------------------------
1240       double precision function gen_psi()
1241       implicit none
1242       integer i
1243       double precision x,ran_number
1244       include 'COMMON.IOUNITS'
1245       include 'COMMON.GEO'
1246       x=0.0D0
1247       do i=1,100
1248         x=ran_number(-pi,pi)
1249         if (dabs(x).gt.angmin) then
1250           gen_psi=x
1251           return
1252         endif
1253       enddo
1254       write (iout,'(a)')'From Gen_Psi: Cannot generate angle increment.'
1255       gen_psi=0.0D0
1256       return
1257       end
1258 c----------------------------------------------------------------------------
1259       subroutine metropolis(n,xcur,xold,ecur,eold,accepted,similar,
1260      &                      enelower)
1261       implicit real*8 (a-h,o-z)
1262       include 'DIMENSIONS'
1263       include 'COMMON.MCM'
1264       include 'COMMON.IOUNITS'
1265       include 'COMMON.VAR'
1266       include 'COMMON.GEO'
1267 crc      include 'COMMON.DEFORM'
1268       double precision ecur,eold,xx,ran_number,bol
1269       double precision xcur(n),xold(n)
1270       double precision ecut1 ,ecut2 ,tola
1271       logical accepted,similar,not_done,enelower
1272       logical lprn
1273       data ecut1 /-1.0D-5/,ecut2 /5.0D-3/,tola/5.0D0/
1274 !      ecut1=-5*enedif
1275 !      ecut2=50*enedif
1276 !      tola=5.0d0
1277 C Set lprn=.true. for debugging.
1278       lprn=.false.
1279       if (lprn)
1280      &write(iout,*)'enedif',enedif,' ecut1',ecut1,' ecut2',ecut2
1281       similar=.false.
1282       enelower=.false.
1283       accepted=.false.
1284 C Check if the conformation is similar.
1285       difene=ecur-eold
1286       reldife=difene/dmax1(dabs(eold),dabs(ecur),1.0D0)
1287       if (lprn) then
1288         write (iout,*) 'Metropolis'
1289         write(iout,*)'ecur,eold,difene,reldife',ecur,eold,difene,reldife
1290       endif
1291 C If energy went down remarkably, we accept the new conformation
1292 C unconditionally.
1293 cjp      if (reldife.lt.ecut1) then
1294       if (difene.lt.ecut1) then
1295         accepted=.true.
1296         EneLower=.true.
1297         if (lprn) write (iout,'(a)')
1298      &   'Conformation accepted, because energy has lowered remarkably.'
1299 !      elseif (reldife.lt.ecut2 .and. dif_ang(nphi,xcur,xold).lt.tola)
1300 cjp      elseif (reldife.lt.ecut2)
1301       elseif (difene.lt.ecut2)
1302      & then
1303 C Reject the conf. if energy has changed insignificantly and there is not
1304 C much change in conformation.
1305         if (lprn)
1306      &   write (iout,'(2a)') 'Conformation rejected, because it is',
1307      &      ' similar to the preceding one.'
1308         accepted=.false.
1309         similar=.true.
1310       else
1311 C Else carry out Metropolis test.
1312         EneLower=.false.
1313         xx=ran_number(0.0D0,1.0D0)
1314         xxh=betbol*difene
1315         if (lprn)
1316      &    write (iout,*) 'betbol=',betbol,' difene=',difene,' xxh=',xxh
1317         if (xxh.gt.50.0D0) then
1318           bol=0.0D0
1319         else
1320           bol=exp(-xxh)
1321         endif
1322         if (lprn) write (iout,*) 'bol=',bol,' xx=',xx
1323         if (bol.gt.xx) then
1324           accepted=.true.
1325           if (lprn) write (iout,'(a)')
1326      &    'Conformation accepted, because it passed Metropolis test.'
1327         else
1328           accepted=.false.
1329           if (lprn) write (iout,'(a)')
1330      & 'Conformation rejected, because it did not pass Metropolis test.'
1331         endif
1332       endif
1333 #ifdef AIX
1334       call flush_(iout)
1335 #else
1336       call flush(iout)
1337 #endif
1338       return
1339       end
1340 c------------------------------------------------------------------------------
1341       integer function conf_comp(x,ene)
1342       implicit real*8 (a-h,o-z)
1343       include 'DIMENSIONS'
1344       include 'COMMON.MCM'
1345       include 'COMMON.VAR'
1346       include 'COMMON.IOUNITS'
1347       include 'COMMON.GEO'
1348       double precision etol , angtol
1349       double precision x(maxvar)
1350       double precision dif_ang,difa
1351       data etol /0.1D0/, angtol /20.0D0/
1352       do ii=nsave,1,-1
1353 c       write (iout,*) 'ii=',ii,'ene=',ene,esave(ii),dabs(ene-esave(ii))
1354         if (dabs(ene-esave(ii)).lt.etol) then
1355           difa=dif_ang(nphi,x,varsave(1,ii))
1356 c         do i=1,nphi
1357 c           write(iout,'(i3,3f8.3)')i,rad2deg*x(i),
1358 c    &          rad2deg*varsave(i,ii)
1359 c         enddo
1360 c         write(iout,*) 'ii=',ii,' difa=',difa,' angtol=',angtol
1361           if (difa.le.angtol) then
1362             if (print_mc.gt.0) then
1363             write (iout,'(a,i5,2(a,1pe15.4))')
1364      &      'Current conformation matches #',ii,
1365      &      ' in the store array ene=',ene,' esave=',esave(ii)
1366 c           write (*,'(a,i5,a)') 'Current conformation matches #',ii,
1367 c    &      ' in the store array.'
1368             endif ! print_mc.gt.0
1369             if (print_mc.gt.1) then
1370             do i=1,nphi
1371               write(iout,'(i3,3f8.3)')i,rad2deg*x(i),
1372      &            rad2deg*varsave(i,ii)
1373             enddo
1374             endif ! print_mc.gt.1
1375             nrepm=nrepm+1
1376             conf_comp=ii
1377             return
1378           endif
1379         endif
1380       enddo
1381       conf_comp=0
1382       return
1383       end
1384 C----------------------------------------------------------------------------
1385       double precision function dif_ang(n,x,y)
1386       implicit none
1387       integer i,n
1388       double precision x(n),y(n)
1389       double precision w,wa,dif,difa
1390       double precision pinorm
1391       include 'COMMON.GEO'
1392       wa=0.0D0
1393       difa=0.0D0
1394       do i=1,n
1395         dif=dabs(pinorm(y(i)-x(i)))
1396         if (dabs(dif-dwapi).lt.dif) dif=dabs(dif-dwapi)
1397         w=1.0D0-(2.0D0*(i-1)/(n-1)-1.0D0)**2+1.0D0/n
1398         wa=wa+w
1399         difa=difa+dif*dif*w
1400       enddo
1401       dif_ang=rad2deg*dsqrt(difa/wa)
1402       return
1403       end
1404 c--------------------------------------------------------------------------
1405       subroutine add2cache(n1,n2,ncache,nvar,SourceID,CachSrc,
1406      &                     ecur,xcur,ecache,xcache)
1407       implicit none
1408       include 'COMMON.GEO'
1409       include 'COMMON.IOUNITS'
1410       integer n1,n2,ncache,nvar,SourceID,CachSrc(n2)
1411       integer i,ii,j
1412       double precision ecur,xcur(nvar),ecache(n2),xcache(n1,n2)
1413 cd    write (iout,*) 'Enter ADD2CACHE ncache=',ncache ,' ecur',ecur
1414 cd    write (iout,'(10f8.3)') (rad2deg*xcur(i),i=1,nvar)
1415 cd    write (iout,*) 'Old CACHE array:'
1416 cd    do i=1,ncache
1417 cd      write (iout,*) 'i=',i,' ecache=',ecache(i),' CachSrc',CachSrc(i)
1418 cd      write (iout,'(10f8.3)') (rad2deg*xcache(j,i),j=1,nvar)
1419 cd    enddo
1420
1421       i=ncache
1422       do while (i.gt.0 .and. ecur.lt.ecache(i))
1423         i=i-1
1424       enddo
1425       i=i+1
1426 cd    write (iout,*) 'i=',i,' ncache=',ncache
1427       if (ncache.eq.n2) then
1428         write (iout,*) 'Cache dimension exceeded',ncache,n2
1429         write (iout,*) 'Highest-energy conformation will be removed.'
1430         ncache=ncache-1
1431       endif
1432       do ii=ncache,i,-1
1433         ecache(ii+1)=ecache(ii)
1434         CachSrc(ii+1)=CachSrc(ii)
1435         do j=1,nvar
1436           xcache(j,ii+1)=xcache(j,ii)
1437         enddo
1438       enddo
1439       ecache(i)=ecur
1440       CachSrc(i)=SourceID
1441       do j=1,nvar
1442         xcache(j,i)=xcur(j)
1443       enddo
1444       ncache=ncache+1
1445 cd    write (iout,*) 'New CACHE array:'
1446 cd    do i=1,ncache
1447 cd      write (iout,*) 'i=',i,' ecache=',ecache(i),' CachSrc',CachSrc(i)
1448 cd      write (iout,'(10f8.3)') (rad2deg*xcache(j,i),j=1,nvar)
1449 cd    enddo
1450       return
1451       end
1452 c--------------------------------------------------------------------------
1453       subroutine rm_from_cache(i,n1,n2,ncache,nvar,CachSrc,ecache,
1454      &                         xcache)
1455       implicit none
1456       include 'COMMON.GEO'
1457       include 'COMMON.IOUNITS'
1458       integer n1,n2,ncache,nvar,CachSrc(n2)
1459       integer i,ii,j
1460       double precision ecache(n2),xcache(n1,n2)
1461
1462 cd    write (iout,*) 'Enter RM_FROM_CACHE'
1463 cd    write (iout,*) 'Old CACHE array:'
1464 cd    do ii=1,ncache
1465 cd    write (iout,*)'i=',ii,' ecache=',ecache(ii),' CachSrc',CachSrc(ii)
1466 cd      write (iout,'(10f8.3)') (rad2deg*xcache(j,ii),j=1,nvar)
1467 cd    enddo
1468
1469       do ii=i+1,ncache
1470         ecache(ii-1)=ecache(ii)
1471         CachSrc(ii-1)=CachSrc(ii)
1472         do j=1,nvar
1473           xcache(j,ii-1)=xcache(j,ii)
1474         enddo
1475       enddo
1476       ncache=ncache-1
1477 cd    write (iout,*) 'New CACHE array:'
1478 cd    do i=1,ncache
1479 cd      write (iout,*) 'i=',i,' ecache=',ecache(i),' CachSrc',CachSrc(i)
1480 cd      write (iout,'(10f8.3)') (rad2deg*xcache(j,i),j=1,nvar)
1481 cd    enddo
1482       return
1483       end