Ensemble-Averaged Activation Strain Analysis

(1)$\Delta \overline{E} = \Delta \overline{E}_{\text{strain}} + \Delta \overline{E}_{\text{int}}$

Herein we describe an Ensemble-Averaged extension of the activation/strain analysis (ASA; also known as the distortion/interaction model), wherein the ASA is utilized for the analyses of entire molecular dynamics trajectories. The implementation utilizes CHARMM-style forcefields for the calculation of all energy terms.

Note

Throughout this document an overline will be used to distinguish between “normal” and ensemble-averaged quantities: e.g. $$E_{\text{strain}}$$ versus $$\overline{E}_{\text{strain}}$$.

Strain/Distortion

The ensemble averaged strain $$\Delta \overline{E}_{\text{strain}}$$ represents the distortion of all ligands with respect to their equilibrium geometry. Given an MD trajectory with $$m$$ iterations and $$n$$ ligands per quantum dot, the energy is averaged over all $$m$$ MD iterations and summed over all $$n$$ ligands.

The magnitude of this term is determined by all covalent and non-covalent intra-ligand interactions. As this term quantifies the deviation of a ligand from its equilibrium geometry, it is, by definition, always positive.

(2)$\Delta E_{\text{strain}} = E_{\text{lig-pert}} - E_{\text{lig-eq}} \quad \Rightarrow \quad \Delta \overline{E}_{\text{strain}} = \frac{1}{m} \sum_{i=0}^{m} \sum_{j=0}^{n} E_{\text{lig-pert}}(i, j) - E_{\text{lig-eq}}$
(3)$\Delta E_{\text{strain}} = \Delta V_{\text{bond}} + \Delta V_{\text{angle}} + \Delta V_{\text{Urey-Bradley}} + \Delta V_{\text{dihedral}} + \Delta V_{\text{improper}} + \Delta V_{\text{Lennard-Jones}} + \Delta V_{\text{elstat}}$

$$E_{\text{lig-eq}}$$ is herein the total energy of a (single) ligand at its equilibrium geometry, while $$E_{\text{lig-pert}}(i, j)$$ is the total energy of the (perturbed) ligand $$j$$ at MD iteration $$i$$.

Interaction

The ensemble averaged interaction $$\Delta \overline{E}_{\text{int}}$$ represents the mutual interaction between all ligands in a molecule. The interaction is, again, averaged over all MD iterations and summed over all ligand-pairs.

The magnitude of this term is determined by all non-covalent inter-ligand interactions and can be either positive (dominated by Pauli and/or Coulombic repulsion) or negative (dominated by dispersion and/or Coulombic attraction).

(4)$\Delta E_{\text{int}} = \sum_{j=0}^{n} \sum_{k \gt j}^{n} \Delta E_{\text{lig-int}} (j, k) \quad \Rightarrow \quad \Delta \overline{E}_{\text{int}} = \frac{1}{m} \sum_{i=0}^{m} \sum_{j=0}^{n} \sum_{k \gt j}^{n} \Delta E_{\text{lig-int}} (i, j, k)$
(5)$\Delta E_{\text{int}} = \Delta V_{\text{Lennard-Jones}} + \Delta V_{\text{elstat}}$

$$\Delta E_{\text{lig-int}}(i, j, k)$$ represents the pair-wise interactions between ligands $$j$$ and $$k$$ at MD iteration $$i$$. Double counting is avoided by ensuring that $$k > j$$.

Note

In order to avoid the substantial Coulombic repulsion between negatively charged ligands, its parameters are substituted with those from its neutral (i.e. protonated) counterpart. This correction is applied, exclusively, for the calculation of $$\Delta E_{\text{lig-int}}$$.

Total Energy

The total (ensemble-averaged) energy is the sum of $$\Delta \overline{E}_{\text{strain}}$$ and $$\Delta \overline{E}_{\text{int}}$$. Note that the energy is associated with a set of $$n$$ ligands, i.e. the distortion and mutual interaction between all $$n$$ ligands. Division by $$n$$ will thus yield the averaged energy per ligand per MD iteration.

(6)$\Delta \overline{E} = \Delta \overline{E}_{\text{strain}} + \Delta \overline{E}_{\text{int}} = \frac{1}{m} \sum_{i=0}^{m} \Delta E_{\text{strain}}(i) + \Delta E_{\text{int}}(i)$

Examples

An example input script using the Cd68Se55 core and OC(=O)CC ligand.

The activation_strain.md key enables the MD-ASA procedure; activation_strain.use_ff ensures that the user-specified forcefield is used during the construction of the MD trajectory.

path: ...

input_cores:
- Cd68Se55.xyz:
guess_bonds: False

input_ligands:
- OC(=O)CC

optional:
core:
dummy: Cl

ligand:
optimize: True
split: True

qd:
activation_strain:
use_ff: True
md: True
job1: Cp2kJob

forcefield:
charge:
keys: [input, force_eval, mm, forcefield, charge]
Cd: 0.9768
Se: -0.9768
O2D2: -0.4704
C2O3: 0.4524
epsilon:
unit: kjmol
keys: [input, force_eval, mm, forcefield, nonbonded, lennard-jones]
Cd Cd: 0.3101
Se Se: 0.4266
Cd Se: 1.5225
Cd O2D2: 1.8340
Se O2D2: 1.6135
sigma:
unit: nm
keys: [input, force_eval, mm, forcefield, nonbonded, lennard-jones]
Cd Cd: 0.1234
Se Se: 0.4852
Cd Se: 0.2940
Cd O2D2: 0.2471
Se O2D2: 0.3526


activation_strain

optional.qd.activation_strain

All settings related to the activation strain analyses.

Example:

optional:
qd:
activation_strain:
use_ff: True
md: True
iter_start: 500
dump_csv: False

el_scale14: 1.0
lj_scale14: 1.0

distance_upper_bound: "inf"
k: 20
shift_cutoff: True

job1: cp2kjob
s1: ...

forcefield:
...


optional.qd.activation_strain.use_ff
Parameter

Utilize the parameters supplied in the optional.forcefield block.

optional.qd.activation_strain.md
Parameter

Perform an ensemble-averaged activation strain analysis.

If True, perform the analysis along an entire molecular dynamics trajectory. If False, only use a single geometry instead.

optional.qd.activation_strain.iter_start
Parameter

The MD iteration at which the ASA will be started.

All preceding iteration are disgarded, treated as pre-equilibration steps. Note that this refers to the iteration is specified in the .xyz file. For example, if a geometry is written to the .xyz file very 10 iterations (as is the default), then iter_start=500 is equivalent to MD iteration 5000.

optional.qd.activation_strain.dump_csv
Parameter

Dump a set of .csv files containing all potential energies gathered over the course of the MD simulation.

For each quantum dot two files are created in the .../qd/asa/ directory, one containing the potentials over the course of the MD simulation (.qd.csv) and for the optimized ligand (.lig.csv).

optional.qd.activation_strain.el_scale14
Parameter

Scaling factor to apply to all 1,4-nonbonded electrostatic interactions.

Serves the same purpose as the cp2k EI_SCALE14 keyword.

optional.qd.activation_strain.lj_scale14
Parameter

Scaling factor to apply to all 1,4-nonbonded Lennard-Jones interactions.

Serves the same purpose as the cp2k VDW_SCALE14 keyword.

optional.qd.activation_strain.distance_upper_bound
Parameter

Consider only atom-pairs within this distance for calculating inter-ligand interactions.

Units are in Angstrom. Using "inf" will default to the full, untruncated, distance matrix.

optional.qd.activation_strain.k
Parameter

The (maximum) number of to-be considered distances per atom.

Only relevant when distance_upper_bound != "inf".

optional.qd.activation_strain.shift_cutoff
Parameter

Add a constant to all electrostatic and Lennard-Jones potentials such that the potential is zero at the distance upper bound.

Serves the same purpose as the cp2k SHIFT_CUTOFF keyword. Only relevant when distance_upper_bound != "inf".

optional.qd.activation_strain.job1
Parameter

A type object of a Job subclass, used for performing the activation strain analysis.

Should be set to Cp2kJob if activation_strain.md = True.

optional.qd.activation_strain.s1
Parameter
• Default value – See below

s1:
input:
motion:
print:
trajectory:
each:
md: 10
md:
ensemble: NVT
temperature: 300.0
timestep: 1.0
steps: 15000
thermostat:
type: CSVR
csvr:
timecon: 1250

force_eval:
method: FIST
mm:
forcefield:
ei_scale14: 1.0
vdw_scale14: 1.0
ignore_missing_critical_params: ''
parmtype: CHM
parm_file_name: null
do_nonbonded: ''
shift_cutoff: .TRUE.
spline:
emax_spline: 10e10
r0_nb: 0.2
poisson:
periodic: NONE
ewald:
ewald_type: NONE
subsys:
cell:
abc: '[angstrom] 100.0 100.0 100.0'
periodic: NONE
topology:
conn_file_format: PSF
conn_file_name: null
coord_file_format: 'OFF'
center_coordinates:
center_point: 0.0 0.0 0.0

global:
print_level: low
project: cp2k
run_type: MD


The job settings used for calculating the performing the ASA.

Alternatively, a path can be provided to .json or .yaml file containing the job settings.

The default settings above are specifically for the ensemble-averaged ASA (activation_strain.md = True.).