source/unres/src_MD/readrtns.F

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