source/unres/src_MD/readrtns.F

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