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