source/unres/src_MD/readrtns.F

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