To interpret ultrafast dynamics experiments on large molecules, computer simulation is required due to the complex response to the laser field. We present a method capable of efficiently computing the static electronic response of large systems to external electric fields. This is achieved by extending the density-functional tight binding method to include larger basis sets and by multipole expansion of the charge density into electrostatically interacting Gaussian distributions. Polarizabilities for a range of hydrocarbon molecules are computed for a multipole expansion up to quadrupole order, giving excellent agreement with experimental values, with average errors similar to those from density functional theory, but at a small fraction of the cost. We apply the model in conjunction with the polarizable-point-dipoles model to estimate the internal fields in amorphous poly(3-hexylthiophene-2,5-diyl).

## I. INTRODUCTION

Recent advances in the field of ultrafast dynamics have made it possible to characterize the harmonic response of small molecules with sub-femtosecond resolution.^{1} While experimental work is rapidly progressing, it is nonetheless very difficult to obtain information about the dynamics of the excited electrons directly. Computational modelling holds the potential to assist the interpretation of experimental data, although the systems encountered in ultrafast experiments are particularly complex due to the presence of strong, time-dependent electric fields with intensities on the order of 10^{14} W/cm^{2}.

Simulations of such systems require computational methods with a reliable description of electrostatic polarization, as well as being sufficiently computationally efficient to enable demanding time-dependent simulations of molecule in condensed phases. The challenge is hence to develop a theory with predictive capabilities comparable to density-functional theory (DFT), but with significantly less computational effort. We shall in this work consider the electrostatic response of molecules to external fields, with an investigation on the time-dependent response following in a future publication.

Density-functional tight binding (DFTB)^{2} is a particularly promising method, as it offers an ideal trade-off between computational time and accuracy. Current formulations provide a poor description of electronic polarizability, necessitating further improvements to the model. The implementation of a first principles, self-consistent charge multipoles formalism in DFTB has been suggested before.^{3} Here we formulate an extended DFTB approach including self-consistent polarized charges and polarization orbitals. The electrostatic integrals are evaluated analytically for the chosen basis set by using Gaussian expansions, resulting in an internally consistent theory without the need for a fitting procedure. The resulting Gaussian Tight Binding (GTB) model gives molecular polarizabilities in excellent agreement with experimental data for hydrocarbons, having errors of the same order as DFT using the PBE^{4} exchange-correlation functional. Various empirical extensions to tight-binding with a self-consistent description of ion polarizability have previously been formulated and tested, such as the polarizable-ion tight-binding model,^{5,6} or parametrically extended DFTB^{7–9} models. We present a systematically improvable model derived from first-principles within the DFTB formalism.

In Sec. II we outline the GTB model. The accuracy is assessed by the level of agreement with PBE-DFT and experimental data on the example of various hydrocarbons’ mean polarizabilities. Finally we use GTB in conjunction with the polarizable-point-dipoles (PPD) method to estimate the internal electric fields of amorphous bulk poly(3-hexylthiophene-2,5-diyl) under the influence of an external field.

## II. GAUSSIAN TIGHT BINDING

### A. Density functional tight binding

See Refs. 2 and 10–12 for extensive reviews on the density-functional tight-binding method. We start with the total energy expression of the one-electron mean-field system,^{10}

Here *ρ* is the electron number density, *ψ*^{(n)} the eigenstate of level *n* with occupancy *f _{n}*, $ H \u02c6 [ \rho ] $ the Kohn-Sham Hamiltonian,

*E*the ion-ion repulsion, and

_{ion}*E*the double-counting term,

_{dc}*E _{Ha}* is the Hartree energy and

*E*the exchange-correlation functional. We now define a reference density

_{xc}*ρ*

^{0}as the superposition of spherical, neutral atomic densities $ \rho I 0 $, where the subscript

*I*denotes the atomic site. The total density is then given by the sum of the reference density and the density variation,

*δρ*:

*ρ*=

*ρ*

^{0}+

*δρ*. This decomposition is used to expand the total energy in terms of the density variation,

*δρ*, in atomic units,

Here *f _{xc}* is the functional derivative of the exchange-correlation potential. Up to this point the expression is still exact to

*O*(

*δρ*

^{2}). Truncating the series at first order in

*δρ*returns the non-self-consistent Harris-Foulkes energy, with higher orders reintroducing self-consistency. In practice we truncate the series at the second order. We employ the common tight-binding approximations, which are the neglect of three-centre and crystal field terms; see Refs. 11 and 12 for an extensive review. The Hamiltonian matrix elements of the reference density H[

*ρ*

^{0}] are tabulated for a chosen atomic basis set, while the remaining terms dependent on the reference density are approximated by a repulsive pair-potential.

^{13}

### B. Multipole expansion

Our attention is now directed to the second order correction to the Hartree energy,

as it is particularly expensive to evaluate in a linear-combination of atomic orbitals approach due to the involvement of four-centre integrals. There are several established ways to reduce the computational load, for instance, by resolution of the identity approximation or density fitting.^{14} We want to capture the essential physics of electrostatic interaction, while keeping the computational time to a minimum. For this purpose, we have chosen to adopt a simple density fitting scheme based on a multipole expansion. We employ an expansion in real cubic harmonics *K _{lm}* to avoid using a complex basis set,

*l*and degree

*m*, with −

*l*≤

*m*≤

*l*. We have chosen the principal axis to lie in the z-direction, hence

**= (**

*r**K*

_{10}

*K*

_{11}

*K*

_{12})

^{T}= (

*z x y*)

^{T}.

The density variation *δρ* can be written as a superposition of localized atomically centered components $\delta \rho = \u2211 I \delta \rho I $ without error, provided a finite range atomic basis set is used. Without introducing any further approximations, we express these in radial and angular components by a multipole expansion,

We wish to express *a _{Ilm}*(

*r*) by a set of radial functions

*f*(

_{Ilm}*r*) for which there is a straightforward integrable solution of the electrostatic integral (Section II C). We approach this problem by defining the charge density moments

*Q*as the expansion coefficients of the radial functions

_{Ilm}*f*(

_{Ilm}*r*),

where the charge density moments *Q _{Ilm}* are defined in terms of the cubic harmonics,

We find a normalisation condition for *f _{Ilm}* by substituting the variational density (7) into the density moments (8),

The variational density (7) is now substituted into the second order correction of the Hartree energy (4), with *R*_{IJ} = *R*_{J} − *R*_{I},

where the integral can be identified with the Madelung structure constant *B*_{Ilm,Jl′m′}(*R*_{IJ}), leading to the compact expression,

This formalism enables us to express *δE _{Ha}* as a function of the charge density moments

*Q*(8). The total energy can hence be self-consistently minimized with respect to the moments

_{Ilm}*Q*.

_{Ilm}### C. Self-consistent solution

The second order correction to the Hartree energy *δE _{Ha}* depends directly on the electron density

*ρ*via the density moments

*Q*(8). Solving for

_{Ilm}*ρ*will hence require a self-consistent approach. We shall first define the electron density in a linear-combination of atomic orbitals approach.

Let the atom on site *I* have real atomic-type orbitals *ϕ*_{Iα} attached to it. The eigenstate *ψ*^{(n)} then expands into $ \psi ( n ) ( r ) = \u2211 I \alpha c I \alpha ( n ) \varphi I \alpha ( r \u2212 R I ) $, where $ c I \alpha ( n ) $ are the expansion coefficients. Assuming real orbitals, the total electron density is given by

with the density matrix

where *f _{n}* is the electron occupancy of state

*n*. Substituting this definition of the density (12) into (8), we can express the charge moments

*Q*in terms of the density matrix

_{Ilm} where we define *Z _{I}* as the ionic charge of site

*I*and $ M I \alpha J \beta l m $ as the cubic moment integral over the basis functions

*ϕ*

_{Iα}and

*ϕ*

_{Jβ},

Here we expressed the density variation as $\delta \rho I ( r ) = \rho I ( r ) \u2212 \rho I 0 ( r ) $. The atomic reference density $ \rho I 0 ( r ) $ is spherically symmetric, therefore $\u222b \rho I 0 ( r \u2212 R I ) K l m ( r \u2212 R I ) = Z I \delta l 0 \delta m 0 $.

The self-consistent Hamiltonian contribution *δH* is subsequently found by differentiating the Hartree energy correction with respect to the expansion coefficients

giving

We control the level of accuracy in *δE _{Ha}* by varying the order of the multipole expansion. Consider a zero-order expansion, where

*δρ*is approximated as spherically symmetric. The self-consistent solution will add or remove charge at each site to minimize the total energy. This approximation is known as self-charge-consistent tight-binding (SCC-TB).

_{I}^{15}

While the zeroth order approximation enables a basic description of charge transfer, it is still a rather crude picture. Our model allows a flexible choice of the expansion order, allowing us to treat *δE _{Ha}* in a systematically improvable framework. A first-order multipole expansion, for example, considers

*δρ*to be polarizable, leading to an overall improved picture of electrostatic interaction.

_{I}### D. Implementation

We have not yet chosen a form for our radial functions *f _{Ilm}*(

*r*). There are two points to consider for the choice. Firstly, the radial functions will be used to evaluate the integral (10), hence it is beneficial to choose functions for which the integral becomes analytical. Secondly, the quality of approximation improves the closer

*f*(

_{Ilm}*r*) matches the true radial distribution

*a*(

_{Ilm}*r*). Past work in the monopole approximation has shown a single Gaussian function for

*f*(

_{Ilm}*r*) to perform well;

^{16}this motivates us to pursue this approach here,

wherer *k _{Ilm}* is the normalisation constant introduced to satisfy Equation (9),

where *α _{Ilm}* is the inverse square width of the Gaussian charge distribution associated with site

*I*and angular momentum {

*l*,

*m*}. In practice we use the same Gaussian width

*α*over all momenta for each atomic species. This width can be associated with the Hubbard parameter

_{I}*U*of the atomic species as shown by Elstner

_{I}*et al.*

^{2}We follow this approach here, using

While this value can be calculated from first principles, we leave it as an adjustable parameter and use *U _{I}* from the mio-1-1 parameter set,

^{2}averaged over all angular momenta. The GTB model is implemented in the tight-binding code plato

^{17}and will be released in the future.

The computational overhead of the multipole model is dominated by the use of the larger polarizable basis set, slowing down the self-consistent solution by direct diagonalization roughly by a factor of 10. Some DFTB parametrizations of certain atoms already include polarization orbitals by default,^{18,19} and would therefore not be affected by this computational cost.

Our current implementation constructs the multipole moments on-the-fly explicitly from the Gaussian basis set. While this is not a significant expense for single-point calculations, it represents about 3/4 of the computational effort in time-dependent runs with forces (plato implementation). However, we emphasize that the multipole integrals can also be obtained from tabulated Slater-Koster overlap tables of even higher shells. This would almost entirely eliminate the overhead associated with constructing the integrals, and only leave the overhead associated with using the larger basis-set size.

### E. Coupling to external fields

We can easily introduce coupling to an external uniform electric field ** E**, as the dipole moments are already evaluated for the self-consistent minimization. Our total Hamiltonian will then be composed of the Harris-Foulkes tight binding Hamiltonian

*H*

_{0}, the self-consistent contribution

*δH*, and the external field contribution

*δH*,

^{ext} The external field coupling energy is found by integration over the total charge density *n*(** r**),

The monopole moment *Q*_{I00} can be identified as the net charge of the total atomic density *ρ* on site *I*. The contribution *δH ^{ext}* to the self-consistent Hamiltonian

*δH*becomes

We now have all the ingredients ready to implement the GTB model. The Madelung integrals *B*_{IlmJ′l′m′} and density moments *Q _{Ilm}* are computed on the fly, and the tight-binding problem can be solved self-consistently.

### F. Polarizability benchmarking

We seek to identify how accurate our tight binding method is compared to DFT. Our physical observables of interest are molecular polarizabilities and charge moments, and we will show them to be significantly affected by the second order Hartree correction. The molecular multipole moments shall first be formally defined.

Consider a molecule under the influence of a constant uniform electric field ** E**. The electron density will rearrange, forming net molecular dipole

*p*

_{μ}and quadrupole

*Q*

_{μν}moments given by the linear relations,

^{20}

where *α* is the 3 × 3 dipole polarizability matrix, and *a* is the 3 × 3 × 3 dipole-quadrupole polarizability tensor. The static dipole and quadrupole moments are, respectively, given by $ p \mu 0 $ and $ Q \mu \nu 0 $. The moments are computed by integration over the total charge density *n*(** r**),

where $ r \mu \u2032 = r \mu \u2212 r \mu 0 $, with $ r \mu 0 $ being the chosen multipole origin. This definition is similar to that of the atomically expanded charge density moments (8); however, we now compute the total charge moments of the molecule. The dipole moment is independent of origin provided the molecule has no net charge, which is the case here. The quadrupole moments however do depend on the choice of origin if a net dipole moment is present.^{21}

We now turn our attention to generating the reference data. The DFT dipole polarizabilities *α*_{μν} are computed by linear response over a range of basis sets, using the correlation-consistent series cc-pVXZ^{22} and the diffuse-function augmented series aug-cc-pXDZ,^{22} where X=(D,T,Q) stands for double-, triple-, and quadruple-zeta. The DFT dipole polarizability simulations were carried out with the PBE exchange-correlation functional^{4} and the Resolution of the Identity approximation (RI-J)^{14} (with the exception of H_{2}) using the software orca.^{23} The geometries of the benchmark molecules are optimized in PBE-DFT using the Gaussian 09 software^{24} with the cc-pVTZ basis set.

. | cc-pVXZ . | aug-cc-pVXZ . | . | ||||
---|---|---|---|---|---|---|---|

. | D . | T . | Q . | D . | T . | Q . | Expt. . |

Acetylene | 2.39 | 2.92 | 3.20 | 3.47 | 3.56 | 3.58 | 3.49^{25} |

Azulene | 17.0 | 18.38 | 19.11 | 19.85 | 19.89 | 19.89 | 15.52^{26} |

Benzene | 8.69 | 9.53 | 9.98 | 10.48 | 10.52 | 10.52 | 9.9^{27} |

Butadiene | 7.08 | 7.83 | 8.23 | 8.70 | 8.74 | 8.75 | 8.10^{28} |

Butane | 7.14 | 7.88 | 8.14 | 8.30 | 8.33 | 8.33 | 7.69^{25} |

Ethane | 3.68 | 4.16 | 4.32 | 4.45 | 4.48 | 4.48 | 4.23^{25} |

Ethylene | 3.21 | 3.66 | 3.90 | 4.20 | 4.25 | 4.26 | 4.19^{25} |

H_{2} | 0.45 | 0.62 | 0.70 | 0.84 | 0.86 | 0.86 | 0.79^{25} |

Methane | 1.99 | 2.33 | 2.46 | 2.60 | 2.62 | 2.63 | 2.45^{25} |

Naphthalene | 15.71 | 16.98 | 17.65 | 18.31 | 18.35 | 18.40 | 16.6^{27} |

Propane | 5.41 | 6.02 | 6.23 | 6.37 | 6.39 | 6.40 | 5.92^{25} |

Δα | −0.15 | −0.04 | 0.02 | 0.08 | 0.09 | 0.09 | |

|Δα| | 0.17 | 0.08 | 0.07 | 0.08 | 0.09 | 0.09 |

. | cc-pVXZ . | aug-cc-pVXZ . | . | ||||
---|---|---|---|---|---|---|---|

. | D . | T . | Q . | D . | T . | Q . | Expt. . |

Acetylene | 2.39 | 2.92 | 3.20 | 3.47 | 3.56 | 3.58 | 3.49^{25} |

Azulene | 17.0 | 18.38 | 19.11 | 19.85 | 19.89 | 19.89 | 15.52^{26} |

Benzene | 8.69 | 9.53 | 9.98 | 10.48 | 10.52 | 10.52 | 9.9^{27} |

Butadiene | 7.08 | 7.83 | 8.23 | 8.70 | 8.74 | 8.75 | 8.10^{28} |

Butane | 7.14 | 7.88 | 8.14 | 8.30 | 8.33 | 8.33 | 7.69^{25} |

Ethane | 3.68 | 4.16 | 4.32 | 4.45 | 4.48 | 4.48 | 4.23^{25} |

Ethylene | 3.21 | 3.66 | 3.90 | 4.20 | 4.25 | 4.26 | 4.19^{25} |

H_{2} | 0.45 | 0.62 | 0.70 | 0.84 | 0.86 | 0.86 | 0.79^{25} |

Methane | 1.99 | 2.33 | 2.46 | 2.60 | 2.62 | 2.63 | 2.45^{25} |

Naphthalene | 15.71 | 16.98 | 17.65 | 18.31 | 18.35 | 18.40 | 16.6^{27} |

Propane | 5.41 | 6.02 | 6.23 | 6.37 | 6.39 | 6.40 | 5.92^{25} |

Δα | −0.15 | −0.04 | 0.02 | 0.08 | 0.09 | 0.09 | |

|Δα| | 0.17 | 0.08 | 0.07 | 0.08 | 0.09 | 0.09 |

The mean polarizabilities *α _{m}* = Tr(

*α*)/3 are listed in Table I as polarizability volumes

*α*/(4

*π*ε

_{0}) along with experimental values. The errors relative to experiment are given by the mean relative error $\Delta \alpha m =>"\; xmlns="">({\alpha}_{m}^{\mathit{sim}}-{\alpha}_{m}^{\mathit{exp}})/{\alpha}_{m}^{\mathit{exp}}$ and the mean relative absolute error $|\Delta \alpha |=>"\; xmlns="">|({\alpha}_{m}^{\mathit{sim}}-{\alpha}_{m}^{\mathit{exp}})/{\alpha}_{m}^{\mathit{exp}}|$. The addition of diffuse functions has a significant effect on convergence; however, the converged values overestimate polarizability considerably. This is a known issue for pure DFT.

^{29}

Next we apply external fields up to ±0.0039 E_{h}/ea_{0} to each molecule in Gaussian tight-binding using the plato ^{17} software. A linear fit is applied to the variation of dipole moments *q*_{μ} with external field strength *E*_{μ}, and the polarizabilities *α*_{μμ} extracted from the slope. We run the simulations with a minimal basis set (SV) and with the addition of polarization orbitals (SVP). For H atoms we include p orbitals, while for C we include d orbitals. The details of the method used to construct these models will be described in a separate paper. Table II shows the mean polarizabilities and their error relative to experimental data. The experimental polarizabilities are static^{25,26,28} and dynamic.^{27}

. | GTB-SV, order . | GTB-SVP, order . | . | ||||
---|---|---|---|---|---|---|---|

. | 0 . | 1 . | 2 . | 0 . | 1 . | 2 . | Expt. . |

Acetylene | 1.35 | 1.50 | 1.47 | 1.91 | 3.10 | 3.04 | 3.49^{25} |

Azulene | 12.31 | 11.87 | 11.92 | 17.17 | 17.57 | 17.58 | 15.52^{26} |

Benzene | 5.80 | 5.34 | 5.32 | 9.06 | 9.18 | 9.13 | 9.9^{27} |

Butadiene | 4.76 | 4.43 | 4.40 | 7.28 | 7.75 | 7.71 | 8.10^{28} |

Butane | 4.39 | 2.80 | 2.81 | 9.05 | 7.20 | 7.22 | 7.69^{25} |

Ethane | 2.25 | 1.45 | 1.46 | 4.50 | 3.92 | 3.94 | 4.23^{25} |

Ethylene | 1.96 | 1.73 | 1.72 | 3.27 | 3.71 | 3.68 | 4.19^{25} |

H_{2} | 0.18 | 0.18 | 0.18 | 0.20 | 0.90 | 0.90 | 0.79^{25} |

Methane | 1.16 | 0.75 | 0.77 | 2.28 | 2.28 | 2.30 | 2.45^{25} |

Naphthalene | 11.01 | 10.27 | 10.30 | 16.43 | 16.11 | 16.12 | 16.6^{27} |

Propane | 3.31 | 2.12 | 2.13 | 6.68 | 5.54 | 5.56 | 5.92^{25} |

Δα | −0.47 | −0.55 | −0.55 | −0.11 | −0.03 | −0.04 | |

|Δα| | 0.47 | 0.55 | 0.55 | 0.20 | 0.08 | 0.08 |

. | GTB-SV, order . | GTB-SVP, order . | . | ||||
---|---|---|---|---|---|---|---|

. | 0 . | 1 . | 2 . | 0 . | 1 . | 2 . | Expt. . |

Acetylene | 1.35 | 1.50 | 1.47 | 1.91 | 3.10 | 3.04 | 3.49^{25} |

Azulene | 12.31 | 11.87 | 11.92 | 17.17 | 17.57 | 17.58 | 15.52^{26} |

Benzene | 5.80 | 5.34 | 5.32 | 9.06 | 9.18 | 9.13 | 9.9^{27} |

Butadiene | 4.76 | 4.43 | 4.40 | 7.28 | 7.75 | 7.71 | 8.10^{28} |

Butane | 4.39 | 2.80 | 2.81 | 9.05 | 7.20 | 7.22 | 7.69^{25} |

Ethane | 2.25 | 1.45 | 1.46 | 4.50 | 3.92 | 3.94 | 4.23^{25} |

Ethylene | 1.96 | 1.73 | 1.72 | 3.27 | 3.71 | 3.68 | 4.19^{25} |

H_{2} | 0.18 | 0.18 | 0.18 | 0.20 | 0.90 | 0.90 | 0.79^{25} |

Methane | 1.16 | 0.75 | 0.77 | 2.28 | 2.28 | 2.30 | 2.45^{25} |

Naphthalene | 11.01 | 10.27 | 10.30 | 16.43 | 16.11 | 16.12 | 16.6^{27} |

Propane | 3.31 | 2.12 | 2.13 | 6.68 | 5.54 | 5.56 | 5.92^{25} |

Δα | −0.47 | −0.55 | −0.55 | −0.11 | −0.03 | −0.04 | |

|Δα| | 0.47 | 0.55 | 0.55 | 0.20 | 0.08 | 0.08 |

Standard monopole-order tight binding with the minimal basis (SV) severely underestimates polarizabilities by 50% on average, which is corrected by the addition of polarization orbitals (SVP). Including the dipole order improves the agreement to the point where GTB performs comparably to DFT with the cc-pVTZ basis set. Additional quadrupole terms do not notably influence the polarizabilities. The large error in the monopole-order polarizability can be attributed to out-of-plane polarizabilities *α*_{⊥} in planar molecules, for instance, in benzene or azulene. In these cases the polarizability is determined by the atomic polarization, which is not included in the monopole-order.

The same fitting procedure is applied to compute the dipole-quadrupole polarizability volumes *a*_{ξμν}/(4*π*ε_{0}). The DFT dipole-quadrupole polarizability simulations were carried out with Gaussian 09,^{24} also using the finite-difference method. Table III refers to the errors in *a*_{ξμν} between GTB and DFT. Only tensor elements larger than 0.5 Å^{4} are considered for the errors, so as not to include large relative errors from quadrupole moments that vary only little with field strength. As expected, the GTB model with SVP at dipole order performs comparably to DFT with cc-pVTZ.

. | GTB-SV, order . | GTB-SVP, order . | ||||
---|---|---|---|---|---|---|

. | 0 . | 1 . | 2 . | 0 . | 1 . | 2 . |

Errors relative to cc-pVTZ DFT | ||||||

Δα | 0.39 | 0.44 | 0.43 | 0.03 | 0.03 | 0.02 |

|Δα| | 0.40 | 0.53 | 0.55 | 0.33 | 0.06 | 0.08 |

Errors relative to aug-cc-pVTZ DFT | ||||||

Δα | 0.42 | 0.46 | 0.46 | 0.13 | 0.09 | 0.09 |

|Δα| | 0.46 | 0.60 | 0.60 | 0.29 | 0.15 | 0.16 |

. | GTB-SV, order . | GTB-SVP, order . | ||||
---|---|---|---|---|---|---|

. | 0 . | 1 . | 2 . | 0 . | 1 . | 2 . |

Errors relative to cc-pVTZ DFT | ||||||

Δα | 0.39 | 0.44 | 0.43 | 0.03 | 0.03 | 0.02 |

|Δα| | 0.40 | 0.53 | 0.55 | 0.33 | 0.06 | 0.08 |

Errors relative to aug-cc-pVTZ DFT | ||||||

Δα | 0.42 | 0.46 | 0.46 | 0.13 | 0.09 | 0.09 |

|Δα| | 0.46 | 0.60 | 0.60 | 0.29 | 0.15 | 0.16 |

The benchmarks indicate that GTB describes electronic response to external fields with an accuracy comparable to PBE-DFT. This is particularly satisfying due to the simplicity of multipole expansion employed in this method. Inclusion of the quadrupole terms does not generally improve the description of the dipole and dipole-quadrupole polarizability, hence we shall proceed with order-1 GTB for the next set of simulations.

## III. INDUCED FIELDS IN AMORPHOUS P3HT

The performance benefits of GTB make it an efficient choice for computing the properties of large molecules. We demonstrate this on the example of an amorphous bulk polymer under the influence of an external field. This system is especially challenging to model with electronic structure methods due to the polymers’ large molecular size as well as their close proximity to one another. We aim to estimate the total electric fields at each polymer site for a given point in time, using solely input from isolated strand GTB calculations in conjunction with structural data from classical molecular dynamics. This will enable future time-dependent simulations of molecular fragments in external fields to take into account the electrostatic screening by the neglected surrounding material.

With future applications in mind, we choose the polymer poly(3-hexylthiophene-2,5-diyl), P3HT, see Figure 1. The polymer strands considered here are 20 monomers long, and as such contain 500 atoms each. 160 P3HT strands were randomly packed using packmol ^{30} in a simulation box with a density well below the experimentally reported density of P3HT. This large cell of 160 P3HT strands was equilibrated using molecular dynamics (MD) simulations. P3HT is described using the force field developed by Moreno *et al.*^{31} MD simulations were performed using the gromacs-4.6.5 package.^{32–35} Periodic boundary conditions were applied in all directions. The equilibration protocol is as follows. First, the system was equilibrated successively in the NVT and NPT ensemble at 600 K to create an amorphous P3HT melt. Then, the system was quenched in steps of 50 K and equilibrated at each temperature, including the temperature of interest here. The densities extracted for each temperature are comparable with experimental data and a change of slope arises around 200 K, which corresponds to the reported glass transition of P3HT. Finally, an acquisition run was performed for 30 ns in the NPT ensemble. The complete details of the MD simulations can be found elsewhere.^{36} Snapshots of the system have been extracted every 1000 ps. Strands crossing the periodic boundaries have been reconstructed to form whole chains. We have shown previously that these simulations can reproduce quasi-elastic neutron scattering (QENS) data for P3HT at various temperatures.^{36} The experimental signal is likely to be dominated by the amorphous phase of P3HT, therefore the resulting structures are considered representatives of the amorphous system.

The polarizability of every whole chain in isolation is computed in dipole-order GTB with the SVP basis set by applying external fields of ±0.0005 E_{h}/ea_{0}, as outlined in Section II F. The polarizabilities, static dipole moments (Figure 2), and geometries act as input for a damped Polarizable-Point-Dipoles (PPD) model,^{37–39} which approximates the polymer strands as classical polarizable particles. This enables the efficient calculation of the induced fields at any point of the system. We shall briefly outline the PPD model as used in this approach.

### A. Polarizable-point-dipoles model

Consider an interacting polymer strand on site *i* to have a mean ionic position *R*_{i}, a dipole moment $ p i = p i s + p i i $ given by the sum of static $ p i s $ and induced $ p i i $ dipole moments, and an anisotropic polarizability tensor $ \alpha \u02c6 i $. The strands are considered as charge neutral, hence only intermolecular dipole interactions are taken into account. The total electrostatic energy *U _{ele}* is then given as

The first term is the interaction between total dipole moments on different strands; the second term is the energy associated with polarizing each strand; the final term is the interaction of the total dipoles with the external field. The sum is taken over all interacting dipoles in the system. Here $ T \u02c6 i j $ is the 3 × 3 electrostatic dipole-dipole interaction tensor, defined by the second-order derivative of the damping function *s*_{0}(** r**),

^{39}

which depends on the normalized charge distribution *n _{i}* of the polarizable particles,

We shall recast (26) into a more convenient form, defining the matrices

where T_{ii} = 0. The site-dependent vectors are further written into single columns, e.g., $P= p 1 p 2 \u2026 p N T $. The total electrostatic energy is then given as

The total dipole moments P are found by the stationary condition ∂_{P}*U _{ele}* = 0, which follows from minimizing

*U*, and which gives

_{ele} The super-polarizability matrix^{37} of the simulation cell is hence identified as $ \alpha \u02c6 t o t = \u2211 i j J \u2212 1 i j $, from which the dielectric constant can be derived. We can now formulate an expression for the site-dependent internal fields. Consider again each dipole to be given by a static and induced component. The induced component on site *i* is determined by the total internal field on site *i*, following the relation $ p i i = \alpha \u02c6 i E i t o t $. Thus the internal fields are given as E^{tot} = AP^{i}:

where 𝕀 is the NxN identity matrix. The total fields hence depend linearly on the applied external field, with an initial value determined by the static dipoles.

The choice of *n _{i}*(

**) needs to be addressed next. The PPD model makes the assumption that the charge distributions of the interacting particles are of small spatial extent compared to their distance to one another. Violation of this condition can lead to a break-down of the model; this is called the “polarization catastrophe.”**

*r*^{37}In this case the J-matrix is not positive definite, which means

*U*does not have a minimum solution.

_{ele}^{37}The damping function

*s*

_{0}(

**) acts to reduce the overlap between close particles by assigning a localized charge distribution**

*r**n*(

_{i}**) to them, subsequently lowering their contribution to the electrostatic energy. The optimal choice of the damping function has been, and still is, the extensive subject of research.**

*r*^{37,40,41}In this work we choose a Gaussian PPD model

^{38,39,42}for its ease of implementation and its simple reciprocal representation for periodic boundary conditions,

P3HT polymer strands have a particularly anisotropic shape, which makes this system very sensitive to a breakdown of the dipole approximation. We shall therefore test two models that avoid this problem.

#### 1. 1G model

Each strand is described by a single Gaussian-damped polarizable particle (Figure 3). A very generous smearing of at least $\kappa =0.001/ a 0 2 $ was required to keep the J-matrix positive definite. This corresponds to a Gaussian full-width half-maximum (FWHM) of 52.7 *a*_{0}, or 30% of the simulation box length.

#### 2. NG model

Each polymer strand is split into *N*_{ν} individually polarizable beads, with the strand polarizability $ \alpha \u02c6 i $ and static dipole $ p i s $ evenly divided amongst the beads (Figure 4). Beads within the same strand do not interact with each other, though they do interact with dipoles on other strands. The Gaussians are centered on the sulfur atoms at *R*_{ν}, as these follow the polymer backbone well. Using this approach, a damping of $\kappa =0.0022/ a 0 2 $ was found to give stable results, which corresponds to a FWHM of 35.5 *a*_{0}, or 20% of the simulation box length.

### B. Simulation of internal fields

The damped PPD models are run with xyz-periodic boundary conditions, evaluating the $ T \u02c6 i j $ matrix elements by Ewald summation with “tinfoil” surface conditions (ε_{Surface} → ∞).^{38} The dielectric constant ε_{r} is found as 3.32 for 1G and 3.41 for NG. This is consistent with ε_{r} ≈ 3 chosen for P3HT models^{43} and ε_{r} ≈ 2.9 for film measurements,^{44} considering the simple form of the electrostatic model. DFTB has a tendency to overestimate polarizabilities of large conjugated chain molecules, similar to DFT.^{45}

The total internal field equation (32) can furthermore be cast into self- and externally induced components

where

is the self-induced internal field produced by the system’s static dipoles, and

is defined as the screening matrix, which determines both magnitude and direction of the externally induced field, $ \gamma \u02c6 i E ext $. The screening matrix is generally asymmetric and anisotropic. We present histograms (Figure 5) of the angular-averaged screening

which measures the screening of the magnitude of an external unit field averaged over all angular directions, and the angular screening anisotropy

which is the standard-deviation of the angular-averaged screening and as such measures the anisotropy of the screening with respect to the direction of an external unit field.

We find the screening properties to vary considerably over the polymer sites (Figure 5). The 1G model yields more narrow distributions which is consistent with the use of wide Gaussian densities. Interactions between close strands are damped strongly, subsequently less site-dependence is expected. The NG model in contrast retains more interaction with the immediate environment of a strand, hence the screening properties follow a more irregular distribution. This is no surprise as the isolated strands show very strong variations in polarizability (Figure 2). A key result is that both expectation values, $E "> \gamma \u02c6 i \Omega $ and $E \sigma \gamma \u02c6 i $, are very similar between the 1G and NG models. The statistically averaged screening behavior is thus not particularly sensitive to the choice of the model; the NG model however differentiates site-resolved properties better.

For illustrative purposes we select a representative P3HT site *j* with $ "> \gamma \u02c6 j \mathrm{\Omega}$ and *σ*_{γj} = 0.08 close to the expectation values. A strong external field of 0.1 (in E_{h}/ea_{0}) leads to an angular-averaged externally induced field of 0.027, with directionally dependent variations of *σ* = 0.08. The self-induced fields in this case are negligible, as their magnitude reaches at most 0.001 (Figure 5). The total induced field at the representative P3HT site is listed in Table IV over various field strengths in $ e \u02c6 x $ direction for the NG model. A key observation is that the internal field direction does not fully align with the external field direction, but rather converges to a site-dependent constant angle for large fields. This effect is a direct consequence of the electrostatic coupling between polymer strands and is reflected in the anisotropy of the screening-matrices. Figure 6 shows the self- and externally induced fields for a slice through the simulation cell.

$ E ext x $ . | $ E j tot $ . | $\u2221 e \u02c6 x , E j tot $ . | ||
---|---|---|---|---|

0.000 | −0.000 17 | 0.000 34 | −0.000 19 | 114.00 |

0.001 | 0.000 10 | 0.000 33 | −0.000 17 | 75.43 |

0.005 | 0.001 18 | 0.000 28 | −0.000 08 | 13.98 |

0.010 | 0.002 54 | 0.000 22 | 0.000 20 | 5.10 |

0.050 | 0.013 39 | −0.000 23 | 0.000 88 | 3.88 |

0.100 | 0.026 95 | −0.000 80 | 0.001 94 | 4.47 |

$ E ext x $ . | $ E j tot $ . | $\u2221 e \u02c6 x , E j tot $ . | ||
---|---|---|---|---|

0.000 | −0.000 17 | 0.000 34 | −0.000 19 | 114.00 |

0.001 | 0.000 10 | 0.000 33 | −0.000 17 | 75.43 |

0.005 | 0.001 18 | 0.000 28 | −0.000 08 | 13.98 |

0.010 | 0.002 54 | 0.000 22 | 0.000 20 | 5.10 |

0.050 | 0.013 39 | −0.000 23 | 0.000 88 | 3.88 |

0.100 | 0.026 95 | −0.000 80 | 0.001 94 | 4.47 |

In summary, the approach taken here gives us an estimate on the internal field of a molecule in the midst of an extended system in a static external field. The site-resolved internal fields are found by using solely MD structural data and GTB data of isolated strands as an input. The screening matrices and the self-induced fields are an intuitive way to characterize the intermolecular interactions, allowing a basic description of the electrostatic environment in isolated strand simulations. More accurate values can likely be found by using a coarse-grained charge-equilibration model,^{46} though its parametrisation would require considerably more effort than the straightforward NG model outlined in this work. The PPD model can be further extended to include time-dependent external fields by using frequency-dependent polarizabilities from GTB; this would enable simulation of the extended system in conditions similar to those encountered in laser experiments.

### C. Error analysis

We would like to emphasize that LDA and GGA-DFT are known to severely overestimate the polarizability of extended conjugated systems, such as polythiophene,^{47} due to the self-interaction error.^{47} GTB, as an approximate PBE-DFT method, inherits this limitation.^{48} We demonstrate this by computing the polarizabilities of a P3HT strand for chain lengths of varying monomer count *n*, where *n* = 20 corresponds to a full P3HT chain. The chain segments are terminated with hydrogen. Our methods of choice are MP2 in the RIJCOSX approximation^{49} and PBE-DFT using the 6-31G* basis-set,^{50} and GTB-1 using the SVP basis-set. We monitor the relative error *r*_{α} in the isotropic polarizabilities of GTB and DFT to those of MP2 over monomer count *n*.

It is found that while GTB starts with a relative error of 11% for *n* = 1, the error linearly increases up to 25% for *n* = 10. We extrapolate the GTB polarizability error by a linear fit to 44% for *n* = 20. The same behavior is seen for PBE-DFT, where *r*_{α} = 6% for *n* = 1, 33% for *n* = 10, and extrapolate *r*_{α} = 62% for *n* = 20. This error is less severe than previously reported values^{47} of 100% for polythiophene at *n* = 8, possibly due to the non-planar structure of the polymer strand used here and the addition of propane side-chains.

Next we assess the effect of this error on the computed bulk properties of P3HT. We dampen the GTB strand polarizabilities by a factor of 2/3 (≈1/1.44). The 1G (NG) model now predicts a dielectric constant of 2.98 (3.00), and an angular-averaged mean screening of 0.26 (0.25). We note that the screening remains mostly unchanged. A smaller Gaussian density smearing can be used in the PPD model as a consequence of the polarizability reduction, which increases the dielectric constant and therefore largely negates the previous changes. This makes clear the need for future research on efficient self-interaction correction methods^{48,51,52} for tight-binding in conjunction with the electrostatic multipoles model. We refer to the recently developed long-range corrected DFTB method^{48} as a possible solution.

## IV. CONCLUSION

We have presented the GTB method which is a polarizable-ion tight-binding theory systematically derived from DFT. The GTB method improves upon standard self-consistent TB by a more accurate estimation of the second-order variation in the Hartree energy. This is achieved by an on-the-fly multipole expansion of the atomic densities into Gaussian-type functions, and only requires a single adjustable parameter for the atomic hardness per atomic species. GTB computes molecular polarizabilities with PBE-DFT accuracy, but with considerably less computational time due to the usual tight-binding approximations.

We have applied GTB in conjunction with classical molecular dynamics and the polarizable-point-dipoles model to give an estimate of the internal fields in bulk amorphous P3HT as a function of an applied external field. We found a strong site-dependence in the field screening strength, which is consistent with the large variety of isolated strand polarizabilities displayed. The field screening serves as a useful input in future dynamical GTB simulations to correct for the missing electrostatic environment.

The GTB method delivers reliable electrostatic polarizabilities for small systems; however as in PBE-DFT, the polarizability of conjugated extended systems is strongly overestimated. This presents the need to implement approximate self-interaction corrections in GTB for a reliably description of electrostatics on all length-scales.

## Acknowledgments

Effort sponsored by the Air Force Office of Scientific Research, Air Force Material Command, USAF, under grant number FA8655-12-1-2105. The U.S. Government is authorized to reproduce and distribute reprints for Governmental purpose notwithstanding any copyright notation thereon. M.B. was supported by the CDT in Theory and Simulation of Materials at Imperial College London funded by EPSRC Grant No. EP/G036888/1, the Thomas Young Centre under Grant No. TYC-101. A.A.Y.G. was supported by the EPSRC Supergen hub (Grant No. EP/J017361) and EPSRC Grant No. EP/K029843. Useful conversations with Aaron Thong and Peter Haynes are gratefully acknowledged. The simulation input and output files for this article can be accessed on Zenodo (DOI:10.5281/zenodo.154180) under the Creative Commons Attribution 4.0 license.