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