Molecular Dynamics¶

quadrantChart
    title 
    x-axis Classical Mechanic Forces --> Quantum Electronic Structure
    y-axis Quantum Nuclear Effects --> Classical Motion of Nucleus
    quadrant-1 ab-initio MD
    quadrant-2 Classical MD
    quadrant-3 Classical PIMD
    quadrant-4 ab-initio PIMD

Introduction to ab-initio/Born-Oppenheimer Molecular Dynamics¶

The Born-Oppenheimer approximation is a fundamental quantum physical approximation applied to nuclei-electron systems which atoms and molecules are made of. The approximation allows to separate the nuclei and electron wavefuctions \(\Psi_n\) and \(\Psi_e\) to be analysed separately due to the very high mass difference between a nucleus of an atom and an electron (\(m_n >> m_e\)) ¹:

\[\Psi(r,R)_{total} = \Psi_e(r) \Psi_n(R)\]

And as a reminder: The wavefunction squared \(|\Psi|^2\) is the probability for an electron being at the position r or a nucleus for being at the position R in 3 dimensional space. Thus, r and R are 3 column matrices with the xyz coordinates of all electrons and nuclei being conventionally stored row-wise. Here an example of the row-wise xyz positioning of the atoms (nuclei as the center point) of a water molecule (excerpt of a Quantum Espresso input file):

ATOMIC_POSITIONS angstrom
O 2.525856153 1.482855778 8.928868428
H 2.241743628 2.141988100 9.587993065
H 3.499400219 1.513156121 8.959138507

The time-dependent Schrödinger equation is:

\[i \hbar \frac{d}{dt} \ket{\Psi(r,R,t)} = H \ket{\Psi(r,R,t)}\]

To track the motion of the nuclei and electrons via Newton's equations of motion over time one has to solve the Schrödinger equation with fixed nuclei in space, meaning the nuclei variables are defined as parameters and do not vary for each time step.

This separatation via the Born-Oppenheimer approximation leads to the following electronic and nuclear terms which add up to the total energy of the molecular system ²:

\[E = E_n + E_{el}\]

The nuclear Schrödinger equation which is solved to retrieve the nuclear Energy \(E_n\) is:

\[ \hat{H} \Psi(R) = E_{n} \Psi(R)\]

Analogous to the retrieval of \(E_n\), the electronic Schrödinger equation is solved for obtaining \(E_e\):

\[ \hat{H} \Psi(r) = E_{el}(n) \Psi_n(r)\]

The wavefunction still depends on the nucleus position and during a molecular dynamics calculation the nucleii positions are fixed for each time step. To analyse the dynamic behaviour of Born-Oppenheimer approximated molecular systems one can solve the electronic Schrödinger equation in combination the classical mechanics equations by Newton. The Newton equations describe the motion of the nuclei and are:

\[ M_I R_{I} = - \nabla_I min_{\psi_i} \braket{\Psi_0 | H_e^{HF} | \Psi_0} \]

\[ 0 = - H^{HF}_e \psi_i + \sum_j \Lambda_{ij} \psi_i \]

The set of equations for calculating the Hartree-Fock \(H_e^{HF}\) energy could also be replaced with the Kohn-Sham effective one-particle Hamiltonian with the set of equations derived from Hohenberg-Kohn-Sham density functional theory. The set of equations are solved for the ground state within the self-consistent field (SCF) scheme for each increasing trajectory time step with fixed nuclear positions (R).

By also having the Born-Oppenheimer approximation included in the system of equations and solved via a finite difference scheme such as the Verlet algorithm, one can take into account quantum effects such as zero point motion, quantum tunneling or quantum fluctuations.

For Car-Parrinello MD (CPMD) the SCF solution is obtained by solving the Khon-Sham equations only for the initial nuclear configuration (time step t=0) and the wavefunction for the following steps is propagated along the atoms' trajectories via the Lagrangian scheme. CPMD also applies the Born-Oppenheimer approximation and this can lead to confusion since in the literature one has the established nomenclature of BOMD and CPMD.

flowchart TD
    A[Time-dependent Schrödinger equation] -->|"Born-Oppenheimer (BO) approximation"| B("BO <br> Molecular Dynamics (MD) equations")
    B --> C{"classical or <br> quantum mechanically treated nuclei"}
    C -->|Classical treated nuclei| D["Born-Oppenheimer MD (BOMD) <br> Car-Parrinello MD (CPMD)"]
    C -->|Quantum mechanical treated nuclei| F["Path Integral MD (PIMD)"]

How to send off a task plane wave DFT molecular dynamics task to the Cebule engine¶

Under the hood the Quantum Espresso package is deployed for running BOMD and CPMD calculations.

You can send a Quantum Espresso input file in the following way to Cebule:

# Submit the contents of "....in" file as a new task
with open("....in", "r") as file:
    task = session.cebule.create_task("My Task", TaskType.BORN_OPPENHEIMER_MD, file.read())
    print(task)

Get the status of the task via its ID¶

task = session.cebule.task_status(task.id)
print(task)

Get the results of the last task¶

result = session.cebule.get_result(task.id, "stdout")
print(result.decode("utf-8"))

Cancel a running task by its ID¶

task = session.cebule.cancel_task(task.id)
print(task)

Get the available compute time for the user¶

quota = session.cebule.get_quota()
print(quota)

The whole code repository including the markdown documentation and example notebooks of the MQS Python SDK can be found here: https://gitlab.com/mqsdk/python-sdk BOMD and CPMD examples can be found in the notebooks folder.

Force Field Molecular Dynamics¶

Force Field Molecular Dynamics (FFMD) simulates molecular systems by numerically integrating Newton's equations of motion using classical force fields to model atomic interactions. This approach enables resource-efficient studies of dynamic behavior, equilibrium properties, and structural characteristics of complex molecular systems.

Key Components¶

OpenFF SAGE Force Field (v2.2.1): A high-quality general-purpose molecular mechanics force field developed by the Open Force Field Initiative. It provides parameters for a wide range of organic molecules and has been extensively validated against experimental and quantum mechanical data in the SAGE force field paper. SAGE is the second-generation force field from OpenFF, featuring comprehensively refit Lennard-Jones parameters against condensed phase data and improved valence parameters.
Langevin Integrator: An enhanced molecular dynamics integrator that incorporates random forces and frictional terms to mimic the effect of solvent molecules, maintaining constant temperature through stochastic collisions.
Monte Carlo Barostat: Controls pressure by periodically attempting random changes to the simulation cell volume, accepting or rejecting changes based on the Metropolis criterion to maintain a target pressure.

Simulation Ensembles¶

NPT (constant Number, Pressure, Temperature): Simulates a system with fixed number of particles, pressure, and temperature. This ensemble is particularly useful for condensed phase systems and is the default in our implementation.
NVT (constant Number, Volume, Temperature): Simulates a system with fixed number of particles, volume, and temperature. This ensemble is useful for systems where volume changes are not expected.

Cebule SDK TaskType: FORCE_FIELD_MD¶

Performs molecular dynamics simulation with OpenFF SAGE force field using OpenMM
Inputs:
smiles_primary: str - SMILES of primary molecule (the molecule to calculate average forces on)
copies_primary: int - Number of copies of primary molecule. The output forces will be the average forces across copies, reduced to the length of a single primary molecule.
smiles_list_secondary: List[str] - SMILES of secondary molecules (e.g., solvents, ions)
copies_list_secondary: List[int] - Number of copies of each secondary molecule
temperature: float = 293.15 - Simulation temperature in Kelvin
box_length_nm: float = 3.5 - Simulation periodic box length in nanometers
time_fs: int = 100000 - Simulation time in femtoseconds
Cebule max_processors: Controls the number of CPU cores used for the simulation
Output: List of average forces (magnitude) for each atom in the primary molecule after equilibration (from the last 1/3 of the run). See Atom Order which defines the order of atoms in this outputted forces list.

Example:

task_ffmd = session.cebule.create_task("ffmd", 
                                       TaskType.FORCE_FIELD_MD, 
                                       smiles_primary="CCO", 
                                       copies_primary=10, 
                                       smiles_list_secondary=["O", "[Na+].[Cl-]"], 
                                       copies_list_secondary=[1000, 30], 
                                       time_fs=100000,
                                       max_processors=64)

The simulation places the specified molecules in a periodic box, equilibrates the system, and returns the average forces on the primary molecule's atoms after equilibration. This can be useful for understanding molecular behavior in solutions, examining solvent effects, and studying molecular interactions.

J. Grotendorst, S. Blügel, Forschungszentrum Jülich, and D. Marx. Computational Nanoscience: Do it Yourself!: Winter School, 14 - 22 February 2006, Forschungszentrum Jülich, Germany ; Lecture Notes. NIC series. NIC, 2006. ISBN 9783000173509. URL: https://books.google.dk/books?id=XZpKBAAACAAJ. ↩
Jan Halborg Jensen. Molecular modeling basics. CRC Press, 2010. ISBN 978-1-4200-7526-7. ↩