 Open Source Biology & Genetics Interest Group

Open source scripts, reports, and preprints for in vitro biology, genetics, bioinformatics, crispr, and other biotech applications.

# Shell molecular dynamics — GROMACS 2020.4 documentation

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}_Aright).]

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}_Aright)]

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}_Aright) + k_b left( mathbf{x}_S(n+1) – mathbf{x}_Aright) = 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.] 