MolecularDynamics

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

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. If a trajectory is given, then trajectory_interval and trajectory_object_id should be None.

  • 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 interval at which information, such as time, energy, temperature, etc. is written to the log output.
    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.

  • hook_functions (HookFunctions) – An optional HookFunctions object, which contains user-defined functions or a list of functions which will be called just before (pre-step) or just after (post-step) the forces evaluation. It can also contain measurement hooks, which are called at the end of each MD step after all constraints have been applied.
    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 simulation. 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.

  • trajectory_interval (int) – The resolution used in saving steps to a trajectory file. A value of 1 results in all steps being saved; a value of 2 results in every second step being saved; etc. If configuration is a trajectory (i.e. this is a restart MD simulation) this parameter will use the same value as was used in the previous trajectory.
    Default: The same value as log_interval.

  • trajectory_object_id (str | None) – The object id of the trajectory written to trajectory_filename. If a value of None is given, then an object id will be chosen automatically.
    Default: None.

  • scf_initial_guess (UsePreviousConfiguration | DIISExtrapolation) –

    Determines how the state of the system at each MD step is used as an initial guess for the next MD step, in the case when a self-consistent calculation is performed at each step. UsePreviousConfiguration means that the density matrix (or valence density/wave functions in a plane-wave calculation) of the previous configuration will be used as an initial guess for the next SCF calculation. For LCAO and semi-empirical calculations, it is possible and recommended to use DIISExtrapolation, which calculates an improved initial guess using the DIIS extrapolation method for the previous density matrices.

    This parameter has no effect when the calculator attached to configuration is a force-field calculator.
    Default: UsePreviousConfiguration for PlaneWaveCalculator, DIISExtrapolation(4) for LCAOCalculator and SemiEmpiricalCalculator.

  • pre_step_hook (callable | list of callable | None) – A 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. The return status is ignored. Unhandled exceptions will terminate a Molecular Dynamics evaluation. If a list is given, the functions will be called in the given order.
    Default: None.

  • post_step_hook (callable | list of callable | None) – A 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. Optional arguments include, forces, local_forces (for distributed MD), stress, trajectory (the MD trajectory), temperature, pressure, potential_energy, and/or kinetic_energy. The return status is ignored. Unhandled exceptions will terminate a Molecular Dynamics evaluation. If a list is given, the functions will be called in the given order.
    Default: None.

  • measurement_hook (callable | list of callable | None) – A user-defined function or a list of functions which will be called at the end of each MD step after all constraints have been applied. The signature of the function requires the arguments step, time, configuration. Optional arguments include forces, local_forces (for distributed MD), stress, trajectory (the MD trajectory), temperature, pressure, potential_energy, and/or kinetic_energy. The return value should be a dictionary that maps string keys to values. The values may be numbers, numpy arrays, or PhysicalQuantities. These values are stored on the MDTrajectory and can be accessed using the measurement method. Unhandled exceptions will terminate a Molecular Dynamics evaluation. If a list is given, the functions will be called in the given order.
    Default: None.

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 measurement 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, volume):

    # Convert volume to units of Ang**3.
    volume = volume.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,
    measurement_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.

  • When continuing a simulation, an existing MD trajectory object is used to initialize the configuration of the next simulation. In this case, the initial_velocities argument of the MD method should be set to ConfigurationVelocities.

  • 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.

    The measurement_hook is invoked last. It may return a dictionary that describes the data that should be saved to the MDTrajectory. The keys of the dictionary should be strings and the values should be numeric (e.g. numpy array, PhysicalQuantity, integer, or float).

    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.

    Measurement hooks may accept additional arguments beyond what is available to pre- and post-step hook functions:

    • temperature, the instantaneous kinetic temperature

    • pressure, the instantaneous pressure tensor

    • potential_energy, the current total potential energy of the system

    • kinetic_energy, the current total kinetic energy of the system

    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.

  • The inclusion of a homogeneous external electric field in bulk periodic DFT calculations with periodic boundary conditions is supported via the object ElectricFieldConstraint. For other geometries (i.e. slab or molecules), it is also possible to add an electric field by means of metallic regions in the Poisson solver (see Poisson solvers).

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).