source/unres/src_MD/src/readrtns_CSA.F

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