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