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