From c7170a4296b4d48d80a3e44c8e385bbb826c2ad3 Mon Sep 17 00:00:00 2001
From: athomps <athomps@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Thu, 15 Oct 2015 17:39:43 +0000
Subject: [PATCH] Added Vashishta potential
git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@14117 f3b2605a-c512-4ea7-a41b-209d697bcdaa
---
doc/Section_commands.txt | 1 +
doc/Section_example.txt | 1 +
doc/pair_style.txt | 1 +
doc/pair_vashishta.html | 211 +++++++++
doc/pair_vashishta.txt | 204 +++++++++
examples/README | 1 +
examples/vashishta/InP.vashishta | 38 ++
examples/vashishta/SiO.1990.vashishta | 41 ++
examples/vashishta/data.quartz | 26 ++
examples/vashishta/in.indiumphosphide | 75 ++++
examples/vashishta/in.sio2 | 28 ++
potentials/InP.vashishta | 38 ++
potentials/README | 1 +
potentials/SiC.vashishta | 41 ++
potentials/SiO.1990.vashishta | 41 ++
potentials/SiO.1994.vashishta | 38 ++
potentials/SiO.1997.vashishta | 37 ++
src/MANYBODY/pair_vashishta.cpp | 619 ++++++++++++++++++++++++++
src/MANYBODY/pair_vashishta.h | 122 +++++
19 files changed, 1564 insertions(+)
create mode 100644 doc/pair_vashishta.html
create mode 100644 doc/pair_vashishta.txt
create mode 100755 examples/vashishta/InP.vashishta
create mode 100755 examples/vashishta/SiO.1990.vashishta
create mode 100644 examples/vashishta/data.quartz
create mode 100644 examples/vashishta/in.indiumphosphide
create mode 100644 examples/vashishta/in.sio2
create mode 100755 potentials/InP.vashishta
create mode 100755 potentials/SiC.vashishta
create mode 100755 potentials/SiO.1990.vashishta
create mode 100755 potentials/SiO.1994.vashishta
create mode 100755 potentials/SiO.1997.vashishta
create mode 100755 src/MANYBODY/pair_vashishta.cpp
create mode 100755 src/MANYBODY/pair_vashishta.h
diff --git a/doc/Section_commands.txt b/doc/Section_commands.txt
index fcdc83034d..0553cddfed 100644
--- a/doc/Section_commands.txt
+++ b/doc/Section_commands.txt
@@ -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)
diff --git a/doc/Section_example.txt b/doc/Section_example.txt
index 3cbd54a967..8cda58fb16 100644
--- a/doc/Section_example.txt
+++ b/doc/Section_example.txt
@@ -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:
diff --git a/doc/pair_style.txt b/doc/pair_style.txt
index 760d35b600..1f26b4d522 100644
--- a/doc/pair_style.txt
+++ b/doc/pair_style.txt
@@ -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
diff --git a/doc/pair_vashishta.html b/doc/pair_vashishta.html
new file mode 100644
index 0000000000..f2b42588e4
--- /dev/null
+++ b/doc/pair_vashishta.html
@@ -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>
diff --git a/doc/pair_vashishta.txt b/doc/pair_vashishta.txt
new file mode 100644
index 0000000000..053db937a1
--- /dev/null
+++ b/doc/pair_vashishta.txt
@@ -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
+
diff --git a/examples/README b/examples/README
index c2a9023108..765f88a73a 100644
--- a/examples/README
+++ b/examples/README
@@ -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
diff --git a/examples/vashishta/InP.vashishta b/examples/vashishta/InP.vashishta
new file mode 100755
index 0000000000..9fefd4ef19
--- /dev/null
+++ b/examples/vashishta/InP.vashishta
@@ -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
diff --git a/examples/vashishta/SiO.1990.vashishta b/examples/vashishta/SiO.1990.vashishta
new file mode 100755
index 0000000000..f0f9ab354f
--- /dev/null
+++ b/examples/vashishta/SiO.1990.vashishta
@@ -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
diff --git a/examples/vashishta/data.quartz b/examples/vashishta/data.quartz
new file mode 100644
index 0000000000..6c06201ff0
--- /dev/null
+++ b/examples/vashishta/data.quartz
@@ -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
diff --git a/examples/vashishta/in.indiumphosphide b/examples/vashishta/in.indiumphosphide
new file mode 100644
index 0000000000..41a87345d9
--- /dev/null
+++ b/examples/vashishta/in.indiumphosphide
@@ -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
diff --git a/examples/vashishta/in.sio2 b/examples/vashishta/in.sio2
new file mode 100644
index 0000000000..885eb61a7f
--- /dev/null
+++ b/examples/vashishta/in.sio2
@@ -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
+
diff --git a/potentials/InP.vashishta b/potentials/InP.vashishta
new file mode 100755
index 0000000000..9fefd4ef19
--- /dev/null
+++ b/potentials/InP.vashishta
@@ -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
diff --git a/potentials/README b/potentials/README
index 967ec41b18..d4dba824ca 100644
--- a/potentials/README
+++ b/potentials/README
@@ -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
diff --git a/potentials/SiC.vashishta b/potentials/SiC.vashishta
new file mode 100755
index 0000000000..4444fa24e5
--- /dev/null
+++ b/potentials/SiC.vashishta
@@ -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
+
diff --git a/potentials/SiO.1990.vashishta b/potentials/SiO.1990.vashishta
new file mode 100755
index 0000000000..f0f9ab354f
--- /dev/null
+++ b/potentials/SiO.1990.vashishta
@@ -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
diff --git a/potentials/SiO.1994.vashishta b/potentials/SiO.1994.vashishta
new file mode 100755
index 0000000000..d5e23fd1a8
--- /dev/null
+++ b/potentials/SiO.1994.vashishta
@@ -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
diff --git a/potentials/SiO.1997.vashishta b/potentials/SiO.1997.vashishta
new file mode 100755
index 0000000000..eeadf4d6c0
--- /dev/null
+++ b/potentials/SiO.1997.vashishta
@@ -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
diff --git a/src/MANYBODY/pair_vashishta.cpp b/src/MANYBODY/pair_vashishta.cpp
new file mode 100755
index 0000000000..3a5be84492
--- /dev/null
+++ b/src/MANYBODY/pair_vashishta.cpp
@@ -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(¶ms[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(¶ms[ijparam],¶ms[ikparam],¶ms[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;
+}
diff --git a/src/MANYBODY/pair_vashishta.h b/src/MANYBODY/pair_vashishta.h
new file mode 100755
index 0000000000..540c42d054
--- /dev/null
+++ b/src/MANYBODY/pair_vashishta.h
@@ -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.
+
+*/
--
GitLab