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