diff --git a/doc/pair_tersoff.html b/doc/pair_tersoff.html index e38c67bd95..f3ac21bf29 100644 --- a/doc/pair_tersoff.html +++ b/doc/pair_tersoff.html @@ -212,6 +212,6 @@ Condens. Matter, 15, 5649(2003).

-

(Tersoff_2) J. Tersoff, Phys Rev B, 39, 5566 (1989) +

(Tersoff_2) J. Tersoff, Phys Rev B, 39, 5566 (1989); errata (PRB 41, 3248)

diff --git a/doc/pair_tersoff.txt b/doc/pair_tersoff.txt index 6024af80c4..600eee4ac1 100644 --- a/doc/pair_tersoff.txt +++ b/doc/pair_tersoff.txt @@ -206,4 +206,4 @@ appropriate units if your simulation doesn't use "metal" units. Condens. Matter, 15, 5649(2003). :link(Tersoff_2) -[(Tersoff_2)] J. Tersoff, Phys Rev B, 39, 5566 (1989) +[(Tersoff_2)] J. Tersoff, Phys Rev B, 39, 5566 (1989); errata (PRB 41, 3248) diff --git a/examples/ELASTIC/in.elastic b/examples/ELASTIC/in.elastic index 8ce78717ad..4eb09b7d2f 100644 --- a/examples/ELASTIC/in.elastic +++ b/examples/ELASTIC/in.elastic @@ -1,5 +1,7 @@ # Compute elastic constant tensor for a crystal # +# Written by Aidan Thompson (Sandia, athomps@sandia.gov) +# # This script uses the following three include files. # # init.mod (must be modified for different crystal structures) diff --git a/potentials/SiC.tersoff b/potentials/SiC.tersoff index 9a06932846..e42942b752 100644 --- a/potentials/SiC.tersoff +++ b/potentials/SiC.tersoff @@ -1,6 +1,7 @@ # Si and C mixture, parameterized for Tersoff potential # this file is from Rutuparna.Narulkar @ okstate.edu # values are from Phys Rev B, 39, 5566-5568 (1989) +# and errata (PRB 41, 3248) # Tersoff parameters for various elements and mixtures # multiple entries can be added to this file, LAMMPS reads the ones it needs diff --git a/potentials/SiC_1989.tersoff b/potentials/SiC_1989.tersoff new file mode 100644 index 0000000000..7e73697485 --- /dev/null +++ b/potentials/SiC_1989.tersoff @@ -0,0 +1,28 @@ +# Si and C mixture, parameterized for Tersoff potential +# this file is from Saurav Goel - sg258@hw.ac.uk +# J. Tersoff, PRB, 39, 5566 (1989) + errata (PRB 41, 3248) + +# Tersoff parameters for various elements and mixtures +# multiple entries can be added to this file, LAMMPS reads the ones it needs +# these entries are in LAMMPS "metal" units: +# A,B = eV; lambda1,lambda2,lambda3 = 1/Angstroms; R,D = Angstroms +# other quantities are unitless + +# format of a single entry (one or more lines): +# element 1, element 2, element 3, +# m, gamma, lambda3, c, d, costheta0, n, +# beta, lambda2, B, R, D, lambda1, A + +# Tersoff's nomenclature - - - c d h n beta mu B R = (R+S)/2 D = (S-R)/2 Lambda A + +C C C 3 1 0 38049 4.3484 -0.57058 0.72751 1.5724E-07 2.2119 346.74 1.95 0.15 3.4879 1393.6 +Si Si Si 3 1 0 100390 16.217 -0.59825 0.78734 1.1E-06 1.7322 471.18 2.85 0.15 2.4799 1830.8 + + +C Si Si 3 1 0 38049 4.3484 -0.57058 0.72751 1.5724E-07 1.97205 395.145 2.35726 0.1527 2.9839 1597.311 +C Si C 3 1 0 38049 4.3484 -0.57058 0.72751 0 0 0 1.95 0.15 0 0 +C C Si 3 1 0 38049 4.3484 -0.57058 0.72751 0 0 0 2.357260 0.15271 0 0 + +Si C C 3 1 0 100390 16.217 -0.59825 0.78734 1.1E-06 1.97205 395.145 2.35726 0.1527 2.9839 1597.31114 +Si Si C 3 1 0 100390 16.217 -0.59825 0.78734 0 0 0 2.35726 0.15271 0 0 +Si C Si 3 1 0 100390 16.217 -0.59825 0.78734 0 0 0 2.85 0.15 0 0 diff --git a/potentials/SiC_1990.tersoff b/potentials/SiC_1990.tersoff new file mode 100644 index 0000000000..2c511d79a7 --- /dev/null +++ b/potentials/SiC_1990.tersoff @@ -0,0 +1,30 @@ +# Si and C mixture, parameterized for Tersoff potential +# this file is from Saurav Goel - sg258@hw.ac.uk +# J. Tersoff, Phys. Rev. Lett. 64, 1757 (1990). + +# Tersoff parameters for various elements and mixtures +# multiple entries can be added to this file, LAMMPS reads the ones it needs +# these entries are in LAMMPS "metal" units: +# A,B = eV; lambda1,lambda2,lambda3 = 1/Angstroms; R,D = Angstroms +# other quantities are unitless + +# format of a single entry (one or more lines): +# element 1, element 2, element 3, +# m, gamma, lambda3, c, d, costheta0, n, +# beta, lambda2, B, R, D, lambda1, A + +# Tersoff's nomenclature - - - c d h n beta mu B R = (R+S)/2 D = (S-R)/2 Lambda A + +C C C 3 1 0 19981 7.034 -0.33953 0.99054 4.1612E-06 2.3064 389.63 2.5 0 3.4653 1544.8 +Si Si Si 3 1 0 100390 16.217 -0.59825 0.78734 1.1E-06 1.7322 471.18 2.5 0 2.4799 1830.8 + + +C Si Si 3 1 0 19981 7.034 -0.33953 0.99054 4.1612E-06 2.0193 427.269 2.5 0 2.9726 1681.7312 +C Si C 3 1 0 19981 7.034 -0.33953 0.99054 0 0 0 2.5 0 0 0 +C C Si 3 1 0 19981 7.034 -0.33953 0.99054 0 0 0 2.5 0 0 0 + +Si C C 3 1 0 100390 16.217 -0.59825 0.78734 1.1E-06 2.0193 427.269 2.5 0 2.9726 1681.7312 +Si Si C 3 1 0 100390 16.217 -0.59825 0.78734 0 0 0 2.5 0 0 0 +Si C Si 3 1 0 100390 16.217 -0.59825 0.78734 0 0 0 2.5 0 0 0 + + diff --git a/potentials/SiC_1994.tersoff b/potentials/SiC_1994.tersoff new file mode 100644 index 0000000000..f020d59db9 --- /dev/null +++ b/potentials/SiC_1994.tersoff @@ -0,0 +1,28 @@ +# Si and C mixture, parameterized for Tersoff potential +# this file is from Saurav Goel - sg258@hw.ac.uk +#J. Tersoff, Phys. Rev. B49, 16349 (1994). + +# Tersoff parameters for various elements and mixtures +# multiple entries can be added to this file, LAMMPS reads the ones it needs +# these entries are in LAMMPS "metal" units: +# A,B = eV; lambda1,lambda2,lambda3 = 1/Angstroms; R,D = Angstroms +# other quantities are unitless + +# format of a single entry (one or more lines): +# element 1, element 2, element 3, +# m, gamma, lambda3, c, d, costheta0, n, +# beta, lambda2, B, R, D, lambda1, A + +# Tersoff's nomenclature - - - c d h n beta mu B R = (R+S)/2 D = (S-R)/2 Lambda A + +C C C 3 1 0 19981 7.034 -0.33953 0.99054 4.1612E-06 2.3064 389.63 1.95 0.15 3.4653 1544.8 +Si Si Si 3 1 0 100390 16.217 -0.59825 0.78734 1.1E-06 1.7322 471.18 2.85 0.15 2.4799 1830.8 + + +C Si Si 3 1 0 19981 7.034 -0.33953 0.99054 4.1612E-06 2.0193 432.1540 2.357260 0.15271 2.9726 1681.731 +C Si C 3 1 0 19981 7.034 -0.33953 0.99054 0 0 0 1.95 0.15 0 0 +C C Si 3 1 0 19981 7.034 -0.33953 0.99054 0 0 0 2.35726042 0.15271 0 0 + +Si C C 3 1 0 100390 16.217 -0.59825 0.78734 1.1E-06 2.0193 432.154 2.357260 0.152719 2.9726 1681.73125 +Si Si C 3 1 0 100390 16.217 -0.59825 0.78734 0 0 0 2.35726 0.152719 0 0 +Si C Si 3 1 0 100390 16.217 -0.59825 0.78734 0 0 0 2.85 0.15 0 0 diff --git a/src/ASPHERE/fix_nve_asphere.h b/src/ASPHERE/fix_nve_asphere.h index b535a47a14..d9a659a8b4 100755 --- a/src/ASPHERE/fix_nve_asphere.h +++ b/src/ASPHERE/fix_nve_asphere.h @@ -34,9 +34,6 @@ class FixNVEAsphere : public FixNVE { private: double dtq; class AtomVecEllipsoid *avec; - - void richardson(double *, double *, double *); - void omega_from_mq(double *, double *, double *, double *); }; } diff --git a/src/COLLOID/pair_lubricate.cpp b/src/COLLOID/pair_lubricate.cpp index 6d8dec2690..6932620cfe 100644 --- a/src/COLLOID/pair_lubricate.cpp +++ b/src/COLLOID/pair_lubricate.cpp @@ -385,15 +385,15 @@ void PairLubricate::init_style() int irequest = neighbor->request(this); - // require that atom radii are identical within each type require - // monodisperse system with same radii for all types + // require that atom radii are identical within each type + // require monodisperse system with same radii for all types double rad,radtype; for (int i = 1; i <= atom->ntypes; i++) { if (!atom->radius_consistency(i,radtype)) - error->all("Pair colloid requires atoms with same type have same radius"); + error->all("Pair lubricate requires monodisperse particles"); if (i > 1 && radtype != rad) - error->all("Pair colloid requires mono-disperse particles"); + error->all("Pair lubricate requires monodisperse particles"); rad = radtype; } } diff --git a/src/COLLOID/pair_yukawa_colloid.cpp b/src/COLLOID/pair_yukawa_colloid.cpp index 46f66624ee..734be8608e 100644 --- a/src/COLLOID/pair_yukawa_colloid.cpp +++ b/src/COLLOID/pair_yukawa_colloid.cpp @@ -12,7 +12,7 @@ ------------------------------------------------------------------------- */ /* ---------------------------------------------------------------------- - Contributing authors: Randy Schunk (SNL) + Contributing authors: Randy Schunk (Sandia) ------------------------------------------------------------------------- */ #include "math.h" diff --git a/src/atom_vec_ellipsoid.cpp b/src/atom_vec_ellipsoid.cpp index b8e711f142..929e10e09b 100755 --- a/src/atom_vec_ellipsoid.cpp +++ b/src/atom_vec_ellipsoid.cpp @@ -125,7 +125,6 @@ void AtomVecEllipsoid::grow_bonus() /* ---------------------------------------------------------------------- copy atom I info to atom J - if delflag and atom J has bonus data, then delete it ------------------------------------------------------------------------- */ void AtomVecEllipsoid::copy(int i, int j, int delflag) @@ -146,12 +145,17 @@ void AtomVecEllipsoid::copy(int i, int j, int delflag) angmom[j][1] = angmom[i][1]; angmom[j][2] = angmom[i][2]; + // if delflag and atom J has bonus data, then delete it + if (delflag && ellipsoid[j] >= 0) { copy_bonus(nlocal_bonus-1,ellipsoid[j]); nlocal_bonus--; } + + // if atom I has bonus data and not deleting I, repoint I's bonus to J + + if (ellipsoid[i] >= 0 && i != j) bonus[ellipsoid[i]].ilocal = j; ellipsoid[j] = ellipsoid[i]; - if (ellipsoid[j] >= 0) bonus[ellipsoid[j]].ilocal = j; if (atom->nextra_grow) for (int iextra = 0; iextra < atom->nextra_grow; iextra++) @@ -165,24 +169,23 @@ void AtomVecEllipsoid::copy(int i, int j, int delflag) void AtomVecEllipsoid::copy_bonus(int i, int j) { - double *ishape = bonus[i].shape; - double *iquat = bonus[i].quat; - double *jshape = bonus[j].shape; - double *jquat = bonus[j].quat; - jshape[0] = ishape[0]; - jshape[1] = ishape[1]; - jshape[2] = ishape[2]; - jquat[0] = iquat[0]; - jquat[1] = iquat[1]; - jquat[2] = iquat[2]; - jquat[3] = iquat[3]; - int m = bonus[i].ilocal; - bonus[j].ilocal = m; - ellipsoid[m] = j; + memcpy(&bonus[j],&bonus[i],sizeof(Bonus)); + ellipsoid[bonus[j].ilocal] = j; +} + +/* ---------------------------------------------------------------------- + clear ghost info in bonus data + called before ghosts are recommunicated in comm and irregular +------------------------------------------------------------------------- */ + +void AtomVecEllipsoid::clear_bonus() +{ + nghost_bonus = 0; } /* ---------------------------------------------------------------------- set shape values in bonus data for particle I + oriented aligned with xyz axes this may create or delete entry in bonus data ------------------------------------------------------------------------- */ @@ -215,20 +218,10 @@ void AtomVecEllipsoid::set_shape(int i, } } -/* ---------------------------------------------------------------------- - clear ghost info in bonus data - called before ghosts are recommunicated in comm and irregular -------------------------------------------------------------------------- */ - -void AtomVecEllipsoid::clear_bonus() -{ - nghost_bonus = 0; -} - /* ---------------------------------------------------------------------- */ int AtomVecEllipsoid::pack_comm(int n, int *list, double *buf, - int pbc_flag, int *pbc) + int pbc_flag, int *pbc) { int i,j,m; double dx,dy,dz; @@ -767,7 +760,7 @@ void AtomVecEllipsoid::unpack_border(int n, int first, double *buf) type[i] = static_cast (buf[m++]); mask[i] = static_cast (buf[m++]); ellipsoid[i] = static_cast (buf[m++]); - if (ellipsoid[i] < 0) ellipsoid[i] = 0; + if (ellipsoid[i] == 0) ellipsoid[i] = -1; else { j = nlocal_bonus + nghost_bonus; if (j == nmax_bonus) grow_bonus(); @@ -805,9 +798,8 @@ void AtomVecEllipsoid::unpack_border_vel(int n, int first, double *buf) type[i] = static_cast (buf[m++]); mask[i] = static_cast (buf[m++]); ellipsoid[i] = static_cast (buf[m++]); - if (ellipsoid[i] < 0) ellipsoid[i] = 0; + if (ellipsoid[i] == 0) ellipsoid[i] = -1; else { - j = nlocal_bonus + nghost_bonus; if (j == nmax_bonus) grow_bonus(); shape = bonus[j].shape; quat = bonus[j].quat; @@ -842,7 +834,7 @@ int AtomVecEllipsoid::unpack_border_hybrid(int n, int first, double *buf) last = first + n; for (i = first; i < last; i++) { ellipsoid[i] = static_cast (buf[m++]); - if (ellipsoid[i] < 0) ellipsoid[i] = 0; + if (ellipsoid[i] == 0) ellipsoid[i] = -1; else { j = nlocal_bonus + nghost_bonus; if (j == nmax_bonus) grow_bonus(); @@ -891,13 +883,15 @@ int AtomVecEllipsoid::pack_exchange(int i, double *buf) else { buf[m++] = 1; int j = ellipsoid[i]; - buf[m++] = bonus[j].shape[0]; - buf[m++] = bonus[j].shape[1]; - buf[m++] = bonus[j].shape[2]; - buf[m++] = bonus[j].quat[0]; - buf[m++] = bonus[j].quat[1]; - buf[m++] = bonus[j].quat[2]; - buf[m++] = bonus[j].quat[3]; + double *shape = bonus[j].shape; + double *quat = bonus[j].quat; + buf[m++] = shape[0]; + buf[m++] = shape[1]; + buf[m++] = shape[2]; + buf[m++] = quat[0]; + buf[m++] = quat[1]; + buf[m++] = quat[2]; + buf[m++] = quat[3]; } if (atom->nextra_grow) @@ -933,7 +927,8 @@ int AtomVecEllipsoid::unpack_exchange(double *buf) angmom[nlocal][2] = buf[m++]; ellipsoid[nlocal] = static_cast (buf[m++]); - if (ellipsoid[nlocal]) { + if (ellipsoid[nlocal] == 0) ellipsoid[nlocal] = -1; + else { if (nlocal_bonus == nmax_bonus) grow_bonus(); double *shape = bonus[nlocal_bonus].shape; double *quat = bonus[nlocal_bonus].quat; @@ -1057,7 +1052,8 @@ int AtomVecEllipsoid::unpack_restart(double *buf) angmom[nlocal][2] = buf[m++]; ellipsoid[nlocal] = static_cast (buf[m++]); - if (ellipsoid[nlocal]) { + if (ellipsoid[nlocal] == 0) ellipsoid[nlocal] = -1; + else { if (nlocal_bonus == nmax_bonus) grow_bonus(); double *shape = bonus[nlocal_bonus].shape; double *quat = bonus[nlocal_bonus].quat; diff --git a/src/fix_nh.cpp b/src/fix_nh.cpp index 9f5d84b4c8..cf4d5535fb 100644 --- a/src/fix_nh.cpp +++ b/src/fix_nh.cpp @@ -19,6 +19,7 @@ #include "stdlib.h" #include "math.h" #include "fix_nh.h" +#include "math_extra.h" #include "atom.h" #include "force.h" #include "comm.h" @@ -31,7 +32,6 @@ #include "domain.h" #include "memory.h" #include "error.h" -#include "math_extra.h" using namespace LAMMPS_NS; diff --git a/src/fix_store_state.cpp b/src/fix_store_state.cpp index bbae452f96..54c94ad050 100644 --- a/src/fix_store_state.cpp +++ b/src/fix_store_state.cpp @@ -468,7 +468,7 @@ double FixStoreState::memory_usage() void FixStoreState::grow_arrays(int nmax) { - memory->grow(values,nmax,nvalues,"fix_store:values"); + memory->grow(values,nmax,nvalues,"store/state:values"); if (nvalues == 1) { if (nmax) vector_atom = &values[0][0]; else vector_atom = NULL; diff --git a/src/math_extra.cpp b/src/math_extra.cpp index 535a1e576b..c4318e8bbe 100644 --- a/src/math_extra.cpp +++ b/src/math_extra.cpp @@ -428,8 +428,8 @@ void quat_to_mat_trans(const double *quat, double mat[3][3]) /* ---------------------------------------------------------------------- compute space-frame inertia tensor of an ellipsoid - quat = orientiation quaternion of ellipsoid radii = 3 radii of ellipsoid + quat = orientiation quaternion of ellipsoid return symmetric inertia tensor as 6-vector in Voigt notation ------------------------------------------------------------------------- */ @@ -454,6 +454,36 @@ void inertia_ellipsoid(double *radii, double *quat, double mass, inertia[5] = tensor[0][1]; } +/* ---------------------------------------------------------------------- + compute space-frame inertia tensor of a line segment in 2d + length = length of line + theta = orientiation of line + return symmetric inertia tensor as 6-vector in Voigt notation +------------------------------------------------------------------------- */ + +void inertia_line(double length, double theta, double mass, double *inertia) +{ + double p[3][3],ptrans[3][3],itemp[3][3],tensor[3][3]; + double q[4],idiag[3]; + + q[0] = cos(0.5*theta); + q[1] = q[2] = 0.0; + q[3] = sin(0.5*theta); + MathExtra::quat_to_mat(q,p); + MathExtra::quat_to_mat_trans(q,ptrans); + idiag[0] = 0.0; + idiag[1] = 1.0/12.0 * mass * length*length; + idiag[2] = 1.0/12.0 * mass * length*length; + MathExtra::diag_times3(idiag,ptrans,itemp); + MathExtra::times3(p,itemp,tensor); + inertia[0] = tensor[0][0]; + inertia[1] = tensor[1][1]; + inertia[2] = tensor[2][2]; + inertia[3] = tensor[1][2]; + inertia[4] = tensor[0][2]; + inertia[5] = tensor[0][1]; +} + /* ---------------------------------------------------------------------- compute space-frame inertia tensor of a triangle v0,v1,v2 = 3 vertices of triangle diff --git a/src/math_extra.h b/src/math_extra.h index c9a117fe03..3ca98f8f12 100755 --- a/src/math_extra.h +++ b/src/math_extra.h @@ -114,6 +114,8 @@ namespace MathExtra { void inertia_ellipsoid(double *shape, double *quat, double mass, double *inertia); + void inertia_line(double length, double theta, double mass, + double *inertia); void inertia_triangle(double *v0, double *v1, double *v2, double mass, double *inertia); } diff --git a/src/read_data.cpp b/src/read_data.cpp index 7eef7dc4b4..dd3099ab87 100644 --- a/src/read_data.cpp +++ b/src/read_data.cpp @@ -41,7 +41,7 @@ using namespace LAMMPS_NS; #define DELTA 4 // must be 2 or larger // customize for new sections -#define NSECTIONS 23 // change when add to header::section_keywords +#define NSECTIONS 21 // change when add to header::section_keywords #define MIN(a,b) ((a) < (b) ? (a) : (b)) #define MAX(a,b) ((a) > (b) ? (a) : (b))