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