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