Code for Performing Semi-Classical Quantum Dynamics

The methods that are implimented in this code are : PLDM, spin-PLDM, MFE (Mean-Field Ehrenfest) and SQC. The present code works for slurm based High-Performance Computing Cluster (HPCC) as well as on personal computers.

Usage

Step 1

Create a folder and git clone this repository.

git clone https://github.com/arkajitmandal/SemiClassical-NAMD

Step 2

Code up the model system in a python file inside the "Model" folder and name it 'modelName.py'.

The 'modelName.py' should look like:

import numpy as np

class parameters():
   # some parameters
   # Nuclear Timestep (a.u.)
   dtN = 2

   # Number of Nuclear Steps 
   # length of simulation : dtN x Nsteps
   NSteps = 600  
   
   # Number trajectories per cpu
   NTraj = 100

   # Electronic Timestep (a.u.)
   # Please use even number
   dtE = dtN/40

   # Mass of nuclear particles in a.u.
   # is a vector of length NR : number 
   # of nuclear DOF
   M = np.array([1836.0])

   # Initial electronic state
   initState = 0
   
   # Save data every nskip steps
   nskip = 5

def Hel(R):
    # Diabatic potential energies Vij(R) in a.u., 
    # return Vij(R) : NxN Matrix 
    # (N  = number of electronic states)
    # R : Vector of length NR
    # (NR = number of nuclear degrees of freedom)
    # See details below 
    return Vij

def dHel0(R):
    # State independent gradient 
    # return dV0(R) : Vector of length NR
    # (NR = number of nuclear degrees of freedom) 
    # See details below
    return dV0

def dHel(R):
    # Gradient of the diabatic
    # potential energy matrix elements
    # Vij(R) in a.u.  
    # return dVij(R) : Array of dimention (N x N x NR)
    # See details below
    return dVij

def initR():
    # Provide initial values of R, P in a.u.
    # R, P : Both are vectors of length NR
    # (NR = number of nuclear degrees of freedom) 
    # See details below
    return R, P

def initHel0(R):
    #------ This part will be only for NRPMD----------
    #-------while running condensed phase system------
    #R : is a 2D array of dimensionality ndof,nb
         #where, nb is the bead and ndof is the number of dofs
    #M : mass of the particle
    #ω = frequency of the particle
    #R0 = initial positon of the photo-excitation 
    # see details below
    return  np.sum(0.5 *M* ω**2 * (R-R0)**2.0)

You can find several examples of model files inside the "Model" folder. I will explain each parts of this file in more detain in a section below.

Step 3 (simple | serial; on your pc/mac/linux)

Prepare an input file (name it : 'whateverInput.txt'):

Model                = tully2
Method               = pldm-focused 

• Model : The right hand side of the first line, tully2, tells the code to look for tully2.py inside the folder "Model". If you name your model file as modelName.py then you should write 'Model = modelName' (without the '.py' part).
• Method : Written as, method-methodOption. Select a quantum dynamics method. The available methods are :
  • mfe : Mean-Field Ehrenfest Approach. Kind of worst approach you can think of.
  • pldm-focused : Partial Linearized Density Matrix (PLDM) [1] with focused initial conditions. Should be similar to mfe. Maybe slightly better.
  • pldm-sampled : Partial Linearized Density Matrix (PLDM) [1] with sampled initial conditions or the original PLDM approach. Most of the time works well, sometimes does not. Very good if your potentials are Hermonic (like Spin-Boson systems)
  • spinpldm-all: The Spin-Mapping PLDM [2] approach with full sampling. Often better than PLDM. Reliable but slighly slow. If your initial electronic state is a pure state |i⟩⟨i| (you could start from a super position state, but you have to hack into this code to do that) use spinpldm-half to get the same result but much faster (by half).
  • spinpldm-half: The Spin-Mapping PLDM approach, but with our in-house approximation. Works perfectly if starting with an initial electronic state that is a pure state |i⟩⟨i| (you could start from a super position state, but you have to hack into this code to do that).
  • spinpldm-focused: The Spin-Mapping PLDM approach, approximated. Good for short-time calculation and to get a general trend for longer time.
  • sqc-square: The Symmetric Quasi-Classical Approach, with square window [3]. Better than MFE. Cannot use it for more than several electronic states.
  • sqc-triangle: The Symmetric Quasi-Classical Approach, with triangle window [4]. Better than sqc-square.
  • zpesqc-triangle: The zero-point energy corrected Symmetric Quasi-Classical Approach [5], with triangle window. As good as spin-PLDM or better.
  • zpesqc-square: The zero-point energy corrected Symmetric Quasi-Classical Approach [5], with square window. Slightly worse than zpesqc-triangle.
  • nrpmd-{n} : The non-adiabatic ring polymer molecular dynamics[6,7] framework for accurately captures nuclear quantum effects while predicting efficient short-time and reliable longer time dynamics. Good for electron/charge transfer dynamics. Here {"n"} represents the number of beads.

The output file containing population dynamics is 'method-methodOption-modelName.txt', for the above input file it would be:

pldm-focused-tully2.txt

Step 3 (advanced | parallel jobs; slurm on HPCC)

Prepare an input file (name it : 'whateverInput.txt') for slurm submission in your computing cluster:

Model                = tully2
Method               = pldm-focused 

System               = slurm
Nodes                = 2
Cpus                 = 24
Partition            = action    

For first two lines see previous section.

Last four lines provide additional commands for slurm submission. For adding additional slurm ('#SBATCH') command, add them in the preamble of the 'parallel.py'. The default preamble looks like:

#!/usr/bin/env python
#SBATCH -o output.log
#SBATCH -t 1:00:00

Please dont add lines like "#SBATCH -N 1" or "#SBATCH --ntasks-per-node=24" in the preamble as they are declared in the input file.

The output file containing population dynamics is 'method-methodOption-modelName.txt', for the above input file it would be:

pldm-focused-tully2.txt

Step 4

Run the this code with python3.

python3 run.py whateverInput.txt

Where 'whateverInput.txt' is the input file described above. If your inputfile is named 'input.txt' then you could also just run,

python3 run.py

Details of Model Hamiltonian

In all of the approaches coded up here, the nuclear DOF {R,P} are evolved classically (their equation motion evolves under a classical like force) and the electronic DOF are described with the diabatic electronic states {|i⟩}.

A molecular Hamiltonian in the diabatic representation is written as:

where P_k is the momentum for the k th nuclear degrees of freedom with mass M_k. Further, V₀({R_k}) and V_ij({R_k}) are the state-independent and state-dependent part of the electronic Hamiltonian H_el({R_k}) in the diabatic basis {|i⟩}. That is: ⟨i| Ĥ - ∑_k P²_k/2M_k |j⟩ = V_ij({R_k}) + V₀({R_k})δ_ij. Of course most of times, we wave our hands, and make up models that describe V_ij({R_k}) with some analytical functions of {R_k}. If you know the analytical form of V_ij({R_k}) you can write a model file: modelName.py.

One can always, set V₀({R_k}) = 0, and instead redefine V_ij({R_k}) ⟶ V_ij({R_k}) + V₀({R_k})δ_ij and they should be equivalent in principle. However, some of the semiclassical approaches (pldm-sampled, sqc-square and sqc-triangle) produce results that depend on how one separates the state-independent and state-dependent parts of the gradient of the electronic Hamiltonian. This is because, this separation provides state dependent and independent gradients : ∇_kV₀({R_k}) and ∇_kV_ij({R_k}) which will provide different forces on nuclear particle for different choices of separations in some of these approximate quantum dynamics approaches. The nuclear forces computed in all of these approaches assumes this general form:

where Λ_ij vaguely resembles the electornic density matrix elements. For example, in MFE, Λ_ij = c_i*c_j. For methods that have ∑_iΛ_ii = 1 (like MFE) for individual trajectories this separation of state-dependent and independent does not matter. For other's as I said before, it does. In my experience, the more you can put in the independent part the better.

When such separation is not aparent, one can separate state-independent and state-dependent parts in the following manner. Consider the following molecular Hamiltonian with no aparent state-independent part of the electronic Hamiltonian:

One can define a state-independent part as:

and consequently the new state-dependent part (with matrix elements V'_ij({R_k}) ) becomes:

With this choice one gets the following state dependent and independent part of gradients: ∇_kV₀({R_k}) and ∇_kV'_ij({R_k}).

Details of a model file ('modelName.py')

Hel(R)

In the Hel(R) function inside the 'modelName.py' one have to define NxN matrix elements of the state-dependent electronic part of the Hamiltonian. Here you will code up V_ij({R_k}).

dHel(R)

In the dHel(R) function inside the 'modelName.py' one have to define NxNxNR matrix elements of the state-dependent gradient electronic part of the Hamiltonian. Here you will code up ∇_kV_ij({R_k}).

dHel0(R)

In the dHel0(R) function inside the 'modelName.py' one have to define a array of length NR describing the state-independent gradient of electronic part of the Hamiltonian. Here you will code up ∇_kV₀({R_k}).

NOTE: You dont need to code up V₀({R_k}) as such a zero-point energy shift does not affect dynamics.

initR()

Sample R, P from a wigner distribution.

---

to be continued...

References

---

[1] Huo and Coker J. Chem. Phys. 135, 201101 (2011)
[2] Mannouch and Richardson J. Chem. Phys. 153, 194109 (2020)
[3] Cotton and Miller J. Chem. Phys. 139, 234112 (2013)
[4] Cotton and Miller J. Chem. Phys. 145, 144108 (2016)
[5] Cotton and Miller J. Chem. Phys. 150, 194110 (2019)
[6] S. N. Chowdhury and P.Huo J. Chem. Phys. 147, 214109 (2017)
[7] S. N. Chowdhury and P.Huo J. Chem. Phys. 150, 244102 (2019)

email: arkajitmandal@gmail.com