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).

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 1014 W/cm2.

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 PBE4 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 DFTB7–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.

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 

E = n bands f n ψ ( n ) | H ˆ [ ρ ] | ψ ( n ) + E dc [ ρ ] + E ion .
(1)

Here ρ is the electron number density, ψ(n) the eigenstate of level n with occupancy fn, H ˆ [ ρ ] the Kohn-Sham Hamiltonian, Eion the ion-ion repulsion, and Edc the double-counting term,

E dc [ ρ ] = E H a [ ρ ] + E x c [ ρ ] ( E x c [ ρ ] ρ ) ρ ( r ) d r .
(2)

EHa is the Hartree energy and Exc the exchange-correlation functional. We now define a reference density ρ0 as the superposition of spherical, neutral atomic densities ρ 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,

E = n f n ψ ( n ) | H ˆ [ ρ 0 ] | ψ ( n ) + E dc [ ρ 0 ] + E ion + 1 2 δ ρ ( r ) δ ρ ( r ) 1 r r f x c d r d r + O ( δ ρ 3 ) .
(3)

Here fxc 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 

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

δ E H a = 1 2 δ ρ ( r ) δ ρ ( r ) r r d r d r ,
(4)

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 Klm to avoid using a complex basis set,

(5)
K l , 2 m 1 = 2 π 2 l + 1 r l ( 1 ) m Y l m + Y l m ,
K l , 2 m = i 2 π 2 l + 1 r l ( 1 ) m Y l m Y l m ,
where Y l m are spherical harmonics of order l and degree m, with −lml. We have chosen the principal axis to lie in the z-direction, hence r = (K10K11K12)T = (z x y)T.

The density variation δρ can be written as a superposition of localized atomically centered components δ ρ = I δ ρ 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,

δ ρ I ( r ) = l m a I l m ( r R I ) K l m ( r R I ) .
(6)

We wish to express aIlm(r) by a set of radial functions fIlm(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 QIlm as the expansion coefficients of the radial functions fIlm(r),

δ ρ I ( r ) = l m Q I l m f I l m ( r R I ) K l m ( r R I ) ,
(7)

where the charge density moments QIlm are defined in terms of the cubic harmonics,

Q I l m = δ ρ I ( r ) K l m ( r ) d r .
(8)

We find a normalisation condition for fIlm by substituting the variational density (7) into the density moments (8),

f I l m ( r ) | K l m ( r ) | 2 d r = 1 .
(9)

The variational density (7) is now substituted into the second order correction of the Hartree energy (4), with RIJ = RJRI,

δ E H a = 1 2 I l m J l m Q I l m Q J l m f I l m ( r ) K l m ( r ) f J l m ( r R I J ) K l m ( r R I J ) r r d r d r ,
(10)

where the integral can be identified with the Madelung structure constant BIlm,Jlm(RIJ), leading to the compact expression,

δ E H a = 1 2 I l m J l m Q I l m Q J l m B I l m , J l m ( R I J ) .
(11)

This formalism enables us to express δEHa as a function of the charge density moments QIlm(8). The total energy can hence be self-consistently minimized with respect to the moments QIlm.

The second order correction to the Hartree energy δEHa depends directly on the electron density ρ via the density moments QIlm(8). Solving for ρ 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 ϕ attached to it. The eigenstate ψ(n) then expands into ψ ( n ) ( r ) = I α c I α ( n ) ϕ I α ( r R I ) , where c I α ( n ) are the expansion coefficients. Assuming real orbitals, the total electron density is given by

ρ ( r ) = I α J β ρ I α J β ϕ I α ( r R I ) ϕ J β ( r R J ) ,
(12)

with the density matrix

ρ I α J β = n f n c I α ( n ) c J β ( n ) ,
(13)

where fn is the electron occupancy of state n. Substituting this definition of the density (12) into (8), we can express the charge moments QIlm in terms of the density matrix

Q I l m = Z I δ l 0 δ m 0 + α J β ρ I α J β M I α J β l m ,
(14)

where we define ZI as the ionic charge of site I and M I α J β l m as the cubic moment integral over the basis functions ϕ and ϕ,

M I α J β l m ( R I J ) = ϕ I α ( r ) ϕ J β ( r R I J ) K l m ( r ) d r .
(15)

Here we expressed the density variation as δ ρ I ( r ) = ρ I ( r ) ρ I 0 ( r ) . The atomic reference density ρ I 0 ( r ) is spherically symmetric, therefore ρ I 0 ( r R I ) K l m ( r R I ) = Z I δ l 0 δ m 0 .

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

J β δ H I α J β c J β ( n ) = ( δ E H a ) c I α ( n ) ,
(16)

giving

δ H I α J β = 1 4 K l m l m Q K l m ( M J β I α l m B J l m K l m + B K l m J l m + M I α J β l m B I l m K l m + B K l m I l m ) .
(17)

We control the level of accuracy in δEHa by varying the order of the multipole expansion. Consider a zero-order expansion, where δρI 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).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 δEHa in a systematically improvable framework. A first-order multipole expansion, for example, considers δρI to be polarizable, leading to an overall improved picture of electrostatic interaction.

We have not yet chosen a form for our radial functions fIlm(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 fIlm(r) matches the true radial distribution aIlm(r). Past work in the monopole approximation has shown a single Gaussian function for fIlm(r) to perform well;16 this motivates us to pursue this approach here,

f I l m ( r ) = k I l m e α I l m r 2
(18)

wherer kIlm is the normalisation constant introduced to satisfy Equation (9),

k I 00 = α I 00 π 3 / 2 , k I 1 m = 2 α I 1 m α I 1 m π 3 / 2 , k I 2 m = 4 3 α I 2 m 2 α I 2 m π 3 / 2 ,
(19)

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 αI over all momenta for each atomic species. This width can be associated with the Hubbard parameter UI of the atomic species as shown by Elstner et al.2 We follow this approach here, using

α I = π 8 U I 2 .
(20)

While this value can be calculated from first principles, we leave it as an adjustable parameter and use UI from the mio-1-1 parameter set,2 averaged over all angular momenta. The GTB model is implemented in the tight-binding code plato17 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.

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 H0, the self-consistent contribution δH, and the external field contribution δHext,

H ˆ = H ˆ 0 + δ H ˆ + δ H ˆ e x t .
(21)

The external field coupling energy is found by integration over the total charge density n(r),

E e x t = E ( r r 0 ) n ( r ) d r = I ν ( E ) ν ( Q I 00 R I r 0 ν + Q I 1 ν ) .
(22)

The monopole moment QI00 can be identified as the net charge of the total atomic density ρ on site I. The contribution δHext to the self-consistent Hamiltonian δH becomes

δ H I α J β e x t = ν ( E ) ν ( M I α J β 00 R I r 0 ν + M I α J β 1 ν ) .
(23)

We now have all the ingredients ready to implement the GTB model. The Madelung integrals BIlmJlm and density moments QIlm are computed on the fly, and the tight-binding problem can be solved self-consistently.

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 

p μ = p μ 0 + ν α μ ν E ν , Q μ ν = Q μ ν 0 + ξ a ξ μ ν E ξ ,
(24)

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 μ 0 and Q μ ν 0 . The moments are computed by integration over the total charge density n(r),

p μ = n ( r ) r μ d r , Q μ ν = 1 2 n ( r ) 3 r μ r ν δ μ ν r 2 d r ,
(25)

where r μ = r μ r μ 0 , with r μ 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-pVXZ22 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 functional4 and the Resolution of the Identity approximation (RI-J)14 (with the exception of H2) using the software orca.23 The geometries of the benchmark molecules are optimized in PBE-DFT using the Gaussian 09 software24 with the cc-pVTZ basis set.

TABLE I.

Basis set convergence of mean polarizability volumes αm = Tr(α)/3 in Å3 for various hydrocarbons in PBE-DFT. The mean relative error and mean relative absolute errors to experiment are given as Δα and |Δα|.

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.4925  
Azulene  17.0  18.38  19.11  19.85  19.89  19.89  15.5226  
Benzene  8.69  9.53  9.98  10.48  10.52  10.52  9.927  
Butadiene  7.08  7.83  8.23  8.70  8.74  8.75  8.1028  
Butane  7.14  7.88  8.14  8.30  8.33  8.33  7.6925  
Ethane  3.68  4.16  4.32  4.45  4.48  4.48  4.2325  
Ethylene  3.21  3.66  3.90  4.20  4.25  4.26  4.1925  
H2  0.45  0.62  0.70  0.84  0.86  0.86  0.7925  
Methane  1.99  2.33  2.46  2.60  2.62  2.63  2.4525  
Naphthalene  15.71  16.98  17.65  18.31  18.35  18.40  16.627  
Propane  5.41  6.02  6.23  6.37  6.39  6.40  5.9225  
Δα  −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.4925  
Azulene  17.0  18.38  19.11  19.85  19.89  19.89  15.5226  
Benzene  8.69  9.53  9.98  10.48  10.52  10.52  9.927  
Butadiene  7.08  7.83  8.23  8.70  8.74  8.75  8.1028  
Butane  7.14  7.88  8.14  8.30  8.33  8.33  7.6925  
Ethane  3.68  4.16  4.32  4.45  4.48  4.48  4.2325  
Ethylene  3.21  3.66  3.90  4.20  4.25  4.26  4.1925  
H2  0.45  0.62  0.70  0.84  0.86  0.86  0.7925  
Methane  1.99  2.33  2.46  2.60  2.62  2.63  2.4525  
Naphthalene  15.71  16.98  17.65  18.31  18.35  18.40  16.627  
Propane  5.41  6.02  6.23  6.37  6.39  6.40  5.9225  
Δα  −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 Δ α m = ( α m sim α m exp ) / α m exp and the mean relative absolute error | Δ α | = | ( α m sim α m exp ) / α m 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 Eh/ea0 to each molecule in Gaussian tight-binding using the plato17 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 static25,26,28 and dynamic.27 

TABLE II.

Gaussian tight binding calculations of mean polarizability volumes αm = Tr(α)/3 in Å3 for single-valence (SV) and polarizable (SVP) basis sets, as well as various orders of multipole expansion. The mean relative error and mean relative absolute errors to experiment are given as Δα and |Δα|.

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.4925  
Azulene  12.31  11.87  11.92  17.17  17.57  17.58  15.5226  
Benzene  5.80  5.34  5.32  9.06  9.18  9.13  9.927  
Butadiene  4.76  4.43  4.40  7.28  7.75  7.71  8.1028  
Butane  4.39  2.80  2.81  9.05  7.20  7.22  7.6925  
Ethane  2.25  1.45  1.46  4.50  3.92  3.94  4.2325  
Ethylene  1.96  1.73  1.72  3.27  3.71  3.68  4.1925  
H2  0.18  0.18  0.18  0.20  0.90  0.90  0.7925  
Methane  1.16  0.75  0.77  2.28  2.28  2.30  2.4525  
Naphthalene  11.01  10.27  10.30  16.43  16.11  16.12  16.627  
Propane  3.31  2.12  2.13  6.68  5.54  5.56  5.9225  
Δα  −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.4925  
Azulene  12.31  11.87  11.92  17.17  17.57  17.58  15.5226  
Benzene  5.80  5.34  5.32  9.06  9.18  9.13  9.927  
Butadiene  4.76  4.43  4.40  7.28  7.75  7.71  8.1028  
Butane  4.39  2.80  2.81  9.05  7.20  7.22  7.6925  
Ethane  2.25  1.45  1.46  4.50  3.92  3.94  4.2325  
Ethylene  1.96  1.73  1.72  3.27  3.71  3.68  4.1925  
H2  0.18  0.18  0.18  0.20  0.90  0.90  0.7925  
Methane  1.16  0.75  0.77  2.28  2.28  2.30  2.4525  
Naphthalene  11.01  10.27  10.30  16.43  16.11  16.12  16.627  
Propane  3.31  2.12  2.13  6.68  5.54  5.56  5.9225  
Δα  −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.

TABLE III.

Errors of GTB quadrupole-dipole polarizabilities aξμν to DFT for single-valence (SV) and polarizable (SVP) basis sets, as well as various orders of multipole expansion. The mean relative error and mean relative absolute errors are given as Δα and |Δα|.

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.

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.

FIG. 1.

A P3HT dimer. The polymers in this work were chosen to be 20 monomers long (n = 10), giving them a mass of 3.3 kDa and an end-to-end length of about 70 Å. The simulation cell is cubic with a cell length of 95.42 Å and a density of 1017 kg/m3.

FIG. 1.

A P3HT dimer. The polymers in this work were chosen to be 20 monomers long (n = 10), giving them a mass of 3.3 kDa and an end-to-end length of about 70 Å. The simulation cell is cubic with a cell length of 95.42 Å and a density of 1017 kg/m3.

Close modal

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 packmol30 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.

FIG. 2.

Electrostatic properties of isolated P3HT strands from the amorphous phase: histograms of mean polarizability volume αm (a) and polarizability anisotropy Δα (b) in 103 Å3, and static dipole moments p i s (c) in ea0. Here Δα is defined by 2Δα2 = (αxx  −  αyy)2 + (αyy  −  αzz)2 + (αzz  −  αxx)2.

FIG. 2.

Electrostatic properties of isolated P3HT strands from the amorphous phase: histograms of mean polarizability volume αm (a) and polarizability anisotropy Δα (b) in 103 Å3, and static dipole moments p i s (c) in ea0. Here Δα is defined by 2Δα2 = (αxx  −  αyy)2 + (αyy  −  αzz)2 + (αzz  −  αxx)2.

Close modal

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 Eh/ea0, 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.

Consider an interacting polymer strand on site i to have a mean ionic position Ri, 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 α ˆ i . The strands are considered as charge neutral, hence only intermolecular dipole interactions are taken into account. The total electrostatic energy Uele is then given as

U e l e = U d , d + U d , p o l + U d , e x t = i , j > i p i T ˆ i j p j + 1 2 i p i i α ˆ i 1 p i i + i p i . E e x t .
(26)

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 ˆ i j is the 3 × 3 electrostatic dipole-dipole interaction tensor, defined by the second-order derivative of the damping function s0(r),39 

T ˆ i j μ ν = r μ r ν s 0 ( r ) / r | r = R i j ,
(27)

which depends on the normalized charge distribution ni of the polarizable particles,

s 0 ( R i j ) = R i j n i ( r R i ) n j ( r R j ) r r d r d r .
(28)

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

A i j = α ˆ i 1 δ i j , T i j = T ˆ i j , J i j = A i j T i j ,
(29)

where Tii = 0. The site-dependent vectors are further written into single columns, e.g., P = p 1 p 2 p N T . The total electrostatic energy is then given as

U e l e = 1 2 P T P + 1 2 P i A P i + P E ext .
(30)

The total dipole moments P are found by the stationary condition ∂PUele = 0, which follows from minimizing Uele, and which gives

P = J 1 A P s + J 1 E e x t .
(31)

The super-polarizability matrix37 of the simulation cell is hence identified as α ˆ t o t = i j J 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 = α ˆ i E i t o t . Thus the internal fields are given as Etot = APi:

E tot = A J 1 A + I P s + A J 1 E ext ,
(32)

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 ni(r) 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.”37 In this case the J-matrix is not positive definite, which means Uele does not have a minimum solution.37 The damping function s0(r) acts to reduce the overlap between close particles by assigning a localized charge distribution ni(r) 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.37,40,41 In this work we choose a Gaussian PPD model38,39,42 for its ease of implementation and its simple reciprocal representation for periodic boundary conditions,

n i ( r ) = κ π 3 / 2 e κ r 2 .
(33)

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 κ = 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 a0, or 30% of the simulation box length.

FIG. 3.

Coordinate system in the 1G model. Each P3HT strand is described as one polarizable particle. The charge distribution of a strand is approximated by a diffuse Gaussian (patterned area).

FIG. 3.

Coordinate system in the 1G model. Each P3HT strand is described as one polarizable particle. The charge distribution of a strand is approximated by a diffuse Gaussian (patterned area).

Close modal

2. NG model

Each polymer strand is split into Nν individually polarizable beads, with the strand polarizability α ˆ 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 κ = 0 . 0022 / a 0 2 was found to give stable results, which corresponds to a FWHM of 35.5 a0, or 20% of the simulation box length.

FIG. 4.

Coordinate system in the NG model. The P3HT strands are split into beads centered on the sulfur atoms (orange). The charge distribution of a bead is approximated by a diffuse Gaussian (patterned area). The polarizability and static dipole of a strand are evenly split amongst the beads.

FIG. 4.

Coordinate system in the NG model. The P3HT strands are split into beads centered on the sulfur atoms (orange). The charge distribution of a bead is approximated by a diffuse Gaussian (patterned area). The polarizability and static dipole of a strand are evenly split amongst the beads.

Close modal

The damped PPD models are run with xyz-periodic boundary conditions, evaluating the T ˆ 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 models43 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

E i tot = E i self + γ ˆ i E ext ,
(34)

where

E i self = j A J 1 A + I i j p j s
(35)

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

γ ˆ i = j A J 1 i j
(36)

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

γ ˆ i Ω = 1 / 4 π γ ˆ i . e ˆ r ( θ , ϕ ) d Ω ,
(37)

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

σ γ ˆ i = 1 / 4 π γ ˆ i . e ˆ r ( θ , ϕ ) 2 d Ω γ i ˆ Ω 2 1 / 2 ,
(38)

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.

FIG. 5.

Internal field properties of amorphous P3HT in the 1G (blue, solid fill) and NG (red, patterned fill) models: histograms of angular-averaged screening γ ˆ i Ω (a), angular screening anisotropy σ γ ˆ i (b), and the self-induced fields E i s (c) in 10−3 Eh/ea0. The expectation values (vertical lines) in the 1G (NG) model are E γ ˆ i Ω = 0 . 29 (0.27), E σ γ ˆ i = 0 . 08 (0.09), and E E i s = 0 . 16 (0.3) 10−3 Eh/ea0.

FIG. 5.

Internal field properties of amorphous P3HT in the 1G (blue, solid fill) and NG (red, patterned fill) models: histograms of angular-averaged screening γ ˆ i Ω (a), angular screening anisotropy σ γ ˆ i (b), and the self-induced fields E i s (c) in 10−3 Eh/ea0. The expectation values (vertical lines) in the 1G (NG) model are E γ ˆ i Ω = 0 . 29 (0.27), E σ γ ˆ i = 0 . 08 (0.09), and E E i s = 0 . 16 (0.3) 10−3 Eh/ea0.

Close modal

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 γ ˆ i Ω and E σ γ ˆ 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 γ ˆ j Ω = 0 . 27 and σγj = 0.08 close to the expectation values. A strong external field of 0.1 (in Eh/ea0) 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 ˆ 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.

TABLE IV.

Total induced fields in Eh/ea0 of an example strand site j in amorphous P3HT for various external fields in e ˆ x direction. The angle between external and total internal field e ˆ x , E j t o t gradually converges to a limit of 5.10. The screening factor for this direction is found as 0.27.

E ext x E j tot e ˆ 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 e ˆ 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.

FIG. 6.

Internal fields of amorphous P3HT with applied external field (red arrows) of E ext = 0 . 002 E h / ea 0 e ˆ x and without field (black arrows) in the NG model. The dots represent the mean position of P3HT strands. The largest internal field E i tot = 0 . 000 94 E h / ea 0 is marked for reference. A slice of 20 Å width is shown for improved visibility.

FIG. 6.

Internal fields of amorphous P3HT with applied external field (red arrows) of E ext = 0 . 002 E h / ea 0 e ˆ x and without field (black arrows) in the NG model. The dots represent the mean position of P3HT strands. The largest internal field E i tot = 0 . 000 94 E h / ea 0 is marked for reference. A slice of 20 Å width is shown for improved visibility.

Close modal

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 approximation49 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 values47 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 methods48,51,52 for tight-binding in conjunction with the electrostatic multipoles model. We refer to the recently developed long-range corrected DFTB method48 as a possible solution.

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.

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.

1.
S.
Baker
,
J. S.
Robinson
,
C.
Haworth
,
H.
Teng
,
R.
Smith
,
C.
Chirilă
,
M.
Lein
,
J.
Tisch
, and
J.
Marangos
, “
Probing proton dynamics in molecules on an attosecond time scale
,”
Science
312
,
424
427
(
2006
).
2.
M.
Elstner
,
D.
Porezag
,
G.
Jungnickel
,
J.
Elsner
,
M.
Haugk
,
T.
Frauenheim
,
S.
Suhai
, and
G.
Seifert
, “
Self-consistent-charge density-functional tight-binding method for simulations of complex materials properties
,”
Phys. Rev. B
58
,
7260
(
1998
).
3.
Z.
Bodrog
and
B.
Aradi
, “
Possible improvements to the self-consistent-charges density-functional tight-binding method within the second order
,”
Phys. Status Solidi B
249
,
259
269
(
2012
).
4.
J. P.
Perdew
,
K.
Burke
, and
M.
Ernzerhof
, “
Generalized gradient approximation made simple
,”
Phys. Rev. Lett.
77
,
3865
(
1996
).
5.
M.
Finnis
,
A.
Paxton
,
M.
Methfesselt
, and
M.
van Schilfgaarde
, “
Self-consistent tight-binding approximation including polarisable ions
,” in
MRS Proceedings
(
Cambridge University Press
,
1997
), Vol.
491
, p.
265
.
6.
A.
Paxton
,
T.
Todorov
, and
A.
Elena
, “
Ring currents in azulene
,”
Chem. Phys. Lett.
483
,
154
158
(
2009
).
7.
S.
Kaminski
,
T. J.
Giese
,
M.
Gaus
,
D. M.
York
, and
M.
Elstner
, “
Extended polarization in third-order SCC-DFTB from chemical-potential equalization
,”
J. Phys. Chem. A
116
,
9131
9141
(
2012
).
8.
S.
Kaminski
,
M.
Gaus
, and
M.
Elstner
, “
Improved electronic properties from third-order SCC-DFTB with cost efficient post-scf extensions
,”
J. Phys. Chem. A
116
,
11927
11937
(
2012
).
9.
A. S.
Christensen
,
M.
Elstner
, and
Q.
Cui
, “
Improving intermolecular interactions in DFTB3 using extended polarization from chemical-potential equalization
,”
J. Chem. Phys.
143
,
084123
(
2015
).
10.
W. M. C.
Foulkes
and
R.
Haydock
, “
Tight-binding models and density-functional theory
,”
Phys. Rev. B
39
,
12520
(
1989
).
11.
C.
Goringe
,
D.
Bowler
 et al, “
Tight-binding modelling of materials
,”
Rep. Prog. Phys.
60
,
1447
(
1997
).
12.
A.
Horsfield
and
A.
Bratkovsky
, “
Ab initio tight binding
,”
J. Phys.: Condens. Matter
12
,
R1
(
2000
).
13.
A. P.
Sutton
,
M. W.
Finnis
,
D. G.
Pettifor
, and
Y.
Ohta
, “
The tight-binding bond model
,”
J. Phys. C: Solid State Phys.
21
,
35
(
1988
).
14.
C.-K.
Skylaris
,
L.
Gagliardi
,
N. C.
Handy
,
A. G.
Ioannou
,
S.
Spencer
, and
A.
Willetts
, “
On the resolution of identity coulomb energy approximation in density functional theory
,”
J. Mol. Struct.: THEOCHEM
501
,
229
239
(
2000
).
15.
M.
Elstner
, “
The SCC-DFTB method and its application to biological systems
,”
Theor. Chem. Acc.
116
,
316
325
(
2006
).
16.
P.
Soin
,
A.
Horsfield
, and
D.
Nguyen-Manh
, “
Efficient self-consistency for magnetic tight binding
,”
Comput. Phys. Commun.
182
,
1350
1360
(
2011
).
17.
S. D.
Kenny
and
A. P.
Horsfield
, “
Plato: A localised orbital based density functional theory code
,”
Comput. Phys. Commun.
180
,
2616
2621
(
2009
).
18.
T. A.
Niehaus
,
M.
Elstner
,
T.
Frauenheim
, and
S.
Suhai
, “
Application of an approximate density-functional method to sulfur containing compounds
,”
J. Mol. Struct.: THEOCHEM
541
,
185
194
(
2001
).
19.
A.
Fihey
,
C.
Hettich
,
J.
Touzeau
,
F.
Maurel
,
A.
Perrier
,
C.
Köhler
,
B.
Aradi
, and
T.
Frauenheim
, “
SCC-DFTB parameters for simulating hybrid gold-thiolates compounds
,”
J. Comput. Chem.
36
,
2075
2087
(
2015
).
20.
R. E.
Raab
and
O. L.
De Lange
,
Multipole Theory in Electromagnetism: Classical, Quantum, and Symmetry Aspects, with Applications
(
Oxford University Press
,
USA
,
2005
).
21.
J. D.
Jackson
,
Classical Electrodynamics
(
Wiley
,
1999
).
22.
T. H.
Dunning
, Jr.
, “
Gaussian basis sets for use in correlated molecular calculations. I. The atoms boron through neon and hydrogen
,”
J. Chem. Phys.
90
,
1007
1023
(
1989
).
23.
F.
Neese
, “
The orca program system
,”
Wiley Interdiscip. Rev.: Comput. Mol. Sci.
2
,
73
78
(
2012
).
24.
M.
Frisch
,
G.
Trucks
,
H. B.
Schlegel
,
G.
Scuseria
,
M.
Robb
,
J.
Cheeseman
,
G.
Scalmani
,
V.
Barone
,
B.
Mennucci
,
G.
Petersson
 et al, gaussian 09, Revision a.02, Gaussian, Inc., Wallingford, CT, 2009.
25.
T. N.
Olney
,
N.
Cann
,
G.
Cooper
, and
C.
Brion
, “
Absolute scale determination for photoabsorption spectra and the calculation of molecular properties using dipole sum-rules
,”
Chem. Phys.
223
,
59
98
(
1997
).
26.
W.
Baumann
, “
Dipole moments and polarizabilities of azulene in its ground state and first and second excited singlet state
,”
Chem. Phys.
20
,
17
24
(
1977
).
27.
C. L.
Cheng
,
D.
Murthy
, and
G.
Ritchie
, “
Molecular polarizabilities from the cotton-mouton effect
,”
Aust. J. Chem.
25
,
1301
1305
(
1972
).
28.
G.
Maroulis
,
C.
Makris
,
U.
Hohm
, and
U.
Wachsmuth
, “
Determination of the complete polarizability tensor of 1,3-butadiene by combination of refractive index and light scattering measurements and accurate quantum chemical ab initio calculations
,”
J. Phys. Chem. A
103
,
4359
4367
(
1999
).
29.
B.
Champagne
,
E. A.
Perpete
,
S. J.
van Gisbergen
,
E.-J.
Baerends
,
J. G.
Snijders
,
C.
Soubra-Ghaoui
,
K. A.
Robins
, and
B.
Kirtman
, “
Assessment of conventional density functional schemes for computing the polarizabilities and hyperpolarizabilities of conjugated oligomers: An ab initio investigation of polyacetylene chains
,”
J. Chem. Phys.
109
,
10489
10498
(
1998
).
30.
L.
Martínez
,
R.
Andrade
,
E. G.
Birgin
, and
J. M.
Martínez
, “
Packmol: A package for building initial configurations for molecular dynamics simulations
,”
J. Comput. Chem.
30
,
2157
2164
(
2009
).
31.
M.
Moreno
,
M.
Casalegno
,
G.
Raos
,
S. V.
Meille
, and
R.
Po
, “
Molecular modeling of crystalline alkylthiophene oligomers and polymers
,”
J. Phys. Chem. B
114
,
1591
1602
(
2010
).
32.
H. J.
Berendsen
,
D.
Van Der Spoel
, and
R.
Van Drunen
, “
GROMACS: A message-passing parallel molecular dynamics implementation
,”
Comput. Phys. Commun.
91
,
43
56
(
1995
).
33.
S.
Pronk
,
S.
Páll
,
R.
Schulz
,
P.
Larsson
,
P.
Bjelkmar
,
R.
Apostolov
,
M. R.
Shirts
,
J. C.
Smith
,
P. M.
Kasson
,
D.
Van Der Spoel
 et al, “
GROMACS 4.5: A high-throughput and highly parallel open source molecular simulation toolkit
,”
Bioinformatics
29
,
845
(
2013
).
34.
D.
Van Der Spoel
,
E.
Lindahl
,
B.
Hess
,
G.
Groenhof
,
A. E.
Mark
, and
H. J.
Berendsen
, “
GROMACS: Fast, flexible, and free
,”
J. Comput. Chem.
26
,
1701
1718
(
2005
).
35.
B.
Hess
,
C.
Kutzner
,
D.
Van Der Spoel
, and
E.
Lindahl
, “
GROMACS 4: Algorithms for highly efficient, load-balanced, and scalable molecular simulation
,”
J. Chem. Theory Comput.
4
,
435
447
(
2008
).
36.
A. A.
Guilbert
,
A.
Urbina
,
J.
Abad
,
C.
Díaz-Paniagua
,
F.
Batallán
,
T.
Seydel
,
M.
Zbiri
,
V.
García-Sakai
, and
J.
Nelson
, “
Temperature-dependent dynamics of polyalkylthiophene conjugated polymers: A combined neutron scattering and simulation study
,”
Chem. Mater.
27
,
7652
7661
(
2015
).
37.
B. T.
Thole
, “
Molecular polarizabilities calculated with a modified dipole interaction
,”
Chem. Phys.
59
,
341
350
(
1981
).
38.
M. P.
Allen
and
D. J.
Tildesley
,
Computer Simulation of Liquids
(
Oxford University Press
,
1989
).
39.
J.
Sala
,
E.
Guàrdia
, and
M.
Masia
, “
The polarizable point dipoles method with electrostatic damping: Implementation on a model system
,”
J. Chem. Phys.
133
,
234101
(
2010
).
40.
L.
Jensen
,
P.-O.
Åstrand
,
A.
Osted
,
J.
Kongsted
, and
K. V.
Mikkelsen
, “
Polarizability of molecular clusters as calculated by a dipole interaction model
,”
J. Chem. Phys.
116
,
4001
4010
(
2002
).
41.
I.
Harczuk
,
O.
Vahtras
, and
H.
Ågren
, “
Hyperpolarizabilities of extended molecular mechanical systems
,”
Phys. Chem. Chem. Phys.
18
,
8710
8722
(
2016
).
42.
D.
Elking
,
T.
Darden
, and
R. J.
Woods
, “
Gaussian induced dipole polarization model
,”
J. Comput. Chem.
28
,
1261
1274
(
2007
).
43.
C.
Tanase
,
E.
Meijer
,
P.
Blom
, and
D.
De Leeuw
, “
Unification of the hole transport in polymeric field-effect transistors and light-emitting diodes
,”
Phys. Rev. Lett.
91
,
216601
(
2003
).
44.
A. J.
Morfa
,
T. M.
Barnes
,
A. J.
Ferguson
,
D. H.
Levi
,
G.
Rumbles
,
K. L.
Rowlen
, and
J.
van de Lagemaat
, “
Optical characterization of pristine poly (3-hexyl thiophene) films
,”
J. Polym. Sci., B: Polym. Phys.
49
,
186
194
(
2011
).
45.
T. A.
Niehaus
and
F.
Della Sala
, “
Range separated functionals in the density functional based tight-binding method: Formalism
,”
Phys. Status Solidi B
249
,
237
244
(
2012
).
46.
A. K.
Rappe
and
W. A.
Goddard
III
, “
Charge equilibration for molecular dynamics simulations
,”
J. Phys. Chem.
95
,
3358
3363
(
1991
).
47.
P.
Mori-Sánchez
,
Q.
Wu
, and
W.
Yang
, “
Accurate polymer polarizabilities with exact exchange density-functional theory
,”
J. Chem. Phys.
119
,
11001
11004
(
2003
).
48.
V.
Lutsker
,
B.
Aradi
, and
T. A.
Niehaus
, “
Implementation and benchmark of a long-range corrected functional in the density functional based tight-binding method
,”
J. Chem. Phys.
143
,
184107
(
2015
).
49.
S.
Kossmann
and
F.
Neese
, “
Efficient structure optimization with second-order many-body perturbation theory: The RIJCOSX-MP2 method
,”
J. Chem. Theory Comput.
6
,
2325
2338
(
2010
).
50.
W. J.
Hehre
,
R.
Ditchfield
, and
J. A.
Pople
, “
Selfconsistent molecular orbital methods. XII. Further extensions of gaussiantype basis sets for use in molecular orbital studies of organic molecules
,”
J. Chem. Phys.
56
,
2257
2261
(
1972
).
51.
J.
Messud
,
P. M.
Dinh
,
P.-G.
Reinhard
, and
E.
Suraud
, “
The generalized SIC-OEP formalism and the generalized SIC-Slater approximation (stationary and time-dependent cases)
,”
Ann. Phys.
523
,
270
290
(
2011
).
52.
P.
Klüpfel
,
P. M.
Dinh
,
P.-G.
Reinhard
, and
E.
Suraud
, “
Koopmans condition in self-interaction-corrected density-functional theory
,”
Phys. Rev. A
88
,
052501
(
2013
).