MolecularDynamics

MolecularDynamics(configuration, constraints=None, trajectory_filename=None, steps=None, log_interval=None, method=None, xyz_filename=None, pre_step_hook=None, post_step_hook=None, write_velocities=True, write_forces=True, write_stresses=None, domain_decomposition_pattern=<class 'NL.ComputerScienceUtilities.NLFlag._NLFlag.Automatic'>)

Function for performing a molecular dynamics simulation.

Parameters:
  • configuration (MoleculeConfiguration | BulkConfiguration | DeviceConfiguration | SurfaceConfiguration | MDTrajectory) – The initial configuration or a trajectory. If a trajectory is given the last configuration will be used as the initial configuration for this MD run and the trajectory will be appended to.
  • constraints (list of ints | list of BaseConstraint) – The list of atomic indices, denoting fixed atoms, or constraint objects, such as RigidBody.
    Default: []
  • trajectory_filename (str | None) – The filename of the file to be used for storing the trajectory, or None if no trajectory should be written. A trajectory filename should not be given if configuration is a trajectory.
    Default: None
  • steps (int) – The number of time-steps to take in simulation.
    Default: 50
  • log_interval (int) – The resolution used in saving steps to a trajectory file ad writing the log file, where a value of 1 results in all steps being saved, and e.g. a value of 2 resulting in every other step being discarded. If configuration is a trajectory (i.e. this is a restart MD calculation) this parameter will be the same value as it was in the previous trajectory.
    Default: 1
  • method (BaseMDmethod) – The molecular dynamics method used for the simulation.
    Default: NVEVelocityVerlet
  • xyz_filename (str) – The name of the file to be used for storing the xyz trajectory, or None if no xyz-trajectory should be written.
    Default: None
  • pre_step_hook (function | list of functions | None) – An optional user-defined function or a list of functions which will be called just before the forces evaluation. The signature of the function requires the arguments (step, time, configuration, forces, stress). The return status is ignored. Unhandled exceptions will abort the Molecular Dynamics evaluation. If a list is given the functions will be called in the given order.
    Default: None
  • post_step_hook (function | list of functions | None) – An optional user-defined function or a list of functions which will be called just after the forces evaluation. The signature of the function requires the arguments (step, time, configuration, forces, stress). The return status is ignored. Unhandled exceptions will abort the Molecular Dynamics evaluation. If a list is given the functions will be called in the given order.
    Default: None
  • write_velocities (bool) –

    Write the velocities to the trajectory file every log_interval steps.

    If configuration is a trajectory (i.e. this is a restart MD calculation) this parameter will be the same value as it was in the previous trajectory.
    Default: True

  • write_forces (bool) –

    Write the forces to the trajectory file every log_interval steps.

    If configuration is a trajectory (i.e. this is a restart MD calculation) this parameter will be the same value as it was in the previous trajectory.
    Default: True

  • write_stresses (bool) –

    Write the stress to the trajectory file every log_interval steps. A value of None means that write_stresses is True by default for NPT methods and False otherwise. This is to avoid the additional work of calculating the stress when it is not needed.

    If configuration is a trajectory (i.e. this is a restart MD calculation) this parameter will be the same value as it was in the previous trajectory.
    Default: None

  • domain_decomposition_pattern (list of type int | Automatic | None) – The pattern how the domains should be arranged in a parallel simulations. E.g. [1, 2, 4] means 1 domain in A-, 2 in B-, and 4 in C-direction. If Automatic domain decomposition is used, then the simulation cell will be divided into domains whose edges are as close together in length as possible. If None is given then domain decomposition will be disabled.
    Default: Automatic
Returns:

The molecular dynamics trajectory.

Return type:

MDTrajectory

Usage Examples

Perform a molecular dynamics run of 500 steps of a bulk silicon crystal, using the Stillinger-Weber potential:

# -------------------------------------------------------------
# Bulk configuration
# -------------------------------------------------------------

# Set up lattice
vector_a = [5.4306, 0.0, 0.0]*Angstrom
vector_b = [0.0, 5.4306, 0.0]*Angstrom
vector_c = [0.0, 0.0, 5.4306]*Angstrom
lattice = UnitCell(vector_a, vector_b, vector_c)

# Define elements
elements = [Silicon, Silicon, Silicon, Silicon, Silicon, Silicon, Silicon,
            Silicon]

# Define coordinates
fractional_coordinates = [[0.0 ,  0.0 ,  0.0 ],
                          [0.25,  0.25,  0.25],
                          [0.5 ,  0.5 ,  0.0 ],
                          [0.75,  0.75,  0.25],
                          [0.5 ,  0.0 ,  0.5 ],
                          [0.75,  0.25,  0.75],
                          [0.0 ,  0.5 ,  0.5 ],
                          [0.25,  0.75,  0.75]]

# Set up configuration
bulk_configuration = BulkConfiguration(
    bravais_lattice=lattice,
    elements=elements,
    fractional_coordinates=fractional_coordinates
    )

# -------------------------------------------------------------
# Calculator
# -------------------------------------------------------------

potentialSet = StillingerWeber_Si_1985()
calculator = TremoloXCalculator(parameters=potentialSet)

bulk_configuration.setCalculator(calculator)
bulk_configuration.update()

# -------------------------------------------------------------
# Molecular Dynamics
# -------------------------------------------------------------

initial_velocity = MaxwellBoltzmannDistribution(
    temperature=300.0*Kelvin,
    remove_center_of_mass_momentum=True
)

method = NVEVelocityVerlet(
    time_step=1*femtoSecond,
    initial_velocity=initial_velocity
)

md_trajectory = MolecularDynamics(
    bulk_configuration,
    trajectory_filename='trajectory.nc',
    steps=500,
    log_interval=50,
    method=method
)

bulk_configuration = md_trajectory.lastImage()

# Save the final configuration
nlsave('final_configuration.nc', bulk_configuration)

# Get the list of the kinetic energies of the snapshots
kinetic_energies = md_trajectory.kineticEnergies()

# Get the list of the potential energies of the snapshots
potential_energies = md_trajectory.potentialEnergies()


md_example1.py

Usage example of the hook functionality to print the pressure and the cell volume to the log file at each step (e.g. in an NPT simulation).

# Define a hook function that prints the instantaneous cell volume and the pressure at each MD step
def print_cell_and_stress(step, time, configuration, forces, stress):

    # Get the observables from the configuration
    volume = configuration.bravaisLattice().unitCellVolume().inUnitsOf(Ang**3)

    # Calculate the hydrostatic pressure from the stress tensor
    pressure = -(stress[0, 0] + stress[1, 1] + stress [2, 2]) / 3.0

    # Print the output to the log file
    print("| MD step %i:  Pressure : %f bar    Volume:  %f Angstrom**3" % (step, pressure, volume))

# Run the MD simulation
md_trajectory = MolecularDynamics(
    bulk_configuration,
    constraints=[],
    trajectory_filename='trajectory.nc',
    steps=200,
    post_step_hook=print_cell_and_stress,
    log_interval=50,
    method=NPTMartynaTobiasKlein()
)


md_example2-full.py

Example of a post-step hook function for MD simulations using a domain decomposition (see Notes below). Here the z-coordinate of the first atom (index 0) is tethered to a fixed position using a harmonic spring potential. Note, that this example also works in serial simulations.

def parallelBiasForceHook(step, time, configuration, local_forces):
    # Add a bias force to the force acting on the first atom
    # in the z-direction.
    # First, check if the first atom is in the local domain.
    if 0 in configuration.localIndices():
        # Get the local index of the atom with first atom.
        local_index = numpy.where(configuration.localIndices() == 0)[0]

        # Get the local coordinates and calculate the bias force.
        local_positions = configuration.localCartesianCoordinates()
        bias_force = -0.5*eV/Ang**2*(local_positions[local_index, 2] - 2.0*Ang)

        # Add the bias force to the local forces vector.
        local_forces[local_index, 2] += bias_force

md_parallel_hook_example.py

Notes

  • The MolecularDynamics function returns an MDTrajectory object. After the simulation has completed, you can extract various properties of the stored snapshots, such as coordinates, forces, or kinetic energies, by using the MDTrajectory class methods.

  • You can perform custom operations on the configuration during the simulation, by using the pre_step_hook or post_step_hook functionality.

    The pre-step-hook one is invoked immediately before the force and stress calculations, which allows you to modify atomic positions or the cell vectors.

    The post-step-hook is invoked after the force and stress calculation, which means that you can modify the forces and stress, e.g. by adding a bias potential.

    Suitable functions must at least have the signature (step, time, configuration) where step denotes the integration step number, time is the current simulation time, and configuration is the current configuration. Optionally, additional arguments forces, local_forces, and stress can be passed. forces is the PhysicalQuantity array holding the current (global) forces, and stress is the stress tensor. local_forces are the forces on the atoms in local domain when run in distributed mode (see notes below), otherwise the array is identical to forces. Note, that you need to modify these two arrays inplace, as the return status of the function is ignored. You can use the class methods of the configuration object (see e.g. BulkConfiguration) to perform operations on the atoms.

    Moreover, it is possible to give a list of hook functions, which will be called sequentially at each MD step.

  • Constraints are specified by adding constraint objects, such as FixAtomConstraints, RigidBody, or FixCenterOfMass, to the constraints list.

    As an alternative to FixAtomConstraints, you can still use the old way of specifying the atomic indices of the fixed atoms in the constraints list. When using a thermostat to control the temperature, the reduced number of degrees of freedom due to the constraints is automatically taken into account.

Notes for Parallel MD Simulations

MolecularDynamics behaves differently when run in parallel with different types of calculators.

For DFT calculators (e.g. LCAOCalculator) running a MD simulation using more than one MPI process will parallelize over the energy, force and stress calculations as in a normal DFT calculation on a single configuration. The MD integration will be carried using an identical copy of the configuration on each process. In this case, when running MD with custom scripts (e.g. hook functions), you need to ensure that every MPI process follows the exact same behavior (e.g. when using random numbers, make sure to use a fixed seed).

For force field calculations (i.e. TremoloXCalculator), it is also possible to parallelize over multiple MPI processes. When more than one MPI process is used (and distributed mode is enabled, see below), the MD simulation will be executed using a domain decomposition. This means that each MPI process will only contain a sub-region of the configuration (a domain). The number of atoms in these local domains changes as atoms move across the domain boundaries during the MD simulation. This mode makes most sense for large configurations (> 10 000 atoms) as the domain size must be larger than the largest cutoff in the potential. The distributed mode will automatically be triggered if all of the following conditions are fulfilled:

  • There is more than one MPI process
  • The number of atoms in the configuration is larger than 1000 atoms per MPI process
  • All potential functions in the TremoloXCalculator support MPI

When these conditions are not met, the simulation will fall back to the serial mode. This means that each MPI process will run the same simulation and no parallel speedup will occur. For smaller system sizes, it may therefore be more efficient to use a single MPI process but many OpenMP threads instead.

When using hook functions in MD simulations, which are run with a domain decomposition, it is important to note that there are a few pitfalls, which can decrease the performance:

  • Using the forces argument in the signature will trigger gathering of all forces on all processes each time the hook function is called. Instead you should use local_forces if possible.
  • Calling configuration.cartesianCoordinates() or configuration.velocities() will trigger gathering of all positions respectively velocities on all processes. Instead, you should try to use configuration.localCartesianCoordinates() or configuration.localVelocities() if possible. These methods will return the arrays with the only the positions or velocities of the current atoms in the local domain (see the third usage example above).