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
 860       do i=1,nres
 861         thetabound(1,i)=0
 862         thetabound(2,i)=pi
 863       enddo
 864 C begin reading theta constrains this is quartic constrains allowing to
 865 C have smooth second derivative
 866       if (with_theta_constr) then
 867 C with_theta_constr is keyword allowing for occurance of theta constrains
 868       read (inp,*) ntheta_constr
 869 C ntheta_constr is the number of theta constrains
 870       if (ntheta_constr.gt.0) then
 871 C        read (inp,*) ftors
 872         read (inp,*) (itheta_constr(i),theta_constr0(i),
 873      &  theta_drange(i),for_thet_constr(i),
 874      &  i=1,ntheta_constr)
 875 C the above code reads from 1 to ntheta_constr
 876 C itheta_constr(i) residue i for which is theta_constr
 877 C theta_constr0 the global minimum value
 878 C theta_drange is range for which there is no energy penalty
 879 C for_thet_constr is the force constant for quartic energy penalty
 880 C E=k*x**4
 881         if(me.eq.king.or..not.out1file)then
 882          write (iout,*)
 883      &   'There are',ntheta_constr,' constraints on phi angles.'
 884          do i=1,ntheta_constr
 885           write (iout,'(i5,3f8.3)') itheta_constr(i),theta_constr0(i),
 886      &    theta_drange(i),
 887      &    for_thet_constr(i)
 888          enddo
 889         endif
 890         do i=1,nthet_constr
 891           theta_constr0(i)=deg2rad*theta_constr0(i)
 892           theta_drange(i)=deg2rad*theta_drange(i)
 893         enddo
 894 C        if(me.eq.king.or..not.out1file)
 895 C     &   write (iout,*) 'FTORS',ftors
 896         do i=1,ntheta_constr
 897           ii = itheta_constr(i)
 898           thetabound(1,ii) = phi0(i)-drange(i)
 899           thetabound(2,ii) = phi0(i)+drange(i)
 900         enddo
 901       endif ! ntheta_constr.gt.0
 902       endif! with_theta_constr
 903 C
 904 C      with_dihed_constr = index(controlcard,"WITH_DIHED_CONSTR").gt.0
 905 C      write (iout,*) "with_dihed_constr ",with_dihed_constr
 906       nnt=1
 907 #ifdef MPI
 908       if (me.eq.king) then
 909 #endif
 910        write (iout,'(a)') 'Boundaries in phi angle sampling:'
 911        do i=1,nres
 912          write (iout,'(a3,i5,2f10.1)')
 913      &   restyp(itype(i)),i,phibound(1,i)*rad2deg,phibound(2,i)*rad2deg
 914        enddo
 915 #ifdef MP
 916       endif
 917 #endif
 918       nct=nres
 919 cd      print *,'NNT=',NNT,' NCT=',NCT
 920       if (itype(1).eq.ntyp1) nnt=2
 921       if (itype(nres).eq.ntyp1) nct=nct-1
 922       if (pdbref) then
 923         if(me.eq.king.or..not.out1file)
 924      &   write (iout,'(a,i3)') 'nsup=',nsup
 925         nstart_seq=nnt
 926         if (nsup.le.(nct-nnt+1)) then
 927           do i=0,nct-nnt+1-nsup
 928             if (seq_comp(itype(nnt+i),itype_pdb(nstart_sup),nsup)) then
 929               nstart_seq=nnt+i
 930               goto 111
 931             endif
 932           enddo
 933           write (iout,'(a)')
 934      &            'Error - sequences to be superposed do not match.'
 935           stop
 936         else
 937           do i=0,nsup-(nct-nnt+1)
 938             if (seq_comp(itype(nnt),itype_pdb(nstart_sup+i),nct-nnt+1))
 939      &      then
 940               nstart_sup=nstart_sup+i
 941               nsup=nct-nnt+1
 942               goto 111
 943             endif
 944           enddo
 945           write (iout,'(a)')
 946      &            'Error - sequences to be superposed do not match.'
 947         endif
 948   111   continue
 949         if (nsup.eq.0) nsup=nct-nnt
 950         if (nstart_sup.eq.0) nstart_sup=nnt
 951         if (nstart_seq.eq.0) nstart_seq=nnt
 952         if(me.eq.king.or..not.out1file)
 953      &   write (iout,*) 'nsup=',nsup,' nstart_sup=',nstart_sup,
 954      &                 ' nstart_seq=',nstart_seq
 955       endif
 956 c--- Zscore rms -------
 957       if (nz_start.eq.0) nz_start=nnt
 958       if (nz_end.eq.0 .and. nsup.gt.0) then
 959         nz_end=nnt+nsup-1
 960       else if (nz_end.eq.0) then
 961         nz_end=nct
 962       endif
 963       if(me.eq.king.or..not.out1file)then
 964        write (iout,*) 'NZ_START=',nz_start,' NZ_END=',nz_end
 965        write (iout,*) 'IZ_SC=',iz_sc
 966       endif
 967 c----------------------
 968       call init_int_table
 969       if (refstr) then
 970         if (.not.pdbref) then
 971           call read_angles(inp,*38)
 972           goto 39
 973    38     write (iout,'(a)') 'Error reading reference structure.'
 974 #ifdef MPI
 975           call MPI_Finalize(MPI_COMM_WORLD,IERROR)
 976           stop 'Error reading reference structure'
 977 #endif
 978    39     call chainbuild
 979           call setup_var
 980 czscore          call geom_to_var(nvar,coord_exp_zs(1,1))
 981           nstart_sup=nnt
 982           nstart_seq=nnt
 983           nsup=nct-nnt+1
 984           kkk=1
 985           do i=1,2*nres
 986             do j=1,3
 987               cref(j,i,kkk)=c(j,i)
 988             enddo
 989           enddo
 990           call contact(.true.,ncont_ref,icont_ref,co)
 991         endif
 992         endif
 993         print *, "A TU"
 994         write (iout,*) "constr_dist",constr_dist,nstart_sup,nsup
 995         call flush(iout)
 996         if (constr_dist.gt.0) call read_dist_constr
 997         write (iout,*) "After read_dist_constr nhpb",nhpb
 998         call hpb_partition
 999         if(me.eq.king.or..not.out1file)
1000      &   write (iout,*) 'Contact order:',co
1001         if (pdbref) then
1002         if(me.eq.king.or..not.out1file)
1003      &   write (2,*) 'Shifting contacts:',nstart_seq,nstart_sup
1004         do i=1,ncont_ref
1005           do j=1,2
1006             icont_ref(j,i)=icont_ref(j,i)+nstart_seq-nstart_sup
1007           enddo
1008           if(me.eq.king.or..not.out1file)
1009      &     write (2,*) i,' ',restyp(itype(icont_ref(1,i))),' ',
1010      &     icont_ref(1,i),' ',
1011      &     restyp(itype(icont_ref(2,i))),' ',icont_ref(2,i)
1012         enddo
1013         endif
1014 C      endif
1015       if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1016      &    .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1017      &    modecalc.ne.10) then
1018 C If input structure hasn't been supplied from the PDB file read or generate
1019 C initial geometry.
1020         if (iranconf.eq.0 .and. .not. extconf) then
1021           if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1022      &     write (iout,'(a)') 'Initial geometry will be read in.'
1023           if (read_cart) then
1024             read(inp,'(8f10.5)',end=36,err=36)
1025      &       ((c(l,k),l=1,3),k=1,nres),
1026      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
1027             write (iout,*) "Exit READ_CART"
1028             write (iout,'(8f10.5)')
1029      &       ((c(l,k),l=1,3),k=1,nres),
1030      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
1031             call int_from_cart1(.true.)
1032             write (iout,*) "Finish INT_TO_CART"
1033             do i=1,nres-1
1034               do j=1,3
1035                 dc(j,i)=c(j,i+1)-c(j,i)
1036                 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1037               enddo
1038             enddo
1039             do i=nnt,nct
1040               if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1041                 do j=1,3
1042                   dc(j,i+nres)=c(j,i+nres)-c(j,i)
1043                   dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1044                 enddo
1045               endif
1046             enddo
1047             return
1048           else
1049             call read_angles(inp,*36)
1050           endif
1051           goto 37
1052    36     write (iout,'(a)') 'Error reading angle file.'
1053 #ifdef MPI
1054           call mpi_finalize( MPI_COMM_WORLD,IERR )
1055 #endif
1056           stop 'Error reading angle file.'
1057    37     continue
1058         else if (extconf) then
1059          if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1060      &    write (iout,'(a)') 'Extended chain initial geometry.'
1061          do i=3,nres
1062           theta(i)=90d0*deg2rad
1063          enddo
1064          do i=4,nres
1065           phi(i)=180d0*deg2rad
1066          enddo
1067          do i=2,nres-1
1068           alph(i)=110d0*deg2rad
1069          enddo
1070          do i=2,nres-1
1071           omeg(i)=-120d0*deg2rad
1072           if (itype(i).le.0) omeg(i)=-omeg(i)
1073          enddo
1074         else
1075           if(me.eq.king.or..not.out1file)
1076      &     write (iout,'(a)') 'Random-generated initial geometry.'
1077
1078
1079 #ifdef MPI
1080           if (me.eq.king  .or. fg_rank.eq.0 .and. (
1081      &           modecalc.eq.12 .or. modecalc.eq.14) ) then
1082 #endif
1083             do itrial=1,100
1084               itmp=1
1085               call gen_rand_conf(itmp,*30)
1086               goto 40
1087    30         write (iout,*) 'Failed to generate random conformation',
1088      &          ', itrial=',itrial
1089               write (*,*) 'Processor:',me,
1090      &          ' Failed to generate random conformation',
1091      &          ' itrial=',itrial
1092               call intout
1093
1094 #ifdef AIX
1095               call flush_(iout)
1096 #else
1097               call flush(iout)
1098 #endif
1099             enddo
1100             write (iout,'(a,i3,a)') 'Processor:',me,
1101      &        ' error in generating random conformation.'
1102             write (*,'(a,i3,a)') 'Processor:',me,
1103      &        ' error in generating random conformation.'
1104             call flush(iout)
1105 #ifdef MPI
1106             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1107    40       continue
1108           endif
1109 #else
1110           do itrial=1,100
1111             itmp=1
1112             call gen_rand_conf(itmp,*30)
1113             goto 40
1114    30       write (iout,*) 'Failed to generate random conformation',
1115      &        ', itrial=',itrial
1116             write (*,*) 'Failed to generate random conformation',
1117      &        ', itrial=',itrial
1118           enddo
1119           write (iout,'(a,i3,a)') 'Processor:',me,
1120      &      ' error in generating random conformation.'
1121           write (*,'(a,i3,a)') 'Processor:',me,
1122      &      ' error in generating random conformation.'
1123           stop
1124    40     continue
1125 #endif
1126         endif
1127       elseif (modecalc.eq.4) then
1128         read (inp,'(a)') intinname
1129         open (intin,file=intinname,status='old',err=333)
1130         if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1131      &  write (iout,'(a)') 'intinname',intinname
1132         write (*,'(a)') 'Processor',myrank,' intinname',intinname
1133         goto 334
1134   333   write (iout,'(2a)') 'Error opening angle file ',intinname
1135 #ifdef MPI
1136         call MPI_Finalize(MPI_COMM_WORLD,IERR)
1137 #endif
1138         stop 'Error opening angle file.'
1139   334   continue
1140
1141       endif
1142 C Generate distance constraints, if the PDB structure is to be regularized.
1143       if (nthread.gt.0) then
1144         call read_threadbase
1145       endif
1146       call setup_var
1147       if (me.eq.king .or. .not. out1file)
1148      & call intout
1149       if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1150         write (iout,'(/a,i3,a)')
1151      &  'The chain contains',ns,' disulfide-bridging cysteines.'
1152         write (iout,'(20i4)') (iss(i),i=1,ns)
1153        if (dyn_ss) then
1154           write(iout,*)"Running with dynamic disulfide-bond formation"
1155        else
1156         write (iout,'(/a/)') 'Pre-formed links are:'
1157         do i=1,nss
1158           i1=ihpb(i)-nres
1159           i2=jhpb(i)-nres
1160           it1=itype(i1)
1161           it2=itype(i2)
1162           write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1163      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1164      &    ebr,forcon(i)
1165         enddo
1166         write (iout,'(a)')
1167        endif
1168       endif
1169       if (ns.gt.0.and.dyn_ss) then
1170           do i=nss+1,nhpb
1171             ihpb(i-nss)=ihpb(i)
1172             jhpb(i-nss)=jhpb(i)
1173             forcon(i-nss)=forcon(i)
1174             dhpb(i-nss)=dhpb(i)
1175           enddo
1176           nhpb=nhpb-nss
1177           nss=0
1178           call hpb_partition
1179           do i=1,ns
1180             dyn_ss_mask(iss(i))=.true.
1181           enddo
1182       endif
1183       if (i2ndstr.gt.0) call secstrp2dihc
1184 c      call geom_to_var(nvar,x)
1185 c      call etotal(energia(0))
1186 c      call enerprint(energia(0))
1187 c      call briefout(0,etot)
1188 c      stop
1189 cd    write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1190 cd    write (iout,'(a)') 'Variable list:'
1191 cd    write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1192 #ifdef MPI
1193       if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1194      &  write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1195      &  'Processor',myrank,': end reading molecular data.'
1196 #endif
1197       print *,"A TU?"
1198       return
1199       end
1200 c--------------------------------------------------------------------------
1201       logical function seq_comp(itypea,itypeb,length)
1202       implicit none
1203       integer length,itypea(length),itypeb(length)
1204       integer i
1205       do i=1,length
1206         if (itypea(i).ne.itypeb(i)) then
1207           seq_comp=.false.
1208           return
1209         endif
1210       enddo
1211       seq_comp=.true.
1212       return
1213       end
1214 c-----------------------------------------------------------------------------
1215       subroutine read_bridge
1216 C Read information about disulfide bridges.
1217       implicit real*8 (a-h,o-z)
1218       include 'DIMENSIONS'
1219 #ifdef MPI
1220       include 'mpif.h'
1221 #endif
1222       include 'COMMON.IOUNITS'
1223       include 'COMMON.GEO'
1224       include 'COMMON.VAR'
1225       include 'COMMON.INTERACT'
1226       include 'COMMON.LOCAL'
1227       include 'COMMON.NAMES'
1228       include 'COMMON.CHAIN'
1229       include 'COMMON.FFIELD'
1230       include 'COMMON.SBRIDGE'
1231       include 'COMMON.HEADER'
1232       include 'COMMON.CONTROL'
1233       include 'COMMON.DBASE'
1234       include 'COMMON.THREAD'
1235       include 'COMMON.TIME1'
1236       include 'COMMON.SETUP'
1237 C Read bridging residues.
1238       read (inp,*) ns,(iss(i),i=1,ns)
1239       print *,'ns=',ns
1240       if(me.eq.king.or..not.out1file)
1241      &  write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1242 C Check whether the specified bridging residues are cystines.
1243       do i=1,ns
1244         if (itype(iss(i)).ne.1) then
1245           if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1246      &   'Do you REALLY think that the residue ',
1247      &    restyp(itype(iss(i))),i,
1248      &   ' can form a disulfide bridge?!!!'
1249           write (*,'(2a,i3,a)')
1250      &   'Do you REALLY think that the residue ',
1251      &    restyp(itype(iss(i))),i,
1252      &   ' can form a disulfide bridge?!!!'
1253 #ifdef MPI
1254          call MPI_Finalize(MPI_COMM_WORLD,ierror)
1255          stop
1256 #endif
1257         endif
1258       enddo
1259 C Read preformed bridges.
1260       if (ns.gt.0) then
1261       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1262       if(fg_rank.eq.0)
1263      & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1264       if (nss.gt.0) then
1265         nhpb=nss
1266 C Check if the residues involved in bridges are in the specified list of
1267 C bridging residues.
1268         do i=1,nss
1269           do j=1,i-1
1270             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1271      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1272               write (iout,'(a,i3,a)') 'Disulfide pair',i,
1273      &      ' contains residues present in other pairs.'
1274               write (*,'(a,i3,a)') 'Disulfide pair',i,
1275      &      ' contains residues present in other pairs.'
1276 #ifdef MPI
1277               call MPI_Finalize(MPI_COMM_WORLD,ierror)
1278               stop
1279 #endif
1280             endif
1281           enddo
1282           do j=1,ns
1283             if (ihpb(i).eq.iss(j)) goto 10
1284           enddo
1285           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1286    10     continue
1287           do j=1,ns
1288             if (jhpb(i).eq.iss(j)) goto 20
1289           enddo
1290           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1291    20     continue
1292           dhpb(i)=dbr
1293           forcon(i)=fbr
1294         enddo
1295         do i=1,nss
1296           ihpb(i)=ihpb(i)+nres
1297           jhpb(i)=jhpb(i)+nres
1298         enddo
1299       endif
1300       endif
1301       return
1302       end
1303 c----------------------------------------------------------------------------
1304       subroutine read_x(kanal,*)
1305       implicit real*8 (a-h,o-z)
1306       include 'DIMENSIONS'
1307       include 'COMMON.GEO'
1308       include 'COMMON.VAR'
1309       include 'COMMON.CHAIN'
1310       include 'COMMON.IOUNITS'
1311       include 'COMMON.CONTROL'
1312       include 'COMMON.LOCAL'
1313       include 'COMMON.INTERACT'
1314 c Read coordinates from input
1315 c
1316       read(kanal,'(8f10.5)',end=10,err=10)
1317      &  ((c(l,k),l=1,3),k=1,nres),
1318      &  ((c(l,k+nres),l=1,3),k=nnt,nct)
1319       do j=1,3
1320         c(j,nres+1)=c(j,1)
1321         c(j,2*nres)=c(j,nres)
1322       enddo
1323       call int_from_cart1(.false.)
1324       do i=1,nres-1
1325         do j=1,3
1326           dc(j,i)=c(j,i+1)-c(j,i)
1327           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1328         enddo
1329       enddo
1330       do i=nnt,nct
1331         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1332           do j=1,3
1333             dc(j,i+nres)=c(j,i+nres)-c(j,i)
1334             dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1335           enddo
1336         endif
1337       enddo
1338
1339       return
1340    10 return1
1341       end
1342 c----------------------------------------------------------------------------
1343       subroutine read_threadbase
1344       implicit real*8 (a-h,o-z)
1345       include 'DIMENSIONS'
1346       include 'COMMON.IOUNITS'
1347       include 'COMMON.GEO'
1348       include 'COMMON.VAR'
1349       include 'COMMON.INTERACT'
1350       include 'COMMON.LOCAL'
1351       include 'COMMON.NAMES'
1352       include 'COMMON.CHAIN'
1353       include 'COMMON.FFIELD'
1354       include 'COMMON.SBRIDGE'
1355       include 'COMMON.HEADER'
1356       include 'COMMON.CONTROL'
1357       include 'COMMON.DBASE'
1358       include 'COMMON.THREAD'
1359       include 'COMMON.TIME1'
1360 C Read pattern database for threading.
1361       read (icbase,*) nseq
1362       do i=1,nseq
1363         read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1364      &   nres_base(2,i),nres_base(3,i)
1365         read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1366      &   nres_base(1,i))
1367 c       write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1368 c    &   nres_base(2,i),nres_base(3,i)
1369 c       write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1370 c    &   nres_base(1,i))
1371       enddo
1372       close (icbase)
1373       if (weidis.eq.0.0D0) weidis=0.1D0
1374       do i=nnt,nct
1375         do j=i+2,nct
1376           nhpb=nhpb+1
1377           ihpb(nhpb)=i
1378           jhpb(nhpb)=j
1379           forcon(nhpb)=weidis
1380         enddo
1381       enddo
1382       read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1383       write (iout,'(a,i5)') 'nexcl: ',nexcl
1384       write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1385       return
1386       end
1387 c------------------------------------------------------------------------------
1388       subroutine setup_var
1389       implicit real*8 (a-h,o-z)
1390       include 'DIMENSIONS'
1391       include 'COMMON.IOUNITS'
1392       include 'COMMON.GEO'
1393       include 'COMMON.VAR'
1394       include 'COMMON.INTERACT'
1395       include 'COMMON.LOCAL'
1396       include 'COMMON.NAMES'
1397       include 'COMMON.CHAIN'
1398       include 'COMMON.FFIELD'
1399       include 'COMMON.SBRIDGE'
1400       include 'COMMON.HEADER'
1401       include 'COMMON.CONTROL'
1402       include 'COMMON.DBASE'
1403       include 'COMMON.THREAD'
1404       include 'COMMON.TIME1'
1405 C Set up variable list.
1406       ntheta=nres-2
1407       nphi=nres-3
1408       nvar=ntheta+nphi
1409       nside=0
1410       do i=2,nres-1
1411         if (itype(i).ne.10 .and. itype(i).ne.ntyp1) then
1412           nside=nside+1
1413           ialph(i,1)=nvar+nside
1414           ialph(nside,2)=i
1415         endif
1416       enddo
1417       if (indphi.gt.0) then
1418         nvar=nphi
1419       else if (indback.gt.0) then
1420         nvar=nphi+ntheta
1421       else
1422         nvar=nvar+2*nside
1423       endif
1424 cd    write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1425       return
1426       end
1427 c----------------------------------------------------------------------------
1428       subroutine gen_dist_constr
1429 C Generate CA distance constraints.
1430       implicit real*8 (a-h,o-z)
1431       include 'DIMENSIONS'
1432       include 'COMMON.IOUNITS'
1433       include 'COMMON.GEO'
1434       include 'COMMON.VAR'
1435       include 'COMMON.INTERACT'
1436       include 'COMMON.LOCAL'
1437       include 'COMMON.NAMES'
1438       include 'COMMON.CHAIN'
1439       include 'COMMON.FFIELD'
1440       include 'COMMON.SBRIDGE'
1441       include 'COMMON.HEADER'
1442       include 'COMMON.CONTROL'
1443       include 'COMMON.DBASE'
1444       include 'COMMON.THREAD'
1445       include 'COMMON.TIME1'
1446       dimension itype_pdb(maxres)
1447       common /pizda/ itype_pdb
1448       character*2 iden
1449 cd      print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1450 cd      write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1451 cd     & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1452 cd     & ' nsup',nsup
1453       do i=nstart_sup,nstart_sup+nsup-1
1454 cd      write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1455 cd     &    ' seq_pdb', restyp(itype_pdb(i))
1456         do j=i+2,nstart_sup+nsup-1
1457           nhpb=nhpb+1
1458           ihpb(nhpb)=i+nstart_seq-nstart_sup
1459           jhpb(nhpb)=j+nstart_seq-nstart_sup
1460           forcon(nhpb)=weidis
1461           dhpb(nhpb)=dist(i,j)
1462         enddo
1463       enddo
1464 cd      write (iout,'(a)') 'Distance constraints:'
1465 cd      do i=nss+1,nhpb
1466 cd        ii=ihpb(i)
1467 cd        jj=jhpb(i)
1468 cd        iden='CA'
1469 cd        if (ii.gt.nres) then
1470 cd          iden='SC'
1471 cd          ii=ii-nres
1472 cd          jj=jj-nres
1473 cd        endif
1474 cd        write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1475 cd     &  restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1476 cd     &  dhpb(i),forcon(i)
1477 cd      enddo
1478       return
1479       end
1480 c----------------------------------------------------------------------------
1481       subroutine map_read
1482       implicit real*8 (a-h,o-z)
1483       include 'DIMENSIONS'
1484       include 'COMMON.MAP'
1485       include 'COMMON.IOUNITS'
1486       character*3 angid(4) /'THE','PHI','ALP','OME'/
1487       character*80 mapcard,ucase
1488       do imap=1,nmap
1489         read (inp,'(a)') mapcard
1490         mapcard=ucase(mapcard)
1491         if (index(mapcard,'PHI').gt.0) then
1492           kang(imap)=1
1493         else if (index(mapcard,'THE').gt.0) then
1494           kang(imap)=2
1495         else if (index(mapcard,'ALP').gt.0) then
1496           kang(imap)=3
1497         else if (index(mapcard,'OME').gt.0) then
1498           kang(imap)=4
1499         else
1500           write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1501           stop 'Error - illegal variable spec in MAP card.'
1502         endif
1503         call readi (mapcard,'RES1',res1(imap),0)
1504         call readi (mapcard,'RES2',res2(imap),0)
1505         if (res1(imap).eq.0) then
1506           res1(imap)=res2(imap)
1507         else if (res2(imap).eq.0) then
1508           res2(imap)=res1(imap)
1509         endif
1510         if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1511           write (iout,'(a)')
1512      &    'Error - illegal definition of variable group in MAP.'
1513           stop 'Error - illegal definition of variable group in MAP.'
1514         endif
1515         call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1516         call reada(mapcard,'TO',ang_to(imap),0.0D0)
1517         call readi(mapcard,'NSTEP',nstep(imap),0)
1518         if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1519           write (iout,'(a)')
1520      &     'Illegal boundary and/or step size specification in MAP.'
1521           stop 'Illegal boundary and/or step size specification in MAP.'
1522         endif
1523       enddo ! imap
1524       return
1525       end
1526 c----------------------------------------------------------------------------
1527       subroutine csaread
1528       implicit real*8 (a-h,o-z)
1529       include 'DIMENSIONS'
1530       include 'COMMON.IOUNITS'
1531       include 'COMMON.GEO'
1532       include 'COMMON.CSA'
1533       include 'COMMON.BANK'
1534       include 'COMMON.CONTROL'
1535       character*80 ucase
1536       character*620 mcmcard
1537       call card_concat(mcmcard)
1538
1539       call readi(mcmcard,'NCONF',nconf,50)
1540       call readi(mcmcard,'NADD',nadd,0)
1541       call readi(mcmcard,'JSTART',jstart,1)
1542       call readi(mcmcard,'JEND',jend,1)
1543       call readi(mcmcard,'NSTMAX',nstmax,500000)
1544       call readi(mcmcard,'N0',n0,1)
1545       call readi(mcmcard,'N1',n1,6)
1546       call readi(mcmcard,'N2',n2,4)
1547       call readi(mcmcard,'N3',n3,0)
1548       call readi(mcmcard,'N4',n4,0)
1549       call readi(mcmcard,'N5',n5,0)
1550       call readi(mcmcard,'N6',n6,10)
1551       call readi(mcmcard,'N7',n7,0)
1552       call readi(mcmcard,'N8',n8,0)
1553       call readi(mcmcard,'N9',n9,0)
1554       call readi(mcmcard,'N14',n14,0)
1555       call readi(mcmcard,'N15',n15,0)
1556       call readi(mcmcard,'N16',n16,0)
1557       call readi(mcmcard,'N17',n17,0)
1558       call readi(mcmcard,'N18',n18,0)
1559
1560       vdisulf=(index(mcmcard,'DYNSS').gt.0)
1561
1562       call readi(mcmcard,'NDIFF',ndiff,2)
1563       call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1564       call readi(mcmcard,'IS1',is1,1)
1565       call readi(mcmcard,'IS2',is2,8)
1566       call readi(mcmcard,'NRAN0',nran0,4)
1567       call readi(mcmcard,'NRAN1',nran1,2)
1568       call readi(mcmcard,'IRR',irr,1)
1569       call readi(mcmcard,'NSEED',nseed,20)
1570       call readi(mcmcard,'NTOTAL',ntotal,10000)
1571       call reada(mcmcard,'CUT1',cut1,2.0d0)
1572       call reada(mcmcard,'CUT2',cut2,5.0d0)
1573       call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1574       call readi(mcmcard,'ICMAX',icmax,3)
1575       call readi(mcmcard,'IRESTART',irestart,0)
1576 c!bankt      call readi(mcmcard,'NBANKTM',ntbankm,0)
1577       ntbankm=0
1578 c!bankt
1579       call reada(mcmcard,'DELE',dele,20.0d0)
1580       call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1581       call readi(mcmcard,'IREF',iref,0)
1582       call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1583       call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1584       call readi(mcmcard,'NCONF_IN',nconf_in,0)
1585       call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1586       write (iout,*) "NCONF_IN",nconf_in
1587       return
1588       end
1589 c----------------------------------------------------------------------------
1590 cfmc      subroutine mcmfread
1591 cfmc      implicit real*8 (a-h,o-z)
1592 cfmc      include 'DIMENSIONS'
1593 cfmc      include 'COMMON.MCMF'
1594 cfmc      include 'COMMON.IOUNITS'
1595 cfmc      include 'COMMON.GEO'
1596 cfmc      character*80 ucase
1597 cfmc      character*620 mcmcard
1598 cfmc      call card_concat(mcmcard)
1599 cfmc
1600 cfmc      call readi(mcmcard,'MAXRANT',maxrant,1000)
1601 cfmc      write(iout,*)'MAXRANT=',maxrant
1602 cfmc      call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1603 cfmc      write(iout,*)'MAXFAM=',maxfam
1604 cfmc      call readi(mcmcard,'NNET1',nnet1,5)
1605 cfmc      write(iout,*)'NNET1=',nnet1
1606 cfmc      call readi(mcmcard,'NNET2',nnet2,4)
1607 cfmc      write(iout,*)'NNET2=',nnet2
1608 cfmc      call readi(mcmcard,'NNET3',nnet3,4)
1609 cfmc      write(iout,*)'NNET3=',nnet3
1610 cfmc      call readi(mcmcard,'ILASTT',ilastt,0)
1611 cfmc      write(iout,*)'ILASTT=',ilastt
1612 cfmc      call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1613 cfmc      write(iout,*)'MAXSTR=',maxstr
1614 cfmc      maxstr_f=maxstr/maxfam
1615 cfmc      write(iout,*)'MAXSTR_F=',maxstr_f
1616 cfmc      call readi(mcmcard,'NMCMF',nmcmf,10)
1617 cfmc      write(iout,*)'NMCMF=',nmcmf
1618 cfmc      call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1619 cfmc      write(iout,*)'IFOCUS=',ifocus
1620 cfmc      call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1621 cfmc      write(iout,*)'NLOCMCMF=',nlocmcmf
1622 cfmc      call readi(mcmcard,'INTPRT',intprt,1000)
1623 cfmc      write(iout,*)'INTPRT=',intprt
1624 cfmc      call readi(mcmcard,'IPRT',iprt,100)
1625 cfmc      write(iout,*)'IPRT=',iprt
1626 cfmc      call readi(mcmcard,'IMAXTR',imaxtr,100)
1627 cfmc      write(iout,*)'IMAXTR=',imaxtr
1628 cfmc      call readi(mcmcard,'MAXEVEN',maxeven,1000)
1629 cfmc      write(iout,*)'MAXEVEN=',maxeven
1630 cfmc      call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1631 cfmc      write(iout,*)'MAXEVEN1=',maxeven1
1632 cfmc      call readi(mcmcard,'INIMIN',inimin,200)
1633 cfmc      write(iout,*)'INIMIN=',inimin
1634 cfmc      call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1635 cfmc      write(iout,*)'NSTEPMCMF=',nstepmcmf
1636 cfmc      call readi(mcmcard,'NTHREAD',nthread,5)
1637 cfmc      write(iout,*)'NTHREAD=',nthread
1638 cfmc      call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1639 cfmc      write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1640 cfmc      call readi(mcmcard,'MAXPERT',maxpert,9)
1641 cfmc      write(iout,*)'MAXPERT=',maxpert
1642 cfmc      call readi(mcmcard,'IRMSD',irmsd,1)
1643 cfmc      write(iout,*)'IRMSD=',irmsd
1644 cfmc      call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1645 cfmc      write(iout,*)'DENEMIN=',denemin
1646 cfmc      call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1647 cfmc      write(iout,*)'RCUT1S=',rcut1s
1648 cfmc      call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1649 cfmc      write(iout,*)'RCUT1E=',rcut1e
1650 cfmc      call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1651 cfmc      write(iout,*)'RCUT2S=',rcut2s
1652 cfmc      call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1653 cfmc      write(iout,*)'RCUT2E=',rcut2e
1654 cfmc      call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1655 cfmc      write(iout,*)'DPERT1=',d_pert1
1656 cfmc      call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1657 cfmc      write(iout,*)'DPERT1A=',d_pert1a
1658 cfmc      call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1659 cfmc      write(iout,*)'DPERT2=',d_pert2
1660 cfmc      call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1661 cfmc      write(iout,*)'DPERT2A=',d_pert2a
1662 cfmc      call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1663 cfmc      write(iout,*)'DPERT2B=',d_pert2b
1664 cfmc      call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1665 cfmc      write(iout,*)'DPERT2C=',d_pert2c
1666 cfmc      d_pert1=deg2rad*d_pert1
1667 cfmc      d_pert1a=deg2rad*d_pert1a
1668 cfmc      d_pert2=deg2rad*d_pert2
1669 cfmc      d_pert2a=deg2rad*d_pert2a
1670 cfmc      d_pert2b=deg2rad*d_pert2b
1671 cfmc      d_pert2c=deg2rad*d_pert2c
1672 cfmc      call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1673 cfmc      write(iout,*)'KT_MCMF1=',kt_mcmf1
1674 cfmc      call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1675 cfmc      write(iout,*)'KT_MCMF2=',kt_mcmf2
1676 cfmc      call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1677 cfmc      write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1678 cfmc      call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1679 cfmc      write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1680 cfmc      call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1681 cfmc      write(iout,*)'RCUTINI=',rcutini
1682 cfmc      call reada(mcmcard,'GRAT',grat,0.5D0)
1683 cfmc      write(iout,*)'GRAT=',grat
1684 cfmc      call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1685 cfmc      write(iout,*)'BIAS_MCMF=',bias_mcmf
1686 cfmc
1687 cfmc      return
1688 cfmc      end
1689 c----------------------------------------------------------------------------
1690       subroutine mcmread
1691       implicit real*8 (a-h,o-z)
1692       include 'DIMENSIONS'
1693       include 'COMMON.MCM'
1694       include 'COMMON.MCE'
1695       include 'COMMON.IOUNITS'
1696       character*80 ucase
1697       character*320 mcmcard
1698       call card_concat(mcmcard)
1699       call readi(mcmcard,'MAXACC',maxacc,100)
1700       call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1701       call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1702       call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1703       call readi(mcmcard,'MAXREPM',maxrepm,200)
1704       call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1705       call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1706       call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1707       call reada(mcmcard,'E_UP',e_up,5.0D0)
1708       call reada(mcmcard,'DELTE',delte,0.1D0)
1709       call readi(mcmcard,'NSWEEP',nsweep,5)
1710       call readi(mcmcard,'NSTEPH',nsteph,0)
1711       call readi(mcmcard,'NSTEPC',nstepc,0)
1712       call reada(mcmcard,'TMIN',tmin,298.0D0)
1713       call reada(mcmcard,'TMAX',tmax,298.0D0)
1714       call readi(mcmcard,'NWINDOW',nwindow,0)
1715       call readi(mcmcard,'PRINT_MC',print_mc,0)
1716       print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1717       print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1718       ent_read=(index(mcmcard,'ENT_READ').gt.0)
1719       call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1720       call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1721       call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1722       call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1723       call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1724       if (nwindow.gt.0) then
1725         read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1726         do i=1,nwindow
1727           winlen(i)=winend(i)-winstart(i)+1
1728         enddo
1729       endif
1730       if (tmax.lt.tmin) tmax=tmin
1731       if (tmax.eq.tmin) then
1732         nstepc=0
1733         nsteph=0
1734       endif
1735       if (nstepc.gt.0 .and. nsteph.gt.0) then
1736         tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1737         tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1738       endif
1739 C Probabilities of different move types
1740       sumpro_type(0)=0.0D0
1741       call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1742       call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1743       sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1744       call reada(mcmcard,'THETA'     ,sumpro_type(3),0.0d0)
1745       sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1746       call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1747       sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1748       do i=1,MaxMoveType
1749         print *,'i',i,' sumprotype',sumpro_type(i)
1750         sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1751         print *,'i',i,' sumprotype',sumpro_type(i)
1752       enddo
1753       return
1754       end
1755 c----------------------------------------------------------------------------
1756       subroutine read_minim
1757       implicit real*8 (a-h,o-z)
1758       include 'DIMENSIONS'
1759       include 'COMMON.MINIM'
1760       include 'COMMON.IOUNITS'
1761       character*80 ucase
1762       character*320 minimcard
1763       call card_concat(minimcard)
1764       call readi(minimcard,'MAXMIN',maxmin,2000)
1765       call readi(minimcard,'MAXFUN',maxfun,5000)
1766       call readi(minimcard,'MINMIN',minmin,maxmin)
1767       call readi(minimcard,'MINFUN',minfun,maxmin)
1768       call reada(minimcard,'TOLF',tolf,1.0D-2)
1769       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1770       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1771       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1772       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1773       write (iout,'(/80(1h*)/20x,a/80(1h*))')
1774      &         'Options in energy minimization:'
1775       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1776      & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1777      & 'MinMin:',MinMin,' MinFun:',MinFun,
1778      & ' TolF:',TolF,' RTolF:',RTolF
1779       return
1780       end
1781 c----------------------------------------------------------------------------
1782       subroutine read_angles(kanal,*)
1783       implicit real*8 (a-h,o-z)
1784       include 'DIMENSIONS'
1785       include 'COMMON.GEO'
1786       include 'COMMON.VAR'
1787       include 'COMMON.CHAIN'
1788       include 'COMMON.IOUNITS'
1789       include 'COMMON.CONTROL'
1790 c Read angles from input
1791 c
1792        read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1793        read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1794        read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1795        read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1796
1797        do i=1,nres
1798 c 9/7/01 avoid 180 deg valence angle
1799         if (theta(i).gt.179.99d0) theta(i)=179.99d0
1800 c
1801         theta(i)=deg2rad*theta(i)
1802         phi(i)=deg2rad*phi(i)
1803         alph(i)=deg2rad*alph(i)
1804         omeg(i)=deg2rad*omeg(i)
1805        enddo
1806       return
1807    10 return1
1808       end
1809 c----------------------------------------------------------------------------
1810       subroutine reada(rekord,lancuch,wartosc,default)
1811       implicit none
1812       character*(*) rekord,lancuch
1813       double precision wartosc,default
1814       integer ilen,iread
1815       external ilen
1816       iread=index(rekord,lancuch)
1817       if (iread.eq.0) then
1818         wartosc=default
1819         return
1820       endif
1821       iread=iread+ilen(lancuch)+1
1822       read (rekord(iread:),*,err=10,end=10) wartosc
1823       return
1824   10  wartosc=default
1825       return
1826       end
1827 c----------------------------------------------------------------------------
1828       subroutine readi(rekord,lancuch,wartosc,default)
1829       implicit none
1830       character*(*) rekord,lancuch
1831       integer wartosc,default
1832       integer ilen,iread
1833       external ilen
1834       iread=index(rekord,lancuch)
1835       if (iread.eq.0) then
1836         wartosc=default
1837         return
1838       endif
1839       iread=iread+ilen(lancuch)+1
1840       read (rekord(iread:),*,err=10,end=10) wartosc
1841       return
1842   10  wartosc=default
1843       return
1844       end
1845 c----------------------------------------------------------------------------
1846       subroutine multreadi(rekord,lancuch,tablica,dim,default)
1847       implicit none
1848       integer dim,i
1849       integer tablica(dim),default
1850       character*(*) rekord,lancuch
1851       character*80 aux
1852       integer ilen,iread
1853       external ilen
1854       do i=1,dim
1855         tablica(i)=default
1856       enddo
1857       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1858       if (iread.eq.0) return
1859       iread=iread+ilen(lancuch)+1
1860       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1861    10 return
1862       end
1863 c----------------------------------------------------------------------------
1864       subroutine multreada(rekord,lancuch,tablica,dim,default)
1865       implicit none
1866       integer dim,i
1867       double precision tablica(dim),default
1868       character*(*) rekord,lancuch
1869       character*80 aux
1870       integer ilen,iread
1871       external ilen
1872       do i=1,dim
1873         tablica(i)=default
1874       enddo
1875       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1876       if (iread.eq.0) return
1877       iread=iread+ilen(lancuch)+1
1878       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
1879    10 return
1880       end
1881 c----------------------------------------------------------------------------
1882       subroutine openunits
1883       implicit real*8 (a-h,o-z)
1884       include 'DIMENSIONS'
1885 #ifdef MPI
1886       include 'mpif.h'
1887       character*16 form,nodename
1888       integer nodelen
1889 #endif
1890       include 'COMMON.SETUP'
1891       include 'COMMON.IOUNITS'
1892       include 'COMMON.MD'
1893       include 'COMMON.CONTROL'
1894       integer lenpre,lenpot,ilen,lentmp
1895       external ilen
1896       character*3 out1file_text,ucase
1897       character*3 ll
1898       external ucase
1899 c      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
1900       call getenv_loc("PREFIX",prefix)
1901       pref_orig = prefix
1902       call getenv_loc("POT",pot)
1903       call getenv_loc("DIRTMP",tmpdir)
1904       call getenv_loc("CURDIR",curdir)
1905       call getenv_loc("OUT1FILE",out1file_text)
1906 c      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
1907       out1file_text=ucase(out1file_text)
1908       if (out1file_text(1:1).eq."Y") then
1909         out1file=.true.
1910       else
1911         out1file=fg_rank.gt.0
1912       endif
1913       lenpre=ilen(prefix)
1914       lenpot=ilen(pot)
1915       lentmp=ilen(tmpdir)
1916       if (lentmp.gt.0) then
1917           write (*,'(80(1h!))')
1918           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
1919           write (*,'(80(1h!))')
1920           write (*,*)"All output files will be on node /tmp directory."
1921 #ifdef MPI
1922         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
1923         if (me.eq.king) then
1924           write (*,*) "The master node is ",nodename
1925         else if (fg_rank.eq.0) then
1926           write (*,*) "I am the CG slave node ",nodename
1927         else
1928           write (*,*) "I am the FG slave node ",nodename
1929         endif
1930 #endif
1931         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
1932         lenpre = lentmp+lenpre+1
1933       endif
1934       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
1935 C Get the names and open the input files
1936 #if defined(WINIFL) || defined(WINPGI)
1937       open(1,file=pref_orig(:ilen(pref_orig))//
1938      &  '.inp',status='old',readonly,shared)
1939        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1940 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1941 C Get parameter filenames and open the parameter files.
1942       call getenv_loc('BONDPAR',bondname)
1943       open (ibond,file=bondname,status='old',readonly,shared)
1944       call getenv_loc('THETPAR',thetname)
1945       open (ithep,file=thetname,status='old',readonly,shared)
1946       call getenv_loc('ROTPAR',rotname)
1947       open (irotam,file=rotname,status='old',readonly,shared)
1948       call getenv_loc('TORPAR',torname)
1949       open (itorp,file=torname,status='old',readonly,shared)
1950       call getenv_loc('TORDPAR',tordname)
1951       open (itordp,file=tordname,status='old',readonly,shared)
1952       call getenv_loc('FOURIER',fouriername)
1953       open (ifourier,file=fouriername,status='old',readonly,shared)
1954       call getenv_loc('ELEPAR',elename)
1955       open (ielep,file=elename,status='old',readonly,shared)
1956       call getenv_loc('SIDEPAR',sidename)
1957       open (isidep,file=sidename,status='old',readonly,shared)
1958 #elif (defined CRAY) || (defined AIX)
1959       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
1960      &  action='read')
1961 c      print *,"Processor",myrank," opened file 1"
1962       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1963 c      print *,"Processor",myrank," opened file 9"
1964 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1965 C Get parameter filenames and open the parameter files.
1966       call getenv_loc('BONDPAR',bondname)
1967       open (ibond,file=bondname,status='old',action='read')
1968 c      print *,"Processor",myrank," opened file IBOND"
1969       call getenv_loc('THETPAR',thetname)
1970       open (ithep,file=thetname,status='old',action='read')
1971 c      print *,"Processor",myrank," opened file ITHEP"
1972       call getenv_loc('ROTPAR',rotname)
1973       open (irotam,file=rotname,status='old',action='read')
1974 c      print *,"Processor",myrank," opened file IROTAM"
1975       call getenv_loc('TORPAR',torname)
1976       open (itorp,file=torname,status='old',action='read')
1977 c      print *,"Processor",myrank," opened file ITORP"
1978       call getenv_loc('TORDPAR',tordname)
1979       open (itordp,file=tordname,status='old',action='read')
1980 c      print *,"Processor",myrank," opened file ITORDP"
1981       call getenv_loc('SCCORPAR',sccorname)
1982       open (isccor,file=sccorname,status='old',action='read')
1983 c      print *,"Processor",myrank," opened file ISCCOR"
1984       call getenv_loc('FOURIER',fouriername)
1985       open (ifourier,file=fouriername,status='old',action='read')
1986 c      print *,"Processor",myrank," opened file IFOURIER"
1987       call getenv_loc('ELEPAR',elename)
1988       open (ielep,file=elename,status='old',action='read')
1989 c      print *,"Processor",myrank," opened file IELEP"
1990       call getenv_loc('SIDEPAR',sidename)
1991       open (isidep,file=sidename,status='old',action='read')
1992 c      print *,"Processor",myrank," opened file ISIDEP"
1993 c      print *,"Processor",myrank," opened parameter files"
1994 #elif (defined G77)
1995       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
1996       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
1997 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
1998 C Get parameter filenames and open the parameter files.
1999       call getenv_loc('BONDPAR',bondname)
2000       open (ibond,file=bondname,status='old')
2001       call getenv_loc('THETPAR',thetname)
2002       open (ithep,file=thetname,status='old')
2003       call getenv_loc('ROTPAR',rotname)
2004       open (irotam,file=rotname,status='old')
2005       call getenv_loc('TORPAR',torname)
2006       open (itorp,file=torname,status='old')
2007       call getenv_loc('TORDPAR',tordname)
2008       open (itordp,file=tordname,status='old')
2009       call getenv_loc('SCCORPAR',sccorname)
2010       open (isccor,file=sccorname,status='old')
2011       call getenv_loc('FOURIER',fouriername)
2012       open (ifourier,file=fouriername,status='old')
2013       call getenv_loc('ELEPAR',elename)
2014       open (ielep,file=elename,status='old')
2015       call getenv_loc('SIDEPAR',sidename)
2016       open (isidep,file=sidename,status='old')
2017 #else
2018       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2019      &  readonly)
2020        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2021 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2022 C Get parameter filenames and open the parameter files.
2023       call getenv_loc('BONDPAR',bondname)
2024       open (ibond,file=bondname,status='old',readonly)
2025       call getenv_loc('THETPAR',thetname)
2026       open (ithep,file=thetname,status='old',readonly)
2027       call getenv_loc('ROTPAR',rotname)
2028       open (irotam,file=rotname,status='old',readonly)
2029       call getenv_loc('TORPAR',torname)
2030       open (itorp,file=torname,status='old',readonly)
2031       call getenv_loc('TORDPAR',tordname)
2032       open (itordp,file=tordname,status='old',readonly)
2033       call getenv_loc('SCCORPAR',sccorname)
2034       open (isccor,file=sccorname,status='old',readonly)
2035 #ifndef CRYST_THETA
2036       call getenv_loc('THETPARPDB',thetname_pdb)
2037       print *,"thetname_pdb ",thetname_pdb
2038       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2039       print *,ithep_pdb," opened"
2040 #endif
2041       call getenv_loc('FOURIER',fouriername)
2042       open (ifourier,file=fouriername,status='old',readonly)
2043       call getenv_loc('ELEPAR',elename)
2044       open (ielep,file=elename,status='old',readonly)
2045       call getenv_loc('SIDEPAR',sidename)
2046       open (isidep,file=sidename,status='old',readonly)
2047 #ifndef CRYST_SC
2048       call getenv_loc('ROTPARPDB',rotname_pdb)
2049       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2050 #endif
2051 #endif
2052 #ifndef OLDSCP
2053 C
2054 C 8/9/01 In the newest version SCp interaction constants are read from a file
2055 C Use -DOLDSCP to use hard-coded constants instead.
2056 C
2057       call getenv_loc('SCPPAR',scpname)
2058 #if defined(WINIFL) || defined(WINPGI)
2059       open (iscpp,file=scpname,status='old',readonly,shared)
2060 #elif (defined CRAY)  || (defined AIX)
2061       open (iscpp,file=scpname,status='old',action='read')
2062 #elif (defined G77)
2063       open (iscpp,file=scpname,status='old')
2064 #else
2065       open (iscpp,file=scpname,status='old',readonly)
2066 #endif
2067 #endif
2068       call getenv_loc('PATTERN',patname)
2069 #if defined(WINIFL) || defined(WINPGI)
2070       open (icbase,file=patname,status='old',readonly,shared)
2071 #elif (defined CRAY)  || (defined AIX)
2072       open (icbase,file=patname,status='old',action='read')
2073 #elif (defined G77)
2074       open (icbase,file=patname,status='old')
2075 #else
2076       open (icbase,file=patname,status='old',readonly)
2077 #endif
2078 #ifdef MPI
2079 C Open output file only for CG processes
2080 c      print *,"Processor",myrank," fg_rank",fg_rank
2081       if (fg_rank.eq.0) then
2082
2083       if (nodes.eq.1) then
2084         npos=3
2085       else
2086         npos = dlog10(dfloat(nodes-1))+1
2087       endif
2088       if (npos.lt.3) npos=3
2089       write (liczba,'(i1)') npos
2090       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2091      &  //')'
2092       write (liczba,form) me
2093       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2094      &  liczba(:ilen(liczba))
2095       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2096      &  //'.int'
2097       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2098      &  //'.pdb'
2099       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2100      &  liczba(:ilen(liczba))//'.mol2'
2101       statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2102      &  liczba(:ilen(liczba))//'.stat'
2103       if (lentmp.gt.0)
2104      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2105      &      //liczba(:ilen(liczba))//'.stat')
2106       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2107      &  //'.rst'
2108       if(usampl) then
2109           qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2110      & liczba(:ilen(liczba))//'.const'
2111       endif
2112
2113       endif
2114 #else
2115       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2116       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2117       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2118       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2119       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2120       if (lentmp.gt.0)
2121      &  call copy_to_tmp(pref_orig(:ile(pref_orig))//'_'//pot(:lenpot)//
2122      &    //'.stat')
2123       rest2name=prefix(:ilen(prefix))//'.rst'
2124       if(usampl) then
2125          qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2126       endif
2127 #endif
2128 #if defined(AIX) || defined(PGI)
2129       if (me.eq.king .or. .not. out1file)
2130      &   open(iout,file=outname,status='unknown')
2131 #ifdef DEBUG
2132       if (fg_rank.gt.0) then
2133         write (liczba,'(i3.3)') myrank/nfgtasks
2134         write (ll,'(bz,i3.3)') fg_rank
2135         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2136      &   status='unknown')
2137       endif
2138 #endif
2139       if(me.eq.king) then
2140        open(igeom,file=intname,status='unknown',position='append')
2141        open(ipdb,file=pdbname,status='unknown')
2142        open(imol2,file=mol2name,status='unknown')
2143        open(istat,file=statname,status='unknown',position='append')
2144       else
2145 c1out       open(iout,file=outname,status='unknown')
2146       endif
2147 #else
2148       if (me.eq.king .or. .not.out1file)
2149      &    open(iout,file=outname,status='unknown')
2150 #ifdef DEBUG
2151       if (fg_rank.gt.0) then
2152         write (liczba,'(i3.3)') myrank/nfgtasks
2153         write (ll,'(bz,i3.3)') fg_rank
2154         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2155      &   status='unknown')
2156       endif
2157 #endif
2158       if(me.eq.king) then
2159        open(igeom,file=intname,status='unknown',access='append')
2160        open(ipdb,file=pdbname,status='unknown')
2161        open(imol2,file=mol2name,status='unknown')
2162        open(istat,file=statname,status='unknown',access='append')
2163       else
2164 c1out       open(iout,file=outname,status='unknown')
2165       endif
2166 #endif
2167       csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2168       csa_seed=prefix(:lenpre)//'.CSA.seed'
2169       csa_history=prefix(:lenpre)//'.CSA.history'
2170       csa_bank=prefix(:lenpre)//'.CSA.bank'
2171       csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2172       csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2173       csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2174 c!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2175       csa_int=prefix(:lenpre)//'.int'
2176       csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2177       csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2178       csa_in=prefix(:lenpre)//'.CSA.in'
2179 c      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2180 C Write file names
2181       if (me.eq.king)then
2182       write (iout,'(80(1h-))')
2183       write (iout,'(30x,a)') "FILE ASSIGNMENT"
2184       write (iout,'(80(1h-))')
2185       write (iout,*) "Input file                      : ",
2186      &  pref_orig(:ilen(pref_orig))//'.inp'
2187       write (iout,*) "Output file                     : ",
2188      &  outname(:ilen(outname))
2189       write (iout,*)
2190       write (iout,*) "Sidechain potential file        : ",
2191      &  sidename(:ilen(sidename))
2192 #ifndef OLDSCP
2193       write (iout,*) "SCp potential file              : ",
2194      &  scpname(:ilen(scpname))
2195 #endif
2196       write (iout,*) "Electrostatic potential file    : ",
2197      &  elename(:ilen(elename))
2198       write (iout,*) "Cumulant coefficient file       : ",
2199      &  fouriername(:ilen(fouriername))
2200       write (iout,*) "Torsional parameter file        : ",
2201      &  torname(:ilen(torname))
2202       write (iout,*) "Double torsional parameter file : ",
2203      &  tordname(:ilen(tordname))
2204       write (iout,*) "SCCOR parameter file : ",
2205      &  sccorname(:ilen(sccorname))
2206       write (iout,*) "Bond & inertia constant file    : ",
2207      &  bondname(:ilen(bondname))
2208       write (iout,*) "Bending parameter file          : ",
2209      &  thetname(:ilen(thetname))
2210       write (iout,*) "Rotamer parameter file          : ",
2211      &  rotname(:ilen(rotname))
2212       write (iout,*) "Thetpdb parameter file          : ",
2213      &  thetname_pdb(:ilen(thetname_pdb))
2214       write (iout,*) "Threading database              : ",
2215      &  patname(:ilen(patname))
2216       if (lentmp.ne.0)
2217      &write (iout,*)" DIRTMP                          : ",
2218      &  tmpdir(:lentmp)
2219       write (iout,'(80(1h-))')
2220       endif
2221       return
2222       end
2223 c----------------------------------------------------------------------------
2224       subroutine card_concat(card)
2225       implicit real*8 (a-h,o-z)
2226       include 'DIMENSIONS'
2227       include 'COMMON.IOUNITS'
2228       character*(*) card
2229       character*80 karta,ucase
2230       external ilen
2231       read (inp,'(a)') karta
2232       karta=ucase(karta)
2233       card=' '
2234       do while (karta(80:80).eq.'&')
2235         card=card(:ilen(card)+1)//karta(:79)
2236         read (inp,'(a)') karta
2237         karta=ucase(karta)
2238       enddo
2239       card=card(:ilen(card)+1)//karta
2240       return
2241       end
2242 c----------------------------------------------------------------------------------
2243       subroutine readrst
2244       implicit real*8 (a-h,o-z)
2245       include 'DIMENSIONS'
2246       include 'COMMON.CHAIN'
2247       include 'COMMON.IOUNITS'
2248       include 'COMMON.MD'
2249       open(irest2,file=rest2name,status='unknown')
2250       read(irest2,*) totT,EK,potE,totE,t_bath
2251       do i=1,2*nres
2252          read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2253       enddo
2254       do i=1,2*nres
2255          read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2256       enddo
2257       if(usampl) then
2258              read (irest2,*) iset
2259       endif
2260       close(irest2)
2261       return
2262       end
2263 c---------------------------------------------------------------------------------
2264       subroutine read_fragments
2265       implicit real*8 (a-h,o-z)
2266       include 'DIMENSIONS'
2267 #ifdef MPI
2268       include 'mpif.h'
2269 #endif
2270       include 'COMMON.SETUP'
2271       include 'COMMON.CHAIN'
2272       include 'COMMON.IOUNITS'
2273       include 'COMMON.MD'
2274       include 'COMMON.CONTROL'
2275       read(inp,*) nset,nfrag,npair,nfrag_back
2276       if(me.eq.king.or..not.out1file)
2277      & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2278      &  " nfrag_back",nfrag_back
2279       do iset=1,nset
2280          read(inp,*) mset(iset)
2281        do i=1,nfrag
2282          read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2283      &     qinfrag(i,iset)
2284          if(me.eq.king.or..not.out1file)
2285      &    write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2286      &     ifrag(2,i,iset), qinfrag(i,iset)
2287        enddo
2288        do i=1,npair
2289         read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2290      &    qinpair(i,iset)
2291         if(me.eq.king.or..not.out1file)
2292      &   write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2293      &    ipair(2,i,iset), qinpair(i,iset)
2294        enddo
2295        do i=1,nfrag_back
2296         read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2297      &     wfrag_back(3,i,iset),
2298      &     ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2299         if(me.eq.king.or..not.out1file)
2300      &   write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2301      &   wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2302        enddo
2303       enddo
2304       return
2305       end
2306 c-------------------------------------------------------------------------------
2307       subroutine read_dist_constr
2308       implicit real*8 (a-h,o-z)
2309       include 'DIMENSIONS'
2310 #ifdef MPI
2311       include 'mpif.h'
2312 #endif
2313       include 'COMMON.SETUP'
2314       include 'COMMON.CONTROL'
2315       include 'COMMON.CHAIN'
2316       include 'COMMON.IOUNITS'
2317       include 'COMMON.SBRIDGE'
2318       integer ifrag_(2,100),ipair_(2,100)
2319       double precision wfrag_(100),wpair_(100)
2320       character*500 controlcard
2321       print *, "WCHODZE"
2322       write (iout,*) "Calling read_dist_constr"
2323 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2324 c      call flush(iout)
2325       call card_concat(controlcard)
2326       call readi(controlcard,"NFRAG",nfrag_,0)
2327       call readi(controlcard,"NPAIR",npair_,0)
2328       call readi(controlcard,"NDIST",ndist_,0)
2329       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2330       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2331       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2332       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2333       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2334 c      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2335 c      write (iout,*) "IFRAG"
2336 c      do i=1,nfrag_
2337 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2338 c      enddo
2339 c      write (iout,*) "IPAIR"
2340 c      do i=1,npair_
2341 c        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2342 c      enddo
2343       call flush(iout)
2344       do i=1,nfrag_
2345         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2346         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2347      &    ifrag_(2,i)=nstart_sup+nsup-1
2348 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2349         call flush(iout)
2350         if (wfrag_(i).gt.0.0d0) then
2351         do j=ifrag_(1,i),ifrag_(2,i)-1
2352           do k=j+1,ifrag_(2,i)
2353 c            write (iout,*) "j",j," k",k
2354             ddjk=dist(j,k)
2355             if (constr_dist.eq.1) then
2356             nhpb=nhpb+1
2357             ihpb(nhpb)=j
2358             jhpb(nhpb)=k
2359               dhpb(nhpb)=ddjk
2360             forcon(nhpb)=wfrag_(i)
2361             else if (constr_dist.eq.2) then
2362               if (ddjk.le.dist_cut) then
2363                 nhpb=nhpb+1
2364                 ihpb(nhpb)=j
2365                 jhpb(nhpb)=k
2366                 dhpb(nhpb)=ddjk
2367                 forcon(nhpb)=wfrag_(i)
2368               endif
2369             else
2370               nhpb=nhpb+1
2371               ihpb(nhpb)=j
2372               jhpb(nhpb)=k
2373               dhpb(nhpb)=ddjk
2374               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2375             endif
2376 #ifdef MPI
2377             if (.not.out1file .or. me.eq.king)
2378      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2379      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2380 #else
2381             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2382      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2383 #endif
2384           enddo
2385         enddo
2386         endif
2387       enddo
2388       do i=1,npair_
2389         if (wpair_(i).gt.0.0d0) then
2390         ii = ipair_(1,i)
2391         jj = ipair_(2,i)
2392         if (ii.gt.jj) then
2393           itemp=ii
2394           ii=jj
2395           jj=itemp
2396         endif
2397         do j=ifrag_(1,ii),ifrag_(2,ii)
2398           do k=ifrag_(1,jj),ifrag_(2,jj)
2399             nhpb=nhpb+1
2400             ihpb(nhpb)=j
2401             jhpb(nhpb)=k
2402             forcon(nhpb)=wpair_(i)
2403             dhpb(nhpb)=dist(j,k)
2404 #ifdef MPI
2405             if (.not.out1file .or. me.eq.king)
2406      &      write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2407      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2408 #else
2409             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2410      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2411 #endif
2412           enddo
2413         enddo
2414         endif
2415       enddo
2416       print *,ndist_
2417       do i=1,ndist_
2418         if (constr_dist.eq.11) then
2419         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2420      &     ibecarb(i),forcon(nhpb+1),fordepth(nhpb+1)
2421         fordepth(nhpb+1)=fordepth(nhpb+1)/forcon(nhpb+1)
2422         else
2423 C        print *,"in else"
2424         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2425      &     ibecarb(i),forcon(nhpb+1)
2426         endif
2427         if (forcon(nhpb+1).gt.0.0d0) then
2428           nhpb=nhpb+1
2429           if (ibecarb(i).gt.0) then
2430             ihpb(i)=ihpb(i)+nres
2431             jhpb(i)=jhpb(i)+nres
2432           endif
2433           if (dhpb(nhpb).eq.0.0d0)
2434      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2435         endif
2436 C        read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),forcon(nhpb+1)
2437 C        if (forcon(nhpb+1).gt.0.0d0) then
2438 C          nhpb=nhpb+1
2439 C          dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2440 #ifdef MPI
2441           if (.not.out1file .or. me.eq.king)
2442      &    write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2443      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2444 #else
2445           write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2446      &     nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2447 #endif
2448
2449       enddo
2450       call flush(iout)
2451       return
2452       end
2453 c-------------------------------------------------------------------------------
2454 #ifdef WINIFL
2455       subroutine flush(iu)
2456       return
2457       end
2458 #endif
2459 #ifdef AIX
2460       subroutine flush(iu)
2461       call flush_(iu)
2462       return
2463       end
2464 #endif
2465 c------------------------------------------------------------------------------
2466       subroutine copy_to_tmp(source)
2467       include "DIMENSIONS"
2468       include "COMMON.IOUNITS"
2469       character*(*) source
2470       character* 256 tmpfile
2471       integer ilen
2472       external ilen
2473       logical ex
2474       tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
2475       inquire(file=tmpfile,exist=ex)
2476       if (ex) then
2477         write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
2478      &   " to temporary directory..."
2479         write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
2480         call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
2481       endif
2482       return
2483       end
2484 c------------------------------------------------------------------------------
2485       subroutine move_from_tmp(source)
2486       include "DIMENSIONS"
2487       include "COMMON.IOUNITS"
2488       character*(*) source
2489       integer ilen
2490       external ilen
2491       write (*,*) "Moving ",source(:ilen(source)),
2492      & " from temporary directory to working directory"
2493       write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
2494       call system("/bin/mv "//source(:ilen(source))//" "//curdir)
2495       return
2496       end
2497 c------------------------------------------------------------------------------
2498       subroutine random_init(seed)
2499 C
2500 C Initialize random number generator
2501 C
2502       implicit real*8 (a-h,o-z)
2503       include 'DIMENSIONS'
2504 #ifdef MPI
2505       include 'mpif.h'
2506       logical OKRandom, prng_restart
2507       real*8  r1
2508       integer iseed_array(4)
2509 #endif
2510       include 'COMMON.IOUNITS'
2511       include 'COMMON.TIME1'
2512       include 'COMMON.THREAD'
2513       include 'COMMON.SBRIDGE'
2514       include 'COMMON.CONTROL'
2515       include 'COMMON.MCM'
2516       include 'COMMON.MAP'
2517       include 'COMMON.HEADER'
2518       include 'COMMON.CSA'
2519       include 'COMMON.CHAIN'
2520       include 'COMMON.MUCA'
2521       include 'COMMON.MD'
2522       include 'COMMON.FFIELD'
2523       include 'COMMON.SETUP'
2524       iseed=-dint(dabs(seed))
2525       if (iseed.eq.0) then
2526         write (iout,'(/80(1h*)/20x,a/80(1h*))')
2527      &    'Random seed undefined. The program will stop.'
2528         write (*,'(/80(1h*)/20x,a/80(1h*))')
2529      &    'Random seed undefined. The program will stop.'
2530 #ifdef MPI
2531         call mpi_finalize(mpi_comm_world,ierr)
2532 #endif
2533         stop 'Bad random seed.'
2534       endif
2535 #ifdef MPI
2536       if (fg_rank.eq.0) then
2537       seed=seed*(me+1)+1
2538 #ifdef AMD64
2539       if(me.eq.king)
2540      &  write (iout,*) 'MPI: node= ', me, ' iseed= ',iseed
2541       OKRandom = prng_restart(me,iseed)
2542 #else
2543       do i=1,4
2544        tmp=65536.0d0**(4-i)
2545        iseed_array(i) = dint(seed/tmp)
2546        seed=seed-iseed_array(i)*tmp
2547       enddo
2548       if(me.eq.king)
2549      & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
2550      &                 (iseed_array(i),i=1,4)
2551       write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
2552      &                 (iseed_array(i),i=1,4)
2553       OKRandom = prng_restart(me,iseed_array)
2554 #endif
2555       if (OKRandom) then
2556 c        r1 = prng_next(me)
2557         r1=ran_number(0.0D0,1.0D0)
2558         if(me.eq.king)
2559      &   write (iout,*) 'ran_num',r1
2560         if (r1.lt.0.0d0) OKRandom=.false.
2561       endif
2562       if (.not.OKRandom) then
2563         write (iout,*) 'PRNG IS NOT WORKING!!!'
2564         print *,'PRNG IS NOT WORKING!!!'
2565         if (me.eq.0) then
2566          call flush(iout)
2567          call mpi_abort(mpi_comm_world,error_msg,ierr)
2568          stop
2569         else
2570          write (iout,*) 'too many processors for parallel prng'
2571          write (*,*) 'too many processors for parallel prng'
2572          call flush(iout)
2573          stop
2574         endif
2575       endif
2576       endif
2577 #else
2578       call vrndst(iseed)
2579       write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
2580 #endif
2581       return
2582       end