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