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

Development of an anisotropic polarizable model for the all-atom AMOEBA force field†

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

Abstract

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.

1. Introduction

π-Conjugated molecules and polymers have garnered significant interest as active functional materials in both crystalline and thin-film electronic and electro-optic applications.1–3 In various electro-optic applications, the rates of charge-carrier transport crucially determine device performance, primarily influenced by key parameters: (1) the strength of electronic coupling between adjacent molecules, (2) the reorganization energy arising from charge transfer processes between sites, and (3) the degree of disorder in site energies within the solid phase.4–6 For weak interaction-based organic molecules, all these parameters are intricately linked to intermolecular interactions, including electrostatic and inductive interactions.

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 scheme29 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 AMOBEA42 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 method40,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.

2. Theory and computational details

2.1. Induced dipole model

In induced-dipole-based polarizable force field, the induced dipole is solved by the following equation
 
μindi = αi(Fpermi + Findi) (1)
where F represents the applied electric field, α denotes the atomic polarizability, the superscript perm and ind refer to the i-th atomic feeling electric fields respectively coming from the contributions of permanent (perm) atomic multipoles and induced (ind) dipole contributions and the subscript i represents the atomic index. In the AMOEBA polarizable force field,28 atomic polarizability (α) adopts the isotropic polarizability. Thus, we call this polarizable model as an isotropic polarizable model. When α adopts the second-order tensor form, then eqn (1) becomes an anisotropic polarizable model. In the standard AMOEBA force field, the direct electric field (Fpermi) and mutual electric field (Findi) are defined as the gradient of electrostatic potentials including the contributions of permanent multipoles and induced dipoles. Thus,
 
image file: d4cp01568a-t1.tif(2a)
 
image file: d4cp01568a-t2.tif(2b)
where ∇ means the gradient operator and Vi(Mj) and image file: d4cp01568a-t3.tif represent the electrostatic potentials from the contributions of atomic permanent multipoles (Mj) and mutual induced dipole image file: d4cp01568a-t4.tif, respectively. For the detailed expression of the electrostatic potential, please see Stone's textbook notations.50

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

 
image file: d4cp01568a-t5.tif(3)
where α−1i means the inverse of atomic polarizability.

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 Thole31 or Applequist29 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,
 
image file: d4cp01568a-t6.tif(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.2. Electrostatic and van del Waals interactions

For clarity, the electrostatic and van der Waals (vdW) interactions strictly adhere to the AMOEBA functional form, specifically the use of atomic multipoles (AMPs) and the buffered 14-7 vdW formula.51 The comprehensive vdW functional form and parameters for aromatic C and H are directly obtained from the literature by Ponder and Ren.28,34 AMPs are derived from quantum mechanics (QM) calculated electron densities using the GDMA program developed by Prof. Stone.52 To account for AMP-dependent electrostatics, various density functional theory (DFT) functionals, including GGA PBE, meta-GGA M06 and M06-L, hybrid functionals B3LYP, range-separation hybrid functionals CAM-B3LYP and ωB97X-D, and double hybrid B2PLYP functional, are employed to calculate the electron densities of naphthalene and anthracene.

2.3. Intramolecular force field

Intramolecular bond stretching, angle bending, and torsional rotation follow the CHARMM functional form,8 which is detailed in our previous works.23,25 Utilizing the adiabatic connection scheme, the local vibrations of bond stretching and angle bending are correlated with QM-based normal vibrations, and the force constants of these local vibrations are acquired through the LocalModes program.53,54 Balanced bond lengths and angles are directly determined from QM optimized molecular structures. Torsional force constants adhere to the default parameters provided by AMOEBA, obtained from the valence program within the Tinker package.34

2.4. Refinement and prediction of crystal structure

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 method61 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 as62
 
image file: d4cp01568a-t7.tif(5)
where ri 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 (RMSDn) 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 RMSD20 is defined as
 
image file: d4cp01568a-t8.tif(6)
where K the number of atoms in molecular clusters, di 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 RMSD20 differs slightly from the RMSD20 definition in Mercury,65 where the comparison involves clusters of 20 molecules.
Table 1 Interaction potentials for the three force fields
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. Computational details

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
 
image file: d4cp01568a-t9.tif(7)
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. Fy 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 R3 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.

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 (αglobali) 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.


image file: d4cp01568a-f1.tif
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,

 
image file: d4cp01568a-t10.tif(8)
where ma is the atomic mass and a represents the atoms of this molecule.

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

 
image file: d4cp01568a-t11.tif(9)
where C is the 3 × 3 modal matrix, C−1 is the inverse matrix of C, I1/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.

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

 
αglobali = CαstiC−1 (10)
where αsti 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.

(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

3. Results and discussions

3.1. QM-derived atomic polarizability

The QM-derived atomic polarizabilities of naphthalene are listed in Table 2. And the atomic polarizabilities for other oligoacenes including anthracene, tetracene and pentacene are listed in Table S3 (ESI†). We define the numerical isotropic atomic polarizability as αnumiso = (αxx + αyy + αzz)/3. As naphthalene belongs to the D2h point group, there are five inequivalent atoms i.e. 1-C, 2-C, 3-C, 6-H and 7-H as shown in Fig. 1. Their numerical isotropic polarizabilities (αnumiso) are 1.333, 1.422, 1.504, 0.384 and 0.343 in 4πε0 Å3 unit, respectively. Under the molecular standard coordinates shown in Fig. 1,the diagonal elements (αxx, αyy and αzz) for the equivalent C atoms such as 1-C, 5-C, 10-C and 14-C have identical values, while the off-diagonal elements of equivalent atoms have either the same or opposite values depending on the symmetry. The sum of diagonal atomic polarizabilities is equal to the corresponding diagonal elements of the molecular polarizability as shown in Table 3, while the sum of off-diagonal elements of atomic polarizabilities equals zero. This indicates that the molecular polarizability tensor is accurately decomposed by our numerical atomic polarizabilities.
Table 2 QM derived atomic polarizability of naphthalene using eqn (7). The unit of polarizability is in 4πε0 Å3. The x, y and z directions correspond to the long, short and normal axes of naphthalene. Atomic type is the same as the atomic index as shown in Fig. 1
Type αxx αxy αxz αyz αyy αzz αnumiso
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

Table 3 Molecular polarizability tensors and their anisotropy (k) of oligoacene obtained from QM calculations at the CAM-B3LYP/6-311+G(d,p) level, isotropic atomic polarizability (ISO-polar) and anisotropic atomic polarizability (ANISO-polar). The unit of polarizability (α) 4πε0 Å3. The x, y and z are defined under the molecular standard coordinates i.e. the molecular x (short), y (long) and z (normal) axes, respectively. The numbers in parentheses are relative errors of the anisotropy of molecular polarizability obtained from ISO-polar and QM methods
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

 
image file: d4cp01568a-t12.tif(11)
where αxx, αyy and αzz respectively are the three principal molecular polarizability elements and αiso is the isotropic molecular polarizability. From eqn (11), it is evident that k reaches its maximum value, i.e., 1 for linear molecules, and k attains its minimum value, i.e., 0 for isotropic single atoms. For most molecules, k falls within the range from 0 to 1. Analyzing the k values of oligoacenes obtained from QM calculations in Table 3, we observe that the anisotropy of molecular polarizability follows the order: naphthalene < anthracene < tetracene < pentacene, with the longest molecule, pentacene, exhibiting the most significant anisotropy. The ANISO-polar model accurately reproduces both the QM-calculated molecular polarizability and its corresponding anisotropy parameters k. However, ISO-polar significantly underestimates the anisotropy of oligoacenes, especially for the longest molecule, pentacene, with the largest relative error of approximately 32.6%. Based on the molecular polarizability tensor, it can be concluded that the ANISO-polar model provides more accurate electron polarization behaviors than ISO-polar, particularly for molecules with the extensive π-electron conjugation, such as pentacene.

3.2. Refinement of experimental crystal structures

After parameterizing the three types of force fields, namely NON-polar, ISO-polar, and ANISO-polar models, using the method outlined above, we separately applied these models to optimize the experimental crystal structures of naphthalene (NAPHTA06) and anthracene (ANTCEN01). The atomic multipole moments (AMPs) were derived using seven DFT functionals to calculate electron wavefunctions. For the anisotropic polarizable model, we uniformly adopted the numerical atomic polarizabilities coupled with respective density functionals at the 6-311+G(d,p) level. The molecular packing similarity (RMSD20) between the refined and experimental structures is illustrated in Fig. 2. Detailed lattice parameters and similarity information of the refined crystal structures are provided in Table S1 (ESI†). Overall, regardless of the choice of functionals, NON-polar models (orange bars) exhibited the largest RMSD20, followed by ISO-polar models (green bars), with ANISO-polar models (purple bars) showing the smallest RMSD20. This suggests that the ANISO-polar model slightly outperforms the others in optimizing molecular crystals.
image file: d4cp01568a-f2.tif
Fig. 2 Molecular packing similarity (RMSD20) 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 RMSD20 values for naphthalene (0.49 Å) and anthracene (0.49 Å) are derived from the double hybrid functional B2PLYP, while the largest RMSD20 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.


image file: d4cp01568a-f3.tif
Fig. 3 The electrostatic potentials (ESP) in kcal mol−1 on the molecular 0.001 iso-electron density surface obtained from seven DFT functional calculations. AMP/CAM-B3LYP refers to the ESP calculated using atomic multipoles (AMP) derived from CAM-B3LYP density. QM-AMP refers to the difference in ESP between QM/CAM-B3LYP and AMP/CAM-B3LYP.

ISO-polar models (green bars) demonstrate smaller RMSD20 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 (RMSD20). 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.

3.3. Prediction of crystal structures

Based on CAM-B3LYP calculations, the force field parameters for the anisotropic polarizable models of oligoacenes, ranging from naphthalene to pentacene, are determined using the previously mentioned methods. The molecular crystal structures are predicted using the USPEX program coupled with the xtalmin crystal optimizer. The crystal candidates’ total energies are utilized to rank the predicted crystals, where the most stable crystal with the lowest energy is designated as Rank 1, the subsequent crystal with the second lowest energy as Rank 2, and so forth. The low-energy regions of USPEX-generated crystals with three models are depicted in Fig. 4. All the predicted crystal structures, whether using the NON-polar, ISO-polar, or ANISO-polar models, the matched structures of naphthalene (NAPHTA06), anthracene (ANTCEN01), tetracene (TETCEN), and pentacene (PENCEN) are positioned at very low total energy levels in their respective energy landscapes. This means that for oligoacene systems, AMP-based electrostatics have already captured the key factors needed to predict crystal structures. The performance of CSP using the NON-polar model is very similar to the results of the ISO-polar and ANISO-polar models, possibly due to the weak electrostatic interactions in non-polar oligoacenes, leading to weaker inductive interactions. Besides the slight improvement in molecular packing similarity, as seen from RMSD20 with black numbers shown in Fig. 4, the polarizable models stabilize the total energy of the matched crystals more effectively by considering the polarization effect, as seen from the total energy with red numbers shown in Fig. 4.
image file: d4cp01568a-f4.tif
Fig. 4 The energy landscapes of crystal structure predictions of (a) naphthalene (the first row), (b) anthracene (the second row), (c) tetracene (the third row) and (d) pentacene (the forth row) using NON-polar (left column), ISO-polar (middle column) and ANISO-polar (right column) models, respectively. The matched structures are labeled by an arrow. The number in black is the molecular packing similarity (RMSD20 unit in Å) between matched structures and experimental structures and the number in red is the total energy of the matched structure in kcal mol−1.

A careful comparison of the molecular packing similarity (RMSD20) between the matched and experimental crystals, as shown in Fig. 4, reveals that for all oligoacene crystals, the ANISO-polar model has the smallest RMSD20, followed by the ISO-polar model, and the NON-polar model has the largest RMSD20. 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.


image file: d4cp01568a-f5.tif
Fig. 5 The overlaps of the predicted crystal structures using ISO-polar (green) and ANISO-polar(yellow) models with the experimental (red) structure. RMSD20 values between the predicted and experimental structures are shown in parentheses.

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.

4. Conclusions

Linear response theory coupled with atoms in molecules (AIMALL) theory can effectively obtain the distributed atomic polarizability tensor, which can exactly reproduce the QM-based molecular polarizability. In this way, the anisotropic atomic polarizability can be derived from QM-based calculations. Utilizing the QM-derived atomic polarizability tensors, an anisotropic polarizable model is developed to describe polarization interactions. This model is integrated into the Tinker program, extending the capabilities of the AMOEBA polarizable model. By leveraging atomic multipoles derived from seven DFT functionals, non-polar, anisotropic, and isotropic polarizable models are individually employed to refine the experimental crystal structures of naphthalene and anthracene. The results indicate that non-polar atomic multipoles from different functionals exhibit similar performance in optimizing crystal structures. Polarizable models enhance crystal refinement compared to non-polar models, with the anisotropic polarizable model demonstrating superior performance in optimizing crystals. Furthermore, both anisotropic and isotropic models are applied to predict the molecular crystal structures of oligoacenes using the USPEX program. The lowest energy USPEX-generated crystals with both polarizable models exhibit herringbone packing structures, successfully predicting the molecular crystal structures. For naphthalene, anthracene, and tetracene, the anisotropic model slightly improves the packing similarity between predicted and experimental structures. Notably, for pentacene, which has the largest anisotropy of molecular polarizability, the anisotropic polarizable model significantly enhances the packing similarity due to its more accurate polarizable interactions.

Data availability

The data supporting this article have been included as part of the ESI.†

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

SW Yin acknowledges the National Science Foundation of China (Grant No. 22273054 and 22073060), and the Ministry of Education, China 111 Project (Grant No. B14041) for financial support.

References

  1. 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 .
  2. J. Zaumseil and H. Sirringhaus, Electron and Ambipolar Transport in Organic Field-Effect Transistors, Chem. Rev., 2007, 107(4), 1296–1323 CrossRef CAS PubMed .
  3. C. L. Wang, H. L. Dong, L. Jiang and W. P. Hu, Organic semiconductor crystals, Chem. Soc. Rev., 2018, 47(2), 422–500 RSC .
  4. 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 .
  5. 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 .
  6. 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 .
  7. 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 .
  8. 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 .
  9. H. Sun, COMPASS:[thin space (1/6-em)] 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 .
  10. G. M. Day, W. D. S. Motherwell and W. Jones, Beyond the Isotropic Atom Model in Crystal Structure Prediction of Rigid Molecules:[thin space (1/6-em)] Atomic Multipoles versus Point Charges, Cryst. Growth Des., 2005, 5(3), 1023–1033 CrossRef CAS .
  11. 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 .
  12. 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 .
  13. 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 .
  14. 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 .
  15. 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 .
  16. Y. Nagata, Polarizable Atomistic Calculation of Site Energy Disorder in Amorphous Alq3, ChemPhysChem, 2010, 11(2), 474–479 CrossRef CAS PubMed .
  17. 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 .
  18. 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 .
  19. 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 .
  20. R. A. Marcus, On the Theory of Oxidation-Reduction Reactions Involving Electron Transfer. I, J. Chem. Phys., 1956, 24(5), 966–978 CrossRef CAS .
  21. R. A. Marcus, Electron transfer reactions in chemistry. Theory and experiment, Rev. Mod. Phys., 1993, 65(3), 599–610 CrossRef CAS .
  22. N. S. Hush, Homogeneous and heterogeneous optical and thermal electron transfer, Electrochim. Acta, 1968, 13(5), 1005–1023 CrossRef CAS .
  23. 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 .
  24. 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 .
  25. 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 .
  26. A. K. Rappe and W. A. Goddard III, Charge equilibration for molecular dynamics simulations, J. Phys. Chem., 1991, 95(8), 3358–3363 CrossRef CAS .
  27. 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 .
  28. 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 .
  29. 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 .
  30. 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 .
  31. B. T. Thole, Molecular polarizabilities calculated with a modified dipole interaction, Chem. Phys., 1981, 59(3), 341–350 CrossRef CAS .
  32. 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 .
  33. 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 .
  34. 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 .
  35. 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 .
  36. A. J. Stone, Distributed polarizabilities, Mol. Phys., 1985, 56(5), 1065–1082 CrossRef CAS .
  37. 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 .
  38. 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 .
  39. 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 .
  40. 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 .
  41. M. Swart and P. T. van Duijnen, DRF90: a polarizable force field, Mol. Simul., 2006, 32(6), 471–484 CrossRef CAS .
  42. 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 .
  43. 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 .
  44. 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 .
  45. 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 .
  46. A. J. Misquitta and A. J. Stone, Accurate Induction Energies for Small Organic Molecules:[thin space (1/6-em)] 1. Theory, J. Chem. Theory Comput., 2008, 4(1), 7–18 CrossRef CAS PubMed .
  47. 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 .
  48. 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 .
  49. J. W. Ponder, Tinker, version 8.1, https://dasher.wustl.edu/tinker/, 2011 Search PubMed .
  50. A. J. Stone, The Theory of Intermolecular Forces. Oxford University press: Great Clarendon Street, 2013 Search PubMed .
  51. 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 .
  52. 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 .
  53. 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 .
  54. 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 .
  55. 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 .
  56. 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 .
  57. 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 .
  58. 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 .
  59. 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 .
  60. 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 .
  61. D. C. Liu and J. Nocedal, On the limited memory BFGS method for large scale optimization, Math. Prog., 1989, 45(1), 503–528 CrossRef .
  62. 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 .
  63. J. Rohlicek and E. Skorepova, CrystalCMP: automatic comparison of molecular structures, J. Appl. Crystallogr., 2020, 53(3), 841–847 CrossRef CAS PubMed .
  64. 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 .
  65. 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 .
  66. 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 .
  67. 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 .
  68. 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 .
  69. 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 .
  70. C. R. L. Sueur and A. J. Stone, Practical schemes for distributed polarizabilities, Mol. Phys., 1993, 78(5), 1267–1291 CrossRef .
  71. T. A. Keith, AIMAll Version 16.10.31, TK Gristmill Software, Overland Park KS, USA: 2016 Search PubMed .
  72. 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 .
  73. 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