In this paper, we present CTRAMER (Charge-Transfer RAtes from Molecular dynamics, Electronic structure, and Rate theory)—an open-source software package for calculating interfacial charge-transfer (CT) rate constants in organic photovoltaic (OPV) materials based on ab initio calculations and molecular dynamics simulations. The software is based on identifying representative donor/acceptor geometries within interfacial structures obtained from molecular dynamics simulation of donor/acceptor blends and calculating the corresponding Fermi's golden rule CT rate constants within the framework of the linearized-semiclassical approximation. While the methods used are well established, the integration of these state-of-the-art tools originating from different disciplines to study photoinduced CT processes with explicit treatment of the environment, in our opinion, makes this package unique and innovative. The software also provides tools for investigating other observables of interest. After outlining the features and implementation details, the usage and performance of the software are demonstrated with results from an example OPV system.

The low production and environmental costs, as well as improved plasticity and synthetic tunability, of organic materials in comparison to their inorganic counterparts motivate the development of photovoltaic devices based on organic materials [so called organic photovoltaics (OPV)].1–13 Within OPV devices, photoexcitation of the donor material leads to the formation of excitons. The excitons then diffuse to the donor/acceptor (D/A) interface, where charge transfer (CT), namely, electron transfer from the donor to the acceptor, occurs.1–3,14 This is followed by charge separation, where the electron and hole diffuse away from the interface within the donor and acceptor layers, respectively. Thus, a better understanding of the correlation between the effect of varying interfacial D/A pair geometries on CT rate constants is required to improve OPV performance.1,6,15–18

In this paper, we introduce a new software package, ctramer (Charge Transfer RAtes from Molecular dynamics, Electronic structure, and Rate theory), that provides computational tools for correlating interfacial CT rates with the underlying interfacial structure. The approach that has been benchmarked and employed by us in previous work1,19 combines state-of-the-art electronic structure calculations and molecular dynamics (MD) simulations to compute representative interfacial D/A geometries and the corresponding CT rate constants. The CT rate constants are calculated within the framework of Fermi’s golden rule (FGR) and based on the linearized-semiclassical (LSC) approximation.19–25 Support for other levels of CT theory is planned to be added in future versions. Each of the methods used in ctramer has been chosen due to being well-studied and performing well in benchmarks. It is the combination of these state-of-the-art methods from different fields that makes ctramer unique.

It should be noted that the FGR/LSC framework currently used for calculating CT rate constants in ctramer is based on treating the environment of the D/A pair at the molecular level, as opposed to treating it as a polarizable continuum or a harmonic bath.1,21 Such resolution is required to account for the heterogeneity of the solid state environment and the distribution of D/A geometries and CT rates it can give rise to.1,15,20,27–30 It should also be noted that the molecular models we use are parameterized based on inputs obtained from electronic structure calculations.

Our primary goal in this paper is to introduce a general-purpose software package based on the computational framework summarize its available features in it, and demonstrate its applicability and scope by presenting the results from an example OPV system where boron subphthalocyanine (SubPC) serves as the donor and fullerene (C60) serves as the acceptor.2,4,6,12,31–37 This framework could be broadly applied to other materials.

The overall workflow for ctramer is outlined in Fig. 1. The algorithm is initiated by the molecular coordinates of the donor and acceptor [Fig. 1(a)]. The output corresponds to CT rate constants for different interfacial D/A geometries, which can be used to correlate the interfacial CT rate constants to the interfacial structure.

FIG. 1.

The overall workflow for ctramer (illustrated on the SubPC/C60 system). (a)/(b) corresponds to the C60/SubPC molecules, respectively, where the blue dashed enclosures represent electronic-structure calculations done for each molecule separately. (c) A condensed-phase system constructed using the molecular coordinates from (a) and (b) and the Mulliken charges obtained therein. (d) An equilibrated snapshot of the (c) system after equilibration. (e) The potential of mean force for the system in (d), where a coordinate in (R,θ) space corresponds to a SubPC/C60 pair. (f) A SubPC/C60 pair selected from (e), with the blue dashed enclosure representing electronic-structure calculations done for the pair. (g) Scatterplot of characteristics for excited states calculated in (f). (h) Distribution of fluctuations in the energy gap of transitions between states selected from those calculated in (f). (i) Charge-transfer rate-constant densities for the configuration from (f)–(h) (in this case, hollow) and two others. The numbers for each model vary according to the population of each configuration.

FIG. 1.

The overall workflow for ctramer (illustrated on the SubPC/C60 system). (a)/(b) corresponds to the C60/SubPC molecules, respectively, where the blue dashed enclosures represent electronic-structure calculations done for each molecule separately. (c) A condensed-phase system constructed using the molecular coordinates from (a) and (b) and the Mulliken charges obtained therein. (d) An equilibrated snapshot of the (c) system after equilibration. (e) The potential of mean force for the system in (d), where a coordinate in (R,θ) space corresponds to a SubPC/C60 pair. (f) A SubPC/C60 pair selected from (e), with the blue dashed enclosure representing electronic-structure calculations done for the pair. (g) Scatterplot of characteristics for excited states calculated in (f). (h) Distribution of fluctuations in the energy gap of transitions between states selected from those calculated in (f). (i) Charge-transfer rate-constant densities for the configuration from (f)–(h) (in this case, hollow) and two others. The numbers for each model vary according to the population of each configuration.

Close modal

ctramer consists of five modules as shown in Fig. 2. The modules address the different scales needed for describing the CT process occurring at a D/A core that is affected by its molecular environment. The modular nature of ctramer allows the user to exclude or replace steps as needed. To initiate the calculation, the individual donor and acceptor molecules [see Figs. 1(a) and 1(b)] atomic coordinates are obtained from the literature or experiment. In Module 1, the atomic charges of the donor and acceptor molecules are calculated. (Alternatively, the charges can be assigned using published force fields.) Module 2 combines multiple donor and acceptor molecules [Fig. 1(c)] and equilibrates the overall system using the charges from Module 1 to a pre-assigned temperature and pressure [Fig. 1(d)]. This is followed by a determination of the distribution of interfacial D/A pair geometries in the form of multiple parallel trajectories [Fig. 1(e)]. Module 3 performs electronic state calculations on selected representative interfacial D/A pairs that correspond to different classes of geometries [Fig. 1(f)]. From these calculations, excited donor and acceptor states are identified using preset criteria of excited state properties [oscillator strength (OS) and CT characteristics] [Fig. 1(g)]. Here, the relevant transitions and the coupling between the donor and acceptor states are obtained. Module 4 uses MD simulations to calculate the fluctuations in the D/A energy gap in the condensed phase [Fig. 1(h)]. Combining these fluctuations with CT characteristics from Module 3, LSC-based E-FGR CT rates [Fig. 1(i)] are obtained by Module 5 at the chosen level of theory.

FIG. 2.

A flowchart representing information flow during each section of the software. Each box represents a process, while the arrows represent data within the software. Input parameters are the number of representative structures from MD (N1), excited states from Q-Chem (N2), and excitonic (N3)/CT (N4) states for which transitions are selected.

FIG. 2.

A flowchart representing information flow during each section of the software. Each box represents a process, while the arrows represent data within the software. Input parameters are the number of representative structures from MD (N1), excited states from Q-Chem (N2), and excitonic (N3)/CT (N4) states for which transitions are selected.

Close modal

Figure 3 illustrates the actual scripts that comprise each module in ctramer. The electronic structure calculations and the MD simulations reported in this manuscript are based on Q-Chem 4.038 and AMBER12,39 respectively. All the scripts are available on GitHub40 (https://github.com/ctramer/ctramer). Below, we provide additional information on the procedures for establishing the different parameters under each module.

FIG. 3.

The execution of the scripts in order (available on GitHub: https://github.com/ctramer/ctramer).

FIG. 3.

The execution of the scripts in order (available on GitHub: https://github.com/ctramer/ctramer).

Close modal

Time-dependent density-functional theory (TD-DFT) calculations are used in two modules of ctramer. First, in Module 1, TD-DFT is used on individual molecules to obtain partial charges for use in the MD simulations of Module 2. In Module 3, TD-DFT is used on selected interfacial D/A molecule pairs to obtain excited states and their partial charges, oscillator strength (OS), and relative energy.41 These electronic structure protocols are benchmarked against experimental measured excitation energies27,42,43 (including those of charge transfer states) and measured rates.2,44 We note that, in the condensed phase, as molecules tend to neighbor several molecules, each molecule can be involved in multiple D/A pairs. We hypothesize that, on the macroscopic scale, CT follows the paths involving the D/A pairs with the fastest rate constants on the microscopic scale.

Important excited states are selected and then classified as donor and/or acceptor states. ctramer calculates the electronic coupling coefficients between electronic states using the fragment-charge-difference (FCD) method.45 The choice of functional and basis set can be customized for the system under study based on the literature and experimentation.

For the results presented in this paper, ctramer uses the 6-31G* basis set46 and Baer–Neuhauser–Livshits (BNL) range-separated hybrid functional47,48 for TDDFT calculations. A γ value of 0.167 bohr−1, tuned to an optimal value for the on-top geometry based on the J2(γ) scheme,49 is used for all the geometries.

The calculated excited states are classified as follows: First, the charge of the donor molecule, QD, is used to classify states either as non-CT (QD < 0.25 e) or CT (QD > 0.25 e). Second, states with a significant OS are referred to as light-absorbing states or bright.50 CT states with negligible OS are referred to as dark (dCT), and those with a significant OS are addressed as bright CT (bCT). Non-CT states with significant OS are referred to as excitonic (EX). As ctramer studies direct CT from photoexcited (bright) states to CT states, states that are both dark and non-CT are not addressed. However, extensions to consider additional processes can be developed in future versions. The states are then named as EXn, bCTn, or dCTn. The index n refers to the rank of a state’s energy, from the lowest to the highest, within EX, bCT, and dCT states at the same geometry. Classification thresholds can be customized within ctramer.

Finally, using these classified states, transitions from a bright state to a CT state are selected. The user selects the maximum number of donor and acceptor states between which transition rates are calculated. The computational cost is linear with the number of donor states but is largely unaffected by the number of acceptor states (due to the method of calculating rates detailed in Subsection II E). As part of these electronic structure calculations, electronic coupling coefficients are computed for each transition.

Module 2 focuses on creating a condensed-phase system from the coordinates and charges established in Module 1, equilibrating it to a given temperature and pressure and then analyzing the distribution of structures. From this distribution, representative structures are selected for further analysis along with their surrounding molecules. In simulations during Module 4, an interfacial molecular pair is kept fixed for use in calculating energy gaps between excited states. These energy gaps are used in Module 5 to calculate rate constants.

ctramer uses Packmol51 to construct condensed-phase OPV systems. Packmol places a user-defined number of molecules into a region of space, the size of which is also set by the user. ctramer by default uses six square layers into each of which it places 25 molecules. Each layer consists of only one type of molecule. Space is placed between separate types of layers to simulate the fabrication procedure. Too much or too little space between molecules can cause the system to not equilibrate correctly.

The MD simulations of ctramer utilize the generalized AMBER force field (GAFF),39,52,53 augmented by other force fields (FFs) as needed. For example, interactions involving the boron atom that is found at the center of SubPC are not given in the GAFF and are here taken from Refs. 54 and 55. The MD simulations are performed using the amber12 program sander.39ctramer fits the simulation box to the constructed system and then applies periodic boundary conditions. Both equilibration and production MD simulations have a time step of 2.0 fs. The SHAKE algorithm56 is used to constrain bonds involving hydrogen. Neighbor list updates, real-space Coulomb interactions, and van der Waals interactions utilize a threshold of 10.0 Å. ctramer calculates electrostatic interactions with the particle-mesh Ewald method.57,58

By default, a hybrid algorithm is used to perform energy minimization. The default settings for this algorithm, used for the results in this paper, are listed here. The conjugated-gradient method is used for 500 time steps, followed by 4500 time steps of the steepest-descent method. The minimized system is then heated to 298.15 K using the canonical ensemble (NVT) gradually over a period of 10 ps. ctramer then fixes the system at this temperature for 1.0 ns and equilibrates at a constant pressure of 1.0 bar using a NPT ensemble with 1.0 ps as the pressure relaxation time. The simulation box is then refit to the equilibrated system with periodic boundary conditions being kept. A planned future update is for ctramer to simulate automatically until equilibration is reached.

Production runs in Module 2 use the equilibrated system with the adjusted box size in a NVT ensemble. As with all other parameters, the number and length of production runs can be specified by the user. For the results shown in this paper, one run of 20 ns was used. In our tests on ctramer, multiple production runs can aid in sampling multiple local minima.

Production runs in Module 4 use the entire system in the time step in Module 2 from which a representative geometry was chosen. The atoms of the representative pair are fixed in place using a harmonic potential. By default, ctramer uses a force constant of 50 kcal/mol Å2. However, as the energies calculated by ctramer do not include the restraint potential, the size of the force constant will not affect the results.39,59 The goal of fixing the atoms is to keep the molecules in the representative geometry. These fixed atoms are assigned with Mulliken partial charges, from Module 3, of a donor state for that representative pair geometry. The rest of the molecules in the system have ground state charges from Module 1 and are unrestrained. After equilibration, production runs are performed using the same parameters as Module 2. The effect on uncertainty in rate constants in Module 5 by changing the number and length of production runs in Module 4 is discussed in Sec. II E.

Module 5 can evaluate CT rates at different levels. Below, transition rate constants are calculated following a Marcus-level linearized semi-classical (LSC)1,19,21,22 approximation to Fermi’s golden rule, where the donor-to-acceptor transition rate constant (kM) is given by

kM=|ΓDA|22πσU2expU22σU2.
(1)

Here, UR=VDRVAR is the D/A energy gap as a function of the nuclear coordinates, R, where VD and VA are the potential-energy surfaces of the donor and acceptor states, respectively. U and σU are the first and second moments of UR, respectively. This rate constant is derived from Fermi’s golden rule by the LSC approximation, as detailed in Refs. 19, 21, 23, and 24, but is referred to as the “Marcus-level” as the reorganization energy, Er, the reaction free energy, ΔE, and the activation energy, Ea, can be calculated from U and σU: Er=σU2/(2kBT), ΔE=ErU, and Ea=kBTU2/(2σU2). These parameters allow for the analysis of rates calculated using Marcus theory.60–63 The electronic excitation energies are from the Module 3 electronic structure calculations.

UR is determined for each transition using the classically sampled trajectories from Module 4 as detailed in Sec. II D. First, the potential-energy surface from MD at a given time step, VαM, is calculated for each state, with α denoting the donor or the acceptor, by recalculating the energy of the entire system. To avoid double counting of potential energy by the electronic structure and MD methods, VαM is corrected to Vα=VαM+ Wα, where Wα is the difference in the single-point energy of each state between calculations by the electronic structure and MD methods. U is then determined using the difference of VD and VA. Finally, U and σU are obtained using the moments of NR MD trajectories of length LR. We here use NR = 40 runs and LR = 40 ns, which can be adjusted for the desired accuracy. The uncertainty for rate constants decreases approximately as 1NRLR.

As shown in Fig. 1(h) and Fig. S1, the distribution of UR can be significantly non-Gaussian for many transitions. We attribute this to ctramer accounting for how the condensed phase’s heterogeneity can lead to multiple local energy minima. As the LSC approximation uses a Gaussian distribution to model the probability density at UR=0, ctramer resolves these energy minima by best fitting the probability density function of UR to a sum of Gaussian distributions chosen by least-squares regression. ctramer increases the number of distributions, beginning at 1, until the 95% confidence interval for σU includes non-positive numbers or a maximum of three is reached. These settings can be customized within ctramer.

While kM measures the transition rate between states, we multiply it by the amount of the corresponding charge transferred, (ΔQD), to measure the rate of CT (kC) for a transition as follows:

kC=ΔQDkM.
(2)

Both kC and kM are obtained in the context of a single transition for one structure. For aid in comparison between structures, these CT rates are summed over all the identified transitions, t, for a given representative structure i, to give a structure-level CT rate constant (KiC) as follows:

KiC=tktC.
(3)

A system-level CT density (ωC) can then be established by averaging KC over the area of the D/A interface in the simulation as follows:

ωC=iKiCni×1A,
(4)

where ni is the number of pairs represented by the structure i from Eq. (3) and A is the approximate area of the interface (which is by default calculated using ctramer but can replaced).

ctramer is used to analyze the CT rate in the interface of SubPC/C60, pair of donor–acceptor organic materials used in model OPV studies. ctramer can be used for many materials, but these results are presented here as an example use of the software. Data are available at https://github.com/ctramer/ctramer.40 

The coordinates of the optimized SubPC and C60 molecules are provided in GitHub40 (https://github.com/ctramer/ctramer), while the references used to determine the atomic charges are listed in Table S1.

The multilayered OPV system is represented using six alternating layers of 25 SubPC or C60 molecules each. A large ensemble of interfacial D/A pairs is then obtained, where a interfacial D/A pair is defined as 1 SubPC molecule and 1 C60 molecule where the minimum distance between any atom from separate molecules is less than 5 Å.

Analysis of the different interfacial pairs is aided by two order parameters [see Fig. 4(b)]: first, the distance (R) between the SubPC boron atom and the C60 center of mass and second, the angle (θ) between the vector from the C60 center of mass to the SubPC boron atom and the vector from the SubPC boron atom to the SubPC chlorine atom.

FIG. 4.

The potential of mean force (a) for the SubPC/C60 pair on a Rθ coordinate system. The color is scaled by kBT. Also shown are the representative structures (on top, hollow, and edge) of the SubPC/C60 pair at each of the three major geometries and the percentage of interfacial pairs in that geometry. R and θ are defined in (b), where the yellow bead corresponds to the C60 center of mass and the red and blue beads correspond to the SubPC boron and SubPC chlorine atoms, respectively. R is the distance between the SubPC boron atom and the C60 center of mass. θ is the angle between vectors V1 and V2.

FIG. 4.

The potential of mean force (a) for the SubPC/C60 pair on a Rθ coordinate system. The color is scaled by kBT. Also shown are the representative structures (on top, hollow, and edge) of the SubPC/C60 pair at each of the three major geometries and the percentage of interfacial pairs in that geometry. R and θ are defined in (b), where the yellow bead corresponds to the C60 center of mass and the red and blue beads correspond to the SubPC boron and SubPC chlorine atoms, respectively. R is the distance between the SubPC boron atom and the C60 center of mass. θ is the angle between vectors V1 and V2.

Close modal

Inspection of the potential of mean force (PMF) using a (R,θ) coordinate system [shown in Fig. 4(a)] reveals two pronounced regions of low energy centered at ∼(7.5 Å, 0°) and (7.5 Å, 120°), which correspond to D/A pair geometries identified as on top and hollow in previous gas-phase and mean-field studies1,2,32,49,64 [examples of which are shown in Fig. 4(a)]. The approximate percentage of sampled D/A pairs corresponding to these geometries [shown in Fig. 4(a)] is calculated using D/A pairs where R < 8.5 Å and θ < 38° for on top and D/A pairs where R < 9 Å and 95° < θ < 160° for hollow.

However, the majority of the D/A pairs in the PMF do not correspond to either on top or hollow. Instead, these pairs correspond to a geometry, noted as the edge, that was recently identified using a condensed-phase analysis.1 While most of the edge pairs fall within the range of 10 Å < R < 14 Å and 40° < θ < 100°, they have a large variance in R and θ [one example structure is shown in Fig. 4(a)]. The edge geometry consists of pairs where only the edge of a SubPC arm is close to the C60 molecule. Most of the edge pairs fall within the following range on the PMF: 10 Å < R < 14 Å and 40° < θ < 100°.

Among the ensemble structures from each basin (e.g., on-top, hollow, or edge ensembles), the most probable structure is selected as the representative D/A pair geometry for analysis in Module 2.

Electronic structure calculations are performed on the interfacial D/A geometries selected from Module 2. The important parameters are (1) the charge of the D/A pair’s donor molecule (QD), (2) the energy for the D/A pair while isolated in the gas-phase (Egas), and (3) the pair’s oscillator strength (OS).50 These three parameters are shown by the x axis, y axis, and color bar, respectively, in Fig. 5.

FIG. 5.

Scatterplot of the excitation energy in the gas phase (Egas) vs the charge of the donor molecule (QD) for the 25 lowest excited states of each representative geometry. Dots are colored according to their oscillator strength (OS). States used for calculations in this paper are shown by large squares and labeled with their names, while other states are denoted by smaller circles.

FIG. 5.

Scatterplot of the excitation energy in the gas phase (Egas) vs the charge of the donor molecule (QD) for the 25 lowest excited states of each representative geometry. Dots are colored according to their oscillator strength (OS). States used for calculations in this paper are shown by large squares and labeled with their names, while other states are denoted by smaller circles.

Close modal

Transitions studied here are from the lowest-energy bright EX state to each of the dark CT states, with one exception. In the on-top geometry, a bright CT state is considered as a donor state in addition to the lowest-energy EX state and also as an acceptor state. The OS of the bright states for the on-top geometry is about half that of the EX1 state for the hollow and edge geometries, as shown in Fig. 5 and listed in Table S2. Hollow and edge each also have another EX state that is similar in both OS and Egas but is not used for these results. However, the dCT states for on top range from approximately equal to slightly lower in energy than its bright states, while all the dCT states for hollow and edge are much higher in energy than the bright states. The exact values shown in Fig. 5 are recorded in Table S2. EX and CT states as labeled in Fig. 5 are used for analysis in Module 4.

Next, we calculate the electronic-population-transfer rate constant, kM, and CT rate constant, kC, for transitions from a donor state (EX1 or bCT1) to an acceptor state (bCT1, dCT1, dCT2, or dCT3). The required inputs to calculate kM, as shown in Eq. (1), are the D/A energy-gap first and second moments, the coupling coefficient, ΓDA, and the excitation-energy correction, Wα (see Table I and Table S3). The relative energies between states are changed in MD simulations due to the molecular environment. Clearly, the CT states are expected to be stabilized by the condensed-phase polarizable environment more than the localized excitations.26 The rate constant, kC (see Table I), is the product of kM and ΔQD, the difference in charge of the donor molecule between the donor and acceptor states.

TABLE I.

Interfacial charge-transfer rates for representative geometries.

GeometryTransitionΓDA (meV)ΔQD (e)〈U〉 (meV)kC (nA)
On-top EX1 → dCT1 4.03 0.706 217 ± 2 180 ± 10 
EX1 → dCT2 24.46 0.616 148 ± 11 520 ± 70 
EX1 → dCT3 5.47 0.825 −168 ± 1 220 ± 10 
EX1 → bCT1 25.82 0.314 114 ± 5 510 ± 40 
bCT1 → dCT1 74.16 0.391 123 ± 2 24 ± 2 
bCT1 → dCT2 72.27 0.302 61.6 ± 0.1 1010 ± 20 
bCT1 → dCT3 21.13 0.511 −240 ± 12 0.014 ± 0.008 
Hollow EX1 → dCT1 1.85 0.943 −743 ± 14 0.006 ± 0.002 
EX1 → dCT2 20.21 0.959 −837 ± 14 0.04 ± 0.02 
EX1→ dCT3 15.53 0.905 −936 ± 30 0.006 ± 0.004 
Edge EX1 → dCT1 10.30 0.977 −354 ± 23 65 ± 5 
EX1 → dCT2 14.02 0.781 −432 ± 15 0.32 ± 0.07 
EX1 → dCT3 17.22 0.953 −481 ± 20 3.1 ± 0.6 
GeometryTransitionΓDA (meV)ΔQD (e)〈U〉 (meV)kC (nA)
On-top EX1 → dCT1 4.03 0.706 217 ± 2 180 ± 10 
EX1 → dCT2 24.46 0.616 148 ± 11 520 ± 70 
EX1 → dCT3 5.47 0.825 −168 ± 1 220 ± 10 
EX1 → bCT1 25.82 0.314 114 ± 5 510 ± 40 
bCT1 → dCT1 74.16 0.391 123 ± 2 24 ± 2 
bCT1 → dCT2 72.27 0.302 61.6 ± 0.1 1010 ± 20 
bCT1 → dCT3 21.13 0.511 −240 ± 12 0.014 ± 0.008 
Hollow EX1 → dCT1 1.85 0.943 −743 ± 14 0.006 ± 0.002 
EX1 → dCT2 20.21 0.959 −837 ± 14 0.04 ± 0.02 
EX1→ dCT3 15.53 0.905 −936 ± 30 0.006 ± 0.004 
Edge EX1 → dCT1 10.30 0.977 −354 ± 23 65 ± 5 
EX1 → dCT2 14.02 0.781 −432 ± 15 0.32 ± 0.07 
EX1 → dCT3 17.22 0.953 −481 ± 20 3.1 ± 0.6 

In the on-top geometry, the ΓDA values for the bCT1 → dCT1 and bCT1 → dCT2 transitions are at least double the size of any other ΓDA values considered. On the other hand, as transitions involving bCT1 use a CT state as the donor state, the ΔQD value for these transitions is significantly less than that of transitions with an EX donor state. As a result, the values of kC for transitions with either EX1 or bCT1 as the donor state in the on-top geometry are within an order of magnitude of each other. Additionally, the values of kM are comparable to those from a non-condensed-phase analysis based on optimized geometries.2 

For both the hollow and edge geometries, smaller ΓDA values for transitions from EX1 to dCT states than those in the on-top geometry lead to smaller rate constants. However, the kC values for the edge geometry are larger overall than those of the hollow geometry. This difference between the edge and hollow can be traced back to the fact that the transitions in the edge geometry correspond to a much smaller U value than those of the hollow geometry (see Table I). The kM values for transitions in the hollow geometry are observed to be significantly smaller than those previously reported in a non-condensed-phase analysis based on optimized geometries (see Table S4).2 Marcus theory parameters are available for all transitions in the supplementary material (Table S5).

We have described the software package ctramer for the analysis of CT rates based on electronic structure calculations, MD simulations, and rate theory. ctramer is a unique combination of well-established methods from different disciplines, which allows for a more high-quality study of photoinduced CT processes with explicit description of the molecular environment. The customizable features, software architecture, and guidelines for usage were discussed. Additionally, the scientific justification behind the different approaches, as well as example results, were provided.

ctramer will continue to be actively developed and supported, as it will remain a long-term focus of the authors. Additions planned for the immediate future are automated tuning of the parameters for electronic structure calculations, enabling different levels of theory, and machine-learning clustering methods for selecting representative structures. Other future goals include various extensions to improve the accessibility and computational efficiency of the software.

See the supplementary material for additional information regarding the results: parameters used for MD not included in GAFF, a table of excited state properties, energy correction terms for each state, information on non-Gaussian energy gap distributions, and Marcus theory parameters.

E.G., B.D.D., and M.S.C. acknowledge support from the Department of Energy (DOE), Basic Energy Sciences through the Chemical Sciences, Geosciences and Biosciences Division (Grant No. DE-SC0016501). X.S. acknowledges support from the National Natural Science Foundation of China (Grant No. 21903054) and the Program for Eastern Young Scholar at Shanghai Institutions of Higher Learning. Computational resources are provided by the National Energy Research Scientific Computing Center (NERSC), a U.S. Department of Energy Office of Science User Facility operated under Contract No. DE-AC02-05CH11231, and the uHPC cluster managed by the University of Houston and acquired through NSF Award No. OAC 1531814.

The data that support the findings of this study are openly available in GitHub.40 

1.
J.
Tinnin
,
S.
Bhandari
,
P.
Zhang
,
H.
Aksu
,
B.
Maiti
,
E.
Geva
,
B. D.
Dunietz
,
X.
Sun
, and
M. S.
Cheung
,
Phys. Rev. Appl.
13
(
5
),
054075
(
2020
).
2.
M. H.
Lee
,
E.
Geva
, and
B. D.
Dunietz
,
J. Phys. Chem. C
118
(
18
),
9780
(
2014
).
3.
Y.
Zhao
and
W.
Liang
,
Chem. Soc. Rev.
41
(
3
),
1075
(
2012
).
4.
J.
Kim
and
S.
Yim
,
Appl. Phys. Lett.
99
(
19
),
193303
(
2011
).
5.
H.
Tian
,
Z.
Yu
,
A.
Hagfeldt
,
L.
Kloo
, and
L.
Sun
,
J. Am. Chem. Soc.
133
(
24
),
9413
(
2011
).
6.
X.
Tong
,
B. E.
Lassiter
, and
S. R.
Forrest
,
Org. Electron.
11
(
4
),
705
(
2010
).
7.
A.
Hagfeldt
,
G.
Boschloo
,
L.
Sun
,
L.
Kloo
, and
H.
Pettersson
,
Chem. Rev.
110
(
11
),
6595
(
2010
).
8.
T.
Wang
and
J.-L.
Brédas
,
Matter
2
(
1
),
119
(
2020
).
9.
J.
Kurpiers
,
T.
Ferron
,
S.
Roland
,
M.
Jakoby
,
T.
Thiede
,
F.
Jaiser
,
S.
Albrecht
,
S.
Janietz
,
B. A.
Collins
,
I. A.
Howard
, and
D.
Neher
,
Nat. Commun.
9
(
1
),
2038
(
2018
).
10.
S. M.
Falke
,
C. A.
Rozzi
,
D.
Brida
,
M.
Maiuri
,
M.
Amato
,
E.
Sommer
,
A.
De Sio
,
A.
Rubio
,
G.
Cerullo
,
E.
Molinari
, and
C.
Lienau
,
Science
344
(
6187
),
1001
(
2014
).
11.
L.
Benatto
,
C. F. N.
Marchiori
,
T.
Talka
,
M.
Aramini
,
N. A. D.
Yamamoto
,
S.
Huotari
,
L. S.
Roman
, and
M.
Koehler
,
Thin Solid Films
697
,
137827
(
2020
).
12.
N. C.
Giebink
,
G. P.
Wiederrecht
,
M. R.
Wasielewski
, and
S. R.
Forrest
,
Phys. Rev. B
82
(
15
),
155305
(
2010
).
13.
N. C.
Giebink
,
G. P.
Wiederrecht
,
M. R.
Wasielewski
, and
S. R.
Forrest
,
Phys. Rev. B
83
(
19
),
195326
(
2011
).
14.
M. H.
Lee
,
B. D.
Dunietz
, and
E.
Geva
,
J. Phys. Chem. Lett.
5
(
21
),
3810
(
2014
).
15.
D.
Balamurugan
,
A. J. A.
Aquino
,
F.
de Dios
,
L.
Flores
,
H.
Lischka
, and
M. S.
Cheung
,
J. Phys. Chem. B
117
(
40
),
12065
(
2013
).
16.
J.-L.
Brédas
,
D.
Beljonne
,
V.
Coropceanu
, and
J.
Cornil
,
Chem. Rev.
104
(
11
),
4971
(
2004
).
17.
V.
D’Innocenzo
,
G.
Grancini
,
M. J. P.
Alcocer
,
A. R. S.
Kandada
,
S. D.
Stranks
,
M. M. M.
Lee
,
G.
Lanzani
,
H. J.
Snaith
, and
A.
Petrozza
,
Nat. Commun.
5
(
1
),
3586
(
2014
).
18.
C. W.
Tang
and
S. A.
VanSlyke
,
Appl. Phys. Lett.
51
,
913
(
1987
).
19.
X.
Sun
and
E.
Geva
,
J. Chem. Theory Comput.
12
(
6
),
2926
(
2016
).
20.
X.
Sun
,
P.
Zhang
,
Y.
Lai
,
K. L.
Williams
,
M. S.
Cheung
,
B. D.
Dunietz
, and
E.
Geva
,
J. Phys. Chem. C
122
(
21
),
11288
(
2018
).
21.
X.
Sun
and
E.
Geva
,
J. Phys. Chem. A
120
(
19
),
2976
(
2016
).
22.
Z.
Hu
,
Z.
Tong
,
M. S.
Cheung
,
B. D.
Dunietz
,
E.
Geva
, and
X.
Sun
,
J. Phys. Chem. B
124
(
43
),
9579
(
2020
).
23.
X.
Sun
and
E.
Geva
,
J. Chem. Phys.
144
(
24
),
244105
(
2016
).
24.
X.
Sun
and
E.
Geva
,
J. Chem. Phys.
145
(
6
),
064109
(
2016
).
25.
Z.
Tong
,
X.
Gao
,
M. S.
Cheung
,
B. D.
Dunietz
,
E.
Geva
, and
X.
Sun
,
J. Chem. Phys.
153
(
4
),
044105
(
2020
).
26.
A. K.
Manna
,
D.
Balamurugan
,
M. S.
Cheung
, and
B. D.
Dunietz
,
J. Phys. Chem. Lett.
6
(
7
),
1231
(
2015
).
27.
H.
Aksu
,
A.
Schubert
,
E.
Geva
, and
B. D.
Dunietz
,
J. Phys. Chem. B
123
(
42
),
8970
(
2019
).
28.
M. A.
Steffen
,
K.
Lao
, and
S. G.
Boxer
,
Science
264
(
5160
),
810
(
1994
).
29.
A.
Niedringhaus
,
V. R.
Policht
,
R.
Sechrist
,
A.
Konar
,
P. D.
Laible
,
D. F.
Bocian
,
D.
Holten
,
C.
Kirmaier
, and
J. P.
Ogilvie
,
Proc. Natl. Acad. Sci. U. S. A.
115
(
14
),
3563
(
2018
).
30.
J.
Koepke
,
E.-M.
Krammer
,
A. R.
Klingen
,
P.
Sebban
,
G. M.
Ullmann
, and
G.
Fritzsch
,
J. Mol. Biol.
371
(
2
),
396
(
2007
).
31.
K. L.
Mutolo
,
E. I.
Mayo
,
B. P.
Rand
,
S. R.
Forrest
, and
M. E.
Thompson
,
J. Am. Chem. Soc.
128
(
25
),
8108
(
2006
).
32.
N.
Ilyas
,
S. S.
Harivyasi
,
P.
Zahl
,
R.
Cortes
,
O. T.
Hofmann
,
P.
Sutter
,
E.
Zojer
, and
O. L. A.
Monti
,
J. Phys. Chem. C
120
(
13
),
7113
(
2016
).
33.
G.
D’Avino
,
L.
Muccioli
,
F.
Castet
,
C.
Poelking
,
D.
Andrienko
,
Z. G.
Soos
,
J.
Cornil
, and
D.
Beljonne
,
J. Phys.: Condens. Matter
28
(
43
),
433002
(
2016
).
34.
D. S.
Josey
,
S. R.
Nyikos
,
R. K.
Garner
,
A.
Dovijarski
,
J. S.
Castrucci
,
J. M.
Wang
,
G. J.
Evans
, and
T. P.
Bender
,
ACS Energy Lett.
2
(
3
),
726
(
2017
).
35.
R.
Pandey
,
A. A.
Gunawan
,
K. A.
Mkhoyan
, and
R. J.
Holmes
,
Adv. Funct. Mater.
22
(
3
),
617
(
2012
).
36.
D. E.
Wilcox
,
M. H.
Lee
,
M. E.
Sykes
,
A.
Niedringhaus
,
E.
Geva
,
B. D.
Dunietz
,
M.
Shtein
, and
J. P.
Ogilvie
,
J. Phys. Chem. Lett.
6
(
3
),
569
(
2015
).
37.
M. H.
Lee
,
E.
Geva
, and
B. D.
Dunietz
,
J. Phys. Chem. A
120
(
19
),
2970
(
2016
).
38.
Y.
Shao
,
L. F.
Molnar
,
Y.
Jung
,
J.
Kussmann
,
C.
Ochsenfeld
,
S. T.
Brown
,
A. T. B.
Gilbert
,
L. V.
Slipchenko
,
S. V.
Levchenko
,
D. P.
O’Neill
,
R. A.
DiStasio
, Jr.
,
R. C.
Lochan
,
T.
Wang
,
G. J. O.
Beran
,
N. A.
Besley
,
J. M.
Herbert
,
C.
Yeh Lin
,
T.
Van Voorhis
,
S.
Hung Chien
,
A.
Sodt
,
R. P.
Steele
,
V. A.
Rassolov
,
P. E.
Maslen
,
P. P.
Korambath
,
R. D.
Adamson
,
B.
Austin
,
J.
Baker
,
E. F. C.
Byrd
,
H.
Dachsel
,
R. J.
Doerksen
,
A.
Dreuw
,
B. D.
Dunietz
,
A. D.
Dutoi
,
T. R.
Furlani
,
S. R.
Gwaltney
,
A.
Heyden
,
S.
Hirata
,
C.-P.
Hsu
,
G.
Kedziora
,
R. Z.
Khalliulin
,
P.
Klunzinger
,
A. M.
Lee
,
M. S.
Lee
,
W.
Liang
,
I.
Lotan
,
N.
Nair
,
B.
Peters
,
E. I.
Proynov
,
P. A.
Pieniazek
,
Y.
Min Rhee
,
J.
Ritchie
,
E.
Rosta
,
C.
David Sherrill
,
A. C.
Simmonett
,
J. E.
Subotnik
,
H.
Lee Woodcock
 III
,
W.
Zhang
,
A. T.
Bell
,
A. K.
Chakraborty
,
D. M.
Chipman
,
F. J.
Keil
,
A.
Warshel
,
W. J.
Hehre
,
H. F.
Schaefer
 III
,
J.
Kong
,
A. I.
Krylov
,
P. M. W.
Gill
, and
M.
Head-Gordon
,
Phys. Chem. Chem. Phys.
8
(
27
),
3172
(
2006
).
39.
D. A.
Case
,
T. A.
Darden
,
I. T. E.
Cheatham
,
C. L.
Simmerling
,
J.
Wang
,
R. E.
Duke
,
R.
Luo
,
R. C.
Walker
,
W.
Zhang
,
K. M.
Merz
,
B.
Roberts
,
S.
Hayik
,
A.
Roitberg
,
G.
Seabra
,
J.
Swails
,
A. W.
Götz
,
I.
Kolossváry
,
K. F.
Wong
,
F.
Paesani
,
J.
Vanicek
,
R. M.
Wolf
,
J.
Liu
,
X.
Wu
,
S. R.
Brozell
,
T.
Steinbrecher
,
H.
Gohlke
,
Q.
Cai
,
X.
Ye
,
J.
Wang
,
M.-J.
Hsieh
,
G.
Cui
,
D. R.
Roe
,
D. H.
Mathews
,
M. G.
Seetin
,
R.
Salomon-Ferrer
,
C.
Sagui
,
V.
Babin
,
T.
Luchko
,
S.
Gusarov
,
A.
Kovalenko
, and
P. A.
Kollman
,
AMBER12
(
University of California
,
San Francisco
,
2012
).
40.
J.
Tinnin
,
H.
Aksu
,
Z.
Tong
,
P.
Zhang
,
E.
Geva
,
B. D.
Dunietz
,
X.
Sun
, and
M. S.
Cheung
,
2021
, ctramer SubPC/C60 Data, Github, https://github.com/ctramer/ctramer.
41.
E.
Runge
and
E. K. U.
Gross
,
Phys. Rev. Lett.
52
(
12
),
997
(
1984
).
42.
K.
Begam
,
S.
Bhandari
,
B.
Maiti
, and
B. D.
Dunietz
,
J. Chem. Theory Comput.
16
(
5
),
3287
(
2020
).
43.
S.
Bhandari
and
B. D.
Dunietz
,
J. Chem. Theory Comput.
15
(
8
),
4305
(
2019
).
44.
M. H.
Lee
,
B. D.
Dunietz
, and
E.
Geva
,
J. Phys. Chem. C
117
(
44
),
23391
(
2013
).
45.
A. A.
Voityuk
and
N.
Rösch
,
J. Chem. Phys.
117
(
12
),
5607
(
2002
).
46.
V. A.
Rassolov
,
J. A.
Pople
,
M. A.
Ratner
, and
T. L.
Windus
,
J. Chem. Phys.
109
(
4
),
1223
(
1998
).
47.
R.
Baer
and
D.
Neuhauser
,
Phys. Rev. Lett.
94
(
4
),
043002
(
2005
).
48.
E.
Livshits
and
R.
Baer
,
Phys. Chem. Chem. Phys.
9
(
23
),
2932
(
2007
).
49.
T.
Stein
,
L.
Kronik
, and
R.
Baer
,
J. Am. Chem. Soc.
131
(
8
),
002818
2820
(
2009
).
50.
J. W.
Robinson
,
Atomic Spectroscopy
, 2nd ed. (
Marcel Dekker
,
New York
,
1996
).
51.
L.
Martínez
,
R.
Andrade
,
E. G.
Birgin
, and
J. M.
Martínez
,
J. Comput. Chem.
30
(
13
),
2157
(
2009
).
52.
J.
Wang
,
W.
Wang
,
P. A.
Kollman
, and
D. A.
Case
,
J. Mol. Graphics
25
(
2
),
247
(
2006
).
53.
J.
Wang
,
R. M.
Wolf
,
J. W.
Caldwell
,
P. A.
Kollman
, and
D. A.
Case
,
J. Comput. Chem.
25
(
9
),
1157
(
2004
).
54.
C. G.
Claessens
,
D.
González-Rodríguez
, and
T.
Torres
,
Chem. Rev.
102
(
3
),
835
(
2002
).
55.
X.
Wu
,
Z.
Liu
,
S.
Huang
, and
W.
Wang
,
Phys. Chem. Chem. Phys.
7
(
14
),
2771
(
2005
).
56.
G.
Ciccotti
and
J. P.
Ryckaert
,
Comput. Phys. Rep.
4
(
6
),
346
(
1986
).
57.
T.
Darden
,
D.
York
, and
L.
Pedersen
,
J. Chem. Phys.
98
(
12
),
10089
(
1993
).
58.
D. M.
York
,
T. A.
Darden
, and
L. G.
Pedersen
,
J. Chem. Phys.
99
(
10
),
8345
(
1993
).
59.
J.
Han
,
P.
Zhang
,
H.
Aksu
,
B.
Maiti
,
X.
Sun
,
E.
Geva
,
B. D.
Dunietz
, and
M. S.
Cheung
,
J. Chem. Theory Comput.
16
(
10
),
6481
(
2020
).
60.
P. F.
Barbara
,
T. J.
Meyer
, and
M. A.
Ratner
,
J. Phys. Chem.
100
(
31
),
13148
(
1996
).
61.
R. A.
Marcus
,
J. Chem. Phys.
24
(
5
),
966
(
1956
).
62.
R. A.
Marcus
,
J. Chem. Phys.
24
(
5
),
979
(
1956
).
63.
R. A.
Marcus
,
Rev. Mod. Phys.
65
(
3
),
599
(
1993
).
64.
C.
Xu
,
G.
Li
, and
S.
Zheng
,
J. Photochem. Photobiol., A
403
,
112852
(
2020
).

Supplementary Material