From 5cd856c97f977ddba1a25dbdbf4c2bb4d137b286 Mon Sep 17 00:00:00 2001 From: Steve Plimpton Date: Tue, 17 Jan 2017 09:02:56 -0700 Subject: [PATCH] fix spring doc page update --- doc/src/fix_spring.txt | 6 +- examples/USER/misc/basal/in.basal | 326 +++++++++++++++--------------- src/atom_vec_atomic.cpp | 2 +- src/atom_vec_body.cpp | 2 +- src/atom_vec_charge.cpp | 2 +- src/atom_vec_ellipsoid.cpp | 2 +- src/atom_vec_hybrid.cpp | 2 +- src/atom_vec_line.cpp | 2 +- src/atom_vec_sphere.cpp | 2 +- src/atom_vec_tri.cpp | 2 +- src/comm_brick.cpp | 12 +- src/domain.cpp | 34 ++++ src/domain.h | 1 + src/math_extra.h | 13 ++ src/npair.cpp | 2 +- src/region.cpp | 6 +- 16 files changed, 233 insertions(+), 183 deletions(-) diff --git a/doc/src/fix_spring.txt b/doc/src/fix_spring.txt index 1d0bd47148..5f94f4cdae 100644 --- a/doc/src/fix_spring.txt +++ b/doc/src/fix_spring.txt @@ -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:] diff --git a/examples/USER/misc/basal/in.basal b/examples/USER/misc/basal/in.basal index fa30c38b74..7fe4c5d2a0 100644 --- a/examples/USER/misc/basal/in.basal +++ b/examples/USER/misc/basal/in.basal @@ -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" diff --git a/src/atom_vec_atomic.cpp b/src/atom_vec_atomic.cpp index c29e04ea87..eda1a33152 100644 --- a/src/atom_vec_atomic.cpp +++ b/src/atom_vec_atomic.cpp @@ -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"); diff --git a/src/atom_vec_body.cpp b/src/atom_vec_body.cpp index 86d3ed8720..ca080ff0b6 100644 --- a/src/atom_vec_body.cpp +++ b/src/atom_vec_body.cpp @@ -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"); diff --git a/src/atom_vec_charge.cpp b/src/atom_vec_charge.cpp index 08c3186a45..a93a29662b 100644 --- a/src/atom_vec_charge.cpp +++ b/src/atom_vec_charge.cpp @@ -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"); diff --git a/src/atom_vec_ellipsoid.cpp b/src/atom_vec_ellipsoid.cpp index 4d1dc01c07..858b89d62b 100644 --- a/src/atom_vec_ellipsoid.cpp +++ b/src/atom_vec_ellipsoid.cpp @@ -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"); diff --git a/src/atom_vec_hybrid.cpp b/src/atom_vec_hybrid.cpp index 7d34931b44..54bd78a83c 100644 --- a/src/atom_vec_hybrid.cpp +++ b/src/atom_vec_hybrid.cpp @@ -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 diff --git a/src/atom_vec_line.cpp b/src/atom_vec_line.cpp index 0e534577f3..3451285370 100644 --- a/src/atom_vec_line.cpp +++ b/src/atom_vec_line.cpp @@ -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"); diff --git a/src/atom_vec_sphere.cpp b/src/atom_vec_sphere.cpp index 7bf4d40825..a72704b4c6 100644 --- a/src/atom_vec_sphere.cpp +++ b/src/atom_vec_sphere.cpp @@ -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"); diff --git a/src/atom_vec_tri.cpp b/src/atom_vec_tri.cpp index 8ffc39cec3..eb87e75b18 100644 --- a/src/atom_vec_tri.cpp +++ b/src/atom_vec_tri.cpp @@ -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"); diff --git a/src/comm_brick.cpp b/src/comm_brick.cpp index 289b11782c..d6cbed40af 100644 --- a/src/comm_brick.cpp +++ b/src/comm_brick.cpp @@ -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; diff --git a/src/domain.cpp b/src/domain.cpp index f627048cfa..52ac9d3d1b 100644 --- a/src/domain.cpp +++ b/src/domain.cpp @@ -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 diff --git a/src/domain.h b/src/domain.h index b8bf1657ce..22e3191231 100644 --- a/src/domain.h +++ b/src/domain.h @@ -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 &); diff --git a/src/math_extra.h b/src/math_extra.h index 75dd492845..a67acce3c6 100644 --- a/src/math_extra.h +++ b/src/math_extra.h @@ -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 ------------------------------------------------------------------------- */ diff --git a/src/npair.cpp b/src/npair.cpp index 6ea4e62550..1fd115d00c 100644 --- a/src/npair.cpp +++ b/src/npair.cpp @@ -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; diff --git a/src/region.cpp b/src/region.cpp index e109b7fd6b..e69fdc79d5 100644 --- a/src/region.cpp +++ b/src/region.cpp @@ -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];