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