+Syntax: +
+pair_style vashishta ++
Examples: +
+pair_style vashishta +pair_coeff * * SiC.vashishta Si C ++
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, +SiC Vashishta2007, +and InP Branicio2009. +
+The potential for the energy U of a system of atoms is +
+
+where we follow the notation used in 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: +
+See the pair_coeff 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 ++
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: +
+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. +
+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 +shift, table, and tail options. +
+This pair style does not write its information to binary restart +files, 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 command. It does not support the +inner, middle, outer keywords. +
+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 for more info. +
+This pair style requires the newton setting to be "on" +for pair interactions. +
+The Vashishta potential files provided with LAMMPS (see the +potentials directory) are parameterized for metal units. +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: +
+ +Default: none +
+(Vashishta1990) P. Vashishta, R. K. Kalia, J. P. Rino, Phys. Rev. B 41, 12197 (1990). +
+ + +(Vashishta2007) P. Vashishta, R. K. Kalia, A. Nakano, J. P. Rino. J. Appl. Phys. 101, 103515 (2007). +
+ + +(Branicio2009) Branicio, Rino, Gan and Tsuzuki, J. Phys Condensed Matter 21 (2009) 095002 +
+ diff --git a/doc/pair_vashishta.txt b/doc/pair_vashishta.txt new file mode 100644 index 000000000..053db937a --- /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 c2a902310..765f88a73 100644 --- a/examples/README +++ b/examples/README @@ -1,169 +1,170 @@ LAMMPS example problems There are 3 flavors of sub-directories in this file, each with sample problems you can run with LAMMPS. lower-case directories = simple test problems for LAMMPS and its packages upper-case directories = more complex problems USER directory with its own sub-directories = tests for USER packages Each is discussed below. ------------------------------------------ Lower-case directories Each of these sub-directories contains a sample problem you can run with LAMMPS. Most are 2d models so that they run quickly, requiring a few seconds to a few minutes to run on a desktop machine. Each problem has an input script (in.*) and produces a log file (log.*) and (optionally) a dump file (dump.*) or image files (image.*) or movie (movie.mpg) when it runs. Some use a data file (data.*) of initial coordinates as additional input. Some require that you install one or more optional LAMMPS packages. A few sample log file outputs on different machines and different numbers of processors are included in the directories to compare your answers to. E.g. a log file like log.crack.date.foo.P means it ran on P processors of machine "foo" with the dated version of LAMMPS. Note that these problems should get statistically similar answers when run on different machines or different numbers of processors, but not identical answers to those in the log of dump files included here. See the Errors section of the LAMMPS documentation for more discussion. Most of the example input scripts have commented-out lines that produce dump snapshots of the running simulation in any of 3 formats. If you uncomment the dump command in the input script, a text dump file will be produced, which can be animated by various visualization programs (see http://lammps.sandia.gov/viz.html) such as VMD or AtomEye. It can also be animated using the xmovie tool described in the Additional Tools section of the LAMMPS documentation. If you uncomment the dump image command in the input script, and assuming you have built LAMMPS with a JPG library, JPG snapshot images will be produced when the simulation runs. They can be quickly post-processed into a movie using commands described on the dump image doc page. If you uncomment the dump movie command in the input script, and assuming you have built LAMMPS with the FFMPEG library, an MPG movie will be produced when the simulation runs. The movie file can be played using various viewers, such as mplayer or QuickTime. Animations of many of these examples can be viewed on the Movies section of the LAMMPS WWW Site. These are the sample problems and their output in the various sub-directories: accelerate: use of all the various accelerator packages balance: dynamic load balancing, 2d system body: body particles, 2d system colloid: big colloid particles in a small particle solvent, 2d system coreshell: adiabatic core/shell model comb: models using the COMB potential crack: crack propagation in a 2d solid deposit: deposition of atoms and molecules onto a 3d substrate dipole: point dipolar particles, 2d system dreiding: methanol via Dreiding FF eim: NaCl using the EIM potential ellipse: ellipsoidal particles in spherical solvent, 2d system flow: Couette and Poiseuille flow in a 2d channel friction: frictional contact of spherical asperities between 2d surfaces hugoniostat: Hugoniostat shock dynamics indent: spherical indenter into a 2d solid kim: use of potentials in Knowledge Base for Interatomic Models (KIM) meam: MEAM test for SiC and shear (same as shear examples) melt: rapid melt of 3d LJ system micelle: self-assembly of small lipid-like molecules into 2d bilayers min: energy minimization of 2d LJ melt msst: MSST shock dynamics nb3b: use of nonbonded 3-body harmonic pair style neb: nudged elastic band (NEB) calculation for barrier finding nemd: non-equilibrium MD of 2d sheared system obstacle: flow around two voids in a 2d channel peptide: dynamics of a small solvated peptide chain (5-mer) peri: Peridynamic model of cylinder impacted by indenter pour: pouring of granular particles into a 3d box, then chute flow prd: parallel replica dynamics of vacancy diffusion in bulk Si python: use of PYTHON package to invoke Python code from input script qeq: use of QEQ pacakge for charge equilibration reax: RDX and TATB models using the ReaxFF rigid: rigid bodies modeled as independent or coupled shear: sideways shear applied to 2d solid, with and without a void snap: use of SNAP potential for Ta 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 LAMMPS executable "lmp_mpi" with all the packages needed by all the examples, with the exception of the accelerate sub-directory. See the accelerate/README for Make.py commands suitable for its example scripts. cd src Make.py -j 16 -p none std no-lib reax meam poems reaxc orig -a lib-all mpi Here is how you might run and visualize one of the sample problems: cd indent cp ../../src/lmp_mpi . # copy LAMMPS executable to this dir lmp_mpi < in.indent # run the problem Running the simulation produces the files {dump.indent} and {log.lammps}. You can visualize the dump file as follows: ../../tools/xmovie/xmovie -scale dump.indent If you uncomment the dump image line(s) in the input script a series of JPG images will be produced by the run. These can be viewed individually or turned into a movie or animated by tools like ImageMagick or QuickTime or various Windows-based tools. See the dump image doc page for more details. E.g. this Imagemagick command would create a GIF file suitable for viewing in a browser. % convert -loop 1 *.jpg foo.gif ------------------------------------------ Upper-case directories The ASPHERE directory has examples of how to model aspherical particles with or without solvent, in 3 styles LAMMPS provides. Namely point ellipsoids, rigid bodies, and generalized aspherical bodies built from line/triangle surface facets in 2d/3d. See the ASPHERE/README file to get started. The COUPLE directory has examples of how to use LAMMPS as a library, either by itself or in tandem with another code or library. See the COUPLE/README file to get started. The ELASTIC directory has an example script for computing elastic constants at zero temperature, using an Si example. See the ELASTIC/in.elastic file for more info. The ELASTIC_T directory has an example script for computing elastic constants at finite temperature, using an Si example. See the ELASTIC_T/in.elastic file for more info. The KAPPA directory has example scripts for computing the thermal conductivity (kappa) of a LJ liquid using 4 different methods. See the KAPPA/README file for more info. The MC directory has an example script for using LAMMPS as an energy-evaluation engine in a iterative Monte Carlo energy-relaxation loop. The USER directory contains subdirectories of user-provided example scripts for ser packages. See the README files in those directories for more info. See the doc/Section_start.html file for more info about installing and building user packages. The VISCOSITY directory has example scripts for computing the viscosity of a LJ liquid using 4 different methods. See the VISCOSITY/README file for more info. diff --git a/examples/vashishta/InP.vashishta b/examples/vashishta/InP.vashishta new file mode 100755 index 000000000..9fefd4ef1 --- /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 000000000..f0f9ab354 --- /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 000000000..6c06201ff --- /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 000000000..41a87345d --- /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 000000000..885eb61a7 --- /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 000000000..9fefd4ef1 --- /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 967ec41b1..d4dba824c 100644 --- a/potentials/README +++ b/potentials/README @@ -1,102 +1,103 @@ This directory contains potential files for different elements and alloys, as used by LAMMPS for various pair styles. See the description of the "pair_style" and "pair_coeff" commands for details of the file formats and the various styles in LAMMPS that read these files. IMPORTANT NOTE: These files are provided primarily to demonstrate the different types of interatomic potentials that LAMMPS supports. Each file has a header line with a date for when it was added to the LAMMPS distribution. Also a citation and contact info for the person who contributed it to LAMMPS (if we remember who that is). This info is not meant to "guarantee" that the potential is correct. I.e. that the contributor transcribed the info from the paper correctly or that the paper itself had no errors. In many cases (but not all), we or other LAMMPS users have confirmed that when the potential file is used with the current version of LAMMPS, it reproduces results in the cited publication. In some cases, this accuracy check may require other parameters not contained in the potential file to be specified as part of a LAMMPS simulation, e.g. a distance cutoff. Also, for particular materials and applications modeled with a pair style coded in LAMMPS, a different potential file may be more suitable than the one provided here. For best results when choosing a potential, you should do a thorough search of published literature and on-line databases such as the Interatomic Potentials Repository Project (NIST) or the Knowledgebase of Interatomic Models (KIM). Whatever potential you choose for your application, you should verify that you have defined it and are using it correctly in LAMMPS, by comparing with published results for that potential. 2nd IMPORTANT NOTE: The DATE field in the first line of each of these files is printed to the screen and log file when it is read by a LAMMPS input script. If an updated or corrected version of the same potential file is later added to the LAMMPS distribution, then a new DATE will be added to the file. This means you can "diff" an old and new log file and see that the potential file changed, which could affect your simulation results. A small amount of metadata is included in the first line of each file in order to track the provenance of each file. The metadata is indicated by a keyword followed by white space, followed by the metadata, followed by whitespace. The metadata is intended to be straightforward and human-readable, while still conforming to a standard format. DATE: Format is "yyyy-mm-dd". This indicates the date of a significant change to the file. Multiple entries can appear in reverse chronological order. As described above, the first of these will be printed to the screen and log file when it is read by a LAMMPS input script. CONTRIBUTOR: Format is "name[, email address]". This indicates the person who contributed the file and/or who is best able to provide more details about its provenance. CITATION: Format is "surname[[, surname] and surname], Publication abbreviation with spaces and no periods, volume, page[-page], (year)" COMMENT: This one is optional and is used to hold any other text that can not go elsewhere. If the first line of the file is always skipped by the file reader, then the first line should begin with the DATE keyword. If the file format supports comment lines, then the first line should be a comment line with the metadata e.g. "# DATE: 2010-01-01..." If the first line of the file is required to begin with data, then the metadata will be appended to the first line e.g. "7 DATE: 2010-01-01..." The prefix of each file indicates the element(s) it is parameterized for. An additional lower-case identification tag may be appended. Si = Silicon SiC = Silicon and Carbon Au_u3 = Gold universal 3 The suffix of each file indicates the pair style it is used with: adp ADP angular dependent potential airebo AI-REBO and REBO potentials bop.table BOP potential, tabulated form cdeam concentration-dependent EAM comb COMB potential comb3 COMB3 potential eam embedded atom method (EAM) single element, DYNAMO funcfl format eam.alloy EAM multi-element alloy, DYNAMO setfl format eam.fs Finnis-Sinclair EAM format (single element or alloy) edip EDIP potential for silicon-based materials eim embedded-ion method (EIM) potential lcbop LCBOP long-range bond-order potential meam modified EAM (MEAM) library and individual elements/alloys meam.spline modified EAM (MEAM) spline potential meam.sw.spline modified EAM (MEAM) Stillinger-Weber spline potential nb3b.harmonic nonbonded 3-body harmonic potential poly polymorphic 3-body potential reax ReaxFF potential (see README.reax for more info) snap SNAP potential snapcoeff SNAP potential snapparam SNAP potential streitz Coulombic portion of Streitz-Mintmire potential 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 000000000..4444fa24e --- /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 000000000..f0f9ab354 --- /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 000000000..d5e23fd1a --- /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 000000000..eeadf4d6c --- /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 000000000..3a5be8449 --- /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 000000000..540c42d05 --- /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. + +*/