ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Added Vashishta potential

Browse Source 
git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@14117 f3b2605a-c512-4ea7-a41b-209d697bcdaa
This commit is contained in:
athomps
2015-10-15 17:39:43 +00:00
parent edf114fd25
commit c7170a4296
19 changed files with 1564 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
doc/Section_commands.txt
View File

@ -887,6 +887,7 @@ KOKKOS, o = USER-OMP, t = OPT.
"tip4p/cut (o)"_pair_coul.html,
"tip4p/long (o)"_pair_coul.html,
"tri/lj (o)"_pair_tri_lj.html,
"vashishta"_pair_vashishta.html,
"yukawa (go)"_pair_yukawa.html,
"yukawa/colloid (go)"_pair_yukawa_colloid.html,
"zbl (go)"_pair_zbl.html :tb(c=4,ea=c)

1
doc/Section_example.txt
View File

@ -80,6 +80,7 @@ snap: NVE dynamics for BCC tantalum crystal using SNAP potential
srd: stochastic rotation dynamics (SRD) particles as solvent
tad: temperature-accelerated dynamics of vacancy diffusion in bulk Si
tri: triangular particles in rigid bodies :tb(s=:)
vashishta: models using the Vashishta potential
Here is how you might run and visualize one of the sample problems:

1
doc/pair_style.txt
View File

@ -200,6 +200,7 @@ in the pair section of "this page"_Section_commands.html#cmd_5.
"pair_style tip4p/cut"_pair_coul.html - Coulomb for TIP4P water w/out LJ
"pair_style tip4p/long"_pair_coul.html - long-range Coulombics for TIP4P water w/out LJ
"pair_style tri/lj"_pair_tri_lj.html - LJ potential between triangles
"pair_style vashishta"_pair_vahishta.html - Vashishta 2-body and 3-body potential
"pair_style yukawa"_pair_yukawa.html - Yukawa potential
"pair_style yukawa/colloid"_pair_yukawa_colloid.html - screened Yukawa potential for finite-size particles
"pair_style zbl"_pair_zbl.html - Ziegler-Biersack-Littmark potential :ul

211
doc/pair_vashishta.html Normal file
View File

@ -0,0 +1,211 @@
<HTML>
<CENTER><A HREF = "http://lammps.sandia.gov">LAMMPS WWW Site</A> - <A HREF = "Manual.html">LAMMPS Documentation</A> - <A HREF = "Section_commands.html#comm">LAMMPS Commands</A>
</CENTER>
<HR>
<H3>pair_style vashishta command
</H3>
<P><B>Syntax:</B>
</P>
<PRE>pair_style vashishta
</PRE>
<P><B>Examples:</B>
</P>
<PRE>pair_style vashishta
pair_coeff * * SiC.vashishta Si C
</PRE>
<P><B>Description:</B>
</P>
<P>The <I>vashishta</I> style computes the combined 2-body and 3-body
family of potentials developed in the group of Vashishta and
co-workers. By combining repulsive, screened Coulombic,
screened charge-dipole, and dispersion interactions with a
bond-angle energy based on the Stillinger-Weber potential,
this potential has been used to describe a variety of inorganic
compounds, including SiO2 <A HREF = "#Vashishta1990">Vashishta1990</A>,
SiC <A HREF = "#Vashishta2007">Vashishta2007</A>,
and InP <A HREF = "#Branicio2009">Branicio2009</A>.
</P>
<P>The potential for the energy U of a system of atoms is
</P>
<CENTER><IMG SRC = "Eqs/pair_vashishta.jpg">
</CENTER>
<P>where we follow the notation used in <A HREF = "#Branicio2009">Branicio2009</A>.
U2 is a two-body term and U3 is a three-body term. The
summation over two-body terms is over all neighbors J within
a cutoff distance = <I>rc</I>. The twobody terms are shifted and
tilted by a linear function so that the energy and force are
both zero at <I>rc</I>. The summation over three-body terms
is over all neighbors J and K within a cut-off distance = <I>r0</I>,
where the exponential screening function becomes zero.
</P>
<P>Only a single pair_coeff command is used with the <I>vashishta</I> style which
specifies a Vashishta potential file with parameters for all
needed elements. These are mapped to LAMMPS atom types by specifying
N additional arguments after the filename in the pair_coeff command,
where N is the number of LAMMPS atom types:
</P>
<UL><LI>filename
<LI>N element names = mapping of Vashishta elements to atom types
</UL>
<P>See the <A HREF = "pair_coeff.html">pair_coeff</A> doc page for alternate ways
to specify the path for the potential file.
</P>
<P>As an example, imagine a file SiC.vashishta has parameters for
Si and C. If your LAMMPS simulation has 4 atoms types and you want
the 1st 3 to be Si, and the 4th to be C, you would use the following
pair_coeff command:
</P>
<PRE>pair_coeff * * SiC.vashishta Si Si Si C
</PRE>
<P>The 1st 2 arguments must be * * so as to span all LAMMPS atom types.
The first three Si arguments map LAMMPS atom types 1,2,3 to the Si
element in the file. The final C argument maps LAMMPS atom type 4
to the C element in the file. If a mapping value is specified as
NULL, the mapping is not performed. This can be used when a <I>vashishta</I>
potential is used as part of the <I>hybrid</I> pair style. The NULL values
are placeholders for atom types that will be used with other
potentials.
</P>
<P>Vashishta files in the <I>potentials</I> directory of the LAMMPS
distribution have a ".vashishta" suffix. Lines that are not blank or
comments (starting with #) define parameters for a triplet of
elements. The parameters in a single entry correspond to the two-body
and three-body coefficients in the formulae above:
</P>
<UL><LI>element 1 (the center atom in a 3-body interaction)
<LI>element 2
<LI>element 3
<LI>H (energy units)
<LI>eta
<LI>Zi (electron charge units)
<LI>Zj (electron charge units)
<LI>lambda1 (distance units)
<LI>D (energy units)
<LI>lambda4 (distance units)
<LI>W (energy units)
<LI>rc (distance units)
<LI>B (energy units)
<LI>gamma
<LI>r0 (distance units)
<LI>C
<LI>costheta0
</UL>
<P>The non-annotated parameters are unitless.
The Vashishta potential file must contain entries for all the
elements listed in the pair_coeff command. It can also contain
entries for additional elements not being used in a particular
simulation; LAMMPS ignores those entries.
For a single-element simulation, only a single entry is required
(e.g. SiSiSi). For a two-element simulation, the file must contain 8
entries (for SiSiSi, SiSiC, SiCSi, SiCC, CSiSi, CSiC, CCSi, CCC), that
specify parameters for all permutations of the two elements
interacting in three-body configurations. Thus for 3 elements, 27
entries would be required, etc.
</P>
<P>Depending on the particular version of the Vashishta potential,
the values of these parameters may be keyed to the identities of
zero, one, two, or three elements.
In order to make the input file format unambiguous, general,
and simple to code,
LAMMPS uses a slightly confusing method for specifying parameters.
All parameters are divided into two classes: two-body and three-body.
Two-body and three-body parameters are handled differently,
as described below.
The two-body parameters are H, eta, lambda1, D, lambda4, W, rc, gamma, and r0.
They appear in the above formulae with two subscripts.
The parameters Zi and Zj are also classified as two-body parameters,
even though they only have 1 subscript.
The three-body parameters are B, C, costheta0.
They appear in the above formulae with three subscripts.
Two-body and three-body parameters are handled differently,
as described below.
</P>
<P>The first element in each entry is the center atom
in a three-body interaction, while the second and third elements
are two neighbor atoms. Three-body parameters for a central atom I
and two neighbors J and K are taken from the IJK entry.
Note that even though three-body parameters do not depend on the order of
J and K, LAMMPS stores three-body parameters for both IJK and IKJ.
The user must ensure that these values are equal.
Two-body parameters for an atom I interacting with atom J are taken from
the IJJ entry, where the 2nd and 3rd
elements are the same. Thus the two-body parameters
for Si interacting with C come from the SiCC entry. Note that even
though two-body parameters (except possibly gamma and r0 in U3)
do not depend on the order of the two elements,
LAMMPS will get the Si-C value from the SiCC entry
and the C-Si value from the CSiSi entry. The user must ensure
that these values are equal. Two-body parameters appearing
in entries where the 2nd and 3rd elements are different are
stored but never used. It is good practice to enter zero for
these values. Note that the three-body function U3 above
contains the two-body parameters gamma and r0. So U3 for a
central C atom bonded to an Si atom and a second C atom
will take three-body parameters from the CSiC entry, but
two-body parameters from the CCC and CSiSi entries.
</P>
<HR>
<P><B>Mixing, shift, table, tail correction, restart, rRESPA info</B>:
</P>
<P>For atom type pairs I,J and I != J, where types I and J correspond to
two different element types, mixing is performed by LAMMPS as
described above from values in the potential file.
</P>
<P>This pair style does not support the <A HREF = "pair_modify.html">pair_modify</A>
shift, table, and tail options.
</P>
<P>This pair style does not write its information to <A HREF = "restart.html">binary restart
files</A>, since it is stored in potential files. Thus, you
need to re-specify the pair_style and pair_coeff commands in an input
script that reads a restart file.
</P>
<P>This pair style can only be used via the <I>pair</I> keyword of the
<A HREF = "run_style.html">run_style respa</A> command. It does not support the
<I>inner</I>, <I>middle</I>, <I>outer</I> keywords.
</P>
<HR>
<P><B>Restrictions:</B>
</P>
<P>This pair style is part of the MANYBODY package. It is only enabled
if LAMMPS was built with that package (which it is by default). See
the <A HREF = "Section_start.html#start_3">Making LAMMPS</A> section for more info.
</P>
<P>This pair style requires the <A HREF = "newton.html">newton</A> setting to be "on"
for pair interactions.
</P>
<P>The Vashishta potential files provided with LAMMPS (see the
potentials directory) are parameterized for metal <A HREF = "units.html">units</A>.
You can use the Vashishta potential with any LAMMPS units, but you would need
to create your own Vashishta potential file with coefficients listed in the
appropriate units if your simulation doesn't use "metal" units.
</P>
<P><B>Related commands:</B>
</P>
<P><A HREF = "pair_coeff.html">pair_coeff</A>
</P>
<P><B>Default:</B> none
</P>
<HR>
<A NAME = "Vashishta1990"></A>
<P><B>(Vashishta1990)</B> P. Vashishta, R. K. Kalia, J. P. Rino, Phys. Rev. B 41, 12197 (1990).
</P>
<A NAME = "Vashishta2007"></A>
<P><B>(Vashishta2007)</B> P. Vashishta, R. K. Kalia, A. Nakano, J. P. Rino. J. Appl. Phys. 101, 103515 (2007).
</P>
<A NAME = "Branicio2009"></A>
<P><B>(Branicio2009)</B> Branicio, Rino, Gan and Tsuzuki, J. Phys Condensed Matter 21 (2009) 095002
</P>
</HTML>

204
doc/pair_vashishta.txt Normal file
View File

@ -0,0 +1,204 @@
"LAMMPS WWW Site"_lws - "LAMMPS Documentation"_ld - "LAMMPS Commands"_lc :c
:link(lws,http://lammps.sandia.gov)
:link(ld,Manual.html)
:link(lc,Section_commands.html#comm)
:line
pair_style vashishta command :h3
[Syntax:]
pair_style vashishta :pre
[Examples:]
pair_style vashishta
pair_coeff * * SiC.vashishta Si C :pre
[Description:]
The {vashishta} style computes the combined 2-body and 3-body
family of potentials developed in the group of Vashishta and
co-workers. By combining repulsive, screened Coulombic,
screened charge-dipole, and dispersion interactions with a
bond-angle energy based on the Stillinger-Weber potential,
this potential has been used to describe a variety of inorganic
compounds, including SiO2 "Vashishta1990"_#Vashishta1990,
SiC "Vashishta2007"_#Vashishta2007,
and InP "Branicio2009"_#Branicio2009.
The potential for the energy U of a system of atoms is
:c,image(Eqs/pair_vashishta.jpg)
where we follow the notation used in "Branicio2009"_#Branicio2009.
U2 is a two-body term and U3 is a three-body term. The
summation over two-body terms is over all neighbors J within
a cutoff distance = {rc}. The twobody terms are shifted and
tilted by a linear function so that the energy and force are
both zero at {rc}. The summation over three-body terms
is over all neighbors J and K within a cut-off distance = {r0},
where the exponential screening function becomes zero.
Only a single pair_coeff command is used with the {vashishta} style which
specifies a Vashishta potential file with parameters for all
needed elements. These are mapped to LAMMPS atom types by specifying
N additional arguments after the filename in the pair_coeff command,
where N is the number of LAMMPS atom types:
filename
N element names = mapping of Vashishta elements to atom types :ul
See the "pair_coeff"_pair_coeff.html doc page for alternate ways
to specify the path for the potential file.
As an example, imagine a file SiC.vashishta has parameters for
Si and C. If your LAMMPS simulation has 4 atoms types and you want
the 1st 3 to be Si, and the 4th to be C, you would use the following
pair_coeff command:
pair_coeff * * SiC.vashishta Si Si Si C :pre
The 1st 2 arguments must be * * so as to span all LAMMPS atom types.
The first three Si arguments map LAMMPS atom types 1,2,3 to the Si
element in the file. The final C argument maps LAMMPS atom type 4
to the C element in the file. If a mapping value is specified as
NULL, the mapping is not performed. This can be used when a {vashishta}
potential is used as part of the {hybrid} pair style. The NULL values
are placeholders for atom types that will be used with other
potentials.
Vashishta files in the {potentials} directory of the LAMMPS
distribution have a ".vashishta" suffix. Lines that are not blank or
comments (starting with #) define parameters for a triplet of
elements. The parameters in a single entry correspond to the two-body
and three-body coefficients in the formulae above:
element 1 (the center atom in a 3-body interaction)
element 2
element 3
H (energy units)
eta
Zi (electron charge units)
Zj (electron charge units)
lambda1 (distance units)
D (energy units)
lambda4 (distance units)
W (energy units)
rc (distance units)
B (energy units)
gamma
r0 (distance units)
C
costheta0 :ul
The non-annotated parameters are unitless.
The Vashishta potential file must contain entries for all the
elements listed in the pair_coeff command. It can also contain
entries for additional elements not being used in a particular
simulation; LAMMPS ignores those entries.
For a single-element simulation, only a single entry is required
(e.g. SiSiSi). For a two-element simulation, the file must contain 8
entries (for SiSiSi, SiSiC, SiCSi, SiCC, CSiSi, CSiC, CCSi, CCC), that
specify parameters for all permutations of the two elements
interacting in three-body configurations. Thus for 3 elements, 27
entries would be required, etc.
Depending on the particular version of the Vashishta potential,
the values of these parameters may be keyed to the identities of
zero, one, two, or three elements.
In order to make the input file format unambiguous, general,
and simple to code,
LAMMPS uses a slightly confusing method for specifying parameters.
All parameters are divided into two classes: two-body and three-body.
Two-body and three-body parameters are handled differently,
as described below.
The two-body parameters are H, eta, lambda1, D, lambda4, W, rc, gamma, and r0.
They appear in the above formulae with two subscripts.
The parameters Zi and Zj are also classified as two-body parameters,
even though they only have 1 subscript.
The three-body parameters are B, C, costheta0.
They appear in the above formulae with three subscripts.
Two-body and three-body parameters are handled differently,
as described below.
The first element in each entry is the center atom
in a three-body interaction, while the second and third elements
are two neighbor atoms. Three-body parameters for a central atom I
and two neighbors J and K are taken from the IJK entry.
Note that even though three-body parameters do not depend on the order of
J and K, LAMMPS stores three-body parameters for both IJK and IKJ.
The user must ensure that these values are equal.
Two-body parameters for an atom I interacting with atom J are taken from
the IJJ entry, where the 2nd and 3rd
elements are the same. Thus the two-body parameters
for Si interacting with C come from the SiCC entry. Note that even
though two-body parameters (except possibly gamma and r0 in U3)
do not depend on the order of the two elements,
LAMMPS will get the Si-C value from the SiCC entry
and the C-Si value from the CSiSi entry. The user must ensure
that these values are equal. Two-body parameters appearing
in entries where the 2nd and 3rd elements are different are
stored but never used. It is good practice to enter zero for
these values. Note that the three-body function U3 above
contains the two-body parameters gamma and r0. So U3 for a
central C atom bonded to an Si atom and a second C atom
will take three-body parameters from the CSiC entry, but
two-body parameters from the CCC and CSiSi entries.
:line
[Mixing, shift, table, tail correction, restart, rRESPA info]:
For atom type pairs I,J and I != J, where types I and J correspond to
two different element types, mixing is performed by LAMMPS as
described above from values in the potential file.
This pair style does not support the "pair_modify"_pair_modify.html
shift, table, and tail options.
This pair style does not write its information to "binary restart
files"_restart.html, since it is stored in potential files. Thus, you
need to re-specify the pair_style and pair_coeff commands in an input
script that reads a restart file.
This pair style can only be used via the {pair} keyword of the
"run_style respa"_run_style.html command. It does not support the
{inner}, {middle}, {outer} keywords.
:line
[Restrictions:]
This pair style is part of the MANYBODY package. It is only enabled
if LAMMPS was built with that package (which it is by default). See
the "Making LAMMPS"_Section_start.html#start_3 section for more info.
This pair style requires the "newton"_newton.html setting to be "on"
for pair interactions.
The Vashishta potential files provided with LAMMPS (see the
potentials directory) are parameterized for metal "units"_units.html.
You can use the Vashishta potential with any LAMMPS units, but you would need
to create your own Vashishta potential file with coefficients listed in the
appropriate units if your simulation doesn't use "metal" units.
[Related commands:]
"pair_coeff"_pair_coeff.html
[Default:] none
:line
:link(Vashishta1990)
[(Vashishta1990)] P. Vashishta, R. K. Kalia, J. P. Rino, Phys. Rev. B 41, 12197 (1990).
:link(Vashishta2007)
[(Vashishta2007)] P. Vashishta, R. K. Kalia, A. Nakano, J. P. Rino. J. Appl. Phys. 101, 103515 (2007).
:link(Branicio2009)
[(Branicio2009)] Branicio, Rino, Gan and Tsuzuki, J. Phys Condensed Matter 21 (2009) 095002

1
examples/README
View File

@ -98,6 +98,7 @@ srd: stochastic rotation dynamics (SRD) particles as solvent
snap: NVE dynamics for BCC tantalum crystal using SNAP potential
streitz: Streitz-Mintmire potential for Al2O3
tad: temperature-accelerated dynamics of vacancy diffusion in bulk Si
vashishta: models using the Vashishta potential
voronoi: test of Voronoi tesselation in compute voronoi/atom
Here is a src/Make.py command which will perform a parallel build of a

38
examples/vashishta/InP.vashishta Executable file
View File

@ -0,0 +1,38 @@
# DATE: 2015-10-14 CONTRIBUTOR: Aidan Thompson, athomps@sandia.gov CITATION: Branicio, Rino, Gan and Tsuzuki, J. Phys Condensed Matter 21 (2009) 095002
#
# Vashishta potential file for InP, Branicio, Rino, Gan and Tsuzuki,
# J. Phys Condensed Matter 21 (2009) 095002
#
# These entries are in LAMMPS "metal" units:
# H = eV*Angstroms^eta; Zi, Zj = |e| (e = electronic charge);
# lambda1, lambda4, rc, r0, gamma = Angstroms;
# D = eV*Angstroms^4; W = eV*Angstroms^6; B = eV;
# other quantities are unitless
# element1 element2 element3
# H eta Zi Zj lambda1 D lambda4
# W rc B gamma r0 C cos(theta)
In In In 273.584 7 -1.21 -1.21 4.5 0.0 2.75
0.0 6.0 0.0 0.0 0.0 0.0 0.0
P P P 1813.06 7 1.21 1.21 4.5 52.7067 2.75
0.0 6.0 0.0 0.0 0.0 0.0 -0.333333333333
In P P 4847.09 9 1.21 -1.21 4.5 26.3533 2.75
270.105 6.0 4.34967 1.0 3.55 7.0 -0.333333333333
P In In 4847.09 9 1.21 -1.21 4.5 26.3533 2.75
270.105 6.0 4.34967 1.0 3.55 7.0 -0.333333333333
In In P 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0
In P In 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0
P In P 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0
P P In 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0

41
examples/vashishta/SiO.1990.vashishta Executable file
View File

@ -0,0 +1,41 @@
# DATE: 2015-10-14 CONTRIBUTOR: Aidan Thompson, athomps@sandia.gov CITATION: P. Vashishta, R. K. Kalia, J. P. Rino, and I. Ebbsjo, Phys. Rev. B 41, 12197 (1990).
#
# Vashishta potential file for SiO2, P. Vashishta, R. K. Kalia, J. P. Rino, and I. Ebbsjo,
# Phys. Rev. B 41, 12197 (1990).
#
# These parameters, some inferred indirectly, give a good
# match to the energy-volume curve for alpha-quartz in Fig. 2 of the paper.
#
# These entries are in LAMMPS "metal" units:
# H = eV*Angstroms^eta; Zi, Zj = |e| (e = electronic charge);
# lambda1, lambda4, rc, r0, gamma = Angstroms;
# D = eV*Angstroms^4; W = eV*Angstroms^6; B = eV;
# other quantities are unitless
#
# element1 element2 element3
# H eta Zi Zj lambda1 D lambda4
# W rc B gamma r0 C cos(theta)
Si Si Si 0.82023 11 1.6 1.6 999 0.0 4.43
0.0 10.0 0.0 0.0 0.0 0.0 0.0
O O O 743.848 7 -0.8 -0.8 999 22.1179 4.43
0.0 10.0 0.0 0.0 0.0 0.0 0.0
O Si Si 163.859 9 -0.8 1.6 999 44.2357 4.43
0.0 10.0 20.146 1.0 2.60 0.0 -0.77714596
Si O O 163.859 9 1.6 -0.8 999 44.2357 4.43
0.0 10.0 5.0365 1.0 2.60 0.0 -0.333333333333
Si O Si 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0
Si Si O 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0
O Si O 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0
O O Si 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0

26
examples/vashishta/data.quartz Normal file
View File

@ -0,0 +1,26 @@
# SiO2 alpha quartz
9 atoms
2 atom types
0 4.913400 xlo xhi
0 4.255129 ylo yhi
0 5.405200 zlo zhi
-2.456700 0.0 0.0 xy xz yz
Masses
1 28.0855
2 15.9994
Atoms
1 1 2.308807 0.000000 3.603467
2 1 -1.154403 1.999485 1.801733
3 1 -1.154403 -1.999485 0.000000
4 2 1.375998 1.140800 4.245244
5 2 -1.675961 0.621249 7.848711
6 2 0.299963 -1.762049 6.046977
7 2 0.299963 1.762049 -4.245244
8 2 -1.675961 -0.621249 -0.64177
9 2 1.375998 -1.140800 -2.443511

75
examples/vashishta/in.indiumphosphide Normal file
View File

@ -0,0 +1,75 @@
# calculate the energy volume curve for InP zincblende
# define volume range and filename
variable ndelta equal 100
variable volatom_min equal 20.0
variable volatom_max equal 29.0
variable evsvolfile string evsvol.dat
# set up cell
units metal
boundary p p p
# setup loop variables for box volume
variable amin equal ${volatom_min}^(1/3)*2
variable delta equal (${volatom_max}-${volatom_min})/${ndelta}
variable scale equal (${delta}/v_volatom+1)^(1/3)
# set up 8 atom InP zincblende unit cell
lattice diamond ${amin}
region box prism &
0 1 &
0 1 &
0 1 &
0 0 0
create_box 2 box
create_atoms 1 box &
basis 5 2 &
basis 6 2 &
basis 7 2 &
basis 8 2
mass 1 114.76
mass 2 30.98
# choose potential
pair_style vashishta
pair_coeff * * InP.vashishta In P
# setup neighbor style
neighbor 1.0 nsq
neigh_modify once no every 1 delay 0 check yes
# setup output
thermo_style custom step temp pe press vol
thermo_modify norm no
variable volatom equal vol/atoms
variable eatom equal pe/atoms
print "# Volume [A^3/atom] Energy [eV/atom]" file ${evsvolfile}
# loop over range of volumes
label loop
variable i loop ${ndelta}
change_box all x scale ${scale} y scale ${scale} z scale ${scale} remap
# calculate energy
# no energy minimization needed for zincblende
run 0
print "${volatom} ${eatom}" append ${evsvolfile}
next i
jump SELF loop

28
examples/vashishta/in.sio2 Normal file
View File

@ -0,0 +1,28 @@
# test Vashishta potential for quartz
units metal
boundary p p p
atom_style atomic
read_data data.quartz
replicate 4 4 4
velocity all create 2000.0 277387 mom yes
displace_atoms all move 0.05 0.9 0.4 units box
pair_style vashishta
pair_coeff * * SiO.1990.vashishta Si O
neighbor 0.3 bin
neigh_modify delay 10
fix 1 all nve
thermo 10
timestep 0.001
#dump 1 all cfg 10 *.cfg mass type xs ys zs vx vy vz fx fy fz
#dump_modify 1 element Si O
run 100

38
potentials/InP.vashishta Executable file
View File

@ -0,0 +1,38 @@
# DATE: 2015-10-14 CONTRIBUTOR: Aidan Thompson, athomps@sandia.gov CITATION: Branicio, Rino, Gan and Tsuzuki, J. Phys Condensed Matter 21 (2009) 095002
#
# Vashishta potential file for InP, Branicio, Rino, Gan and Tsuzuki,
# J. Phys Condensed Matter 21 (2009) 095002
#
# These entries are in LAMMPS "metal" units:
# H = eV*Angstroms^eta; Zi, Zj = |e| (e = electronic charge);
# lambda1, lambda4, rc, r0, gamma = Angstroms;
# D = eV*Angstroms^4; W = eV*Angstroms^6; B = eV;
# other quantities are unitless
# element1 element2 element3
# H eta Zi Zj lambda1 D lambda4
# W rc B gamma r0 C cos(theta)
In In In 273.584 7 -1.21 -1.21 4.5 0.0 2.75
0.0 6.0 0.0 0.0 0.0 0.0 0.0
P P P 1813.06 7 1.21 1.21 4.5 52.7067 2.75
0.0 6.0 0.0 0.0 0.0 0.0 -0.333333333333
In P P 4847.09 9 1.21 -1.21 4.5 26.3533 2.75
270.105 6.0 4.34967 1.0 3.55 7.0 -0.333333333333
P In In 4847.09 9 1.21 -1.21 4.5 26.3533 2.75
270.105 6.0 4.34967 1.0 3.55 7.0 -0.333333333333
In In P 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0
In P In 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0
P In P 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0
P P In 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0

1
potentials/README
View File

@ -100,3 +100,4 @@ sw Stillinger-Weber potential
tersoff Tersoff potential
tersoff.mod modified Tersoff potential
tersoff.zbl Tersoff with ZBL core
vashishta Vashishta 2-body and 3-body potential

41
potentials/SiC.vashishta Executable file
View File

@ -0,0 +1,41 @@
# DATE: 2015-10-14 CONTRIBUTOR: Aidan Thompson, athomps@sandia.gov CITATION: P. Vashishta, R. K. Kalia, A. Nakano, and J. P. Rino. J. Appl. Phys. 101, 103515 (2007).
#
# Vashishta potential file for SiC, P. Vashishta, R. K. Kalia, A. Nakano,
# and J. P. Rino. J. Appl. Phys. 101, 103515 (2007).
#
# These entries are in LAMMPS "metal" units:
# H = eV*Angstroms^eta; Zi, Zj = |e| (e = electronic charge);
# lambda1, lambda4, rc, r0, gamma = Angstroms;
# D = eV*Angstroms^4; W = eV*Angstroms^6; B = eV;
# other quantities are unitless
#
# Note: Value of D here equals D/2 in paper
#
# element1 element2 element3
# H eta Zi Zj lambda1 D lambda4
# W rc B gamma r0 C cos(theta)
C C C 471.74538 7 -1.201 -1.201 5.0 0.0 3.0
0.0 7.35 0.0 0.0 0.0 0.0 0.0
Si Si Si 23.67291 7 1.201 1.201 5.0 15.575 3.0
0.0 7.35 0.0 0.0 0.0 0.0 0.0
C Si Si 447.09026 9 -1.201 1.201 5.0 7.7874 3.0
61.4694 7.35 9.003 1.0 2.90 5.0 -0.333333333333
Si C C 447.09026 9 1.201 -1.201 5.0 7.7874 3.0
61.4694 7.35 9.003 1.0 2.90 5.0 -0.333333333333
C C Si 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0
C Si C 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0
Si C Si 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0
Si Si C 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0

41
potentials/SiO.1990.vashishta Executable file
View File

@ -0,0 +1,41 @@
# DATE: 2015-10-14 CONTRIBUTOR: Aidan Thompson, athomps@sandia.gov CITATION: P. Vashishta, R. K. Kalia, J. P. Rino, and I. Ebbsjo, Phys. Rev. B 41, 12197 (1990).
#
# Vashishta potential file for SiO2, P. Vashishta, R. K. Kalia, J. P. Rino, and I. Ebbsjo,
# Phys. Rev. B 41, 12197 (1990).
#
# These parameters, some inferred indirectly, give a good
# match to the energy-volume curve for alpha-quartz in Fig. 2 of the paper.
#
# These entries are in LAMMPS "metal" units:
# H = eV*Angstroms^eta; Zi, Zj = |e| (e = electronic charge);
# lambda1, lambda4, rc, r0, gamma = Angstroms;
# D = eV*Angstroms^4; W = eV*Angstroms^6; B = eV;
# other quantities are unitless
#
# element1 element2 element3
# H eta Zi Zj lambda1 D lambda4
# W rc B gamma r0 C cos(theta)
Si Si Si 0.82023 11 1.6 1.6 999 0.0 4.43
0.0 10.0 0.0 0.0 0.0 0.0 0.0
O O O 743.848 7 -0.8 -0.8 999 22.1179 4.43
0.0 10.0 0.0 0.0 0.0 0.0 0.0
O Si Si 163.859 9 -0.8 1.6 999 44.2357 4.43
0.0 10.0 20.146 1.0 2.60 0.0 -0.77714596
Si O O 163.859 9 1.6 -0.8 999 44.2357 4.43
0.0 10.0 5.0365 1.0 2.60 0.0 -0.333333333333
Si O Si 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0
Si Si O 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0
O Si O 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0
O O Si 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0

38
potentials/SiO.1994.vashishta Executable file
View File

@ -0,0 +1,38 @@
# DATE: 2015-10-14 CONTRIBUTOR: Aidan Thompson, athomps@sandia.gov CITATION: Nakano, Kalia, and Vashishta, J. Non-Crystal. Solids, v. 171, p. 157 (1994)
#
# Vashishta potential file for SiO2, Nakano, Kalia, and Vashishta,
# J. Non-Crystal. Solids, v. 171, p. 157 (1994)
#
# These entries are in LAMMPS "metal" units:
# H = eV*Angstroms^eta; Zi, Zj = |e| (e = electronic charge);
# lambda1, lambda4, rc, r0, gamma = Angstroms;
# D = eV*Angstroms^4; W = eV*Angstroms^6; B = eV;
# other quantities are unitless
#
# element1 element2 element3
# H eta Zi Zj lambda1 D lambda4
# W rc B gamma r0 C cos(theta)
Si Si Si 0.80603 11 1.76 1.76 4.43 0.0 2.5
0.0 5.5 0.0 0.0 0.0 0.0 0.0
O O O 730.17 7 -0.88 -0.88 4.43 26.7447 2.5
0.0 5.5 0.0 0.0 0.0 0.0 0.0
O Si Si 160.849 9 -0.88 1.76 4.43 53.4894 2.5
0.0 5.5 19.972 1.0 2.60 0.0 -0.77714596
Si O O 160.849 9 1.76 -0.88 4.43 53.4894 2.5
0.0 5.5 4.993 1.0 2.60 0.0 -0.333333333333
Si O Si 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0
Si Si O 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0
O Si O 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0
O O Si 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0

37
potentials/SiO.1997.vashishta Executable file
View File

@ -0,0 +1,37 @@
# DATE: 2015-10-14 CONTRIBUTOR: Aidan Thompson, athomps@sandia.gov CITATION: Broughton, Meli,Vashishta,and Kalia, Phys Rev B, v. 56, p. 611 (1997)
#
# Vashishta potential file for SiO2, Broughton, Meli,Vashishta,and Kalia, Phys Rev B, v. 56, p. 611 (1997)
#
# These entries are in LAMMPS "metal" units:
# H = eV*Angstroms^eta; Zi, Zj = |e| (e = electronic charge);
# lambda1, lambda4, rc, r0, gamma = Angstroms;
# D = eV*Angstroms^4; W = eV*Angstroms^6; B = eV;
# other quantities are unitless
#
# element1 element2 element3
# H eta Zi Zj lambda1 D lambda4
# W rc B gamma r0 C cos(theta)
Si Si Si 0.15500 11 0.7872 0.7872 4.43 0.0 2.5
0.0 5.5 0.0 0.0 0.0 0.0 0.0
O O O 140.38 7 -0.3936 -0.3936 4.43 5.3504 2.5
0.0 5.5 0.0 0.0 0.0 0.0 0.0
O Si Si 30.923 9 -0.3936 0.7872 4.43 10.701 2.5
0.0 5.5 19.972 1.0 2.60 0.0 -0.80593
Si O O 30.923 9 0.7872 -0.3936 4.43 10.701 2.5
0.0 5.5 4.993 1.0 2.60 0.0 -0.333333333333
Si O Si 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0
Si Si O 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0
O Si O 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0
O O Si 0.0 0.0 0.0 0.0 0.0 0.0 0.0
0.0 0.0 0.0 0.0 0.0 0.0 0.0

619
src/MANYBODY/pair_vashishta.cpp Executable file
View File

@ -0,0 +1,619 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing author: Aidan Thompson (SNL)
------------------------------------------------------------------------- */
#include "math.h"
#include "stdio.h"
#include "stdlib.h"
#include "string.h"
#include "pair_vashishta.h"
#include "atom.h"
#include "neighbor.h"
#include "neigh_request.h"
#include "force.h"
#include "comm.h"
#include "memory.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "memory.h"
#include "error.h"
using namespace LAMMPS_NS;
//using namespace MathConst;
#define MAXLINE 1024
#define DELTA 4
/* ---------------------------------------------------------------------- */
PairVashishta::PairVashishta(LAMMPS *lmp) : Pair(lmp)
{
single_enable = 0;
restartinfo = 0;
one_coeff = 1;
manybody_flag = 1;
nelements = 0;
elements = NULL;
nparams = maxparam = 0;
params = NULL;
elem2param = NULL;
}
/* ----------------------------------------------------------------------
check if allocated, since class can be destructed when incomplete
------------------------------------------------------------------------- */
PairVashishta::~PairVashishta()
{
if (elements)
for (int i = 0; i < nelements; i++) delete [] elements[i];
delete [] elements;
memory->destroy(params);
memory->destroy(elem2param);
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
delete [] map;
}
}
/* ---------------------------------------------------------------------- */
void PairVashishta::compute(int eflag, int vflag)
{
int i,j,k,ii,jj,kk,inum,jnum,jnumm1;
int itype,jtype,ktype,ijparam,ikparam,ijkparam;
tagint itag,jtag;
double xtmp,ytmp,ztmp,delx,dely,delz,evdwl,fpair;
double rsq,rsq1,rsq2;
double delr1[3],delr2[3],fj[3],fk[3];
int *ilist,*jlist,*numneigh,**firstneigh;
evdwl = 0.0;
if (eflag || vflag) ev_setup(eflag,vflag);
else evflag = vflag_fdotr = 0;
double **x = atom->x;
double **f = atom->f;
tagint *tag = atom->tag;
int *type = atom->type;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// loop over full neighbor list of my atoms
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
itag = tag[i];
itype = map[type[i]];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
// two-body interactions, skip half of them
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtag = tag[j];
if (itag > jtag) {
if ((itag+jtag) % 2 == 0) continue;
} else if (itag < jtag) {
if ((itag+jtag) % 2 == 1) continue;
} else {
if (x[j][2] < ztmp) continue;
if (x[j][2] == ztmp && x[j][1] < ytmp) continue;
if (x[j][2] == ztmp && x[j][1] == ytmp && x[j][0] < xtmp) continue;
}
jtype = map[type[j]];
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
ijparam = elem2param[itype][jtype][jtype];
if (rsq > params[ijparam].cutsq) continue;
twobody(&params[ijparam],rsq,fpair,eflag,evdwl);
f[i][0] += delx*fpair;
f[i][1] += dely*fpair;
f[i][2] += delz*fpair;
f[j][0] -= delx*fpair;
f[j][1] -= dely*fpair;
f[j][2] -= delz*fpair;
if (evflag) ev_tally(i,j,nlocal,newton_pair,
evdwl,0.0,fpair,delx,dely,delz);
}
jnumm1 = jnum - 1;
for (jj = 0; jj < jnumm1; jj++) {
j = jlist[jj];
j &= NEIGHMASK;
jtype = map[type[j]];
ijparam = elem2param[itype][jtype][jtype];
delr1[0] = x[j][0] - xtmp;
delr1[1] = x[j][1] - ytmp;
delr1[2] = x[j][2] - ztmp;
rsq1 = delr1[0]*delr1[0] + delr1[1]*delr1[1] + delr1[2]*delr1[2];
if (rsq1 >= params[ijparam].cutsq2) continue;
for (kk = jj+1; kk < jnum; kk++) {
k = jlist[kk];
k &= NEIGHMASK;
ktype = map[type[k]];
ikparam = elem2param[itype][ktype][ktype];
ijkparam = elem2param[itype][jtype][ktype];
delr2[0] = x[k][0] - xtmp;
delr2[1] = x[k][1] - ytmp;
delr2[2] = x[k][2] - ztmp;
rsq2 = delr2[0]*delr2[0] + delr2[1]*delr2[1] + delr2[2]*delr2[2];
if (rsq2 >= params[ikparam].cutsq2) continue;
threebody(&params[ijparam],&params[ikparam],&params[ijkparam],
rsq1,rsq2,delr1,delr2,fj,fk,eflag,evdwl);
f[i][0] -= fj[0] + fk[0];
f[i][1] -= fj[1] + fk[1];
f[i][2] -= fj[2] + fk[2];
f[j][0] += fj[0];
f[j][1] += fj[1];
f[j][2] += fj[2];
f[k][0] += fk[0];
f[k][1] += fk[1];
f[k][2] += fk[2];
if (evflag) ev_tally3(i,j,k,evdwl,0.0,fj,fk,delr1,delr2);
}
}
}
if (vflag_fdotr) virial_fdotr_compute();
}
/* ---------------------------------------------------------------------- */
void PairVashishta::allocate()
{
allocated = 1;
int n = atom->ntypes;
memory->create(setflag,n+1,n+1,"pair:setflag");
memory->create(cutsq,n+1,n+1,"pair:cutsq");
map = new int[n+1];
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairVashishta::settings(int narg, char **arg)
{
if (narg != 0) error->all(FLERR,"Illegal pair_style command");
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairVashishta::coeff(int narg, char **arg)
{
int i,j,n;
if (!allocated) allocate();
if (narg != 3 + atom->ntypes)
error->all(FLERR,"Incorrect args for pair coefficients");
// insure I,J args are * *
if (strcmp(arg[0],"*") != 0 || strcmp(arg[1],"*") != 0)
error->all(FLERR,"Incorrect args for pair coefficients");
// read args that map atom types to elements in potential file
// map[i] = which element the Ith atom type is, -1 if NULL
// nelements = # of unique elements
// elements = list of element names
if (elements) {
for (i = 0; i < nelements; i++) delete [] elements[i];
delete [] elements;
}
elements = new char*[atom->ntypes];
for (i = 0; i < atom->ntypes; i++) elements[i] = NULL;
nelements = 0;
for (i = 3; i < narg; i++) {
if (strcmp(arg[i],"NULL") == 0) {
map[i-2] = -1;
continue;
}
for (j = 0; j < nelements; j++)
if (strcmp(arg[i],elements[j]) == 0) break;
map[i-2] = j;
if (j == nelements) {
n = strlen(arg[i]) + 1;
elements[j] = new char[n];
strcpy(elements[j],arg[i]);
nelements++;
}
}
// read potential file and initialize potential parameters
read_file(arg[2]);
setup();
// clear setflag since coeff() called once with I,J = * *
n = atom->ntypes;
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
setflag[i][j] = 0;
// set setflag i,j for type pairs where both are mapped to elements
int count = 0;
for (int i = 1; i <= n; i++)
for (int j = i; j <= n; j++)
if (map[i] >= 0 && map[j] >= 0) {
setflag[i][j] = 1;
count++;
}
if (count == 0) error->all(FLERR,"Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
init specific to this pair style
------------------------------------------------------------------------- */
void PairVashishta::init_style()
{
if (atom->tag_enable == 0)
error->all(FLERR,"Pair style Vashishta requires atom IDs");
if (force->newton_pair == 0)
error->all(FLERR,"Pair style Vashishta requires newton pair on");
// need a full neighbor list
int irequest = neighbor->request(this);
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->full = 1;
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairVashishta::init_one(int i, int j)
{
if (setflag[i][j] == 0) error->all(FLERR,"All pair coeffs are not set");
return cutmax;
}
/* ---------------------------------------------------------------------- */
void PairVashishta::read_file(char *file)
{
int params_per_line = 17;
char **words = new char*[params_per_line+1];
memory->sfree(params);
params = NULL;
nparams = maxparam = 0;
// open file on proc 0
FILE *fp;
if (comm->me == 0) {
fp = force->open_potential(file);
if (fp == NULL) {
char str[128];
sprintf(str,"Cannot open Vashishta potential file %s",file);
error->one(FLERR,str);
}
}
// read each set of params from potential file
// one set of params can span multiple lines
// store params if all 3 element tags are in element list
int n,nwords,ielement,jelement,kelement;
char line[MAXLINE],*ptr;
int eof = 0;
while (1) {
if (comm->me == 0) {
ptr = fgets(line,MAXLINE,fp);
if (ptr == NULL) {
eof = 1;
fclose(fp);
} else n = strlen(line) + 1;
}
MPI_Bcast(&eof,1,MPI_INT,0,world);
if (eof) break;
MPI_Bcast(&n,1,MPI_INT,0,world);
MPI_Bcast(line,n,MPI_CHAR,0,world);
// strip comment, skip line if blank
if ((ptr = strchr(line,'#'))) *ptr = '\0';
nwords = atom->count_words(line);
if (nwords == 0) continue;
// concatenate additional lines until have params_per_line words
while (nwords < params_per_line) {
n = strlen(line);
if (comm->me == 0) {
ptr = fgets(&line[n],MAXLINE-n,fp);
if (ptr == NULL) {
eof = 1;
fclose(fp);
} else n = strlen(line) + 1;
}
MPI_Bcast(&eof,1,MPI_INT,0,world);
if (eof) break;
MPI_Bcast(&n,1,MPI_INT,0,world);
MPI_Bcast(line,n,MPI_CHAR,0,world);
if ((ptr = strchr(line,'#'))) *ptr = '\0';
nwords = atom->count_words(line);
}
if (nwords != params_per_line)
error->all(FLERR,"Incorrect format in Vashishta potential file");
// words = ptrs to all words in line
nwords = 0;
words[nwords++] = strtok(line," \t\n\r\f");
while ((words[nwords++] = strtok(NULL," \t\n\r\f"))) continue;
// ielement,jelement,kelement = 1st args
// if all 3 args are in element list, then parse this line
// else skip to next entry in file
for (ielement = 0; ielement < nelements; ielement++)
if (strcmp(words[0],elements[ielement]) == 0) break;
if (ielement == nelements) continue;
for (jelement = 0; jelement < nelements; jelement++)
if (strcmp(words[1],elements[jelement]) == 0) break;
if (jelement == nelements) continue;
for (kelement = 0; kelement < nelements; kelement++)
if (strcmp(words[2],elements[kelement]) == 0) break;
if (kelement == nelements) continue;
// load up parameter settings and error check their values
if (nparams == maxparam) {
maxparam += DELTA;
params = (Param *) memory->srealloc(params,maxparam*sizeof(Param),
"pair:params");
}
params[nparams].ielement = ielement;
params[nparams].jelement = jelement;
params[nparams].kelement = kelement;
params[nparams].bigh = atof(words[3]);
params[nparams].eta = atof(words[4]);
params[nparams].zi = atof(words[5]);
params[nparams].zj = atof(words[6]);
params[nparams].lambda1 = atof(words[7]);
params[nparams].bigd = atof(words[8]);
params[nparams].lambda4 = atof(words[9]);
params[nparams].bigw = atof(words[10]);
params[nparams].cut = atof(words[11]);
params[nparams].bigb = atof(words[12]);
params[nparams].gamma = atof(words[13]);
params[nparams].r0 = atof(words[14]);
params[nparams].bigc = atof(words[15]);
params[nparams].costheta = atof(words[16]);
if (params[nparams].bigb < 0.0 || params[nparams].gamma < 0.0 ||
params[nparams].r0 < 0.0 || params[nparams].bigc < 0.0 ||
params[nparams].bigh < 0.0 || params[nparams].eta < 0.0 ||
params[nparams].lambda1 < 0.0 || params[nparams].bigd < 0.0 ||
params[nparams].lambda4 < 0.0 || params[nparams].bigw < 0.0 ||
params[nparams].cut < 0.0)
error->all(FLERR,"Illegal Vashishta parameter");
nparams++;
}
delete [] words;
}
/* ---------------------------------------------------------------------- */
void PairVashishta::setup()
{
int i,j,k,m,n;
// set elem2param for all triplet combinations
// must be a single exact match to lines read from file
// do not allow for ACB in place of ABC
memory->destroy(elem2param);
memory->create(elem2param,nelements,nelements,nelements,"pair:elem2param");
for (i = 0; i < nelements; i++)
for (j = 0; j < nelements; j++)
for (k = 0; k < nelements; k++) {
n = -1;
for (m = 0; m < nparams; m++) {
if (i == params[m].ielement && j == params[m].jelement &&
k == params[m].kelement) {
if (n >= 0) error->all(FLERR,"Potential file has duplicate entry");
n = m;
}
}
if (n < 0) error->all(FLERR,"Potential file is missing an entry");
elem2param[i][j][k] = n;
}
// compute parameter values derived from inputs
// set cutsq using shortcut to reduce neighbor list for accelerated
// calculations. cut must remain unchanged as it is a potential parameter
for (m = 0; m < nparams; m++) {
params[m].cutsq = params[m].cut * params[m].cut;
params[m].cutsq2 = params[m].r0 * params[m].r0;
params[m].lam1inv = 1.0/params[m].lambda1;
params[m].lam4inv = 1.0/params[m].lambda4;
params[m].zizj = params[m].zi*params[m].zj * force->qqr2e;
// Note that bigd does not have 1/2 factor
params[m].mbigd = params[m].bigd;
params[m].heta = params[m].bigh*params[m].eta;
params[m].big2b = 2.0*params[m].bigb;
params[m].big6w = 6.0*params[m].bigw;
params[m].rcinv = 1.0/params[m].cut;
params[m].rc2inv = 1.0/params[m].cutsq;
params[m].rc4inv = params[m].rc2inv*params[m].rc2inv;
params[m].rc6inv = params[m].rc2inv*params[m].rc4inv;
params[m].rceta = pow(params[m].rcinv,params[m].eta);
params[m].lam1rc = params[m].cut*params[m].lam1inv;
params[m].lam4rc = params[m].cut*params[m].lam4inv;
params[m].vrcc2 = params[m].zizj*params[m].rcinv *
exp(-params[m].lam1rc);
params[m].vrcc3 = params[m].mbigd*params[m].rc4inv *
exp(-params[m].lam4rc);
params[m].vrc = params[m].bigh*params[m].rceta +
params[m].vrcc2 - params[m].vrcc3 -
params[m].bigw*params[m].rc6inv;
params[m].dvrc1 = params[m].heta*params[m].rceta*params[m].rcinv;
params[m].dvrc2 = params[m].vrcc2 *
(params[m].rcinv+params[m].lam1inv);
params[m].dvrc3 = params[m].vrcc3 *
(4.0*params[m].rcinv+params[m].lam4inv);
params[m].dvrc4 = params[m].big6w*params[m].rc6inv *
params[m].rcinv;
params[m].dvrc = params[m].dvrc3 + params[m].dvrc4 -
params[m].dvrc1 - params[m].dvrc2;
params[m].c5 = params[m].cut*params[m].dvrc - params[m].vrc;
}
// set cutmax to max of all params
cutmax = 0.0;
for (m = 0; m < nparams; m++) {
if (params[m].cut > cutmax) cutmax = params[m].cut;
if (params[m].r0 > cutmax) cutmax = params[m].r0;
}
}
/* ---------------------------------------------------------------------- */
void PairVashishta::twobody(Param *param, double rsq, double &fforce,
int eflag, double &eng)
{
double r,rinvsq,r4inv,r6inv,reta,lam1r,lam4r;
double vc2,vc3,vo2;
r = sqrt(rsq);
rinvsq = 1.0/rsq;
r4inv = rinvsq*rinvsq;
r6inv = rinvsq*r4inv;
reta = pow(r,-param->eta);
lam1r = r*param->lam1inv;
lam4r = r*param->lam4inv;
vc2 = param->zizj*exp(-lam1r)/r;
vc3 = param->mbigd*r4inv*exp(-lam4r);
vo2 = param->bigh*reta + vc2 - vc3 - param->bigw*r6inv;
fforce = (param->dvrc*r - (4.0*vc3 + lam4r*vc3+param->big6w*r6inv-
param->heta*reta - vc2-lam1r*vc2)) * rinvsq;
if (eflag) eng = vo2 - r*param->dvrc + param->c5;
}
/* ---------------------------------------------------------------------- */
void PairVashishta::threebody(Param *paramij, Param *paramik, Param *paramijk,
double rsq1, double rsq2,
double *delr1, double *delr2,
double *fj, double *fk, int eflag, double &eng)
{
double r1,rinvsq1,rainv1,gsrainv1,gsrainvsq1,expgsrainv1;
double r2,rinvsq2,rainv2,gsrainv2,gsrainvsq2,expgsrainv2;
double rinv12,cs,delcs,delcssq,facexp,facrad,frad1,frad2,pcsinv,pcsinvsq,pcs;
double facang,facang12,csfacang,csfac1,csfac2;
r1 = sqrt(rsq1);
rinvsq1 = 1.0/rsq1;
rainv1 = 1.0/(r1 - paramij->r0);
gsrainv1 = paramij->gamma * rainv1;
gsrainvsq1 = gsrainv1*rainv1/r1;
expgsrainv1 = exp(gsrainv1);
r2 = sqrt(rsq2);
rinvsq2 = 1.0/rsq2;
rainv2 = 1.0/(r2 - paramik->r0);
gsrainv2 = paramik->gamma * rainv2;
gsrainvsq2 = gsrainv2*rainv2/r2;
expgsrainv2 = exp(gsrainv2);
rinv12 = 1.0/(r1*r2);
cs = (delr1[0]*delr2[0] + delr1[1]*delr2[1] + delr1[2]*delr2[2]) * rinv12;
delcs = cs - paramijk->costheta;
delcssq = delcs*delcs;
pcsinv = paramijk->bigc*delcssq + 1.0;
pcsinvsq = pcsinv*pcsinv;
pcs = delcssq/pcsinv;
facexp = expgsrainv1*expgsrainv2;
facrad = paramijk->bigb * facexp * pcs;
frad1 = facrad*gsrainvsq1;
frad2 = facrad*gsrainvsq2;
facang = paramijk->big2b * facexp * delcs/pcsinvsq;
facang12 = rinv12*facang;
csfacang = cs*facang;
csfac1 = rinvsq1*csfacang;
fj[0] = delr1[0]*(frad1+csfac1)-delr2[0]*facang12;
fj[1] = delr1[1]*(frad1+csfac1)-delr2[1]*facang12;
fj[2] = delr1[2]*(frad1+csfac1)-delr2[2]*facang12;
csfac2 = rinvsq2*csfacang;
fk[0] = delr2[0]*(frad2+csfac2)-delr1[0]*facang12;
fk[1] = delr2[1]*(frad2+csfac2)-delr1[1]*facang12;
fk[2] = delr2[2]*(frad2+csfac2)-delr1[2]*facang12;
if (eflag) eng = facrad;
}

122
src/MANYBODY/pair_vashishta.h Executable file
View File

@ -0,0 +1,122 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
http://lammps.sandia.gov, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
PairStyle(vashishta,PairVashishta)
#else
#ifndef LMP_PAIR_Vashishta_H
#define LMP_PAIR_Vashishta_H
#include "pair.h"
namespace LAMMPS_NS {
class PairVashishta : public Pair {
public:
PairVashishta(class LAMMPS *);
virtual ~PairVashishta();
virtual void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
virtual double init_one(int, int);
virtual void init_style();
protected:
struct Param {
double bigb,gamma,r0,bigc,costheta;
double bigh,eta,zi,zj;
double lambda1,bigd,mbigd,lambda4,bigw,cut;
double lam1inv,lam4inv,zizj,heta,big2b,big6w;
double rcinv,rc2inv,rc4inv,rc6inv,rceta;
double cutsq2,cutsq;
double lam1rc,lam4rc,vrcc2,vrcc3,vrc;
double dvrc1,dvrc2,dvrc3,dvrc4,dvrc,c5;
int ielement,jelement,kelement;
};
double cutmax; // max cutoff for all elements
int nelements; // # of unique elements
char **elements; // names of unique elements
int ***elem2param; // mapping from element triplets to parameters
int *map; // mapping from atom types to elements
int nparams; // # of stored parameter sets
int maxparam; // max # of parameter sets
Param *params; // parameter set for an I-J-K interaction
virtual void allocate();
void read_file(char *);
void setup();
void twobody(Param *, double, double &, int, double &);
void threebody(Param *, Param *, Param *, double, double, double *, double *,
double *, double *, int, double &);
};
}
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
E: Pair style Vashishta requires atom IDs
This is a requirement to use the Vashishta potential.
E: Pair style Vashishta requires newton pair on
See the newton command. This is a restriction to use the Vashishta
potential.
E: All pair coeffs are not set
All pair coefficients must be set in the data file or by the
pair_coeff command before running a simulation.
E: Cannot open Vashishta potential file %s
The specified Vashishta potential file cannot be opened. Check that the path
and name are correct.
E: Incorrect format in Vashishta potential file
Incorrect number of words per line in the potential file.
E: Illegal Vashishta parameter
One or more of the coefficients defined in the potential file is
invalid.
E: Potential file has duplicate entry
The potential file for a Vashishta or Tersoff potential has more than
one entry for the same 3 ordered elements.
E: Potential file is missing an entry
The potential file for a Vashishta or Tersoff potential does not have a
needed entry.
*/