[gmx-users] Check for Vacuum Bubbles & Density

[Johnny Lu]

Hi.

The NPT simulation has 9K TIP3P water with 173 amino acid residues and uses
Amber99SB-ILDN force field, and use Berendsen barostat and V-rescale
thermostat. I plan to run NVT with the same pressure an temperature after
this.

The average density is 1015 kg/m^3 over 35ns. The density was near 1015
over more than 300 ns, in previous NPT simulations.

    Statistics over 34744941 steps [ 0.0000 through 34744.9400 ps ], 5 data
sets
    All statistics are over 6948989 points

    Energy                      Average   Err.Est.       RMSD  Tot-Drift

-------------------------------------------------------------------------------
    Total Energy                -307197        210     880.41   -1408.29
(kJ/mol)
    Temperature                 299.993     0.0037    1.70813  0.0273197
(K)
    Pressure                    1.00317      0.009    138.426 -0.0381886
(bar)
    Volume                      298.567      0.049   0.371817  -0.314958
(nm^3)
    Density                      1014.8       0.17    1.26377    1.07078
(kg/m^3)

This is higher than the 982.3 value in the mail-list [
https://www.mail-archive.com/gmx-users@gromacs.org/msg29327.html]. Why the
density is different?

I used a 1.0 nm cut off, following the force-field paper: Improved
side-chain torsion potentials for the Amber ff99SB protein force field [
http://www.ncbi.nlm.nih.gov/pmc/articles/PMC2970904/], but I also added
dispersion correction. Should I not use dispersion correction?

I look at the trajectory in VMD with periodic images turn on. But there
were many water that overlap each other, and I can not see any vacuum
bubble clearly  (if there is a vacuum bubble) [
http://archive.ambermd.org/201103/0431.html]. Is there a better way to
check for vacuum bubbles?

I didn't know a long NVT equilibration at 300 K, right after minimization,
can make vacuum bubbles. The NVT equilibration was 300 ps.

I also didn't know that I should run berendsen barostat before
Parrinello-Rahman, so those 300+ ns NPT simulations before this one all
uses
Parrinello-Rahman barostat.

If there is no vacuum bubbles, I think I was just lucky.

Below is the mdp file I used for the NPT simulation:
    title           = NPT production
    ; Run parameters
    integrator      = md-vv         ; leap-frog integrator
    nsteps          = 36000000      ; 1 ns
    dt              = 0.001 ; 1 fs

    ; Output control
    nstenergy       = 10           ; save energies every 1.0 ps

    nstlog          = 5000          ; update log file every 0.2 ps
    ; Bond parameters
    continuation    = yes           ; Restarting after NVT
    constraint_algorithm = lincs    ; holonomic constraints
    constraints     = H-bonds       ; all bonds (even heavy atom-H bonds)
constrained
    lincs_iter      = 2             ; accuracy of LINCS
    lincs_order     = 4             ; also related to accuracy
    ; Neighborsearching
    ns_type         = grid          ; search neighboring grid cells
    nstlist         = 10            ; 10 fs
    rlist           = 1.0           ; short-range neighborlist cutoff (in
nm)
    rcoulomb        = 1.0           ; short-range electrostatic cutoff (in
nm)
    rvdw            = 1.0           ; short-range van der Waals cutoff (in
nm)
    ; Electrostatics
    coulombtype     = PME           ; Particle Mesh Ewald for long-range
electrostatics
    pme_order       = 4             ; cubic interpolation
    fourierspacing  = 0.1           ; grid spacing for FFT
    ; Temperature coupling is on
    tcoupl          = V-rescale     ; modified Berendsen thermostat
    tc-grps         = system        ; single thermostat
    tau_t           = 1.0           ; time constant, in ps
    ref_t           = 300           ; reference temperature, one for each
group, in K
    ; Pressure coupling is on
    pcoupl          = berendsen     ; Pressure coupling on in NPT
    pcoupltype      = isotropic     ; uniform scaling of box vectors
    tau_p           = 2.0           ; time constant, in ps
    ref_p           = 1.0           ; reference pressure, in bar
    compressibility = 4.5e-5        ; isothermal compressibility of water,
bar^-1
    refcoord_scaling = com
    ; Periodic boundary conditions
    pbc             = xyz           ; 3-D PBC
    ; Dispersion correction
    DispCorr        = EnerPres      ; account for cut-off vdW scheme
    ; Velocity generation
    gen_vel         = no            ; Velocity generation is off
    cutoff-scheme = Verlet

Thank you again.