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(iabs(itype(i)))
 934           vbld_inv(i+nres)=dsc_inv(iabs(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           if (itype(i).le.0) omeg(i)=-omeg(i)
1155          enddo
1156         else
1157           if(me.eq.king.or..not.out1file)
1158      &     write (iout,'(a)') 'Random-generated initial geometry.'
1159
1160
1161 #ifdef MPI
1162           if (me.eq.king  .or. fg_rank.eq.0 .and. (
1163      &           modecalc.eq.12 .or. modecalc.eq.14) ) then
1164 #endif
1165             do itrial=1,100
1166               itmp=1
1167               call gen_rand_conf(itmp,*30)
1168               goto 40
1169    30         write (iout,*) 'Failed to generate random conformation',
1170      &          ', itrial=',itrial
1171               write (*,*) 'Processor:',me,
1172      &          ' Failed to generate random conformation',
1173      &          ' itrial=',itrial
1174               call intout
1175
1176 #ifdef AIX
1177               call flush_(iout)
1178 #else
1179               call flush(iout)
1180 #endif
1181             enddo
1182             write (iout,'(a,i3,a)') 'Processor:',me,
1183      &        ' error in generating random conformation.'
1184             write (*,'(a,i3,a)') 'Processor:',me,
1185      &        ' error in generating random conformation.'
1186             call flush(iout)
1187 #ifdef MPI
1188             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1189    40       continue
1190           endif
1191 #else
1192    40     continue
1193 #endif
1194         endif
1195       elseif (modecalc.eq.4) then
1196         read (inp,'(a)') intinname
1197         open (intin,file=intinname,status='old',err=333)
1198         if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1199      &  write (iout,'(a)') 'intinname',intinname
1200         write (*,'(a)') 'Processor',myrank,' intinname',intinname
1201         goto 334
1202   333   write (iout,'(2a)') 'Error opening angle file ',intinname
1203 #ifdef MPI
1204         call MPI_Finalize(MPI_COMM_WORLD,IERR)
1205 #endif
1206         stop 'Error opening angle file.'
1207   334   continue
1208
1209       endif
1210 C Generate distance constraints, if the PDB structure is to be regularized.
1211       if (nthread.gt.0) then
1212         call read_threadbase
1213       endif
1214       call setup_var
1215       if (me.eq.king .or. .not. out1file)
1216      & call intout
1217       if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1218         write (iout,'(/a,i3,a)')
1219      &  'The chain contains',ns,' disulfide-bridging cysteines.'
1220         write (iout,'(20i4)') (iss(i),i=1,ns)
1221         write (iout,'(/a/)') 'Pre-formed links are:'
1222         do i=1,nss
1223           i1=ihpb(i)-nres
1224           i2=jhpb(i)-nres
1225           it1=itype(i1)
1226           it2=itype(i2)
1227           if (me.eq.king.or..not.out1file)
1228      &    write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1229      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1230      &    ebr,forcon(i)
1231         enddo
1232         write (iout,'(a)')
1233       endif
1234       if (i2ndstr.gt.0) call secstrp2dihc
1235 c      call geom_to_var(nvar,x)
1236 c      call etotal(energia(0))
1237 c      call enerprint(energia(0))
1238 c      call briefout(0,etot)
1239 c      stop
1240 cd    write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1241 cd    write (iout,'(a)') 'Variable list:'
1242 cd    write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1243 #ifdef MPI
1244       if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1245      &  write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1246      &  'Processor',myrank,': end reading molecular data.'
1247 #endif
1248       return
1249       end
1250 c--------------------------------------------------------------------------
1251       logical function seq_comp(itypea,itypeb,length)
1252       implicit none
1253       integer length,itypea(length),itypeb(length)
1254       integer i
1255       do i=1,length
1256         if (itypea(i).ne.itypeb(i)) then
1257           seq_comp=.false.
1258           return
1259         endif
1260       enddo
1261       seq_comp=.true.
1262       return
1263       end
1264 c-----------------------------------------------------------------------------
1265       subroutine read_bridge
1266 C Read information about disulfide bridges.
1267       implicit real*8 (a-h,o-z)
1268       include 'DIMENSIONS'
1269 #ifdef MPI
1270       include 'mpif.h'
1271 #endif
1272       include 'COMMON.IOUNITS'
1273       include 'COMMON.GEO'
1274       include 'COMMON.VAR'
1275       include 'COMMON.INTERACT'
1276       include 'COMMON.LOCAL'
1277       include 'COMMON.NAMES'
1278       include 'COMMON.CHAIN'
1279       include 'COMMON.FFIELD'
1280       include 'COMMON.SBRIDGE'
1281       include 'COMMON.HEADER'
1282       include 'COMMON.CONTROL'
1283       include 'COMMON.DBASE'
1284       include 'COMMON.THREAD'
1285       include 'COMMON.TIME1'
1286       include 'COMMON.SETUP'
1287 C Read bridging residues.
1288       read (inp,*) ns,(iss(i),i=1,ns)
1289       print *,'ns=',ns
1290       if(me.eq.king.or..not.out1file)
1291      &  write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1292 C Check whether the specified bridging residues are cystines.
1293       do i=1,ns
1294         if (itype(iss(i)).ne.1) then
1295           if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1296      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
1297      &   ' can form a disulfide bridge?!!!'
1298           write (*,'(2a,i3,a)')
1299      &   'Do you REALLY think that the residue ',restyp(iss(i)),i,
1300      &   ' can form a disulfide bridge?!!!'
1301 #ifdef MPI
1302          call MPI_Finalize(MPI_COMM_WORLD,ierror)
1303          stop
1304 #endif
1305         endif
1306       enddo
1307 C Read preformed bridges.
1308       if (ns.gt.0) then
1309       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1310       write (iout,*) 'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1311       if (nss.gt.0) then
1312         nhpb=nss
1313 C Check if the residues involved in bridges are in the specified list of
1314 C bridging residues.
1315         do i=1,nss
1316           do j=1,i-1
1317             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1318      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1319               write (iout,'(a,i3,a)') 'Disulfide pair',i,
1320      &      ' contains residues present in other pairs.'
1321               write (*,'(a,i3,a)') 'Disulfide pair',i,
1322      &      ' contains residues present in other pairs.'
1323 #ifdef MPI
1324               call MPI_Finalize(MPI_COMM_WORLD,ierror)
1325               stop
1326 #endif
1327             endif
1328           enddo
1329           do j=1,ns
1330             if (ihpb(i).eq.iss(j)) goto 10
1331           enddo
1332           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1333    10     continue
1334           do j=1,ns
1335             if (jhpb(i).eq.iss(j)) goto 20
1336           enddo
1337           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1338    20     continue
1339           dhpb(i)=dbr
1340           forcon(i)=fbr
1341         enddo
1342         do i=1,nss
1343           ihpb(i)=ihpb(i)+nres
1344           jhpb(i)=jhpb(i)+nres
1345         enddo
1346       endif
1347       endif
1348       return
1349       end
1350 c----------------------------------------------------------------------------
1351       subroutine read_x(kanal,*)
1352       implicit real*8 (a-h,o-z)
1353       include 'DIMENSIONS'
1354       include 'COMMON.GEO'
1355       include 'COMMON.VAR'
1356       include 'COMMON.CHAIN'
1357       include 'COMMON.IOUNITS'
1358       include 'COMMON.CONTROL'
1359       include 'COMMON.LOCAL'
1360       include 'COMMON.INTERACT'
1361 c Read coordinates from input
1362 c
1363       read(kanal,'(8f10.5)',end=10,err=10)
1364      &  ((c(l,k),l=1,3),k=1,nres),
1365      &  ((c(l,k+nres),l=1,3),k=nnt,nct)
1366       do j=1,3
1367         c(j,nres+1)=c(j,1)
1368         c(j,2*nres)=c(j,nres)
1369       enddo
1370       call int_from_cart1(.false.)
1371       do i=1,nres-1
1372         do j=1,3
1373           dc(j,i)=c(j,i+1)-c(j,i)
1374           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1375         enddo
1376       enddo
1377       do i=nnt,nct
1378         if (itype(i).ne.10) then
1379           do j=1,3
1380             dc(j,i+nres)=c(j,i+nres)-c(j,i)
1381             dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1382           enddo
1383         endif
1384       enddo
1385
1386       return
1387    10 return1
1388       end
1389 c----------------------------------------------------------------------------
1390       subroutine read_threadbase
1391       implicit real*8 (a-h,o-z)
1392       include 'DIMENSIONS'
1393       include 'COMMON.IOUNITS'
1394       include 'COMMON.GEO'
1395       include 'COMMON.VAR'
1396       include 'COMMON.INTERACT'
1397       include 'COMMON.LOCAL'
1398       include 'COMMON.NAMES'
1399       include 'COMMON.CHAIN'
1400       include 'COMMON.FFIELD'
1401       include 'COMMON.SBRIDGE'
1402       include 'COMMON.HEADER'
1403       include 'COMMON.CONTROL'
1404       include 'COMMON.DBASE'
1405       include 'COMMON.THREAD'
1406       include 'COMMON.TIME1'
1407 C Read pattern database for threading.
1408       read (icbase,*) nseq
1409       do i=1,nseq
1410         read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1411      &   nres_base(2,i),nres_base(3,i)
1412         read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1413      &   nres_base(1,i))
1414 c       write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1415 c    &   nres_base(2,i),nres_base(3,i)
1416 c       write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1417 c    &   nres_base(1,i))
1418       enddo
1419       close (icbase)
1420       if (weidis.eq.0.0D0) weidis=0.1D0
1421       do i=nnt,nct
1422         do j=i+2,nct
1423           nhpb=nhpb+1
1424           ihpb(nhpb)=i
1425           jhpb(nhpb)=j
1426           forcon(nhpb)=weidis
1427         enddo
1428       enddo
1429       read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1430       write (iout,'(a,i5)') 'nexcl: ',nexcl
1431       write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1432       return
1433       end
1434 c------------------------------------------------------------------------------
1435       subroutine setup_var
1436       implicit real*8 (a-h,o-z)
1437       include 'DIMENSIONS'
1438       include 'COMMON.IOUNITS'
1439       include 'COMMON.GEO'
1440       include 'COMMON.VAR'
1441       include 'COMMON.INTERACT'
1442       include 'COMMON.LOCAL'
1443       include 'COMMON.NAMES'
1444       include 'COMMON.CHAIN'
1445       include 'COMMON.FFIELD'
1446       include 'COMMON.SBRIDGE'
1447       include 'COMMON.HEADER'
1448       include 'COMMON.CONTROL'
1449       include 'COMMON.DBASE'
1450       include 'COMMON.THREAD'
1451       include 'COMMON.TIME1'
1452 C Set up variable list.
1453       ntheta=nres-2
1454       nphi=nres-3
1455       nvar=ntheta+nphi
1456       nside=0
1457       do i=2,nres-1
1458         if (itype(i).ne.10) then
1459           nside=nside+1
1460           ialph(i,1)=nvar+nside
1461           ialph(nside,2)=i
1462         endif
1463       enddo
1464       if (indphi.gt.0) then
1465         nvar=nphi
1466       else if (indback.gt.0) then
1467         nvar=nphi+ntheta
1468       else
1469         nvar=nvar+2*nside
1470       endif
1471 cd    write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1472       return
1473       end
1474 c----------------------------------------------------------------------------
1475       subroutine gen_dist_constr
1476 C Generate CA distance constraints.
1477       implicit real*8 (a-h,o-z)
1478       include 'DIMENSIONS'
1479       include 'COMMON.IOUNITS'
1480       include 'COMMON.GEO'
1481       include 'COMMON.VAR'
1482       include 'COMMON.INTERACT'
1483       include 'COMMON.LOCAL'
1484       include 'COMMON.NAMES'
1485       include 'COMMON.CHAIN'
1486       include 'COMMON.FFIELD'
1487       include 'COMMON.SBRIDGE'
1488       include 'COMMON.HEADER'
1489       include 'COMMON.CONTROL'
1490       include 'COMMON.DBASE'
1491       include 'COMMON.THREAD'
1492       include 'COMMON.TIME1'
1493       dimension itype_pdb(maxres)
1494       common /pizda/ itype_pdb
1495       character*2 iden
1496 cd      print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1497 cd      write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1498 cd     & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1499 cd     & ' nsup',nsup
1500       do i=nstart_sup,nstart_sup+nsup-1
1501 cd      write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1502 cd     &    ' seq_pdb', restyp(itype_pdb(i))
1503         do j=i+2,nstart_sup+nsup-1
1504           nhpb=nhpb+1
1505           ihpb(nhpb)=i+nstart_seq-nstart_sup
1506           jhpb(nhpb)=j+nstart_seq-nstart_sup
1507           forcon(nhpb)=weidis
1508           dhpb(nhpb)=dist(i,j)
1509         enddo
1510       enddo
1511 cd      write (iout,'(a)') 'Distance constraints:'
1512 cd      do i=nss+1,nhpb
1513 cd        ii=ihpb(i)
1514 cd        jj=jhpb(i)
1515 cd        iden='CA'
1516 cd        if (ii.gt.nres) then
1517 cd          iden='SC'
1518 cd          ii=ii-nres
1519 cd          jj=jj-nres
1520 cd        endif
1521 cd        write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1522 cd     &  restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1523 cd     &  dhpb(i),forcon(i)
1524 cd      enddo
1525       return
1526       end
1527 c----------------------------------------------------------------------------
1528       subroutine map_read
1529       implicit real*8 (a-h,o-z)
1530       include 'DIMENSIONS'
1531       include 'COMMON.MAP'
1532       include 'COMMON.IOUNITS'
1533       character*3 angid(4) /'THE','PHI','ALP','OME'/
1534       character*80 mapcard,ucase
1535       do imap=1,nmap
1536         read (inp,'(a)') mapcard
1537         mapcard=ucase(mapcard)
1538         if (index(mapcard,'PHI').gt.0) then
1539           kang(imap)=1
1540         else if (index(mapcard,'THE').gt.0) then
1541           kang(imap)=2
1542         else if (index(mapcard,'ALP').gt.0) then
1543           kang(imap)=3
1544         else if (index(mapcard,'OME').gt.0) then
1545           kang(imap)=4
1546         else
1547           write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1548           stop 'Error - illegal variable spec in MAP card.'
1549         endif
1550         call readi (mapcard,'RES1',res1(imap),0)
1551         call readi (mapcard,'RES2',res2(imap),0)
1552         if (res1(imap).eq.0) then
1553           res1(imap)=res2(imap)
1554         else if (res2(imap).eq.0) then
1555           res2(imap)=res1(imap)
1556         endif
1557         if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1558           write (iout,'(a)')
1559      &    'Error - illegal definition of variable group in MAP.'
1560           stop 'Error - illegal definition of variable group in MAP.'
1561         endif
1562         call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1563         call reada(mapcard,'TO',ang_to(imap),0.0D0)
1564         call readi(mapcard,'NSTEP',nstep(imap),0)
1565         if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1566           write (iout,'(a)')
1567      &     'Illegal boundary and/or step size specification in MAP.'
1568           stop 'Illegal boundary and/or step size specification in MAP.'
1569         endif
1570       enddo ! imap
1571       return
1572       end
1573 c----------------------------------------------------------------------------
1574 csa      subroutine csaread
1575 csa      implicit real*8 (a-h,o-z)
1576 csa      include 'DIMENSIONS'
1577 csa      include 'COMMON.IOUNITS'
1578 csa      include 'COMMON.GEO'
1579 csa      include 'COMMON.CSA'
1580 csa      include 'COMMON.BANK'
1581 csa      include 'COMMON.CONTROL'
1582 csa      character*80 ucase
1583 csa      character*620 mcmcard
1584 csa      call card_concat(mcmcard)
1585 csa
1586 csa      call readi(mcmcard,'NCONF',nconf,50)
1587 csa      call readi(mcmcard,'NADD',nadd,0)
1588 csa      call readi(mcmcard,'JSTART',jstart,1)
1589 csa      call readi(mcmcard,'JEND',jend,1)
1590 csa      call readi(mcmcard,'NSTMAX',nstmax,500000)
1591 csa      call readi(mcmcard,'N0',n0,1)
1592 csa      call readi(mcmcard,'N1',n1,6)
1593 csa      call readi(mcmcard,'N2',n2,4)
1594 csa      call readi(mcmcard,'N3',n3,0)
1595 csa      call readi(mcmcard,'N4',n4,0)
1596 csa      call readi(mcmcard,'N5',n5,0)
1597 csa      call readi(mcmcard,'N6',n6,10)
1598 csa      call readi(mcmcard,'N7',n7,0)
1599 csa      call readi(mcmcard,'N8',n8,0)
1600 csa      call readi(mcmcard,'N9',n9,0)
1601 csa      call readi(mcmcard,'N14',n14,0)
1602 csa      call readi(mcmcard,'N15',n15,0)
1603 csa      call readi(mcmcard,'N16',n16,0)
1604 csa      call readi(mcmcard,'N17',n17,0)
1605 csa      call readi(mcmcard,'N18',n18,0)
1606 csa
1607 csa      vdisulf=(index(mcmcard,'DYNSS').gt.0)
1608 csa
1609 csa      call readi(mcmcard,'NDIFF',ndiff,2)
1610 csa      call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1611 csa      call readi(mcmcard,'IS1',is1,1)
1612 csa      call readi(mcmcard,'IS2',is2,8)
1613 csa      call readi(mcmcard,'NRAN0',nran0,4)
1614 csa      call readi(mcmcard,'NRAN1',nran1,2)
1615 csa      call readi(mcmcard,'IRR',irr,1)
1616 csa      call readi(mcmcard,'NSEED',nseed,20)
1617 csa      call readi(mcmcard,'NTOTAL',ntotal,10000)
1618 csa      call reada(mcmcard,'CUT1',cut1,2.0d0)
1619 csa      call reada(mcmcard,'CUT2',cut2,5.0d0)
1620 csa      call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1621 csa      call readi(mcmcard,'ICMAX',icmax,3)
1622 csa      call readi(mcmcard,'IRESTART',irestart,0)
1623 csac!bankt      call readi(mcmcard,'NBANKTM',ntbankm,0)
1624 csa      ntbankm=0
1625 csac!bankt
1626 csa      call reada(mcmcard,'DELE',dele,20.0d0)
1627 csa      call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1628 csa      call readi(mcmcard,'IREF',iref,0)
1629 csa      call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1630 csa      call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1631 csa      call readi(mcmcard,'NCONF_IN',nconf_in,0)
1632 csa      call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1633 csa      write (iout,*) "NCONF_IN",nconf_in
1634 csa      return
1635 csa      end
1636 c----------------------------------------------------------------------------
1637 cfmc      subroutine mcmfread
1638 cfmc      implicit real*8 (a-h,o-z)
1639 cfmc      include 'DIMENSIONS'
1640 cfmc      include 'COMMON.MCMF'
1641 cfmc      include 'COMMON.IOUNITS'
1642 cfmc      include 'COMMON.GEO'
1643 cfmc      character*80 ucase
1644 cfmc      character*620 mcmcard
1645 cfmc      call card_concat(mcmcard)
1646 cfmc
1647 cfmc      call readi(mcmcard,'MAXRANT',maxrant,1000)
1648 cfmc      write(iout,*)'MAXRANT=',maxrant
1649 cfmc      call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1650 cfmc      write(iout,*)'MAXFAM=',maxfam
1651 cfmc      call readi(mcmcard,'NNET1',nnet1,5)
1652 cfmc      write(iout,*)'NNET1=',nnet1
1653 cfmc      call readi(mcmcard,'NNET2',nnet2,4)
1654 cfmc      write(iout,*)'NNET2=',nnet2
1655 cfmc      call readi(mcmcard,'NNET3',nnet3,4)
1656 cfmc      write(iout,*)'NNET3=',nnet3
1657 cfmc      call readi(mcmcard,'ILASTT',ilastt,0)
1658 cfmc      write(iout,*)'ILASTT=',ilastt
1659 cfmc      call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1660 cfmc      write(iout,*)'MAXSTR=',maxstr
1661 cfmc      maxstr_f=maxstr/maxfam
1662 cfmc      write(iout,*)'MAXSTR_F=',maxstr_f
1663 cfmc      call readi(mcmcard,'NMCMF',nmcmf,10)
1664 cfmc      write(iout,*)'NMCMF=',nmcmf
1665 cfmc      call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1666 cfmc      write(iout,*)'IFOCUS=',ifocus
1667 cfmc      call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1668 cfmc      write(iout,*)'NLOCMCMF=',nlocmcmf
1669 cfmc      call readi(mcmcard,'INTPRT',intprt,1000)
1670 cfmc      write(iout,*)'INTPRT=',intprt
1671 cfmc      call readi(mcmcard,'IPRT',iprt,100)
1672 cfmc      write(iout,*)'IPRT=',iprt
1673 cfmc      call readi(mcmcard,'IMAXTR',imaxtr,100)
1674 cfmc      write(iout,*)'IMAXTR=',imaxtr
1675 cfmc      call readi(mcmcard,'MAXEVEN',maxeven,1000)
1676 cfmc      write(iout,*)'MAXEVEN=',maxeven
1677 cfmc      call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1678 cfmc      write(iout,*)'MAXEVEN1=',maxeven1
1679 cfmc      call readi(mcmcard,'INIMIN',inimin,200)
1680 cfmc      write(iout,*)'INIMIN=',inimin
1681 cfmc      call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1682 cfmc      write(iout,*)'NSTEPMCMF=',nstepmcmf
1683 cfmc      call readi(mcmcard,'NTHREAD',nthread,5)
1684 cfmc      write(iout,*)'NTHREAD=',nthread
1685 cfmc      call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1686 cfmc      write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1687 cfmc      call readi(mcmcard,'MAXPERT',maxpert,9)
1688 cfmc      write(iout,*)'MAXPERT=',maxpert
1689 cfmc      call readi(mcmcard,'IRMSD',irmsd,1)
1690 cfmc      write(iout,*)'IRMSD=',irmsd
1691 cfmc      call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1692 cfmc      write(iout,*)'DENEMIN=',denemin
1693 cfmc      call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1694 cfmc      write(iout,*)'RCUT1S=',rcut1s
1695 cfmc      call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1696 cfmc      write(iout,*)'RCUT1E=',rcut1e
1697 cfmc      call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1698 cfmc      write(iout,*)'RCUT2S=',rcut2s
1699 cfmc      call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1700 cfmc      write(iout,*)'RCUT2E=',rcut2e
1701 cfmc      call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1702 cfmc      write(iout,*)'DPERT1=',d_pert1
1703 cfmc      call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1704 cfmc      write(iout,*)'DPERT1A=',d_pert1a
1705 cfmc      call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1706 cfmc      write(iout,*)'DPERT2=',d_pert2
1707 cfmc      call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1708 cfmc      write(iout,*)'DPERT2A=',d_pert2a
1709 cfmc      call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1710 cfmc      write(iout,*)'DPERT2B=',d_pert2b
1711 cfmc      call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1712 cfmc      write(iout,*)'DPERT2C=',d_pert2c
1713 cfmc      d_pert1=deg2rad*d_pert1
1714 cfmc      d_pert1a=deg2rad*d_pert1a
1715 cfmc      d_pert2=deg2rad*d_pert2
1716 cfmc      d_pert2a=deg2rad*d_pert2a
1717 cfmc      d_pert2b=deg2rad*d_pert2b
1718 cfmc      d_pert2c=deg2rad*d_pert2c
1719 cfmc      call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1720 cfmc      write(iout,*)'KT_MCMF1=',kt_mcmf1
1721 cfmc      call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1722 cfmc      write(iout,*)'KT_MCMF2=',kt_mcmf2
1723 cfmc      call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1724 cfmc      write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1725 cfmc      call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1726 cfmc      write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1727 cfmc      call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1728 cfmc      write(iout,*)'RCUTINI=',rcutini
1729 cfmc      call reada(mcmcard,'GRAT',grat,0.5D0)
1730 cfmc      write(iout,*)'GRAT=',grat
1731 cfmc      call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1732 cfmc      write(iout,*)'BIAS_MCMF=',bias_mcmf
1733 cfmc
1734 cfmc      return
1735 cfmc      end
1736 c----------------------------------------------------------------------------
1737       subroutine mcmread
1738       implicit real*8 (a-h,o-z)
1739       include 'DIMENSIONS'
1740       include 'COMMON.MCM'
1741       include 'COMMON.MCE'
1742       include 'COMMON.IOUNITS'
1743       character*80 ucase
1744       character*320 mcmcard
1745       call card_concat(mcmcard)
1746       call readi(mcmcard,'MAXACC',maxacc,100)
1747       call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1748       call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1749       call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1750       call readi(mcmcard,'MAXREPM',maxrepm,200)
1751       call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1752       call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1753       call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1754       call reada(mcmcard,'E_UP',e_up,5.0D0)
1755       call reada(mcmcard,'DELTE',delte,0.1D0)
1756       call readi(mcmcard,'NSWEEP',nsweep,5)
1757       call readi(mcmcard,'NSTEPH',nsteph,0)
1758       call readi(mcmcard,'NSTEPC',nstepc,0)
1759       call reada(mcmcard,'TMIN',tmin,298.0D0)
1760       call reada(mcmcard,'TMAX',tmax,298.0D0)
1761       call readi(mcmcard,'NWINDOW',nwindow,0)
1762       call readi(mcmcard,'PRINT_MC',print_mc,0)
1763       print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1764       print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1765       ent_read=(index(mcmcard,'ENT_READ').gt.0)
1766       call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1767       call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1768       call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1769       call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1770       call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1771       if (nwindow.gt.0) then
1772         read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1773         do i=1,nwindow
1774           winlen(i)=winend(i)-winstart(i)+1
1775         enddo
1776       endif
1777       if (tmax.lt.tmin) tmax=tmin
1778       if (tmax.eq.tmin) then
1779         nstepc=0
1780         nsteph=0
1781       endif
1782       if (nstepc.gt.0 .and. nsteph.gt.0) then
1783         tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1784         tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1785       endif
1786 C Probabilities of different move types
1787       sumpro_type(0)=0.0D0
1788       call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1789       call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1790       sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1791       call reada(mcmcard,'THETA'     ,sumpro_type(3),0.0d0)
1792       sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1793       call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1794       sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1795       do i=1,MaxMoveType
1796         print *,'i',i,' sumprotype',sumpro_type(i)
1797         sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1798         print *,'i',i,' sumprotype',sumpro_type(i)
1799       enddo
1800       return
1801       end
1802 c----------------------------------------------------------------------------
1803       subroutine read_minim
1804       implicit real*8 (a-h,o-z)
1805       include 'DIMENSIONS'
1806       include 'COMMON.MINIM'
1807       include 'COMMON.IOUNITS'
1808       character*80 ucase
1809       character*320 minimcard
1810       call card_concat(minimcard)
1811       call readi(minimcard,'MAXMIN',maxmin,2000)
1812       call readi(minimcard,'MAXFUN',maxfun,5000)
1813       call readi(minimcard,'MINMIN',minmin,maxmin)
1814       call readi(minimcard,'MINFUN',minfun,maxmin)
1815       call reada(minimcard,'TOLF',tolf,1.0D-2)
1816       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1817       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1818       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1819       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1820       write (iout,'(/80(1h*)/20x,a/80(1h*))')
1821      &         'Options in energy minimization:'
1822       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1823      & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1824      & 'MinMin:',MinMin,' MinFun:',MinFun,
1825      & ' TolF:',TolF,' RTolF:',RTolF
1826       return
1827       end
1828 c----------------------------------------------------------------------------
1829       subroutine read_angles(kanal,*)
1830       implicit real*8 (a-h,o-z)
1831       include 'DIMENSIONS'
1832       include 'COMMON.GEO'
1833       include 'COMMON.VAR'
1834       include 'COMMON.CHAIN'
1835       include 'COMMON.IOUNITS'
1836       include 'COMMON.CONTROL'
1837 c Read angles from input
1838 c
1839        read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1840        read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1841        read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1842        read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1843
1844        do i=1,nres
1845 c 9/7/01 avoid 180 deg valence angle
1846         if (theta(i).gt.179.99d0) theta(i)=179.99d0
1847 c
1848         theta(i)=deg2rad*theta(i)
1849         phi(i)=deg2rad*phi(i)
1850         alph(i)=deg2rad*alph(i)
1851         omeg(i)=deg2rad*omeg(i)
1852        enddo
1853       return
1854    10 return1
1855       end
1856 c----------------------------------------------------------------------------
1857       subroutine reada(rekord,lancuch,wartosc,default)
1858       implicit none
1859       character*(*) rekord,lancuch
1860       double precision wartosc,default
1861       integer ilen,iread
1862       external ilen
1863       iread=index(rekord,lancuch)
1864       if (iread.eq.0) then
1865         wartosc=default
1866         return
1867       endif
1868       iread=iread+ilen(lancuch)+1
1869       read (rekord(iread:),*,err=10,end=10) wartosc
1870       return
1871   10  wartosc=default
1872       return
1873       end
1874 c----------------------------------------------------------------------------
1875       subroutine readi(rekord,lancuch,wartosc,default)
1876       implicit none
1877       character*(*) rekord,lancuch
1878       integer wartosc,default
1879       integer ilen,iread
1880       external ilen
1881       iread=index(rekord,lancuch)
1882       if (iread.eq.0) then
1883         wartosc=default
1884         return
1885       endif
1886       iread=iread+ilen(lancuch)+1
1887       read (rekord(iread:),*,err=10,end=10) wartosc
1888       return
1889   10  wartosc=default
1890       return
1891       end
1892 c----------------------------------------------------------------------------
1893       subroutine multreadi(rekord,lancuch,tablica,dim,default)
1894       implicit none
1895       integer dim,i
1896       integer tablica(dim),default
1897       character*(*) rekord,lancuch
1898       character*80 aux
1899       integer ilen,iread
1900       external ilen
1901       do i=1,dim
1902         tablica(i)=default
1903       enddo
1904       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1905       if (iread.eq.0) return
1906       iread=iread+ilen(lancuch)+1
1907       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1908    10 return
1909       end
1910 c----------------------------------------------------------------------------
1911       subroutine multreada(rekord,lancuch,tablica,dim,default)
1912       implicit none
1913       integer dim,i
1914       double precision tablica(dim),default
1915       character*(*) rekord,lancuch
1916       character*80 aux
1917       integer ilen,iread
1918       external ilen
1919       do i=1,dim
1920         tablica(i)=default
1921       enddo
1922       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1923       if (iread.eq.0) return
1924       iread=iread+ilen(lancuch)+1
1925       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1926    10 return
1927       end
1928 c----------------------------------------------------------------------------
1929       subroutine openunits
1930       implicit real*8 (a-h,o-z)
1931       include 'DIMENSIONS'
1932 #ifdef MPI
1933       include 'mpif.h'
1934       character*16 form,nodename
1935       integer nodelen
1936 #endif
1937       include 'COMMON.SETUP'
1938       include 'COMMON.IOUNITS'
1939       include 'COMMON.MD'
1940       include 'COMMON.CONTROL'
1941       integer lenpre,lenpot,ilen,lentmp
1942       external ilen
1943       character*3 out1file_text,ucase
1944       character*3 ll
1945       external ucase
1946 c      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1947       call getenv_loc("PREFIX",prefix)
1948       pref_orig = prefix
1949       call getenv_loc("POT",pot)
1950       call getenv_loc("DIRTMP",tmpdir)
1951       call getenv_loc("CURDIR",curdir)
1952       call getenv_loc("OUT1FILE",out1file_text)
1953 c      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1954       out1file_text=ucase(out1file_text)
1955       if (out1file_text(1:1).eq."Y") then
1956         out1file=.true.
1957       else
1958         out1file=fg_rank.gt.0
1959       endif
1960       lenpre=ilen(prefix)
1961       lenpot=ilen(pot)
1962       lentmp=ilen(tmpdir)
1963       if (lentmp.gt.0) then
1964           write (*,'(80(1h!))')
1965           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
1966           write (*,'(80(1h!))')
1967           write (*,*)"All output files will be on node /tmp directory."
1968 #ifdef MPI
1969         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1970         if (me.eq.king) then
1971           write (*,*) "The master node is ",nodename
1972         else if (fg_rank.eq.0) then
1973           write (*,*) "I am the CG slave node ",nodename
1974         else
1975           write (*,*) "I am the FG slave node ",nodename
1976         endif
1977 #endif
1978         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1979         lenpre = lentmp+lenpre+1
1980       endif
1981       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1982 C Get the names and open the input files
1983 #if defined(WINIFL) || defined(WINPGI)
1984       open(1,file=pref_orig(:ilen(pref_orig))//
1985      &  '.inp',status='old',readonly,shared)
1986        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1987 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1988 C Get parameter filenames and open the parameter files.
1989       call getenv_loc('BONDPAR',bondname)
1990       open (ibond,file=bondname,status='old',readonly,shared)
1991       call getenv_loc('THETPAR',thetname)
1992       open (ithep,file=thetname,status='old',readonly,shared)
1993 #ifndef CRYST_THETA
1994       call getenv_loc('THETPARPDB',thetname_pdb)
1995       open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
1996 #endif
1997       call getenv_loc('ROTPAR',rotname)
1998       open (irotam,file=rotname,status='old',readonly,shared)
1999 #ifndef CRYST_SC
2000       call getenv_loc('ROTPARPDB',rotname_pdb)
2001       open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
2002 #endif
2003       call getenv_loc('TORPAR',torname)
2004       open (itorp,file=torname,status='old',readonly,shared)
2005       call getenv_loc('TORDPAR',tordname)
2006       open (itordp,file=tordname,status='old',readonly,shared)
2007       call getenv_loc('FOURIER',fouriername)
2008       open (ifourier,file=fouriername,status='old',readonly,shared)
2009       call getenv_loc('ELEPAR',elename)
2010       open (ielep,file=elename,status='old',readonly,shared)
2011       call getenv_loc('SIDEPAR',sidename)
2012       open (isidep,file=sidename,status='old',readonly,shared)
2013 #elif (defined CRAY) || (defined AIX)
2014       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2015      &  action='read')
2016 c      print *,"Processor",myrank," opened file 1"
2017       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2018 c      print *,"Processor",myrank," opened file 9"
2019 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2020 C Get parameter filenames and open the parameter files.
2021       call getenv_loc('BONDPAR',bondname)
2022       open (ibond,file=bondname,status='old',action='read')
2023 c      print *,"Processor",myrank," opened file IBOND"
2024       call getenv_loc('THETPAR',thetname)
2025       open (ithep,file=thetname,status='old',action='read')
2026 c      print *,"Processor",myrank," opened file ITHEP"
2027 #ifndef CRYST_THETA
2028       call getenv_loc('THETPARPDB',thetname_pdb)
2029       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2030 #endif
2031       call getenv_loc('ROTPAR',rotname)
2032       open (irotam,file=rotname,status='old',action='read')
2033 c      print *,"Processor",myrank," opened file IROTAM"
2034 #ifndef CRYST_SC
2035       call getenv_loc('ROTPARPDB',rotname_pdb)
2036       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2037 #endif
2038       call getenv_loc('TORPAR',torname)
2039       open (itorp,file=torname,status='old',action='read')
2040 c      print *,"Processor",myrank," opened file ITORP"
2041       call getenv_loc('TORDPAR',tordname)
2042       open (itordp,file=tordname,status='old',action='read')
2043 c      print *,"Processor",myrank," opened file ITORDP"
2044       call getenv_loc('SCCORPAR',sccorname)
2045       open (isccor,file=sccorname,status='old',action='read')
2046 c      print *,"Processor",myrank," opened file ISCCOR"
2047       call getenv_loc('FOURIER',fouriername)
2048       open (ifourier,file=fouriername,status='old',action='read')
2049 c      print *,"Processor",myrank," opened file IFOURIER"
2050       call getenv_loc('ELEPAR',elename)
2051       open (ielep,file=elename,status='old',action='read')
2052 c      print *,"Processor",myrank," opened file IELEP"
2053       call getenv_loc('SIDEPAR',sidename)
2054       open (isidep,file=sidename,status='old',action='read')
2055 c      print *,"Processor",myrank," opened file ISIDEP"
2056 c      print *,"Processor",myrank," opened parameter files"
2057 #elif (defined G77)
2058       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2059       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2060 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2061 C Get parameter filenames and open the parameter files.
2062       call getenv_loc('BONDPAR',bondname)
2063       open (ibond,file=bondname,status='old')
2064       call getenv_loc('THETPAR',thetname)
2065       open (ithep,file=thetname,status='old')
2066 #ifndef CRYST_THETA
2067       call getenv_loc('THETPARPDB',thetname_pdb)
2068       open (ithep_pdb,file=thetname_pdb,status='old')
2069 #endif
2070       call getenv_loc('ROTPAR',rotname)
2071       open (irotam,file=rotname,status='old')
2072 #ifndef CRYST_SC
2073       call getenv_loc('ROTPARPDB',rotname_pdb)
2074       open (irotam_pdb,file=rotname_pdb,status='old')
2075 #endif
2076       call getenv_loc('TORPAR',torname)
2077       open (itorp,file=torname,status='old')
2078       call getenv_loc('TORDPAR',tordname)
2079       open (itordp,file=tordname,status='old')
2080       call getenv_loc('SCCORPAR',sccorname)
2081       open (isccor,file=sccorname,status='old')
2082       call getenv_loc('FOURIER',fouriername)
2083       open (ifourier,file=fouriername,status='old')
2084       call getenv_loc('ELEPAR',elename)
2085       open (ielep,file=elename,status='old')
2086       call getenv_loc('SIDEPAR',sidename)
2087       open (isidep,file=sidename,status='old')
2088 #else
2089       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2090      &action='read')
2091        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2092 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2093 C Get parameter filenames and open the parameter files.
2094       call getenv_loc('BONDPAR',bondname)
2095       open (ibond,file=bondname,status='old',action='read')
2096       call getenv_loc('THETPAR',thetname)
2097       open (ithep,file=thetname,status='old',action='read')
2098 #ifndef CRYST_THETA
2099       call getenv_loc('THETPARPDB',thetname_pdb)
2100       print *,"thetname_pdb ",thetname_pdb
2101       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2102       print *,ithep_pdb," opened"
2103 #endif
2104       call getenv_loc('ROTPAR',rotname)
2105       open (irotam,file=rotname,status='old',action='read')
2106 #ifndef CRYST_SC
2107       call getenv_loc('ROTPARPDB',rotname_pdb)
2108       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2109 #endif
2110       call getenv_loc('TORPAR',torname)
2111       open (itorp,file=torname,status='old',action='read')
2112       call getenv_loc('TORDPAR',tordname)
2113       open (itordp,file=tordname,status='old',action='read')
2114       call getenv_loc('SCCORPAR',sccorname)
2115       open (isccor,file=sccorname,status='old',action='read')
2116       call getenv_loc('FOURIER',fouriername)
2117       open (ifourier,file=fouriername,status='old',action='read')
2118       call getenv_loc('ELEPAR',elename)
2119       open (ielep,file=elename,status='old',action='read')
2120       call getenv_loc('SIDEPAR',sidename)
2121       open (isidep,file=sidename,status='old',action='read')
2122 #endif
2123 #ifndef OLDSCP
2124 C
2125 C 8/9/01 In the newest version SCp interaction constants are read from a file
2126 C Use -DOLDSCP to use hard-coded constants instead.
2127 C
2128       call getenv_loc('SCPPAR',scpname)
2129 #if defined(WINIFL) || defined(WINPGI)
2130       open (iscpp,file=scpname,status='old',readonly,shared)
2131 #elif (defined CRAY)  || (defined AIX)
2132       open (iscpp,file=scpname,status='old',action='read')
2133 #elif (defined G77)
2134       open (iscpp,file=scpname,status='old')
2135 #else
2136       open (iscpp,file=scpname,status='old',action='read')
2137 #endif
2138 #endif
2139       call getenv_loc('PATTERN',patname)
2140 #if defined(WINIFL) || defined(WINPGI)
2141       open (icbase,file=patname,status='old',readonly,shared)
2142 #elif (defined CRAY)  || (defined AIX)
2143       open (icbase,file=patname,status='old',action='read')
2144 #elif (defined G77)
2145       open (icbase,file=patname,status='old')
2146 #else
2147       open (icbase,file=patname,status='old',action='read')
2148 #endif
2149 #ifdef MPI
2150 C Open output file only for CG processes
2151 c      print *,"Processor",myrank," fg_rank",fg_rank
2152       if (fg_rank.eq.0) then
2153
2154       if (nodes.eq.1) then
2155         npos=3
2156       else
2157         npos = dlog10(dfloat(nodes-1))+1
2158       endif
2159       if (npos.lt.3) npos=3
2160       write (liczba,'(i1)') npos
2161       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2162      &  //')'
2163       write (liczba,form) me
2164       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2165      &  liczba(:ilen(liczba))
2166       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2167      &  //'.int'
2168       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2169      &  //'.pdb'
2170       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2171      &  liczba(:ilen(liczba))//'.mol2'
2172       statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2173      &  liczba(:ilen(liczba))//'.stat'
2174       if (lentmp.gt.0)
2175      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2176      &      //liczba(:ilen(liczba))//'.stat')
2177       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2178      &  //'.rst'
2179       if(usampl) then
2180           qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2181      & liczba(:ilen(liczba))//'.const'
2182       endif
2183
2184       endif
2185 #else
2186       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2187       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2188       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2189       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2190       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2191       if (lentmp.gt.0)
2192      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2193      &    //'.stat')
2194       rest2name=prefix(:ilen(prefix))//'.rst'
2195       if(usampl) then
2196          qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2197       endif
2198 #endif
2199 #if defined(AIX) || defined(PGI)
2200       if (me.eq.king .or. .not. out1file)
2201      &   open(iout,file=outname,status='unknown')
2202 c#define DEBUG
2203 #ifdef DEBUG
2204       if (fg_rank.gt.0) then
2205         write (liczba,'(i3.3)') myrank/nfgtasks
2206         write (ll,'(bz,i3.3)') fg_rank
2207         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2208      &   status='unknown')
2209       endif
2210 #endif
2211 c#undef DEBUG
2212       if(me.eq.king) then
2213        open(igeom,file=intname,status='unknown',position='append')
2214        open(ipdb,file=pdbname,status='unknown')
2215        open(imol2,file=mol2name,status='unknown')
2216        open(istat,file=statname,status='unknown',position='append')
2217       else
2218 c1out       open(iout,file=outname,status='unknown')
2219       endif
2220 #else
2221       if (me.eq.king .or. .not.out1file)
2222      &    open(iout,file=outname,status='unknown')
2223 c#define DEBUG
2224 #ifdef DEBUG
2225       if (fg_rank.gt.0) then
2226         print "Processor",fg_rank," opening output file"
2227         write (liczba,'(i3.3)') myrank/nfgtasks
2228         write (ll,'(bz,i3.3)') fg_rank
2229         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2230      &   status='unknown')
2231       endif
2232 #endif
2233 c#undef DEBUG
2234       if(me.eq.king) then
2235        open(igeom,file=intname,status='unknown',access='append')
2236        open(ipdb,file=pdbname,status='unknown')
2237        open(imol2,file=mol2name,status='unknown')
2238        open(istat,file=statname,status='unknown',access='append')
2239       else
2240 c1out       open(iout,file=outname,status='unknown')
2241       endif
2242 #endif
2243 csa      csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2244 csa      csa_seed=prefix(:lenpre)//'.CSA.seed'
2245 csa      csa_history=prefix(:lenpre)//'.CSA.history'
2246 csa      csa_bank=prefix(:lenpre)//'.CSA.bank'
2247 csa      csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2248 csa      csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2249 csa      csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2250 csac!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2251 csa      csa_int=prefix(:lenpre)//'.int'
2252 csa      csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2253 csa      csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2254 csa      csa_in=prefix(:lenpre)//'.CSA.in'
2255 c      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2256 C Write file names
2257       if (me.eq.king)then
2258       write (iout,'(80(1h-))')
2259       write (iout,'(30x,a)') "FILE ASSIGNMENT"
2260       write (iout,'(80(1h-))')
2261       write (iout,*) "Input file                      : ",
2262      &  pref_orig(:ilen(pref_orig))//'.inp'
2263       write (iout,*) "Output file                     : ",
2264      &  outname(:ilen(outname))
2265       write (iout,*)
2266       write (iout,*) "Sidechain potential file        : ",
2267      &  sidename(:ilen(sidename))
2268 #ifndef OLDSCP
2269       write (iout,*) "SCp potential file              : ",
2270      &  scpname(:ilen(scpname))
2271 #endif
2272       write (iout,*) "Electrostatic potential file    : ",
2273      &  elename(:ilen(elename))
2274       write (iout,*) "Cumulant coefficient file       : ",
2275      &  fouriername(:ilen(fouriername))
2276       write (iout,*) "Torsional parameter file        : ",
2277      &  torname(:ilen(torname))
2278       write (iout,*) "Double torsional parameter file : ",
2279      &  tordname(:ilen(tordname))
2280       write (iout,*) "SCCOR parameter file : ",
2281      &  sccorname(:ilen(sccorname))
2282       write (iout,*) "Bond & inertia constant file    : ",
2283      &  bondname(:ilen(bondname))
2284       write (iout,*) "Bending parameter file          : ",
2285      &  thetname(:ilen(thetname))
2286       write (iout,*) "Rotamer parameter file          : ",
2287      &  rotname(:ilen(rotname))
2288       write (iout,*) "Threading database              : ",
2289      &  patname(:ilen(patname))
2290       if (lentmp.ne.0)
2291      &write (iout,*)" DIRTMP                          : ",
2292      &  tmpdir(:lentmp)
2293       write (iout,'(80(1h-))')
2294       endif
2295       return
2296       end
2297 c----------------------------------------------------------------------------
2298       subroutine card_concat(card)
2299       implicit real*8 (a-h,o-z)
2300       include 'DIMENSIONS'
2301       include 'COMMON.IOUNITS'
2302       character*(*) card
2303       character*80 karta,ucase
2304       external ilen
2305       read (inp,'(a)') karta
2306       karta=ucase(karta)
2307       card=' '
2308       do while (karta(80:80).eq.'&')
2309         card=card(:ilen(card)+1)//karta(:79)
2310         read (inp,'(a)') karta
2311         karta=ucase(karta)
2312       enddo
2313       card=card(:ilen(card)+1)//karta
2314       return
2315       end
2316 c----------------------------------------------------------------------------------
2317       subroutine readrst
2318       implicit real*8 (a-h,o-z)
2319       include 'DIMENSIONS'
2320       include 'COMMON.CHAIN'
2321       include 'COMMON.IOUNITS'
2322       include 'COMMON.MD'
2323       open(irest2,file=rest2name,status='unknown')
2324       read(irest2,*) totT,EK,potE,totE,t_bath
2325       do i=1,2*nres
2326          read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2327       enddo
2328       do i=1,2*nres
2329          read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2330       enddo
2331       if(usampl) then
2332              read (irest2,*) iset
2333       endif
2334       close(irest2)
2335       return
2336       end
2337 c---------------------------------------------------------------------------------
2338       subroutine read_fragments
2339       implicit real*8 (a-h,o-z)
2340       include 'DIMENSIONS'
2341 #ifdef MPI
2342       include 'mpif.h'
2343 #endif
2344       include 'COMMON.SETUP'
2345       include 'COMMON.CHAIN'
2346       include 'COMMON.IOUNITS'
2347       include 'COMMON.MD'
2348       include 'COMMON.CONTROL'
2349       read(inp,*) nset,nfrag,npair,nfrag_back
2350       if(me.eq.king.or..not.out1file)
2351      & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2352      &  " nfrag_back",nfrag_back
2353       do iset=1,nset
2354          read(inp,*) mset(iset)
2355        do i=1,nfrag
2356          read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2357      &     qinfrag(i,iset)
2358          if(me.eq.king.or..not.out1file)
2359      &    write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2360      &     ifrag(2,i,iset), qinfrag(i,iset)
2361        enddo
2362        do i=1,npair
2363         read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2364      &    qinpair(i,iset)
2365         if(me.eq.king.or..not.out1file)
2366      &   write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2367      &    ipair(2,i,iset), qinpair(i,iset)
2368        enddo
2369        do i=1,nfrag_back
2370         read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2371      &     wfrag_back(3,i,iset),
2372      &     ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2373         if(me.eq.king.or..not.out1file)
2374      &   write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2375      &   wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2376        enddo
2377       enddo
2378       return
2379       end
2380 c-------------------------------------------------------------------------------
2381       subroutine read_dist_constr
2382       implicit real*8 (a-h,o-z)
2383       include 'DIMENSIONS'
2384 #ifdef MPI
2385       include 'mpif.h'
2386 #endif
2387       include 'COMMON.SETUP'
2388       include 'COMMON.CONTROL'
2389       include 'COMMON.CHAIN'
2390       include 'COMMON.IOUNITS'
2391       include 'COMMON.SBRIDGE'
2392       integer ifrag_(2,100),ipair_(2,100)
2393       double precision wfrag_(100),wpair_(100)
2394       character*500 controlcard
2395 c      write (iout,*) "Calling read_dist_constr"
2396 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2397 c      call flush(iout)
2398       call card_concat(controlcard)
2399       call readi(controlcard,"NFRAG",nfrag_,0)
2400       call readi(controlcard,"NPAIR",npair_,0)
2401       call readi(controlcard,"NDIST",ndist_,0)
2402       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2403       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2404       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2405       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2406       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2407 c      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2408 c      write (iout,*) "IFRAG"
2409 c      do i=1,nfrag_
2410 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2411 c      enddo
2412 c      write (iout,*) "IPAIR"
2413 c      do i=1,npair_
2414 c        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2415 c      enddo
2416       if (.not.refstr .and. nfrag.gt.0) then
2417         write (iout,*)
2418      &  "ERROR: no reference structure to compute distance restraints"
2419         write (iout,*)
2420      &  "Restraints must be specified explicitly (NDIST=number)"
2421         stop
2422       endif
2423       if (nfrag.lt.2 .and. npair.gt.0) then
2424         write (iout,*) "ERROR: Less than 2 fragments specified",
2425      &   " but distance restraints between pairs requested"
2426         stop
2427       endif
2428       call flush(iout)
2429       do i=1,nfrag_
2430         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2431         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2432      &    ifrag_(2,i)=nstart_sup+nsup-1
2433 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2434         call flush(iout)
2435         if (wfrag_(i).gt.0.0d0) then
2436         do j=ifrag_(1,i),ifrag_(2,i)-1
2437           do k=j+1,ifrag_(2,i)
2438             write (iout,*) "j",j," k",k
2439             ddjk=dist(j,k)
2440             if (constr_dist.eq.1) then
2441             nhpb=nhpb+1
2442             ihpb(nhpb)=j
2443             jhpb(nhpb)=k
2444               dhpb(nhpb)=ddjk
2445             forcon(nhpb)=wfrag_(i)
2446             else if (constr_dist.eq.2) then
2447               if (ddjk.le.dist_cut) then
2448                 nhpb=nhpb+1
2449                 ihpb(nhpb)=j
2450                 jhpb(nhpb)=k
2451                 dhpb(nhpb)=ddjk
2452                 forcon(nhpb)=wfrag_(i)
2453               endif
2454             else
2455               nhpb=nhpb+1
2456               ihpb(nhpb)=j
2457               jhpb(nhpb)=k
2458               dhpb(nhpb)=ddjk
2459               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2460             endif
2461 #ifdef MPI
2462             if (.not.out1file .or. me.eq.king)
2463      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2464      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2465 #else
2466             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2467      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2468 #endif
2469           enddo
2470         enddo
2471         endif
2472       enddo
2473       do i=1,npair_
2474         if (wpair_(i).gt.0.0d0) then
2475         ii = ipair_(1,i)
2476         jj = ipair_(2,i)
2477         if (ii.gt.jj) then
2478           itemp=ii
2479           ii=jj
2480           jj=itemp
2481         endif
2482         do j=ifrag_(1,ii),ifrag_(2,ii)
2483           do k=ifrag_(1,jj),ifrag_(2,jj)
2484             nhpb=nhpb+1
2485             ihpb(nhpb)=j
2486             jhpb(nhpb)=k
2487             forcon(nhpb)=wpair_(i)
2488             dhpb(nhpb)=dist(j,k)
2489 #ifdef MPI
2490             if (.not.out1file .or. me.eq.king)
2491      &      write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2492      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2493 #else
2494             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2495      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2496 #endif
2497           enddo
2498         enddo
2499         endif
2500       enddo
2501       do i=1,ndist_
2502         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2503      &     ibecarb(i),forcon(nhpb+1)
2504         if (forcon(nhpb+1).gt.0.0d0) then
2505           nhpb=nhpb+1
2506           if (ibecarb(i).gt.0) then
2507             ihpb(i)=ihpb(i)+nres
2508             jhpb(i)=jhpb(i)+nres
2509           endif
2510           if (dhpb(nhpb).eq.0.0d0)
2511      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2512         endif
2513       enddo
2514 #ifdef MPI
2515       if (.not.out1file .or. me.eq.king) then
2516 #endif
2517       do i=1,nhpb
2518           write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
2519      &     i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
2520       enddo
2521       call flush(iout)
2522 #ifdef MPI
2523       endif
2524 #endif
2525       return
2526       end
2527 c-------------------------------------------------------------------------------
2528 #ifdef WINIFL
2529       subroutine flush(iu)
2530       return
2531       end
2532 #endif
2533 #ifdef AIX
2534       subroutine flush(iu)
2535       call flush_(iu)
2536       return
2537       end
2538 #endif
2539 c------------------------------------------------------------------------------
2540       subroutine copy_to_tmp(source)
2541       include "DIMENSIONS"
2542       include "COMMON.IOUNITS"
2543       character*(*) source
2544       character* 256 tmpfile
2545       integer ilen
2546       external ilen
2547       logical ex
2548       tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2549       inquire(file=tmpfile,exist=ex)
2550       if (ex) then
2551         write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2552      &   " to temporary directory..."
2553         write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2554         call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2555       endif
2556       return
2557       end
2558 c------------------------------------------------------------------------------
2559       subroutine move_from_tmp(source)
2560       include "DIMENSIONS"
2561       include "COMMON.IOUNITS"
2562       character*(*) source
2563       integer ilen
2564       external ilen
2565       write (*,*) "Moving ",source(:ilen(source)),
2566      & " from temporary directory to working directory"
2567       write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2568       call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2569       return
2570       end
2571 c------------------------------------------------------------------------------
2572       subroutine random_init(seed)
2573 C
2574 C Initialize random number generator
2575 C
2576       implicit real*8 (a-h,o-z)
2577       include 'DIMENSIONS'
2578 #ifdef AMD64
2579       integer*8 iseedi8
2580 #endif
2581 #ifdef MPI
2582       include 'mpif.h'
2583       logical OKRandom, prng_restart
2584       real*8  r1
2585       integer iseed_array(4)
2586 #endif
2587       include 'COMMON.IOUNITS'
2588       include 'COMMON.TIME1'
2589       include 'COMMON.THREAD'
2590       include 'COMMON.SBRIDGE'
2591       include 'COMMON.CONTROL'
2592       include 'COMMON.MCM'
2593       include 'COMMON.MAP'
2594       include 'COMMON.HEADER'
2595 csa      include 'COMMON.CSA'
2596       include 'COMMON.CHAIN'
2597       include 'COMMON.MUCA'
2598       include 'COMMON.MD'
2599       include 'COMMON.FFIELD'
2600       include 'COMMON.SETUP'
2601       iseed=-dint(dabs(seed))
2602       if (iseed.eq.0) then
2603         write (iout,'(/80(1h*)/20x,a/80(1h*))')
2604      &    'Random seed undefined. The program will stop.'
2605         write (*,'(/80(1h*)/20x,a/80(1h*))')
2606      &    'Random seed undefined. The program will stop.'
2607 #ifdef MPI
2608         call mpi_finalize(mpi_comm_world,ierr)
2609 #endif
2610         stop 'Bad random seed.'
2611       endif
2612 #ifdef MPI
2613       if (fg_rank.eq.0) then
2614       seed=seed*(me+1)+1
2615 #ifdef AMD64
2616       iseedi8=dint(seed)
2617       if(me.eq.king .or. .not. out1file)
2618      &  write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
2619       write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
2620       OKRandom = prng_restart(me,iseedi8)
2621 #else
2622       do i=1,4
2623        tmp=65536.0d0**(4-i)
2624        iseed_array(i) = dint(seed/tmp)
2625        seed=seed-iseed_array(i)*tmp
2626       enddo
2627       if(me.eq.king .or. .not. out1file)
2628      & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2629      &                 (iseed_array(i),i=1,4)
2630       write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2631      &                 (iseed_array(i),i=1,4)
2632       OKRandom = prng_restart(me,iseed_array)
2633 #endif
2634       if (OKRandom) then
2635         r1=ran_number(0.0D0,1.0D0)
2636         if(me.eq.king .or. .not. out1file)
2637      &   write (iout,*) 'ran_num',r1
2638         if (r1.lt.0.0d0) OKRandom=.false.
2639       endif
2640       if (.not.OKRandom) then
2641         write (iout,*) 'PRNG IS NOT WORKING!!!'
2642         print *,'PRNG IS NOT WORKING!!!'
2643         if (me.eq.0) then
2644          call flush(iout)
2645          call mpi_abort(mpi_comm_world,error_msg,ierr)
2646          stop
2647         else
2648          write (iout,*) 'too many processors for parallel prng'
2649          write (*,*) 'too many processors for parallel prng'
2650          call flush(iout)
2651          stop
2652         endif
2653       endif
2654       endif
2655 #else
2656       call vrndst(iseed)
2657       write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
2658 #endif
2659       return
2660       end