source/unres/src_MD/readrtns.F

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