Molecular Dynamics¶

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)"]

Introduction to 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_e\) and \(\Psi_n\) 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 [@Jensen2010]:

\[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.

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.

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