Shell molecular dynamics¶

GROMACS can simulate polarizability using the shell model of Dick and Overhauser 43. In such models a shell particle representing the electronic degrees of freedom is attached to a nucleus by a spring. The potential energy is minimized with respect to the shell position at every step of the simulation (see below). Successful applications of shell models in GROMACS have been published for $$N_2$$ 44 and water45.

Optimization of the shell positions¶

The force $$\mathbf{F}$$$$_S$$ on a shell particle $$S$$ can be decomposed into two components

(1)$\mathbf{F}_S ~=~ \mathbf{F}_{bond} + \mathbf{F}_{nb}$

where $$\mathbf{F}_{bond}$$ denotes the component representing the polarization energy, usually represented by a harmonic potential and $$\mathbf{F}_{nb}$$ is the sum of Coulomb and van der Waals interactions. If we assume that $$\mathbf{F}_{nb}$$ is almost constant we can analytically derive the optimal position of the shell, i.e. where $$\mathbf{F}_S$$ = 0. If we have the shell S connected to atom A we have

(2)$\mathbf{F}_{bond} ~=~ k_b \left( \mathbf{x}_S - \mathbf{x}_A\right).$

In an iterative solver, we have positions $$\mathbf{x}_S(n)$$ where $$n$$ is the iteration count. We now have at iteration $$n$$

(3)$\mathbf{F}_{nb} ~=~ \mathbf{F}_S - k_b \left( \mathbf{x}_S(n) - \mathbf{x}_A\right)$

and the optimal position for the shells $$x_S(n+1)$$ thus follows from

(4)$\mathbf{F}_S - k_b \left( \mathbf{x}_S(n) - \mathbf{x}_A\right) + k_b \left( \mathbf{x}_S(n+1) - \mathbf{x}_A\right) = 0$

if we write

(5)$\Delta \mathbf{x}_S = \mathbf{x}_S(n+1) - \mathbf{x}_S(n)$

we finally obtain

(6)$\Delta \mathbf{x}_S = \mathbf{F}_S/k_b$

which then yields the algorithm to compute the next trial in the optimization of shell positions

(7)$\mathbf{x}_S(n+1) ~=~ \mathbf{x}_S(n) + \mathbf{F}_S/k_b.$