A LAMMPS package for molecular dynamics under extensional flow fields
Note: As of 10/23/2017, this package has been included in LAMMPS as the USER-UEF package. Installation and usage details can be found here.
UEF is a LAMMPS package for non-equilibrium molecular dynamics (NEMD) under diagonal flow fields, including uniaxial and biaxial flow. With this package, simulations under extensional flow may be carried out for an indefinite amount of time. It is an implementation of the boundary conditions developed by Matthew Dobson, and also uses numerical lattice reduction as was proposed by Thomas Hunt. The lattice reduction algorithm is from Igor Semaev. The package is intended for simulations of homogeneous flows, and integrates the SLLOD equations of motion.
Authored by:
David Nicholson
Massachusetts Institute of Technology
Initial commit: Oct 19, 2016
Support provided via issues and/or email.
- Installation
- Usage
- fix nvt/uef
- fix npt/uef
- compute temp/uef
- compute pressure/uef
- dump cfg/uef
- Implementation Details
- Examples
- Error and Warning Messages
- Citing the UEF package
The implementation has been tested in the Jul. 30, 2016 stable version of LAMMPS. Older versions may not be compatible.
The UEF package is compiled within LAMMPS, and doesn't require any additional libraries or modification to the LAMMPS source code.
To install the package from a LAMMPS distribution located at lammps/
- Clone this repository and copy the
USER-UEF/
subdirectory into thelammps/src/
directory. - Move into
lammps/src/
then runmake yes-USER-UEF
. - Compile LAMMPS. See the LAMMPS documentation for compilation instructions.
The following commands will install the package in a fresh version of the current stable LAMMPS package to the current directory.
wget http://lammps.sandia.gov/tars/lammps-stable.tar.gz
tar -xvf lammps-stable.tar.gz
git clone https://github.com/RutledgeGroupMIT/UEF.git
cp -r UEF/USER-UEF/ lammps-*/src
cd lammps-*/src/
make yes-user-uef
[compile LAMMPS (e.g. make mpi)]
The package defines fix nvt/uef
and fix npt/uef
for constant volume and stress-controlled simulations, compute pressure/uef
and compute temp/uef
to compute the pressure and kinetic energy tensors, and dump cfg/uef
for outputting properly oriented atomic coordinates.
fix ID all nvt/uef temp Tstart Tstop Tdamp erate eps_x eps_y keyword values ...
- ID = name for the fix
- Tstart,Tstop = external temperature at start/end of run
- Tdamp = temperature damping parameter (time units)
- eps_x = strain rate in x dimension 1/(time units)
- eps_y = strain rate in y dimension 1/(time units)
Additional keywords: - strain = initial level of strain (default="0 0"). Use of this keyword is not recommended, but may be necessary when resuming a run from a data file. This keyword should be leftunse when restart files are used.
- The following additional keywords from
fix nvt
can be used with this fix: tchain, tloop, drag
- Uniaxial flow
fix f1 all nvt/uef temp 400 400 100 erate 0.00001 -0.000005
- Biaxial flow
fix f2 all nvt/uef temp 400 400 100 erate 0.000005 0.000005
The applied flow field must be traceless, and therefore only the xx and yy components of the rate-of-deformation tensor need to be specified ( by eps_x
and eps_y
respectively). The zz component of the rate-of-deformation tensor is always -(eps_x
+ eps_y
).
Due to requirements of the boundary conditions, when the strain keyword is unset, or set to zero, the initial simulation box must be cubic and have style triclinic. If the box is initially of type ortho, use the command change box all triclinic
before invoking the fix.
This fix integrates the SLLOD equations of motion, which lead to an instability in the center of mass velocity under extension. A fix momentum
should be used to regularly reset the linear momentum. Additionally, this fix stores the peculiar velocity of each atom, defined as the velocity relative to the streaming velocity. This is in contrast to the LAMMPS fix nvt/sllod
command, which uses a lab-frame velocity.
This fix defines a compute pressure/uef
and compute temp/uef
that can be accessed at c_ID_press
and c_ID_temp
respectively for scalar values, or c_ID_press[i]
and c_ID_temp[i]
for the pressure and kinetic energy tensors.
When this fix is applied, any orientation-dependent vector or tensor-valued quantities computed, except for the tensors from compute pressure/uef
/compute temp/uef
and coordinates from dump cfg/uef
, will not be in the same coordinate system as the flow field. See the implementation details for further information.
This fix can be used with write_restart
and read_restart
, run_style respa
, and fix modify
, however custom pressure and temperature computes must be of type pressure/uef
and temp/uef
. When resuming from restart files, you may need to use box tilt large
since LAMMPS does not always agree that the simulation box is fully reduced.
fix ID all npt/uef temp Tstart Tend Tdamp erate eps_x eps_y keyword value...
- ID = name for the fix
- Tstart,Tstop = external temperature at start/end of run
- Tdamp = temperature damping parameter (time units)
- Pstart,Pstop = external pressure at start/end of run
- Pdamp = pressure damping parameter (time units)
- eps_x = strain rate in x dimension 1/(time units)
- eps_y = strain rate in y dimension 1/(time units)
Additional keywords: - x or y or z = Pstart Pstop Pdamp
- iso = Pstart Pstop Pdamp
- strain = initial level of strain (default=0). Use of this keyword is not recommended, but may be necessary when resuming a run with data file. This keyword is not needed when restart files are used.
- ext = x or y or z or xy or xz or yz or xyz (default=xyz). These are "external" dimensions used in pressure control. For example, for uniaxial extension in the z direction, x and y correspond to free surfaces. The setting xy will only control (P_xx+P_yy)/2 to the target external pressure.
- The following additional keywords from
fix nvt
can be used with this fix: couple, tchain, pchain, tloop, ploop, drag
- Uniaxial flow
fix f1 all npt/uef temp 400 400 300 iso 1 1 3000 erate 0.00001 -0.000005 ext yz
- Biaxial flow
fix f2 all npt/uef temp 400 400 300 z 1 1 3000 erate 0.000005 0.000005
The usage notes for fix nvt/uef
apply to fix npt/uef
as well. There are two ways to control the pressure using this fix. The first method involves using theext
keyword along with the iso
pressure style. With this method, the pressure is controlled by scaling the simulation box isotropically to achieve the average stress in the directions specified by ext
.
For example, this command will control the hydrostatic pressure under uniaxial tension:
fix f1 all npt/uef temp 0.7 0.7 0.5 iso 1 1 5 erate -0.5 -0.5 ext xyz
This command will control the average stress in compression directions under uniaxial tension:
fix f1 all npt/uef temp 0.7 0.7 0.5 iso 1 1 5 erate -0.5 -0.5 ext xy
The second method involves setting the normal stresses using the x
, y
, and/or z
keywords. When using this method, the same pressure must be specified via Pstart
and Pstop
for all dimensions controlled. Any choice of pressure conditions that would cause LAMMPS to compute a deviatoric stress are not permissible and will result in an error. Additionally, all dimensions with controlled stress must have the same applied strain rate. The ext
keyword must be set to the default value (xyz
) when using this method.
For example, the following commands will work:
fix f1 all npt/uef temp 0.7 0.7 0.5 x 1 1 5 y 1 1 5 erate -0.5 -0.5
fix f1 all npt/uef temp 0.7 0.7 0.5 z 1 1 5 erate 0.5 0.5
The following commands will not work:
fix f1 all npt/uef temp 0.7 0.7 0.5 x 1 1 5 z 1 1 5 erate -0.5 -0.5
fix f1 all npt/uef temp 0.7 0.7 0.5 x 1 1 5 z 2 2 5 erate 0.5 0.5
compute ID all temp/uef
- ID = name for the compute
compute c1 all temp/uef
This compute requires a fix nvt/uef
or fix npt/uef
. It computes the kinetic energy tensor in the reference frame of the flow field.
See the documentation for compute temp
for further details on output.
compute ID all pressure/uef temp-ID
- ID = name for the compute
- temp-ID = ID of compute that calculates temperature
Additional keywords: - The following additional keywords from
compute pressure
may be used with this fix: ke or pair or bond or angle or dihedral or improper or kspace or fix or virial
compute c1 all pressure/uef c_1_temp
This compute requires afix nvt/uef
or fix npt/uef
. It computes the pressure tensor in the reference frame of the flow field.
The pressure tensor computed from compute pressure/uef
is only accurate if its temperature compute, specified by temp-ID
, is a compute temp/uef
.
See the documentation for compute pressure
for further details on output.
dump ID all cfg/uef N file mass type xs ys zs keyword value
- N = dump every this many timesteps
- file = name of file to write dump info to
Additional keywords: - See the documentation for
dump cfg
for additional keywords.
dump d1 all cfg/uef 100 dump.*.cfg mass type xs ys zs
This command requires a fix nvt/uef
or fix npt/uef
. It outputs the atomic positions in the reference frame of the flow field. Only the positions are in the proper reference frame; if the atomic velocities are specified as an output, for example, they will not be in the flow field reference frame.
See the dump cfg
documentation for further information on writing trajectories with cfg files.
The simulation box used in the boundary conditions developed by Hunt and Dobson does not have a consistent alignment relative to the applied flow field. This can be seen in the video of the simulation box under uniaxial extension at the top of the page. LAMMPS utilizes an upper-triangular simulation box, making it impossible to express the evolving simulation box in the same coordinate system as the flow field. The UEF package keeps track of two coordinate systems: the flow frame, and the upper triangular LAMMPS frame. The coordinate systems are related to each other through the QR decomposition, as is illustrated in the image below.
During most molecular dynamics operations, the system is represented in the LAMMPS frame. Only when the positions and velocities are updated is the system rotated to the flow frame, and it is rotated back to the LAMMPS frame immediately afterwards. For this reason, all vector-valued quantities (except for the tensors from compute pressure/uef
and compute temp/uef
) will be computed in the LAMMPS frame. Rotationally invariant scalar quantities like the temperature and hydrostatic pressure, on the other hand, will be computed correctly. Additionally, the system is in the LAMMPS frame during all of the output steps, and therefore trajectory files made using the dump
command will be in the LAMMPS frame unless the dump cfg/uef
command is used.
Sample input scripts and data files for uniaxial and biaxial extension of a WCA fluid can be found in the examples/
subdirectory. Please see examples/README
for further information.
The methods described with inherit error/warning messages from fix npt/nvt
, compute pressure
, and dump cfg
. Additional error messages associated with this package are below:
- "Illegal fix nvt/npt/uef command" - The fix was not called correctly. Check the usage notes.
- "Keyword erate must be set for fix nvt/npt/uef command" - Self-explanatory.
- "Simulation box must be triclinic for fix/nvt/npt/uef" - Self-explanatory.
- "Only normal stresses can be controlled with fix/nvt/npt/uef" - The keywords xy xz and yz cannot be used for pressure control.
- "The ext keyword may only be used with iso pressure control" - Self-explanatory
- "All controlled stresses must have the same value in fix/nvt/npt/uef" - Stress control is only possible when the stress specified for each dimension is the same
- "Dimensions with controlled stresses must have same strain rate in fix/nvt/npt/uef" - Stress-controlled dimensions with the same strain rate must have the same target stress.
- "Can't use another fix which changes box shape with fix/nvt/npt/uef" - The
fix npt/nvt/uef
command must have full control over the box shape. You cannot use a simultaneousfix deform
command, for example. - "Pressure ID for fix/nvt/uef doesn't exist" - The
compute pressure
introduced viafix_modify
does not exist - "Using fix nvt/npt/uef without a compute pressure/uef" - The
compute pressure
introduced viafix_modify
is not acompute pressure/uef
. - "Initial box is not close enough to the expected uef box" - The initial box does not correspond to the shape required by the value of the strain keyword. If the default strain value of zero was used, the initial box is not cubic.
- "Can't use compute pressure/uef without defining a fix nvt/npt/uef" - Self-explanatory.
- "Can't use compute pressure/uef without defining a fix nvt/npt/uef" - Self-explanatory.
- "Can't use compute/temp without defining a fix nvt/npt/uef" - Self-explanatory.
- "Temperature control must be used with fix nvt/uef" - Self-explanatory.
- "Pressure control can't be used with fix nvt/uef" - Self-explanatory.
- "Temperature control must be used with fix npt/uef" - Self-explanatory.
- "Pressure control must be used with fix npt/uef" - Self-explanatory.
Users of this package are encouraged to cite the following article in scientific publications:
D.A. Nicholson, G.C. Rutledge, "Molecular simulation of flow-enhanced nucleation in n-eicosane melts under steady shear and uniaxial extension", J. Chem Phys., 2016, 145 (24), http://aip.scitation.org/doi/full/10.1063/1.4972894.