This commit is contained in:
Mike Brown
2011-04-28 15:33:57 -04:00
16 changed files with 168 additions and 54 deletions

View File

@ -212,6 +212,6 @@ Condens. Matter, 15, 5649(2003).
</P>
<A NAME = "Tersoff_2"></A>
<P><B>(Tersoff_2)</B> J. Tersoff, Phys Rev B, 39, 5566 (1989)
<P><B>(Tersoff_2)</B> J. Tersoff, Phys Rev B, 39, 5566 (1989); errata (PRB 41, 3248)
</P>
</HTML>

View File

@ -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)

View File

@ -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)

View File

@ -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

View File

@ -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

View File

@ -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

View File

@ -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

View File

@ -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 *);
};
}

View File

@ -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;
}
}

View File

@ -12,7 +12,7 @@
------------------------------------------------------------------------- */
/* ----------------------------------------------------------------------
Contributing authors: Randy Schunk (SNL)
Contributing authors: Randy Schunk (Sandia)
------------------------------------------------------------------------- */
#include "math.h"

View File

@ -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<int> (buf[m++]);
mask[i] = static_cast<int> (buf[m++]);
ellipsoid[i] = static_cast<int> (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<int> (buf[m++]);
mask[i] = static_cast<int> (buf[m++]);
ellipsoid[i] = static_cast<int> (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<int> (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<int> (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<int> (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;

View File

@ -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;

View File

@ -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;

View File

@ -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

View File

@ -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);
}

View File

@ -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))