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 N2 44 and water45.

Optimization of the shell positions#

The force FS on a shell particle S can be decomposed into two components

(93)#FS = Fbond+Fnb

where Fbond denotes the component representing the polarization energy, usually represented by a harmonic potential and Fnb is the sum of Coulomb and van der Waals interactions. If we assume that Fnb is almost constant we can analytically derive the optimal position of the shell, i.e. where FS=0. If we have the shell S connected to atom A we have

(94)#Fbond = kb(xSxA).

In an iterative solver, we have positions xS(n) where n is the iteration count. We now have at iteration n

(95)#Fnb = FSkb(xS(n)xA)

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

(96)#FSkb(xS(n)xA)+kb(xS(n+1)xA)=0

if we write

(97)#ΔxS=xS(n+1)xS(n)

we finally obtain

(98)#ΔxS=FS/kb

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

(99)#xS(n+1) = xS(n)+FS/kb.