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 C Assign initial virtual bond lengths
1001         do i=2,nres
1002           vbld(i)=vbl
1003           vbld_inv(i)=vblinv
1004         enddo
1005         do i=2,nres-1
1006           vbld(i+nres)=dsc(itype(i))
1007           vbld_inv(i+nres)=dsc_inv(itype(i))
1008 c          write (iout,*) "i",i," itype",itype(i),
1009 c     &      " dsc",dsc(itype(i))," vbld",vbld(i),vbld(i+nres)
1010         enddo
1011       endif
1012 c      print *,nres
1013 c      print '(20i4)',(itype(i),i=1,nres)
1014       do i=1,nres
1015 #ifdef PROCOR
1016         if (itype(i).eq.21 .or. itype(i+1).eq.21) then
1017 #else
1018         if (itype(i).eq.21) then
1019 #endif
1020           itel(i)=0
1021 #ifdef PROCOR
1022         else if (itype(i+1).ne.20) then
1023 #else
1024         else if (itype(i).ne.20) then
1025 #endif
1026           itel(i)=1
1027         else
1028           itel(i)=2
1029         endif
1030       enddo
1031       if(me.eq.king.or..not.out1file)then
1032        write (iout,*) "ITEL"
1033        do i=1,nres-1
1034          write (iout,*) i,itype(i),itel(i)
1035        enddo
1036        print *,'Call Read_Bridge.'
1037       endif
1038       call read_bridge
1039 C 8/13/98 Set limits to generating the dihedral angles
1040       do i=1,nres
1041         phibound(1,i)=-pi
1042         phibound(2,i)=pi
1043       enddo
1044       read (inp,*) ndih_constr
1045       if (ndih_constr.gt.0) then
1046         read (inp,*) ftors
1047         read (inp,*) (idih_constr(i),phi0(i),drange(i),i=1,ndih_constr)
1048         if(me.eq.king.or..not.out1file)then
1049          write (iout,*)
1050      &   'There are',ndih_constr,' constraints on phi angles.'
1051          do i=1,ndih_constr
1052           write (iout,'(i5,2f8.3)') idih_constr(i),phi0(i),drange(i)
1053          enddo
1054         endif
1055         do i=1,ndih_constr
1056           phi0(i)=deg2rad*phi0(i)
1057           drange(i)=deg2rad*drange(i)
1058         enddo
1059         if(me.eq.king.or..not.out1file)
1060      &   write (iout,*) 'FTORS',ftors
1061         do i=1,ndih_constr
1062           ii = idih_constr(i)
1063           phibound(1,ii) = phi0(i)-drange(i)
1064           phibound(2,ii) = phi0(i)+drange(i)
1065         enddo
1066       endif
1067       nnt=1
1068 #ifdef MPI
1069       if (me.eq.king) then
1070 #endif
1071        write (iout,'(a)') 'Boundaries in phi angle sampling:'
1072        do i=1,nres
1073          write (iout,'(a3,i5,2f10.1)')
1074      &   restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
1075        enddo
1076 #ifdef MP
1077       endif
1078 #endif
1079       nct=nres
1080 cd      print *,'NNT=',NNT,' NCT=',NCT
1081       if (itype(1).eq.21) nnt=2
1082       if (itype(nres).eq.21) nct=nct-1
1083
1084 C     Bartek:READ init_vars
1085 C     Initialize variables!
1086 C     Juyong:READ read_info
1087 C     READ fragment information!!
1088 C     both routines should be in dfa.F file!!
1089
1090 #ifdef DFA
1091       if (.not. (wdfa_dist.eq.0.0 .and. wdfa_tor.eq.0.0 .and.
1092      &            wdfa_nei.eq.0.0 .and. wdfa_beta.eq.0.0)) then
1093        call init_dfa_vars
1094        print*, 'init_dfa_vars finished!'
1095        call read_dfa_info
1096        print*, 'read_dfa_info finished!'
1097       endif
1098 #endif
1099 C
1100       if (pdbref) then
1101         if(me.eq.king.or..not.out1file)
1102      &   write (iout,'(a,i3)') 'nsup=',nsup
1103         nstart_seq=nnt
1104         if (nsup.le.(nct-nnt+1)) then
1105           do i=0,nct-nnt+1-nsup
1106             if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
1107               nstart_seq=nnt+i
1108               goto 111
1109             endif
1110           enddo
1111           write (iout,'(a)')
1112      &            'Error - sequences to be superposed do not match.'
1113           stop
1114         else
1115           do i=0,nsup-(nct-nnt+1)
1116             if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
1117      &      then
1118               nstart_sup=nstart_sup+i
1119               nsup=nct-nnt+1
1120               goto 111
1121             endif
1122           enddo
1123           write (iout,'(a)')
1124      &            'Error - sequences to be superposed do not match.'
1125         endif
1126   111   continue
1127         if (nsup.eq.0) nsup=nct-nnt
1128         if (nstart_sup.eq.0) nstart_sup=nnt
1129         if (nstart_seq.eq.0) nstart_seq=nnt
1130         if(me.eq.king.or..not.out1file)
1131      &   write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
1132      &                 ' nstart_seq=',nstart_seq
1133       endif
1134 c--- Zscore rms -------
1135       if (nz_start.eq.0) nz_start=nnt
1136       if (nz_end.eq.0 .and. nsup.gt.0) then
1137         nz_end=nnt+nsup-1
1138       else if (nz_end.eq.0) then
1139         nz_end=nct
1140       endif
1141       if(me.eq.king.or..not.out1file)then
1142        write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
1143        write (iout,*) 'IZ_SC=',iz_sc
1144       endif
1145 c----------------------
1146       call init_int_table
1147       if (refstr) then
1148         if (.not.pdbref) then
1149           call read_angles(inp,*38)
1150           goto 39
1151    38     write (iout,'(a)') 'Error reading reference structure.'
1152 #ifdef MPI
1153           call MPI_Finalize(MPI_COMM_WORLD,IERROR)
1154           stop 'Error reading reference structure'
1155 #endif
1156    39     call chainbuild
1157           call setup_var
1158 czscore          call geom_to_var(nvar,coord_exp_zs(1,1))
1159           nstart_sup=nnt
1160           nstart_seq=nnt
1161           nsup=nct-nnt+1
1162           do i=1,2*nres
1163             do j=1,3
1164               cref(j,i)=c(j,i)
1165             enddo
1166           enddo
1167           call contact(.true.,ncont_ref,icont_ref,co)
1168         endif
1169         if(me.eq.king.or..not.out1file)
1170      &   write (iout,*) 'Contact order:',co
1171         if (pdbref) then
1172         if(me.eq.king.or..not.out1file)
1173      &   write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1174         do i=1,ncont_ref
1175           do j=1,2
1176             icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1177           enddo
1178           if(me.eq.king.or..not.out1file)
1179      &     write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1180      &     icont_ref(1,i),' ',
1181      &     restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1182         enddo
1183         endif
1184       endif
1185 c        write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
1186       if (constr_dist.gt.0) then
1187         call read_dist_constr
1188       endif
1189
1190
1191       if (constr_homology.gt.0) then
1192         call read_constr_homology
1193         if (indpdb.gt.0 .or. pdbref) then
1194           do i=1,2*nres
1195             do j=1,3
1196               c(j,i)=crefjlee(j,i)
1197               cref(j,i)=crefjlee(j,i)
1198             enddo
1199           enddo
1200         endif
1201 #ifdef DEBUG
1202         write (iout,*) "Array C"
1203         do i=1,nres
1204           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(c(j,i),j=1,3),
1205      &      (c(j,i+nres),j=1,3)
1206         enddo
1207         write (iout,*) "Array Cref"
1208         do i=1,nres
1209           write (iout,'(i5,3f8.3,5x,3f8.3)') i,(cref(j,i),j=1,3),
1210      &      (cref(j,i+nres),j=1,3)
1211         enddo
1212 #endif
1213        call int_from_cart1(.false.)
1214        call sc_loc_geom(.false.)
1215        do i=1,nres
1216          thetaref(i)=theta(i)
1217          phiref(i)=phi(i)
1218        enddo
1219        do i=1,nres-1
1220          do j=1,3
1221            dc(j,i)=c(j,i+1)-c(j,i)
1222            dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1223          enddo
1224        enddo
1225        do i=2,nres-1
1226          do j=1,3
1227            dc(j,i+nres)=c(j,i+nres)-c(j,i)
1228            dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1229          enddo
1230        enddo
1231       else
1232         homol_nset=0
1233       endif
1234
1235
1236       if (nhpb.gt.0) call hpb_partition
1237 c      write (iout,*) "After read_dist_constr nhpb",nhpb
1238 c      call flush(iout)
1239       if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1240      &    .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1241      &    modecalc.ne.10) then
1242 C If input structure hasn't been supplied from the PDB file read or generate
1243 C initial geometry.
1244         if (iranconf.eq.0 .and. .not. extconf) then
1245           if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1246      &     write (iout,'(a)') 'Initial geometry will be read in.'
1247           if (read_cart) then
1248             read (inp,*) time,potE,uconst,t_bath,
1249      &       nss,(ihpb(j),jhpb(j),j=1,nss), nn, (qfrag(i),i=1,nn)
1250             read(inp,'(8f10.5)',end=36,err=36)
1251      &       ((c(l,k),l=1,3),k=1,nres),
1252      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
1253             call int_from_cart1(.false.)
1254             do i=1,nres-1
1255               do j=1,3
1256                 dc(j,i)=c(j,i+1)-c(j,i)
1257                 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1258               enddo
1259             enddo
1260             do i=nnt,nct
1261               if (itype(i).ne.10) then
1262                 do j=1,3
1263                   dc(j,i+nres)=c(j,i+nres)-c(j,i)
1264                   dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1265                 enddo
1266               endif
1267             enddo
1268 c            return
1269           else
1270             call read_angles(inp,*36)
1271           endif
1272           goto 37
1273    36     write (iout,'(a)') 'Error reading angle file.'
1274 #ifdef MPI
1275           call mpi_finalize( MPI_COMM_WORLD,IERR )
1276 #endif
1277           stop 'Error reading angle file.'
1278    37     continue
1279         else if (extconf) then
1280          if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1281      &    write (iout,'(a)') 'Extended chain initial geometry.'
1282          do i=3,nres
1283           theta(i)=90d0*deg2rad
1284          enddo
1285          do i=4,nres
1286           phi(i)=180d0*deg2rad
1287          enddo
1288          do i=2,nres-1
1289           alph(i)=110d0*deg2rad
1290          enddo
1291          do i=2,nres-1
1292           omeg(i)=-120d0*deg2rad
1293          enddo
1294         else
1295           if(me.eq.king.or..not.out1file)
1296      &     write (iout,'(a)') 'Random-generated initial geometry.'
1297
1298
1299 #ifdef MPI
1300           if (me.eq.king  .or. fg_rank.eq.0 .and. (
1301      &           modecalc.eq.12 .or. modecalc.eq.14) ) then
1302 #endif
1303             do itrial=1,100
1304               itmp=1
1305               call gen_rand_conf(itmp,*30)
1306               goto 40
1307    30         write (iout,*) 'Failed to generate random conformation',
1308      &          ', itrial=',itrial
1309               write (*,*) 'Processor:',me,
1310      &          ' Failed to generate random conformation',
1311      &          ' itrial=',itrial
1312               call intout
1313
1314 #ifdef AIX
1315               call flush_(iout)
1316 #else
1317               call flush(iout)
1318 #endif
1319             enddo
1320             write (iout,'(a,i3,a)') 'Processor:',me,
1321      &        ' error in generating random conformation.'
1322             write (*,'(a,i3,a)') 'Processor:',me,
1323      &        ' error in generating random conformation.'
1324             call flush(iout)
1325 #ifdef MPI
1326             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1327    40       continue
1328           endif
1329 #else
1330    40     continue
1331 #endif
1332         endif
1333       elseif (modecalc.eq.4) then
1334         read (inp,'(a)') intinname
1335         open (intin,file=intinname,status='old',err=333)
1336         if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1337      &  write (iout,'(a)') 'intinname',intinname
1338         write (*,'(a)') 'Processor',myrank,' intinname',intinname
1339         goto 334
1340   333   write (iout,'(2a)') 'Error opening angle file ',intinname
1341 #ifdef MPI
1342         call MPI_Finalize(MPI_COMM_WORLD,IERR)
1343 #endif
1344         stop 'Error opening angle file.'
1345   334   continue
1346
1347       endif
1348 C Generate distance constraints, if the PDB structure is to be regularized.
1349       if (nthread.gt.0) then
1350         call read_threadbase
1351       endif
1352       write (iout,*) "READRTNS: Calling setup_var"
1353       call flush(iout)
1354       call setup_var
1355       if (me.eq.king .or. .not. out1file)
1356      & call intout
1357       if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1358         write (iout,'(/a,i3,a)')
1359      &  'The chain contains',ns,' disulfide-bridging cysteines.'
1360         write (iout,'(20i4)') (iss(i),i=1,ns)
1361        if (dyn_ss) then
1362           write(iout,*)"Running with dynamic disulfide-bond formation"
1363        else
1364         write (iout,'(/a/)') 'Pre-formed links are:'
1365         do i=1,nss
1366           i1=ihpb(i)-nres
1367           i2=jhpb(i)-nres
1368           it1=itype(i1)
1369           it2=itype(i2)
1370           write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1371      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1372      &    ebr,forcon(i)
1373         enddo
1374         write (iout,'(a)')
1375        endif
1376       endif
1377       if (ns.gt.0.and.dyn_ss) then
1378           do i=nss+1,nhpb
1379             ihpb(i-nss)=ihpb(i)
1380             jhpb(i-nss)=jhpb(i)
1381             forcon(i-nss)=forcon(i)
1382             dhpb(i-nss)=dhpb(i)
1383           enddo
1384           nhpb=nhpb-nss
1385           nss=0
1386           call hpb_partition
1387           do i=1,ns
1388             dyn_ss_mask(iss(i))=.true.
1389           enddo
1390       endif
1391       if (i2ndstr.gt.0) call secstrp2dihc
1392 c      call geom_to_var(nvar,x)
1393 c      call etotal(energia(0))
1394 c      call enerprint(energia(0))
1395 c      call briefout(0,etot)
1396 c      stop
1397 cd    write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1398 cd    write (iout,'(a)') 'Variable list:'
1399 cd    write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1400 #ifdef MPI
1401       if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1402      &  write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1403      &  'Processor',myrank,': end reading molecular data.'
1404 #endif
1405       return
1406       end
1407 c--------------------------------------------------------------------------
1408       logical function seq_comp(itypea,itypeb,length)
1409       implicit none
1410       integer length,itypea(length),itypeb(length)
1411       integer i
1412       do i=1,length
1413         if (itypea(i).ne.itypeb(i)) then
1414           seq_comp=.false.
1415           return
1416         endif
1417       enddo
1418       seq_comp=.true.
1419       return
1420       end
1421 c-----------------------------------------------------------------------------
1422       subroutine read_bridge
1423 C Read information about disulfide bridges.
1424       implicit real*8 (a-h,o-z)
1425       include 'DIMENSIONS'
1426 #ifdef MPI
1427       include 'mpif.h'
1428 #endif
1429       include 'COMMON.IOUNITS'
1430       include 'COMMON.GEO'
1431       include 'COMMON.VAR'
1432       include 'COMMON.INTERACT'
1433       include 'COMMON.LOCAL'
1434       include 'COMMON.NAMES'
1435       include 'COMMON.CHAIN'
1436       include 'COMMON.FFIELD'
1437       include 'COMMON.SBRIDGE'
1438       include 'COMMON.HEADER'
1439       include 'COMMON.CONTROL'
1440       include 'COMMON.DBASE'
1441       include 'COMMON.THREAD'
1442       include 'COMMON.TIME1'
1443       include 'COMMON.SETUP'
1444 C Read bridging residues.
1445       read (inp,*) ns,(iss(i),i=1,ns)
1446       print *,'ns=',ns
1447       if(me.eq.king.or..not.out1file)
1448      &  write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1449 C Check whether the specified bridging residues are cystines.
1450       do i=1,ns
1451         if (itype(iss(i)).ne.1) then
1452           if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1453      &   'Do you REALLY think that the residue ',
1454      &    restyp(itype(iss(i))),i,
1455      &   ' can form a disulfide bridge?!!!'
1456           write (*,'(2a,i3,a)')
1457      &   'Do you REALLY think that the residue ',
1458      &    restyp(itype(iss(i))),i,
1459      &   ' can form a disulfide bridge?!!!'
1460 #ifdef MPI
1461          call MPI_Finalize(MPI_COMM_WORLD,ierror)
1462          stop
1463 #endif
1464         endif
1465       enddo
1466 C Read preformed bridges.
1467       if (ns.gt.0) then
1468       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1469       if(fg_rank.eq.0)
1470      & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1471       if (nss.gt.0) then
1472         nhpb=nss
1473 C Check if the residues involved in bridges are in the specified list of
1474 C bridging residues.
1475         do i=1,nss
1476           do j=1,i-1
1477             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1478      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1479               write (iout,'(a,i3,a)') 'Disulfide pair',i,
1480      &      ' contains residues present in other pairs.'
1481               write (*,'(a,i3,a)') 'Disulfide pair',i,
1482      &      ' contains residues present in other pairs.'
1483 #ifdef MPI
1484               call MPI_Finalize(MPI_COMM_WORLD,ierror)
1485               stop
1486 #endif
1487             endif
1488           enddo
1489           do j=1,ns
1490             if (ihpb(i).eq.iss(j)) goto 10
1491           enddo
1492           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1493    10     continue
1494           do j=1,ns
1495             if (jhpb(i).eq.iss(j)) goto 20
1496           enddo
1497           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1498    20     continue
1499           dhpb(i)=dbr
1500           forcon(i)=fbr
1501         enddo
1502         do i=1,nss
1503           ihpb(i)=ihpb(i)+nres
1504           jhpb(i)=jhpb(i)+nres
1505         enddo
1506       endif
1507       endif
1508       return
1509       end
1510 c----------------------------------------------------------------------------
1511       subroutine read_x(kanal,*)
1512       implicit real*8 (a-h,o-z)
1513       include 'DIMENSIONS'
1514       include 'COMMON.GEO'
1515       include 'COMMON.VAR'
1516       include 'COMMON.CHAIN'
1517       include 'COMMON.IOUNITS'
1518       include 'COMMON.CONTROL'
1519       include 'COMMON.LOCAL'
1520       include 'COMMON.INTERACT'
1521 c Read coordinates from input
1522 c
1523       read(kanal,'(8f10.5)',end=10,err=10)
1524      &  ((c(l,k),l=1,3),k=1,nres),
1525      &  ((c(l,k+nres),l=1,3),k=nnt,nct)
1526       do j=1,3
1527         c(j,nres+1)=c(j,1)
1528         c(j,2*nres)=c(j,nres)
1529       enddo
1530       call int_from_cart1(.false.)
1531       do i=1,nres-1
1532         do j=1,3
1533           dc(j,i)=c(j,i+1)-c(j,i)
1534           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1535         enddo
1536       enddo
1537       do i=nnt,nct
1538         if (itype(i).ne.10) then
1539           do j=1,3
1540             dc(j,i+nres)=c(j,i+nres)-c(j,i)
1541             dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1542           enddo
1543         endif
1544       enddo
1545
1546       return
1547    10 return1
1548       end
1549 c----------------------------------------------------------------------------
1550       subroutine read_threadbase
1551       implicit real*8 (a-h,o-z)
1552       include 'DIMENSIONS'
1553       include 'COMMON.IOUNITS'
1554       include 'COMMON.GEO'
1555       include 'COMMON.VAR'
1556       include 'COMMON.INTERACT'
1557       include 'COMMON.LOCAL'
1558       include 'COMMON.NAMES'
1559       include 'COMMON.CHAIN'
1560       include 'COMMON.FFIELD'
1561       include 'COMMON.SBRIDGE'
1562       include 'COMMON.HEADER'
1563       include 'COMMON.CONTROL'
1564       include 'COMMON.DBASE'
1565       include 'COMMON.THREAD'
1566       include 'COMMON.TIME1'
1567 C Read pattern database for threading.
1568       read (icbase,*) nseq
1569       do i=1,nseq
1570         read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1571      &   nres_base(2,i),nres_base(3,i)
1572         read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1573      &   nres_base(1,i))
1574 c       write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1575 c    &   nres_base(2,i),nres_base(3,i)
1576 c       write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1577 c    &   nres_base(1,i))
1578       enddo
1579       close (icbase)
1580       if (weidis.eq.0.0D0) weidis=0.1D0
1581       do i=nnt,nct
1582         do j=i+2,nct
1583           nhpb=nhpb+1
1584           ihpb(nhpb)=i
1585           jhpb(nhpb)=j
1586           forcon(nhpb)=weidis
1587         enddo
1588       enddo
1589       read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1590       write (iout,'(a,i5)') 'nexcl: ',nexcl
1591       write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1592       return
1593       end
1594 c------------------------------------------------------------------------------
1595       subroutine setup_var
1596       implicit real*8 (a-h,o-z)
1597       include 'DIMENSIONS'
1598       include 'COMMON.IOUNITS'
1599       include 'COMMON.GEO'
1600       include 'COMMON.VAR'
1601       include 'COMMON.INTERACT'
1602       include 'COMMON.LOCAL'
1603       include 'COMMON.NAMES'
1604       include 'COMMON.CHAIN'
1605       include 'COMMON.FFIELD'
1606       include 'COMMON.SBRIDGE'
1607       include 'COMMON.HEADER'
1608       include 'COMMON.CONTROL'
1609       include 'COMMON.DBASE'
1610       include 'COMMON.THREAD'
1611       include 'COMMON.TIME1'
1612 C Set up variable list.
1613       ntheta=nres-2
1614       nphi=nres-3
1615       nvar=ntheta+nphi
1616       nside=0
1617       do i=2,nres-1
1618         if (itype(i).ne.10) then
1619           nside=nside+1
1620           ialph(i,1)=nvar+nside
1621           ialph(nside,2)=i
1622         endif
1623       enddo
1624       if (indphi.gt.0) then
1625         nvar=nphi
1626       else if (indback.gt.0) then
1627         nvar=nphi+ntheta
1628       else
1629         nvar=nvar+2*nside
1630       endif
1631 cd    write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1632       return
1633       end
1634 c----------------------------------------------------------------------------
1635       subroutine gen_dist_constr
1636 C Generate CA distance constraints.
1637       implicit real*8 (a-h,o-z)
1638       include 'DIMENSIONS'
1639       include 'COMMON.IOUNITS'
1640       include 'COMMON.GEO'
1641       include 'COMMON.VAR'
1642       include 'COMMON.INTERACT'
1643       include 'COMMON.LOCAL'
1644       include 'COMMON.NAMES'
1645       include 'COMMON.CHAIN'
1646       include 'COMMON.FFIELD'
1647       include 'COMMON.SBRIDGE'
1648       include 'COMMON.HEADER'
1649       include 'COMMON.CONTROL'
1650       include 'COMMON.DBASE'
1651       include 'COMMON.THREAD'
1652       include 'COMMON.TIME1'
1653       dimension itype_pdb(maxres)
1654       common /pizda/ itype_pdb
1655       character*2 iden
1656 cd      print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1657 cd      write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1658 cd     & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1659 cd     & ' nsup',nsup
1660       do i=nstart_sup,nstart_sup+nsup-1
1661 cd      write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1662 cd     &    ' seq_pdb', restyp(itype_pdb(i))
1663         do j=i+2,nstart_sup+nsup-1
1664           nhpb=nhpb+1
1665           ihpb(nhpb)=i+nstart_seq-nstart_sup
1666           jhpb(nhpb)=j+nstart_seq-nstart_sup
1667           forcon(nhpb)=weidis
1668           dhpb(nhpb)=dist(i,j)
1669         enddo
1670       enddo
1671 cd      write (iout,'(a)') 'Distance constraints:'
1672 cd      do i=nss+1,nhpb
1673 cd        ii=ihpb(i)
1674 cd        jj=jhpb(i)
1675 cd        iden='CA'
1676 cd        if (ii.gt.nres) then
1677 cd          iden='SC'
1678 cd          ii=ii-nres
1679 cd          jj=jj-nres
1680 cd        endif
1681 cd        write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1682 cd     &  restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1683 cd     &  dhpb(i),forcon(i)
1684 cd      enddo
1685       return
1686       end
1687 c----------------------------------------------------------------------------
1688       subroutine map_read
1689       implicit real*8 (a-h,o-z)
1690       include 'DIMENSIONS'
1691       include 'COMMON.MAP'
1692       include 'COMMON.IOUNITS'
1693       character*3 angid(4) /'THE','PHI','ALP','OME'/
1694       character*80 mapcard,ucase
1695       do imap=1,nmap
1696         read (inp,'(a)') mapcard
1697         mapcard=ucase(mapcard)
1698         if (index(mapcard,'PHI').gt.0) then
1699           kang(imap)=1
1700         else if (index(mapcard,'THE').gt.0) then
1701           kang(imap)=2
1702         else if (index(mapcard,'ALP').gt.0) then
1703           kang(imap)=3
1704         else if (index(mapcard,'OME').gt.0) then
1705           kang(imap)=4
1706         else
1707           write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1708           stop 'Error - illegal variable spec in MAP card.'
1709         endif
1710         call readi (mapcard,'RES1',res1(imap),0)
1711         call readi (mapcard,'RES2',res2(imap),0)
1712         if (res1(imap).eq.0) then
1713           res1(imap)=res2(imap)
1714         else if (res2(imap).eq.0) then
1715           res2(imap)=res1(imap)
1716         endif
1717         if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1718           write (iout,'(a)')
1719      &    'Error - illegal definition of variable group in MAP.'
1720           stop 'Error - illegal definition of variable group in MAP.'
1721         endif
1722         call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1723         call reada(mapcard,'TO',ang_to(imap),0.0D0)
1724         call readi(mapcard,'NSTEP',nstep(imap),0)
1725         if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1726           write (iout,'(a)')
1727      &     'Illegal boundary and/or step size specification in MAP.'
1728           stop 'Illegal boundary and/or step size specification in MAP.'
1729         endif
1730       enddo ! imap
1731       return
1732       end
1733 c----------------------------------------------------------------------------
1734 csa      subroutine csaread
1735 csa      implicit real*8 (a-h,o-z)
1736 csa      include 'DIMENSIONS'
1737 csa      include 'COMMON.IOUNITS'
1738 csa      include 'COMMON.GEO'
1739 csa      include 'COMMON.CSA'
1740 csa      include 'COMMON.BANK'
1741 csa      include 'COMMON.CONTROL'
1742 csa      character*80 ucase
1743 csa      character*620 mcmcard
1744 csa      call card_concat(mcmcard)
1745 csa
1746 csa      call readi(mcmcard,'NCONF',nconf,50)
1747 csa      call readi(mcmcard,'NADD',nadd,0)
1748 csa      call readi(mcmcard,'JSTART',jstart,1)
1749 csa      call readi(mcmcard,'JEND',jend,1)
1750 csa      call readi(mcmcard,'NSTMAX',nstmax,500000)
1751 csa      call readi(mcmcard,'N0',n0,1)
1752 csa      call readi(mcmcard,'N1',n1,6)
1753 csa      call readi(mcmcard,'N2',n2,4)
1754 csa      call readi(mcmcard,'N3',n3,0)
1755 csa      call readi(mcmcard,'N4',n4,0)
1756 csa      call readi(mcmcard,'N5',n5,0)
1757 csa      call readi(mcmcard,'N6',n6,10)
1758 csa      call readi(mcmcard,'N7',n7,0)
1759 csa      call readi(mcmcard,'N8',n8,0)
1760 csa      call readi(mcmcard,'N9',n9,0)
1761 csa      call readi(mcmcard,'N14',n14,0)
1762 csa      call readi(mcmcard,'N15',n15,0)
1763 csa      call readi(mcmcard,'N16',n16,0)
1764 csa      call readi(mcmcard,'N17',n17,0)
1765 csa      call readi(mcmcard,'N18',n18,0)
1766 csa
1767 csa      vdisulf=(index(mcmcard,'DYNSS').gt.0)
1768 csa
1769 csa      call readi(mcmcard,'NDIFF',ndiff,2)
1770 csa      call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1771 csa      call readi(mcmcard,'IS1',is1,1)
1772 csa      call readi(mcmcard,'IS2',is2,8)
1773 csa      call readi(mcmcard,'NRAN0',nran0,4)
1774 csa      call readi(mcmcard,'NRAN1',nran1,2)
1775 csa      call readi(mcmcard,'IRR',irr,1)
1776 csa      call readi(mcmcard,'NSEED',nseed,20)
1777 csa      call readi(mcmcard,'NTOTAL',ntotal,10000)
1778 csa      call reada(mcmcard,'CUT1',cut1,2.0d0)
1779 csa      call reada(mcmcard,'CUT2',cut2,5.0d0)
1780 csa      call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1781 csa      call readi(mcmcard,'ICMAX',icmax,3)
1782 csa      call readi(mcmcard,'IRESTART',irestart,0)
1783 csac!bankt      call readi(mcmcard,'NBANKTM',ntbankm,0)
1784 csa      ntbankm=0
1785 csac!bankt
1786 csa      call reada(mcmcard,'DELE',dele,20.0d0)
1787 csa      call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1788 csa      call readi(mcmcard,'IREF',iref,0)
1789 csa      call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1790 csa      call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1791 csa      call readi(mcmcard,'NCONF_IN',nconf_in,0)
1792 csa      call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1793 csa      write (iout,*) "NCONF_IN",nconf_in
1794 csa      return
1795 csa      end
1796 c----------------------------------------------------------------------------
1797 cfmc      subroutine mcmfread
1798 cfmc      implicit real*8 (a-h,o-z)
1799 cfmc      include 'DIMENSIONS'
1800 cfmc      include 'COMMON.MCMF'
1801 cfmc      include 'COMMON.IOUNITS'
1802 cfmc      include 'COMMON.GEO'
1803 cfmc      character*80 ucase
1804 cfmc      character*620 mcmcard
1805 cfmc      call card_concat(mcmcard)
1806 cfmc
1807 cfmc      call readi(mcmcard,'MAXRANT',maxrant,1000)
1808 cfmc      write(iout,*)'MAXRANT=',maxrant
1809 cfmc      call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1810 cfmc      write(iout,*)'MAXFAM=',maxfam
1811 cfmc      call readi(mcmcard,'NNET1',nnet1,5)
1812 cfmc      write(iout,*)'NNET1=',nnet1
1813 cfmc      call readi(mcmcard,'NNET2',nnet2,4)
1814 cfmc      write(iout,*)'NNET2=',nnet2
1815 cfmc      call readi(mcmcard,'NNET3',nnet3,4)
1816 cfmc      write(iout,*)'NNET3=',nnet3
1817 cfmc      call readi(mcmcard,'ILASTT',ilastt,0)
1818 cfmc      write(iout,*)'ILASTT=',ilastt
1819 cfmc      call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1820 cfmc      write(iout,*)'MAXSTR=',maxstr
1821 cfmc      maxstr_f=maxstr/maxfam
1822 cfmc      write(iout,*)'MAXSTR_F=',maxstr_f
1823 cfmc      call readi(mcmcard,'NMCMF',nmcmf,10)
1824 cfmc      write(iout,*)'NMCMF=',nmcmf
1825 cfmc      call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1826 cfmc      write(iout,*)'IFOCUS=',ifocus
1827 cfmc      call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1828 cfmc      write(iout,*)'NLOCMCMF=',nlocmcmf
1829 cfmc      call readi(mcmcard,'INTPRT',intprt,1000)
1830 cfmc      write(iout,*)'INTPRT=',intprt
1831 cfmc      call readi(mcmcard,'IPRT',iprt,100)
1832 cfmc      write(iout,*)'IPRT=',iprt
1833 cfmc      call readi(mcmcard,'IMAXTR',imaxtr,100)
1834 cfmc      write(iout,*)'IMAXTR=',imaxtr
1835 cfmc      call readi(mcmcard,'MAXEVEN',maxeven,1000)
1836 cfmc      write(iout,*)'MAXEVEN=',maxeven
1837 cfmc      call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1838 cfmc      write(iout,*)'MAXEVEN1=',maxeven1
1839 cfmc      call readi(mcmcard,'INIMIN',inimin,200)
1840 cfmc      write(iout,*)'INIMIN=',inimin
1841 cfmc      call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1842 cfmc      write(iout,*)'NSTEPMCMF=',nstepmcmf
1843 cfmc      call readi(mcmcard,'NTHREAD',nthread,5)
1844 cfmc      write(iout,*)'NTHREAD=',nthread
1845 cfmc      call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1846 cfmc      write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1847 cfmc      call readi(mcmcard,'MAXPERT',maxpert,9)
1848 cfmc      write(iout,*)'MAXPERT=',maxpert
1849 cfmc      call readi(mcmcard,'IRMSD',irmsd,1)
1850 cfmc      write(iout,*)'IRMSD=',irmsd
1851 cfmc      call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1852 cfmc      write(iout,*)'DENEMIN=',denemin
1853 cfmc      call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1854 cfmc      write(iout,*)'RCUT1S=',rcut1s
1855 cfmc      call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1856 cfmc      write(iout,*)'RCUT1E=',rcut1e
1857 cfmc      call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1858 cfmc      write(iout,*)'RCUT2S=',rcut2s
1859 cfmc      call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1860 cfmc      write(iout,*)'RCUT2E=',rcut2e
1861 cfmc      call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1862 cfmc      write(iout,*)'DPERT1=',d_pert1
1863 cfmc      call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1864 cfmc      write(iout,*)'DPERT1A=',d_pert1a
1865 cfmc      call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1866 cfmc      write(iout,*)'DPERT2=',d_pert2
1867 cfmc      call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1868 cfmc      write(iout,*)'DPERT2A=',d_pert2a
1869 cfmc      call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1870 cfmc      write(iout,*)'DPERT2B=',d_pert2b
1871 cfmc      call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1872 cfmc      write(iout,*)'DPERT2C=',d_pert2c
1873 cfmc      d_pert1=deg2rad*d_pert1
1874 cfmc      d_pert1a=deg2rad*d_pert1a
1875 cfmc      d_pert2=deg2rad*d_pert2
1876 cfmc      d_pert2a=deg2rad*d_pert2a
1877 cfmc      d_pert2b=deg2rad*d_pert2b
1878 cfmc      d_pert2c=deg2rad*d_pert2c
1879 cfmc      call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1880 cfmc      write(iout,*)'KT_MCMF1=',kt_mcmf1
1881 cfmc      call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1882 cfmc      write(iout,*)'KT_MCMF2=',kt_mcmf2
1883 cfmc      call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1884 cfmc      write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1885 cfmc      call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1886 cfmc      write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1887 cfmc      call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1888 cfmc      write(iout,*)'RCUTINI=',rcutini
1889 cfmc      call reada(mcmcard,'GRAT',grat,0.5D0)
1890 cfmc      write(iout,*)'GRAT=',grat
1891 cfmc      call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1892 cfmc      write(iout,*)'BIAS_MCMF=',bias_mcmf
1893 cfmc
1894 cfmc      return
1895 cfmc      end
1896 c----------------------------------------------------------------------------
1897       subroutine mcmread
1898       implicit real*8 (a-h,o-z)
1899       include 'DIMENSIONS'
1900       include 'COMMON.MCM'
1901       include 'COMMON.MCE'
1902       include 'COMMON.IOUNITS'
1903       character*80 ucase
1904       character*320 mcmcard
1905       call card_concat(mcmcard)
1906       call readi(mcmcard,'MAXACC',maxacc,100)
1907       call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1908       call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1909       call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1910       call readi(mcmcard,'MAXREPM',maxrepm,200)
1911       call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1912       call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1913       call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1914       call reada(mcmcard,'E_UP',e_up,5.0D0)
1915       call reada(mcmcard,'DELTE',delte,0.1D0)
1916       call readi(mcmcard,'NSWEEP',nsweep,5)
1917       call readi(mcmcard,'NSTEPH',nsteph,0)
1918       call readi(mcmcard,'NSTEPC',nstepc,0)
1919       call reada(mcmcard,'TMIN',tmin,298.0D0)
1920       call reada(mcmcard,'TMAX',tmax,298.0D0)
1921       call readi(mcmcard,'NWINDOW',nwindow,0)
1922       call readi(mcmcard,'PRINT_MC',print_mc,0)
1923       print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1924       print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1925       ent_read=(index(mcmcard,'ENT_READ').gt.0)
1926       call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1927       call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1928       call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1929       call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1930       call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1931       if (nwindow.gt.0) then
1932         read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1933         do i=1,nwindow
1934           winlen(i)=winend(i)-winstart(i)+1
1935         enddo
1936       endif
1937       if (tmax.lt.tmin) tmax=tmin
1938       if (tmax.eq.tmin) then
1939         nstepc=0
1940         nsteph=0
1941       endif
1942       if (nstepc.gt.0 .and. nsteph.gt.0) then
1943         tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1944         tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1945       endif
1946 C Probabilities of different move types
1947       sumpro_type(0)=0.0D0
1948       call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1949       call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1950       sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1951       call reada(mcmcard,'THETA'     ,sumpro_type(3),0.0d0)
1952       sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1953       call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1954       sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1955       do i=1,MaxMoveType
1956         print *,'i',i,' sumprotype',sumpro_type(i)
1957         sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1958         print *,'i',i,' sumprotype',sumpro_type(i)
1959       enddo
1960       return
1961       end
1962 c----------------------------------------------------------------------------
1963       subroutine read_minim
1964       implicit real*8 (a-h,o-z)
1965       include 'DIMENSIONS'
1966       include 'COMMON.MINIM'
1967       include 'COMMON.IOUNITS'
1968       character*80 ucase
1969       character*320 minimcard
1970       call card_concat(minimcard)
1971       call readi(minimcard,'MAXMIN',maxmin,2000)
1972       call readi(minimcard,'MAXFUN',maxfun,5000)
1973       call readi(minimcard,'MINMIN',minmin,maxmin)
1974       call readi(minimcard,'MINFUN',minfun,maxmin)
1975       call reada(minimcard,'TOLF',tolf,1.0D-2)
1976       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1977       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1978       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1979       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1980       write (iout,'(/80(1h*)/20x,a/80(1h*))')
1981      &         'Options in energy minimization:'
1982       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1983      & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1984      & 'MinMin:',MinMin,' MinFun:',MinFun,
1985      & ' TolF:',TolF,' RTolF:',RTolF
1986       return
1987       end
1988 c----------------------------------------------------------------------------
1989       subroutine read_angles(kanal,*)
1990       implicit real*8 (a-h,o-z)
1991       include 'DIMENSIONS'
1992       include 'COMMON.GEO'
1993       include 'COMMON.VAR'
1994       include 'COMMON.CHAIN'
1995       include 'COMMON.IOUNITS'
1996       include 'COMMON.CONTROL'
1997 c Read angles from input
1998 c
1999        read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
2000        read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
2001        read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
2002        read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
2003
2004        do i=1,nres
2005 c 9/7/01 avoid 180 deg valence angle
2006         if (theta(i).gt.179.99d0) theta(i)=179.99d0
2007 c
2008         theta(i)=deg2rad*theta(i)
2009         phi(i)=deg2rad*phi(i)
2010         alph(i)=deg2rad*alph(i)
2011         omeg(i)=deg2rad*omeg(i)
2012        enddo
2013       return
2014    10 return1
2015       end
2016 c----------------------------------------------------------------------------
2017       subroutine reada(rekord,lancuch,wartosc,default)
2018       implicit none
2019       character*(*) rekord,lancuch
2020       double precision wartosc,default
2021       integer ilen,iread
2022       external ilen
2023       iread=index(rekord,lancuch)
2024       if (iread.eq.0) then
2025         wartosc=default
2026         return
2027       endif
2028       iread=iread+ilen(lancuch)+1
2029       read (rekord(iread:),*,err=10,end=10) wartosc
2030       return
2031   10  wartosc=default
2032       return
2033       end
2034 c----------------------------------------------------------------------------
2035       subroutine readi(rekord,lancuch,wartosc,default)
2036       implicit none
2037       character*(*) rekord,lancuch
2038       integer wartosc,default
2039       integer ilen,iread
2040       external ilen
2041       iread=index(rekord,lancuch)
2042       if (iread.eq.0) then
2043         wartosc=default
2044         return
2045       endif
2046       iread=iread+ilen(lancuch)+1
2047       read (rekord(iread:),*,err=10,end=10) wartosc
2048       return
2049   10  wartosc=default
2050       return
2051       end
2052 c----------------------------------------------------------------------------
2053       subroutine multreadi(rekord,lancuch,tablica,dim,default)
2054       implicit none
2055       integer dim,i
2056       integer tablica(dim),default
2057       character*(*) rekord,lancuch
2058       character*80 aux
2059       integer ilen,iread
2060       external ilen
2061       do i=1,dim
2062         tablica(i)=default
2063       enddo
2064       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2065       if (iread.eq.0) return
2066       iread=iread+ilen(lancuch)+1
2067       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2068    10 return
2069       end
2070 c----------------------------------------------------------------------------
2071       subroutine multreada(rekord,lancuch,tablica,dim,default)
2072       implicit none
2073       integer dim,i
2074       double precision tablica(dim),default
2075       character*(*) rekord,lancuch
2076       character*80 aux
2077       integer ilen,iread
2078       external ilen
2079       do i=1,dim
2080         tablica(i)=default
2081       enddo
2082       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2083       if (iread.eq.0) return
2084       iread=iread+ilen(lancuch)+1
2085       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2086    10 return
2087       end
2088 c----------------------------------------------------------------------------
2089       subroutine openunits
2090       implicit real*8 (a-h,o-z)
2091       include 'DIMENSIONS'
2092 #ifdef MPI
2093       include 'mpif.h'
2094       character*16 form,nodename
2095       integer nodelen
2096 #endif
2097       include 'COMMON.SETUP'
2098       include 'COMMON.IOUNITS'
2099       include 'COMMON.MD'
2100       include 'COMMON.CONTROL'
2101       integer lenpre,lenpot,ilen,lentmp
2102       external ilen
2103       character*3 out1file_text,ucase
2104       character*3 ll
2105       external ucase
2106 c      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
2107       call getenv_loc("PREFIX",prefix)
2108       pref_orig = prefix
2109       call getenv_loc("POT",pot)
2110       call getenv_loc("DIRTMP",tmpdir)
2111       call getenv_loc("CURDIR",curdir)
2112       call getenv_loc("OUT1FILE",out1file_text)
2113 c      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
2114       out1file_text=ucase(out1file_text)
2115       if (out1file_text(1:1).eq."Y") then
2116         out1file=.true.
2117       else
2118         out1file=fg_rank.gt.0
2119       endif
2120       lenpre=ilen(prefix)
2121       lenpot=ilen(pot)
2122       lentmp=ilen(tmpdir)
2123       if (lentmp.gt.0) then
2124           write (*,'(80(1h!))')
2125           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
2126           write (*,'(80(1h!))')
2127           write (*,*)"All output files will be on node /tmp directory."
2128 #ifdef MPI
2129         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2130         if (me.eq.king) then
2131           write (*,*) "The master node is ",nodename
2132         else if (fg_rank.eq.0) then
2133           write (*,*) "I am the CG slave node ",nodename
2134         else
2135           write (*,*) "I am the FG slave node ",nodename
2136         endif
2137 #endif
2138         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2139         lenpre = lentmp+lenpre+1
2140       endif
2141       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2142 C Get the names and open the input files
2143 #if defined(WINIFL) || defined(WINPGI)
2144       open(1,file=pref_orig(:ilen(pref_orig))//
2145      &  '.inp',status='old',readonly,shared)
2146        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2147 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2148 C Get parameter filenames and open the parameter files.
2149       call getenv_loc('BONDPAR',bondname)
2150       open (ibond,file=bondname,status='old',readonly,shared)
2151       call getenv_loc('THETPAR',thetname)
2152       open (ithep,file=thetname,status='old',readonly,shared)
2153 #ifndef CRYST_THETA
2154       call getenv_loc('THETPARPDB',thetname_pdb)
2155       open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
2156 #endif
2157       call getenv_loc('ROTPAR',rotname)
2158       open (irotam,file=rotname,status='old',readonly,shared)
2159 #ifndef CRYST_SC
2160       call getenv_loc('ROTPARPDB',rotname_pdb)
2161       open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
2162 #endif
2163       call getenv_loc('TORPAR',torname)
2164       open (itorp,file=torname,status='old',readonly,shared)
2165       call getenv_loc('TORDPAR',tordname)
2166       open (itordp,file=tordname,status='old',readonly,shared)
2167       call getenv_loc('FOURIER',fouriername)
2168       open (ifourier,file=fouriername,status='old',readonly,shared)
2169       call getenv_loc('ELEPAR',elename)
2170       open (ielep,file=elename,status='old',readonly,shared)
2171       call getenv_loc('SIDEPAR',sidename)
2172       open (isidep,file=sidename,status='old',readonly,shared)
2173 #elif (defined CRAY) || (defined AIX)
2174       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2175      &  action='read')
2176 c      print *,"Processor",myrank," opened file 1"
2177       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2178 c      print *,"Processor",myrank," opened file 9"
2179 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2180 C Get parameter filenames and open the parameter files.
2181       call getenv_loc('BONDPAR',bondname)
2182       open (ibond,file=bondname,status='old',action='read')
2183 c      print *,"Processor",myrank," opened file IBOND"
2184       call getenv_loc('THETPAR',thetname)
2185       open (ithep,file=thetname,status='old',action='read')
2186 c      print *,"Processor",myrank," opened file ITHEP"
2187 #ifndef CRYST_THETA
2188       call getenv_loc('THETPARPDB',thetname_pdb)
2189       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2190 #endif
2191       call getenv_loc('ROTPAR',rotname)
2192       open (irotam,file=rotname,status='old',action='read')
2193 c      print *,"Processor",myrank," opened file IROTAM"
2194 #ifndef CRYST_SC
2195       call getenv_loc('ROTPARPDB',rotname_pdb)
2196       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2197 #endif
2198       call getenv_loc('TORPAR',torname)
2199       open (itorp,file=torname,status='old',action='read')
2200 c      print *,"Processor",myrank," opened file ITORP"
2201       call getenv_loc('TORDPAR',tordname)
2202       open (itordp,file=tordname,status='old',action='read')
2203 c      print *,"Processor",myrank," opened file ITORDP"
2204       call getenv_loc('SCCORPAR',sccorname)
2205       open (isccor,file=sccorname,status='old',action='read')
2206 c      print *,"Processor",myrank," opened file ISCCOR"
2207       call getenv_loc('FOURIER',fouriername)
2208       open (ifourier,file=fouriername,status='old',action='read')
2209 c      print *,"Processor",myrank," opened file IFOURIER"
2210       call getenv_loc('ELEPAR',elename)
2211       open (ielep,file=elename,status='old',action='read')
2212 c      print *,"Processor",myrank," opened file IELEP"
2213       call getenv_loc('SIDEPAR',sidename)
2214       open (isidep,file=sidename,status='old',action='read')
2215 c      print *,"Processor",myrank," opened file ISIDEP"
2216 c      print *,"Processor",myrank," opened parameter files"
2217 #elif (defined G77)
2218       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2219       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2220 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2221 C Get parameter filenames and open the parameter files.
2222       call getenv_loc('BONDPAR',bondname)
2223       open (ibond,file=bondname,status='old')
2224       call getenv_loc('THETPAR',thetname)
2225       open (ithep,file=thetname,status='old')
2226 #ifndef CRYST_THETA
2227       call getenv_loc('THETPARPDB',thetname_pdb)
2228       open (ithep_pdb,file=thetname_pdb,status='old')
2229 #endif
2230       call getenv_loc('ROTPAR',rotname)
2231       open (irotam,file=rotname,status='old')
2232 #ifndef CRYST_SC
2233       call getenv_loc('ROTPARPDB',rotname_pdb)
2234       open (irotam_pdb,file=rotname_pdb,status='old')
2235 #endif
2236       call getenv_loc('TORPAR',torname)
2237       open (itorp,file=torname,status='old')
2238       call getenv_loc('TORDPAR',tordname)
2239       open (itordp,file=tordname,status='old')
2240       call getenv_loc('SCCORPAR',sccorname)
2241       open (isccor,file=sccorname,status='old')
2242       call getenv_loc('FOURIER',fouriername)
2243       open (ifourier,file=fouriername,status='old')
2244       call getenv_loc('ELEPAR',elename)
2245       open (ielep,file=elename,status='old')
2246       call getenv_loc('SIDEPAR',sidename)
2247       open (isidep,file=sidename,status='old')
2248 #else
2249       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2250      &action='read')
2251        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2252 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2253 C Get parameter filenames and open the parameter files.
2254       call getenv_loc('BONDPAR',bondname)
2255       open (ibond,file=bondname,status='old',action='read')
2256       call getenv_loc('THETPAR',thetname)
2257       open (ithep,file=thetname,status='old',action='read')
2258 #ifndef CRYST_THETA
2259       call getenv_loc('THETPARPDB',thetname_pdb)
2260       print *,"thetname_pdb ",thetname_pdb
2261       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2262       print *,ithep_pdb," opened"
2263 #endif
2264       call getenv_loc('ROTPAR',rotname)
2265       open (irotam,file=rotname,status='old',action='read')
2266 #ifndef CRYST_SC
2267       call getenv_loc('ROTPARPDB',rotname_pdb)
2268       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2269 #endif
2270       call getenv_loc('TORPAR',torname)
2271       open (itorp,file=torname,status='old',action='read')
2272       call getenv_loc('TORDPAR',tordname)
2273       open (itordp,file=tordname,status='old',action='read')
2274       call getenv_loc('SCCORPAR',sccorname)
2275       open (isccor,file=sccorname,status='old',action='read')
2276       call getenv_loc('FOURIER',fouriername)
2277       open (ifourier,file=fouriername,status='old',action='read')
2278       call getenv_loc('ELEPAR',elename)
2279       open (ielep,file=elename,status='old',action='read')
2280       call getenv_loc('SIDEPAR',sidename)
2281       open (isidep,file=sidename,status='old',action='read')
2282 #endif
2283 #ifndef OLDSCP
2284 C
2285 C 8/9/01 In the newest version SCp interaction constants are read from a file
2286 C Use -DOLDSCP to use hard-coded constants instead.
2287 C
2288       call getenv_loc('SCPPAR',scpname)
2289 #if defined(WINIFL) || defined(WINPGI)
2290       open (iscpp,file=scpname,status='old',readonly,shared)
2291 #elif (defined CRAY)  || (defined AIX)
2292       open (iscpp,file=scpname,status='old',action='read')
2293 #elif (defined G77)
2294       open (iscpp,file=scpname,status='old')
2295 #else
2296       open (iscpp,file=scpname,status='old',action='read')
2297 #endif
2298 #endif
2299       call getenv_loc('PATTERN',patname)
2300 #if defined(WINIFL) || defined(WINPGI)
2301       open (icbase,file=patname,status='old',readonly,shared)
2302 #elif (defined CRAY)  || (defined AIX)
2303       open (icbase,file=patname,status='old',action='read')
2304 #elif (defined G77)
2305       open (icbase,file=patname,status='old')
2306 #else
2307       open (icbase,file=patname,status='old',action='read')
2308 #endif
2309 #ifdef MPI
2310 C Open output file only for CG processes
2311 c      print *,"Processor",myrank," fg_rank",fg_rank
2312       if (fg_rank.eq.0) then
2313
2314       if (nodes.eq.1) then
2315         npos=3
2316       else
2317         npos = dlog10(dfloat(nodes-1))+1
2318       endif
2319       if (npos.lt.3) npos=3
2320       write (liczba,'(i1)') npos
2321       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2322      &  //')'
2323       write (liczba,form) me
2324       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2325      &  liczba(:ilen(liczba))
2326       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2327      &  //'.int'
2328       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2329      &  //'.pdb'
2330       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2331      &  liczba(:ilen(liczba))//'.mol2'
2332       cartname=prefix(:lenpre)//'_'//pot(:lenpot)//
2333      &  liczba(:ilen(liczba))//'.x'
2334       statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2335      &  liczba(:ilen(liczba))//'.stat'
2336       if (lentmp.gt.0)
2337      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2338      &      //liczba(:ilen(liczba))//'.stat')
2339       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2340      &  //'.rst'
2341       if(usampl) then
2342           qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2343      & liczba(:ilen(liczba))//'.const'
2344       endif
2345
2346       endif
2347 #else
2348       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2349       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2350       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2351       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2352       cartname=prefix(:lenpre)//'_'//pot(:lenpot)//'.x'
2353       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2354       if (lentmp.gt.0)
2355      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2356      &    //'.stat')
2357       rest2name=prefix(:ilen(prefix))//'.rst'
2358       if(usampl) then
2359          qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2360       endif
2361 #endif
2362 #if defined(AIX) || defined(PGI)
2363       if (me.eq.king .or. .not. out1file)
2364      &   open(iout,file=outname,status='unknown')
2365 c#define DEBUG
2366 #ifdef DEBUG
2367       if (fg_rank.gt.0) then
2368         write (liczba,'(i3.3)') myrank/nfgtasks
2369         write (ll,'(bz,i3.3)') fg_rank
2370         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2371      &   status='unknown')
2372       endif
2373 #endif
2374 c#undef DEBUG
2375       if(me.eq.king) then
2376        open(igeom,file=intname,status='unknown',position='append')
2377        open(ipdb,file=pdbname,status='unknown')
2378        open(imol2,file=mol2name,status='unknown')
2379        open(istat,file=statname,status='unknown',position='append')
2380       else
2381 c1out       open(iout,file=outname,status='unknown')
2382       endif
2383 #else
2384       if (me.eq.king .or. .not.out1file)
2385      &    open(iout,file=outname,status='unknown')
2386 c#define DEBUG
2387 #ifdef DEBUG
2388       if (fg_rank.gt.0) then
2389         print "Processor",fg_rank," opening output file"
2390         write (liczba,'(i3.3)') myrank/nfgtasks
2391         write (ll,'(bz,i3.3)') fg_rank
2392         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2393      &   status='unknown')
2394       endif
2395 #endif
2396 c#undef DEBUG
2397       if(me.eq.king) then
2398        open(igeom,file=intname,status='unknown',access='append')
2399        open(ipdb,file=pdbname,status='unknown')
2400        open(imol2,file=mol2name,status='unknown')
2401        open(istat,file=statname,status='unknown',access='append')
2402       else
2403 c1out       open(iout,file=outname,status='unknown')
2404       endif
2405 #endif
2406 csa      csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2407 csa      csa_seed=prefix(:lenpre)//'.CSA.seed'
2408 csa      csa_history=prefix(:lenpre)//'.CSA.history'
2409 csa      csa_bank=prefix(:lenpre)//'.CSA.bank'
2410 csa      csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2411 csa      csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2412 csa      csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2413 csac!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2414 csa      csa_int=prefix(:lenpre)//'.int'
2415 csa      csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2416 csa      csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2417 csa      csa_in=prefix(:lenpre)//'.CSA.in'
2418 c      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2419 C Write file names
2420       if (me.eq.king)then
2421       write (iout,'(80(1h-))')
2422       write (iout,'(30x,a)') "FILE ASSIGNMENT"
2423       write (iout,'(80(1h-))')
2424       write (iout,*) "Input file                      : ",
2425      &  pref_orig(:ilen(pref_orig))//'.inp'
2426       write (iout,*) "Output file                     : ",
2427      &  outname(:ilen(outname))
2428       write (iout,*)
2429       write (iout,*) "Sidechain potential file        : ",
2430      &  sidename(:ilen(sidename))
2431 #ifndef OLDSCP
2432       write (iout,*) "SCp potential file              : ",
2433      &  scpname(:ilen(scpname))
2434 #endif
2435       write (iout,*) "Electrostatic potential file    : ",
2436      &  elename(:ilen(elename))
2437       write (iout,*) "Cumulant coefficient file       : ",
2438      &  fouriername(:ilen(fouriername))
2439       write (iout,*) "Torsional parameter file        : ",
2440      &  torname(:ilen(torname))
2441       write (iout,*) "Double torsional parameter file : ",
2442      &  tordname(:ilen(tordname))
2443       write (iout,*) "SCCOR parameter file : ",
2444      &  sccorname(:ilen(sccorname))
2445       write (iout,*) "Bond & inertia constant file    : ",
2446      &  bondname(:ilen(bondname))
2447       write (iout,*) "Bending parameter file          : ",
2448      &  thetname(:ilen(thetname))
2449       write (iout,*) "Rotamer parameter file          : ",
2450      &  rotname(:ilen(rotname))
2451       write (iout,*) "Threading database              : ",
2452      &  patname(:ilen(patname))
2453       if (lentmp.ne.0)
2454      &write (iout,*)" DIRTMP                          : ",
2455      &  tmpdir(:lentmp)
2456       write (iout,'(80(1h-))')
2457       endif
2458       return
2459       end
2460 c----------------------------------------------------------------------------
2461       subroutine card_concat(card)
2462       implicit real*8 (a-h,o-z)
2463       include 'DIMENSIONS'
2464       include 'COMMON.IOUNITS'
2465       character*(*) card
2466       character*80 karta,ucase
2467       external ilen
2468       read (inp,'(a)') karta
2469       karta=ucase(karta)
2470       card=' '
2471       do while (karta(80:80).eq.'&')
2472         card=card(:ilen(card)+1)//karta(:79)
2473         read (inp,'(a)') karta
2474         karta=ucase(karta)
2475       enddo
2476       card=card(:ilen(card)+1)//karta
2477       return
2478       end
2479 c----------------------------------------------------------------------------------
2480       subroutine readrst
2481       implicit real*8 (a-h,o-z)
2482       include 'DIMENSIONS'
2483       include 'COMMON.CHAIN'
2484       include 'COMMON.IOUNITS'
2485       include 'COMMON.MD'
2486       include 'COMMON.CONTROL'
2487       open(irest2,file=rest2name,status='unknown')
2488       read(irest2,*) totT,EK,potE,totE,t_bath
2489       do i=1,2*nres
2490          read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2491       enddo
2492       do i=1,2*nres
2493          read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2494       enddo
2495       if(usampl.or.homol_nset.gt.1) then
2496              read (irest2,*) iset
2497       endif
2498       close(irest2)
2499       return
2500       end
2501 c---------------------------------------------------------------------------------
2502       subroutine read_fragments
2503       implicit real*8 (a-h,o-z)
2504       include 'DIMENSIONS'
2505 #ifdef MPI
2506       include 'mpif.h'
2507 #endif
2508       include 'COMMON.SETUP'
2509       include 'COMMON.CHAIN'
2510       include 'COMMON.IOUNITS'
2511       include 'COMMON.MD'
2512       include 'COMMON.CONTROL'
2513       integer iset1
2514       read(inp,*) nset,nfrag,npair,nfrag_back
2515       if(me.eq.king.or..not.out1file)
2516      & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2517      &  " nfrag_back",nfrag_back
2518       do iset1=1,nset
2519          read(inp,*) mset(iset1)
2520        do i=1,nfrag
2521          read(inp,*) wfrag(i,iset1),ifrag(1,i,iset1),ifrag(2,i,iset1),
2522      &     qinfrag(i,iset1)
2523          if(me.eq.king.or..not.out1file)
2524      &    write(iout,*) "R ",i,wfrag(i,iset1),ifrag(1,i,iset1),
2525      &     ifrag(2,i,iset1), qinfrag(i,iset1)
2526        enddo
2527        do i=1,npair
2528         read(inp,*) wpair(i,iset1),ipair(1,i,iset1),ipair(2,i,iset1),
2529      &    qinpair(i,iset1)
2530         if(me.eq.king.or..not.out1file)
2531      &   write(iout,*) "R ",i,wpair(i,iset1),ipair(1,i,iset1),
2532      &    ipair(2,i,iset1), qinpair(i,iset1)
2533        enddo
2534        do i=1,nfrag_back
2535         read(inp,*) wfrag_back(1,i,iset1),wfrag_back(2,i,iset1),
2536      &     wfrag_back(3,i,iset1),
2537      &     ifrag_back(1,i,iset1),ifrag_back(2,i,iset1)
2538         if(me.eq.king.or..not.out1file)
2539      &   write(iout,*) "A",i,wfrag_back(1,i,iset1),
2540      &   wfrag_back(2,i,iset1),
2541      &   wfrag_back(3,i,iset1),ifrag_back(1,i,iset1),
2542      &   ifrag_back(2,i,iset1)
2543        enddo
2544       enddo
2545       return
2546       end
2547 c-------------------------------------------------------------------------------
2548       subroutine read_dist_constr
2549       implicit real*8 (a-h,o-z)
2550       include 'DIMENSIONS'
2551 #ifdef MPI
2552       include 'mpif.h'
2553 #endif
2554       include 'COMMON.SETUP'
2555       include 'COMMON.CONTROL'
2556       include 'COMMON.CHAIN'
2557       include 'COMMON.IOUNITS'
2558       include 'COMMON.SBRIDGE'
2559       integer ifrag_(2,100),ipair_(2,100)
2560       double precision wfrag_(100),wpair_(100)
2561       character*500 controlcard
2562 c      write (iout,*) "Calling read_dist_constr"
2563 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2564 c      call flush(iout)
2565       call card_concat(controlcard)
2566       call readi(controlcard,"NFRAG",nfrag_,0)
2567       call readi(controlcard,"NPAIR",npair_,0)
2568       call readi(controlcard,"NDIST",ndist_,0)
2569       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2570       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2571       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2572       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2573       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2574 c      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2575 c      write (iout,*) "IFRAG"
2576 c      do i=1,nfrag_
2577 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2578 c      enddo
2579 c      write (iout,*) "IPAIR"
2580 c      do i=1,npair_
2581 c        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2582 c      enddo
2583       if (.not.refstr .and. nfrag.gt.0) then
2584         write (iout,*)
2585      &  "ERROR: no reference structure to compute distance restraints"
2586         write (iout,*)
2587      &  "Restraints must be specified explicitly (NDIST=number)"
2588         stop
2589       endif
2590       if (nfrag.lt.2 .and. npair.gt.0) then
2591         write (iout,*) "ERROR: Less than 2 fragments specified",
2592      &   " but distance restraints between pairs requested"
2593         stop
2594       endif
2595       call flush(iout)
2596       do i=1,nfrag_
2597         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2598         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2599      &    ifrag_(2,i)=nstart_sup+nsup-1
2600 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2601         call flush(iout)
2602         if (wfrag_(i).gt.0.0d0) then
2603         do j=ifrag_(1,i),ifrag_(2,i)-1
2604           do k=j+1,ifrag_(2,i)
2605 c            write (iout,*) "j",j," k",k
2606             ddjk=dist(j,k)
2607             if (constr_dist.eq.1) then
2608             nhpb=nhpb+1
2609             ihpb(nhpb)=j
2610             jhpb(nhpb)=k
2611               dhpb(nhpb)=ddjk
2612             forcon(nhpb)=wfrag_(i)
2613             else if (constr_dist.eq.2) then
2614               if (ddjk.le.dist_cut) then
2615                 nhpb=nhpb+1
2616                 ihpb(nhpb)=j
2617                 jhpb(nhpb)=k
2618                 dhpb(nhpb)=ddjk
2619                 forcon(nhpb)=wfrag_(i)
2620               endif
2621             else
2622               nhpb=nhpb+1
2623               ihpb(nhpb)=j
2624               jhpb(nhpb)=k
2625               dhpb(nhpb)=ddjk
2626               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2627             endif
2628 #ifdef MPI
2629             if (.not.out1file .or. me.eq.king)
2630      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2631      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2632 #else
2633             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2634      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2635 #endif
2636           enddo
2637         enddo
2638         endif
2639       enddo
2640       do i=1,npair_
2641         if (wpair_(i).gt.0.0d0) then
2642         ii = ipair_(1,i)
2643         jj = ipair_(2,i)
2644         if (ii.gt.jj) then
2645           itemp=ii
2646           ii=jj
2647           jj=itemp
2648         endif
2649         do j=ifrag_(1,ii),ifrag_(2,ii)
2650           do k=ifrag_(1,jj),ifrag_(2,jj)
2651             nhpb=nhpb+1
2652             ihpb(nhpb)=j
2653             jhpb(nhpb)=k
2654             forcon(nhpb)=wpair_(i)
2655             dhpb(nhpb)=dist(j,k)
2656 #ifdef MPI
2657             if (.not.out1file .or. me.eq.king)
2658      &      write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2659      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2660 #else
2661             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2662      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2663 #endif
2664           enddo
2665         enddo
2666         endif
2667       enddo
2668       do i=1,ndist_
2669         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2670      &     ibecarb(i),forcon(nhpb+1)
2671         if (forcon(nhpb+1).gt.0.0d0) then
2672           nhpb=nhpb+1
2673           if (ibecarb(i).gt.0) then
2674             ihpb(i)=ihpb(i)+nres
2675             jhpb(i)=jhpb(i)+nres
2676           endif
2677           if (dhpb(nhpb).eq.0.0d0)
2678      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2679         endif
2680       enddo
2681 #ifdef MPI
2682       if (.not.out1file .or. me.eq.king) then
2683 #endif
2684       do i=1,nhpb
2685           write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
2686      &     i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
2687       enddo
2688       call flush(iout)
2689 #ifdef MPI
2690       endif
2691 #endif
2692       return
2693       end
2694 c-------------------------------------------------------------------------------
2695
2696       subroutine read_constr_homology
2697
2698       include 'DIMENSIONS'
2699 #ifdef MPI
2700       include 'mpif.h'
2701 #endif
2702       include 'COMMON.SETUP'
2703       include 'COMMON.CONTROL'
2704       include 'COMMON.CHAIN'
2705       include 'COMMON.IOUNITS'
2706       include 'COMMON.MD'
2707       include 'COMMON.GEO'
2708       include 'COMMON.INTERACT'
2709 c
2710 c For new homol impl
2711 c
2712       include 'COMMON.VAR'
2713 c
2714
2715 c     double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
2716 c    &                 dist_cut
2717 c     common /przechowalnia/ odl_temp(maxres,maxres,max_template),
2718 c    &    sigma_odl_temp(maxres,maxres,max_template)
2719       character*2 kic2
2720       character*24 model_ki_dist, model_ki_angle
2721       character*500 controlcard
2722       integer ki, i, j, k, l, ii_in_use(maxdim),i_tmp,idomain_tmp
2723       logical lprn /.true./
2724       integer ilen
2725       external ilen
2726 c
2727 c     FP - Nov. 2014 Temporary specifications for new vars
2728 c
2729       double precision rescore_tmp,x12,y12,z12,rescore2_tmp
2730       double precision, dimension (max_template,maxres) :: rescore
2731       double precision, dimension (max_template,maxres) :: rescore2
2732       character*24 pdbfile,tpl_k_rescore
2733 c -----------------------------------------------------------------
2734 c Reading multiple PDB ref structures and calculation of retraints
2735 c not using pre-computed ones stored in files model_ki_{dist,angle}
2736 c FP (Nov., 2014)
2737 c -----------------------------------------------------------------
2738 c
2739 c
2740 c Alternative: reading from input
2741       call card_concat(controlcard)
2742       call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
2743       call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
2744       call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
2745       call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
2746       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
2747
2748       call readi(controlcard,"HOMOL_NSET",homol_nset,1)
2749       read2sigma=(index(controlcard,'READ2SIGMA').gt.0)
2750       if (homol_nset.gt.1)then
2751          call card_concat(controlcard)
2752          read(controlcard,*) (waga_homology(i),i=1,homol_nset)
2753          if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
2754           write(iout,*) "iset homology_weight "
2755           do i=1,homol_nset
2756            write(iout,*) i,waga_homology(i)
2757           enddo
2758          endif
2759          iset=mod(kolor,homol_nset)+1
2760       else
2761        iset=1
2762        waga_homology(1)=1.0
2763       endif
2764
2765 cd      write (iout,*) "nnt",nnt," nct",nct
2766 cd      call flush(iout)
2767
2768
2769       lim_odl=0
2770       lim_dih=0
2771 c
2772 c  New
2773 c
2774       lim_theta=0
2775       lim_xx=0
2776 c
2777       write(iout,*) 'nnt=',nnt,'nct=',nct
2778 c
2779       do i = nnt,nct
2780         do k=1,constr_homology
2781          idomain(k,i)=0
2782         enddo
2783       enddo
2784
2785       ii=0
2786       do i = nnt,nct-2
2787         do j=i+2,nct
2788         ii=ii+1
2789         ii_in_use(ii)=0
2790         enddo
2791       enddo
2792
2793       do k=1,constr_homology
2794
2795         read(inp,'(a)') pdbfile
2796 c  Next stament causes error upon compilation (?)
2797 c       if(me.eq.king.or. .not. out1file)
2798 c         write (iout,'(2a)') 'PDB data will be read from file ',
2799 c    &   pdbfile(:ilen(pdbfile))
2800          write (iout,'(a,5x,a)') 'HOMOL: Opening PDB file',
2801      &  pdbfile(:ilen(pdbfile))
2802         open(ipdbin,file=pdbfile,status='old',err=33)
2803         goto 34
2804   33    write (iout,'(a,5x,a)') 'Error opening PDB file',
2805      &  pdbfile(:ilen(pdbfile))
2806         stop
2807   34    continue
2808 c        print *,'Begin reading pdb data'
2809 c
2810 c Files containing res sim or local scores (former containing sigmas)
2811 c
2812
2813         write(kic2,'(bz,i2.2)') k
2814
2815         tpl_k_rescore="template"//kic2//".sco"
2816
2817         unres_pdb=.false.
2818         if (read2sigma) then
2819           call readpdb_template(k)
2820         else
2821           call readpdb
2822         endif
2823 c
2824 c     Distance restraints
2825 c
2826 c          ... --> odl(k,ii)
2827 C Copy the coordinates from reference coordinates (?)
2828         do i=1,2*nres
2829           do j=1,3
2830             c(j,i)=cref(j,i)
2831 c           write (iout,*) "c(",j,i,") =",c(j,i)
2832           enddo
2833         enddo
2834 c
2835 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
2836 c
2837 c         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
2838           open (ientin,file=tpl_k_rescore,status='old')
2839           if (nnt.gt.1) rescore(k,1)=0.0d0
2840           do irec=nnt,maxdim ! loop for reading res sim
2841             if (read2sigma) then
2842              read (ientin,*,end=1401) i_tmp,rescore2_tmp,rescore_tmp,
2843      &                                idomain_tmp
2844              i_tmp=i_tmp+nnt-1
2845              idomain(k,i_tmp)=idomain_tmp
2846              rescore(k,i_tmp)=0.5d0*(rescore_tmp+0.5d0)
2847              rescore2(k,i_tmp)=0.5d0*(rescore2_tmp+0.5d0)
2848             else
2849              idomain(k,irec)=1
2850              read (ientin,*,end=1401) rescore_tmp
2851
2852 c           rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
2853              rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
2854 c           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
2855             endif
2856           enddo
2857  1401   continue
2858         close (ientin)
2859         if (waga_dist.ne.0.0d0) then
2860           ii=0
2861           do i = nnt,nct-2
2862             do j=i+2,nct
2863
2864             if (idomain(k,i).eq.idomain(k,j).and.idomain(k,i).ne.0) then
2865
2866               ii=ii+1
2867               ii_in_use(ii)=1
2868               l_homo(k,ii)=.true.
2869
2870 c             write (iout,*) "k",k
2871 c             write (iout,*) "i",i," j",j," constr_homology",
2872 c    &                       constr_homology
2873               ires_homo(ii)=i
2874               jres_homo(ii)=j
2875 c
2876 c Attempt to replace dist(i,j) by its definition in ...
2877 c
2878               x12=c(1,i)-c(1,j)
2879               y12=c(2,i)-c(2,j)
2880               z12=c(3,i)-c(3,j)
2881               distal=dsqrt(x12*x12+y12*y12+z12*z12)
2882               odl(k,ii)=distal
2883               if (read2sigma) then
2884                 sigma_odl(k,ii)=0
2885                 do ik=i,j
2886                  sigma_odl(k,ii)=sigma_odl(k,ii)+rescore2(k,ik)
2887                 enddo
2888                 sigma_odl(k,ii)=sigma_odl(k,ii)/(j-i+1)
2889                 if (odl(k,ii).gt.dist_cut) sigma_odl(k,ii) =
2890      &        sigma_odl(k,ii)*dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2891               else
2892                 if (odl(k,ii).le.dist_cut) then
2893                  sigma_odl(k,ii)=rescore(k,i)+rescore(k,j) ! other exprs possible
2894 c                sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
2895                 else
2896 #ifdef OLDSIGMA
2897                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error
2898      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
2899 #else
2900                  sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error
2901      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2-0.5d0)
2902 #endif
2903                 endif
2904               endif
2905 c   Following expr replaced by a positive exp argument
2906 c             sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
2907 c    &                      dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
2908
2909 c             sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)*
2910 c    &                      dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
2911 c
2912               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) ! rescore ~ error
2913 c             sigma_odl(k,ii)=sigma_odl(k,ii)*sigma_odl(k,ii)
2914 c
2915 c             sigma_odl(k,ii)=sigma_odl_temp(i,j,k)* ! new var read from file (?)
2916 c    &                        sigma_odl_temp(i,j,k)  ! not inverse because of use of res. similarity
2917             else
2918               ii=ii+1
2919               l_homo(k,ii)=.false.
2920             endif
2921             enddo
2922           enddo
2923         lim_odl=ii
2924         endif
2925 c
2926 c     Theta, dihedral and SC retraints
2927 c
2928         if (waga_angle.gt.0.0d0) then
2929 c         open (ientin,file=tpl_k_sigma_dih,status='old')
2930 c         do irec=1,maxres-3 ! loop for reading sigma_dih
2931 c            read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
2932 c            if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
2933 c            sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
2934 c    &                            sigma_dih(k,i+nnt-1)
2935 c         enddo
2936 c1402   continue
2937 c         close (ientin)
2938           do i = nnt+3,nct
2939             if (idomain(k,i).eq.0) then
2940                sigma_dih(k,i)=0.0
2941                cycle
2942             endif
2943             dih(k,i)=phiref(i) ! right?
2944 c           read (ientin,*) sigma_dih(k,i) ! original variant
2945 c             write (iout,*) "dih(",k,i,") =",dih(k,i)
2946 c             write(iout,*) "rescore(",k,i,") =",rescore(k,i),
2947 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
2948 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
2949 c    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
2950
2951             sigma_dih(k,i)=(rescore(k,i)+rescore(k,i-1)+
2952      &                     rescore(k,i-2)+rescore(k,i-3))/4.0
2953 c            if (read2sigma) sigma_dih(k,i)=sigma_dih(k,i)/4.0
2954 c           write (iout,*) "Raw sigmas for dihedral angle restraints"
2955 c           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
2956 c           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
2957 c                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
2958 c   Instead of res sim other local measure of b/b str reliability possible
2959             sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
2960 c           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
2961           enddo
2962           lim_dih=nct-nnt-2
2963         endif
2964
2965         if (waga_theta.gt.0.0d0) then
2966 c         open (ientin,file=tpl_k_sigma_theta,status='old')
2967 c         do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
2968 c            read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
2969 c            sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
2970 c    &                              sigma_theta(k,i+nnt-1)
2971 c         enddo
2972 c1403   continue
2973 c         close (ientin)
2974
2975           do i = nnt+2,nct ! right? without parallel.
2976 c         do i = i=1,nres ! alternative for bounds acc to readpdb?
2977 c         do i=ithet_start,ithet_end ! with FG parallel.
2978              if (idomain(k,i).eq.0) then
2979               sigma_theta(k,i)=0.0
2980               cycle
2981              endif
2982              thetatpl(k,i)=thetaref(i)
2983 c            write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
2984 c            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
2985 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
2986 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
2987 c            read (ientin,*) sigma_theta(k,i) ! 1st variant
2988              sigma_theta(k,i)=(rescore(k,i)+rescore(k,i-1)+
2989      &                        rescore(k,i-2))/3.0
2990 c             if (read2sigma) sigma_theta(k,i)=sigma_theta(k,i)/3.0
2991              sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
2992
2993 c            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
2994 c                             rescore(k,i-2) !  right expression ?
2995 c            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
2996           enddo
2997         endif
2998         lim_theta=nct-nnt-1
2999
3000         if (waga_d.gt.0.0d0) then
3001 c       open (ientin,file=tpl_k_sigma_d,status='old')
3002 c         do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
3003 c            read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
3004 c            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
3005 c    &                          sigma_d(k,i+nnt-1)
3006 c         enddo
3007 c1404   continue
3008
3009           do i = nnt,nct ! right? without parallel.
3010 c         do i=2,nres-1 ! alternative for bounds acc to readpdb?
3011 c         do i=loc_start,loc_end ! with FG parallel.
3012                if (itype(i).eq.10) cycle
3013                if (idomain(k,i).eq.0 ) then
3014                   sigma_d(k,i)=0.0
3015                   cycle
3016                endif
3017                xxtpl(k,i)=xxref(i)
3018                yytpl(k,i)=yyref(i)
3019                zztpl(k,i)=zzref(i)
3020 c              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
3021 c              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
3022 c              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
3023 c              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
3024                sigma_d(k,i)=rescore(k,i) !  right expression ?
3025                sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
3026
3027 c              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
3028 c              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
3029 c              read (ientin,*) sigma_d(k,i) ! 1st variant
3030                if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right?
3031           enddo
3032           lim_xx=nct-nnt+1
3033         endif
3034       enddo
3035 c
3036 c remove distance restraints not used in any model from the list
3037 c shift data in all arrays
3038 c
3039       if (waga_dist.ne.0.0d0) then
3040         ii=0
3041         do i=nnt,nct-2
3042          do j=i+2,nct
3043           ii=ii+1
3044           if (ii_in_use(ii).eq.0) then
3045              do ki=ii,lim_odl-1
3046               ires_homo(ki)=ires_homo(ki+1)
3047               jres_homo(ki)=jres_homo(ki+1)
3048               ii_in_use(ki)=ii_in_use(ki+1)
3049               do k=1,constr_homology
3050                odl(k,ki)=odl(k,ki+1)
3051                sigma_odl(k,ki)=sigma_odl(k,ki+1)
3052                l_homo(k,ki)=l_homo(k,ki+1)
3053               enddo
3054              enddo
3055              ii=ii-1
3056              lim_odl=lim_odl-1
3057           endif
3058          enddo
3059         enddo
3060       endif
3061       if (constr_homology.gt.0) call homology_partition
3062       if (constr_homology.gt.0) call init_int_table
3063 cd      write (iout,*) "homology_partition: lim_theta= ",lim_theta,
3064 cd     & "lim_xx=",lim_xx
3065 c     write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
3066 c     write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
3067 c
3068 c Print restraints
3069 c
3070       if (.not.lprn) return
3071 cd      write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
3072       if(me.eq.king .or. .not. out1file .and. fg_rank.eq.0) then
3073        write (iout,*) "Distance restraints from templates"
3074        do ii=1,lim_odl
3075        write(iout,'(3i5,10(2f8.2,1x,l1,4x))')
3076      &  ii,ires_homo(ii),jres_homo(ii),
3077      &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),l_homo(ki,ii),
3078      &  ki=1,constr_homology)
3079        enddo
3080        write (iout,*) "Dihedral angle restraints from templates"
3081        do i=nnt+3,lim_dih
3082         write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*dih(ki,i),
3083      &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
3084        enddo
3085        write (iout,*) "Virtual-bond angle restraints from templates"
3086        do i=nnt+2,lim_theta
3087         write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*thetatpl(ki,i),
3088      &      rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
3089        enddo
3090        write (iout,*) "SC restraints from templates"
3091        do i=nnt,lim_xx
3092         write(iout,'(i5,10(4f8.2,4x))') i,
3093      &  (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
3094      &   1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
3095        enddo
3096       endif
3097 c -----------------------------------------------------------------
3098       return
3099       end
3100 c----------------------------------------------------------------------
3101
3102 #ifdef WINIFL
3103       subroutine flush(iu)
3104       return
3105       end
3106 #endif
3107 #ifdef AIX
3108       subroutine flush(iu)
3109       call flush_(iu)
3110       return
3111       end
3112 #endif
3113
3114 c------------------------------------------------------------------------------
3115       subroutine copy_to_tmp(source)
3116       include "DIMENSIONS"
3117       include "COMMON.IOUNITS"
3118       character*(*) source
3119       character* 256 tmpfile
3120       integer ilen
3121       external ilen
3122       logical ex
3123       tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
3124       inquire(file=tmpfile,exist=ex)
3125       if (ex) then
3126         write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
3127      &   " to temporary directory..."
3128         write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
3129         call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
3130       endif
3131       return
3132       end
3133 c------------------------------------------------------------------------------
3134       subroutine move_from_tmp(source)
3135       include "DIMENSIONS"
3136       include "COMMON.IOUNITS"
3137       character*(*) source
3138       integer ilen
3139       external ilen
3140       write (*,*) "Moving ",source(:ilen(source)),
3141      & " from temporary directory to working directory"
3142       write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
3143       call system("/bin/mv "//source(:ilen(source))//" "//curdir)
3144       return
3145       end
3146 c------------------------------------------------------------------------------
3147       subroutine random_init(seed)
3148 C
3149 C Initialize random number generator
3150 C
3151       implicit real*8 (a-h,o-z)
3152       include 'DIMENSIONS'
3153 #ifdef AMD64
3154       integer*8 iseedi8
3155 #endif
3156 #ifdef MPI
3157       include 'mpif.h'
3158       logical OKRandom, prng_restart
3159       real*8  r1
3160       integer iseed_array(4)
3161 #endif
3162       include 'COMMON.IOUNITS'
3163       include 'COMMON.TIME1'
3164       include 'COMMON.THREAD'
3165       include 'COMMON.SBRIDGE'
3166       include 'COMMON.CONTROL'
3167       include 'COMMON.MCM'
3168       include 'COMMON.MAP'
3169       include 'COMMON.HEADER'
3170 csa      include 'COMMON.CSA'
3171       include 'COMMON.CHAIN'
3172       include 'COMMON.MUCA'
3173       include 'COMMON.MD'
3174       include 'COMMON.FFIELD'
3175       include 'COMMON.SETUP'
3176       iseed=-dint(dabs(seed))
3177       if (iseed.eq.0) then
3178         write (iout,'(/80(1h*)/20x,a/80(1h*))')
3179      &    'Random seed undefined. The program will stop.'
3180         write (*,'(/80(1h*)/20x,a/80(1h*))')
3181      &    'Random seed undefined. The program will stop.'
3182 #ifdef MPI
3183         call mpi_finalize(mpi_comm_world,ierr)
3184 #endif
3185         stop 'Bad random seed.'
3186       endif
3187 #ifdef MPI
3188       if (fg_rank.eq.0) then
3189       seed=seed*(me+1)+1
3190 #ifdef AMD64
3191       iseedi8=dint(seed)
3192       if(me.eq.king .or. .not. out1file)
3193      &  write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
3194       write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
3195       OKRandom = prng_restart(me,iseedi8)
3196 #else
3197       do i=1,4
3198        tmp=65536.0d0**(4-i)
3199        iseed_array(i) = dint(seed/tmp)
3200        seed=seed-iseed_array(i)*tmp
3201       enddo
3202       if(me.eq.king .or. .not. out1file)
3203      & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
3204      &                 (iseed_array(i),i=1,4)
3205       write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
3206      &                 (iseed_array(i),i=1,4)
3207       OKRandom = prng_restart(me,iseed_array)
3208 #endif
3209       if (OKRandom) then
3210         r1=ran_number(0.0D0,1.0D0)
3211         if(me.eq.king .or. .not. out1file)
3212      &   write (iout,*) 'ran_num',r1
3213         if (r1.lt.0.0d0) OKRandom=.false.
3214       endif
3215       if (.not.OKRandom) then
3216         write (iout,*) 'PRNG IS NOT WORKING!!!'
3217         print *,'PRNG IS NOT WORKING!!!'
3218         if (me.eq.0) then
3219          call flush(iout)
3220          call mpi_abort(mpi_comm_world,error_msg,ierr)
3221          stop
3222         else
3223          write (iout,*) 'too many processors for parallel prng'
3224          write (*,*) 'too many processors for parallel prng'
3225          call flush(iout)
3226          stop
3227         endif
3228       endif
3229       endif
3230 #else
3231       call vrndst(iseed)
3232       write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
3233 #endif
3234       return
3235       end