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