Pressure Profile Calculations

NAMD supports the calculation of lateral pressure profiles as a function of the z-coordinate in the system. The algorithm is essentially that described by Lindahl and Edholm. The simulation space is partitioned into slabs, and the virial due to the interaction between two particles is distributed across the slabs containing the particles as well as the slabs that lie between the particles (taking periodic boundary conditions into account). The diagonal components of the pressure tensor for each slab are recorded in the NAMD output file. The units of pressure are the same as in the regular NAMD pressure output; i.e., bar.

The total virial contains contributions from three components: kinetic energy, bonded interactions, and nonbonded interactions. The kinetic and bonded components are easily calculated every timestep, and thus may be computed during a normal simulation run. The nonbonded component, however, adds significant overhead and has not been implemented for PME calculations. The calculation of the nonbonded contribution should be performed offline, using the saved frames of the trajectory file and a long nonbonded cutoff.

Pressure profile calculations may be performed in either constant volume or constant pressure conditions. If constant pressure is enabled, the slabs thickness will be rescaled along with the unit cell; the dcdUnitCell option will also be switched on so that unit cell information is stored in the trajectory file.

Periodic boundary conditions must also be enabled.

pressureProfile $<$ compute pressure profile $>$
Acceptable Values: on or off
Default Value: off
Description: When active, NAMD will compute kinetic and bonded contributions to the pressure profile. Results will be recorded in the NAMD output file in lines with the format PRESSUREPROFILE: ts Axx Ayy Azz Bxx Byy Bzz ... , where ts is the timestep, followed by the three diagonal components of the pressure tensor in the first slab (the slab with lowest z), then the next lowest slab, and so forth. NAMD will also output the pressure profile averaged over all the steps since the last output. The format of this line is PRESSUREPROFILEAVG: ts Axx Ayy Azz ... ; i.e., exactly as for the instantaneous pressure. It is recommended to use the averaged pressure profile instead of the instantaneous value as this will give better statistics and may prevent artifacts when using multiple timestepping. pressureProfileSlabs $<$ Number of slabs in the spatial partition $>$
Acceptable Values: Positive integer
Default Value: 10
Description: NAMD divides the entire periodic cell into horizontal slabs of equal thickness; pressureProfileSlabs specifies the number of such slabs.

pressureProfileFreq $<$ How often to output pressure profile data $>$
Acceptable Values: Positive integer
Default Value: 1
Description: Specifies the number of timesteps between output of pressure profile data.

pressureProfileNonbonded $<$ Compute nonbonded pressure profile contribution $>$
Acceptable Values: on or off
Default Value: off
Description: When enabled, only the nonbonded contribution to the pressure profile will be computed. For trajectory analysis the recommended way to use this option is to use the NAMD Tcl scripting interface as described in Sec. 2.2.2 to run for 0 steps, so that NAMD prints the pressure profile without performing any dynamics.

Here is an example snippet from a NAMD input that can be used to compute the nonbonded component of the pressure profile. It assumes that the coordinates were saved in three dcd files (pp01.dcd, pp02.dcd, and pp03.dcd) every 500 timesteps. The pressureProfileSlabs must be the same as was used for the calculation of the bonded part of the pressure. Note the use of no switching function; this way the calculated electrostatic energy will be closest to what would be calculated by PME.

switching	     off
cutoff		     20.0
pairlistdist   20.5

pressureProfile	on
pressureProfileSlabs 60 
pressureProfileFreq 100 
pressureProfileNonbonded on

# Assume that coordinates were written to the dcd files every 500 timesteps.
set ts 500
foreach dcd { pp01.dcd pp02.dcd pp03.dcd } {
  coorfile open dcd $dcd
  while { [coorfile read] != -1 } {
    firstTimestep $ts 
    run 0 
    incr ts 500
  }
  coorfile close
}