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