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