source/unres/src_MD-M/readrtns_CSA.F

   1       subroutine readrtns
   2       implicit real*8 (a-h,o-z)
   3       include 'DIMENSIONS'
   4 #ifdef MPI
   5       include 'mpif.h'
   6 #endif
   7       include 'COMMON.SETUP'
   8       include 'COMMON.CONTROL'
   9       include 'COMMON.SBRIDGE'
  10       include 'COMMON.IOUNITS'
  11       logical file_exist
  12 C Read force-field parameters except weights
  13       call parmread
  14 C Read job setup parameters
  15       call read_control
  16 C Read control parameters for energy minimzation if required
  17       if (minim) call read_minim
  18 C Read MCM control parameters if required
  19       if (modecalc.eq.3 .or. modecalc.eq.6) call mcmread
  20 C Read MD control parameters if reqjuired
  21       if (modecalc.eq.12) call read_MDpar
  22 C Read MREMD control parameters if required
  23       if (modecalc.eq.14) then
  24          call read_MDpar
  25          call read_REMDpar
  26       endif
  27 C Read MUCA control parameters if required
  28       if (lmuca) call read_muca
  29 C Read CSA control parameters if required (from fort.40 if exists
  30 C otherwise from general input file)
  31       if (modecalc.eq.8) then
  32        inquire (file="fort.40",exist=file_exist)
  33        if (.not.file_exist) call csaread
  34       endif
  35 cfmc      if (modecalc.eq.10) call mcmfread
  36 C Read molecule information, molecule geometry, energy-term weights, and
  37 C restraints if requested
  38       call molread
  39 C Print restraint information
  40 #ifdef MPI
  41       if (.not. out1file .or. me.eq.king) then
  42 #endif
  43       if (nhpb.gt.nss)
  44      &write (iout,'(a,i5,a)') "The following",nhpb-nss,
  45      & " distance constraints have been imposed"
  46       do i=nss+1,nhpb
  47         write (iout,'(3i6,f10.5)') i-nss,ihpb(i),jhpb(i),forcon(i)
  48       enddo
  49 #ifdef MPI
  50       endif
  51 #endif
  52 c      print *,"Processor",myrank," leaves READRTNS"
  53       return
  54       end
  55 C-------------------------------------------------------------------------------
  56       subroutine read_control
  57 C
  58 C Read contorl data
  59 C
  60       implicit real*8 (a-h,o-z)
  61       include 'DIMENSIONS'
  62 #ifdef MP
  63       include 'mpif.h'
  64       logical OKRandom, prng_restart
  65       real*8  r1
  66 #endif
  67       include 'COMMON.IOUNITS'
  68       include 'COMMON.TIME1'
  69       include 'COMMON.THREAD'
  70       include 'COMMON.SBRIDGE'
  71       include 'COMMON.CONTROL'
  72       include 'COMMON.MCM'
  73       include 'COMMON.MAP'
  74       include 'COMMON.HEADER'
  75       include 'COMMON.CSA'
  76       include 'COMMON.CHAIN'
  77       include 'COMMON.MUCA'
  78       include 'COMMON.MD'
  79       include 'COMMON.FFIELD'
  80       include 'COMMON.INTERACT'
  81       include 'COMMON.SETUP'
  82       COMMON /MACHSW/ KDIAG,ICORFL,IXDR
  83       character*8 diagmeth(0:3) /'Library','EVVRSP','Givens','Jacobi'/
  84       character*80 ucase
  85       character*320 controlcard
  86
  87       nglob_csa=0
  88       eglob_csa=1d99
  89       nmin_csa=0
  90       read (INP,'(a)') titel
  91       call card_concat(controlcard)
  92 c      out1file=index(controlcard,'OUT1FILE').gt.0 .or. fg_rank.gt.0
  93 c      print *,"Processor",me," fg_rank",fg_rank," out1file",out1file
  94       call reada(controlcard,'SEED',seed,0.0D0)
  95       call random_init(seed)
  96 C Set up the time limit (caution! The time must be input in minutes!)
  97       read_cart=index(controlcard,'READ_CART').gt.0
  98       call readi(controlcard,'CONSTR_DIST',constr_dist,0)
  99       call readi(controlcard,'SYM',symetr,1)
 100       call reada(controlcard,'TIMLIM',timlim,960.0D0) ! default 16 hours
 101       unres_pdb = index(controlcard,'UNRES_PDB') .gt. 0
 102       call reada(controlcard,'SAFETY',safety,30.0D0) ! default 30 minutes
 103       call reada(controlcard,'RMSDBC',rmsdbc,3.0D0)
 104       call reada(controlcard,'RMSDBC1',rmsdbc1,0.5D0)
 105       call reada(controlcard,'RMSDBC1MAX',rmsdbc1max,1.5D0)
 106       call reada(controlcard,'RMSDBCM',rmsdbcm,3.0D0)
 107       call reada(controlcard,'DRMS',drms,0.1D0)
 108       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
 109        write (iout,'(a,f10.1)')'RMSDBC = ',rmsdbc
 110        write (iout,'(a,f10.1)')'RMSDBC1 = ',rmsdbc1
 111        write (iout,'(a,f10.1)')'RMSDBC1MAX = ',rmsdbc1max
 112        write (iout,'(a,f10.1)')'DRMS    = ',drms
 113        write (iout,'(a,f10.1)')'RMSDBCM = ',rmsdbcm
 114        write (iout,'(a,f10.1)') 'Time limit (min):',timlim
 115       endif
 116       call readi(controlcard,'NZ_START',nz_start,0)
 117       call readi(controlcard,'NZ_END',nz_end,0)
 118       call readi(controlcard,'IZ_SC',iz_sc,0)
 119       timlim=60.0D0*timlim
 120       safety = 60.0d0*safety
 121       timem=timlim
 122       modecalc=0
 123       call reada(controlcard,"T_BATH",t_bath,300.0d0)
 124       minim=(index(controlcard,'MINIMIZE').gt.0)
 125       dccart=(index(controlcard,'CART').gt.0)
 126       overlapsc=(index(controlcard,'OVERLAP').gt.0)
 127       overlapsc=.not.overlapsc
 128       searchsc=(index(controlcard,'NOSEARCHSC').gt.0)
 129       searchsc=.not.searchsc
 130       sideadd=(index(controlcard,'SIDEADD').gt.0)
 131       energy_dec=(index(controlcard,'ENERGY_DEC').gt.0)
 132       outpdb=(index(controlcard,'PDBOUT').gt.0)
 133       outmol2=(index(controlcard,'MOL2OUT').gt.0)
 134       pdbref=(index(controlcard,'PDBREF').gt.0)
 135       refstr=pdbref .or. (index(controlcard,'REFSTR').gt.0)
 136       indpdb=index(controlcard,'PDBSTART')
 137       extconf=(index(controlcard,'EXTCONF').gt.0)
 138       call readi(controlcard,'IPRINT',iprint,0)
 139       call readi(controlcard,'MAXGEN',maxgen,10000)
 140       call readi(controlcard,'MAXOVERLAP',maxoverlap,1000)
 141       call readi(controlcard,"KDIAG",kdiag,0)
 142       call readi(controlcard,"RESCALE_MODE",rescale_mode,2)
 143       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0)
 144      & write (iout,*) "RESCALE_MODE",rescale_mode
 145       split_ene=index(controlcard,'SPLIT_ENE').gt.0
 146       if (index(controlcard,'REGULAR').gt.0.0D0) then
 147         call reada(controlcard,'WEIDIS',weidis,0.1D0)
 148         modecalc=1
 149         refstr=.true.
 150       endif
 151       if (index(controlcard,'CHECKGRAD').gt.0) then
 152         modecalc=5
 153         if (index(controlcard,'CART').gt.0) then
 154           icheckgrad=1
 155         elseif (index(controlcard,'CARINT').gt.0) then
 156           icheckgrad=2
 157         else
 158           icheckgrad=3
 159         endif
 160       elseif (index(controlcard,'THREAD').gt.0) then
 161         modecalc=2
 162         call readi(controlcard,'THREAD',nthread,0)
 163         if (nthread.gt.0) then
 164           call reada(controlcard,'WEIDIS',weidis,0.1D0)
 165         else
 166           if (fg_rank.eq.0)
 167      &    write (iout,'(a)')'A number has to follow the THREAD keyword.'
 168           stop 'Error termination in Read_Control.'
 169         endif
 170       else if (index(controlcard,'MCMA').gt.0) then
 171         modecalc=3
 172       else if (index(controlcard,'MCEE').gt.0) then
 173         modecalc=6
 174       else if (index(controlcard,'MULTCONF').gt.0) then
 175         modecalc=4
 176       else if (index(controlcard,'MAP').gt.0) then
 177         modecalc=7
 178         call readi(controlcard,'MAP',nmap,0)
 179       else if (index(controlcard,'CSA').gt.0) then
 180         modecalc=8
 181 crc      else if (index(controlcard,'ZSCORE').gt.0) then
 182 crc
 183 crc  ZSCORE is rm from UNRES, modecalc=9 is available
 184 crc
 185 crc        modecalc=9
 186 cfcm      else if (index(controlcard,'MCMF').gt.0) then
 187 cfmc        modecalc=10
 188       else if (index(controlcard,'SOFTREG').gt.0) then
 189         modecalc=11
 190       else if (index(controlcard,'CHECK_BOND').gt.0) then
 191         modecalc=-1
 192       else if (index(controlcard,'TEST').gt.0) then
 193         modecalc=-2
 194       else if (index(controlcard,'MD').gt.0) then
 195         modecalc=12
 196       else if (index(controlcard,'RE ').gt.0) then
 197         modecalc=14
 198       endif
 199
 200       lmuca=index(controlcard,'MUCA').gt.0
 201       call readi(controlcard,'MUCADYN',mucadyn,0)
 202       call readi(controlcard,'MUCASMOOTH',muca_smooth,0)
 203       if (lmuca .and. (me.eq.king .or. .not.out1file ))
 204      & then
 205        write (iout,*) 'MUCADYN=',mucadyn
 206        write (iout,*) 'MUCASMOOTH=',muca_smooth
 207       endif
 208
 209       iscode=index(controlcard,'ONE_LETTER')
 210       indphi=index(controlcard,'PHI')
 211       indback=index(controlcard,'BACK')
 212       iranconf=index(controlcard,'RAND_CONF')
 213       i2ndstr=index(controlcard,'USE_SEC_PRED')
 214       gradout=index(controlcard,'GRADOUT').gt.0
 215       gnorm_check=index(controlcard,'GNORM_CHECK').gt.0
 216       call reada(controlcard,'DISTCHAINMAX',distchainmax,5.0d0)
 217       if (me.eq.king .or. .not.out1file )
 218      &  write (iout,*) "DISTCHAINMAX",distchainmax
 219
 220       if(me.eq.king.or..not.out1file)
 221      & write (iout,'(2a)') diagmeth(kdiag),
 222      &  ' routine used to diagonalize matrices.'
 223       return
 224       end
 225 c--------------------------------------------------------------------------
 226       subroutine read_REMDpar
 227 C
 228 C Read REMD settings
 229 C
 230       implicit real*8 (a-h,o-z)
 231       include 'DIMENSIONS'
 232       include 'COMMON.IOUNITS'
 233       include 'COMMON.TIME1'
 234       include 'COMMON.MD'
 235 #ifndef LANG0
 236       include 'COMMON.LANGEVIN'
 237 #else
 238       include 'COMMON.LANGEVIN.lang0'
 239 #endif
 240       include 'COMMON.INTERACT'
 241       include 'COMMON.NAMES'
 242       include 'COMMON.GEO'
 243       include 'COMMON.REMD'
 244       include 'COMMON.CONTROL'
 245       include 'COMMON.SETUP'
 246       character*80 ucase
 247       character*320 controlcard
 248       character*3200 controlcard1
 249       integer iremd_m_total
 250
 251       if(me.eq.king.or..not.out1file)
 252      & write (iout,*) "REMD setup"
 253
 254       call card_concat(controlcard)
 255       call readi(controlcard,"NREP",nrep,3)
 256       call readi(controlcard,"NSTEX",nstex,1000)
 257       call reada(controlcard,"RETMIN",retmin,10.0d0)
 258       call reada(controlcard,"RETMAX",retmax,1000.0d0)
 259       mremdsync=(index(controlcard,'SYNC').gt.0)
 260       call readi(controlcard,"NSYN",i_sync_step,100)
 261       restart1file=(index(controlcard,'REST1FILE').gt.0)
 262       traj1file=(index(controlcard,'TRAJ1FILE').gt.0)
 263       call readi(controlcard,"TRAJCACHE",max_cache_traj_use,1)
 264       if(max_cache_traj_use.gt.max_cache_traj)
 265      &           max_cache_traj_use=max_cache_traj
 266       if(me.eq.king.or..not.out1file) then
 267 cd       if (traj1file) then
 268 crc caching is in testing - NTWX is not ignored
 269 cd        write (iout,*) "NTWX value is ignored"
 270 cd        write (iout,*) "  trajectory is stored to one file by master"
 271 cd        write (iout,*) "  before exchange at NSTEX intervals"
 272 cd       endif
 273        write (iout,*) "NREP= ",nrep
 274        write (iout,*) "NSTEX= ",nstex
 275        write (iout,*) "SYNC= ",mremdsync
 276        write (iout,*) "NSYN= ",i_sync_step
 277        write (iout,*) "TRAJCACHE= ",max_cache_traj_use
 278       endif
 279       remd_tlist=.false.
 280       if (index(controlcard,'TLIST').gt.0) then
 281          remd_tlist=.true.
 282          call card_concat(controlcard1)
 283          read(controlcard1,*) (remd_t(i),i=1,nrep)
 284          if(me.eq.king.or..not.out1file)
 285      &    write (iout,*)'tlist',(remd_t(i),i=1,nrep)
 286       endif
 287       remd_mlist=.false.
 288       if (index(controlcard,'MLIST').gt.0) then
 289          remd_mlist=.true.
 290          call card_concat(controlcard1)
 291          read(controlcard1,*) (remd_m(i),i=1,nrep)
 292          if(me.eq.king.or..not.out1file) then
 293           write (iout,*)'mlist',(remd_m(i),i=1,nrep)
 294           iremd_m_total=0
 295           do i=1,nrep
 296            iremd_m_total=iremd_m_total+remd_m(i)
 297           enddo
 298            write (iout,*) 'Total number of replicas ',iremd_m_total
 299           endif
 300          endif
 301       if(me.eq.king.or..not.out1file)
 302      &   write (iout,'(/30(1h=),a,29(1h=)/)') " End of REMD run setup "
 303       return
 304       end
 305 c--------------------------------------------------------------------------
 306       subroutine read_MDpar
 307 C
 308 C Read MD settings
 309 C
 310       implicit real*8 (a-h,o-z)
 311       include 'DIMENSIONS'
 312       include 'COMMON.IOUNITS'
 313       include 'COMMON.TIME1'
 314       include 'COMMON.MD'
 315 #ifndef LANG0
 316       include 'COMMON.LANGEVIN'
 317 #else
 318       include 'COMMON.LANGEVIN.lang0'
 319 #endif
 320       include 'COMMON.INTERACT'
 321       include 'COMMON.NAMES'
 322       include 'COMMON.GEO'
 323       include 'COMMON.SETUP'
 324       include 'COMMON.CONTROL'
 325       include 'COMMON.SPLITELE'
 326       character*80 ucase
 327       character*320 controlcard
 328
 329       call card_concat(controlcard)
 330       call readi(controlcard,"NSTEP",n_timestep,1000000)
 331       call readi(controlcard,"NTWE",ntwe,100)
 332       call readi(controlcard,"NTWX",ntwx,1000)
 333       call reada(controlcard,"DT",d_time,1.0d-1)
 334       call reada(controlcard,"DVMAX",dvmax,2.0d1)
 335       call reada(controlcard,"DAMAX",damax,1.0d1)
 336       call reada(controlcard,"EDRIFTMAX",edriftmax,1.0d+1)
 337       call readi(controlcard,"LANG",lang,0)
 338       RESPA = index(controlcard,"RESPA") .gt. 0
 339       call readi(controlcard,"NTIME_SPLIT",ntime_split,1)
 340       ntime_split0=ntime_split
 341       call readi(controlcard,"MAXTIME_SPLIT",maxtime_split,64)
 342       ntime_split0=ntime_split
 343       call reada(controlcard,"R_CUT",r_cut,2.0d0)
 344       call reada(controlcard,"LAMBDA",rlamb,0.3d0)
 345       rest = index(controlcard,"REST").gt.0
 346       tbf = index(controlcard,"TBF").gt.0
 347       usampl = index(controlcard,"USAMPL").gt.0
 348
 349       mdpdb = index(controlcard,"MDPDB").gt.0
 350       call reada(controlcard,"T_BATH",t_bath,300.0d0)
 351       call reada(controlcard,"TAU_BATH",tau_bath,1.0d-1)
 352       call reada(controlcard,"EQ_TIME",eq_time,1.0d+4)
 353       call readi(controlcard,"RESET_MOMENT",count_reset_moment,1000)
 354       if (count_reset_moment.eq.0) count_reset_moment=1000000000
 355       call readi(controlcard,"RESET_VEL",count_reset_vel,1000)
 356       reset_moment=lang.eq.0 .and. tbf .and. count_reset_moment.gt.0
 357       reset_vel=lang.eq.0 .and. tbf .and. count_reset_vel.gt.0
 358       if (count_reset_vel.eq.0) count_reset_vel=1000000000
 359       large = index(controlcard,"LARGE").gt.0
 360       print_compon = index(controlcard,"PRINT_COMPON").gt.0
 361       rattle = index(controlcard,"RATTLE").gt.0
 362 c  if performing umbrella sampling, fragments constrained are read from the fragment file
 363       nset=0
 364       if(usampl) then
 365         call read_fragments
 366       endif
 367
 368       if(me.eq.king.or..not.out1file) then
 369        write (iout,*)
 370        write (iout,'(27(1h=),a26,27(1h=))') " Parameters of the MD run "
 371        write (iout,*)
 372        write (iout,'(a)') "The units are:"
 373        write (iout,'(a)') "positions: angstrom, time: 48.9 fs"
 374        write (iout,'(2a)') "velocity: angstrom/(48.9 fs),",
 375      &  " acceleration: angstrom/(48.9 fs)**2"
 376        write (iout,'(a)') "energy: kcal/mol, temperature: K"
 377        write (iout,*)
 378        write (iout,'(a60,i10)') "Number of time steps:",n_timestep
 379        write (iout,'(a60,f10.5,a)')
 380      &  "Initial time step of numerical integration:",d_time,
 381      &  " natural units"
 382        write (iout,'(60x,f10.5,a)') d_time*48.9," fs"
 383        if (RESPA) then
 384         write (iout,'(2a,i4,a)')
 385      &    "A-MTS algorithm used; initial time step for fast-varying",
 386      &    " short-range forces split into",ntime_split," steps."
 387         write (iout,'(a,f5.2,a,f5.2)') "Short-range force cutoff",
 388      &   r_cut," lambda",rlamb
 389        endif
 390        write (iout,'(2a,f10.5)')
 391      &  "Maximum acceleration threshold to reduce the time step",
 392      &  "/increase split number:",damax
 393        write (iout,'(2a,f10.5)')
 394      &  "Maximum predicted energy drift to reduce the timestep",
 395      &  "/increase split number:",edriftmax
 396        write (iout,'(a60,f10.5)')
 397      & "Maximum velocity threshold to reduce velocities:",dvmax
 398        write (iout,'(a60,i10)') "Frequency of property output:",ntwe
 399        write (iout,'(a60,i10)') "Frequency of coordinate output:",ntwx
 400        if (rattle) write (iout,'(a60)')
 401      &  "Rattle algorithm used to constrain the virtual bonds"
 402       endif
 403       reset_fricmat=1000
 404       if (lang.gt.0) then
 405         call reada(controlcard,"ETAWAT",etawat,0.8904d0)
 406         call reada(controlcard,"RWAT",rwat,1.4d0)
 407         call reada(controlcard,"SCAL_FRIC",scal_fric,2.0d-2)
 408         surfarea=index(controlcard,"SURFAREA").gt.0
 409         call readi(controlcard,"RESET_FRICMAT",reset_fricmat,1000)
 410         if(me.eq.king.or..not.out1file)then
 411          write (iout,'(/a,$)') "Langevin dynamics calculation"
 412          if (lang.eq.1) then
 413           write (iout,'(a/)')
 414      &      " with direct integration of Langevin equations"
 415          else if (lang.eq.2) then
 416           write (iout,'(a/)') " with TINKER stochasic MD integrator"
 417          else if (lang.eq.3) then
 418           write (iout,'(a/)') " with Ciccotti's stochasic MD integrator"
 419          else if (lang.eq.4) then
 420           write (iout,'(a/)') " in overdamped mode"
 421          else
 422           write (iout,'(//a,i5)')
 423      &      "=========== ERROR: Unknown Langevin dynamics mode:",lang
 424           stop
 425          endif
 426          write (iout,'(a60,f10.5)') "Temperature:",t_bath
 427          write (iout,'(a60,f10.5)') "Viscosity of the solvent:",etawat
 428          write (iout,'(a60,f10.5)') "Radius of solvent molecule:",rwat
 429          write (iout,'(a60,f10.5)')
 430      &   "Scaling factor of the friction forces:",scal_fric
 431          if (surfarea) write (iout,'(2a,i10,a)')
 432      &     "Friction coefficients will be scaled by solvent-accessible",
 433      &     " surface area every",reset_fricmat," steps."
 434         endif
 435 c Calculate friction coefficients and bounds of stochastic forces
 436         eta=6*pi*cPoise*etawat
 437         if(me.eq.king.or..not.out1file)
 438      &   write(iout,'(a60,f10.5)')"Eta of the solvent in natural units:"
 439      &   ,eta
 440         gamp=scal_fric*(pstok+rwat)*eta
 441         stdfp=dsqrt(2*Rb*t_bath/d_time)
 442         do i=1,ntyp
 443           gamsc(i)=scal_fric*(restok(i)+rwat)*eta
 444           stdfsc(i)=dsqrt(2*Rb*t_bath/d_time)
 445         enddo
 446         if(me.eq.king.or..not.out1file)then
 447          write (iout,'(/2a/)')
 448      &   "Radii of site types and friction coefficients and std's of",
 449      &   " stochastic forces of fully exposed sites"
 450          write (iout,'(a5,f5.2,2f10.5)')'p',pstok,gamp,stdfp*dsqrt(gamp)
 451          do i=1,ntyp
 452           write (iout,'(a5,f5.2,2f10.5)') restyp(i),restok(i),
 453      &     gamsc(i),stdfsc(i)*dsqrt(gamsc(i))
 454          enddo
 455         endif
 456       else if (tbf) then
 457         if(me.eq.king.or..not.out1file)then
 458          write (iout,'(a)') "Berendsen bath calculation"
 459          write (iout,'(a60,f10.5)') "Temperature:",t_bath
 460          write (iout,'(a60,f10.5)') "Coupling constant (tau):",tau_bath
 461          if (reset_moment)
 462      &   write (iout,'(a,i10,a)') "Momenta will be reset at zero every",
 463      &   count_reset_moment," steps"
 464          if (reset_vel)
 465      &    write (iout,'(a,i10,a)')
 466      &    "Velocities will be reset at random every",count_reset_vel,
 467      &   " steps"
 468         endif
 469       else
 470         if(me.eq.king.or..not.out1file)
 471      &   write (iout,'(a31)') "Microcanonical mode calculation"
 472       endif
 473       if(me.eq.king.or..not.out1file)then
 474        if (rest) write (iout,'(/a/)') "===== Calculation restarted ===="
 475        if (usampl) then
 476           write(iout,*) "MD running with constraints."
 477           write(iout,*) "Equilibration time ", eq_time, " mtus."
 478           write(iout,*) "Constraining ", nfrag," fragments."
 479           write(iout,*) "Length of each fragment, weight and q0:"
 480           do iset=1,nset
 481            write (iout,*) "Set of restraints #",iset
 482            do i=1,nfrag
 483               write(iout,'(2i5,f8.1,f7.4)') ifrag(1,i,iset),
 484      &           ifrag(2,i,iset),wfrag(i,iset),qinfrag(i,iset)
 485            enddo
 486            write(iout,*) "constraints between ", npair, "fragments."
 487            write(iout,*) "constraint pairs, weights and q0:"
 488            do i=1,npair
 489             write(iout,'(2i5,f8.1,f7.4)') ipair(1,i,iset),
 490      &             ipair(2,i,iset),wpair(i,iset),qinpair(i,iset)
 491            enddo
 492            write(iout,*) "angle constraints within ", nfrag_back,
 493      &      "backbone fragments."
 494            write(iout,*) "fragment, weights:"
 495            do i=1,nfrag_back
 496             write(iout,'(2i5,3f8.1)') ifrag_back(1,i,iset),
 497      &         ifrag_back(2,i,iset),wfrag_back(1,i,iset),
 498      &         wfrag_back(2,i,iset),wfrag_back(3,i,iset)
 499            enddo
 500           enddo
 501         iset=mod(kolor,nset)+1
 502        endif
 503       endif
 504       if(me.eq.king.or..not.out1file)
 505      & write (iout,'(/30(1h=),a,29(1h=)/)') " End of MD run setup "
 506       return
 507       end
 508 c------------------------------------------------------------------------------
 509       subroutine molread
 510 C
 511 C Read molecular data.
 512 C
 513       implicit real*8 (a-h,o-z)
 514       include 'DIMENSIONS'
 515 #ifdef MPI
 516       include 'mpif.h'
 517       integer error_msg
 518 #endif
 519       include 'COMMON.IOUNITS'
 520       include 'COMMON.GEO'
 521       include 'COMMON.VAR'
 522       include 'COMMON.INTERACT'
 523       include 'COMMON.LOCAL'
 524       include 'COMMON.NAMES'
 525       include 'COMMON.CHAIN'
 526       include 'COMMON.FFIELD'
 527       include 'COMMON.SBRIDGE'
 528       include 'COMMON.HEADER'
 529       include 'COMMON.CONTROL'
 530       include 'COMMON.DBASE'
 531       include 'COMMON.THREAD'
 532       include 'COMMON.CONTACTS'
 533       include 'COMMON.TORCNSTR'
 534       include 'COMMON.TIME1'
 535       include 'COMMON.BOUNDS'
 536       include 'COMMON.MD'
 537       include 'COMMON.SETUP'
 538       character*4 sequence(maxres)
 539       integer rescode
 540       double precision x(maxvar)
 541       character*256 pdbfile
 542       character*320 weightcard
 543       character*80 weightcard_t,ucase
 544       dimension itype_pdb(maxres)
 545       common /pizda/ itype_pdb
 546       logical seq_comp,fail
 547       double precision energia(0:n_ene)
 548       integer ilen
 549       external ilen
 550 C
 551 C Body
 552 C
 553 C Read weights of the subsequent energy terms.
 554        call card_concat(weightcard)
 555        call reada(weightcard,'WLONG',wlong,1.0D0)
 556        call reada(weightcard,'WSC',wsc,wlong)
 557        call reada(weightcard,'WSCP',wscp,wlong)
 558        call reada(weightcard,'WELEC',welec,1.0D0)
 559        call reada(weightcard,'WVDWPP',wvdwpp,welec)
 560        call reada(weightcard,'WEL_LOC',wel_loc,1.0D0)
 561        call reada(weightcard,'WCORR4',wcorr4,0.0D0)
 562        call reada(weightcard,'WCORR5',wcorr5,0.0D0)
 563        call reada(weightcard,'WCORR6',wcorr6,0.0D0)
 564        call reada(weightcard,'WTURN3',wturn3,1.0D0)
 565        call reada(weightcard,'WTURN4',wturn4,1.0D0)
 566        call reada(weightcard,'WTURN6',wturn6,1.0D0)
 567        call reada(weightcard,'WSCCOR',wsccor,1.0D0)
 568        call reada(weightcard,'WSTRAIN',wstrain,1.0D0)
 569        call reada(weightcard,'WBOND',wbond,1.0D0)
 570        call reada(weightcard,'WTOR',wtor,1.0D0)
 571        call reada(weightcard,'WTORD',wtor_d,1.0D0)
 572        call reada(weightcard,'WANG',wang,1.0D0)
 573        call reada(weightcard,'WSCLOC',wscloc,1.0D0)
 574        call reada(weightcard,'SCAL14',scal14,0.4D0)
 575        call reada(weightcard,'SCALSCP',scalscp,1.0d0)
 576        call reada(weightcard,'CUTOFF',cutoff_corr,7.0d0)
 577        call reada(weightcard,'DELT_CORR',delt_corr,0.5d0)
 578        call reada(weightcard,'TEMP0',temp0,300.0d0)
 579        if (index(weightcard,'SOFT').gt.0) ipot=6
 580 C 12/1/95 Added weight for the multi-body term WCORR
 581        call reada(weightcard,'WCORRH',wcorr,1.0D0)
 582        if (wcorr4.gt.0.0d0) wcorr=wcorr4
 583        weights(1)=wsc
 584        weights(2)=wscp
 585        weights(3)=welec
 586        weights(4)=wcorr
 587        weights(5)=wcorr5
 588        weights(6)=wcorr6
 589        weights(7)=wel_loc
 590        weights(8)=wturn3
 591        weights(9)=wturn4
 592        weights(10)=wturn6
 593        weights(11)=wang
 594        weights(12)=wscloc
 595        weights(13)=wtor
 596        weights(14)=wtor_d
 597        weights(15)=wstrain
 598        weights(16)=wvdwpp
 599        weights(17)=wbond
 600        weights(18)=scal14
 601        weights(21)=wsccor
 602       if(me.eq.king.or..not.out1file)
 603      & write (iout,10) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
 604      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
 605      &  wturn4,wturn6
 606    10 format (/'Energy-term weights (unscaled):'//
 607      & 'WSCC=   ',f10.6,' (SC-SC)'/
 608      & 'WSCP=   ',f10.6,' (SC-p)'/
 609      & 'WELEC=  ',f10.6,' (p-p electr)'/
 610      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
 611      & 'WBOND=  ',f10.6,' (stretching)'/
 612      & 'WANG=   ',f10.6,' (bending)'/
 613      & 'WSCLOC= ',f10.6,' (SC local)'/
 614      & 'WTOR=   ',f10.6,' (torsional)'/
 615      & 'WTORD=  ',f10.6,' (double torsional)'/
 616      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
 617      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
 618      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
 619      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
 620      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
 621      & 'WSCCOR= ',f10.6,' (back-scloc correlation)'/
 622      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
 623      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
 624      & 'WTURN6= ',f10.6,' (turns, 6th order)')
 625       if(me.eq.king.or..not.out1file)then
 626        if (wcorr4.gt.0.0d0) then
 627         write (iout,'(/2a/)') 'Local-electrostatic type correlation ',
 628      &   'between contact pairs of peptide groups'
 629         write (iout,'(2(a,f5.3/))')
 630      &  'Cutoff on 4-6th order correlation terms: ',cutoff_corr,
 631      &  'Range of quenching the correlation terms:',2*delt_corr
 632        else if (wcorr.gt.0.0d0) then
 633         write (iout,'(/2a/)') 'Hydrogen-bonding correlation ',
 634      &   'between contact pairs of peptide groups'
 635        endif
 636        write (iout,'(a,f8.3)')
 637      &  'Scaling factor of 1,4 SC-p interactions:',scal14
 638        write (iout,'(a,f8.3)')
 639      &  'General scaling factor of SC-p interactions:',scalscp
 640       endif
 641       r0_corr=cutoff_corr-delt_corr
 642       do i=1,ntyp
 643         aad(i,1)=scalscp*aad(i,1)
 644         aad(i,2)=scalscp*aad(i,2)
 645         bad(i,1)=scalscp*bad(i,1)
 646         bad(i,2)=scalscp*bad(i,2)
 647       enddo
 648       call rescale_weights(t_bath)
 649       if(me.eq.king.or..not.out1file)
 650      & write (iout,22) wsc,wscp,welec,wvdwpp,wbond,wang,wscloc,wtor,
 651      &  wtor_d,wstrain,wel_loc,wcorr,wcorr5,wcorr6,wsccor,wturn3,
 652      &  wturn4,wturn6
 653    22 format (/'Energy-term weights (scaled):'//
 654      & 'WSCC=   ',f10.6,' (SC-SC)'/
 655      & 'WSCP=   ',f10.6,' (SC-p)'/
 656      & 'WELEC=  ',f10.6,' (p-p electr)'/
 657      & 'WVDWPP= ',f10.6,' (p-p VDW)'/
 658      & 'WBOND=  ',f10.6,' (stretching)'/
 659      & 'WANG=   ',f10.6,' (bending)'/
 660      & 'WSCLOC= ',f10.6,' (SC local)'/
 661      & 'WTOR=   ',f10.6,' (torsional)'/
 662      & 'WTORD=  ',f10.6,' (double torsional)'/
 663      & 'WSTRAIN=',f10.6,' (SS bridges & dist. cnstr.)'/
 664      & 'WEL_LOC=',f10.6,' (multi-body 3-rd order)'/
 665      & 'WCORR4= ',f10.6,' (multi-body 4th order)'/
 666      & 'WCORR5= ',f10.6,' (multi-body 5th order)'/
 667      & 'WCORR6= ',f10.6,' (multi-body 6th order)'/
 668      & 'WSCCOR= ',f10.6,' (back-scloc correlatkion)'/
 669      & 'WTURN3= ',f10.6,' (turns, 3rd order)'/
 670      & 'WTURN4= ',f10.6,' (turns, 4th order)'/
 671      & 'WTURN6= ',f10.6,' (turns, 6th order)')
 672       if(me.eq.king.or..not.out1file)
 673      & write (iout,*) "Reference temperature for weights calculation:",
 674      &  temp0
 675       call reada(weightcard,"D0CM",d0cm,3.78d0)
 676       call reada(weightcard,"AKCM",akcm,15.1d0)
 677       call reada(weightcard,"AKTH",akth,11.0d0)
 678       call reada(weightcard,"AKCT",akct,12.0d0)
 679       call reada(weightcard,"V1SS",v1ss,-1.08d0)
 680       call reada(weightcard,"V2SS",v2ss,7.61d0)
 681       call reada(weightcard,"V3SS",v3ss,13.7d0)
 682       call reada(weightcard,"EBR",ebr,-5.50D0)
 683       dyn_ss=(index(weightcard,'DYN_SS').gt.0)
 684       do i=1,maxres
 685         dyn_ss_mask(i)=.false.
 686       enddo
 687       do i=1,maxres-1
 688         do j=i+1,maxres
 689           dyn_ssbond_ij(i,j)=1.0d300
 690         enddo
 691       enddo
 692       call reada(weightcard,"HT",Ht,0.0D0)
 693       if (dyn_ss) then
 694         ss_depth=ebr/wsc-0.25*eps(1,1)
 695         Ht=Ht/wsc-0.25*eps(1,1)
 696         akcm=akcm*wstrain/wsc
 697         akth=akth*wstrain/wsc
 698         akct=akct*wstrain/wsc
 699         v1ss=v1ss*wstrain/wsc
 700         v2ss=v2ss*wstrain/wsc
 701         v3ss=v3ss*wstrain/wsc
 702       else
 703         ss_depth=ebr/wstrain-0.25*eps(1,1)*wsc/wstrain
 704       endif
 705
 706       if(me.eq.king.or..not.out1file) then
 707        write (iout,*) "Parameters of the SS-bond potential:"
 708        write (iout,*) "D0CM",d0cm," AKCM",akcm," AKTH",akth,
 709      & " AKCT",akct
 710        write (iout,*) "V1SS",v1ss," V2SS",v2ss," V3SS",v3ss
 711        write (iout,*) "EBR",ebr," SS_DEPTH",ss_depth
 712        write (iout,*)" HT",Ht
 713        print *,'indpdb=',indpdb,' pdbref=',pdbref
 714       endif
 715       if (indpdb.gt.0 .or. pdbref) then
 716         read(inp,'(a)') pdbfile
 717         if(me.eq.king.or..not.out1file)
 718      &   write (iout,'(2a)') 'PDB data will be read from file ',
 719      &   pdbfile(:ilen(pdbfile))
 720         open(ipdbin,file=pdbfile,status='old',err=33)
 721         goto 34
 722   33    write (iout,'(a)') 'Error opening PDB file.'
 723         stop
 724   34    continue
 725 c        write (iout,*) 'Begin reading pdb data'
 726 c        call flush(iout)
 727         call readpdb
 728 c        write (iout,*) 'Finished reading pdb data'
 729 c        call flush(iout)
 730         if(me.eq.king.or..not.out1file)
 731      &   write (iout,'(a,i3,a,i3)')'nsup=',nsup,
 732      &   ' nstart_sup=',nstart_sup
 733         do i=1,nres
 734           itype_pdb(i)=itype(i)
 735         enddo
 736         close (ipdbin)
 737         nnt=nstart_sup
 738         nct=nstart_sup+nsup-1
 739         call contact(.false.,ncont_ref,icont_ref,co)
 740
 741         if (sideadd) then
 742          if(me.eq.king.or..not.out1file)
 743      &    write(iout,*)'Adding sidechains'
 744          maxsi=1000
 745          do i=2,nres-1
 746           iti=itype(i)
 747           if (iti.ne.10 .and. itype(i).ne.ntyp1) then
 748             nsi=0
 749             fail=.true.
 750             do while (fail.and.nsi.le.maxsi)
 751               call gen_side(iti,theta(i+1),alph(i),omeg(i),fail)
 752               nsi=nsi+1
 753             enddo
 754             if(fail) write(iout,*)'Adding sidechain failed for res ',
 755      &              i,' after ',nsi,' trials'
 756           endif
 757          enddo
 758         endif
 759       endif
 760       if (indpdb.eq.0) then
 761 C Read sequence if not taken from the pdb file.
 762         read (inp,*) nres
 763 c        print *,'nres=',nres
 764         if (iscode.gt.0) then
 765           read (inp,'(80a1)') (sequence(i)(1:1),i=1,nres)
 766         else
 767           read (inp,'(20(1x,a3))') (sequence(i),i=1,nres)
 768         endif
 769 C Convert sequence to numeric code
 770         do i=1,nres
 771           itype(i)=rescode(i,sequence(i),iscode)
 772         enddo
 773 C Assign initial virtual bond lengths
 774         do i=2,nres
 775           vbld(i)=vbl
 776           vbld_inv(i)=vblinv
 777         enddo
 778         do i=2,nres-1
 779           vbld(i+nres)=dsc(iabs(itype(i)))
 780           vbld_inv(i+nres)=dsc_inv(iabs(itype(i)))
 781 c          write (iout,*) "i",i," itype",itype(i),
 782 c     &      " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
 783         enddo
 784       endif
 785 c      print *,nres
 786 c      print '(20i4)',(itype(i),i=1,nres)
 787       do i=1,nres
 788 #ifdef PROCOR
 789         if (itype(i).eq.ntyp1 .or. itype(i+1).eq.ntyp1) then
 790 #else
 791         if (itype(i).eq.ntyp1) then
 792 #endif
 793           itel(i)=0
 794 #ifdef PROCOR
 795         else if (iabs(itype(i+1)).ne.20) then
 796 #else
 797         else if (iabs(itype(i)).ne.20) then
 798 #endif
 799           itel(i)=1
 800         else
 801           itel(i)=2
 802         endif
 803       enddo
 804       if(me.eq.king.or..not.out1file)then
 805        write (iout,*) "ITEL"
 806        do i=1,nres-1
 807          write (iout,*) i,itype(i),itel(i)
 808        enddo
 809        print *,'Call Read_Bridge.'
 810       endif
 811       call read_bridge
 812 C 8/13/98 Set limits to generating the dihedral angles
 813       do i=1,nres
 814         phibound(1,i)=-pi
 815         phibound(2,i)=pi
 816       enddo
 817       read (inp,*) ndih_constr
 818       if (ndih_constr.gt.0) then
 819         read (inp,*) ftors
 820         read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
 821         if(me.eq.king.or..not.out1file)then
 822          write (iout,*)
 823      &   'There are',ndih_constr,' constraints on phi angles.'
 824          do i=1,ndih_constr
 825           write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
 826          enddo
 827         endif
 828         do i=1,ndih_constr
 829           phi0(i)=deg2rad*phi0(i)
 830           drange(i)=deg2rad*drange(i)
 831         enddo
 832         if(me.eq.king.or..not.out1file)
 833      &   write (iout,*) 'FTORS',ftors
 834         do i=1,ndih_constr
 835           ii = idih_constr(i)
 836           phibound(1,ii) = phi0(i)-drange(i)
 837           phibound(2,ii) = phi0(i)+drange(i)
 838         enddo
 839       endif
 840       nnt=1
 841 #ifdef MPI
 842       if (me.eq.king) then
 843 #endif
 844        write (iout,'(a)') 'Boundaries in phi angle sampling:'
 845        do i=1,nres
 846          write (iout,'(a3,i5,2f10.1)')
 847      &   restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
 848        enddo
 849 #ifdef MP
 850       endif
 851 #endif
 852       nct=nres
 853 cd      print *,'NNT=',NNT,' NCT=',NCT
 854       if (itype(1).eq.ntyp1) nnt=2
 855       if (itype(nres).eq.ntyp1) nct=nct-1
 856       if (pdbref) then
 857         if(me.eq.king.or..not.out1file)
 858      &   write (iout,'(a,i3)') 'nsup=',nsup
 859         nstart_seq=nnt
 860         if (nsup.le.(nct-nnt+1)) then
 861           do i=0,nct-nnt+1-nsup
 862             if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
 863               nstart_seq=nnt+i
 864               goto 111
 865             endif
 866           enddo
 867           write (iout,'(a)')
 868      &            'Error - sequences to be superposed do not match.'
 869           stop
 870         else
 871           do i=0,nsup-(nct-nnt+1)
 872             if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
 873      &      then
 874               nstart_sup=nstart_sup+i
 875               nsup=nct-nnt+1
 876               goto 111
 877             endif
 878           enddo
 879           write (iout,'(a)')
 880      &            'Error - sequences to be superposed do not match.'
 881         endif
 882   111   continue
 883         if (nsup.eq.0) nsup=nct-nnt
 884         if (nstart_sup.eq.0) nstart_sup=nnt
 885         if (nstart_seq.eq.0) nstart_seq=nnt
 886         if(me.eq.king.or..not.out1file)
 887      &   write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
 888      &                 ' nstart_seq=',nstart_seq
 889       endif
 890 c--- Zscore rms -------
 891       if (nz_start.eq.0) nz_start=nnt
 892       if (nz_end.eq.0 .and. nsup.gt.0) then
 893         nz_end=nnt+nsup-1
 894       else if (nz_end.eq.0) then
 895         nz_end=nct
 896       endif
 897       if(me.eq.king.or..not.out1file)then
 898        write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
 899        write (iout,*) 'IZ_SC=',iz_sc
 900       endif
 901 c----------------------
 902       call init_int_table
 903       if (refstr) then
 904         if (.not.pdbref) then
 905           call read_angles(inp,*38)
 906           goto 39
 907    38     write (iout,'(a)') 'Error reading reference structure.'
 908 #ifdef MPI
 909           call MPI_Finalize(MPI_COMM_WORLD,IERROR)
 910           stop 'Error reading reference structure'
 911 #endif
 912    39     call chainbuild
 913           call setup_var
 914 czscore          call geom_to_var(nvar,coord_exp_zs(1,1))
 915           nstart_sup=nnt
 916           nstart_seq=nnt
 917           nsup=nct-nnt+1
 918           kkk=1
 919           do i=1,2*nres
 920             do j=1,3
 921               cref(j,i,kkk)=c(j,i)
 922             enddo
 923           enddo
 924           call contact(.true.,ncont_ref,icont_ref,co)
 925         endif
 926 c        write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
 927         call flush(iout)
 928         if (constr_dist.gt.0) call read_dist_constr
 929         write (iout,*) "After read_dist_constr nhpb",nhpb
 930         call hpb_partition
 931         if(me.eq.king.or..not.out1file)
 932      &   write (iout,*) 'Contact order:',co
 933         if (pdbref) then
 934         if(me.eq.king.or..not.out1file)
 935      &   write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
 936         do i=1,ncont_ref
 937           do j=1,2
 938             icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
 939           enddo
 940           if(me.eq.king.or..not.out1file)
 941      &     write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
 942      &     icont_ref(1,i),' ',
 943      &     restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
 944         enddo
 945         endif
 946       endif
 947       if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
 948      &    .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
 949      &    modecalc.ne.10) then
 950 C If input structure hasn't been supplied from the PDB file read or generate
 951 C initial geometry.
 952         if (iranconf.eq.0 .and. .not. extconf) then
 953           if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
 954      &     write (iout,'(a)') 'Initial geometry will be read in.'
 955           if (read_cart) then
 956             read(inp,'(8f10.5)',end=36,err=36)
 957      &       ((c(l,k),l=1,3),k=1,nres),
 958      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
 959             write (iout,*) "Exit READ_CART"
 960             write (iout,'(8f10.5)')
 961      &       ((c(l,k),l=1,3),k=1,nres),
 962      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
 963             call int_from_cart1(.true.)
 964             write (iout,*) "Finish INT_TO_CART"
 965             do i=1,nres-1
 966               do j=1,3
 967                 dc(j,i)=c(j,i+1)-c(j,i)
 968                 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
 969               enddo
 970             enddo
 971             do i=nnt,nct
 972               if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
 973                 do j=1,3
 974                   dc(j,i+nres)=c(j,i+nres)-c(j,i)
 975                   dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
 976                 enddo
 977               endif
 978             enddo
 979             return
 980           else
 981             call read_angles(inp,*36)
 982           endif
 983           goto 37
 984    36     write (iout,'(a)') 'Error reading angle file.'
 985 #ifdef MPI
 986           call mpi_finalize( MPI_COMM_WORLD,IERR )
 987 #endif
 988           stop 'Error reading angle file.'
 989    37     continue
 990         else if (extconf) then
 991          if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
 992      &    write (iout,'(a)') 'Extended chain initial geometry.'
 993          do i=3,nres
 994           theta(i)=90d0*deg2rad
 995          enddo
 996          do i=4,nres
 997           phi(i)=180d0*deg2rad
 998          enddo
 999          do i=2,nres-1
1000           alph(i)=110d0*deg2rad
1001          enddo
1002          do i=2,nres-1
1003           omeg(i)=-120d0*deg2rad
1004           if (itype(i).le.0) omeg(i)=-omeg(i)
1005          enddo
1006         else
1007           if(me.eq.king.or..not.out1file)
1008      &     write (iout,'(a)') 'Random-generated initial geometry.'
1009
1010
1011 #ifdef MPI
1012           if (me.eq.king  .or. fg_rank.eq.0 .and. (
1013      &           modecalc.eq.12 .or. modecalc.eq.14) ) then
1014 #endif
1015             do itrial=1,100
1016               itmp=1
1017               call gen_rand_conf(itmp,*30)
1018               goto 40
1019    30         write (iout,*) 'Failed to generate random conformation',
1020      &          ', itrial=',itrial
1021               write (*,*) 'Processor:',me,
1022      &          ' Failed to generate random conformation',
1023      &          ' itrial=',itrial
1024               call intout
1025
1026 #ifdef AIX
1027               call flush_(iout)
1028 #else
1029               call flush(iout)
1030 #endif
1031             enddo
1032             write (iout,'(a,i3,a)') 'Processor:',me,
1033      &        ' error in generating random conformation.'
1034             write (*,'(a,i3,a)') 'Processor:',me,
1035      &        ' error in generating random conformation.'
1036             call flush(iout)
1037 #ifdef MPI
1038             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1039    40       continue
1040           endif
1041 #else
1042           do itrial=1,100
1043             itmp=1
1044             call gen_rand_conf(itmp,*30)
1045             goto 40
1046    30       write (iout,*) 'Failed to generate random conformation',
1047      &        ', itrial=',itrial
1048             write (*,*) 'Failed to generate random conformation',
1049      &        ', itrial=',itrial
1050           enddo
1051           write (iout,'(a,i3,a)') 'Processor:',me,
1052      &      ' error in generating random conformation.'
1053           write (*,'(a,i3,a)') 'Processor:',me,
1054      &      ' error in generating random conformation.'
1055           stop
1056    40     continue
1057 #endif
1058         endif
1059       elseif (modecalc.eq.4) then
1060         read (inp,'(a)') intinname
1061         open (intin,file=intinname,status='old',err=333)
1062         if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1063      &  write (iout,'(a)') 'intinname',intinname
1064         write (*,'(a)') 'Processor',myrank,' intinname',intinname
1065         goto 334
1066   333   write (iout,'(2a)') 'Error opening angle file ',intinname
1067 #ifdef MPI
1068         call MPI_Finalize(MPI_COMM_WORLD,IERR)
1069 #endif
1070         stop 'Error opening angle file.'
1071   334   continue
1072
1073       endif
1074 C Generate distance constraints, if the PDB structure is to be regularized.
1075       if (nthread.gt.0) then
1076         call read_threadbase
1077       endif
1078       call setup_var
1079       if (me.eq.king .or. .not. out1file)
1080      & call intout
1081       if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1082         write (iout,'(/a,i3,a)')
1083      &  'The chain contains',ns,' disulfide-bridging cysteines.'
1084         write (iout,'(20i4)') (iss(i),i=1,ns)
1085        if (dyn_ss) then
1086           write(iout,*)"Running with dynamic disulfide-bond formation"
1087        else
1088         write (iout,'(/a/)') 'Pre-formed links are:'
1089         do i=1,nss
1090           i1=ihpb(i)-nres
1091           i2=jhpb(i)-nres
1092           it1=itype(i1)
1093           it2=itype(i2)
1094           write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1095      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1096      &    ebr,forcon(i)
1097         enddo
1098         write (iout,'(a)')
1099        endif
1100       endif
1101       if (ns.gt.0.and.dyn_ss) then
1102           do i=nss+1,nhpb
1103             ihpb(i-nss)=ihpb(i)
1104             jhpb(i-nss)=jhpb(i)
1105             forcon(i-nss)=forcon(i)
1106             dhpb(i-nss)=dhpb(i)
1107           enddo
1108           nhpb=nhpb-nss
1109           nss=0
1110           call hpb_partition
1111           do i=1,ns
1112             dyn_ss_mask(iss(i))=.true.
1113           enddo
1114       endif
1115       if (i2ndstr.gt.0) call secstrp2dihc
1116 c      call geom_to_var(nvar,x)
1117 c      call etotal(energia(0))
1118 c      call enerprint(energia(0))
1119 c      call briefout(0,etot)
1120 c      stop
1121 cd    write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1122 cd    write (iout,'(a)') 'Variable list:'
1123 cd    write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1124 #ifdef MPI
1125       if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1126      &  write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1127      &  'Processor',myrank,': end reading molecular data.'
1128 #endif
1129       return
1130       end
1131 c--------------------------------------------------------------------------
1132       logical function seq_comp(itypea,itypeb,length)
1133       implicit none
1134       integer length,itypea(length),itypeb(length)
1135       integer i
1136       do i=1,length
1137         if (itypea(i).ne.itypeb(i)) then
1138           seq_comp=.false.
1139           return
1140         endif
1141       enddo
1142       seq_comp=.true.
1143       return
1144       end
1145 c-----------------------------------------------------------------------------
1146       subroutine read_bridge
1147 C Read information about disulfide bridges.
1148       implicit real*8 (a-h,o-z)
1149       include 'DIMENSIONS'
1150 #ifdef MPI
1151       include 'mpif.h'
1152 #endif
1153       include 'COMMON.IOUNITS'
1154       include 'COMMON.GEO'
1155       include 'COMMON.VAR'
1156       include 'COMMON.INTERACT'
1157       include 'COMMON.LOCAL'
1158       include 'COMMON.NAMES'
1159       include 'COMMON.CHAIN'
1160       include 'COMMON.FFIELD'
1161       include 'COMMON.SBRIDGE'
1162       include 'COMMON.HEADER'
1163       include 'COMMON.CONTROL'
1164       include 'COMMON.DBASE'
1165       include 'COMMON.THREAD'
1166       include 'COMMON.TIME1'
1167       include 'COMMON.SETUP'
1168 C Read bridging residues.
1169       read (inp,*) ns,(iss(i),i=1,ns)
1170       print *,'ns=',ns
1171       if(me.eq.king.or..not.out1file)
1172      &  write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1173 C Check whether the specified bridging residues are cystines.
1174       do i=1,ns
1175         if (itype(iss(i)).ne.1) then
1176           if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1177      &   'Do you REALLY think that the residue ',
1178      &    restyp(itype(iss(i))),i,
1179      &   ' can form a disulfide bridge?!!!'
1180           write (*,'(2a,i3,a)')
1181      &   'Do you REALLY think that the residue ',
1182      &    restyp(itype(iss(i))),i,
1183      &   ' can form a disulfide bridge?!!!'
1184 #ifdef MPI
1185          call MPI_Finalize(MPI_COMM_WORLD,ierror)
1186          stop
1187 #endif
1188         endif
1189       enddo
1190 C Read preformed bridges.
1191       if (ns.gt.0) then
1192       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1193       if(fg_rank.eq.0)
1194      & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1195       if (nss.gt.0) then
1196         nhpb=nss
1197 C Check if the residues involved in bridges are in the specified list of
1198 C bridging residues.
1199         do i=1,nss
1200           do j=1,i-1
1201             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1202      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1203               write (iout,'(a,i3,a)') 'Disulfide pair',i,
1204      &      ' contains residues present in other pairs.'
1205               write (*,'(a,i3,a)') 'Disulfide pair',i,
1206      &      ' contains residues present in other pairs.'
1207 #ifdef MPI
1208               call MPI_Finalize(MPI_COMM_WORLD,ierror)
1209               stop
1210 #endif
1211             endif
1212           enddo
1213           do j=1,ns
1214             if (ihpb(i).eq.iss(j)) goto 10
1215           enddo
1216           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1217    10     continue
1218           do j=1,ns
1219             if (jhpb(i).eq.iss(j)) goto 20
1220           enddo
1221           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1222    20     continue
1223           dhpb(i)=dbr
1224           forcon(i)=fbr
1225         enddo
1226         do i=1,nss
1227           ihpb(i)=ihpb(i)+nres
1228           jhpb(i)=jhpb(i)+nres
1229         enddo
1230       endif
1231       endif
1232       return
1233       end
1234 c----------------------------------------------------------------------------
1235       subroutine read_x(kanal,*)
1236       implicit real*8 (a-h,o-z)
1237       include 'DIMENSIONS'
1238       include 'COMMON.GEO'
1239       include 'COMMON.VAR'
1240       include 'COMMON.CHAIN'
1241       include 'COMMON.IOUNITS'
1242       include 'COMMON.CONTROL'
1243       include 'COMMON.LOCAL'
1244       include 'COMMON.INTERACT'
1245 c Read coordinates from input
1246 c
1247       read(kanal,'(8f10.5)',end=10,err=10)
1248      &  ((c(l,k),l=1,3),k=1,nres),
1249      &  ((c(l,k+nres),l=1,3),k=nnt,nct)
1250       do j=1,3
1251         c(j,nres+1)=c(j,1)
1252         c(j,2*nres)=c(j,nres)
1253       enddo
1254       call int_from_cart1(.false.)
1255       do i=1,nres-1
1256         do j=1,3
1257           dc(j,i)=c(j,i+1)-c(j,i)
1258           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1259         enddo
1260       enddo
1261       do i=nnt,nct
1262         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1263           do j=1,3
1264             dc(j,i+nres)=c(j,i+nres)-c(j,i)
1265             dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1266           enddo
1267         endif
1268       enddo
1269
1270       return
1271    10 return1
1272       end
1273 c----------------------------------------------------------------------------
1274       subroutine read_threadbase
1275       implicit real*8 (a-h,o-z)
1276       include 'DIMENSIONS'
1277       include 'COMMON.IOUNITS'
1278       include 'COMMON.GEO'
1279       include 'COMMON.VAR'
1280       include 'COMMON.INTERACT'
1281       include 'COMMON.LOCAL'
1282       include 'COMMON.NAMES'
1283       include 'COMMON.CHAIN'
1284       include 'COMMON.FFIELD'
1285       include 'COMMON.SBRIDGE'
1286       include 'COMMON.HEADER'
1287       include 'COMMON.CONTROL'
1288       include 'COMMON.DBASE'
1289       include 'COMMON.THREAD'
1290       include 'COMMON.TIME1'
1291 C Read pattern database for threading.
1292       read (icbase,*) nseq
1293       do i=1,nseq
1294         read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1295      &   nres_base(2,i),nres_base(3,i)
1296         read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1297      &   nres_base(1,i))
1298 c       write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1299 c    &   nres_base(2,i),nres_base(3,i)
1300 c       write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1301 c    &   nres_base(1,i))
1302       enddo
1303       close (icbase)
1304       if (weidis.eq.0.0D0) weidis=0.1D0
1305       do i=nnt,nct
1306         do j=i+2,nct
1307           nhpb=nhpb+1
1308           ihpb(nhpb)=i
1309           jhpb(nhpb)=j
1310           forcon(nhpb)=weidis
1311         enddo
1312       enddo
1313       read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1314       write (iout,'(a,i5)') 'nexcl: ',nexcl
1315       write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1316       return
1317       end
1318 c------------------------------------------------------------------------------
1319       subroutine setup_var
1320       implicit real*8 (a-h,o-z)
1321       include 'DIMENSIONS'
1322       include 'COMMON.IOUNITS'
1323       include 'COMMON.GEO'
1324       include 'COMMON.VAR'
1325       include 'COMMON.INTERACT'
1326       include 'COMMON.LOCAL'
1327       include 'COMMON.NAMES'
1328       include 'COMMON.CHAIN'
1329       include 'COMMON.FFIELD'
1330       include 'COMMON.SBRIDGE'
1331       include 'COMMON.HEADER'
1332       include 'COMMON.CONTROL'
1333       include 'COMMON.DBASE'
1334       include 'COMMON.THREAD'
1335       include 'COMMON.TIME1'
1336 C Set up variable list.
1337       ntheta=nres-2
1338       nphi=nres-3
1339       nvar=ntheta+nphi
1340       nside=0
1341       do i=2,nres-1
1342         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1343           nside=nside+1
1344           ialph(i,1)=nvar+nside
1345           ialph(nside,2)=i
1346         endif
1347       enddo
1348       if (indphi.gt.0) then
1349         nvar=nphi
1350       else if (indback.gt.0) then
1351         nvar=nphi+ntheta
1352       else
1353         nvar=nvar+2*nside
1354       endif
1355 cd    write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1356       return
1357       end
1358 c----------------------------------------------------------------------------
1359       subroutine gen_dist_constr
1360 C Generate CA distance constraints.
1361       implicit real*8 (a-h,o-z)
1362       include 'DIMENSIONS'
1363       include 'COMMON.IOUNITS'
1364       include 'COMMON.GEO'
1365       include 'COMMON.VAR'
1366       include 'COMMON.INTERACT'
1367       include 'COMMON.LOCAL'
1368       include 'COMMON.NAMES'
1369       include 'COMMON.CHAIN'
1370       include 'COMMON.FFIELD'
1371       include 'COMMON.SBRIDGE'
1372       include 'COMMON.HEADER'
1373       include 'COMMON.CONTROL'
1374       include 'COMMON.DBASE'
1375       include 'COMMON.THREAD'
1376       include 'COMMON.TIME1'
1377       dimension itype_pdb(maxres)
1378       common /pizda/ itype_pdb
1379       character*2 iden
1380 cd      print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1381 cd      write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1382 cd     & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1383 cd     & ' nsup',nsup
1384       do i=nstart_sup,nstart_sup+nsup-1
1385 cd      write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1386 cd     &    ' seq_pdb', restyp(itype_pdb(i))
1387         do j=i+2,nstart_sup+nsup-1
1388           nhpb=nhpb+1
1389           ihpb(nhpb)=i+nstart_seq-nstart_sup
1390           jhpb(nhpb)=j+nstart_seq-nstart_sup
1391           forcon(nhpb)=weidis
1392           dhpb(nhpb)=dist(i,j)
1393         enddo
1394       enddo
1395 cd      write (iout,'(a)') 'Distance constraints:'
1396 cd      do i=nss+1,nhpb
1397 cd        ii=ihpb(i)
1398 cd        jj=jhpb(i)
1399 cd        iden='CA'
1400 cd        if (ii.gt.nres) then
1401 cd          iden='SC'
1402 cd          ii=ii-nres
1403 cd          jj=jj-nres
1404 cd        endif
1405 cd        write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1406 cd     &  restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1407 cd     &  dhpb(i),forcon(i)
1408 cd      enddo
1409       return
1410       end
1411 c----------------------------------------------------------------------------
1412       subroutine map_read
1413       implicit real*8 (a-h,o-z)
1414       include 'DIMENSIONS'
1415       include 'COMMON.MAP'
1416       include 'COMMON.IOUNITS'
1417       character*3 angid(4) /'THE','PHI','ALP','OME'/
1418       character*80 mapcard,ucase
1419       do imap=1,nmap
1420         read (inp,'(a)') mapcard
1421         mapcard=ucase(mapcard)
1422         if (index(mapcard,'PHI').gt.0) then
1423           kang(imap)=1
1424         else if (index(mapcard,'THE').gt.0) then
1425           kang(imap)=2
1426         else if (index(mapcard,'ALP').gt.0) then
1427           kang(imap)=3
1428         else if (index(mapcard,'OME').gt.0) then
1429           kang(imap)=4
1430         else
1431           write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1432           stop 'Error - illegal variable spec in MAP card.'
1433         endif
1434         call readi (mapcard,'RES1',res1(imap),0)
1435         call readi (mapcard,'RES2',res2(imap),0)
1436         if (res1(imap).eq.0) then
1437           res1(imap)=res2(imap)
1438         else if (res2(imap).eq.0) then
1439           res2(imap)=res1(imap)
1440         endif
1441         if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1442           write (iout,'(a)')
1443      &    'Error - illegal definition of variable group in MAP.'
1444           stop 'Error - illegal definition of variable group in MAP.'
1445         endif
1446         call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1447         call reada(mapcard,'TO',ang_to(imap),0.0D0)
1448         call readi(mapcard,'NSTEP',nstep(imap),0)
1449         if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1450           write (iout,'(a)')
1451      &     'Illegal boundary and/or step size specification in MAP.'
1452           stop 'Illegal boundary and/or step size specification in MAP.'
1453         endif
1454       enddo ! imap
1455       return
1456       end
1457 c----------------------------------------------------------------------------
1458       subroutine csaread
1459       implicit real*8 (a-h,o-z)
1460       include 'DIMENSIONS'
1461       include 'COMMON.IOUNITS'
1462       include 'COMMON.GEO'
1463       include 'COMMON.CSA'
1464       include 'COMMON.BANK'
1465       include 'COMMON.CONTROL'
1466       character*80 ucase
1467       character*620 mcmcard
1468       call card_concat(mcmcard)
1469
1470       call readi(mcmcard,'NCONF',nconf,50)
1471       call readi(mcmcard,'NADD',nadd,0)
1472       call readi(mcmcard,'JSTART',jstart,1)
1473       call readi(mcmcard,'JEND',jend,1)
1474       call readi(mcmcard,'NSTMAX',nstmax,500000)
1475       call readi(mcmcard,'N0',n0,1)
1476       call readi(mcmcard,'N1',n1,6)
1477       call readi(mcmcard,'N2',n2,4)
1478       call readi(mcmcard,'N3',n3,0)
1479       call readi(mcmcard,'N4',n4,0)
1480       call readi(mcmcard,'N5',n5,0)
1481       call readi(mcmcard,'N6',n6,10)
1482       call readi(mcmcard,'N7',n7,0)
1483       call readi(mcmcard,'N8',n8,0)
1484       call readi(mcmcard,'N9',n9,0)
1485       call readi(mcmcard,'N14',n14,0)
1486       call readi(mcmcard,'N15',n15,0)
1487       call readi(mcmcard,'N16',n16,0)
1488       call readi(mcmcard,'N17',n17,0)
1489       call readi(mcmcard,'N18',n18,0)
1490
1491       vdisulf=(index(mcmcard,'DYNSS').gt.0)
1492
1493       call readi(mcmcard,'NDIFF',ndiff,2)
1494       call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1495       call readi(mcmcard,'IS1',is1,1)
1496       call readi(mcmcard,'IS2',is2,8)
1497       call readi(mcmcard,'NRAN0',nran0,4)
1498       call readi(mcmcard,'NRAN1',nran1,2)
1499       call readi(mcmcard,'IRR',irr,1)
1500       call readi(mcmcard,'NSEED',nseed,20)
1501       call readi(mcmcard,'NTOTAL',ntotal,10000)
1502       call reada(mcmcard,'CUT1',cut1,2.0d0)
1503       call reada(mcmcard,'CUT2',cut2,5.0d0)
1504       call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1505       call readi(mcmcard,'ICMAX',icmax,3)
1506       call readi(mcmcard,'IRESTART',irestart,0)
1507 c!bankt      call readi(mcmcard,'NBANKTM',ntbankm,0)
1508       ntbankm=0
1509 c!bankt
1510       call reada(mcmcard,'DELE',dele,20.0d0)
1511       call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1512       call readi(mcmcard,'IREF',iref,0)
1513       call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1514       call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1515       call readi(mcmcard,'NCONF_IN',nconf_in,0)
1516       call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1517       write (iout,*) "NCONF_IN",nconf_in
1518       return
1519       end
1520 c----------------------------------------------------------------------------
1521 cfmc      subroutine mcmfread
1522 cfmc      implicit real*8 (a-h,o-z)
1523 cfmc      include 'DIMENSIONS'
1524 cfmc      include 'COMMON.MCMF'
1525 cfmc      include 'COMMON.IOUNITS'
1526 cfmc      include 'COMMON.GEO'
1527 cfmc      character*80 ucase
1528 cfmc      character*620 mcmcard
1529 cfmc      call card_concat(mcmcard)
1530 cfmc
1531 cfmc      call readi(mcmcard,'MAXRANT',maxrant,1000)
1532 cfmc      write(iout,*)'MAXRANT=',maxrant
1533 cfmc      call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1534 cfmc      write(iout,*)'MAXFAM=',maxfam
1535 cfmc      call readi(mcmcard,'NNET1',nnet1,5)
1536 cfmc      write(iout,*)'NNET1=',nnet1
1537 cfmc      call readi(mcmcard,'NNET2',nnet2,4)
1538 cfmc      write(iout,*)'NNET2=',nnet2
1539 cfmc      call readi(mcmcard,'NNET3',nnet3,4)
1540 cfmc      write(iout,*)'NNET3=',nnet3
1541 cfmc      call readi(mcmcard,'ILASTT',ilastt,0)
1542 cfmc      write(iout,*)'ILASTT=',ilastt
1543 cfmc      call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1544 cfmc      write(iout,*)'MAXSTR=',maxstr
1545 cfmc      maxstr_f=maxstr/maxfam
1546 cfmc      write(iout,*)'MAXSTR_F=',maxstr_f
1547 cfmc      call readi(mcmcard,'NMCMF',nmcmf,10)
1548 cfmc      write(iout,*)'NMCMF=',nmcmf
1549 cfmc      call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1550 cfmc      write(iout,*)'IFOCUS=',ifocus
1551 cfmc      call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1552 cfmc      write(iout,*)'NLOCMCMF=',nlocmcmf
1553 cfmc      call readi(mcmcard,'INTPRT',intprt,1000)
1554 cfmc      write(iout,*)'INTPRT=',intprt
1555 cfmc      call readi(mcmcard,'IPRT',iprt,100)
1556 cfmc      write(iout,*)'IPRT=',iprt
1557 cfmc      call readi(mcmcard,'IMAXTR',imaxtr,100)
1558 cfmc      write(iout,*)'IMAXTR=',imaxtr
1559 cfmc      call readi(mcmcard,'MAXEVEN',maxeven,1000)
1560 cfmc      write(iout,*)'MAXEVEN=',maxeven
1561 cfmc      call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1562 cfmc      write(iout,*)'MAXEVEN1=',maxeven1
1563 cfmc      call readi(mcmcard,'INIMIN',inimin,200)
1564 cfmc      write(iout,*)'INIMIN=',inimin
1565 cfmc      call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1566 cfmc      write(iout,*)'NSTEPMCMF=',nstepmcmf
1567 cfmc      call readi(mcmcard,'NTHREAD',nthread,5)
1568 cfmc      write(iout,*)'NTHREAD=',nthread
1569 cfmc      call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1570 cfmc      write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1571 cfmc      call readi(mcmcard,'MAXPERT',maxpert,9)
1572 cfmc      write(iout,*)'MAXPERT=',maxpert
1573 cfmc      call readi(mcmcard,'IRMSD',irmsd,1)
1574 cfmc      write(iout,*)'IRMSD=',irmsd
1575 cfmc      call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1576 cfmc      write(iout,*)'DENEMIN=',denemin
1577 cfmc      call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1578 cfmc      write(iout,*)'RCUT1S=',rcut1s
1579 cfmc      call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1580 cfmc      write(iout,*)'RCUT1E=',rcut1e
1581 cfmc      call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1582 cfmc      write(iout,*)'RCUT2S=',rcut2s
1583 cfmc      call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1584 cfmc      write(iout,*)'RCUT2E=',rcut2e
1585 cfmc      call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1586 cfmc      write(iout,*)'DPERT1=',d_pert1
1587 cfmc      call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1588 cfmc      write(iout,*)'DPERT1A=',d_pert1a
1589 cfmc      call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1590 cfmc      write(iout,*)'DPERT2=',d_pert2
1591 cfmc      call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1592 cfmc      write(iout,*)'DPERT2A=',d_pert2a
1593 cfmc      call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1594 cfmc      write(iout,*)'DPERT2B=',d_pert2b
1595 cfmc      call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1596 cfmc      write(iout,*)'DPERT2C=',d_pert2c
1597 cfmc      d_pert1=deg2rad*d_pert1
1598 cfmc      d_pert1a=deg2rad*d_pert1a
1599 cfmc      d_pert2=deg2rad*d_pert2
1600 cfmc      d_pert2a=deg2rad*d_pert2a
1601 cfmc      d_pert2b=deg2rad*d_pert2b
1602 cfmc      d_pert2c=deg2rad*d_pert2c
1603 cfmc      call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1604 cfmc      write(iout,*)'KT_MCMF1=',kt_mcmf1
1605 cfmc      call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1606 cfmc      write(iout,*)'KT_MCMF2=',kt_mcmf2
1607 cfmc      call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1608 cfmc      write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1609 cfmc      call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1610 cfmc      write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1611 cfmc      call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1612 cfmc      write(iout,*)'RCUTINI=',rcutini
1613 cfmc      call reada(mcmcard,'GRAT',grat,0.5D0)
1614 cfmc      write(iout,*)'GRAT=',grat
1615 cfmc      call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1616 cfmc      write(iout,*)'BIAS_MCMF=',bias_mcmf
1617 cfmc
1618 cfmc      return
1619 cfmc      end
1620 c----------------------------------------------------------------------------
1621       subroutine mcmread
1622       implicit real*8 (a-h,o-z)
1623       include 'DIMENSIONS'
1624       include 'COMMON.MCM'
1625       include 'COMMON.MCE'
1626       include 'COMMON.IOUNITS'
1627       character*80 ucase
1628       character*320 mcmcard
1629       call card_concat(mcmcard)
1630       call readi(mcmcard,'MAXACC',maxacc,100)
1631       call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1632       call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1633       call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1634       call readi(mcmcard,'MAXREPM',maxrepm,200)
1635       call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1636       call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1637       call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1638       call reada(mcmcard,'E_UP',e_up,5.0D0)
1639       call reada(mcmcard,'DELTE',delte,0.1D0)
1640       call readi(mcmcard,'NSWEEP',nsweep,5)
1641       call readi(mcmcard,'NSTEPH',nsteph,0)
1642       call readi(mcmcard,'NSTEPC',nstepc,0)
1643       call reada(mcmcard,'TMIN',tmin,298.0D0)
1644       call reada(mcmcard,'TMAX',tmax,298.0D0)
1645       call readi(mcmcard,'NWINDOW',nwindow,0)
1646       call readi(mcmcard,'PRINT_MC',print_mc,0)
1647       print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1648       print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1649       ent_read=(index(mcmcard,'ENT_READ').gt.0)
1650       call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1651       call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1652       call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1653       call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1654       call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1655       if (nwindow.gt.0) then
1656         read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1657         do i=1,nwindow
1658           winlen(i)=winend(i)-winstart(i)+1
1659         enddo
1660       endif
1661       if (tmax.lt.tmin) tmax=tmin
1662       if (tmax.eq.tmin) then
1663         nstepc=0
1664         nsteph=0
1665       endif
1666       if (nstepc.gt.0 .and. nsteph.gt.0) then
1667         tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1668         tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1669       endif
1670 C Probabilities of different move types
1671       sumpro_type(0)=0.0D0
1672       call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1673       call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1674       sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1675       call reada(mcmcard,'THETA'     ,sumpro_type(3),0.0d0)
1676       sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1677       call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1678       sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1679       do i=1,MaxMoveType
1680         print *,'i',i,' sumprotype',sumpro_type(i)
1681         sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1682         print *,'i',i,' sumprotype',sumpro_type(i)
1683       enddo
1684       return
1685       end
1686 c----------------------------------------------------------------------------
1687       subroutine read_minim
1688       implicit real*8 (a-h,o-z)
1689       include 'DIMENSIONS'
1690       include 'COMMON.MINIM'
1691       include 'COMMON.IOUNITS'
1692       character*80 ucase
1693       character*320 minimcard
1694       call card_concat(minimcard)
1695       call readi(minimcard,'MAXMIN',maxmin,2000)
1696       call readi(minimcard,'MAXFUN',maxfun,5000)
1697       call readi(minimcard,'MINMIN',minmin,maxmin)
1698       call readi(minimcard,'MINFUN',minfun,maxmin)
1699       call reada(minimcard,'TOLF',tolf,1.0D-2)
1700       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1701       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1702       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1703       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1704       write (iout,'(/80(1h*)/20x,a/80(1h*))')
1705      &         'Options in energy minimization:'
1706       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1707      & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1708      & 'MinMin:',MinMin,' MinFun:',MinFun,
1709      & ' TolF:',TolF,' RTolF:',RTolF
1710       return
1711       end
1712 c----------------------------------------------------------------------------
1713       subroutine read_angles(kanal,*)
1714       implicit real*8 (a-h,o-z)
1715       include 'DIMENSIONS'
1716       include 'COMMON.GEO'
1717       include 'COMMON.VAR'
1718       include 'COMMON.CHAIN'
1719       include 'COMMON.IOUNITS'
1720       include 'COMMON.CONTROL'
1721 c Read angles from input
1722 c
1723        read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1724        read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1725        read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1726        read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1727
1728        do i=1,nres
1729 c 9/7/01 avoid 180 deg valence angle
1730         if (theta(i).gt.179.99d0) theta(i)=179.99d0
1731 c
1732         theta(i)=deg2rad*theta(i)
1733         phi(i)=deg2rad*phi(i)
1734         alph(i)=deg2rad*alph(i)
1735         omeg(i)=deg2rad*omeg(i)
1736        enddo
1737       return
1738    10 return1
1739       end
1740 c----------------------------------------------------------------------------
1741       subroutine reada(rekord,lancuch,wartosc,default)
1742       implicit none
1743       character*(*) rekord,lancuch
1744       double precision wartosc,default
1745       integer ilen,iread
1746       external ilen
1747       iread=index(rekord,lancuch)
1748       if (iread.eq.0) then
1749         wartosc=default
1750         return
1751       endif
1752       iread=iread+ilen(lancuch)+1
1753       read (rekord(iread:),*,err=10,end=10) wartosc
1754       return
1755   10  wartosc=default
1756       return
1757       end
1758 c----------------------------------------------------------------------------
1759       subroutine readi(rekord,lancuch,wartosc,default)
1760       implicit none
1761       character*(*) rekord,lancuch
1762       integer wartosc,default
1763       integer ilen,iread
1764       external ilen
1765       iread=index(rekord,lancuch)
1766       if (iread.eq.0) then
1767         wartosc=default
1768         return
1769       endif
1770       iread=iread+ilen(lancuch)+1
1771       read (rekord(iread:),*,err=10,end=10) wartosc
1772       return
1773   10  wartosc=default
1774       return
1775       end
1776 c----------------------------------------------------------------------------
1777       subroutine multreadi(rekord,lancuch,tablica,dim,default)
1778       implicit none
1779       integer dim,i
1780       integer tablica(dim),default
1781       character*(*) rekord,lancuch
1782       character*80 aux
1783       integer ilen,iread
1784       external ilen
1785       do i=1,dim
1786         tablica(i)=default
1787       enddo
1788       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1789       if (iread.eq.0) return
1790       iread=iread+ilen(lancuch)+1
1791       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1792    10 return
1793       end
1794 c----------------------------------------------------------------------------
1795       subroutine multreada(rekord,lancuch,tablica,dim,default)
1796       implicit none
1797       integer dim,i
1798       double precision tablica(dim),default
1799       character*(*) rekord,lancuch
1800       character*80 aux
1801       integer ilen,iread
1802       external ilen
1803       do i=1,dim
1804         tablica(i)=default
1805       enddo
1806       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1807       if (iread.eq.0) return
1808       iread=iread+ilen(lancuch)+1
1809       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1810    10 return
1811       end
1812 c----------------------------------------------------------------------------
1813       subroutine openunits
1814       implicit real*8 (a-h,o-z)
1815       include 'DIMENSIONS'
1816 #ifdef MPI
1817       include 'mpif.h'
1818       character*16 form,nodename
1819       integer nodelen
1820 #endif
1821       include 'COMMON.SETUP'
1822       include 'COMMON.IOUNITS'
1823       include 'COMMON.MD'
1824       include 'COMMON.CONTROL'
1825       integer lenpre,lenpot,ilen,lentmp
1826       external ilen
1827       character*3 out1file_text,ucase
1828       character*3 ll
1829       external ucase
1830 c      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1831       call getenv_loc("PREFIX",prefix)
1832       pref_orig = prefix
1833       call getenv_loc("POT",pot)
1834       call getenv_loc("DIRTMP",tmpdir)
1835       call getenv_loc("CURDIR",curdir)
1836       call getenv_loc("OUT1FILE",out1file_text)
1837 c      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1838       out1file_text=ucase(out1file_text)
1839       if (out1file_text(1:1).eq."Y") then
1840         out1file=.true.
1841       else
1842         out1file=fg_rank.gt.0
1843       endif
1844       lenpre=ilen(prefix)
1845       lenpot=ilen(pot)
1846       lentmp=ilen(tmpdir)
1847       if (lentmp.gt.0) then
1848           write (*,'(80(1h!))')
1849           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
1850           write (*,'(80(1h!))')
1851           write (*,*)"All output files will be on node /tmp directory."
1852 #ifdef MPI
1853         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1854         if (me.eq.king) then
1855           write (*,*) "The master node is ",nodename
1856         else if (fg_rank.eq.0) then
1857           write (*,*) "I am the CG slave node ",nodename
1858         else
1859           write (*,*) "I am the FG slave node ",nodename
1860         endif
1861 #endif
1862         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1863         lenpre = lentmp+lenpre+1
1864       endif
1865       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1866 C Get the names and open the input files
1867 #if defined(WINIFL) || defined(WINPGI)
1868       open(1,file=pref_orig(:ilen(pref_orig))//
1869      &  '.inp',status='old',readonly,shared)
1870        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1871 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1872 C Get parameter filenames and open the parameter files.
1873       call getenv_loc('BONDPAR',bondname)
1874       open (ibond,file=bondname,status='old',readonly,shared)
1875       call getenv_loc('THETPAR',thetname)
1876       open (ithep,file=thetname,status='old',readonly,shared)
1877       call getenv_loc('ROTPAR',rotname)
1878       open (irotam,file=rotname,status='old',readonly,shared)
1879       call getenv_loc('TORPAR',torname)
1880       open (itorp,file=torname,status='old',readonly,shared)
1881       call getenv_loc('TORDPAR',tordname)
1882       open (itordp,file=tordname,status='old',readonly,shared)
1883       call getenv_loc('FOURIER',fouriername)
1884       open (ifourier,file=fouriername,status='old',readonly,shared)
1885       call getenv_loc('ELEPAR',elename)
1886       open (ielep,file=elename,status='old',readonly,shared)
1887       call getenv_loc('SIDEPAR',sidename)
1888       open (isidep,file=sidename,status='old',readonly,shared)
1889 #elif (defined CRAY) || (defined AIX)
1890       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1891      &  action='read')
1892 c      print *,"Processor",myrank," opened file 1"
1893       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1894 c      print *,"Processor",myrank," opened file 9"
1895 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1896 C Get parameter filenames and open the parameter files.
1897       call getenv_loc('BONDPAR',bondname)
1898       open (ibond,file=bondname,status='old',action='read')
1899 c      print *,"Processor",myrank," opened file IBOND"
1900       call getenv_loc('THETPAR',thetname)
1901       open (ithep,file=thetname,status='old',action='read')
1902 c      print *,"Processor",myrank," opened file ITHEP"
1903       call getenv_loc('ROTPAR',rotname)
1904       open (irotam,file=rotname,status='old',action='read')
1905 c      print *,"Processor",myrank," opened file IROTAM"
1906       call getenv_loc('TORPAR',torname)
1907       open (itorp,file=torname,status='old',action='read')
1908 c      print *,"Processor",myrank," opened file ITORP"
1909       call getenv_loc('TORDPAR',tordname)
1910       open (itordp,file=tordname,status='old',action='read')
1911 c      print *,"Processor",myrank," opened file ITORDP"
1912       call getenv_loc('SCCORPAR',sccorname)
1913       open (isccor,file=sccorname,status='old',action='read')
1914 c      print *,"Processor",myrank," opened file ISCCOR"
1915       call getenv_loc('FOURIER',fouriername)
1916       open (ifourier,file=fouriername,status='old',action='read')
1917 c      print *,"Processor",myrank," opened file IFOURIER"
1918       call getenv_loc('ELEPAR',elename)
1919       open (ielep,file=elename,status='old',action='read')
1920 c      print *,"Processor",myrank," opened file IELEP"
1921       call getenv_loc('SIDEPAR',sidename)
1922       open (isidep,file=sidename,status='old',action='read')
1923 c      print *,"Processor",myrank," opened file ISIDEP"
1924 c      print *,"Processor",myrank," opened parameter files"
1925 #elif (defined G77)
1926       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
1927       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1928 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1929 C Get parameter filenames and open the parameter files.
1930       call getenv_loc('BONDPAR',bondname)
1931       open (ibond,file=bondname,status='old')
1932       call getenv_loc('THETPAR',thetname)
1933       open (ithep,file=thetname,status='old')
1934       call getenv_loc('ROTPAR',rotname)
1935       open (irotam,file=rotname,status='old')
1936       call getenv_loc('TORPAR',torname)
1937       open (itorp,file=torname,status='old')
1938       call getenv_loc('TORDPAR',tordname)
1939       open (itordp,file=tordname,status='old')
1940       call getenv_loc('SCCORPAR',sccorname)
1941       open (isccor,file=sccorname,status='old')
1942       call getenv_loc('FOURIER',fouriername)
1943       open (ifourier,file=fouriername,status='old')
1944       call getenv_loc('ELEPAR',elename)
1945       open (ielep,file=elename,status='old')
1946       call getenv_loc('SIDEPAR',sidename)
1947       open (isidep,file=sidename,status='old')
1948 #else
1949       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1950      &  readonly)
1951        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1952 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1953 C Get parameter filenames and open the parameter files.
1954       call getenv_loc('BONDPAR',bondname)
1955       open (ibond,file=bondname,status='old',readonly)
1956       call getenv_loc('THETPAR',thetname)
1957       open (ithep,file=thetname,status='old',readonly)
1958       call getenv_loc('ROTPAR',rotname)
1959       open (irotam,file=rotname,status='old',readonly)
1960       call getenv_loc('TORPAR',torname)
1961       open (itorp,file=torname,status='old',readonly)
1962       call getenv_loc('TORDPAR',tordname)
1963       open (itordp,file=tordname,status='old',readonly)
1964       call getenv_loc('SCCORPAR',sccorname)
1965       open (isccor,file=sccorname,status='old',readonly)
1966 #ifndef CRYST_THETA
1967       call getenv_loc('THETPARPDB',thetname_pdb)
1968       print *,"thetname_pdb ",thetname_pdb
1969       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
1970       print *,ithep_pdb," opened"
1971 #endif
1972       call getenv_loc('FOURIER',fouriername)
1973       open (ifourier,file=fouriername,status='old',readonly)
1974       call getenv_loc('ELEPAR',elename)
1975       open (ielep,file=elename,status='old',readonly)
1976       call getenv_loc('SIDEPAR',sidename)
1977       open (isidep,file=sidename,status='old',readonly)
1978 #ifndef CRYST_SC
1979       call getenv_loc('ROTPARPDB',rotname_pdb)
1980       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
1981 #endif
1982 #endif
1983 #ifndef OLDSCP
1984 C
1985 C 8/9/01 In the newest version SCp interaction constants are read from a file
1986 C Use -DOLDSCP to use hard-coded constants instead.
1987 C
1988       call getenv_loc('SCPPAR',scpname)
1989 #if defined(WINIFL) || defined(WINPGI)
1990       open (iscpp,file=scpname,status='old',readonly,shared)
1991 #elif (defined CRAY)  || (defined AIX)
1992       open (iscpp,file=scpname,status='old',action='read')
1993 #elif (defined G77)
1994       open (iscpp,file=scpname,status='old')
1995 #else
1996       open (iscpp,file=scpname,status='old',readonly)
1997 #endif
1998 #endif
1999       call getenv_loc('PATTERN',patname)
2000 #if defined(WINIFL) || defined(WINPGI)
2001       open (icbase,file=patname,status='old',readonly,shared)
2002 #elif (defined CRAY)  || (defined AIX)
2003       open (icbase,file=patname,status='old',action='read')
2004 #elif (defined G77)
2005       open (icbase,file=patname,status='old')
2006 #else
2007       open (icbase,file=patname,status='old',readonly)
2008 #endif
2009 #ifdef MPI
2010 C Open output file only for CG processes
2011 c      print *,"Processor",myrank," fg_rank",fg_rank
2012       if (fg_rank.eq.0) then
2013
2014       if (nodes.eq.1) then
2015         npos=3
2016       else
2017         npos = dlog10(dfloat(nodes-1))+1
2018       endif
2019       if (npos.lt.3) npos=3
2020       write (liczba,'(i1)') npos
2021       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2022      &  //')'
2023       write (liczba,form) me
2024       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2025      &  liczba(:ilen(liczba))
2026       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2027      &  //'.int'
2028       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2029      &  //'.pdb'
2030       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2031      &  liczba(:ilen(liczba))//'.mol2'
2032       statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2033      &  liczba(:ilen(liczba))//'.stat'
2034       if (lentmp.gt.0)
2035      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2036      &      //liczba(:ilen(liczba))//'.stat')
2037       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2038      &  //'.rst'
2039       if(usampl) then
2040           qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2041      & liczba(:ilen(liczba))//'.const'
2042       endif
2043
2044       endif
2045 #else
2046       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2047       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2048       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2049       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2050       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2051       if (lentmp.gt.0)
2052      &  call copy_to_tmp(pref_orig(:ile(pref_orig))//'_'//pot(:lenpot)//
2053      &    //'.stat')
2054       rest2name=prefix(:ilen(prefix))//'.rst'
2055       if(usampl) then
2056          qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2057       endif
2058 #endif
2059 #if defined(AIX) || defined(PGI)
2060       if (me.eq.king .or. .not. out1file)
2061      &   open(iout,file=outname,status='unknown')
2062 #ifdef DEBUG
2063       if (fg_rank.gt.0) then
2064         write (liczba,'(i3.3)') myrank/nfgtasks
2065         write (ll,'(bz,i3.3)') fg_rank
2066         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2067      &   status='unknown')
2068       endif
2069 #endif
2070       if(me.eq.king) then
2071        open(igeom,file=intname,status='unknown',position='append')
2072        open(ipdb,file=pdbname,status='unknown')
2073        open(imol2,file=mol2name,status='unknown')
2074        open(istat,file=statname,status='unknown',position='append')
2075       else
2076 c1out       open(iout,file=outname,status='unknown')
2077       endif
2078 #else
2079       if (me.eq.king .or. .not.out1file)
2080      &    open(iout,file=outname,status='unknown')
2081 #ifdef DEBUG
2082       if (fg_rank.gt.0) then
2083         write (liczba,'(i3.3)') myrank/nfgtasks
2084         write (ll,'(bz,i3.3)') fg_rank
2085         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2086      &   status='unknown')
2087       endif
2088 #endif
2089       if(me.eq.king) then
2090        open(igeom,file=intname,status='unknown',access='append')
2091        open(ipdb,file=pdbname,status='unknown')
2092        open(imol2,file=mol2name,status='unknown')
2093        open(istat,file=statname,status='unknown',access='append')
2094       else
2095 c1out       open(iout,file=outname,status='unknown')
2096       endif
2097 #endif
2098       csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2099       csa_seed=prefix(:lenpre)//'.CSA.seed'
2100       csa_history=prefix(:lenpre)//'.CSA.history'
2101       csa_bank=prefix(:lenpre)//'.CSA.bank'
2102       csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2103       csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2104       csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2105 c!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2106       csa_int=prefix(:lenpre)//'.int'
2107       csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2108       csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2109       csa_in=prefix(:lenpre)//'.CSA.in'
2110 c      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2111 C Write file names
2112       if (me.eq.king)then
2113       write (iout,'(80(1h-))')
2114       write (iout,'(30x,a)') "FILE ASSIGNMENT"
2115       write (iout,'(80(1h-))')
2116       write (iout,*) "Input file                      : ",
2117      &  pref_orig(:ilen(pref_orig))//'.inp'
2118       write (iout,*) "Output file                     : ",
2119      &  outname(:ilen(outname))
2120       write (iout,*)
2121       write (iout,*) "Sidechain potential file        : ",
2122      &  sidename(:ilen(sidename))
2123 #ifndef OLDSCP
2124       write (iout,*) "SCp potential file              : ",
2125      &  scpname(:ilen(scpname))
2126 #endif
2127       write (iout,*) "Electrostatic potential file    : ",
2128      &  elename(:ilen(elename))
2129       write (iout,*) "Cumulant coefficient file       : ",
2130      &  fouriername(:ilen(fouriername))
2131       write (iout,*) "Torsional parameter file        : ",
2132      &  torname(:ilen(torname))
2133       write (iout,*) "Double torsional parameter file : ",
2134      &  tordname(:ilen(tordname))
2135       write (iout,*) "SCCOR parameter file : ",
2136      &  sccorname(:ilen(sccorname))
2137       write (iout,*) "Bond & inertia constant file    : ",
2138      &  bondname(:ilen(bondname))
2139       write (iout,*) "Bending parameter file          : ",
2140      &  thetname(:ilen(thetname))
2141       write (iout,*) "Rotamer parameter file          : ",
2142      &  rotname(:ilen(rotname))
2143       write (iout,*) "Thetpdb parameter file          : ",
2144      &  thetname_pdb(:ilen(thetname_pdb))
2145       write (iout,*) "Threading database              : ",
2146      &  patname(:ilen(patname))
2147       if (lentmp.ne.0)
2148      &write (iout,*)" DIRTMP                          : ",
2149      &  tmpdir(:lentmp)
2150       write (iout,'(80(1h-))')
2151       endif
2152       return
2153       end
2154 c----------------------------------------------------------------------------
2155       subroutine card_concat(card)
2156       implicit real*8 (a-h,o-z)
2157       include 'DIMENSIONS'
2158       include 'COMMON.IOUNITS'
2159       character*(*) card
2160       character*80 karta,ucase
2161       external ilen
2162       read (inp,'(a)') karta
2163       karta=ucase(karta)
2164       card=' '
2165       do while (karta(80:80).eq.'&')
2166         card=card(:ilen(card)+1)//karta(:79)
2167         read (inp,'(a)') karta
2168         karta=ucase(karta)
2169       enddo
2170       card=card(:ilen(card)+1)//karta
2171       return
2172       end
2173 c----------------------------------------------------------------------------------
2174       subroutine readrst
2175       implicit real*8 (a-h,o-z)
2176       include 'DIMENSIONS'
2177       include 'COMMON.CHAIN'
2178       include 'COMMON.IOUNITS'
2179       include 'COMMON.MD'
2180       open(irest2,file=rest2name,status='unknown')
2181       read(irest2,*) totT,EK,potE,totE,t_bath
2182       do i=1,2*nres
2183          read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2184       enddo
2185       do i=1,2*nres
2186          read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2187       enddo
2188       if(usampl) then
2189              read (irest2,*) iset
2190       endif
2191       close(irest2)
2192       return
2193       end
2194 c---------------------------------------------------------------------------------
2195       subroutine read_fragments
2196       implicit real*8 (a-h,o-z)
2197       include 'DIMENSIONS'
2198 #ifdef MPI
2199       include 'mpif.h'
2200 #endif
2201       include 'COMMON.SETUP'
2202       include 'COMMON.CHAIN'
2203       include 'COMMON.IOUNITS'
2204       include 'COMMON.MD'
2205       include 'COMMON.CONTROL'
2206       read(inp,*) nset,nfrag,npair,nfrag_back
2207       if(me.eq.king.or..not.out1file)
2208      & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2209      &  " nfrag_back",nfrag_back
2210       do iset=1,nset
2211          read(inp,*) mset(iset)
2212        do i=1,nfrag
2213          read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2214      &     qinfrag(i,iset)
2215          if(me.eq.king.or..not.out1file)
2216      &    write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2217      &     ifrag(2,i,iset), qinfrag(i,iset)
2218        enddo
2219        do i=1,npair
2220         read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2221      &    qinpair(i,iset)
2222         if(me.eq.king.or..not.out1file)
2223      &   write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2224      &    ipair(2,i,iset), qinpair(i,iset)
2225        enddo
2226        do i=1,nfrag_back
2227         read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2228      &     wfrag_back(3,i,iset),
2229      &     ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2230         if(me.eq.king.or..not.out1file)
2231      &   write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2232      &   wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2233        enddo
2234       enddo
2235       return
2236       end
2237 c-------------------------------------------------------------------------------
2238       subroutine read_dist_constr
2239       implicit real*8 (a-h,o-z)
2240       include 'DIMENSIONS'
2241 #ifdef MPI
2242       include 'mpif.h'
2243 #endif
2244       include 'COMMON.SETUP'
2245       include 'COMMON.CONTROL'
2246       include 'COMMON.CHAIN'
2247       include 'COMMON.IOUNITS'
2248       include 'COMMON.SBRIDGE'
2249       integer ifrag_(2,100),ipair_(2,100)
2250       double precision wfrag_(100),wpair_(100)
2251       character*500 controlcard
2252 c      write (iout,*) "Calling read_dist_constr"
2253 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2254 c      call flush(iout)
2255       call card_concat(controlcard)
2256       call readi(controlcard,"NFRAG",nfrag_,0)
2257       call readi(controlcard,"NPAIR",npair_,0)
2258       call readi(controlcard,"NDIST",ndist_,0)
2259       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2260       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2261       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2262       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2263       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2264 c      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2265 c      write (iout,*) "IFRAG"
2266 c      do i=1,nfrag_
2267 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2268 c      enddo
2269 c      write (iout,*) "IPAIR"
2270 c      do i=1,npair_
2271 c        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2272 c      enddo
2273       call flush(iout)
2274       do i=1,nfrag_
2275         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2276         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2277      &    ifrag_(2,i)=nstart_sup+nsup-1
2278 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2279         call flush(iout)
2280         if (wfrag_(i).gt.0.0d0) then
2281         do j=ifrag_(1,i),ifrag_(2,i)-1
2282           do k=j+1,ifrag_(2,i)
2283 c            write (iout,*) "j",j," k",k
2284             ddjk=dist(j,k)
2285             if (constr_dist.eq.1) then
2286             nhpb=nhpb+1
2287             ihpb(nhpb)=j
2288             jhpb(nhpb)=k
2289               dhpb(nhpb)=ddjk
2290             forcon(nhpb)=wfrag_(i)
2291             else if (constr_dist.eq.2) then
2292               if (ddjk.le.dist_cut) then
2293                 nhpb=nhpb+1
2294                 ihpb(nhpb)=j
2295                 jhpb(nhpb)=k
2296                 dhpb(nhpb)=ddjk
2297                 forcon(nhpb)=wfrag_(i)
2298               endif
2299             else
2300               nhpb=nhpb+1
2301               ihpb(nhpb)=j
2302               jhpb(nhpb)=k
2303               dhpb(nhpb)=ddjk
2304               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2305             endif
2306 #ifdef MPI
2307             if (.not.out1file .or. me.eq.king)
2308      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2309      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2310 #else
2311             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2312      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2313 #endif
2314           enddo
2315         enddo
2316         endif
2317       enddo
2318       do i=1,npair_
2319         if (wpair_(i).gt.0.0d0) then
2320         ii = ipair_(1,i)
2321         jj = ipair_(2,i)
2322         if (ii.gt.jj) then
2323           itemp=ii
2324           ii=jj
2325           jj=itemp
2326         endif
2327         do j=ifrag_(1,ii),ifrag_(2,ii)
2328           do k=ifrag_(1,jj),ifrag_(2,jj)
2329             nhpb=nhpb+1
2330             ihpb(nhpb)=j
2331             jhpb(nhpb)=k
2332             forcon(nhpb)=wpair_(i)
2333             dhpb(nhpb)=dist(j,k)
2334 #ifdef MPI
2335             if (.not.out1file .or. me.eq.king)
2336      &      write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2337      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2338 #else
2339             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2340      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2341 #endif
2342           enddo
2343         enddo
2344         endif
2345       enddo
2346       do i=1,ndist_
2347         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2348         if (forcon(nhpb+1).gt.0.0d0) then
2349           nhpb=nhpb+1
2350           dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2351 #ifdef MPI
2352           if (.not.out1file .or. me.eq.king)
2353      &    write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2354      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2355 #else
2356           write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2357      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2358 #endif
2359         endif
2360       enddo
2361       call flush(iout)
2362       return
2363       end
2364 c-------------------------------------------------------------------------------
2365 #ifdef WINIFL
2366       subroutine flush(iu)
2367       return
2368       end
2369 #endif
2370 #ifdef AIX
2371       subroutine flush(iu)
2372       call flush_(iu)
2373       return
2374       end
2375 #endif
2376 c------------------------------------------------------------------------------
2377       subroutine copy_to_tmp(source)
2378       include "DIMENSIONS"
2379       include "COMMON.IOUNITS"
2380       character*(*) source
2381       character* 256 tmpfile
2382       integer ilen
2383       external ilen
2384       logical ex
2385       tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2386       inquire(file=tmpfile,exist=ex)
2387       if (ex) then
2388         write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2389      &   " to temporary directory..."
2390         write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2391         call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2392       endif
2393       return
2394       end
2395 c------------------------------------------------------------------------------
2396       subroutine move_from_tmp(source)
2397       include "DIMENSIONS"
2398       include "COMMON.IOUNITS"
2399       character*(*) source
2400       integer ilen
2401       external ilen
2402       write (*,*) "Moving ",source(:ilen(source)),
2403      & " from temporary directory to working directory"
2404       write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2405       call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2406       return
2407       end
2408 c------------------------------------------------------------------------------
2409       subroutine random_init(seed)
2410 C
2411 C Initialize random number generator
2412 C
2413       implicit real*8 (a-h,o-z)
2414       include 'DIMENSIONS'
2415 #ifdef MPI
2416       include 'mpif.h'
2417       logical OKRandom, prng_restart
2418       real*8  r1
2419       integer iseed_array(4)
2420 #endif
2421       include 'COMMON.IOUNITS'
2422       include 'COMMON.TIME1'
2423       include 'COMMON.THREAD'
2424       include 'COMMON.SBRIDGE'
2425       include 'COMMON.CONTROL'
2426       include 'COMMON.MCM'
2427       include 'COMMON.MAP'
2428       include 'COMMON.HEADER'
2429       include 'COMMON.CSA'
2430       include 'COMMON.CHAIN'
2431       include 'COMMON.MUCA'
2432       include 'COMMON.MD'
2433       include 'COMMON.FFIELD'
2434       include 'COMMON.SETUP'
2435       iseed=-dint(dabs(seed))
2436       if (iseed.eq.0) then
2437         write (iout,'(/80(1h*)/20x,a/80(1h*))')
2438      &    'Random seed undefined. The program will stop.'
2439         write (*,'(/80(1h*)/20x,a/80(1h*))')
2440      &    'Random seed undefined. The program will stop.'
2441 #ifdef MPI
2442         call mpi_finalize(mpi_comm_world,ierr)
2443 #endif
2444         stop 'Bad random seed.'
2445       endif
2446 #ifdef MPI
2447       if (fg_rank.eq.0) then
2448       seed=seed*(me+1)+1
2449 #ifdef AMD64
2450       if(me.eq.king)
2451      &  write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
2452       OKRandom = prng_restart(me,iseed)
2453 #else
2454       do i=1,4
2455        tmp=65536.0d0**(4-i)
2456        iseed_array(i) = dint(seed/tmp)
2457        seed=seed-iseed_array(i)*tmp
2458       enddo
2459       if(me.eq.king)
2460      & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2461      &                 (iseed_array(i),i=1,4)
2462       write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2463      &                 (iseed_array(i),i=1,4)
2464       OKRandom = prng_restart(me,iseed_array)
2465 #endif
2466       if (OKRandom) then
2467 c        r1 = prng_next(me)
2468         r1=ran_number(0.0D0,1.0D0)
2469         if(me.eq.king)
2470      &   write (iout,*) 'ran_num',r1
2471         if (r1.lt.0.0d0) OKRandom=.false.
2472       endif
2473       if (.not.OKRandom) then
2474         write (iout,*) 'PRNG IS NOT WORKING!!!'
2475         print *,'PRNG IS NOT WORKING!!!'
2476         if (me.eq.0) then
2477          call flush(iout)
2478          call mpi_abort(mpi_comm_world,error_msg,ierr)
2479          stop
2480         else
2481          write (iout,*) 'too many processors for parallel prng'
2482          write (*,*) 'too many processors for parallel prng'
2483          call flush(iout)
2484          stop
2485         endif
2486       endif
2487       endif
2488 #else
2489       call vrndst(iseed)
2490       write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
2491 #endif
2492       return
2493       end