17
$\begingroup$

Often empirical electrostatic models or molecular force fields approximate the molecular polarizability using an additive model, e.g., for N atoms:

$$ \alpha_{mol} = \sum_i^N \alpha_{atom_i} $$

This assumes polarizability is additive — and for many molecules, that assumption works okay. It also assumes isotropic polarizabilities. For atoms, that's a reasonable assumption because they're roughly spherical.

But molecules aren't spherical, so all the elements of the polarizability tensor matter (e.g., the anisotropy).

For a given frame of reference / coordinate system, how do you calculate an approximate $3 \times 3$ polarizability tensor from the additive $\alpha_{atom_i}$ components?

It seems like there should be an easy way (e.g., you calculate the rank-1 dipole moment from the charges and the atomic coordinates) but I haven't found a ready reference.

$\endgroup$
5
  • 2
    $\begingroup$ I believe the Thole damping paper describes how to do this, and the Applequist paper describes in a bit more detail (and is more readable in my opinion). Applequist, J., Carl, J. R., & Fung, K. K. (1972). Atom dipole interaction model for molecular polarizability. Application to polyatomic molecules and determination of atom polarizabilities. Journal of the American Chemical Society, 94(9), 2952-2960. Thole, B. T. (1981). Molecular polarizabilities calculated with a modified dipole interaction. Chemical Physics, 59(3), 341-350. $\endgroup$
    – jheindel
    Feb 7, 2022 at 3:54
  • 1
    $\begingroup$ @jheindel - sounds like a good start towards an answer if you want to write it up in some detail. (Most of these papers are a bit unclear on how the T matrix is supposed to work and how you go from the 3N x 3N matrix down to a 3x3 molecular polarizability tensor.) $\endgroup$ Feb 7, 2022 at 17:53
  • 2
    $\begingroup$ I agree they are unclear, which is why I am hesitant to write an answer :) I've constructed the 3Nx3N matrix quite a few times when implementing various force fields but I've never explicitly constructed the 3x3 molecular polarizability tensor. I'll look into it and see if I can understand it well enough to write an answer. $\endgroup$
    – jheindel
    Feb 7, 2022 at 20:08
  • 2
    $\begingroup$ "But molecules aren't spherical, so all the elements of the polarizability tensor matter (e.g., the anisotropy)." Regarding this sentence... The off-diagonal elements of the polarizability tensor are very low in magnitude. $\endgroup$
    – Pro
    Feb 8, 2022 at 12:27
  • 4
    $\begingroup$ @PrasantaBandyopadhyay I believe that's a question of the coordinate system. If your coordinates happen to align with the eigenvectors of the tensor, then you have zero off-diagonal elements. But, if you rotate the coordinate system by 45 degrees, then the off-diagonal elements will have similar magnitude to the diagonal elements. $\endgroup$ Feb 8, 2022 at 13:16

1 Answer 1

7
+200
$\begingroup$

Luckily, I had a need to do this recently, so I can answer with more confidence than I had when commenting before.

Basically, the procedure is what you described. We need to: (1) compute the interaction of our molecule in an applied electric field, (2) compute the change in dipole moment in this field, (3) compute the polarizability elements by finite difference.

Let's leave aside step (1) for a second since it is slightly different depending on the model of polarization being used.

The central equation needed for steps (2) and (3) is the following:

$$ \mathbf{P}^{ind}=\alpha\mathbf{E} \tag{1} $$

or, to make things simpler, let's only consider applying a field in each direction $w=x,y,z$ separately. Then, we have,

$$ \mathbf{P}^{\textrm{ind}}=\alpha\mathbf{E}_w \tag{2} $$

Then, the induced dipole moment may still have a response in any direction, $q=x,y,z$ due to the anisotropy of the polarizability. Taking the derivative with respect to $\mathbf{E}_w$ lets us write the following:

$$ \frac{dP^{ind}_{q}}{dE_w}=\alpha_{qw} \tag{3} $$

Now that we have an expression for the change of each element of the induced dipole $P^{ind}_{q}$ when placed in a field $E_w$, we can compute each polarizability element by finite difference.

$$ \alpha_{qw} = \frac{P^{ind}_q(E_{+w}) - P^{ind}_q(E_{-w})}{2|E_w|} \tag{4} $$

where $P^{ind}_q(E_{+w})$ is the $q$-component of the induced dipole when experiencing a field coming from the positive-$w$ direction. We then repeat this procedure for each of cartesian axes to get all nine elements of the molecular polarizability.

This procedure really has nothing to do with polarizable force fields and can equally well be used in electronic structure calculations (and it often is).


So, the question now becomes: how do we compute the interaction of a molecule with an applied field in polarizable force field.

I'll address induced dipole force fields and fluctuating charge force fields separately, since the answer is slightly different in each case.

Starting with induced dipole force fields, the working equations are as follows:

$$\tag{5} \begin{bmatrix} \alpha_1^{-1} & -T_{12}^{dd} & \dots & -T_{1N}^{dd} \\ -T_{21}^{dd} & \ddots & \dots & \vdots \\ \vdots & & & \\ -T_{N1}^{dd} & \dots & & \alpha_N^{-1} \\ \end{bmatrix} \begin{bmatrix} \mathbf{P}_1 \\ \mathbf{P}_2 \\ \vdots \\ \mathbf{P}_N \end{bmatrix} = \begin{bmatrix} \mathbf{E}_1^q \\ \mathbf{E}_2^q \\ \vdots \\ \mathbf{E}_N^q \end{bmatrix} $$

$T_{12}^{dd}$ is the dipole-dipole interaction tensor and everything else has the same meaning. For this model, it is rather simple to compute the interaction with the field. We just fill our solution vector with the field applied in each direction and use the finite difference scheme described above.

Something I've never seen discussed in the literature is that you can actually compute the in-molecule atomic polarizabillities in this way and they well generally not equal the atomic polarizabilities used as parameters. This is because the interactions between induced dipoles within a molecule are strongly damped and hence can result in slight changes in the atomic polarizability (after accounting for changes due to orientation).

If you are only interested in the molecular polarzability, then you just do the finite difference scheme on the induced dipole of the molecule which you get by summing the atomic induced dipoles, $$ \mathbf{P}^{ind}_{mol}=\sum_{i=1}^{N_{atom}}\mathbf{P}^{ind}_i \tag{6} $$

If you were using a fluctuating charge model, rather than a dipole-polarization model, the procedure is slightly different.

We need to recognize that the molecular dipole can be computed from the cartesian coordinates and charges as, $$ \mathbf{P}_{mol}=\sum_{i=1}^{N_{atom}}q_i\mathbf{R}_i\tag{7} $$

When charges are allowed to rearrange within a molecule, we get a changes to the charge at each atom, $\delta q$. Then, the total dipole (permanent plus induced dipole) is, $$ \mathbf{P}_{mol}=\sum_{i=1}^{N_{atom}}(q_i+\delta q_i)\mathbf{R}_i\tag{8} $$

The permanent dipole doesn't contribute to the polarizability which is clear from the fact the $q_i$ terms will go to zero when taking the derivative with respect to the field.

I won't write it out in detail here, but in the fluctuating charge polarizabilty model, you end up adding a term to your solution vector which is $\mathbf{R}_i\cdot\mathbf{E}_w$ since multiplying this by $\delta q_i$ gives the interaction of the induced charge with the field.


Hopefully this clears things up. I can add more information if needed. I couldn't find any good references for how to do this in the induced dipole case, so I don't have much to cite but I am confident in the procedure since it's the same as what is done in electronic structure. I found one reference for the fluctuating charge case which I can't seem to locate anymore. I don't know why this hasn't been discussed more in the literature...

Maybe because if you compute the polarizability hyper-surface and compare it with what you get from electronic structure, the comparison is not very favorable even for the best polarizable force fields. Maybe this is a research question that someone (I) could address in the future (present)... :)

$\endgroup$
3
  • $\begingroup$ +10. I've just added equation labels, because even if you're not referring to any specific equations in your answer, someone else might want to say "I used Eq. 2 of this post" instead of saying "I used the second equation in this post". Also if you agree with my edit for the superscript in what's now Eq. 2, please go ahead and make that consistent in the rest of the answer! $\endgroup$ Jan 5, 2023 at 18:32
  • 1
    $\begingroup$ Thanks. Let me look through this and the Applequist paper again. As I mentioned, I basically want to estimate the full anisotropic molecular polarizability tensor from the isotropic atomic polarizabilities. We have a lot of QM-derived polarizabilities, so if you're interested, please drop me an e-mail. $\endgroup$ Jan 5, 2023 at 18:35
  • 1
    $\begingroup$ @GeoffHutchison Right. It shouldn't matter if the atomic polarizabilities are isotropic or anisotropic. Just doing the finite difference of the induced dipole should work. I haven't tried this with anything but relatively small systems though. The anisotropy is induced by the dipole-dipole interactions I'm pretty sure. $\endgroup$
    – jheindel
    Jan 5, 2023 at 22:36

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .