Ensemble-Averaged Activation Strain Analysis
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.
\(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).
\(\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.
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:
Type -
bool
Default value –
False
Utilize the parameters supplied in the
optional.forcefield
block.
- optional.qd.activation_strain.md
- Parameter:
Type -
bool
Default value –
False
Perform an ensemble-averaged activation strain analysis.
If
True
, perform the analysis along an entire molecular dynamics trajectory. IfFalse
, only use a single geometry instead.
- optional.qd.activation_strain.iter_start
- Parameter:
Type -
int
Default value –
500
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:
Type -
bool
Default value –
False
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:
Type -
float
Default value –
1.0
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:
Type -
float
Default value –
1.0
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
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:
Type -
int
Default value –
20
The (maximum) number of to-be considered distances per atom.
Only relevant when
distance_upper_bound != "inf"
.
- optional.qd.activation_strain.shift_cutoff
- Parameter:
Type -
bool
Default value –
True
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
A
type
object of aJob
subclass, used for performing the activation strain analysis.Should be set to
Cp2kJob
ifactivation_strain.md = True
.
- optional.qd.activation_strain.s1
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: MDThe 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
.).