Pressure Control

Constant pressure simulation (and pressure calculation) require periodic boundary conditions. Pressure is controlled by dynamically adjusting the size of the unit cell and rescaling all atomic coordinates (other than those of fixed atoms) during the simulation.

Pressure values in NAMD output are in bar. PRESSURE is the pressure calculated based on individual atoms, while GPRESSURE incorporates hydrogen atoms into the heavier atoms to which they are bonded, producing smaller fluctuations. The TEMPAVG, PRESSAVG, and GPRESSAVG are the average of temperature and pressure values since the previous ENERGY output; for the first step in the simulation they will be identical to TEMP, PRESSURE, and GPRESSURE.

The phenomenological pressure of bulk matter reflects averaging in both space and time of the sum of a large positive term (the kinetic pressure, $nRT/V$ ), and a large cancelling negative term (the static pressure). The instantaneous pressure of a simulation cell as simulated by NAMD will have mean square fluctuations (according to David Case quoting Section 114 of Statistical Physics by Landau and Lifshitz) of $kT/(V \beta)$ , where $\beta$ is the compressibility, which is RMS of roughly 100 bar for a 10,000 atom biomolecular system. Much larger fluctuations are regularly observed in practice.

The instantaneous pressure for a biomolecular system is well defined for ``internal'' forces that are based on particular periodic images of the interacting atoms, conserve momentum, and are translationally invariant. When dealing with externally applied forces such as harmonic constraints, fixed atoms, and various steering forces, NAMD bases its pressure calculation on the relative positions of the affected atoms in the input coordinates and assumes that the net force will average to zero over time. For time periods during with the net force is non-zero, the calculated pressure fluctuations will include a term proportional to the distance to the affected from the user-defined cell origin. A good way to observe these effects and to confirm that pressure for external forces is handled reasonably is to run a constant volume cutoff simulation in a cell that is larger than the molecular system by at least the cutoff distance; the pressure for this isolated system should average to zero over time.

Because NAMD's impluse-basd multiple timestepping system alters the balance between bonded and non-bonded forces from every timestep to an average balance over two steps, the calculated pressure on even and odd steps will be different. The PRESSAVG and GPRESSAVG fields provide the average over the non-printed intermediate steps. If you print energies on every timestep you will see the effect clearly in the PRESSURE field.

The following options affect all pressure control methods.

• useGroupPressure $<$ group or atomic quantities $>$
  Acceptable Values: yes or no
  Default Value: no
  Description: Pressure can be calculated using either the atomic virial and kinetic energy (the default) or a hydrogen-group based pseudo-molecular virial and kinetic energy. The latter fluctuates less and is required in conjunction with rigidBonds (SHAKE).

• useFlexibleCell $<$ anisotropic cell fluctuations $>$
  Acceptable Values: yes or no
  Default Value: no
  Description: NAMD allows the three orthogonal dimensions of the periodic cell to fluctuate independently when this option is enabled.

• useConstantRatio $<$ constant shape in first two cell dimensions $>$
  Acceptable Values: yes or no
  Default Value: no
  Description: When enabled, NAMD keeps the ratio of the unit cell in the x-y plane constant while allowing fluctuations along all axes. The useFlexibleCell option is required for this option.

• useConstantArea $<$ constant area and normal pressure conditions $>$
  Acceptable Values: yes or no
  Default Value: no
  Description: When enabled, NAMD keeps the dimension of the unit cell in the x-y plane constant while allowing fluctuations along the z axis. This is not currently implemented in Berendsen's method.

Berendsen pressure bath coupling

NAMD provides constant pressure simulation using Berendsen's method. The following parameters are used to define the algorithm.

• BerendsenPressure $<$ use Berendsen pressure bath coupling? $>$
  Acceptable Values: on or off
  Default Value: off
  Description: Specifies whether or not Berendsen pressure bath coupling is active. If set to on, then the parameters BerendsenPressureTarget, BerendsenPressureCompressibility and BerendsenPressureRelaxationTime must be set and the parameter BerendsenPressureFreq can optionally be set to control the behavior of this feature.

• BerendsenPressureTarget $<$ target pressure (bar) $>$
  Acceptable Values: positive decimal
  Description: Specifies target pressure for Berendsen's method. A typical value would be 1.01325 bar, atmospheric pressure at sea level.

• BerendsenPressureCompressibility $<$ compressibility (bar$^{-1}$ ) $>$
  Acceptable Values: positive decimal
  Description: Specifies compressibility for Berendsen's method. A typical value would be 4.57E-5 bar$^{-1}$ , corresponding to liquid water. The higher the compressibility, the more volume will be adjusted for a given pressure difference. The compressibility and the relaxation time appear only as a ratio in the dynamics, so a larger compressibility is equivalent to a smaller relaxation time.

• BerendsenPressureRelaxationTime $<$ relaxation time (fs) $>$
  Acceptable Values: positive decimal
  Description: Specifies relaxation time for Berendsen's method. If the instantaneous pressure did not fluctuate randomly during a simulation and the compressibility estimate was exact then the inital pressure would decay exponentially to the target pressure with this time constant. Having a longer relaxation time results in more averaging over pressure measurements and hence smaller fluctuations in the cell volume. A reasonable choice for relaxation time would be 100 fs. The compressibility and the relaxation time appear only as a ratio in the dynamics, so a larger compressibility is equivalent to a smaller relaxation time.

• BerendsenPressureFreq $<$ how often to rescale positions $>$
  Acceptable Values: positive multiple of nonbondedFrequency and fullElectFrequency
  Default Value: nonbondedFrequency or fullElectFrequency if used
  Description: Specifies number of timesteps between position rescalings for Berendsen's method. Primarily to deal with multiple timestepping integrators, but also to reduce cell volume fluctuations, cell rescalings can occur on a longer interval. This could reasonably be between 1 and 20 timesteps, but the relaxation time should be at least ten times larger.

Nosé-Hoover Langevin piston pressure control

NAMD provides constant pressure simulation using a modified Nosé-Hoover method in which Langevin dynamics is used to control fluctuations in the barostat. This method should be combined with a method of temperature control, such as Langevin dynamics, in order to simulate the NPT ensemble.

The Langevin piston Nose-Hoover method in NAMD is a combination of the Nose-Hoover constant pressure method as described in GJ Martyna, DJ Tobias and ML Klein, "Constant pressure molecular dynamics algorithms", J. Chem. Phys 101(5), 1994, with piston fluctuation control implemented using Langevin dynamics as in SE Feller, Y Zhang, RW Pastor and BR Brooks, "Constant pressure molecular dynamics simulation: The Langevin piston method", J. Chem. Phys. 103(11), 1995.

The equations of motion are:

$$r' = p/m + e' r$$

$$p' = F - e' p - g p + R$$

$$V' = 3 V e'$$

$$e'' = 3V/W (P - P_0) - g_e e' + R_e/W$$

$$W = 3 N \tau^2 k T$$

$$<R^2> = 2 m g k T / h$$

$$\tau = {\rm oscillation period}$$

$$<R_e^2> = 2 W g_e k T / h$$

Here, $W$ is the mass of piston, $R$ is noise on atoms, and $R_e$ is the noise on the piston.

The user specifies the desired pressure, oscillation and decay times of the piston, and temperature of the piston. The compressibility of the system is not required. In addition, the user specifies the damping coefficients and temperature of the atoms for Langevin dynamics.

The following parameters are used to define the algorithm.

• LangevinPiston $<$ use Langevin piston pressure control? $>$
  Acceptable Values: on or off
  Default Value: off
  Description: Specifies whether or not Langevin piston pressure control is active. If set to on, then the parameters LangevinPistonTarget, LangevinPistonPeriod, LangevinPistonDecay and LangevinPistonTemp must be set.

• LangevinPistonTarget $<$ target pressure (bar) $>$
  Acceptable Values: positive decimal
  Description: Specifies target pressure for Langevin piston method. A typical value would be 1.01325 bar, atmospheric pressure at sea level.

• LangevinPistonPeriod $<$ oscillation period (fs) $>$
  Acceptable Values: positive decimal
  Description: Specifies barostat oscillation time scale for Langevin piston method. If the instantaneous pressure did not fluctuate randomly during a simulation and the decay time was infinite (no friction) then the cell volume would oscillate with this angular period. Having a longer period results in more averaging over pressure measurements and hence slower fluctuations in the cell volume. A reasonable choice for the piston period would be 200 fs.

• LangevinPistonDecay $<$ damping time scale (fs) $>$
  Acceptable Values: positive decimal
  Description: Specifies barostat damping time scale for Langevin piston method. A value larger than the piston period would result in underdamped dynamics (decaying ringing in the cell volume) while a smaller value approaches exponential decay as in Berendsen's method above. A smaller value also corresponds to larger random forces with increased coupling to the Langevin temperature bath. Typically this would be chosen equal to or smaller than the piston period, such as 100 fs.

• LangevinPistonTemp $<$ noise temperature (K) $>$
  Acceptable Values: positive decimal
  Description: Specifies barostat noise temperature for Langevin piston method. This should be set equal to the target temperature for the chosen method of temperature control.

• SurfaceTensionTarget $<$ Surface tension target (dyn/cm) $>$
  Acceptable Values: decimal
  Default Value: 0.0
  Description: Specifies surface tension target. Must be used with useFlexibleCell and periodic boundary conditions. The pressure specified in LangevinPistonTarget becomes the pressure along the z axis, and surface tension is applied in the x-y plane.

• StrainRate $<$ initial strain rate $>$
  Acceptable Values: decimal triple (x y z)
  Default Value: 0. 0. 0.
  Description: Optionally specifies the initial strain rate for pressure control. Is overridden by value read from file specified with extendedSystem. There is typically no reason to set this parameter.

• ExcludeFromPressure $<$ Should some atoms be excluded from pressure rescaling? $>$
  Acceptable Values: on or off
  Default Value: off
  Description: Specifies whether or not to exclude some atoms from pressure rescaling. The coordinates and velocites of such atoms are not rescaled during constant pressure simulations, though they do contribute to the virial calculation. May be useful for membrane protein simulation. EXPERIMENTAL.

• ExcludeFromPressureFile $<$ File specifying excluded atoms $>$
  Acceptable Values: PDB file
  Default Value: coordinates file
  Description: PDB file with one column specifying which atoms to exclude from pressure rescaling. Specify 1 for excluded and 0 for not excluded.

• ExcludeFromPressureCol $<$ Column in PDB file for specifying excluded atoms $>$
  Acceptable Values: O, B, X, Y, or Z
  Default Value: O
  Description: Specifies which column of the pdb file to check for excluded atoms.

Monte Carlo barostat

NAMD provides constant pressure simulation using Monte Carlo method to control volume fluctuations in the barostat. This method should be combined with a method of temperature control, in order to simulate the NPT ensemble. This feature is only supported in GPU Resident mode.

Please note that before switching to constant pressure simulations, you will need to equilibrate your system using NVT ensemble simulation.

In this method, a trial volume-change, $\Delta V$ , is generated uniformly, within the range of $[-\Delta V_{\rm max}, +\Delta V_{\rm max}]$ , where $\Delta V_{\rm max}$ is the maximum volume-change, adjusted to obtain the target acceptance ratio.

The probability of accepting volume-change from $V$ to $V^{'} = V + \Delta V$ is:

$$P^{acc} = {\rm min} \left[1, \left(\frac{V^{'}}{V}\right)^{N} e^{-\beta \left(\Delta U + P\Delta V - \gamma \Delta A\right)} \right]$$

$$\Delta U = U\left(\textbf{s}^{N}, V^{'}\right) - U\left(\textbf{s}^{N}, V\right)$$

$$\Delta V = V^{'} - V$$

$$\Delta A = {\textrm{change in x-y area}}$$

$$\beta = \frac{1}{k T}$$

$$\textbf{s}^{N} = {\textrm{reduced coordinates}}$$

where, $\gamma$ , $P$ , $T$ , $U\left(\textbf{s}^{N}, V^{'}\right)$ , and $U\left(\textbf{s}^{N}, V\right)$ are the target surface tension, target pressure, system temperature, and total potential energy of the system at new and old configuration, respectively. $N$ denotes the total number of molecule in the system.

The following parameters are used to define the algorithm.

• MonteCarloPressure $<$ use Monte Carlo pressure control? $>$
  Acceptable Values: on or off
  Default Value: off
  Description: Specifies whether or not Monte Carlo pressure control is active. If set to on, then the parameters MonteCarloPressureTarget and MonteCarloTemp must be set. The parameters MonteCarloPressureFreq, MonteCarloMaxVolume, MonteCarloAcceptanceRate, MonteCarloAdjustmentFreq, and SurfaceTensionTarget can optionally be set to control the behavior of this feature.

• MonteCarloPressureTarget $<$ target pressure (bar) $>$
  Acceptable Values: positive decimal
  Description: Specifies target pressure for Monte Carlo barostat. A typical value would be 1.01325 bar, atmospheric pressure at sea level.

• MonteCarloTemp $<$ system temperature (K) $>$
  Acceptable Values: positive decimal
  Description: Specifies temperature for Monte Carlo barostat. This should be set equal to the target temperature for the chosen method of temperature control.

• MonteCarloPressureFreq $<$ how often attempt to change volume? $>$
  Acceptable Values: positive multiple of nonbondedFrequency and fullElectFrequency
  Default Value: 50 * nonbondedFrequency or 50 * fullElectFrequency if used
  Description: Specifies number of timesteps between attempt to change the volume and rescale the coordinates. A typical value for this parameter could be 50.

• MonteCarloMaxVolume $<$ initial maximum volume-change for each axis dimension (Å$^3$ ) $>$
  Acceptable Values: vector of positive decimal and less than $0.3V$
  Default Value: $0.01V$ $0.01V$ $0.01V$
  Description: Optionally specifies the initial maximum volume-change for each axis dimension ( $\Delta V_{\rm max}$ ) in Monte Carlo pressure control. If extendedSystem is also specified, value is overridden. There is typically no reason to set this parameter.

• MonteCarloAcceptanceRate $<$ target volume-change acceptance rate $>$
  Acceptable Values: decimal between 0.2 and 0.8
  Default Value: 0.5
  Description: Specifies the desired acceptance rate for volume-change attempt. The $\Delta V_{\rm max}$ is dynamically adjusted to reach the target acceptance rate. Higher acceptance rate results in small volume fluctuations, while lower acceptance rate results in large volume fluctuations.

• MonteCarloAdjustmentFreq $<$ how often adjust maximum volume-change? $>$
  Acceptable Values: positive integer
  Default Value: 30
  Description: Specifies the number of volume-change attempts, to collect acceptance rate statistics, in order to adjust maximum volume-change ( $\Delta V_{\rm max}$ ). Higher value would provide more accurate statistics, but might require longer simulation to reach equilibration. A typical value of 30 for isotropic and 50 for anisotropic cell fluctuation could be used for this parameter.

• SurfaceTensionTarget $<$ Surface tension target (dyn/cm) $>$
  Acceptable Values: decimal
  Default Value: 0.0
  Description: Specifies surface tension target. Must be used with useFlexibleCell and periodic boundary conditions. The pressure specified in MonteCarloPressureTarget becomes the pressure along the z axis, and surface tension is applied in the x-y plane.