2.22. MedeA LAMMPS: A Powerful Gateway to a Powerful Molecular Dynamics Program
Contents
| download: | pdf |
|---|
2.22.1. Introduction
MedeA LAMMPS provides flexible calculation setup and analysis capabilities to unlock the power of the LAMMPS molecular dynamics (MD) code.
MedeA LAMMPS automates the details of properly formatting molecules, fluids, and solids into the required LAMMPS coordinate, connectivity, forcefield parameter, and command-line formats. It also provides access to the core capabilities of LAMMPS, including minimization, molecular dynamics simulations using the NVE, NVT, and NPT ensembles, and materials properties calculations.
MedeA LAMMPS fully integrates with MedeA Forcefield for advanced forcefield handling and assignment, and any compatible custom forcefield can be used. There are also options for expert LAMMPS users to add any LAMMPS commands to existing protocols, or to prepare completely customized simulations.
After each calculation, MedeA LAMMPS automatically analyzes the block averages and fluctuations of temperature, pressure, density, cell parameters, total energy and all energy components (potential, kinetic, Coulomb and van der Waals), and stress tensor elements, reporting convergence statistics and uncertainties computed according to a method applicable to properties sampled using molecular dynamics or Monte Carlo methods. Additionally, visualization of trajectories of MD simulations and structure optimizations can be performed.
The user interface to LAMMPS is based on flowcharts and those from any previous MedeA LAMMPS calculation can be reused, edited, shared with colleagues, and rerun, even on different systems and compute servers.
Bring up the Flowchart interface by selecting Jobs >> New Job…, and you should see the following:
2.22.2. Variables
A MedeA LAMMPS flowchart usually starts with a Variables stage:
With simulation variables defined in the following way:
If no simulation variables are defined, the LAMMPS stages use the following default variables:
- Time step (tstep): 1 fs
- Temperature (T): 300 K
- Pressure (P): 1 atm
Note
A time step size of 1 fs works well for covalent forcefields with pre-defined bonds, angles, and dihedrals. For reactive potentials without pre-defined bonds, such as Tersoff, ReaxFF, and SNAP, time step sizes from 0.1 fs to 0.4 fs are recommended.
2.22.3. Initialization
Upon opening the LAMMPS stage, you will see the following interface:
The Initialization section in the right-side panel has one available stage: Initialize LAMMPS. The LAMMPS stage must start with the Initialize stage:
The parameters are:
Periodicity: Choose from 3-d periodic for bulk systems and layer perpendicular to Z for slab models.
Hint
Even with 3-D periodic boundary condition, a slab model can effectively be created by having large enough vacuum size (> ~20 \({\mathring{\mathrm{A}}}\))
Nonbond method: For long-range Coulombics, choose from Cutoff for a simple spherical cutoff, PPPM for the particle–particle–particle–mesh long-range electrostatic solver, and Ewald for conventional Ewald. This is a forcefield-sensitive option, meaning this option does not appear for all forcefields.
Nonbond cutoff: The cutoff value for spherical or real-space cutoff. The default value of 9.5 is usually good. This is a forcefield-sensitive option.
Long range precision: Only used with PPPM and Ewald methods. Determines the reciprocal-space convergence. The default value of 0.00001 is usually a good value. This is a forcefield-sensitive option.
Add tail corrections for Van der Waals interactions: The
default is yes. This is a forcefield-sensitive option.Run LAMMPS on: Choose from CPU, OpenMP, and GPU.
- CPU: Simulations will run on the CPU cores and parallelize over MPI ranks.
- OpenMP: Simulations will run on the CPU cores and parallelize over OpenMP threads.
- GPU: Simulations will run on the combination of CPU cores and GPU cards.
Hint
This option is only visible if the
Enable running Job on GPU checkbox from File >>
Preferences… >> Miscellaneous is
checked.Note
MedeA LAMMPS supports NVIDIA GPU cards with compute capabilities (https://developer.nvidia.com/cuda-gpus) from 3.0 to 8.6. This includes the Kepler, Maxwell, Pascal, Volta, Turing, and Ampere series GPU cards. For HPC cluster compute nodes the Tesla Kepler, Pascal, Volta, and Ampere series are recommended. For workstations and desktops, GeForce Titan V, and Quadro GP100 and GV100 are recommended. Both Linux and Windows are supported.
2.22.4. Bias
The Bias section has one available stage: Orientation:
The parameters are:
- Subset: Add bias to this atom pair subset.
- Force: the magnitude of the added bias.
- x, y, and z: the direction of the added bias.
2.22.5. Single Point
The Single Point section has two available stages: Single Point Energy and Single Point Forces. Both of these stages have no adjustable parameters.
2.22.6. Minimization
The Minimization section has one available stage: Minimize.
The Minimize stage has the following parameters:
Relax cell: Choose from No, Isotropically, With x=y, With y=z, With x=z, and Anisotropically. Depending on the choice, there are 0–6 cell parameters to fix (remaining unchanged during minimization), 0–6 stress values to relax each cell parameter to, and maximum volume change:
Method: Choose from Conjugate Gradient (CG), Steepest Descent (SD), and Hessian-free truncated Newton (HFTN). The default is CG.
The CG method is the Polak–Ribiere (PR) version. At each iteration, the force gradient is combined with the previous iteration information to compute a new search direction perpendicular (conjugate) to the previous search direction. The PR variant affects how the direction is chosen and how the CG method is restarted when it ceases to make progress. The PR variant is thought to be the most effective CG choice for most problems.
With the SD method, at each iteration, the search direction is set to the downhill direction corresponding to the force vector (negative gradient of energy). Typically, steepest descent will not converge as quickly as CG but may be more robust in some situations.
With the HFTN algorithm, at each iteration, a quadratic model of the energy potential is solved by a conjugate gradient inner iteration. The Hessian (second derivatives) of the energy is not formed directly but approximated in each conjugate search direction by a finite difference directional derivative. When close to an energy minimum, the algorithm behaves like a Newton method and exhibits a quadratic convergence rate to high accuracy. In most cases, the behavior of HFTN is similar to CG, but it offers an alternative if CG seems to perform poorly.
Warning
When changing cell parameters, only the conjugate gradient and steepest descent methods can be used.
Linesearch: Choose from Fast, Quadratic, or Backtrack. The default is Fast.
Maximum step: Maximum distance (\({\mathring{\mathrm{A}}}\)) to move atoms. The default value of 0.05 is usually good.
Energy tolerance: The default is 0.0 eV or kcal/mol.
Force tolerance: The default is 1.0 eV/\({\mathring{\mathrm{A}}}\) for metallic/NIST and COMB3 forcefields and 1.0 kcal/mol-\({\mathring{\mathrm{A}}}\) for all other forcefields. For numerical stability, the force tolerance value is limited to a minimum of 1e-12 eV/\({\mathring{\mathrm{A}}}\) or kcal/mol-\({\mathring{\mathrm{A}}}\). If a user specifies a smaller value than 1e-12, this value is overwritten with 1e-12 when generating the input files for LAMMPS.
Maximum iterations: Maximum iterations for energy evaluations. The default is 1000.
Maximum force evaluations: Maximum iterations for force evaluations. The default is 10000.
- Write trajectory and the default value for
Every (steps): is 100.
2.22.7. Building and Editing
This section has two available stages: Set Cell and Compress Layer.
2.22.7.1. Set Cell
- Set Cell: Changes the cell dimension by one of the following:
- using density and its Density
- using volume and its Volume
- scaling by a factor and its Factor
All of the above options can be a value or a variable.
2.22.7.2. Compress Layer
The Compress Layer stage performs combined equilibration and compression of layered materials in the simulation cell, ensuring that none of the atoms in the system are outside the desired layer. First, a minimization is performed, then a compression with a set of “indenters” under an NVT ensemble, followed by one last minimization.
The Compress Layer stage has some similar options to the Minimize stage with the following exceptions:
Layer thickness: Set the desired layer thickness with options of Explicit, From target density, or Use initial c dimension:
- Explicit: Enter the desired layer thickness explicitly and the Compress Layer stage adjusts the c dimension accordingly.
- From target density: Enter the desired layer density (excluding vacuum area) and the Compress Layer stage adjusts the c dimension accordingly.
- Use initial c dimension: The final c dimension does not change.
Time: Simulation time.
Time step: Time step size for NVT integration.
Initial Temperature and Final Temperature: Initial and final temperature for NVT integration.
Warning
A Compress Layer stage does not apply periodic boundaries in all 3 directions and must be used with layer perpendicular to Z periodicity and the Cutoff Nonbond method. It should therefore only be used by itself or with other Compress Layer stages. Results may be unpredictable if combined with other stages, such as those found in the Minimize or Dynamics sections that are normally used with 3-dimensional systems.
Warning
The Compress Layer stage may not work correctly when there are frozen atoms in the system. We highly recommend to un-freeze all of the atoms before utilizing a Compress Layer stage.
2.22.8. Dynamics
The Dynamics section has ten stages and the following stages are included in the standard MedeA LAMMPS license:
- Initialize Velocities
- NVE Ensemble
- NVT Ensemble
- NPT Ensemble
While the following stages (modules) require additional licenses:
- Cohesive Energy Density
- Thermal Conductivity
- Viscosity
- Diffusion
- Surface Tension
- Deposition
2.22.8.1. Initialize Velocities
The Initialize Velocities stage thermalizes the cell to a user-defined temperature:
- Initial Temperature: Thermalize the cell to this temperature.
- No net translation: Whether to ensure no net
translational velocity is present in the system.
- No net rotation: Whether to ensure no net
rotational velocity is present in the system.
- Random Seed (positive integer): A random seed for the random number generator. The seed can be assigned manually.
- Assign a random seed automatically: randomly generates
a random seed.
2.22.8.2. NVE Ensemble
The NVE Ensemble stage performs time integration without any alterations to the equations of motion.
The parameters in the Control tab are:
- Time: Duration of the simulation run (defaults to 100 ps).
- Time Step: Time step size employed in solving the equations of motion.
- Sampling: Number of samples employed in performing averaging. This parameter does not affect dynamics.
- Trajectory: Number of trajectory frames saved during the molecular dynamics simulation. This parameter does not affect dynamics.
The Analysis tab allows you to add analyses to the dynamics stage. See The Analysis Tab for more details.
Warning
Using an NVE stage does not automatically guarantee an NVE ensemble. You should examine Job.out and the energy profile to ensure that the total energy is conserved so that an NVE ensemble is achieved.
2.22.8.3. NVT Ensemble
The NVT Ensemble stage performs time integration with a thermostat added to the equations of motion.
The parameters, in addition to those in the NVE Ensemble stage, are:
Initial Temperature and Final Temperature: These two parameters can be the same or different (for establishing cooling or heating). These parameters can be values or variables.
Control: Choose a thermostat algorithm from one of the following:
rescaling: interval (steps), window, amount of rescaling
Langevin: Damping (fs), Random Seed (integer)
Berendsen: Damping (fs)
Nose-Hoover: Damping (fs) and Drag
Hint
The default options of Nose-Hoover and Langevin are the good choices for a thermostat. The Damping parameter is recommended to be 100 times the time step size. For example, for a 1 fs time step, Damping is recommended to be
100while for a 0.2 fs time step, the recommended value is20.
Warning
Using an NVT stage does not automatically guarantee an NVT ensemble. You should examine Job.out and the temperature profile to ensure that temperature is maintained at the defined value so that an NVT ensemble is achieved.
2.22.8.4. NPT Ensemble
The NPT Ensemble stage performs time integration with both a thermostat and a barostat added to the equations of motion.
The parameters, in addition to those in the NVT stage are:
Initial Pressure and Final Pressure: These two parameters can be the same or different (for pressurization or depressurization). These parameters can be values or variables.
Restrict cell motion: Controls how the cell volume and shape are equilibrated and relaxed:
isotropic: Only a, b, and c cell parameters (x, y, and z dimensions) are relaxed and the three components relax to one averaged value.
fixed angles: Also known as anisotropic with a, b, and c cell parameters relaxed independently.
constrained: Allows users to choose to relax any of the six cell parameters (a, b, c, alpha, beta, and gamma) independently.
unconstrained: All six cell parameters relax independently.
Warning
To relax the cell with the constrained or unconstrained options, the cell must be non-orthorhombic or you need to check the box Allow orthorhombic cell angles to relax in the Initialize stage.
Control: Choose a combination of thermostat and barostat from the following options:
Rescaling/Berendsen: Interval (steps), Window, Amount of rescaling, Pressure Damping (fs), Bulk Modulus
Langevin/Berendsen: Damping (fs), Random Seed (integer), Pressure Damping (fs), Bulk Modulus
Berendsen/Berendsen: Damping (fs), Pressure Damping (fs), Bulk Modulus
Nose-Hoover T/Berendsen: Damping (fs), Drag (integer), Pressure Damping (fs), Bulk Modulus
Nose-Hoover T & P: Temperature Damping (fs), Pressure Damping (fs) and Drag (0 to 2’ish)
Hint
The default option of Nose-Hoover T & P is the recommended combination of thermostat and barostat for most systems.
Warning
Using an NPT stage does not automatically guarantee an NPT ensemble. You should examine Job.out and the temperature, stresses, and pressure profiles to ensure that temperature, stresses, and pressure are maintained at the defined value so that an NPT ensemble is achieved.
2.22.8.5. Cohesive Energy Density
The Cohesive Energy Density stage is a separate module available from Materials Design. It performs time integration under an NVT ensemble while making automated extractions and calculations to compute the cohesive energy density.
The parameters, in addition to those in the NVT stage are:
- CED Sampling: Calculate cohesive energy density every this many samples, steps, or amount of time.
- CED Cutoff: The cutoff value to calculate the energy.
Results are written in Job.out.
2.22.8.6. Thermal Conductivity
The Thermal Conducitivity stage is a separate module available from Materials Design and described in section Thermal Conductivity.
2.22.8.7. Viscosity
The Viscosity stage is a separate module available from Materials Design and described in section Viscosity.
2.22.8.8. Diffusion
The Diffusion stage is a separate module available from Materials Design and described in section Diffusion.
2.22.8.9. Surface Tension
The Surface Tension stage is a separate module available from Materials Design. It performs time integration under an NVT ensemble while calculating the surface tension. It has the same parameters as the NVT stage.
2.22.8.10. Deposition
The Deposition stage is a separate module available from Materials Design and described in section MedeA Deposition.
2.22.9. Custom
The Custom Code stage is available with the standard MedeA LAMMPS license. It enables the addition of any custom LAMMPS commands not accessible from the above stages. Variables defined in any previous stages are passed into the Custom Code stage.
Note
Characters in MedeA LAMMPS Custom Code stages with special meaning in Tcl (such as square brackets) will be interpreted in the usual Tcl manner (which can be convenient as this enables access to Flowchart variables in your LAMMPS custom stage). If you wish to avoid such interpretation, you can escape such characters with a preceeding backslash.
2.22.10. The Analysis Tab
The Analysis tab is available in the following stages: NVE, NVT, NPT, Surface Tension, and Deposition. You can add one or more of the following analyses to these stages:
- Distances: The Distance Analysis analyzes the interatomic distance between pairs of atoms defined by a pair subset. Therefore, the subset must be defined for the structure and entered in the Pair subset: box. The pair distances are analyzed every Frequency and plotted using the defined Number of histogram bins. Please see the tutorial Distance Analysis with LAMMPS for more details.
- Distributions: The Distribution Analysis analyzes the spatial distribution of the given Subset along the defined Direction, using spatial bins of the specified width. The Profile can be Number or Density/mass. The Number value means the number is computed for each chunk, i.e., number/chunk. The Density/mass value means the mass density is computed for each chunk, i.e. total-mass/volume/chunk. The Number of intervals defines how many analyses are performed throughout the simulation. Please see the tutorial Deposition of O2 on a Si Surface with Reactive Potentials for more details.
- PairCorrelations: The Pair Correlation Analysis plots the time-averaged pair correlation function (also known as the radial distribution function) between Subset 1 and Subset 2. These subsets must be created for the structure prior to setting up the flowchart, and these subsets must contain only one atom type each. The Number of intervals defines how many analyses are performed throughout the simulation using the defined Number of histogram bins.
- GroupInteractions: The Group Interaction Analysis computes the non-bond interactions between Subset 1 and Subset 2 and the interactions are time-averaged every Frequency. You can choose to include the pair interactions (van der Waals + Coulomb) and kspace contributions (long-range summation) or not. Additional keywords such as boundary and molecule can be added in the Additional Keywords box. See compute group/group for detailed usage of these keywords.
Note
GroupInteractions analysis only supports valence force fields (such as PCFF+).
2.22.11. Example
This flowchart illustrates the efficient use of MedeA LAMMPS flowcharts:
The Variables stage defines the overall parameters for temperature, pressure, and basic time step.
A larger supercell is constructed from the provided molecule.
Two different LAMMPS stages can be used. For instance, an equilibration stage that makes use of a computationally less demanding approach for non-bond interactions (e.g. Cutoff) can be used before a simulation with more computationally intense settings.
| download: | pdf |
|---|