source/unres/src_MD/readrtns.F

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