source/unres/src_MD/readrtns.F

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