fix spring doc page update

This commit is contained in:
Steve Plimpton 2017-01-17 09:02:56 -07:00
parent c31f1e9f22
commit 5cd856c97f
16 changed files with 233 additions and 183 deletions

View File

@ -89,11 +89,7 @@ NOTE: The center of mass of a group of atoms is calculated in
group can straddle a periodic boundary. See the "dump"_dump.html doc
page for a discussion of unwrapped coordinates. It also means that a
spring connecting two groups or a group and the tether point can cross
a periodic boundary and its length be calculated correctly. One
exception is for rigid bodies, which should not be used with the fix
spring command, if the rigid body will cross a periodic boundary.
This is because image flags for rigid bodies are used in a different
way, as explained on the "fix rigid"_fix_rigid.html doc page.
a periodic boundary and its length be calculated correctly.
[Restart, fix_modify, output, run start/stop, minimize info:]

View File

@ -1,163 +1,163 @@
############################################################################
# Input file for investigating twinning nucleation under uniaxial loading with basal plane vector analysis
# Christopher Barrett, March 2013
# This script requires a Mg pair potential file to be in the same directory.
# fname is the file name. It is necessary for loops to work correctly. (See jump command)
variable fname index in.basal
######################################
# POTENTIAL VARIABLES
# lattice parameters and the minimum energy per atom which should be obtained with the current pair potential and homogeneous lattice
variable lx equal 3.181269601
variable b equal sqrt(3)
variable c equal sqrt(8/3)
variable ly equal ${b}*${lx}
variable lz equal ${c}*${lx}
variable pairlocation index almg.liu
variable pairstyle index eam/alloy/opt
######################################
# EQUILIBRATION/DEFORMATION VARIABLES
# eqpress = 10 bar = 1 MPa
# tstep (the timestep) is set to a default value of 0.001 (1 fs)
# seed randomizes the velocity
# srate is the rate of strain in 1/s
# Ndump is the number of timesteps in between each dump of the atom coordinates
variable tstep equal 0.001
variable seed equal 95812384
variable srate equal 1e9
######################################
# INITIALIZATION
units metal
dimension 3
boundary s s s
atom_style atomic
######################################
# ATOM BUILD
atom_modify map array
# lattice custom scale a1 "coordinates of a1" a2 "coordinates of a2" a3 "coordinates of a3" basis "atom1 coordinates" basis "atom2 coordinates" basis "atom3 coordinates" basis "atom4 coordinates" orient x "crystallagraphic orientation of x axis" orient y "crystallagraphic orientation of y axis" z "crystallagraphic orientation of z axis"
lattice custom 3.181269601 a1 1 0 0 a2 0 1.732050808 0 a3 0 0 1.632993162 basis 0.0 0.0 0.0 basis 0.5 0.5 0 basis 0 0.3333333 0.5 basis 0.5 0.833333 0.5 orient x 0 1 1 orient y 1 0 0 orient z 0 1 -1
variable multiple equal 20
variable mx equal "v_lx*v_multiple"
variable my equal "v_ly*v_multiple"
variable mz equal "v_lz*v_multiple"
# the simulation region should be from 0 to a multiple of the periodic boundary in x, y and z.
region whole block 0 ${mz} 0 ${mx} 0 ${my} units box
create_box 2 whole
create_atoms 1 box basis 1 1 basis 2 1 basis 3 1 basis 4 1
region fixed1 block INF INF INF INF INF 10 units box
region fixed2 block INF INF INF INF 100 INF units box
group lower region fixed1
group upper region fixed2
group boundary union upper lower
group mobile subtract all boundary
variable natoms equal "count(all)"
print "# of atoms are: ${natoms}"
######################################
# INTERATOMIC POTENTIAL
pair_style ${pairstyle}
pair_coeff * * ${pairlocation} Mg Mg
######################################
# COMPUTES REQUIRED
compute csym all centro/atom 12
compute eng all pe/atom
compute eatoms all reduce sum c_eng
compute basal all basal/atom
######################################
# MINIMIZATION
# Primarily adjusts the c/a ratio to value predicted by EAM potential
reset_timestep 0
thermo 1
thermo_style custom step pe c_eatoms
min_style cg
minimize 1e-15 1e-15 1000 2000
variable eminimum equal "c_eatoms / count(all)"
print "%%e(it,1)=${eminimum}"
######################################
# EQUILIBRATION
reset_timestep 0
timestep ${tstep}
# atoms are given a random velocity based on a temperature of 100K.
velocity all create 100 ${seed} mom yes rot no
# temperature and pressure are set to 100 and 0
fix 1 all nve
# Set thermo output
thermo 100
thermo_style custom step lx ly lz press pxx pyy pzz pe temp
# Run for at least 2 picosecond (assuming 1 fs timestep)
run 2000
# Loop to run until pressure is below the variable eqpress (defined at beginning of file)
label loopeq
variable eq loop 100
run 250
variable converge equal press
if "${converge} <= 0" then "variable converge equal -press" else "variable converge equal press"
if "${converge} <= 50" then "jump ${fname} breakeq"
next eq
jump ${fname} loopeq
label breakeq
# Store length for strain rate calculations
variable tmp equal "lx"
variable L0 equal ${tmp}
print "Initial Length, L0: ${L0}"
unfix 1
######################################
# DEFORMATION
reset_timestep 0
timestep ${tstep}
# Impose constant strain rate
variable srate1 equal "v_srate / 1.0e10"
velocity upper set 0.0 NULL 0.0 units box
velocity lower set 0.0 NULL 0.0 units box
fix 2 upper setforce 0.0 NULL 0.0
fix 3 lower setforce 0.0 NULL 0.0
fix 1 all nve
# Output strain and stress info to file
# for units metal, pressure is in [bars] = 100 [kPa] = 1/10000 [GPa]
# p2 is in GPa
variable strain equal "(lx - v_L0)/v_L0"
variable p1 equal "v_strain"
variable p2 equal "-pxz/10000"
variable p3 equal "lx"
variable p4 equal "temp"
variable p5 equal "pe"
variable p6 equal "ke"
fix def1 all print 100 "${p1} ${p2} ${p3} ${p4} ${p5} ${p6}" file output.def1.txt screen no
# Dump coordinates to file (for void size calculations)
dump 1 all custom 1000 output.dump.* id x y z c_basal[1] c_basal[2] c_basal[3]
# Display thermo
thermo_style custom step v_strain pxz lx temp pe ke
restart 50000 output.restart
# run deformation for 100000 timesteps (10% strain assuming 1 fs timestep and 1e9/s strainrate)
variable runtime equal 0
label loop
displace_atoms all ramp x 0.0 ${srate1} z 10 100 units box
run 100
variable runtime equal ${runtime}+100
if "${runtime} < 100000" then "jump ${fname} loop"
######################################
# SIMULATION DONE
print "All done"
############################################################################
# Input file for investigating twinning nucleation under uniaxial loading with basal plane vector analysis
# Christopher Barrett, March 2013
# This script requires a Mg pair potential file to be in the same directory.
# fname is the file name. It is necessary for loops to work correctly. (See jump command)
variable fname index in.basal
######################################
# POTENTIAL VARIABLES
# lattice parameters and the minimum energy per atom which should be obtained with the current pair potential and homogeneous lattice
variable lx equal 3.181269601
variable b equal sqrt(3)
variable c equal sqrt(8/3)
variable ly equal ${b}*${lx}
variable lz equal ${c}*${lx}
variable pairlocation index almg.liu
variable pairstyle index eam/alloy/opt
######################################
# EQUILIBRATION/DEFORMATION VARIABLES
# eqpress = 10 bar = 1 MPa
# tstep (the timestep) is set to a default value of 0.001 (1 fs)
# seed randomizes the velocity
# srate is the rate of strain in 1/s
# Ndump is the number of timesteps in between each dump of the atom coordinates
variable tstep equal 0.001
variable seed equal 95812384
variable srate equal 1e9
######################################
# INITIALIZATION
units metal
dimension 3
boundary s s s
atom_style atomic
######################################
# ATOM BUILD
atom_modify map array
# lattice custom scale a1 "coordinates of a1" a2 "coordinates of a2" a3 "coordinates of a3" basis "atom1 coordinates" basis "atom2 coordinates" basis "atom3 coordinates" basis "atom4 coordinates" orient x "crystallagraphic orientation of x axis" orient y "crystallagraphic orientation of y axis" z "crystallagraphic orientation of z axis"
lattice custom 3.181269601 a1 1 0 0 a2 0 1.732050808 0 a3 0 0 1.632993162 basis 0.0 0.0 0.0 basis 0.5 0.5 0 basis 0 0.3333333 0.5 basis 0.5 0.833333 0.5 orient x 0 1 1 orient y 1 0 0 orient z 0 1 -1
variable multiple equal 20
variable mx equal "v_lx*v_multiple"
variable my equal "v_ly*v_multiple"
variable mz equal "v_lz*v_multiple"
# the simulation region should be from 0 to a multiple of the periodic boundary in x, y and z.
region whole block 0 ${mz} 0 ${mx} 0 ${my} units box
create_box 2 whole
create_atoms 1 box basis 1 1 basis 2 1 basis 3 1 basis 4 1
region fixed1 block INF INF INF INF INF 10 units box
region fixed2 block INF INF INF INF 100 INF units box
group lower region fixed1
group upper region fixed2
group boundary union upper lower
group mobile subtract all boundary
variable natoms equal "count(all)"
print "# of atoms are: ${natoms}"
######################################
# INTERATOMIC POTENTIAL
pair_style ${pairstyle}
pair_coeff * * ${pairlocation} Mg Mg
######################################
# COMPUTES REQUIRED
compute csym all centro/atom 12
compute eng all pe/atom
compute eatoms all reduce sum c_eng
compute basal all basal/atom
######################################
# MINIMIZATION
# Primarily adjusts the c/a ratio to value predicted by EAM potential
reset_timestep 0
thermo 1
thermo_style custom step pe c_eatoms
min_style cg
minimize 1e-15 1e-15 1000 2000
variable eminimum equal "c_eatoms / count(all)"
print "%%e(it,1)=${eminimum}"
######################################
# EQUILIBRATION
reset_timestep 0
timestep ${tstep}
# atoms are given a random velocity based on a temperature of 100K.
velocity all create 100 ${seed} mom yes rot no
# temperature and pressure are set to 100 and 0
fix 1 all nve
# Set thermo output
thermo 100
thermo_style custom step lx ly lz press pxx pyy pzz pe temp
# Run for at least 2 picosecond (assuming 1 fs timestep)
run 2000
# Loop to run until pressure is below the variable eqpress (defined at beginning of file)
label loopeq
variable eq loop 100
run 250
variable converge equal press
if "${converge} <= 0" then "variable converge equal -press" else "variable converge equal press"
if "${converge} <= 50" then "jump ${fname} breakeq"
next eq
jump ${fname} loopeq
label breakeq
# Store length for strain rate calculations
variable tmp equal "lx"
variable L0 equal ${tmp}
print "Initial Length, L0: ${L0}"
unfix 1
######################################
# DEFORMATION
reset_timestep 0
timestep ${tstep}
# Impose constant strain rate
variable srate1 equal "v_srate / 1.0e10"
velocity upper set 0.0 NULL 0.0 units box
velocity lower set 0.0 NULL 0.0 units box
fix 2 upper setforce 0.0 NULL 0.0
fix 3 lower setforce 0.0 NULL 0.0
fix 1 all nve
# Output strain and stress info to file
# for units metal, pressure is in [bars] = 100 [kPa] = 1/10000 [GPa]
# p2 is in GPa
variable strain equal "(lx - v_L0)/v_L0"
variable p1 equal "v_strain"
variable p2 equal "-pxz/10000"
variable p3 equal "lx"
variable p4 equal "temp"
variable p5 equal "pe"
variable p6 equal "ke"
fix def1 all print 100 "${p1} ${p2} ${p3} ${p4} ${p5} ${p6}" file output.def1.txt screen no
# Dump coordinates to file (for void size calculations)
dump 1 all custom 1000 output.dump.* id x y z c_basal[1] c_basal[2] c_basal[3]
# Display thermo
thermo_style custom step v_strain pxz lx temp pe ke
restart 50000 output.restart
# run deformation for 100000 timesteps (10% strain assuming 1 fs timestep and 1e9/s strainrate)
variable runtime equal 0
label loop
displace_atoms all ramp x 0.0 ${srate1} z 10 100 units box
run 100
variable runtime equal ${runtime}+100
if "${runtime} < 100000" then "jump ${fname} loop"
######################################
# SIMULATION DONE
print "All done"

View File

@ -51,7 +51,7 @@ void AtomVecAtomic::grow(int n)
if (n == 0) grow_nmax();
else nmax = n;
atom->nmax = nmax;
if (nmax < 0)
if (nmax < 0 || nmax > MAXSMALLINT)
error->one(FLERR,"Per-processor system is too big");
tag = memory->grow(atom->tag,nmax,"atom:tag");

View File

@ -118,7 +118,7 @@ void AtomVecBody::grow(int n)
if (n == 0) grow_nmax();
else nmax = n;
atom->nmax = nmax;
if (nmax < 0)
if (nmax < 0 || nmax > MAXSMALLINT)
error->one(FLERR,"Per-processor system is too big");
tag = memory->grow(atom->tag,nmax,"atom:tag");

View File

@ -53,7 +53,7 @@ void AtomVecCharge::grow(int n)
if (n == 0) grow_nmax();
else nmax = n;
atom->nmax = nmax;
if (nmax < 0)
if (nmax < 0 || nmax > MAXSMALLINT)
error->one(FLERR,"Per-processor system is too big");
tag = memory->grow(atom->tag,nmax,"atom:tag");

View File

@ -72,7 +72,7 @@ void AtomVecEllipsoid::grow(int n)
if (n == 0) grow_nmax();
else nmax = n;
atom->nmax = nmax;
if (nmax < 0)
if (nmax < 0 || nmax > MAXSMALLINT)
error->one(FLERR,"Per-processor system is too big");
tag = memory->grow(atom->tag,nmax,"atom:tag");

View File

@ -144,7 +144,7 @@ void AtomVecHybrid::grow(int n)
if (n == 0) grow_nmax();
else nmax = n;
atom->nmax = nmax;
if (nmax < 0)
if (nmax < 0 || nmax > MAXSMALLINT)
error->one(FLERR,"Per-processor system is too big");
// sub-styles perform all reallocation

View File

@ -83,7 +83,7 @@ void AtomVecLine::grow(int n)
if (n == 0) grow_nmax();
else nmax = n;
atom->nmax = nmax;
if (nmax < 0)
if (nmax < 0 || nmax > MAXSMALLINT)
error->one(FLERR,"Per-processor system is too big");
tag = memory->grow(atom->tag,nmax,"atom:tag");

View File

@ -84,7 +84,7 @@ void AtomVecSphere::grow(int n)
if (n == 0) grow_nmax();
else nmax = n;
atom->nmax = nmax;
if (nmax < 0)
if (nmax < 0 || nmax > MAXSMALLINT)
error->one(FLERR,"Per-processor system is too big");
tag = memory->grow(atom->tag,nmax,"atom:tag");

View File

@ -88,7 +88,7 @@ void AtomVecTri::grow(int n)
if (n == 0) grow_nmax();
else nmax = n;
atom->nmax = nmax;
if (nmax < 0)
if (nmax < 0 || nmax > MAXSMALLINT)
error->one(FLERR,"Per-processor system is too big");
tag = memory->grow(atom->tag,nmax,"atom:tag");

View File

@ -51,10 +51,14 @@ enum{LAYOUT_UNIFORM,LAYOUT_NONUNIFORM,LAYOUT_TILED}; // several files
/* ---------------------------------------------------------------------- */
CommBrick::CommBrick(LAMMPS *lmp) : Comm(lmp),
sendnum(NULL), recvnum(NULL), sendproc(NULL), recvproc(NULL), size_forward_recv(NULL),
size_reverse_send(NULL), size_reverse_recv(NULL), slablo(NULL), slabhi(NULL), multilo(NULL), multihi(NULL),
cutghostmulti(NULL), pbc_flag(NULL), pbc(NULL), firstrecv(NULL), sendlist(NULL), maxsendlist(NULL), buf_send(NULL), buf_recv(NULL)
CommBrick::CommBrick(LAMMPS *lmp) :
Comm(lmp),
sendnum(NULL), recvnum(NULL), sendproc(NULL), recvproc(NULL),
size_forward_recv(NULL),
size_reverse_send(NULL), size_reverse_recv(NULL),
slablo(NULL), slabhi(NULL), multilo(NULL), multihi(NULL),
cutghostmulti(NULL), pbc_flag(NULL), pbc(NULL), firstrecv(NULL),
sendlist(NULL), maxsendlist(NULL), buf_send(NULL), buf_recv(NULL)
{
style = 0;
layout = LAYOUT_UNIFORM;

View File

@ -1102,6 +1102,40 @@ int Domain::closest_image(int i, int j)
return closest;
}
/* ----------------------------------------------------------------------
return local index of atom J or any of its images that is closest to pos
if J is not a valid index like -1, just return it
------------------------------------------------------------------------- */
int Domain::closest_image(double *pos, int j)
{
if (j < 0) return j;
int *sametag = atom->sametag;
double **x = atom->x;
int closest = j;
double delx = pos[0] - x[j][0];
double dely = pos[1] - x[j][1];
double delz = pos[2] - x[j][2];
double rsqmin = delx*delx + dely*dely + delz*delz;
double rsq;
while (sametag[j] >= 0) {
j = sametag[j];
delx = pos[0] - x[j][0];
dely = pos[1] - x[j][1];
delz = pos[2] - x[j][2];
rsq = delx*delx + dely*dely + delz*delz;
if (rsq < rsqmin) {
rsqmin = rsq;
closest = j;
}
}
return closest;
}
/* ----------------------------------------------------------------------
find and return Xj image = periodic image of Xj that is closest to Xi
for triclinic, add/subtract tilt factors in other dims as needed

View File

@ -113,6 +113,7 @@ class Domain : protected Pointers {
void minimum_image(double &, double &, double &);
void minimum_image(double *);
int closest_image(int, int);
int closest_image(double *, int);
void closest_image(const double * const, const double * const,
double * const);
void remap(double *, imageint &);

View File

@ -40,6 +40,7 @@ namespace MathExtra {
inline void sub3(const double *v1, const double *v2, double *ans);
inline double len3(const double *v);
inline double lensq3(const double *v);
inline double distsq3(const double *v1, const double *v2);
inline double dot3(const double *v1, const double *v2);
inline void cross3(const double *v1, const double *v2, double *ans);
@ -265,6 +266,18 @@ inline double MathExtra::lensq3(const double *v)
return v[0]*v[0] + v[1]*v[1] + v[2]*v[2];
}
/* ----------------------------------------------------------------------
ans = distance squared between pts v1 and v2
------------------------------------------------------------------------- */
inline double MathExtra::distsq3(const double *v1, const double *v2)
{
double dx = v1[0] - v2[0];
double dy = v1[1] - v2[1];
double dz = v1[2] - v2[2];
return dx*dx + dy*dy + dz*dz;
}
/* ----------------------------------------------------------------------
dot product of 2 vectors
------------------------------------------------------------------------- */

View File

@ -149,7 +149,7 @@ void NPair::build_setup()
------------------------------------------------------------------------- */
int NPair::exclusion(int i, int j, int itype, int jtype,
int *mask, tagint *molecule) const {
int *mask, tagint *molecule) const {
int m;
if (nex_type && ex_type[itype][jtype]) return 1;

View File

@ -28,8 +28,10 @@ using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
Region::Region(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp),
id(NULL), style(NULL), contact(NULL), list(NULL), xstr(NULL), ystr(NULL), zstr(NULL), tstr(NULL)
Region::Region(LAMMPS *lmp, int narg, char **arg) :
Pointers(lmp),
id(NULL), style(NULL), contact(NULL), list(NULL),
xstr(NULL), ystr(NULL), zstr(NULL), tstr(NULL)
{
int n = strlen(arg[0]) + 1;
id = new char[n];