Building a Simulation
This module is designed to help GPUMD users create a valid simulation for the gpumd executable. The
Simulation class is the main class that facilitates this setup. It will create a run.in file
with a structure that follows the guidlines defined in the gpumd documentation.
Setup the potentials
Minimize system (if needed)
Setup static calculations (cohesive, elastic, phonon)
Run molecular dynamics simulation (if desired)
The Simulation class will also copy potential files to the correct working directory, create the
appropriate xyz.in file, and ensure that the files are consistent with the rules for GPUMD.
At a high level, to build a GPUMD simulation using gpyumd, users first need to instantiate a
Simulation object and then add potentials, static calculations, and runs (with any supported
keywords) to that Simulation object. To add potentials, users must pass
Potential objects to the add_potential() method. To add static
calculations, users must pass Keyword objects for static calculations as defined in the
keywords section of the GPUMD documentation. As of this writing, these are the Minimize,
ComputeCohesive, ComputeElastic, and
ComputePhonon keywords. For the molecular dynamics part of a simulation, users create a series
of Run objects using the add_run() method. Each
Run can be constructed using a set of keywords (non-static calculations) found in the
keyword module.
We will run through two examples here better illustrate how to build a Simulation using gpyumd.
The first example is an equlibrium MD simulation that calculates the thermal conductivity of graphene and demonstrates
how to construct a Run. The second is a phonon calculation of Si, which demonstrates how to add
static calculations.
Example: thermal conductivity of graphene
In this example, we will setup a simulation for an EMD thermal conductivity calculation for graphene using the
Simulation class. First, we will setup the graphene structure (gr) using the
GpumdAtoms class.
# import statements
from ase.build import graphene_nanoribbon
from gpyumd.atoms import GpumdAtoms
from gpyumd.sim import Simulation
import gpyumd.keyword as kwd
# structure creation
gr = GpumdAtoms(graphene_nanoribbon(60, 36, type='armchair', sheet=True, vacuum=3.35/2, C_C=1.44))
gr.euler_rotate(theta=90)
lx, lz, ly = gr.cell.lengths()
gr.cell = gr.cell.new((lx, ly, lz))
gr.center()
gr.pbc = [True, True, False]
gr.set_cutoff(2.1)
gr.set_max_neighbors(3)
With a structure created and finalized, we can begin creating a simulation with:
emd_sim = Simulation(gnr, driver_directory='.')
The driver_directory parameter determines where driver file will be located. This can be important for the path
definitions of potentials, so be sure that it is correct.
Now that we have a Simulation object defined, we want to add a Run. A
Run is simply a group of keywords that are coupled to a specific run keyword in the run.in
file. The first run we want to create is an equilibration and we can do this using the following:
header_comment = "Equilibration" # A comment that will go above the keywords in the run
curr_run = emd_sim.add_run(number_of_steps=1e6, run_name='equilibration', run_header=header_comment)
Note that all of the parameters for the add_run method are optional; however, if you do
not provide the number_of_steps parameter, you will have to add the RunKeyword keyword to
the Run object.
In the equilibration, we want to allow the structure to relax to the appropriate size for the temperature. To do this,
we need to define an integrator suitable for the purpose using the ensemble keyword. In gpyumd, we define the ensemble
keyword using the EnsembleNVE, EnsembleNVT,
EnsembleNPT, and EnsembleHeat classes. In this example, we want
constant atom number, pressure, and temperature, i.e., an NPT ensemble. To handle the extra
complexity for NPT ensembles, we recommend using the class functions for the EnsembleNPT keyword. These class functions
are isotropic(), orthogonal(), and
triclinic(). Below, we initialize the NPT ensemble using the
orthogonal() class method. The keywords for the simulation are then:
pres = 0 # GPa
elast = 53.4059 # GPa
keywords = [
kwd.Velocity(initial_temperature=300),
kwd.EnsembleNPT.orthogonal(method='npt_ber', initial_temperature=300,
final_temperature=300, thermostat_coupling=100,
barostat_coupling=2000, p_xx=pres, p_yy=pres, p_zz=pres,
c_xx=elast, c_yy=elast, c_zz=elast),
kwd.TimeStep(dt_in_fs=1),
kwd.NeighborOff(),
kwd.DumpThermo(1000)
]
for keyword in keywords:
curr_run.add_keyword(keyword)
If you are working in an interactive environment like Jupyter, you can take a quick look at your run by simply running:
curr_run in a cell. Here, this would output:
Run: 'equilibration'
# Equilibration
Velocity(initial_temperature=300)
EnsembleNPT(condition=orthogonal, initial_temperature=300, final_temperature=300, thermostat_coupling=100,
barostat_coupling=2000, p_xx=0, p_yy=0, p_zz=0, c_xx=53.4059, c_yy=53.4059, c_zz=53.4059)
TimeStep(dt_in_fs=1, max_distance_per_step=None)
NeighborOff()
DumpThermo(interval=1000)
run 1000000 steps, dt_in_fs=1 -> 1.0 ns
In the next part of the simulation, the production run, we want to shift to an NVE ensemble and compute the heat current autocorrelation. This can be done with the following:
# Production run
header_comment = "Production"
curr_run = emd_sim.add_run(run_name='production', run_header=header_comment)
# Timestep is propagated from last run
keywords = [
kwd.EnsembleNVE(),
kwd.NeighborOff(),
kwd.ComputeHAC(sample_interval=20, num_corr_steps=50000, output_interval=10),
kwd.RunKeyword(number_of_steps=1e7) # Can add run here instead of during add_run fucntion
]
for keyword in keywords:
curr_run.add_keyword(keyword)
Finally, we need to add the potential to describe our Graphene. Note that we can add this at any point in the
simulation setup. To add the potential, we need to create a Potential keyword object and pass
it to the simulation through the add_potential() method like in the following code:
potential_directory = "path/to/GPUMD/potentials/tersoff"
tersoff_potential = \
kwd.Potential(filename='Graphene_Lindsay_2010_modified.txt', symbols=['C'], directory=potential_directory)
emd_sim.add_potential(tersoff_potential)
The simulation definition is now complete and can be created using the following:
emd_sim.create_simulation(copy_potentials=True)
Note that the copy_potentials parameter determines whether or not you want to copy the potential file to your
simulation directory or if you want to use relative paths to the potential files in the run.in file. After the
create_simulation() method call, the simulation object generates the run.in file and the
xyz.in file and moves the potential to the simulation directory. With the keywords and parameters we used here, the
run.in file’s contents look like:
potential Graphene_Lindsay_2010_modified.txt 0
# Equilibration
velocity 300
time_step 1
ensemble npt_ber 300 300 100 0 0 0 53.4059 53.4059 53.4059 2000
neighbor off
dump_thermo 1000
run 1000000
# Production
ensemble nve
neighbor off
compute_hac 20 50000 10
run 10000000
Example: phonon dispersion of silicon
In this example, we will setup a phonon dispersion calculation of silicon (Si) using the Simulation
class. This example emphasizes the inclusion of static calculations when building a Simulation.
First, we will set up the Si structure. For more information about the structure generation, see The GpumdAtoms Object.
# import statements
from ase.lattice.cubic import Diamond
from ase.build import bulk
from gpyumd.atoms import GpumdAtoms
from gpyumd.sim import Simulation
import gpyumd.keyword as kwd
# Create unit cell
a=5.434
Si_UC = GpumdAtoms(bulk('Si', 'diamond', a=a))
Si_UC.add_basis()
# Create 8 atom diamond structure
Si = Si_UC.repeat([2,2,1])
Si.set_cell([a, a, a])
Si.wrap()
# Complete full supercell
Si = Si.repeat([2,2,2])
Si.set_max_neighbors(4)
Si.set_cutoff(3)
Next, we need to create the relevant extra files needed for the phonon dispersion calculation: kpoints.in and
basis.in. (In the future, these files may also be handled by the Simulation class.)
Si.write_basis()
Si_UC.write_kpoints(path='GXKGL',npoints=400)
We can now create the Simulation object that we want to work with:
phonon_sim = Simulation(Si, driver_directory='.')
We want to calculate the phonon dispersion, so we add the ComputePhonon keyword to the
phonon_sim object using the add_static_calc() method.
phonon_sim.add_static_calc(kwd.ComputePhonon(cutoff=5, displacement=0.005))
Next, we need to define the potential that describes the Si interactions:
potential_directory = "/path/to/GPUMD/potentials/tersoff"
tersoff_potential = \
kwd.Potential(filename='Si_Fan_2019.txt', symbols=['Si'], directory=potential_directory)
phonon_sim.add_potential(tersoff_potential)
Finally, we can generate the run.in file and the xyz.in file and move the potential file to the simulation directory using:
phonon_sim.create_simulation(copy_potentials=True)
The resulting run.in file has the contents:
potential Si_Fan_2019.txt 0
compute_phonon 5 0.005
List of all methods
- class gpyumd.sim.Potentials(gpumd_atoms: GpumdAtoms)[source]
- add_potential(potential: Potential) None[source]
Adds a potential keyword to the current list of potentials.
- Parameters
potential – A potential keyword object
- copy_potentials(directory: str) None[source]
Copies all of the potentials to the selected directory.
- Parameters
directory – Directory to copy potentials to
- class gpyumd.sim.Run(gpumd_atoms: GpumdAtoms, number_of_steps: Optional[int] = None, run_name: Optional[str] = None, run_header: Optional[str] = None)[source]
- Parameters
gpumd_atoms – Atoms for the simulation
number_of_steps – Number of steps to be run in the Run.
run_name – Name of the run
run_header – A comment that will be added to the line above all keywords for this run
- add_keyword(keyword: Keyword, final_check: bool = False) None[source]
- Adds a keyword object to the run. Verifies that the keyword is valid
(to the extent that it can be initially).
- Parameters
keyword – The keyword to add to the run.
final_check – Use only if you know what you’re doing. It is normally used to validate the run when finalized.
- get_output() List[str][source]
Gets the run.in, textual representation of the current run.
- Returns
The Run object in text format of run.in file
- class gpyumd.sim.Simulation(gpumd_atoms: GpumdAtoms, run_directory: Optional[str] = None, driver_directory: Optional[str] = None)[source]
Stores the relevant information for a full ‘gpumd’ executable simulation.
- Parameters
gpumd_atoms – Final structure to be used with the GPUMD simulation.
run_directory – Directory of the simulation.
driver_directory – Directory where the driver input file will be found. Defaults to run_directory.
- add_potential(potential: Potential) None[source]
Adds a potential keyword to the simulation object.
- Parameters
potential – A potential keyword object
- add_run(number_of_steps: Optional[int] = None, run_name: Optional[str] = None, run_header: Optional[str] = None) Run[source]
Adds a new run to a simulation.
- Parameters
number_of_steps – number of steps to run
run_name – A name for the run. Note that runs can share the same name.
run_header – A comment that will be added to the line above all keywords for this run
- Returns
A reference to the new run in the simulation.
- add_static_calc(keyword: Keyword) None[source]
Adds a static calculation to the current simulation object. Currently, only the ‘compute_cohesive’, ‘compute_elastic’, ‘compute_phonon’, and ‘minimize’ keywords are static.
- Parameters
keyword – A keyword object for a static calculation
- create_simulation(copy_potentials: bool = True, use_velocity: bool = False) None[source]
Generates the required files for the gpumd simulation
- Parameters
copy_potentials – Whether or not to copy potentials to the simulation directory. If False, relative path is provided.
use_velocity – Whether or not to add velocities to the xyz.in file