DOI: 10.1039/D4CP01568A
(Paper)
Phys. Chem. Chem. Phys., 2024, Advance Article

Yanyan Yang,
Qianqian Jin and
Shiwei Yin*

Key Laboratory for Macromolecular Science of Shaanxi Province, School of Chemistry and Chemical Engineering, Shaanxi Normal University, Xi’an City 710119, People's Republic of China. E-mail: yin_sw@snnu.edu.cn

Received
17th April 2024
, Accepted 5th August 2024

First published on 8th August 2024

For planar and rigid π-conjugated molecular systems, electrostatic and inductive interactions are pivotal in governing molecular packing structures and electron polarization energies. These electrostatic interactions typically exhibit an anisotropic nature within π-conjugated systems. In this study, we utilize the atoms in molecules (AIM) theory in conjunction with linear response theory to decompose molecular polarizability into distributed atomic polarizability tensors. On the basis of atomic polarizability tensors, we extended an anisotropic polarizable model into the AMOEBA polarizable force field. Both anisotropic and isotropic polarizable models in combination with various density functional theory (DFT)-derived atomic multipoles were applied to optimize the experimental crystals of naphthalene and anthracene. Furthermore, these two types of electrostatic models, coupled with the evolutionary algorithm USPEX program, are utilized to predict the crystal structures of oligoacenes. Our findings demonstrate that the anisotropic polarizable model exhibits superior performance in crystal refinement and crystal structure prediction. This enriched anisotropic polarizable model is seamlessly integrated into the AMOEBA polarizable force field and readily applicable within our modified Tinker program.

In traditional force fields, molecular electrostatic interactions typically employ atom-atom models such as DREIDING,^{7} CHARMM,^{8} COMPASS,^{9} among others, with atomic charges obtained from quantum mechanics (QM)-derived ESP fitting charges. Electrostatic interactions, based on atomic multipoles (AMPs) up to the quadrupole level, are utilized to optimize and rank candidate crystals. This anisotropic electrostatic potential demonstrates excellent performance in predicting molecular crystal structures,^{10–12} significantly impacting electronic couplings and primary transport channels for charge carrier transfer processes.^{13–15} Site energies and their disorder properties serve as key parameters in determining the charge mobilities of amorphous materials.^{16–19} Besides electrostatic interactions, inductive interactions play a crucial role in evaluating site energy and electron polarization energy.^{20–22} Additionally, external reorganization energy is also a critical parameter controlling the rate of charge transfer in organic solids. To explicitly evaluate this, a polarizable force field is necessary to estimate structure relaxation and electron polarization effects during charge transfer processes.^{23–25} Consequently, the development of an accurate polarizable force field, especially tailored for π-conjugated molecules, is imperative for modeling electron transfer processes within their crystalline or amorphous materials.

Currently, numerous polarizable models, such as the fluctuating charge model,^{26} the Drude model,^{27} and the induced dipole model,^{28} are extensively utilized in the literature to describe inductive interactions. Among these, the induced dipole model stands out as the most commonly employed approach for capturing polarization effects. Within linear response theory, the magnitude of induced dipoles is directly proportional to atomic polarizability (α) and the applied electric field. This electric field arises from contributions of both permanent electrostatic and induced dipoles to its surrounding atoms. Consequently, solving for induced dipoles necessitates self-consistent field (SCF) calculations, employing methods such as the Applequist scheme^{29} or Thole's modified model to circumvent the polarization catastrophe resulting from closely situated induced dipoles.^{30,31} This approach ensures the consideration of many-body effects of electron polarization. However, the accurate determination of atomic polarizability is pivotal in evaluating induced dipole moments and essential for the development of accurate polarizable force fields.

There are two widely utilized approaches to determine atomic polarizability. One approach involves fitting experimental or QM-calculated molecular polarizabilities within a training set to derive atomic polarizabilities.^{32,33} The molecular polarizability tensor is obtained from the Applequist scheme or Thole's modified induced dipole model on the basis of isotropic atomic polarizability where the intramolecular induced dipole–dipole interactions need to be considered. This method offers transferable and isotropic atomic polarizabilities, effectively reproducing polarizability tensors for molecules with similar structures to those in the training set. However, if the molecule significantly deviates from the training set, atomic polarizability may prove unsuitable. To address this, new types of molecules should be incorporated into the training set to develop novel atomic polarizabilities. For instance, benzene, naphthalene, and anthracene were added to the training set to establish new aromatic carbon and hydrogen atomic polarizabilities in AMOEBA09.^{34} Another approach involves QM calculations coupled with linear response theory to numerically determine distributed atomic polarizabilities.^{35–38} Although the atomic polarizabilities obtained via QM-based methods are not transferable and depend on the calculation methods, this approach offers several advantages. First, it is based on first principles and is independent of experimental data or training set selection. Secondly, it yields molecule-specific atomic polarizabilities, enabling direct capture of anisotropy from QM calculations, potentially enhancing the description of inductive interactions for specific molecules. Furthermore, QM-based methods can provide state-specific atomic polarizabilities or electrostatic potentials when the molecule is in a specific electronic state, such as an excited or charged state,^{23,39} facilitating explicit assessment of electron polarizations or solvatochromic effects in condensed phases. In addition, QM-derived atomic polarizabilities inherently exist in tensor form, allowing for natural inclusion of anisotropic characteristics.^{35,40}

Most polarizable force fields, such as DRF90,^{41} AMOBEA^{42} and AMBER ff02,^{43} predominantly utilize induced dipole models and isotropic atomic polarizability. Models based on polarizability tensors are rarely reported.^{44,45} Recently, Head-Gordon developed an anisotropic polarizable model (aniso-AMOEBA) for water, utilizing the atomic polarizability tensor of water.^{45} Compared to the isotropic polarizable model AMOEBA03,^{42} aniso-AMOEBA exhibits superior performance in many-body polarization energies, as evidenced by energy decomposition analyses of high-accuracy DFT calculations of potential energy surfaces for water clusters or anion-water trimers.^{45} This suggests that accurate anisotropic atomic polarizabilities indeed improve inductive interactions. In these studies,^{44,45} the atomic polarizability tensors were defined using the Williams–Stone–Misquitta method^{40,46} implemented through Misquitta and Stone's CamCASP program. After that, certain programs such as Orient and DMACRYS need to be applied to transform these spherical tensors into atomic polarizability tensors within the global coordinate frames. Apart from the Williams–Stone–Misquitta method, Macchi and colleagues utilized atoms in molecules (AIM) theory to partition electron density and calculated static atomic polarizability tensors under linear response theory.^{38,47} Using AIM theory coupled with linear response theory, we derived state-specific atomic polarizabilities for oligoacenes and successfully evaluated polarization energy and external reorganization energy during charge transfer processes in organic solids.^{23,25,48} However, anisotropic polarizable models are rarely applied in standard full-atom force fields. Hence, there is a need to develop a more general-purpose anisotropic polarizable model for full-atom force fields to better describe intermolecular interactions.

In this study, we computationally derive the distributed atomic polarizability tensor using AIM theory combined with linear response theory. Utilizing this atomic polarizability tensor, we develop an anisotropic polarizable model and extend the anisotropic AMOEBA polarizable force field to aromatic hydrocarbons, implemented through our modified Tinker program.^{49} As an application, we examine the performance of both anisotropic and isotropic polarizable models on optimizing and predicting crystal structures of rigid aromatic oligoacenes, such as naphthalene, anthracene, tetracene, and pentacene. Our comparative analysis reveals that the more accurate polarizable model, specifically the anisotropic polarizable model, demonstrates superior performance in crystal structure prediction.

μ^{ind}_{i} = α_{i}(F^{perm}_{i} + F^{ind}_{i})
| (1) |

(2a) |

(2b) |

According to eqn (2a) and (2b), eqn (1) can be rewritten as follows

(3) |

To solve for the induced dipole of the i-th atom, we first need to know the induced dipoles of the other atoms. In eqn (3), the direct electric field (the first term on the right side of equation) from the permanent multipoles is fixed, and the mutual electric field (the second term) from the neighboring induced dipoles of the j atoms is updated until the RMS of the induced dipoles is smaller than 0.00001 Debye. Thus, eqn (3) is solved by self-consistent field (SCF) iteration, thereby accounting for the many-body polarization effect.

2.1.1. Isotropic polarizable model for oligoacene. For the isotropic polarizable model, α_{i} adopts the isotropic atomic polarizability in eqn (1) and (3). To account for the anisotropy of molecular polarizability, the induced dipole–dipole interactions from intramolecular dipoles must be included on the basis of Thole^{31} or Applequist^{29} methods by setting the keyword mutual-11-scale as 1.0 in the AMOEBA force field. According to Ponder and Ren's suggestion, the isotropic atomic polarizabilities of aromatic C and H are set as 1.75 and 0.696 4πε_{0} Å^{3}, respectively, which can well describe the polarizability tensors of naphthalene and anthracene molecules.^{34}

2.1.2. Anisotropic polarizable model. When α_{i} adopts the second-order tensor form, that is,

(4) |

Incorporating this into eqn (1), we establish the anisotropic polarizable model. Since the distributed atomic polarizability tensor is derived through molecular quantum mechanics (QM) calculations coupled with a finite field approach, the contribution from intramolecular induced dipole–dipole interaction is already accounted for during the parametrization of atomic polarizability. Therefore, in the anisotropic polarizable model, the intramolecular dipole–dipole interactions should be disregarded by setting the keyword mutual-11-scale as 0.0, and the molecular polarizability tensor should be directly obtained through the summation of distributed atomic polarizability tensors.

2.4.1. Optimization and comparison of crystal structures. The experimental crystal structures of naphthalene (NAPHTA06),^{55} anthracene (ANTCEN01),^{56} tetracene (TETCEN)^{57} and pentacene (PENCEN)^{58} are refined by the xtalmin program without constraining their space groups on the basis of different force fields shown in Table 1. Different from rigid optimization according to rigid body force and torques performed by the DMACRYS program,^{59,60} xtalmin optimized the crystal according to the forces of all atoms and lattice parameters. The forces from inductive interactions are directly included according to the SCF induced dipoles from eqn (3). The limited memory BFGS Quasi-Newton method^{61} is applied to optimize the structure and the convergence criterion is set as 0.01 kcal mol^{−1} Å^{−1} for RMS gradient per atom. For most situations, it takes 2–3 minutes to finish the optimization of a crystal structure using Xeon(R) Gold 6128 CPU at 3.40 GHz. The refined crystals and experimental structures can be compared from two perspectives: crystal similarity and packing similarity. The former is based on atomic positions in crystals to determine similarities of crystals structures. The typical method defines t-RMSD of atoms positions to compare crystal similarity. The t-RMSD is expressed as^{62}

where r_{i} are fractional coordinates of atoms in the crystal structure, G is the transformation matrix from fractional to Cartesian coordinates, N is the number of atoms in crystals, scripts exp/opt are the experimental and refined crystal structures. Here all the H atoms are excluded for comparison because it is hard to detect hydrogen positions using an X-ray diffraction method. Molecular packing similarly compares the similarity of molecular clusters which are generated with the central molecules and their surrounding neighbors in individual crystal structures. Another typical packing similarity (RMSD_{n}) is calculated on the basis of the deviation of centroids of overlapping molecules from each other,^{63,64} subscript n is the number of neighboring molecules in clusters. In this article, we set n as 20, which means 21 molecules are included in molecular clusters. The RMSD_{20} is defined as

where K the number of atoms in molecular clusters, d_{i} is the distance between the closest atoms. Both crystal similarity and packing similarity in this work are calculated using the CrystalCMP program.^{64} Our definition of RMSD_{20} differs slightly from the RMSD_{20} definition in Mercury,^{65} where the comparison involves clusters of 20 molecules.

(5) |

(6) |

Potentials | Electrostatic | vdW | Polarization |
---|---|---|---|

NON-polar | AMP | Buffered 14-7 | None |

ISO-polar | AMP | Buffered 14-7 | Isotropic polarizability |

ANISO-polar | AMP | Buffered 14-7 | Polarizability tensors by eqn (4) |

2.4.2. Crystal structure prediction. The molecular crystal structures are generated using the evolutionary algorithm USPEX code (version10.2),^{66–68} which is coupled with the xtalmin crystal optimizer in the Tinker package. The detailed procedures of crystal structure predictions, please see our previous study.^{12} Three force fields including the non-polarizable AMP model (non-polar), isotropic atomic polarizability-based polarizable model (ISO-polar) and atomic polarizability tensor-based polarizable model (ANISO-polar) are applied to optimize and rank the candidates produced by the USPEX program. The interaction potentials of the three models are summarized in Table 1. The difference between ANISO-polar and ISO-polar models primarily involves two aspects i.e. atomic polarizability (α_{i}) and the definition of neighboring atoms of the i atom in eqn (3). For the ISO-polar model, α_{i} is the isotropic number, and neighboring atoms j’ differ from the i atom and include both intermolecular and intramolecular atoms of the i atom. For the ANISO-polar model, α_{i} is adopted in the tensor form as shown in eqn (4), with neighboring atoms j’ representing only the intermolecular atoms of the i atom.

2.5.1. Parameterization and transformation of atomic polarizability tensors. According to the linear response theory, the static atomic polarizability tensor elements α_{i,xy} can be expressed by the following numerical gradient

where μ^{F/0} is the distributed atomic dipole moment under the electric field F and zero field, respectively. The subscript i is the atomic index. And x, y are the directions of Cartesian axes. F_{y} is the small applied electric field along the y direction. In QTAIM theory,^{69} there have been several proposals to decompose the total molecular polarizability into atomic contributions. This could be obtained by partitioning either the energy or the electron density distribution in R^{3} or in Hilbert space. Stone and his coworker analyzed several methods to partition the molecular polarizability and they concluded that space-partitioned atomic polarizability volumes would be more efficient and stable.^{36,70} Keith later generalized a hard space partitioning of molecular polarizabilities and coded it in the AIMALL program.^{71} In order to get the stable numerical gradient, we set the electric field F to 0.02 a.u. According to the space partition,^{69} the atomic dipole moments can be obtained from QM-calculated electron density with/without applied electric fields by the AIMALL program.^{71} As atomic polarizability, like atomic multipoles, possibly depends on functionals, we adopt the same functional to parameterize both the atomic polarizability tensors and the atomic multipoles. The 6-311+G(d,p) basis set is applied to parameterize the atomic polarizability.

where m_{a} is the atomic mass and a represents the atoms of this molecule.

where C is the 3 × 3 modal matrix, C^{−1} is the inverse matrix of C, I_{1/2/3} represent three principal moments, respectively. And the indices 1, 2 and 3 correspond to the x, y and z axes in our defined molecular standard coordinates.

where α^{st}_{i} represents the atomic polarizability tensor under the molecular standard coordinates. And these global atomic polarizability tensors are directly applied in eqn (3) to compute the induced dipoles.

(7) |

Additionally, atomic polarizability tensors are dependent on the chosen coordinate system, and they need to undergo coordinate conversion when transitioning between different frames of reference. For simplicity, we first define a molecular standard coordinate system to parameterize atomic polarizability. In this frame, the centroid serves as the origin, with x, y, and z aligned along the directions of principal moments of inertia, ordered from the smallest to the largest. For instance, in the case of the naphthalene molecule, x, y, and z correspond respectively to the molecular long (x), short (y), and normal (z) axes, as illustrated in Fig. 1. The atomic polarizability tensors under the standard coordinates can be transformed into global tensors (α^{global}_{i}) in real crystal coordinates using a linear transformation matrix which can diagonalize the moment of the inertia tensor of molecules. The detailed procedure for transforming atomic polarizability tensors is outlined as follows.

Fig. 1 The molecular standard orientation. x, y and z axes are aligned with the molecular long, short and normal axes. |

(1) Moving the centroids of each molecule in crystals to the origin point.

(2) Building the molecular moment of the inertia tensor,

(8) |

(3) Diagonalizing the inertia tensor I matrix to get the linear transformation matrix C, that is,

(9) |

(4) Subsequently, the global atomic polarizability tensors under the global coordinates are derived using the following relation:

α^{global}_{i} = Cα^{st}_{i}C^{−1}
| (10) |

(5) When the molecular Cartesian coordinates are updated during structure optimization, we repeat steps 1 through 4 to obtain updated global atomic polarizability tensors.

It is worth mentioning that if the molecule is flexible and has a lot different conformers, then transformation matrix C would be not unique. In such cases, it is necessary to consider the geometry-dependent atomic polarizability tensors and to utilize local reference coordinates, similar to the coordinate conversion used for atomic multipole moments in AMOBEA. This topic is beyond the scope of this article. Thus, our anisotropic polarizable model is just limited to rigid molecules and not applied to flexible molecules at the stage. In contrast to the rigid molecule optimization performed by the DMACRYS program using anisotropic atomic polarizabilities, the optimization in xtalmin uses both inter and intra-molecular forces to optimize the positions of all atoms and lattice parameters, with the atomic multipoles and polarizabilities fixed. This is a reasonable assumption for these essentially rigid molecules.

2.5.2. QM-derived atomic multipoles (AMPs). The distributed multipole analysis of Gaussian wavefunctions is carried out by the GDMA version 2.2 program to get the distributed AMPs up to the quadrupole level under molecular standard coordinates.^{52} Using the poledit program in the Tinker package, these AMPs in standard coordinates can be transformed into local axes AMPs in the AMOEBA force field. The Tinker program can directly deal with the conversion of these AMPs from local axes to the global axes. As we stated before, electrostatic interaction is determined by AMPs which are related to QM-calculated molecular wavefunctions. To investigate the effect of DFT functionals on AMPs, the wavefunctions of molecules are separately calculated by various DFT functionals to obtain the functional-dependent AMP parameters. These functional-dependent AMPs are applied to refine the molecular crystals. All these QM calculations are performed by the G09 program at the 6-311G(d,p) level.^{72}

Type | α_{xx} |
α_{xy} |
α_{xz} |
α_{yz} |
α_{yy} |
α_{zz} |
α^{num}_{iso} |
---|---|---|---|---|---|---|---|

1-C | 1.711 | −0.248 | 0.000 | 0.000 | 1.422 | 0.866 | 1.333 |

2-C | 2.273 | 0.050 | 0.000 | 0.000 | 1.154 | 0.829 | 1.419 |

3-C | 2.366 | 0.005 | 0.000 | 0.000 | 1.561 | 0.585 | 1.504 |

4-C | 2.273 | −0.050 | 0.000 | 0.000 | 1.154 | 0.829 | 1.419 |

5-C | 1.711 | 0.248 | 0.000 | 0.000 | 1.422 | 0.866 | 1.332 |

6-H | 0.722 | 0.290 | 0.000 | 0.000 | 0.288 | 0.142 | 0.384 |

7-H | 0.147 | 0.073 | 0.000 | 0.000 | 0.741 | 0.140 | 0.343 |

8-H | 0.147 | −0.073 | 0.000 | 0.000 | 0.741 | 0.140 | 0.343 |

9-H | 0.722 | −0.290 | 0.000 | 0.000 | 0.288 | 0.142 | 0.384 |

10-C | 1.711 | −0.248 | 0.000 | 0.000 | 1.422 | 0.866 | 1.333 |

11-C | 2.273 | 0.050 | 0.000 | 0.000 | 1.154 | 0.829 | 1.419 |

12-C | 2.366 | 0.005 | 0.000 | 0.000 | 1.561 | 0.585 | 1.504 |

13-C | 2.273 | −0.050 | 0.000 | 0.000 | 1.154 | 0.829 | 1.419 |

14-C | 1.711 | 0.248 | 0.000 | 0.000 | 1.422 | 0.866 | 1.332 |

15-H | 0.722 | 0.290 | 0.000 | 0.000 | 0.288 | 0.142 | 0.384 |

16-H | 0.147 | 0.073 | 0.000 | 0.000 | 0.741 | 0.140 | 0.343 |

17-H | 0.147 | −0.073 | 0.000 | 0.000 | 0.741 | 0.140 | 0.343 |

18-H | 0.722 | −0.290 | 0.000 | 0.000 | 0.288 | 0.142 | 0.384 |

Σα | 24.14 | 0.0 | 0.0 | 0.0 | 17.54 | 9.08 | 16.92 |

Molecules | Method | α_{xx} |
α_{yy} |
α_{zz} |
α_{iso} |
k |
---|---|---|---|---|---|---|

Naphthalene | QM | 24.20 | 17.60 | 9.07 | 16.96 | 0.258 |

ANISO-polar | 24.14 | 17.54 | 9.08 | 16.92 | 0.258 | |

ISO-polar | 21.88 | 18.50 | 9.77 | 16.72 | 0.216 (−16.4%) | |

Anthracene | QM | 41.50 | 23.80 | 12.20 | 25.77 | 0.330 |

ANISO-polar | 41.37 | 23.80 | 12.00 | 25.72 | 0.332 | |

ISO-polar | 32.86 | 24.59 | 12.86 | 23.43 | 0.248 (−24.9%) | |

Tetracene | QM | 61.80 | 30.10 | 14.80 | 35.57 | 0.390 |

ANISO-polar | 61.82 | 30.06 | 14.81 | 35.56 | 0.390 | |

ISO-polar | 44.62 | 30.52 | 15.86 | 30.33 | 0.274 (−29.6%) | |

Pentacene | QM | 86.50 | 36.80 | 17.70 | 47.00 | 0.436 |

ANISO-polar | 86.58 | 36.82 | 17.69 | 47.03 | 0.437 | |

ISO-polar | 57.05 | 36.50 | 18.93 | 37.49 | 0.294 (−32.6%) |

As we known, the polarizability represents the ability of electron deformation. For aromatic oligoacenes, the delocalized π electrons determine electron deformation, predominantly localized on carbon (C) atoms and rarely on hydrogen (H) atoms. Consequently, the atomic polarizability of C should be significantly greater than that of H atoms. By examining the π-electron densities on various sites, we can deduce that the atomic polarizability of C atoms follows the order 1-C < 2-C < 3-C. Since π electrons primarily extend along the molecular long axis, compared to the short axis, with minimal deformation along the normal axis, the diagonal elements of C atoms exhibit the feature α_{xx} (long axis) > α_{yy} (short axis) ≫ α_{zz} (normal axis). Regarding the off-diagonal elements of C atoms, tensor elements involving the normal axis, such as α_{xz} and α_{yz}, are strictly zero, while α_{xy} is slightly larger than zero.

To assess the accuracy of atomic polarizability parameters, molecular polarizability is typically reproduced by its polarizable model. In Table 3, we present the molecular polarizabilities obtained from isotropic and anisotropic atomic polarizabilities. As mentioned earlier, ISO-polar considers intramolecular induced dipole–dipole interactions by setting intra-mutual-11 as 1.0 to represent the anisotropic molecular polarizability tensor, while ANISO-polar excludes these interactions by setting intra-mutual-11 as 0.0, summing only the distributed atomic polarizability tensors to obtain the molecular polarizability tensor. The isotropic molecular polarizability (α_{iso}) reflects the overall capability of molecular polarization. From Table 3, it is evident that the QM-calculated α_{iso} values are identical to those of the ANISO-polar model. However, ISO-polar fails to simultaneously reproduce these four molecular isotropic polarizabilities. In particular, while ISO-polar adequately represents naphthalene and anthracene, it notably underestimates the QM-calculated results for longer molecules like tetracene and pentacene. As observed previously by Soos, the static polarizability for each CH unit in conjugated polyacetylene significantly increases with the enhancement of π-electron conjugation.^{73} This indicates that the aromatic carbon (C) atomic polarizability strongly correlates with π-electron conjugation, rendering the uniform atomic polarizability, as in the ISO-polar model, unsuitable for aromatic molecules with varying degrees of π-electron conjugation. As mentioned earlier, the uniform aromatic C and H atomic polarizabilities in the AMOEBA force field are parameterized by fitting the molecular polarizability from a database that includes short benzene, naphthalene, and anthracene but excludes longer oligoacenes like tetracene and pentacene.^{34} Consequently, ISO-polar may slightly overestimate the α_{iso} of the short benzene and naphthalene molecules and underestimate the α_{iso} of the long tetracene and pentacene molecules.

In order to evaluate the degree of anisotropy, we define a normalized anisotropy parameter (k) for molecular polarizbility tensors as

(11) |

Fig. 2 Molecular packing similarity (RMSD_{20}) between the refined crystal structure of naphthalene and anthrance and their experimental crystal structures NAPTHA06 and ANTCEN01. Interaction potentials of NON-polar (orange), ISO-polar(green) and ANISO-polar (purple) are listed in Table 1. Atomic multipoles and atomic polarizability tensor for the ANISO-polar model are derived from the respective DFT functional calculations shown on abscissa. |

Considering the influence of functional-dependent atomic multipole moments (AMPs) on NON-polar models, it is observed that the smallest RMSD_{20} values for naphthalene (0.49 Å) and anthracene (0.49 Å) are derived from the double hybrid functional B2PLYP, while the largest RMSD_{20} values for naphthalene (0.52 Å) and anthracene (0.59 Å) are obtained from the meta-GGA M06-L functional. Fig. 3 reveals that despite the differences among the seven functionals, the QM-calculated electrostatic potentials are remarkably similar. The disparity (QM-AMP) between the QM-calculated electrostatic potentials and AMPs on the molecular surface appears negligible (white regions, close to 0.0 kcal mol^{−1}), indicating that AMPs accurately reproduce the QM-calculated molecular electrostatic potentials. Thus, different functional-derived AMPs exhibit comparable performance in optimization.

ISO-polar models (green bars) demonstrate smaller RMSD_{20} values and enhance the refined crystal structures compared to NON-polar models (pink bars) with varying functional-derived AMPs. This underscores the significant role of polarization interactions in optimizing crystal structures. Moreover, ANISO-polar models further improve the refined molecular structures due to their more accurate polarizable model, as evidenced by their ability to precisely reproduce QM-calculated molecular polarizability tensors (as shown in Table 3).

Considering the crystal similarity (t-RMSD) of refined structures of the three force fields depicted in Fig. S1 (ESI†), similar conclusions can be drawn from packing similarity (RMSD_{20}). In other words, ANISO-polar models exhibit the smallest crystal similarity (t-RMSD), NON-polar models demonstrate the largest t-RMSD values, while t-RMSD from ISO-polar models falls between them. Given the consistent results obtained for different DFT-derived atomic multipoles, we uniformly adopt the CAM-B3LYP functional to parameterize the atomic multipoles and polarizability tensors. ISO-polar and ANISO-polar models are subsequently applied to comparatively predict the molecular crystal structures.

A careful comparison of the molecular packing similarity (RMSD_{20}) between the matched and experimental crystals, as shown in Fig. 4, reveals that for all oligoacene crystals, the ANISO-polar model has the smallest RMSD_{20}, followed by the ISO-polar model, and the NON-polar model has the largest RMSD_{20}. This trend is very similar to the results shown in Fig. 2. Comparing the RMSD20 values from the polarizable models, it is evident that the ANISO-polar model slightly improves the molecular packing similarity compared to the ISO-polar model for naphthalene, anthracene, and tetracene, while it significantly enhances the packing similarity for pentacene.

Furthermore, the overlaps of 21 molecular clusters, separately produced according to the predicted and experimental crystal lattices, are shown in Fig. 5. In Fig. 5, the green (ISO-polar) and yellow (ANISO-polar) molecules overlap well with each other and exhibit slight displacement from the red (experimental) molecules for naphthalene, anthracene, and tetracene. However, for pentacene, the yellow (ANISO-polar) molecules overlap with the red (experimental) molecules, while the red (ISO-polar) molecules are notably displaced from the experimental molecules, particularly for the molecules in two blue circles that are distant from the cluster center.

The notable improvement of the ANISO-polar model in crystal structure prediction, particularly for the longest molecule pentacene, may be attributed to two factors. First, pentacene features enhanced π-electron conjugation and the largest anisotropy of molecular polarizability compared to the other three oligoacenes, indicating that ANISO-polar has significantly improved molecular polarizability tensors for pentacene as shown in Table 3. Secondly, the contribution of electrostatic and polarizable interactions to the total energy is relatively small for neutral oligoacene molecules compared to van der Waals interactions. Among them, the electrostatic and polarizable interactions of pentacene are the most notable compared to the other three oligoacenes, as shown in Table S4 (ESI†). Taken together, these factors suggest that the notable improvement in molecular crystal structure prediction by ANISO-polar is primarily due to significant electrostatic interactions in pentacene and substantial enhancement in polarizable interactions facilitated by ANISO-polar models.

- Y. Shirota and H. Kageyama, Charge Carrier Transporting Molecular Materials and Their Applications in Devices, Chem. Rev., 2007, 107(4), 953–1010 CrossRef CAS PubMed .
- J. Zaumseil and H. Sirringhaus, Electron and Ambipolar Transport in Organic Field-Effect Transistors, Chem. Rev., 2007, 107(4), 1296–1323 CrossRef CAS PubMed .
- C. L. Wang, H. L. Dong, L. Jiang and W. P. Hu, Organic semiconductor crystals, Chem. Soc. Rev., 2018, 47(2), 422–500 RSC .
- J. L. Brédas, D. Beljonne, V. Coropceanu and J. Cornil, Charge-Transfer and Energy-Transfer Processes in e-Conjugated Oligomers and Polymers: A Molecular Picture, Chem. Rev., 2004, 104(11), 4971 CrossRef .
- V. Coropceanu, J. Cornil, D. A. da Silva Filho, Y. Olivier, R. Silbey and J.-L. Brédas, Charge Transport in Organic Semiconductors, Chem. Rev., 2007, 107(4), 926–952 CrossRef CAS PubMed .
- H. Oberhofer, K. Reuter and J. Blumberger, Charge Transport in Molecular Materials: An Assessment of Computational Methods, Chem. Rev., 2017, 117(15), 10319–10357 CrossRef CAS PubMed .
- S. L. Mayo, B. D. Olafson and W. A. Goddard, DREIDING: a generic force field for molecular simulations, J. Phys. Chem., 1990, 94(26), 8897–8909 CrossRef CAS .
- A. D. MacKerell Jr., D. Bashford, M. Bellott, R. L. Dunbrack Jr., J. D. Evanseck, M. J. Field, S. Fischer, J. Gao, H. Guo, S. Ha, D. Joseph-McCarthy, L. Kuchnir, K. Kuczera, F. T. K. Lau, C. Mattos, S. Michnick, T. Ngo, D. T. Nguyen, B. Prodhom, W. E. Reiher, B. Roux, M. Schlenkrich, J. C. Smith, R. Stote, J. Straub, M. Watanabe, J. Wiórkiewicz-Kuczera, D. Yin and M. Karplus, All-Atom Empirical Potential for Molecular Modeling and Dynamics Studies of Proteins, J. Phys. Chem. B, 1998, 102(18), 3586–3616 CrossRef PubMed .
- H. Sun, COMPASS: An ab Initio Force-Field Optimized for Condensed-Phase ApplicationsOverview with Details on Alkane and Benzene Compounds. The, J. Phys. Chem. B, 1998, 102(38), 7338–7364 CrossRef CAS .
- G. M. Day, W. D. S. Motherwell and W. Jones, Beyond the Isotropic Atom Model in Crystal Structure Prediction of Rigid Molecules: Atomic Multipoles versus Point Charges, Cryst. Growth Des., 2005, 5(3), 1023–1033 CrossRef CAS .
- G. M. Day and S. L. Price, A Nonempirical Anisotropic Atom−Atom Model Potential for Chlorobenzene Crystals, J. Am. Chem. Soc., 2003, 125(52), 16434–16443 CrossRef CAS PubMed .
- C. Guo, S. Yin and Y. Wang, Influences of electrostatic models on organic crystal structure prediction – a case study of pentacene, CrystEngComm, 2022, 24(48), 8477–8487 RSC .
- Z.-Q. You, Y. Shao and C.-P. Hsu, Calculating electron transfer couplings by the Spin-Flip approach: energy splitting and dynamical correlation effects, Chem. Phys. Lett., 2004, 390(1–3), 116–123 CrossRef CAS .
- L. Wang, G. Nan, X. Yang, Q. Peng, Q. Li and Z. Shuai, Computational methods for design of organic materials with high charge mobility, Chem. Soc. Rev., 2010, 39(2), 423–434 RSC .
- J. L. Brédas, J. P. Calbert, D. A. da Silva Filho and J. Cornil, Organic semiconductors: A theoretical characterization of the basic parameters governing charge transport, Proc. Natl. Acad. Sci. U. S. A., 2002, 99, 5804 CrossRef PubMed .
- Y. Nagata, Polarizable Atomistic Calculation of Site Energy Disorder in Amorphous Alq3, ChemPhysChem, 2010, 11(2), 474–479 CrossRef CAS PubMed .
- Y. Nagata and C. Lennartz, Atomistic simulation on charge mobility of amorphous tris(8-hydroxyquinoline) aluminum (Alq3). Origin of Poole-Frenkel-type behavior, J. Chem. Phys., 2008, 129(3), 034709 CrossRef PubMed .
- F. May, M. Al-Helwi, B. Baumeier, W. Kowalsky, E. Fuchs, C. Lennartz and D. Andrienko, Design Rules for Charge-Transport Efficient Host Materials for Phosphorescent Organic Light-Emitting Diodes, J. Am. Chem. Soc., 2012, 134(33), 13818–13822 CrossRef CAS PubMed .
- A. Fuchs, T. Steinbrecher, M. S. Mommer, Y. Nagata, M. Elstner and C. Lennartz, Molecular origin of differences in hole and electron mobility in amorphous Alq3-a multiscale simulation study, Phys. Chem. Chem. Phys., 2012, 14, 4259–4270 RSC .
- R. A. Marcus, On the Theory of Oxidation-Reduction Reactions Involving Electron Transfer. I, J. Chem. Phys., 1956, 24(5), 966–978 CrossRef CAS .
- R. A. Marcus, Electron transfer reactions in chemistry. Theory and experiment, Rev. Mod. Phys., 1993, 65(3), 599–610 CrossRef CAS .
- N. S. Hush, Homogeneous and heterogeneous optical and thermal electron transfer, Electrochim. Acta, 1968, 13(5), 1005–1023 CrossRef CAS .
- T. Xu, W. Wang and S. Yin, Explicit Method To Evaluate the External Reorganization Energy of Charge-Transfer Reactions in Oligoacene Crystals Using the State-Specific Polarizable Force Field, J. Phys. Chem. A, 2018, 122(45), 8957–8964 CrossRef CAS PubMed .
- D. P. McMahon and A. Troisi, Evaluation of the External Reorganization Energy of Polyacenes, J. Phys. Chem. Lett., 2010, 1(6), 941–946 CrossRef CAS .
- T. Xu, K. Cao, C. Wang and S. Yin, The effect of asymmetric external reorganization energy on electron and hole transport in organic semiconductors, Phys. Chem. Chem. Phys., 2021, 23(28), 15236–15244 RSC .
- A. K. Rappe and W. A. Goddard III, Charge equilibration for molecular dynamics simulations, J. Phys. Chem., 1991, 95(8), 3358–3363 CrossRef CAS .
- J. A. Lemkul, J. Huang, B. Roux and A. D. MacKerell, An Empirical Polarizable Force Field Based on the Classical Drude Oscillator Model: Development History and Recent Applications, Chem. Rev., 2016, 116(9), 4983–5013 CrossRef CAS PubMed .
- J. W. Ponder, C. Wu, P. Ren, V. S. Pande, J. D. Chodera, M. J. Schnieders, I. Haque, D. L. Mobley, D. S. Lambrecht, R. A. DiStasio, M. Head-Gordon, G. N. I. Clark, M. E. Johnson and T. Head-Gordon, Current Status of the AMOEBA Polarizable Force Field, J. Phys. Chem. B, 2010, 114(8), 2549–2564 CrossRef CAS PubMed .
- J. Applequist, J. R. Carl and K.-K. Fung, Atom dipole interaction model for molecular polarizability. Application to polyatomic molecules and determination of atom polarizabilities, J. Am. Chem. Soc., 1972, 94(9), 2952–2960 CrossRef CAS .
- P. Ren and J. W. Ponder, Consistent treatment of inter- and intramolecular polarization in molecular mechanics calculations, J. Comput. Chem., 2002, 23(16), 1497–1506 CrossRef CAS PubMed .
- B. T. Thole, Molecular polarizabilities calculated with a modified dipole interaction, Chem. Phys., 1981, 59(3), 341–350 CrossRef CAS .
- P. T. van Duijnen and M. Swart, Molecular and Atomic Polarizabilities:Thole's Model Revisited, J. Phys. Chem. A, 1998, 102(14), 2399–2407 CrossRef CAS .
- J. Wang, P. Cieplak, J. Li, T. Hou, R. Luo and Y. Duan, Development of Polarizable Models for Molecular Mechanical Calculations I: Parameterization of Atomic Polarizability, J. Phys. Chem. B, 2011, 115(12), 3091–3099 CrossRef CAS PubMed .
- P. Ren, C. Wu and J. W. Ponder, Polarizable Atomic Multipole-Based Molecular Mechanics for Organic Molecules, J. Chem. Theory Comput., 2011, 7(10), 3143–3161 CrossRef CAS PubMed .
- A. J. Misquitta and A. J. Stone, Distributed polarizabilities obtained using a constrained density-fitting algorithm, J. Chem. Phys., 2006, 124(2), 024111 CrossRef PubMed .
- A. J. Stone, Distributed polarizabilities, Mol. Phys., 1985, 56(5), 1065–1082 CrossRef CAS .
- H. Wang and W. Yang, Determining polarizable force fields with electrostatic potentials from quantum mechanical linear response theory, J. Chem. Phys., 2016, 144(22), 224107 CrossRef PubMed .
- L. H. R. D. Santos, A. Krawczuk and P. Macchi, Distributed Atomic Polarizabilities of Amino Acids and their Hydrogen-Bonded Aggregates, J. Phys. Chem. A, 2015, 119, 3285–3298 CrossRef PubMed .
- E. Heid, P. A. Hunt and C. Schröder, Evaluating excited state atomic polarizabilities of chromophores, Phys. Chem. Chem. Phys., 2018, 20(13), 8554–8563 RSC .
- A. J. Misquitta, A. J. Stone and S. L. Price, Accurate Induction Energies for Small Organic Molecules. 2. Development and Testing of Distributed Polarizability Models against SAPT(DFT) Energies, J. Chem. Theory Comput., 2008, 4(1), 19–32 CrossRef CAS PubMed .
- M. Swart and P. T. van Duijnen, DRF90: a polarizable force field, Mol. Simul., 2006, 32(6), 471–484 CrossRef CAS .
- P. Ren and J. W. Ponder, Polarizable Atomic Multipole Water Model for Molecular Mechanics Simulation. The, J. Phys. Chem. B, 2003, 107, 5933–5947 CrossRef CAS .
- Z.-X. Wang, W. Zhang, C. Wu, H. Lei, P. Cieplak and Y. Duan, Strike a balance: Optimization of backbone torsion parameters of AMBER polarizable force field for simulations of proteins and peptides, J. Comput. Chem., 2006, 27(6), 781–790 CrossRef CAS PubMed .
- S. L. Price, M. Leslie, G. W. A. Welch, M. Habgood, L. S. Price, P. G. Karamertzanis and G. M. Day, Modelling organic crystal structures using distributed multipole and polarizability-based model intermolecular potentials, Phys. Chem. Chem. Phys., 2010, 12(30), 8478–8490 RSC .
- A. K. Das, O. N. Demerdash and T. Head-Gordon, Improvements to the AMOEBA Force Field by Introducing Anisotropic Atomic Polarizability of the Water Molecule, J. Chem. Theory Comput., 2018, 14(12), 6722–6733 CrossRef CAS PubMed .
- A. J. Misquitta and A. J. Stone, Accurate Induction Energies for Small Organic Molecules: 1. Theory, J. Chem. Theory Comput., 2008, 4(1), 7–18 CrossRef CAS PubMed .
- A. Krawczuk, D. Pérez and P. Macchi, PolaBer: a program to calculate and visualize distributed atomic polarizabilities based on electron density partitioning, J. Appl. Crystallogr., 2014, 47, 1452–1458 CrossRef CAS .
- T. Xu, W. Wang and S. Yin, Electrostatic Polarization Energies of Charge Carriers in Organic Molecular Crystals: A Comparative Study with Explicit State-Specific Atomic Polarizability Based AMOEBA Force Field and Implicit Solvent Method, J. Chem. Theory Comput., 2018, 14(7), 3728–3739 CrossRef CAS PubMed .
- J. W. Ponder, Tinker, version 8.1, https://dasher.wustl.edu/tinker/, 2011 Search PubMed .
- A. J. Stone, The Theory of Intermolecular Forces. Oxford University press: Great Clarendon Street, 2013 Search PubMed .
- T. A. Halgren, The representation of van der Waals (vdW) interactions in molecular mechanics force fields: potential form, combination rules, and vdW parameters, J. Am. Chem. Soc., 1992, 114(20), 7827–7843 CrossRef CAS .
- A. J. Stone, Distributed Multipole Analysis of Gaussian wavefunctions, GDMA version 2.2.06, https://www-stone.ch.cam.ac.uk/pub/gdma/, 2012 Search PubMed .
- Z. Konkoli and D. Cremer, A new way of analyzing vibrational spectra. III. Characterization of normal vibrational modes in terms of internal vibrational modes, Int. J. Quantum Chem., 1998, 67(1), 29–40 CrossRef CAS .
- W. Zou, R. Kalescky, E. Kraka and D. Cremer, Relating normal vibrational modes to local vibrational modes with the help of an adiabatic connection scheme, J. Chem. Phys., 2012, 137(8), 084114 CrossRef PubMed .
- C. P. Brock and J. D. Dunitz, Temperature-Dependence of Thermal Motion in Crystalline Naphthalene, Acta Crystallogr., Sect. B: Struct. Crystallogr. Cryst. Chem., 1982, B38, 2218–2228 CrossRef CAS .
- C. P. Brock and J. D. Dunitz, Temperature-Dependence of Thermal Motion in Crystalline Anthracene, Acta Crystallogr., Sect. B: Struct. Sci., 1990, B46, 795–806 CrossRef CAS .
- R. B. Campbell, J. M. Robertson and J. Trotter, The crystal structure of hexacene, and a revision of the crystallographic data for tetracene, Acta Crystallogr., 1962, 15(3), 289–290 CrossRef CAS .
- R. B. Campbell, J. M. Robertson and J. Trotter, The crystal and molecular structure of pentacene, Acta Crystallogr., 1961, 14(7), 705–711 CrossRef CAS .
- A. A. Aina, A. J. Misquitta and S. L. Price, A non-empirical intermolecular force-field for trinitrobenzene and its application in crystal structure prediction, J. Chem. Phys., 2021, 154(9), 094123 CrossRef CAS PubMed .
- D. J. Willock, S. L. Price, M. Leslie and C. R. A. Catlow, The relaxation of molecular crystal structures using a distributed multipole electrostatic model, J. Comput. Chem., 1995, 16(5), 628–647 CrossRef CAS .
- D. C. Liu and J. Nocedal, On the limited memory BFGS method for large scale optimization, Math. Prog., 1989, 45(1), 503–528 CrossRef .
- J. V. D. Streek and M. A. Neumann, Validation of experimental molecular crystal structures with dispersion-corrected density functional theory calculations, Acta Crystallogr., Sect. B: Struct. Sci., 2010, 66(5), 544–558 CrossRef PubMed .
- J. Rohlicek and E. Skorepova, CrystalCMP: automatic comparison of molecular structures, J. Appl. Crystallogr., 2020, 53(3), 841–847 CrossRef CAS PubMed .
- J. Rohlíček, E. Skořepová, M. Babor and J. Čejka, CrystalCMP: An easy-to-use tool for fast comparison of molecular packing, J. Appl. Crystallogr., 2016, 49, 2172–2183 CrossRef .
- C. R. Groom, I. J. Bruno, M. P. Lightfoot and S. C. Ward, The Cambridge Structural Database, Acta Crystallogr., Sect. B: Struct. Sci., Cryst. Eng. Mater., 2016, 72(2), 171–179 CrossRef CAS PubMed .
- A. O. Lyakhov, A. R. Oganov, H. T. Stokes and Q. Zhu, New developments in evolutionary structure prediction algorithm USPEX, Comput. Phys. Commun., 2013, 184(4), 1172–1182 CrossRef CAS .
- A. R. Oganov and C. W. Glass, Crystal structure prediction using ab initio evolutionary techniques: Principles and applications, J. Chem. Phys., 2006, 124(24), 244704 CrossRef PubMed .
- M. Pakhnova, I. Kruglov, A. Yanilkin and A. Oganov, The search for stable cocrystals of energetic materials using evolutionary algorithm USPEX, Phys. Chem. Chem. Phys., 2020, 22(29), 16822–16830 RSC .
- R. F. W. Bader, A. Larouche, C. Gatti, M. T. Carroll, P. J. MacDougall and K. B. Wiberg, Properties of atoms in molecules: Dipole moments and transferability of properties, J. Chem. Phys., 1987, 87(2), 1142–1152 CrossRef CAS .
- C. R. L. Sueur and A. J. Stone, Practical schemes for distributed polarizabilities, Mol. Phys., 1993, 78(5), 1267–1291 CrossRef .
- T. A. Keith, AIMAll Version 16.10.31, TK Gristmill Software, Overland Park KS, USA: 2016 Search PubMed .
- M. J. Frisch; G. W. Trucks; H. B. Schlegel; G. E. Scuseria; M. A. Robb; J. R. Cheeseman; G. Scalmani; V. Barone; B. Mennucci; G. A. Petersson; H. Nakatsuji; M. Caricato; X. Li; H. P. Hratchian; A. F. Izmaylov; J. Bloino; G. Zheng; J. L. Sonnenberg; M. Hada; M. Ehara; K. Toyota; R. Fukuda; J. Hasegawa; M. Ishida; T. Nakajima; Y. Honda; O. Kitao; H. Nakai; T. Vreven; J. A. Montgomery; J. E. Peralta; F. Ogliaro; M. Bearpark; J. J. Heyd; E. Brothers; K. N. Kudin; V. N. Staroverov; R. Kobayashi; J. Normand; K. Raghavachari; A. Rendell; J. C. Burant; S. S. Iyengar; J. Tomasi; M. Cossi; N. Rega; J. M. Millam; M. Klene; J. E. Knox; J. B. Cross; V. Bakken; C. Adamo; J. Jaramillo; R. Gomperts; R. E. Stratmann; O. Yazyev; A. J. Austin; R. Cammi; C. Pomelli; J. W. Ochterski; R. L. Martin; K. Morokuma; V. G. Zakrzewski; G. A. Voth; P. Salvador; J. J. Dannenberg; S. Dapprich; A. D. Daniels; O. Farkas; J. B. Foresman; J. V. Ortiz; J. Cioslowski and D. J. Fox, Gaussian 09,Revision C.01, Gaussian, Inc., Wallingford CT, 2009 Search PubMed .
- Z. G. Soos and G. W. Hayden, Static polarizability of interacting π electrons in conjugated polymers, Phys. Rev. B: Condens. Matter Mater. Phys., 1989, 40(5), 3081–3089 CrossRef PubMed .

## Footnote |

† Electronic supplementary information (ESI) available: The lattice parameters of refined experimental crystal structures, utilizing NON-polar, ISO-polar, and ANISO-polar models, are presented in Table S1. Fig. S1 illustrates the crystal structure similarity (t-RMSD) between the refined and experimental structures of naphthalene and anthracene. Table S2 displays the lattice parameters of all matched crystals derived from USPEX-predicted crystals, employing ANISO-polar and ISO-polar models. Table S3 provides the distributed atomic polarizability tensors. Table S4 showcases the results of energy decomposition analysis for the optimized experimental crystals, employing NON-polar, ISO-polar, and ANISO-polar models. Table S5 lists molecular polarizabilities calculated from different functionals. See DOI: https://doi.org/10.1039/d4cp01568a |

This journal is © the Owner Societies 2024 |