source/unres/src_MD/readrtns.F

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