source/unres/src_MD/readrtns.F

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