pair_style lj/BGcleavs3¶

Syntax¶

pair_style lj/BGcleavs3 cutoff1 cutoff2 lambda Dfac

cutoff1 = global internal cut-off
cutoff2 = global external cut-off
lambda = global scaling of the potential
Dfac = multiplicative factor for the derivative of the interaction with respect to lambda

pair_coeff accepts the same arguments as the pair_style. However, keep in mind that parameters in pair coeff have a specific order:

pair_coeff a b lambda Dfac cutoff1 cutoff2

where

a       = atom of type a [mandatory]
b       = atom of type b [mandatory]
lambda  = scaling of the potential
Dfac    = multiplicative factor for the derivative of the interaction with respect to lambda     
cutoff1 = global internal cut-off
cutoff2 = global external cut-off

The arguments in the pair_coeff are all optional except for the atom types. However, if we need to specify a particular argument, we need to specify also all the preceeding ones. If we want to specify a new lambda for a different pair of types we can just write

pair_coeff 1 1 1.0

Instead, in order to specify a new (e.g., external) cutoff we need to rewrite all the preceeding arguments:

pair_coeff 1 1 lambda Dfac cutoff1 cutoff2

Examples¶

pair_style lj/BGcleavs3 2.3 2.5 
pair_coeff 1 1 1.0
pair_coeff 1 2 1.0
pair_coeff 1 3 -1.0

Description¶

This pair style implements the Broughton and Gilmer1 modification to Lennard-Jones potential (see pair lj/BG) to be used in the step3 of the cleaving algorithm. This pair style returns also an array with all the calculated interactions which are needed to calculate the work in the step3. This array can be accessed using the new compute compute paircleav.

The potential implemented in this pair style is

\[\begin{split} U(r_{ln},\lambda) = \begin{cases} \lambda\left[ 4\epsilon \left(\left(\frac{\sigma}{r_{ln}}\right)^{12} -\left(\frac{\sigma}{r_{ln}}\right)^{6} \right)+C_1 \right],\;\mbox{if}\; r_{ln} \leq 2.3\sigma \\ \lambda\left[ C_2\left(\frac{\sigma}{r_{ln}}\right)^{12} + C_3\left(\frac{\sigma}{r_{ln}}\right)^{6} + C_4\left(\frac{r_{ln}}{\sigma}\right)^2 + C_5\right],\;\mbox{if}\; r_{ln} \leq 2.5\sigma \\ 0, \; r_{ln} \geq 2.5\sigma \end{cases} \end{split}\]

where $r_{ln}=|\mathbf{r}_l-\mathbf{r}_n|$ for each couple of atoms $l,n$ in the system, and $C_1, C_2, C_3, C_4, C_5$ are constants hardcorded within the pair style, therefore they not need to be defined. For C_1, C_2, C_3, C_4, C_5$ we used the values reported in Davidchack and Laird2.

Note

$\lambda$ must vary in the (closed) interval [0,1]
There is no internal check that $\lambda$ is changing consistently (i.e., always decreasing or increasing)

The work performed on the system to switch from $\lambda=1$ to $\lambda=0$ is then given by

\[\begin{split} \frac{\partial U(r_{ln},\lambda)}{\partial \lambda} = \begin{cases} \left[ 4\epsilon \left(\left(\frac{\sigma}{r_{ln}}\right)^{12} -\left(\frac{\sigma}{r_{ln}}\right)^{6} \right)+C_1 \right],\;\mbox{if}\; r_{ln} \leq 2.3\sigma \\ \left[ C_2\left(\frac{\sigma}{r_{ln}}\right)^{12} + C_3\left(\frac{\sigma}{r_{ln}}\right)^{6} + C_4\left(\frac{r_{ln}}{\sigma}\right)^2 + C_5\right],\;\mbox{if}\; r_{ln} \leq 2.5\sigma \\ 0, \; r_{ln} \geq 2.5\sigma \end{cases} \end{split}\]

In this case, the the expression of the potential is multiplied by the derivative of $\lambda$ with respect itself, which is equal to 1 (i.e., $\partial \lambda / \partial \lambda=1$). Therefore, the value of Dfrac must be set equal to 1.

Using this style we can represent different switches between different states of the system:

We want to describe the switching between two different states of the system, e.g., $a$ and $a+b$. Let us assume the total Hamiltonian of the system, $H^{T}(\lambda)$, is defined as

\[ H^{T}(\lambda)= \lambda H_{a} + H_{b} \]

for $\lambda \in [0,1]$, where we dropped the dependence on $\mathbf{p}$ and $\mathbf{q}$ for simplicity.

Let’s consider this example:

variable lambda file lambda.dat

pair_style lj/BGNlcleavs3  2.3 2.5 1.0 1.0
pair_coeff 1 1 1.0 1.0
pair_coeff 2 2 1.0 1.0

pair_coeff 1 2 ${lambda} 1.0

The “self” interactions between atoms of the same type remains, but the mixed interactions between atoms of different types are switched-on (if $\lambda$ goes from 0 to 1) or switched-off (if $\lambda$ goes from 1 to 0).

Note

If the switching required is from state $a$ to state $a+b$, then we can use this pair style with a polynomial function of $\lambda$. In this case we need to pass $\lambda^N$ as an argument to the pair style, e.g., with $N=3$:

variable lam file lambda.dat
variable lambda equal ${lam}*${lam}*${lam}

pair_style lj/BGlcleavs3  2.3 2.5 1.0 1.0
pair_coeff 1 1 1.0 1.0
pair_coeff 2 2 1.0 1.0

pair_coeff 1 2 ${lambda} 1.0

However, the work calculated by using the compute paircleav needs to be adjusted by the correct derivative with respect to $\lambda$. E.g., in the case just showed we need to multiply the work by $N \lambda^{N-1}$. For this reason, even if it possible in principle to use this pair style with a polynomial function of $\lambda$ with $N>1$, we suggest for $N>1$ to use the pair_style lj/BGNlcleavs3.

We want to describe the simultaneous switch between two different states of the system, e.g., $a$ and $b$. Let us assume the total Hamiltonian of the system, $H^{T}(\lambda)$, is defined as

\[ H^{T}(\lambda)= \lambda H_{a} + (1-\lambda) H_{b} \]

for $\lambda \in [0,1]$, where we dropped the dependence on $\mathbf{p}$ and $\mathbf{q}$ for simplicity. In this latter case, the command Dfrac needs to be used when we are passing to the pair_style an expression like $(1-\lambda)$. Here, we need to set Dfrac=-1 for some pair coefficients.

variable lambda file lambda.dat
variable minl   equal 1-lambda

pair_style lj/BGcleavs3  2.3 2.5 1.0 1.0
pair_coeff 1 1 1.0 1.0
pair_coeff 2 2 1.0 1.0
pair_coeff 3 3 1.0 1.0

pair_coeff 1 2 ${lambda} 1.0
pair_coeff 1 3 ${minl} -1.0
pair_coeff 2 3 1.0 1.0

The “self” interactions are not modified during the run, as well as the interactions between atoms of type 2 and 3. However, interactions between atoms of type 1 and 2 are switched-on (switched-off) and interactions between atoms of type 1 and 3 are switched-off (switched-on) if $\lambda$ goes from 0 to 1 (from 1 to 0).

Warning

This pair style can be used to simulate a polynomial switching from state $a$ to state $b$ in the form given by:

\[ H^{T}(\lambda)= \lambda^N H_{a} + (1-\lambda^N) H_{b} \]

by defining the polynomial in $\lambda$ outside the pair style, e.g.:

variable lam file lambda.dat
variable N equal 4
variable lambda equal exp($N*log(${lam}))
variable minl   equal 1-${lambda}

pair_style lj/BGcleavs3  2.3 2.5 1.0 1.0
pair_coeff 1 1 1.0 1.0
pair_coeff 2 2 1.0 1.0
pair_coeff 3 3 1.0 1.0

pair_coeff 1 2 ${lambda} 1.0
pair_coeff 1 3 ${minl} -1.0
pair_coeff 2 3 1.0 1.0

However, the energy obtained through the compute paircleav must be multiplied by $N\lambda^{N-1}$ as the pair style does not recognize that $\lambda$ is not a linear parameter. Note that, if a scaling of the form

\[ H^{T}(\lambda)= \lambda^N H_{a} + (1-\lambda)^N H_{b} \]

is needed, while in principle can be also defined with this pair style, we suggest to use the pair_style lj/BGNlcleavs3 instead. The reason is that, with such a scaling the coefficients for $H_{a}$ interactions is going to be $N\lambda^{N-1}$, while the coefficients for $H_{b}$ interactions is given by $-N(1-\lambda)^{N-1}$. The pair_style lj/BGNlcleavs3 is able to automatically account for this.

1: J. Q. Broughton and G. H. Gilmer. Molecular dynamics investigation of the crystal–fluid interface. i. bulk properties. Journal of Chemical Physics, 79(10):5095–5104, 1983.
2: R. L. Davidchack and B. B. Laird. Direct calculation of the crystal–melt interfacial free energies for continuous potentials: application to the lennard-jones system. The Journal of chemical physics, 118(16):7651–7657, 2003.