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       else
1169         homol_nset=0
1170       endif
1171
1172
1173       if (nhpb.gt.0) call hpb_partition
1174 c      write (iout,*) "After read_dist_constr nhpb",nhpb
1175 c      call flush(iout)
1176       if (indpdb.eq.0 .and. modecalc.ne.2 .and. modecalc.ne.4
1177      &    .and. modecalc.ne.8 .and. modecalc.ne.9 .and.
1178      &    modecalc.ne.10) then
1179 C If input structure hasn't been supplied from the PDB file read or generate
1180 C initial geometry.
1181         if (iranconf.eq.0 .and. .not. extconf) then
1182           if(me.eq.king.or..not.out1file .and.fg_rank.eq.0)
1183      &     write (iout,'(a)') 'Initial geometry will be read in.'
1184           if (read_cart) then
1185             read(inp,'(8f10.5)',end=36,err=36)
1186      &       ((c(l,k),l=1,3),k=1,nres),
1187      &       ((c(l,k+nres),l=1,3),k=nnt,nct)
1188             call int_from_cart1(.false.)
1189             do i=1,nres-1
1190               do j=1,3
1191                 dc(j,i)=c(j,i+1)-c(j,i)
1192                 dc_norm(j,i)=dc_norm(j,i)*vbld_inv(i+1)
1193               enddo
1194             enddo
1195             do i=nnt,nct
1196               if (itype(i).ne.10) then
1197                 do j=1,3
1198                   dc(j,i+nres)=c(j,i+nres)-c(j,i)
1199                   dc_norm(j,i+nres)=dc_norm(j,i+nres)*vbld_inv(i+nres)
1200                 enddo
1201               endif
1202             enddo
1203             return
1204           else
1205             call read_angles(inp,*36)
1206           endif
1207           goto 37
1208    36     write (iout,'(a)') 'Error reading angle file.'
1209 #ifdef MPI
1210           call mpi_finalize( MPI_COMM_WORLD,IERR )
1211 #endif
1212           stop 'Error reading angle file.'
1213    37     continue
1214         else if (extconf) then
1215          if(me.eq.king.or..not.out1file .and. fg_rank.eq.0)
1216      &    write (iout,'(a)') 'Extended chain initial geometry.'
1217          do i=3,nres
1218           theta(i)=90d0*deg2rad
1219          enddo
1220          do i=4,nres
1221           phi(i)=180d0*deg2rad
1222          enddo
1223          do i=2,nres-1
1224           alph(i)=110d0*deg2rad
1225          enddo
1226          do i=2,nres-1
1227           omeg(i)=-120d0*deg2rad
1228          enddo
1229         else
1230           if(me.eq.king.or..not.out1file)
1231      &     write (iout,'(a)') 'Random-generated initial geometry.'
1232
1233
1234 #ifdef MPI
1235           if (me.eq.king  .or. fg_rank.eq.0 .and. (
1236      &           modecalc.eq.12 .or. modecalc.eq.14) ) then
1237 #endif
1238             do itrial=1,100
1239               itmp=1
1240               call gen_rand_conf(itmp,*30)
1241               goto 40
1242    30         write (iout,*) 'Failed to generate random conformation',
1243      &          ', itrial=',itrial
1244               write (*,*) 'Processor:',me,
1245      &          ' Failed to generate random conformation',
1246      &          ' itrial=',itrial
1247               call intout
1248
1249 #ifdef AIX
1250               call flush_(iout)
1251 #else
1252               call flush(iout)
1253 #endif
1254             enddo
1255             write (iout,'(a,i3,a)') 'Processor:',me,
1256      &        ' error in generating random conformation.'
1257             write (*,'(a,i3,a)') 'Processor:',me,
1258      &        ' error in generating random conformation.'
1259             call flush(iout)
1260 #ifdef MPI
1261             call MPI_Abort(mpi_comm_world,error_msg,ierrcode)
1262    40       continue
1263           endif
1264 #else
1265    40     continue
1266 #endif
1267         endif
1268       elseif (modecalc.eq.4) then
1269         read (inp,'(a)') intinname
1270         open (intin,file=intinname,status='old',err=333)
1271         if (me.eq.king .or. .not.out1file.and.fg_rank.eq.0)
1272      &  write (iout,'(a)') 'intinname',intinname
1273         write (*,'(a)') 'Processor',myrank,' intinname',intinname
1274         goto 334
1275   333   write (iout,'(2a)') 'Error opening angle file ',intinname
1276 #ifdef MPI
1277         call MPI_Finalize(MPI_COMM_WORLD,IERR)
1278 #endif
1279         stop 'Error opening angle file.'
1280   334   continue
1281
1282       endif
1283 C Generate distance constraints, if the PDB structure is to be regularized.
1284       if (nthread.gt.0) then
1285         call read_threadbase
1286       endif
1287       call setup_var
1288       if (me.eq.king .or. .not. out1file)
1289      & call intout
1290       if (ns.gt.0 .and. (me.eq.king .or. .not.out1file) ) then
1291         write (iout,'(/a,i3,a)')
1292      &  'The chain contains',ns,' disulfide-bridging cysteines.'
1293         write (iout,'(20i4)') (iss(i),i=1,ns)
1294        if (dyn_ss) then
1295           write(iout,*)"Running with dynamic disulfide-bond formation"
1296        else
1297         write (iout,'(/a/)') 'Pre-formed links are:'
1298         do i=1,nss
1299           i1=ihpb(i)-nres
1300           i2=jhpb(i)-nres
1301           it1=itype(i1)
1302           it2=itype(i2)
1303           write (iout,'(2a,i3,3a,i3,a,3f10.3)')
1304      &    restyp(it1),'(',i1,') -- ',restyp(it2),'(',i2,')',dhpb(i),
1305      &    ebr,forcon(i)
1306         enddo
1307         write (iout,'(a)')
1308        endif
1309       endif
1310       if (ns.gt.0.and.dyn_ss) then
1311           do i=nss+1,nhpb
1312             ihpb(i-nss)=ihpb(i)
1313             jhpb(i-nss)=jhpb(i)
1314             forcon(i-nss)=forcon(i)
1315             dhpb(i-nss)=dhpb(i)
1316           enddo
1317           nhpb=nhpb-nss
1318           nss=0
1319           call hpb_partition
1320           do i=1,ns
1321             dyn_ss_mask(iss(i))=.true.
1322           enddo
1323       endif
1324       if (i2ndstr.gt.0) call secstrp2dihc
1325 c      call geom_to_var(nvar,x)
1326 c      call etotal(energia(0))
1327 c      call enerprint(energia(0))
1328 c      call briefout(0,etot)
1329 c      stop
1330 cd    write (iout,'(2(a,i3))') 'NNT',NNT,' NCT',NCT
1331 cd    write (iout,'(a)') 'Variable list:'
1332 cd    write (iout,'(i4,f10.5)') (i,rad2deg*x(i),i=1,nvar)
1333 #ifdef MPI
1334       if (me.eq.king .or. (fg_rank.eq.0 .and. .not.out1file))
1335      &  write (iout,'(//80(1h*)/20x,a,i4,a/80(1h*)//)')
1336      &  'Processor',myrank,': end reading molecular data.'
1337 #endif
1338       return
1339       end
1340 c--------------------------------------------------------------------------
1341       logical function seq_comp(itypea,itypeb,length)
1342       implicit none
1343       integer length,itypea(length),itypeb(length)
1344       integer i
1345       do i=1,length
1346         if (itypea(i).ne.itypeb(i)) then
1347           seq_comp=.false.
1348           return
1349         endif
1350       enddo
1351       seq_comp=.true.
1352       return
1353       end
1354 c-----------------------------------------------------------------------------
1355       subroutine read_bridge
1356 C Read information about disulfide bridges.
1357       implicit real*8 (a-h,o-z)
1358       include 'DIMENSIONS'
1359 #ifdef MPI
1360       include 'mpif.h'
1361 #endif
1362       include 'COMMON.IOUNITS'
1363       include 'COMMON.GEO'
1364       include 'COMMON.VAR'
1365       include 'COMMON.INTERACT'
1366       include 'COMMON.LOCAL'
1367       include 'COMMON.NAMES'
1368       include 'COMMON.CHAIN'
1369       include 'COMMON.FFIELD'
1370       include 'COMMON.SBRIDGE'
1371       include 'COMMON.HEADER'
1372       include 'COMMON.CONTROL'
1373       include 'COMMON.DBASE'
1374       include 'COMMON.THREAD'
1375       include 'COMMON.TIME1'
1376       include 'COMMON.SETUP'
1377 C Read bridging residues.
1378       read (inp,*) ns,(iss(i),i=1,ns)
1379       print *,'ns=',ns
1380       if(me.eq.king.or..not.out1file)
1381      &  write (iout,*) 'ns=',ns,' iss:',(iss(i),i=1,ns)
1382 C Check whether the specified bridging residues are cystines.
1383       do i=1,ns
1384         if (itype(iss(i)).ne.1) then
1385           if (me.eq.king.or..not.out1file) write (iout,'(2a,i3,a)')
1386      &   'Do you REALLY think that the residue ',
1387      &    restyp(itype(iss(i))),i,
1388      &   ' can form a disulfide bridge?!!!'
1389           write (*,'(2a,i3,a)')
1390      &   'Do you REALLY think that the residue ',
1391      &    restyp(itype(iss(i))),i,
1392      &   ' can form a disulfide bridge?!!!'
1393 #ifdef MPI
1394          call MPI_Finalize(MPI_COMM_WORLD,ierror)
1395          stop
1396 #endif
1397         endif
1398       enddo
1399 C Read preformed bridges.
1400       if (ns.gt.0) then
1401       read (inp,*) nss,(ihpb(i),jhpb(i),i=1,nss)
1402       if(fg_rank.eq.0)
1403      & write(iout,*)'nss=',nss,' ihpb,jhpb: ',(ihpb(i),jhpb(i),i=1,nss)
1404       if (nss.gt.0) then
1405         nhpb=nss
1406 C Check if the residues involved in bridges are in the specified list of
1407 C bridging residues.
1408         do i=1,nss
1409           do j=1,i-1
1410             if (ihpb(i).eq.ihpb(j).or.ihpb(i).eq.jhpb(j)
1411      &      .or.jhpb(i).eq.ihpb(j).or.jhpb(i).eq.jhpb(j)) then
1412               write (iout,'(a,i3,a)') 'Disulfide pair',i,
1413      &      ' contains residues present in other pairs.'
1414               write (*,'(a,i3,a)') 'Disulfide pair',i,
1415      &      ' contains residues present in other pairs.'
1416 #ifdef MPI
1417               call MPI_Finalize(MPI_COMM_WORLD,ierror)
1418               stop
1419 #endif
1420             endif
1421           enddo
1422           do j=1,ns
1423             if (ihpb(i).eq.iss(j)) goto 10
1424           enddo
1425           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1426    10     continue
1427           do j=1,ns
1428             if (jhpb(i).eq.iss(j)) goto 20
1429           enddo
1430           write (iout,'(a,i3,a)') 'Pair',i,' contains unknown cystine.'
1431    20     continue
1432           dhpb(i)=dbr
1433           forcon(i)=fbr
1434         enddo
1435         do i=1,nss
1436           ihpb(i)=ihpb(i)+nres
1437           jhpb(i)=jhpb(i)+nres
1438         enddo
1439       endif
1440       endif
1441       return
1442       end
1443 c----------------------------------------------------------------------------
1444       subroutine read_x(kanal,*)
1445       implicit real*8 (a-h,o-z)
1446       include 'DIMENSIONS'
1447       include 'COMMON.GEO'
1448       include 'COMMON.VAR'
1449       include 'COMMON.CHAIN'
1450       include 'COMMON.IOUNITS'
1451       include 'COMMON.CONTROL'
1452       include 'COMMON.LOCAL'
1453       include 'COMMON.INTERACT'
1454 c Read coordinates from input
1455 c
1456       read(kanal,'(8f10.5)',end=10,err=10)
1457      &  ((c(l,k),l=1,3),k=1,nres),
1458      &  ((c(l,k+nres),l=1,3),k=nnt,nct)
1459       do j=1,3
1460         c(j,nres+1)=c(j,1)
1461         c(j,2*nres)=c(j,nres)
1462       enddo
1463       call int_from_cart1(.false.)
1464       do i=1,nres-1
1465         do j=1,3
1466           dc(j,i)=c(j,i+1)-c(j,i)
1467           dc_norm(j,i)=dc(j,i)*vbld_inv(i+1)
1468         enddo
1469       enddo
1470       do i=nnt,nct
1471         if (itype(i).ne.10) then
1472           do j=1,3
1473             dc(j,i+nres)=c(j,i+nres)-c(j,i)
1474             dc_norm(j,i+nres)=dc(j,i+nres)*vbld_inv(i+nres)
1475           enddo
1476         endif
1477       enddo
1478
1479       return
1480    10 return1
1481       end
1482 c----------------------------------------------------------------------------
1483       subroutine read_threadbase
1484       implicit real*8 (a-h,o-z)
1485       include 'DIMENSIONS'
1486       include 'COMMON.IOUNITS'
1487       include 'COMMON.GEO'
1488       include 'COMMON.VAR'
1489       include 'COMMON.INTERACT'
1490       include 'COMMON.LOCAL'
1491       include 'COMMON.NAMES'
1492       include 'COMMON.CHAIN'
1493       include 'COMMON.FFIELD'
1494       include 'COMMON.SBRIDGE'
1495       include 'COMMON.HEADER'
1496       include 'COMMON.CONTROL'
1497       include 'COMMON.DBASE'
1498       include 'COMMON.THREAD'
1499       include 'COMMON.TIME1'
1500 C Read pattern database for threading.
1501       read (icbase,*) nseq
1502       do i=1,nseq
1503         read (icbase,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1504      &   nres_base(2,i),nres_base(3,i)
1505         read (icbase,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1506      &   nres_base(1,i))
1507 c       write (iout,'(i5,2x,a8,2i4)') nres_base(1,i),str_nam(i),
1508 c    &   nres_base(2,i),nres_base(3,i)
1509 c       write (iout,'(9f8.3)') ((cart_base(k,j,i),k=1,3),j=1,
1510 c    &   nres_base(1,i))
1511       enddo
1512       close (icbase)
1513       if (weidis.eq.0.0D0) weidis=0.1D0
1514       do i=nnt,nct
1515         do j=i+2,nct
1516           nhpb=nhpb+1
1517           ihpb(nhpb)=i
1518           jhpb(nhpb)=j
1519           forcon(nhpb)=weidis
1520         enddo
1521       enddo
1522       read (inp,*) nexcl,(iexam(1,i),iexam(2,i),i=1,nexcl)
1523       write (iout,'(a,i5)') 'nexcl: ',nexcl
1524       write (iout,'(2i5)') (iexam(1,i),iexam(2,i),i=1,nexcl)
1525       return
1526       end
1527 c------------------------------------------------------------------------------
1528       subroutine setup_var
1529       implicit real*8 (a-h,o-z)
1530       include 'DIMENSIONS'
1531       include 'COMMON.IOUNITS'
1532       include 'COMMON.GEO'
1533       include 'COMMON.VAR'
1534       include 'COMMON.INTERACT'
1535       include 'COMMON.LOCAL'
1536       include 'COMMON.NAMES'
1537       include 'COMMON.CHAIN'
1538       include 'COMMON.FFIELD'
1539       include 'COMMON.SBRIDGE'
1540       include 'COMMON.HEADER'
1541       include 'COMMON.CONTROL'
1542       include 'COMMON.DBASE'
1543       include 'COMMON.THREAD'
1544       include 'COMMON.TIME1'
1545 C Set up variable list.
1546       ntheta=nres-2
1547       nphi=nres-3
1548       nvar=ntheta+nphi
1549       nside=0
1550       do i=2,nres-1
1551         if (itype(i).ne.10) then
1552           nside=nside+1
1553           ialph(i,1)=nvar+nside
1554           ialph(nside,2)=i
1555         endif
1556       enddo
1557       if (indphi.gt.0) then
1558         nvar=nphi
1559       else if (indback.gt.0) then
1560         nvar=nphi+ntheta
1561       else
1562         nvar=nvar+2*nside
1563       endif
1564 cd    write (iout,'(3i4)') (i,ialph(i,1),ialph(i,2),i=2,nres-1)
1565       return
1566       end
1567 c----------------------------------------------------------------------------
1568       subroutine gen_dist_constr
1569 C Generate CA distance constraints.
1570       implicit real*8 (a-h,o-z)
1571       include 'DIMENSIONS'
1572       include 'COMMON.IOUNITS'
1573       include 'COMMON.GEO'
1574       include 'COMMON.VAR'
1575       include 'COMMON.INTERACT'
1576       include 'COMMON.LOCAL'
1577       include 'COMMON.NAMES'
1578       include 'COMMON.CHAIN'
1579       include 'COMMON.FFIELD'
1580       include 'COMMON.SBRIDGE'
1581       include 'COMMON.HEADER'
1582       include 'COMMON.CONTROL'
1583       include 'COMMON.DBASE'
1584       include 'COMMON.THREAD'
1585       include 'COMMON.TIME1'
1586       dimension itype_pdb(maxres)
1587       common /pizda/ itype_pdb
1588       character*2 iden
1589 cd      print *,'gen_dist_constr: nnt=',nnt,' nct=',nct
1590 cd      write (2,*) 'gen_dist_constr: nnt=',nnt,' nct=',nct,
1591 cd     & ' nstart_sup',nstart_sup,' nstart_seq',nstart_seq,
1592 cd     & ' nsup',nsup
1593       do i=nstart_sup,nstart_sup+nsup-1
1594 cd      write (2,*) 'i',i,' seq ',restyp(itype(i+nstart_seq-nstart_sup)),
1595 cd     &    ' seq_pdb', restyp(itype_pdb(i))
1596         do j=i+2,nstart_sup+nsup-1
1597           nhpb=nhpb+1
1598           ihpb(nhpb)=i+nstart_seq-nstart_sup
1599           jhpb(nhpb)=j+nstart_seq-nstart_sup
1600           forcon(nhpb)=weidis
1601           dhpb(nhpb)=dist(i,j)
1602         enddo
1603       enddo
1604 cd      write (iout,'(a)') 'Distance constraints:'
1605 cd      do i=nss+1,nhpb
1606 cd        ii=ihpb(i)
1607 cd        jj=jhpb(i)
1608 cd        iden='CA'
1609 cd        if (ii.gt.nres) then
1610 cd          iden='SC'
1611 cd          ii=ii-nres
1612 cd          jj=jj-nres
1613 cd        endif
1614 cd        write (iout,'(a,1x,a,i4,3x,a,1x,a,i4,2f10.3)')
1615 cd     &  restyp(itype(ii)),iden,ii,restyp(itype(jj)),iden,jj,
1616 cd     &  dhpb(i),forcon(i)
1617 cd      enddo
1618       return
1619       end
1620 c----------------------------------------------------------------------------
1621       subroutine map_read
1622       implicit real*8 (a-h,o-z)
1623       include 'DIMENSIONS'
1624       include 'COMMON.MAP'
1625       include 'COMMON.IOUNITS'
1626       character*3 angid(4) /'THE','PHI','ALP','OME'/
1627       character*80 mapcard,ucase
1628       do imap=1,nmap
1629         read (inp,'(a)') mapcard
1630         mapcard=ucase(mapcard)
1631         if (index(mapcard,'PHI').gt.0) then
1632           kang(imap)=1
1633         else if (index(mapcard,'THE').gt.0) then
1634           kang(imap)=2
1635         else if (index(mapcard,'ALP').gt.0) then
1636           kang(imap)=3
1637         else if (index(mapcard,'OME').gt.0) then
1638           kang(imap)=4
1639         else
1640           write(iout,'(a)')'Error - illegal variable spec in MAP card.'
1641           stop 'Error - illegal variable spec in MAP card.'
1642         endif
1643         call readi (mapcard,'RES1',res1(imap),0)
1644         call readi (mapcard,'RES2',res2(imap),0)
1645         if (res1(imap).eq.0) then
1646           res1(imap)=res2(imap)
1647         else if (res2(imap).eq.0) then
1648           res2(imap)=res1(imap)
1649         endif
1650         if(res1(imap)*res2(imap).eq.0 .or. res1(imap).gt.res2(imap))then
1651           write (iout,'(a)')
1652      &    'Error - illegal definition of variable group in MAP.'
1653           stop 'Error - illegal definition of variable group in MAP.'
1654         endif
1655         call reada(mapcard,'FROM',ang_from(imap),0.0D0)
1656         call reada(mapcard,'TO',ang_to(imap),0.0D0)
1657         call readi(mapcard,'NSTEP',nstep(imap),0)
1658         if (ang_from(imap).eq.ang_to(imap) .or. nstep(imap).eq.0) then
1659           write (iout,'(a)')
1660      &     'Illegal boundary and/or step size specification in MAP.'
1661           stop 'Illegal boundary and/or step size specification in MAP.'
1662         endif
1663       enddo ! imap
1664       return
1665       end
1666 c----------------------------------------------------------------------------
1667 csa      subroutine csaread
1668 csa      implicit real*8 (a-h,o-z)
1669 csa      include 'DIMENSIONS'
1670 csa      include 'COMMON.IOUNITS'
1671 csa      include 'COMMON.GEO'
1672 csa      include 'COMMON.CSA'
1673 csa      include 'COMMON.BANK'
1674 csa      include 'COMMON.CONTROL'
1675 csa      character*80 ucase
1676 csa      character*620 mcmcard
1677 csa      call card_concat(mcmcard)
1678 csa
1679 csa      call readi(mcmcard,'NCONF',nconf,50)
1680 csa      call readi(mcmcard,'NADD',nadd,0)
1681 csa      call readi(mcmcard,'JSTART',jstart,1)
1682 csa      call readi(mcmcard,'JEND',jend,1)
1683 csa      call readi(mcmcard,'NSTMAX',nstmax,500000)
1684 csa      call readi(mcmcard,'N0',n0,1)
1685 csa      call readi(mcmcard,'N1',n1,6)
1686 csa      call readi(mcmcard,'N2',n2,4)
1687 csa      call readi(mcmcard,'N3',n3,0)
1688 csa      call readi(mcmcard,'N4',n4,0)
1689 csa      call readi(mcmcard,'N5',n5,0)
1690 csa      call readi(mcmcard,'N6',n6,10)
1691 csa      call readi(mcmcard,'N7',n7,0)
1692 csa      call readi(mcmcard,'N8',n8,0)
1693 csa      call readi(mcmcard,'N9',n9,0)
1694 csa      call readi(mcmcard,'N14',n14,0)
1695 csa      call readi(mcmcard,'N15',n15,0)
1696 csa      call readi(mcmcard,'N16',n16,0)
1697 csa      call readi(mcmcard,'N17',n17,0)
1698 csa      call readi(mcmcard,'N18',n18,0)
1699 csa
1700 csa      vdisulf=(index(mcmcard,'DYNSS').gt.0)
1701 csa
1702 csa      call readi(mcmcard,'NDIFF',ndiff,2)
1703 csa      call reada(mcmcard,'DIFFCUT',diffcut,0.0d0)
1704 csa      call readi(mcmcard,'IS1',is1,1)
1705 csa      call readi(mcmcard,'IS2',is2,8)
1706 csa      call readi(mcmcard,'NRAN0',nran0,4)
1707 csa      call readi(mcmcard,'NRAN1',nran1,2)
1708 csa      call readi(mcmcard,'IRR',irr,1)
1709 csa      call readi(mcmcard,'NSEED',nseed,20)
1710 csa      call readi(mcmcard,'NTOTAL',ntotal,10000)
1711 csa      call reada(mcmcard,'CUT1',cut1,2.0d0)
1712 csa      call reada(mcmcard,'CUT2',cut2,5.0d0)
1713 csa      call reada(mcmcard,'ESTOP',estop,-3000.0d0)
1714 csa      call readi(mcmcard,'ICMAX',icmax,3)
1715 csa      call readi(mcmcard,'IRESTART',irestart,0)
1716 csac!bankt      call readi(mcmcard,'NBANKTM',ntbankm,0)
1717 csa      ntbankm=0
1718 csac!bankt
1719 csa      call reada(mcmcard,'DELE',dele,20.0d0)
1720 csa      call reada(mcmcard,'DIFCUT',difcut,720.0d0)
1721 csa      call readi(mcmcard,'IREF',iref,0)
1722 csa      call reada(mcmcard,'RMSCUT',rmscut,4.0d0)
1723 csa      call reada(mcmcard,'PNCCUT',pnccut,0.5d0)
1724 csa      call readi(mcmcard,'NCONF_IN',nconf_in,0)
1725 csa      call reada(mcmcard,'RDIH_BIAS',rdih_bias,0.5d0)
1726 csa      write (iout,*) "NCONF_IN",nconf_in
1727 csa      return
1728 csa      end
1729 c----------------------------------------------------------------------------
1730 cfmc      subroutine mcmfread
1731 cfmc      implicit real*8 (a-h,o-z)
1732 cfmc      include 'DIMENSIONS'
1733 cfmc      include 'COMMON.MCMF'
1734 cfmc      include 'COMMON.IOUNITS'
1735 cfmc      include 'COMMON.GEO'
1736 cfmc      character*80 ucase
1737 cfmc      character*620 mcmcard
1738 cfmc      call card_concat(mcmcard)
1739 cfmc
1740 cfmc      call readi(mcmcard,'MAXRANT',maxrant,1000)
1741 cfmc      write(iout,*)'MAXRANT=',maxrant
1742 cfmc      call readi(mcmcard,'MAXFAM',maxfam,maxfam_p)
1743 cfmc      write(iout,*)'MAXFAM=',maxfam
1744 cfmc      call readi(mcmcard,'NNET1',nnet1,5)
1745 cfmc      write(iout,*)'NNET1=',nnet1
1746 cfmc      call readi(mcmcard,'NNET2',nnet2,4)
1747 cfmc      write(iout,*)'NNET2=',nnet2
1748 cfmc      call readi(mcmcard,'NNET3',nnet3,4)
1749 cfmc      write(iout,*)'NNET3=',nnet3
1750 cfmc      call readi(mcmcard,'ILASTT',ilastt,0)
1751 cfmc      write(iout,*)'ILASTT=',ilastt
1752 cfmc      call readi(mcmcard,'MAXSTR',maxstr,maxstr_mcmf)
1753 cfmc      write(iout,*)'MAXSTR=',maxstr
1754 cfmc      maxstr_f=maxstr/maxfam
1755 cfmc      write(iout,*)'MAXSTR_F=',maxstr_f
1756 cfmc      call readi(mcmcard,'NMCMF',nmcmf,10)
1757 cfmc      write(iout,*)'NMCMF=',nmcmf
1758 cfmc      call readi(mcmcard,'IFOCUS',ifocus,nmcmf)
1759 cfmc      write(iout,*)'IFOCUS=',ifocus
1760 cfmc      call readi(mcmcard,'NLOCMCMF',nlocmcmf,1000)
1761 cfmc      write(iout,*)'NLOCMCMF=',nlocmcmf
1762 cfmc      call readi(mcmcard,'INTPRT',intprt,1000)
1763 cfmc      write(iout,*)'INTPRT=',intprt
1764 cfmc      call readi(mcmcard,'IPRT',iprt,100)
1765 cfmc      write(iout,*)'IPRT=',iprt
1766 cfmc      call readi(mcmcard,'IMAXTR',imaxtr,100)
1767 cfmc      write(iout,*)'IMAXTR=',imaxtr
1768 cfmc      call readi(mcmcard,'MAXEVEN',maxeven,1000)
1769 cfmc      write(iout,*)'MAXEVEN=',maxeven
1770 cfmc      call readi(mcmcard,'MAXEVEN1',maxeven1,3)
1771 cfmc      write(iout,*)'MAXEVEN1=',maxeven1
1772 cfmc      call readi(mcmcard,'INIMIN',inimin,200)
1773 cfmc      write(iout,*)'INIMIN=',inimin
1774 cfmc      call readi(mcmcard,'NSTEPMCMF',nstepmcmf,10)
1775 cfmc      write(iout,*)'NSTEPMCMF=',nstepmcmf
1776 cfmc      call readi(mcmcard,'NTHREAD',nthread,5)
1777 cfmc      write(iout,*)'NTHREAD=',nthread
1778 cfmc      call readi(mcmcard,'MAXSTEPMCMF',maxstepmcmf,2500)
1779 cfmc      write(iout,*)'MAXSTEPMCMF=',maxstepmcmf
1780 cfmc      call readi(mcmcard,'MAXPERT',maxpert,9)
1781 cfmc      write(iout,*)'MAXPERT=',maxpert
1782 cfmc      call readi(mcmcard,'IRMSD',irmsd,1)
1783 cfmc      write(iout,*)'IRMSD=',irmsd
1784 cfmc      call reada(mcmcard,'DENEMIN',denemin,0.01D0)
1785 cfmc      write(iout,*)'DENEMIN=',denemin
1786 cfmc      call reada(mcmcard,'RCUT1S',rcut1s,3.5D0)
1787 cfmc      write(iout,*)'RCUT1S=',rcut1s
1788 cfmc      call reada(mcmcard,'RCUT1E',rcut1e,2.0D0)
1789 cfmc      write(iout,*)'RCUT1E=',rcut1e
1790 cfmc      call reada(mcmcard,'RCUT2S',rcut2s,0.5D0)
1791 cfmc      write(iout,*)'RCUT2S=',rcut2s
1792 cfmc      call reada(mcmcard,'RCUT2E',rcut2e,0.1D0)
1793 cfmc      write(iout,*)'RCUT2E=',rcut2e
1794 cfmc      call reada(mcmcard,'DPERT1',d_pert1,180.0D0)
1795 cfmc      write(iout,*)'DPERT1=',d_pert1
1796 cfmc      call reada(mcmcard,'DPERT1A',d_pert1a,180.0D0)
1797 cfmc      write(iout,*)'DPERT1A=',d_pert1a
1798 cfmc      call reada(mcmcard,'DPERT2',d_pert2,90.0D0)
1799 cfmc      write(iout,*)'DPERT2=',d_pert2
1800 cfmc      call reada(mcmcard,'DPERT2A',d_pert2a,45.0D0)
1801 cfmc      write(iout,*)'DPERT2A=',d_pert2a
1802 cfmc      call reada(mcmcard,'DPERT2B',d_pert2b,90.0D0)
1803 cfmc      write(iout,*)'DPERT2B=',d_pert2b
1804 cfmc      call reada(mcmcard,'DPERT2C',d_pert2c,60.0D0)
1805 cfmc      write(iout,*)'DPERT2C=',d_pert2c
1806 cfmc      d_pert1=deg2rad*d_pert1
1807 cfmc      d_pert1a=deg2rad*d_pert1a
1808 cfmc      d_pert2=deg2rad*d_pert2
1809 cfmc      d_pert2a=deg2rad*d_pert2a
1810 cfmc      d_pert2b=deg2rad*d_pert2b
1811 cfmc      d_pert2c=deg2rad*d_pert2c
1812 cfmc      call reada(mcmcard,'KT_MCMF1',kt_mcmf1,1.0D0)
1813 cfmc      write(iout,*)'KT_MCMF1=',kt_mcmf1
1814 cfmc      call reada(mcmcard,'KT_MCMF2',kt_mcmf2,1.0D0)
1815 cfmc      write(iout,*)'KT_MCMF2=',kt_mcmf2
1816 cfmc      call reada(mcmcard,'DKT_MCMF1',dkt_mcmf1,10.0D0)
1817 cfmc      write(iout,*)'DKT_MCMF1=',dkt_mcmf1
1818 cfmc      call reada(mcmcard,'DKT_MCMF2',dkt_mcmf2,1.0D0)
1819 cfmc      write(iout,*)'DKT_MCMF2=',dkt_mcmf2
1820 cfmc      call reada(mcmcard,'RCUTINI',rcutini,3.5D0)
1821 cfmc      write(iout,*)'RCUTINI=',rcutini
1822 cfmc      call reada(mcmcard,'GRAT',grat,0.5D0)
1823 cfmc      write(iout,*)'GRAT=',grat
1824 cfmc      call reada(mcmcard,'BIAS_MCMF',bias_mcmf,0.0D0)
1825 cfmc      write(iout,*)'BIAS_MCMF=',bias_mcmf
1826 cfmc
1827 cfmc      return
1828 cfmc      end
1829 c----------------------------------------------------------------------------
1830       subroutine mcmread
1831       implicit real*8 (a-h,o-z)
1832       include 'DIMENSIONS'
1833       include 'COMMON.MCM'
1834       include 'COMMON.MCE'
1835       include 'COMMON.IOUNITS'
1836       character*80 ucase
1837       character*320 mcmcard
1838       call card_concat(mcmcard)
1839       call readi(mcmcard,'MAXACC',maxacc,100)
1840       call readi(mcmcard,'MAX_MCM_IT',max_mcm_it,10000)
1841       call readi(mcmcard,'MAXTRIAL',maxtrial,100)
1842       call readi(mcmcard,'MAXTRIAL_ITER',maxtrial_iter,1000)
1843       call readi(mcmcard,'MAXREPM',maxrepm,200)
1844       call reada(mcmcard,'RANFRACT',RanFract,0.5D0)
1845       call reada(mcmcard,'POOL_FRACT',pool_fraction,0.01D0)
1846       call reada(mcmcard,'OVERLAP',overlap_cut,1.0D3)
1847       call reada(mcmcard,'E_UP',e_up,5.0D0)
1848       call reada(mcmcard,'DELTE',delte,0.1D0)
1849       call readi(mcmcard,'NSWEEP',nsweep,5)
1850       call readi(mcmcard,'NSTEPH',nsteph,0)
1851       call readi(mcmcard,'NSTEPC',nstepc,0)
1852       call reada(mcmcard,'TMIN',tmin,298.0D0)
1853       call reada(mcmcard,'TMAX',tmax,298.0D0)
1854       call readi(mcmcard,'NWINDOW',nwindow,0)
1855       call readi(mcmcard,'PRINT_MC',print_mc,0)
1856       print_stat=(index(mcmcard,'NO_PRINT_STAT').le.0)
1857       print_int=(index(mcmcard,'NO_PRINT_INT').le.0)
1858       ent_read=(index(mcmcard,'ENT_READ').gt.0)
1859       call readi(mcmcard,'SAVE_FREQ',save_frequency,1000)
1860       call readi(mcmcard,'MESSAGE_FREQ',message_frequency,1000)
1861       call readi(mcmcard,'POOL_READ_FREQ',pool_read_freq,5000)
1862       call readi(mcmcard,'POOL_SAVE_FREQ',pool_save_freq,1000)
1863       call readi(mcmcard,'PRINT_FREQ',print_freq,1000)
1864       if (nwindow.gt.0) then
1865         read (inp,*) (winstart(i),winend(i),i=1,nwindow)
1866         do i=1,nwindow
1867           winlen(i)=winend(i)-winstart(i)+1
1868         enddo
1869       endif
1870       if (tmax.lt.tmin) tmax=tmin
1871       if (tmax.eq.tmin) then
1872         nstepc=0
1873         nsteph=0
1874       endif
1875       if (nstepc.gt.0 .and. nsteph.gt.0) then
1876         tsteph=(tmax/tmin)**(1.0D0/(nsteph+0.0D0))
1877         tstepc=(tmax/tmin)**(1.0D0/(nstepc+0.0D0))
1878       endif
1879 C Probabilities of different move types
1880       sumpro_type(0)=0.0D0
1881       call reada(mcmcard,'MULTI_BOND',sumpro_type(1),1.0d0)
1882       call reada(mcmcard,'ONE_ANGLE' ,sumpro_type(2),2.0d0)
1883       sumpro_type(2)=sumpro_type(1)+sumpro_type(2)
1884       call reada(mcmcard,'THETA'     ,sumpro_type(3),0.0d0)
1885       sumpro_type(3)=sumpro_type(2)+sumpro_type(3)
1886       call reada(mcmcard,'SIDE_CHAIN',sumpro_type(4),0.5d0)
1887       sumpro_type(4)=sumpro_type(3)+sumpro_type(4)
1888       do i=1,MaxMoveType
1889         print *,'i',i,' sumprotype',sumpro_type(i)
1890         sumpro_type(i)=sumpro_type(i)/sumpro_type(MaxMoveType)
1891         print *,'i',i,' sumprotype',sumpro_type(i)
1892       enddo
1893       return
1894       end
1895 c----------------------------------------------------------------------------
1896       subroutine read_minim
1897       implicit real*8 (a-h,o-z)
1898       include 'DIMENSIONS'
1899       include 'COMMON.MINIM'
1900       include 'COMMON.IOUNITS'
1901       character*80 ucase
1902       character*320 minimcard
1903       call card_concat(minimcard)
1904       call readi(minimcard,'MAXMIN',maxmin,2000)
1905       call readi(minimcard,'MAXFUN',maxfun,5000)
1906       call readi(minimcard,'MINMIN',minmin,maxmin)
1907       call readi(minimcard,'MINFUN',minfun,maxmin)
1908       call reada(minimcard,'TOLF',tolf,1.0D-2)
1909       call reada(minimcard,'RTOLF',rtolf,1.0D-4)
1910       print_min_stat=min0(index(minimcard,'PRINT_MIN_STAT'),1)
1911       print_min_res=min0(index(minimcard,'PRINT_MIN_RES'),1)
1912       print_min_ini=min0(index(minimcard,'PRINT_MIN_INI'),1)
1913       write (iout,'(/80(1h*)/20x,a/80(1h*))')
1914      &         'Options in energy minimization:'
1915       write (iout,'(4(a,i5),a,1pe14.5,a,1pe14.5)')
1916      & 'MaxMin:',MaxMin,' MaxFun:',MaxFun,
1917      & 'MinMin:',MinMin,' MinFun:',MinFun,
1918      & ' TolF:',TolF,' RTolF:',RTolF
1919       return
1920       end
1921 c----------------------------------------------------------------------------
1922       subroutine read_angles(kanal,*)
1923       implicit real*8 (a-h,o-z)
1924       include 'DIMENSIONS'
1925       include 'COMMON.GEO'
1926       include 'COMMON.VAR'
1927       include 'COMMON.CHAIN'
1928       include 'COMMON.IOUNITS'
1929       include 'COMMON.CONTROL'
1930 c Read angles from input
1931 c
1932        read (kanal,*,err=10,end=10) (theta(i),i=3,nres)
1933        read (kanal,*,err=10,end=10) (phi(i),i=4,nres)
1934        read (kanal,*,err=10,end=10) (alph(i),i=2,nres-1)
1935        read (kanal,*,err=10,end=10) (omeg(i),i=2,nres-1)
1936
1937        do i=1,nres
1938 c 9/7/01 avoid 180 deg valence angle
1939         if (theta(i).gt.179.99d0) theta(i)=179.99d0
1940 c
1941         theta(i)=deg2rad*theta(i)
1942         phi(i)=deg2rad*phi(i)
1943         alph(i)=deg2rad*alph(i)
1944         omeg(i)=deg2rad*omeg(i)
1945        enddo
1946       return
1947    10 return1
1948       end
1949 c----------------------------------------------------------------------------
1950       subroutine reada(rekord,lancuch,wartosc,default)
1951       implicit none
1952       character*(*) rekord,lancuch
1953       double precision wartosc,default
1954       integer ilen,iread
1955       external ilen
1956       iread=index(rekord,lancuch)
1957       if (iread.eq.0) then
1958         wartosc=default
1959         return
1960       endif
1961       iread=iread+ilen(lancuch)+1
1962       read (rekord(iread:),*,err=10,end=10) wartosc
1963       return
1964   10  wartosc=default
1965       return
1966       end
1967 c----------------------------------------------------------------------------
1968       subroutine readi(rekord,lancuch,wartosc,default)
1969       implicit none
1970       character*(*) rekord,lancuch
1971       integer wartosc,default
1972       integer ilen,iread
1973       external ilen
1974       iread=index(rekord,lancuch)
1975       if (iread.eq.0) then
1976         wartosc=default
1977         return
1978       endif
1979       iread=iread+ilen(lancuch)+1
1980       read (rekord(iread:),*,err=10,end=10) wartosc
1981       return
1982   10  wartosc=default
1983       return
1984       end
1985 c----------------------------------------------------------------------------
1986       subroutine multreadi(rekord,lancuch,tablica,dim,default)
1987       implicit none
1988       integer dim,i
1989       integer tablica(dim),default
1990       character*(*) rekord,lancuch
1991       character*80 aux
1992       integer ilen,iread
1993       external ilen
1994       do i=1,dim
1995         tablica(i)=default
1996       enddo
1997       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
1998       if (iread.eq.0) return
1999       iread=iread+ilen(lancuch)+1
2000       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2001    10 return
2002       end
2003 c----------------------------------------------------------------------------
2004       subroutine multreada(rekord,lancuch,tablica,dim,default)
2005       implicit none
2006       integer dim,i
2007       double precision tablica(dim),default
2008       character*(*) rekord,lancuch
2009       character*80 aux
2010       integer ilen,iread
2011       external ilen
2012       do i=1,dim
2013         tablica(i)=default
2014       enddo
2015       iread=index(rekord,lancuch(:ilen(lancuch))//"=")
2016       if (iread.eq.0) return
2017       iread=iread+ilen(lancuch)+1
2018       read (rekord(iread:),*,end=10,err=10) (tablica(i),i=1,dim)
2019    10 return
2020       end
2021 c----------------------------------------------------------------------------
2022       subroutine openunits
2023       implicit real*8 (a-h,o-z)
2024       include 'DIMENSIONS'
2025 #ifdef MPI
2026       include 'mpif.h'
2027       character*16 form,nodename
2028       integer nodelen
2029 #endif
2030       include 'COMMON.SETUP'
2031       include 'COMMON.IOUNITS'
2032       include 'COMMON.MD'
2033       include 'COMMON.CONTROL'
2034       integer lenpre,lenpot,ilen,lentmp
2035       external ilen
2036       character*3 out1file_text,ucase
2037       character*3 ll
2038       external ucase
2039 c      print *,"Processor",myrank,"fg_rank",fg_rank," entered openunits"
2040       call getenv_loc("PREFIX",prefix)
2041       pref_orig = prefix
2042       call getenv_loc("POT",pot)
2043       call getenv_loc("DIRTMP",tmpdir)
2044       call getenv_loc("CURDIR",curdir)
2045       call getenv_loc("OUT1FILE",out1file_text)
2046 c      print *,"Processor",myrank,"fg_rank",fg_rank," did GETENV"
2047       out1file_text=ucase(out1file_text)
2048       if (out1file_text(1:1).eq."Y") then
2049         out1file=.true.
2050       else
2051         out1file=fg_rank.gt.0
2052       endif
2053       lenpre=ilen(prefix)
2054       lenpot=ilen(pot)
2055       lentmp=ilen(tmpdir)
2056       if (lentmp.gt.0) then
2057           write (*,'(80(1h!))')
2058           write (*,'(a,19x,a,19x,a)') "!","  A T T E N T I O N  ","!"
2059           write (*,'(80(1h!))')
2060           write (*,*)"All output files will be on node /tmp directory."
2061 #ifdef MPI
2062         call  MPI_GET_PROCESSOR_NAME( nodename, nodelen, IERROR )
2063         if (me.eq.king) then
2064           write (*,*) "The master node is ",nodename
2065         else if (fg_rank.eq.0) then
2066           write (*,*) "I am the CG slave node ",nodename
2067         else
2068           write (*,*) "I am the FG slave node ",nodename
2069         endif
2070 #endif
2071         PREFIX = tmpdir(:lentmp)//'/'//prefix(:lenpre)
2072         lenpre = lentmp+lenpre+1
2073       endif
2074       entname=prefix(:lenpre)//'_'//pot(:lenpot)//'.entr'
2075 C Get the names and open the input files
2076 #if defined(WINIFL) || defined(WINPGI)
2077       open(1,file=pref_orig(:ilen(pref_orig))//
2078      &  '.inp',status='old',readonly,shared)
2079        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2080 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2081 C Get parameter filenames and open the parameter files.
2082       call getenv_loc('BONDPAR',bondname)
2083       open (ibond,file=bondname,status='old',readonly,shared)
2084       call getenv_loc('THETPAR',thetname)
2085       open (ithep,file=thetname,status='old',readonly,shared)
2086 #ifndef CRYST_THETA
2087       call getenv_loc('THETPARPDB',thetname_pdb)
2088       open (ithep_pdb,file=thetname_pdb,status='old',readonly,shared)
2089 #endif
2090       call getenv_loc('ROTPAR',rotname)
2091       open (irotam,file=rotname,status='old',readonly,shared)
2092 #ifndef CRYST_SC
2093       call getenv_loc('ROTPARPDB',rotname_pdb)
2094       open (irotam_pdb,file=rotname_pdb,status='old',readonly,shared)
2095 #endif
2096       call getenv_loc('TORPAR',torname)
2097       open (itorp,file=torname,status='old',readonly,shared)
2098       call getenv_loc('TORDPAR',tordname)
2099       open (itordp,file=tordname,status='old',readonly,shared)
2100       call getenv_loc('FOURIER',fouriername)
2101       open (ifourier,file=fouriername,status='old',readonly,shared)
2102       call getenv_loc('ELEPAR',elename)
2103       open (ielep,file=elename,status='old',readonly,shared)
2104       call getenv_loc('SIDEPAR',sidename)
2105       open (isidep,file=sidename,status='old',readonly,shared)
2106 #elif (defined CRAY) || (defined AIX)
2107       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2108      &  action='read')
2109 c      print *,"Processor",myrank," opened file 1"
2110       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2111 c      print *,"Processor",myrank," opened file 9"
2112 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2113 C Get parameter filenames and open the parameter files.
2114       call getenv_loc('BONDPAR',bondname)
2115       open (ibond,file=bondname,status='old',action='read')
2116 c      print *,"Processor",myrank," opened file IBOND"
2117       call getenv_loc('THETPAR',thetname)
2118       open (ithep,file=thetname,status='old',action='read')
2119 c      print *,"Processor",myrank," opened file ITHEP"
2120 #ifndef CRYST_THETA
2121       call getenv_loc('THETPARPDB',thetname_pdb)
2122       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2123 #endif
2124       call getenv_loc('ROTPAR',rotname)
2125       open (irotam,file=rotname,status='old',action='read')
2126 c      print *,"Processor",myrank," opened file IROTAM"
2127 #ifndef CRYST_SC
2128       call getenv_loc('ROTPARPDB',rotname_pdb)
2129       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2130 #endif
2131       call getenv_loc('TORPAR',torname)
2132       open (itorp,file=torname,status='old',action='read')
2133 c      print *,"Processor",myrank," opened file ITORP"
2134       call getenv_loc('TORDPAR',tordname)
2135       open (itordp,file=tordname,status='old',action='read')
2136 c      print *,"Processor",myrank," opened file ITORDP"
2137       call getenv_loc('SCCORPAR',sccorname)
2138       open (isccor,file=sccorname,status='old',action='read')
2139 c      print *,"Processor",myrank," opened file ISCCOR"
2140       call getenv_loc('FOURIER',fouriername)
2141       open (ifourier,file=fouriername,status='old',action='read')
2142 c      print *,"Processor",myrank," opened file IFOURIER"
2143       call getenv_loc('ELEPAR',elename)
2144       open (ielep,file=elename,status='old',action='read')
2145 c      print *,"Processor",myrank," opened file IELEP"
2146       call getenv_loc('SIDEPAR',sidename)
2147       open (isidep,file=sidename,status='old',action='read')
2148 c      print *,"Processor",myrank," opened file ISIDEP"
2149 c      print *,"Processor",myrank," opened parameter files"
2150 #elif (defined G77)
2151       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old')
2152       open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2153 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2154 C Get parameter filenames and open the parameter files.
2155       call getenv_loc('BONDPAR',bondname)
2156       open (ibond,file=bondname,status='old')
2157       call getenv_loc('THETPAR',thetname)
2158       open (ithep,file=thetname,status='old')
2159 #ifndef CRYST_THETA
2160       call getenv_loc('THETPARPDB',thetname_pdb)
2161       open (ithep_pdb,file=thetname_pdb,status='old')
2162 #endif
2163       call getenv_loc('ROTPAR',rotname)
2164       open (irotam,file=rotname,status='old')
2165 #ifndef CRYST_SC
2166       call getenv_loc('ROTPARPDB',rotname_pdb)
2167       open (irotam_pdb,file=rotname_pdb,status='old')
2168 #endif
2169       call getenv_loc('TORPAR',torname)
2170       open (itorp,file=torname,status='old')
2171       call getenv_loc('TORDPAR',tordname)
2172       open (itordp,file=tordname,status='old')
2173       call getenv_loc('SCCORPAR',sccorname)
2174       open (isccor,file=sccorname,status='old')
2175       call getenv_loc('FOURIER',fouriername)
2176       open (ifourier,file=fouriername,status='old')
2177       call getenv_loc('ELEPAR',elename)
2178       open (ielep,file=elename,status='old')
2179       call getenv_loc('SIDEPAR',sidename)
2180       open (isidep,file=sidename,status='old')
2181 #else
2182       open(1,file=pref_orig(:ilen(pref_orig))//'.inp',status='old',
2183      &action='read')
2184        open (9,file=prefix(:ilen(prefix))//'.intin',status='unknown')
2185 C      open (18,file=prefix(:ilen(prefix))//'.entin',status='unknown')
2186 C Get parameter filenames and open the parameter files.
2187       call getenv_loc('BONDPAR',bondname)
2188       open (ibond,file=bondname,status='old',action='read')
2189       call getenv_loc('THETPAR',thetname)
2190       open (ithep,file=thetname,status='old',action='read')
2191 #ifndef CRYST_THETA
2192       call getenv_loc('THETPARPDB',thetname_pdb)
2193       print *,"thetname_pdb ",thetname_pdb
2194       open (ithep_pdb,file=thetname_pdb,status='old',action='read')
2195       print *,ithep_pdb," opened"
2196 #endif
2197       call getenv_loc('ROTPAR',rotname)
2198       open (irotam,file=rotname,status='old',action='read')
2199 #ifndef CRYST_SC
2200       call getenv_loc('ROTPARPDB',rotname_pdb)
2201       open (irotam_pdb,file=rotname_pdb,status='old',action='read')
2202 #endif
2203       call getenv_loc('TORPAR',torname)
2204       open (itorp,file=torname,status='old',action='read')
2205       call getenv_loc('TORDPAR',tordname)
2206       open (itordp,file=tordname,status='old',action='read')
2207       call getenv_loc('SCCORPAR',sccorname)
2208       open (isccor,file=sccorname,status='old',action='read')
2209       call getenv_loc('FOURIER',fouriername)
2210       open (ifourier,file=fouriername,status='old',action='read')
2211       call getenv_loc('ELEPAR',elename)
2212       open (ielep,file=elename,status='old',action='read')
2213       call getenv_loc('SIDEPAR',sidename)
2214       open (isidep,file=sidename,status='old',action='read')
2215 #endif
2216 #ifndef OLDSCP
2217 C
2218 C 8/9/01 In the newest version SCp interaction constants are read from a file
2219 C Use -DOLDSCP to use hard-coded constants instead.
2220 C
2221       call getenv_loc('SCPPAR',scpname)
2222 #if defined(WINIFL) || defined(WINPGI)
2223       open (iscpp,file=scpname,status='old',readonly,shared)
2224 #elif (defined CRAY)  || (defined AIX)
2225       open (iscpp,file=scpname,status='old',action='read')
2226 #elif (defined G77)
2227       open (iscpp,file=scpname,status='old')
2228 #else
2229       open (iscpp,file=scpname,status='old',action='read')
2230 #endif
2231 #endif
2232       call getenv_loc('PATTERN',patname)
2233 #if defined(WINIFL) || defined(WINPGI)
2234       open (icbase,file=patname,status='old',readonly,shared)
2235 #elif (defined CRAY)  || (defined AIX)
2236       open (icbase,file=patname,status='old',action='read')
2237 #elif (defined G77)
2238       open (icbase,file=patname,status='old')
2239 #else
2240       open (icbase,file=patname,status='old',action='read')
2241 #endif
2242 #ifdef MPI
2243 C Open output file only for CG processes
2244 c      print *,"Processor",myrank," fg_rank",fg_rank
2245       if (fg_rank.eq.0) then
2246
2247       if (nodes.eq.1) then
2248         npos=3
2249       else
2250         npos = dlog10(dfloat(nodes-1))+1
2251       endif
2252       if (npos.lt.3) npos=3
2253       write (liczba,'(i1)') npos
2254       form = '(bz,i'//liczba(:ilen(liczba))//'.'//liczba(:ilen(liczba))
2255      &  //')'
2256       write (liczba,form) me
2257       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)//
2258      &  liczba(:ilen(liczba))
2259       intname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2260      &  //'.int'
2261       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//liczba(:ilen(liczba))
2262      &  //'.pdb'
2263       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//
2264      &  liczba(:ilen(liczba))//'.mol2'
2265       statname=prefix(:lenpre)//'_'//pot(:lenpot)//
2266      &  liczba(:ilen(liczba))//'.stat'
2267       if (lentmp.gt.0)
2268      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2269      &      //liczba(:ilen(liczba))//'.stat')
2270       rest2name=prefix(:ilen(prefix))//"_"//liczba(:ilen(liczba))
2271      &  //'.rst'
2272       if(usampl) then
2273           qname=prefix(:lenpre)//'_'//pot(:lenpot)//
2274      & liczba(:ilen(liczba))//'.const'
2275       endif
2276
2277       endif
2278 #else
2279       outname=prefix(:lenpre)//'.out_'//pot(:lenpot)
2280       intname=prefix(:lenpre)//'_'//pot(:lenpot)//'.int'
2281       pdbname=prefix(:lenpre)//'_'//pot(:lenpot)//'.pdb'
2282       mol2name=prefix(:lenpre)//'_'//pot(:lenpot)//'.mol2'
2283       statname=prefix(:lenpre)//'_'//pot(:lenpot)//'.stat'
2284       if (lentmp.gt.0)
2285      &  call copy_to_tmp(pref_orig(:ilen(pref_orig))//'_'//pot(:lenpot)
2286      &    //'.stat')
2287       rest2name=prefix(:ilen(prefix))//'.rst'
2288       if(usampl) then
2289          qname=prefix(:lenpre)//'_'//pot(:lenpot)//'.const'
2290       endif
2291 #endif
2292 #if defined(AIX) || defined(PGI)
2293       if (me.eq.king .or. .not. out1file)
2294      &   open(iout,file=outname,status='unknown')
2295 c#define DEBUG
2296 #ifdef DEBUG
2297       if (fg_rank.gt.0) then
2298         write (liczba,'(i3.3)') myrank/nfgtasks
2299         write (ll,'(bz,i3.3)') fg_rank
2300         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2301      &   status='unknown')
2302       endif
2303 #endif
2304 c#undef DEBUG
2305       if(me.eq.king) then
2306        open(igeom,file=intname,status='unknown',position='append')
2307        open(ipdb,file=pdbname,status='unknown')
2308        open(imol2,file=mol2name,status='unknown')
2309        open(istat,file=statname,status='unknown',position='append')
2310       else
2311 c1out       open(iout,file=outname,status='unknown')
2312       endif
2313 #else
2314       if (me.eq.king .or. .not.out1file)
2315      &    open(iout,file=outname,status='unknown')
2316 c#define DEBUG
2317 #ifdef DEBUG
2318       if (fg_rank.gt.0) then
2319         print "Processor",fg_rank," opening output file"
2320         write (liczba,'(i3.3)') myrank/nfgtasks
2321         write (ll,'(bz,i3.3)') fg_rank
2322         open(iout,file="debug"//liczba(:ilen(liczba))//"."//ll,
2323      &   status='unknown')
2324       endif
2325 #endif
2326 c#undef DEBUG
2327       if(me.eq.king) then
2328        open(igeom,file=intname,status='unknown',access='append')
2329        open(ipdb,file=pdbname,status='unknown')
2330        open(imol2,file=mol2name,status='unknown')
2331        open(istat,file=statname,status='unknown',access='append')
2332       else
2333 c1out       open(iout,file=outname,status='unknown')
2334       endif
2335 #endif
2336 csa      csa_rbank=prefix(:lenpre)//'.CSA.rbank'
2337 csa      csa_seed=prefix(:lenpre)//'.CSA.seed'
2338 csa      csa_history=prefix(:lenpre)//'.CSA.history'
2339 csa      csa_bank=prefix(:lenpre)//'.CSA.bank'
2340 csa      csa_bank1=prefix(:lenpre)//'.CSA.bank1'
2341 csa      csa_alpha=prefix(:lenpre)//'.CSA.alpha'
2342 csa      csa_alpha1=prefix(:lenpre)//'.CSA.alpha1'
2343 csac!bankt      csa_bankt=prefix(:lenpre)//'.CSA.bankt'
2344 csa      csa_int=prefix(:lenpre)//'.int'
2345 csa      csa_bank_reminimized=prefix(:lenpre)//'.CSA.bank_reminimized'
2346 csa      csa_native_int=prefix(:lenpre)//'.CSA.native.int'
2347 csa      csa_in=prefix(:lenpre)//'.CSA.in'
2348 c      print *,"Processor",myrank,"fg_rank",fg_rank," opened files"
2349 C Write file names
2350       if (me.eq.king)then
2351       write (iout,'(80(1h-))')
2352       write (iout,'(30x,a)') "FILE ASSIGNMENT"
2353       write (iout,'(80(1h-))')
2354       write (iout,*) "Input file                      : ",
2355      &  pref_orig(:ilen(pref_orig))//'.inp'
2356       write (iout,*) "Output file                     : ",
2357      &  outname(:ilen(outname))
2358       write (iout,*)
2359       write (iout,*) "Sidechain potential file        : ",
2360      &  sidename(:ilen(sidename))
2361 #ifndef OLDSCP
2362       write (iout,*) "SCp potential file              : ",
2363      &  scpname(:ilen(scpname))
2364 #endif
2365       write (iout,*) "Electrostatic potential file    : ",
2366      &  elename(:ilen(elename))
2367       write (iout,*) "Cumulant coefficient file       : ",
2368      &  fouriername(:ilen(fouriername))
2369       write (iout,*) "Torsional parameter file        : ",
2370      &  torname(:ilen(torname))
2371       write (iout,*) "Double torsional parameter file : ",
2372      &  tordname(:ilen(tordname))
2373       write (iout,*) "SCCOR parameter file : ",
2374      &  sccorname(:ilen(sccorname))
2375       write (iout,*) "Bond & inertia constant file    : ",
2376      &  bondname(:ilen(bondname))
2377       write (iout,*) "Bending parameter file          : ",
2378      &  thetname(:ilen(thetname))
2379       write (iout,*) "Rotamer parameter file          : ",
2380      &  rotname(:ilen(rotname))
2381       write (iout,*) "Threading database              : ",
2382      &  patname(:ilen(patname))
2383       if (lentmp.ne.0)
2384      &write (iout,*)" DIRTMP                          : ",
2385      &  tmpdir(:lentmp)
2386       write (iout,'(80(1h-))')
2387       endif
2388       return
2389       end
2390 c----------------------------------------------------------------------------
2391       subroutine card_concat(card)
2392       implicit real*8 (a-h,o-z)
2393       include 'DIMENSIONS'
2394       include 'COMMON.IOUNITS'
2395       character*(*) card
2396       character*80 karta,ucase
2397       external ilen
2398       read (inp,'(a)') karta
2399       karta=ucase(karta)
2400       card=' '
2401       do while (karta(80:80).eq.'&')
2402         card=card(:ilen(card)+1)//karta(:79)
2403         read (inp,'(a)') karta
2404         karta=ucase(karta)
2405       enddo
2406       card=card(:ilen(card)+1)//karta
2407       return
2408       end
2409 c----------------------------------------------------------------------------------
2410       subroutine readrst
2411       implicit real*8 (a-h,o-z)
2412       include 'DIMENSIONS'
2413       include 'COMMON.CHAIN'
2414       include 'COMMON.IOUNITS'
2415       include 'COMMON.MD'
2416       include 'COMMON.CONTROL'
2417       open(irest2,file=rest2name,status='unknown')
2418       read(irest2,*) totT,EK,potE,totE,t_bath
2419       do i=1,2*nres
2420          read(irest2,'(3e15.5)') (d_t(j,i),j=1,3)
2421       enddo
2422       do i=1,2*nres
2423          read(irest2,'(3e15.5)') (dc(j,i),j=1,3)
2424       enddo
2425       if(usampl.or.homol_nset.gt.1) then
2426              read (irest2,*) iset
2427       endif
2428       close(irest2)
2429       return
2430       end
2431 c---------------------------------------------------------------------------------
2432       subroutine read_fragments
2433       implicit real*8 (a-h,o-z)
2434       include 'DIMENSIONS'
2435 #ifdef MPI
2436       include 'mpif.h'
2437 #endif
2438       include 'COMMON.SETUP'
2439       include 'COMMON.CHAIN'
2440       include 'COMMON.IOUNITS'
2441       include 'COMMON.MD'
2442       include 'COMMON.CONTROL'
2443       read(inp,*) nset,nfrag,npair,nfrag_back
2444       if(me.eq.king.or..not.out1file)
2445      & write(iout,*) "nset",nset," nfrag",nfrag," npair",npair,
2446      &  " nfrag_back",nfrag_back
2447       do iset=1,nset
2448          read(inp,*) mset(iset)
2449        do i=1,nfrag
2450          read(inp,*) wfrag(i,iset),ifrag(1,i,iset),ifrag(2,i,iset),
2451      &     qinfrag(i,iset)
2452          if(me.eq.king.or..not.out1file)
2453      &    write(iout,*) "R ",i,wfrag(i,iset),ifrag(1,i,iset),
2454      &     ifrag(2,i,iset), qinfrag(i,iset)
2455        enddo
2456        do i=1,npair
2457         read(inp,*) wpair(i,iset),ipair(1,i,iset),ipair(2,i,iset),
2458      &    qinpair(i,iset)
2459         if(me.eq.king.or..not.out1file)
2460      &   write(iout,*) "R ",i,wpair(i,iset),ipair(1,i,iset),
2461      &    ipair(2,i,iset), qinpair(i,iset)
2462        enddo
2463        do i=1,nfrag_back
2464         read(inp,*) wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2465      &     wfrag_back(3,i,iset),
2466      &     ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2467         if(me.eq.king.or..not.out1file)
2468      &   write(iout,*) "A",i,wfrag_back(1,i,iset),wfrag_back(2,i,iset),
2469      &   wfrag_back(3,i,iset),ifrag_back(1,i,iset),ifrag_back(2,i,iset)
2470        enddo
2471       enddo
2472       return
2473       end
2474 c-------------------------------------------------------------------------------
2475       subroutine read_dist_constr
2476       implicit real*8 (a-h,o-z)
2477       include 'DIMENSIONS'
2478 #ifdef MPI
2479       include 'mpif.h'
2480 #endif
2481       include 'COMMON.SETUP'
2482       include 'COMMON.CONTROL'
2483       include 'COMMON.CHAIN'
2484       include 'COMMON.IOUNITS'
2485       include 'COMMON.SBRIDGE'
2486       integer ifrag_(2,100),ipair_(2,100)
2487       double precision wfrag_(100),wpair_(100)
2488       character*500 controlcard
2489 c      write (iout,*) "Calling read_dist_constr"
2490 c      write (iout,*) "nres",nres," nstart_sup",nstart_sup," nsup",nsup
2491 c      call flush(iout)
2492       call card_concat(controlcard)
2493       call readi(controlcard,"NFRAG",nfrag_,0)
2494       call readi(controlcard,"NPAIR",npair_,0)
2495       call readi(controlcard,"NDIST",ndist_,0)
2496       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0)
2497       call multreadi(controlcard,"IFRAG",ifrag_(1,1),2*nfrag_,0)
2498       call multreadi(controlcard,"IPAIR",ipair_(1,1),2*npair_,0)
2499       call multreada(controlcard,"WFRAG",wfrag_(1),nfrag_,0.0d0)
2500       call multreada(controlcard,"WPAIR",wpair_(1),npair_,0.0d0)
2501 c      write (iout,*) "NFRAG",nfrag_," NPAIR",npair_," NDIST",ndist_
2502 c      write (iout,*) "IFRAG"
2503 c      do i=1,nfrag_
2504 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2505 c      enddo
2506 c      write (iout,*) "IPAIR"
2507 c      do i=1,npair_
2508 c        write (iout,*) i,ipair_(1,i),ipair_(2,i),wpair_(i)
2509 c      enddo
2510       if (.not.refstr .and. nfrag.gt.0) then
2511         write (iout,*)
2512      &  "ERROR: no reference structure to compute distance restraints"
2513         write (iout,*)
2514      &  "Restraints must be specified explicitly (NDIST=number)"
2515         stop
2516       endif
2517       if (nfrag.lt.2 .and. npair.gt.0) then
2518         write (iout,*) "ERROR: Less than 2 fragments specified",
2519      &   " but distance restraints between pairs requested"
2520         stop
2521       endif
2522       call flush(iout)
2523       do i=1,nfrag_
2524         if (ifrag_(1,i).lt.nstart_sup) ifrag_(1,i)=nstart_sup
2525         if (ifrag_(2,i).gt.nstart_sup+nsup-1)
2526      &    ifrag_(2,i)=nstart_sup+nsup-1
2527 c        write (iout,*) i,ifrag_(1,i),ifrag_(2,i),wfrag_(i)
2528         call flush(iout)
2529         if (wfrag_(i).gt.0.0d0) then
2530         do j=ifrag_(1,i),ifrag_(2,i)-1
2531           do k=j+1,ifrag_(2,i)
2532 c            write (iout,*) "j",j," k",k
2533             ddjk=dist(j,k)
2534             if (constr_dist.eq.1) then
2535             nhpb=nhpb+1
2536             ihpb(nhpb)=j
2537             jhpb(nhpb)=k
2538               dhpb(nhpb)=ddjk
2539             forcon(nhpb)=wfrag_(i)
2540             else if (constr_dist.eq.2) then
2541               if (ddjk.le.dist_cut) then
2542                 nhpb=nhpb+1
2543                 ihpb(nhpb)=j
2544                 jhpb(nhpb)=k
2545                 dhpb(nhpb)=ddjk
2546                 forcon(nhpb)=wfrag_(i)
2547               endif
2548             else
2549               nhpb=nhpb+1
2550               ihpb(nhpb)=j
2551               jhpb(nhpb)=k
2552               dhpb(nhpb)=ddjk
2553               forcon(nhpb)=wfrag_(i)*dexp(-0.5d0*(ddjk/dist_cut)**2)
2554             endif
2555 #ifdef MPI
2556             if (.not.out1file .or. me.eq.king)
2557      &      write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2558      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2559 #else
2560             write (iout,'(a,3i5,f8.2,1pe12.2)') "+dist.constr ",
2561      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2562 #endif
2563           enddo
2564         enddo
2565         endif
2566       enddo
2567       do i=1,npair_
2568         if (wpair_(i).gt.0.0d0) then
2569         ii = ipair_(1,i)
2570         jj = ipair_(2,i)
2571         if (ii.gt.jj) then
2572           itemp=ii
2573           ii=jj
2574           jj=itemp
2575         endif
2576         do j=ifrag_(1,ii),ifrag_(2,ii)
2577           do k=ifrag_(1,jj),ifrag_(2,jj)
2578             nhpb=nhpb+1
2579             ihpb(nhpb)=j
2580             jhpb(nhpb)=k
2581             forcon(nhpb)=wpair_(i)
2582             dhpb(nhpb)=dist(j,k)
2583 #ifdef MPI
2584             if (.not.out1file .or. me.eq.king)
2585      &      write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2586      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2587 #else
2588             write (iout,'(a,3i5,f8.2,f10.1)') "+dist.constr ",
2589      &       nhpb,ihpb(nhpb),jhpb(nhpb),dhpb(nhpb),forcon(nhpb)
2590 #endif
2591           enddo
2592         enddo
2593         endif
2594       enddo
2595       do i=1,ndist_
2596         read (inp,*) ihpb(nhpb+1),jhpb(nhpb+1),dhpb(i),dhpb1(i),
2597      &     ibecarb(i),forcon(nhpb+1)
2598         if (forcon(nhpb+1).gt.0.0d0) then
2599           nhpb=nhpb+1
2600           if (ibecarb(i).gt.0) then
2601             ihpb(i)=ihpb(i)+nres
2602             jhpb(i)=jhpb(i)+nres
2603           endif
2604           if (dhpb(nhpb).eq.0.0d0)
2605      &       dhpb(nhpb)=dist(ihpb(nhpb),jhpb(nhpb))
2606         endif
2607       enddo
2608 #ifdef MPI
2609       if (.not.out1file .or. me.eq.king) then
2610 #endif
2611       do i=1,nhpb
2612           write (iout,'(a,3i5,2f8.2,i2,f10.1)') "+dist.constr ",
2613      &     i,ihpb(i),jhpb(i),dhpb(i),dhpb1(i),ibecarb(i),forcon(i)
2614       enddo
2615       call flush(iout)
2616 #ifdef MPI
2617       endif
2618 #endif
2619       return
2620       end
2621 c-------------------------------------------------------------------------------
2622
2623       subroutine read_constr_homology
2624
2625       include 'DIMENSIONS'
2626 #ifdef MPI
2627       include 'mpif.h'
2628 #endif
2629       include 'COMMON.SETUP'
2630       include 'COMMON.CONTROL'
2631       include 'COMMON.CHAIN'
2632       include 'COMMON.IOUNITS'
2633       include 'COMMON.MD'
2634       include 'COMMON.GEO'
2635       include 'COMMON.INTERACT'
2636 c
2637 c For new homol impl
2638 c
2639       include 'COMMON.VAR'
2640 c
2641
2642 c     double precision odl_temp,sigma_odl_temp,waga_theta,waga_d,
2643 c    &                 dist_cut
2644 c     common /przechowalnia/ odl_temp(maxres,maxres,max_template),
2645 c    &    sigma_odl_temp(maxres,maxres,max_template)
2646       character*2 kic2
2647       character*24 model_ki_dist, model_ki_angle
2648       character*500 controlcard
2649       integer ki, i, j, k, l
2650       logical lprn /.true./
2651 c
2652 c     FP - Nov. 2014 Temporary specifications for new vars
2653 c
2654       double precision rescore_tmp,x12,y12,z12
2655       double precision, dimension (max_template,maxres) :: rescore
2656       character*24 pdbfile,tpl_k_rescore
2657 c -----------------------------------------------------------------
2658 c Reading multiple PDB ref structures and calculation of retraints
2659 c not using pre-computed ones stored in files model_ki_{dist,angle}
2660 c FP (Nov., 2014)
2661 c -----------------------------------------------------------------
2662 c
2663 c
2664 c Alternative: reading from input
2665       call card_concat(controlcard)
2666       call reada(controlcard,"HOMOL_DIST",waga_dist,1.0d0)
2667       call reada(controlcard,"HOMOL_ANGLE",waga_angle,1.0d0)
2668       call reada(controlcard,"HOMOL_THETA",waga_theta,1.0d0) ! new
2669       call reada(controlcard,"HOMOL_SCD",waga_d,1.0d0) ! new
2670       call reada(controlcard,'DIST_CUT',dist_cut,5.0d0) ! for diff ways of calc sigma
2671
2672       call readi(controlcard,"HOMOL_NSET",homol_nset,1)
2673       if (homol_nset.gt.1)then
2674          call card_concat(controlcard)
2675          read(controlcard,*) (waga_dist1(i),i=1,homol_nset)
2676          call card_concat(controlcard)
2677          read(controlcard,*) (waga_angle1(i),i=1,homol_nset)
2678          write(iout,*) "iset distance_weight angle_weight"
2679          do i=1,homol_nset
2680            write(iout,*) i,waga_dist1(i),waga_angle1(i)
2681          enddo
2682          iset=mod(kolor,homol_nset)+1
2683       else
2684        iset=1
2685 c       call reada(controlcard,"HOMOL_DIST",waga_dist(1),1.0d0)
2686 c       call reada(controlcard,"HOMOL_ANGLE",waga_angle(1),1.0d0)
2687       endif
2688
2689       write (iout,*) "nnt",nnt," nct",nct
2690       call flush(iout)
2691
2692
2693       lim_odl=0
2694       lim_dih=0
2695 c
2696 c  New
2697 c
2698       lim_theta=0
2699       lim_xx=0
2700 c
2701 c  Reading HM global scores (prob not required)
2702 c
2703 c     open (4,file="HMscore")
2704 c     do k=1,constr_homology
2705 c       read (4,*,end=521) hmscore_tmp
2706 c       hmscore(k)=hmscore_tmp ! Another transformation can be used
2707 c       write(*,*) "Model", k, ":", hmscore(k)
2708 c     enddo
2709 c521  continue
2710
2711 c     write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
2712       do k=1,constr_homology
2713
2714         read(inp,'(a)') pdbfile
2715 c  Next stament causes error upon compilation (?)
2716 c       if(me.eq.king.or. .not. out1file)
2717 c         write (iout,'(2a)') 'PDB data will be read from file ',
2718 c    &   pdbfile(:ilen(pdbfile))
2719         open(ipdbin,file=pdbfile,status='old',err=33)
2720         goto 34
2721   33    write (iout,'(a)') 'Error opening PDB file.'
2722         stop
2723   34    continue
2724 c        print *,'Begin reading pdb data'
2725 c
2726 c Files containing res sim or local scores (former containing sigmas)
2727 c
2728
2729         write(kic2,'(bz,i2.2)') k
2730
2731         tpl_k_rescore="template"//kic2//".sco"
2732 c       tpl_k_sigma_odl="template"//kic2//".sigma_odl"
2733 c       tpl_k_sigma_dih="template"//kic2//".sigma_dih"
2734 c       tpl_k_sigma_theta="template"//kic2//".sigma_theta"
2735 c       tpl_k_sigma_d="template"//kic2//".sigma_d"
2736
2737         call readpdb
2738 c
2739 c     Distance restraints
2740 c
2741 c          ... --> odl(k,ii)
2742 C Copy the coordinates from reference coordinates (?)
2743         do i=1,2*nres
2744           do j=1,3
2745             c(j,i)=cref(j,i)
2746 c           write (iout,*) "c(",j,i,") =",c(j,i)
2747           enddo
2748         enddo
2749 c
2750 c From read_dist_constr (commented out 25/11/2014 <-> res sim)
2751 c
2752 c         write(iout,*) "tpl_k_rescore - ",tpl_k_rescore
2753           open (ientin,file=tpl_k_rescore,status='old')
2754           do irec=1,maxdim ! loop for reading res sim
2755             if (irec.eq.1) then
2756                rescore(k,irec)=0.0d0
2757                goto 1301
2758             endif
2759             read (ientin,*,end=1401) rescore_tmp
2760 c           rescore(k,irec)=rescore_tmp+1.0d0 ! to avoid 0 values
2761             rescore(k,irec)=0.5d0*(rescore_tmp+0.5d0) ! alt transf to reduce scores
2762 c           write(iout,*) "rescore(",k,irec,") =",rescore(k,irec)
2763  1301     continue
2764           enddo
2765  1401   continue
2766           close (ientin)
2767 c         open (ientin,file=tpl_k_sigma_odl,status='old')
2768 c         do irec=1,maxdim ! loop for reading sigma_odl
2769 c            read (ientin,*,end=1401) i, j,
2770 c    &                                sigma_odl_temp(i+nnt-1,j+nnt-1,k) ! new variable (?)
2771 c            sigma_odl_temp(j+nnt-1,i+nnt-1,k)= ! which purpose?
2772 c    &       sigma_odl_temp(i+nnt-1,j+nnt-1,k)
2773 c         enddo
2774 c 1401   continue
2775 c         close (ientin)
2776         if (waga_dist.gt.0.0d0) then
2777           ii=0
2778           do i = nnt,nct-2 ! right? without parallel.
2779             do j=i+2,nct ! right?
2780 c         do i = 1,nres ! alternative for bounds as used to set initial values in orig. read_constr_homology
2781 c           do j=i+2,nres ! ibid
2782 c         do i = nnt,nct-2 ! alternative for bounds as used to assign dist restraints in orig. read_constr_homology (s. above)
2783 c           do j=i+2,nct ! ibid
2784               ii=ii+1
2785 c             write (iout,*) "k",k
2786 c             write (iout,*) "i",i," j",j," constr_homology",
2787 c    &                       constr_homology
2788               ires_homo(ii)=i
2789               jres_homo(ii)=j
2790 c
2791 c Attempt to replace dist(i,j) by its definition in ...
2792 c
2793               x12=c(1,i)-c(1,j)
2794               y12=c(2,i)-c(2,j)
2795               z12=c(3,i)-c(3,j)
2796               distal=dsqrt(x12*x12+y12*y12+z12*z12)
2797               odl(k,ii)=distal
2798 c
2799 c             odl(k,ii)=dist(i,j)
2800 c             write (iout,*) "dist(",i,j,") =",dist(i,j)
2801 c             write (iout,*) "distal = ",distal
2802 c             write (iout,*) "odl(",k,ii,") =",odl(k,ii)
2803 c            write(iout,*) "rescore(",k,i,") =",rescore(k,i),
2804 c    &                      "rescore(",k,j,") =",rescore(k,j)
2805 c
2806 c  Calculation of sigma from res sim
2807 c
2808 c             if (odl(k,ii).le.6.0d0) then
2809 c             sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
2810 c    Other functional forms possible depending on odl(k,ii), eg.
2811 c
2812             if (odl(k,ii).le.dist_cut) then
2813               sigma_odl(k,ii)=rescore(k,i)+rescore(k,j) ! other exprs possible
2814 c             sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)
2815             else
2816               sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))* ! sigma ~ rescore ~ error
2817      &                      dexp(0.5d0*(odl(k,ii)/dist_cut)**2)
2818
2819 c   Following expr replaced by a positive exp argument
2820 c             sigma_odl(k,ii)=(rescore(k,i)+rescore(k,j))*
2821 c    &                      dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
2822
2823 c             sigma_odl(k,ii)=hmscore(k)*rescore(k,i)*rescore(k,j)*
2824 c    &                      dexp(-0.5d0*(odl(k,ii)/dist_cut)**2)
2825             endif
2826 c
2827               sigma_odl(k,ii)=1.0d0/(sigma_odl(k,ii)*sigma_odl(k,ii)) ! rescore ~ error
2828 c             sigma_odl(k,ii)=sigma_odl(k,ii)*sigma_odl(k,ii)
2829 c
2830 c             sigma_odl(k,ii)=sigma_odl_temp(i,j,k)* ! new var read from file (?)
2831 c    &                        sigma_odl_temp(i,j,k)  ! not inverse because of use of res. similarity
2832             enddo
2833 c           read (ientin,*) sigma_odl(k,ii) ! 1st variant
2834           enddo
2835 c         lim_odl=ii
2836 c         if (constr_homology.gt.0) call homology_partition
2837         endif
2838 c
2839 c     Theta, dihedral and SC retraints
2840 c
2841         if (waga_angle.gt.0.0d0) then
2842 c         open (ientin,file=tpl_k_sigma_dih,status='old')
2843 c         do irec=1,maxres-3 ! loop for reading sigma_dih
2844 c            read (ientin,*,end=1402) i,j,ki,l,sigma_dih(k,i+nnt-1) ! j,ki,l what for?
2845 c            if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! right?
2846 c            sigma_dih(k,i+nnt-1)=sigma_dih(k,i+nnt-1)* ! not inverse because of use of res. similarity
2847 c    &                            sigma_dih(k,i+nnt-1)
2848 c         enddo
2849 c1402   continue
2850 c         close (ientin)
2851           do i = nnt+3,nct ! right? without parallel.
2852 c         do i=1,nres ! alternative for bounds acc to readpdb?
2853 c         do i=1,nres-3 ! alternative for bounds as used to set initial values in orig. read_constr_homology
2854 c         do i=idihconstr_start_homo,idihconstr_end_homo ! with FG parallel.
2855             dih(k,i)=phiref(i) ! right?
2856 c           read (ientin,*) sigma_dih(k,i) ! original variant
2857 c             write (iout,*) "dih(",k,i,") =",dih(k,i)
2858 c             write(iout,*) "rescore(",k,i,") =",rescore(k,i),
2859 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
2860 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2),
2861 c    &                      "rescore(",k,i-3,") =",rescore(k,i-3)
2862
2863             sigma_dih(k,i)=rescore(k,i)+rescore(k,i-1)+
2864      &                     rescore(k,i-2)+rescore(k,i-3)  !  right expression ?
2865 c
2866 c           write (iout,*) "Raw sigmas for dihedral angle restraints"
2867 c           write (iout,'(i5,10(2f8.2,4x))') i,sigma_dih(k,i)
2868 c           sigma_dih(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
2869 c                          rescore(k,i-2)*rescore(k,i-3)  !  right expression ?
2870 c   Instead of res sim other local measure of b/b str reliability possible
2871             sigma_dih(k,i)=1.0d0/(sigma_dih(k,i)*sigma_dih(k,i))
2872 c           sigma_dih(k,i)=sigma_dih(k,i)*sigma_dih(k,i)
2873             if (i-nnt-2.gt.lim_dih) lim_dih=i-nnt-2 ! right?
2874 c           if (i+nnt-1.gt.lim_dih) lim_dih=i+nnt-1 ! original when readin i from file
2875           enddo
2876         endif
2877
2878         if (waga_theta.gt.0.0d0) then
2879 c         open (ientin,file=tpl_k_sigma_theta,status='old')
2880 c         do irec=1,maxres-2 ! loop for reading sigma_theta, right bounds?
2881 c            read (ientin,*,end=1403) i,j,ki,sigma_theta(k,i+nnt-1) ! j,ki what for?
2882 c            sigma_theta(k,i+nnt-1)=sigma_theta(k,i+nnt-1)* ! not inverse because of use of res. similarity
2883 c    &                              sigma_theta(k,i+nnt-1)
2884 c         enddo
2885 c1403   continue
2886 c         close (ientin)
2887
2888           do i = nnt+2,nct ! right? without parallel.
2889 c         do i = i=1,nres ! alternative for bounds acc to readpdb?
2890 c         do i=ithet_start,ithet_end ! with FG parallel.
2891              thetatpl(k,i)=thetaref(i)
2892 c            write (iout,*) "thetatpl(",k,i,") =",thetatpl(k,i)
2893 c            write(iout,*)  "rescore(",k,i,") =",rescore(k,i),
2894 c    &                      "rescore(",k,i-1,") =",rescore(k,i-1),
2895 c    &                      "rescore(",k,i-2,") =",rescore(k,i-2)
2896 c            read (ientin,*) sigma_theta(k,i) ! 1st variant
2897              sigma_theta(k,i)=rescore(k,i)+rescore(k,i-1)+
2898      &                        rescore(k,i-2) !  right expression ?
2899              sigma_theta(k,i)=1.0d0/(sigma_theta(k,i)*sigma_theta(k,i))
2900
2901 c            sigma_theta(k,i)=hmscore(k)*rescore(k,i)*rescore(k,i-1)*
2902 c                             rescore(k,i-2) !  right expression ?
2903 c            sigma_theta(k,i)=sigma_theta(k,i)*sigma_theta(k,i)
2904              if (i-nnt-1.gt.lim_theta) lim_theta=i-nnt-1 ! right?
2905           enddo
2906         endif
2907
2908         if (waga_d.gt.0.0d0) then
2909 c       open (ientin,file=tpl_k_sigma_d,status='old')
2910 c         do irec=1,maxres-1 ! loop for reading sigma_theta, right bounds?
2911 c            read (ientin,*,end=1404) i,j,sigma_d(k,i+nnt-1) ! j,ki what for?
2912 c            sigma_d(k,i+nnt-1)=sigma_d(k,i+nnt-1)* ! not inverse because of use of res. similarity
2913 c    &                          sigma_d(k,i+nnt-1)
2914 c         enddo
2915 c1404   continue
2916           close (ientin)
2917
2918           do i = nnt,nct ! right? without parallel.
2919 c         do i=2,nres-1 ! alternative for bounds acc to readpdb?
2920 c         do i=loc_start,loc_end ! with FG parallel.
2921              if (itype(i).eq.10) goto 1 ! right?
2922                xxtpl(k,i)=xxref(i)
2923                yytpl(k,i)=yyref(i)
2924                zztpl(k,i)=zzref(i)
2925 c              write (iout,*) "xxtpl(",k,i,") =",xxtpl(k,i)
2926 c              write (iout,*) "yytpl(",k,i,") =",yytpl(k,i)
2927 c              write (iout,*) "zztpl(",k,i,") =",zztpl(k,i)
2928 c              write(iout,*)  "rescore(",k,i,") =",rescore(k,i)
2929                sigma_d(k,i)=rescore(k,i) !  right expression ?
2930                sigma_d(k,i)=1.0d0/(sigma_d(k,i)*sigma_d(k,i))
2931
2932 c              sigma_d(k,i)=hmscore(k)*rescore(k,i) !  right expression ?
2933 c              sigma_d(k,i)=sigma_d(k,i)*sigma_d(k,i)
2934 c              read (ientin,*) sigma_d(k,i) ! 1st variant
2935                if (i-nnt+1.gt.lim_xx) lim_xx=i-nnt+1 ! right?
2936     1     continue
2937           enddo
2938         endif
2939         close(ientin)
2940       enddo
2941       if (waga_dist.gt.0.0d0) lim_odl=ii
2942       if (constr_homology.gt.0) call homology_partition
2943       if (constr_homology.gt.0) call init_int_table
2944       write (iout,*) "homology_partition: lim_theta= ",lim_theta,
2945      & "lim_xx=",lim_xx
2946 c     write (iout,*) "ithet_start =",ithet_start,"ithet_end =",ithet_end
2947 c     write (iout,*) "loc_start =",loc_start,"loc_end =",loc_end
2948 c
2949 c Print restraints
2950 c
2951       if (.not.lprn) return
2952       write(iout,*) "waga_theta",waga_theta,"waga_d",waga_d
2953       write (iout,*) "Distance restraints from templates"
2954       do ii=1,lim_odl
2955        write(iout,'(3i5,10(2f16.2,4x))') ii,ires_homo(ii),jres_homo(ii),
2956      &  (odl(ki,ii),1.0d0/dsqrt(sigma_odl(ki,ii)),ki=1,constr_homology)
2957       enddo
2958       write (iout,*) "Dihedral angle restraints from templates"
2959       do i=nnt+3,lim_dih
2960         write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*dih(ki,i),
2961      &      rad2deg/dsqrt(sigma_dih(ki,i)),ki=1,constr_homology)
2962       enddo
2963       write (iout,*) "Virtual-bond angle restraints from templates"
2964       do i=nnt+2,lim_theta
2965         write (iout,'(i5,10(2f8.2,4x))') i,(rad2deg*thetatpl(ki,i),
2966      &      rad2deg/dsqrt(sigma_theta(ki,i)),ki=1,constr_homology)
2967       enddo
2968       write (iout,*) "SC restraints from templates"
2969       do i=nnt,lim_xx
2970         write(iout,'(i5,10(4f8.2,4x))') i,
2971      &  (xxtpl(ki,i),yytpl(ki,i),zztpl(ki,i),
2972      &   1.0d0/dsqrt(sigma_d(ki,i)),ki=1,constr_homology)
2973       enddo
2974
2975 c -----------------------------------------------------------------
2976       return
2977       end
2978 c----------------------------------------------------------------------
2979
2980 #ifdef WINIFL
2981       subroutine flush(iu)
2982       return
2983       end
2984 #endif
2985 #ifdef AIX
2986       subroutine flush(iu)
2987       call flush_(iu)
2988       return
2989       end
2990 #endif
2991
2992 c------------------------------------------------------------------------------
2993       subroutine copy_to_tmp(source)
2994       include "DIMENSIONS"
2995       include "COMMON.IOUNITS"
2996       character*(*) source
2997       character* 256 tmpfile
2998       integer ilen
2999       external ilen
3000       logical ex
3001       tmpfile=curdir(:ilen(curdir))//"/"//source(:ilen(source))
3002       inquire(file=tmpfile,exist=ex)
3003       if (ex) then
3004         write (*,*) "Copying ",tmpfile(:ilen(tmpfile)),
3005      &   " to temporary directory..."
3006         write (*,*) "/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir
3007         call system("/bin/cp "//tmpfile(:ilen(tmpfile))//" "//tmpdir)
3008       endif
3009       return
3010       end
3011 c------------------------------------------------------------------------------
3012       subroutine move_from_tmp(source)
3013       include "DIMENSIONS"
3014       include "COMMON.IOUNITS"
3015       character*(*) source
3016       integer ilen
3017       external ilen
3018       write (*,*) "Moving ",source(:ilen(source)),
3019      & " from temporary directory to working directory"
3020       write (*,*) "/bin/mv "//source(:ilen(source))//" "//curdir
3021       call system("/bin/mv "//source(:ilen(source))//" "//curdir)
3022       return
3023       end
3024 c------------------------------------------------------------------------------
3025       subroutine random_init(seed)
3026 C
3027 C Initialize random number generator
3028 C
3029       implicit real*8 (a-h,o-z)
3030       include 'DIMENSIONS'
3031 #ifdef AMD64
3032       integer*8 iseedi8
3033 #endif
3034 #ifdef MPI
3035       include 'mpif.h'
3036       logical OKRandom, prng_restart
3037       real*8  r1
3038       integer iseed_array(4)
3039 #endif
3040       include 'COMMON.IOUNITS'
3041       include 'COMMON.TIME1'
3042       include 'COMMON.THREAD'
3043       include 'COMMON.SBRIDGE'
3044       include 'COMMON.CONTROL'
3045       include 'COMMON.MCM'
3046       include 'COMMON.MAP'
3047       include 'COMMON.HEADER'
3048 csa      include 'COMMON.CSA'
3049       include 'COMMON.CHAIN'
3050       include 'COMMON.MUCA'
3051       include 'COMMON.MD'
3052       include 'COMMON.FFIELD'
3053       include 'COMMON.SETUP'
3054       iseed=-dint(dabs(seed))
3055       if (iseed.eq.0) then
3056         write (iout,'(/80(1h*)/20x,a/80(1h*))')
3057      &    'Random seed undefined. The program will stop.'
3058         write (*,'(/80(1h*)/20x,a/80(1h*))')
3059      &    'Random seed undefined. The program will stop.'
3060 #ifdef MPI
3061         call mpi_finalize(mpi_comm_world,ierr)
3062 #endif
3063         stop 'Bad random seed.'
3064       endif
3065 #ifdef MPI
3066       if (fg_rank.eq.0) then
3067       seed=seed*(me+1)+1
3068 #ifdef AMD64
3069       iseedi8=dint(seed)
3070       if(me.eq.king .or. .not. out1file)
3071      &  write (iout,*) 'MPI: node= ', me, ' iseed= ',iseedi8
3072       write (*,*) 'MPI: node= ', me, ' iseed= ',iseedi8
3073       OKRandom = prng_restart(me,iseedi8)
3074 #else
3075       do i=1,4
3076        tmp=65536.0d0**(4-i)
3077        iseed_array(i) = dint(seed/tmp)
3078        seed=seed-iseed_array(i)*tmp
3079       enddo
3080       if(me.eq.king .or. .not. out1file)
3081      & write (iout,*) 'MPI: node= ', me, ' iseed(4)= ',
3082      &                 (iseed_array(i),i=1,4)
3083       write (*,*) 'MPI: node= ',me, ' iseed(4)= ',
3084      &                 (iseed_array(i),i=1,4)
3085       OKRandom = prng_restart(me,iseed_array)
3086 #endif
3087       if (OKRandom) then
3088         r1=ran_number(0.0D0,1.0D0)
3089         if(me.eq.king .or. .not. out1file)
3090      &   write (iout,*) 'ran_num',r1
3091         if (r1.lt.0.0d0) OKRandom=.false.
3092       endif
3093       if (.not.OKRandom) then
3094         write (iout,*) 'PRNG IS NOT WORKING!!!'
3095         print *,'PRNG IS NOT WORKING!!!'
3096         if (me.eq.0) then
3097          call flush(iout)
3098          call mpi_abort(mpi_comm_world,error_msg,ierr)
3099          stop
3100         else
3101          write (iout,*) 'too many processors for parallel prng'
3102          write (*,*) 'too many processors for parallel prng'
3103          call flush(iout)
3104          stop
3105         endif
3106       endif
3107       endif
3108 #else
3109       call vrndst(iseed)
3110       write (iout,*) 'ran_num',ran_number(0.0d0,1.0d0)
3111 #endif
3112       return
3113       end