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