diff --git a/doc/src/Section_howto.txt b/doc/src/Section_howto.txt index 33bba562ed..5e8f676209 100644 --- a/doc/src/Section_howto.txt +++ b/doc/src/Section_howto.txt @@ -1936,18 +1936,22 @@ documentation in the src/library.cpp file for details, including which quantities can be queried by name: void *lammps_extract_global(void *, char *) +void lammps_extract_box(void *, double *, double *, + double *, double *, double *, int *, int *) void *lammps_extract_atom(void *, char *) void *lammps_extract_compute(void *, char *, int, int) void *lammps_extract_fix(void *, char *, int, int, int, int) void *lammps_extract_variable(void *, char *, char *) :pre -int lammps_set_variable(void *, char *, char *) -double lammps_get_thermo(void *, char *) :pre +void lammps_reset_box(void *, double *, double *, double, double, double) +int lammps_set_variable(void *, char *, char *) :pre +double lammps_get_thermo(void *, char *) int lammps_get_natoms(void *) void lammps_gather_atoms(void *, double *) void lammps_scatter_atoms(void *, double *) :pre -void lammps_create_atoms(void *, int, tagint *, int *, double *, double *) :pre +void lammps_create_atoms(void *, int, tagint *, int *, double *, double *, + imageint *, int) :pre The extract functions return a pointer to various global or per-atom quantities stored in LAMMPS or to values calculated by a compute, fix, @@ -1957,10 +1961,16 @@ the other extract functions, the underlying storage may be reallocated as LAMMPS runs, so you need to re-call the function to assure a current pointer or returned value(s). +The lammps_reset_box() function resets the size and shape of the +simulation box, e.g. as part of restoring a previously extracted and +saved state of a simulation. + The lammps_set_variable() function can set an existing string-style variable to a new string value, so that subsequent LAMMPS commands can -access the variable. The lammps_get_thermo() function returns the -current value of a thermo keyword as a double. +access the variable. + +The lammps_get_thermo() function returns the current value of a thermo +keyword as a double precision value. The lammps_get_natoms() function returns the total number of atoms in the system and can be used by the caller to allocate space for the @@ -1973,10 +1983,13 @@ passed by the caller, to each atom owned by individual processors. The lammps_create_atoms() function takes a list of N atoms as input with atom types and coords (required), an optionally atom IDs and -velocities. It uses the coords of each atom to assign it as a new -atom to the processor that owns it. Additional properties for the new -atoms can be assigned via the lammps_scatter_atoms() or -lammps_extract_atom() functions. +velocities and image flags. It uses the coords of each atom to assign +it as a new atom to the processor that owns it. This function is +useful to add atoms to a simulation or (in tandem with +lammps_reset_box()) to restore a previously extracted and saved state +of a simulation. Additional properties for the new atoms can then be +assigned via the lammps_scatter_atoms() or lammps_extract_atom() +functions. The examples/COUPLE and python directories have example C++ and C and Python codes which show how a driver code can link to LAMMPS as a diff --git a/src/KSPACE/ewald_disp.cpp b/src/KSPACE/ewald_disp.cpp index 2b0444c84f..467a748d08 100644 --- a/src/KSPACE/ewald_disp.cpp +++ b/src/KSPACE/ewald_disp.cpp @@ -195,17 +195,6 @@ void EwaldDisp::init() g_ewald = accuracy*sqrt(natoms*(*cutoff)*shape_det(domain->h)) / (2.0*q2); if (g_ewald >= 1.0) g_ewald = (1.35 - 0.15*log(accuracy))/(*cutoff); else g_ewald = sqrt(-log(g_ewald)) / (*cutoff); - } - else if (function[1] || function[2]) { - //Try Newton Solver - //Use old method to get guess - g_ewald = (1.35 - 0.15*log(accuracy))/ *cutoff; - - double g_ewald_new = - NewtonSolve(g_ewald,(*cutoff),natoms,shape_det(domain->h),b2); - if (g_ewald_new > 0.0) g_ewald = g_ewald_new; - else error->warning(FLERR,"Ewald/disp Newton solver failed, " - "using old method to estimate g_ewald"); } else if (function[3]) { //Try Newton Solver //Use old method to get guess @@ -215,6 +204,16 @@ void EwaldDisp::init() if (g_ewald_new > 0.0) g_ewald = g_ewald_new; else error->warning(FLERR,"Ewald/disp Newton solver failed, " "using old method to estimate g_ewald"); + } else if (function[1] || function[2]) { + //Try Newton Solver + //Use old method to get guess + g_ewald = (1.35 - 0.15*log(accuracy))/ *cutoff; + + double g_ewald_new = + NewtonSolve(g_ewald,(*cutoff),natoms,shape_det(domain->h),b2); + if (g_ewald_new > 0.0) g_ewald = g_ewald_new; + else error->warning(FLERR,"Ewald/disp Newton solver failed, " + "using old method to estimate g_ewald"); } } @@ -708,6 +707,8 @@ void EwaldDisp::compute(int eflag, int vflag) compute_virial(); compute_virial_dipole(); compute_virial_peratom(); + + if (slabflag) compute_slabcorr(); } @@ -974,7 +975,6 @@ void EwaldDisp::compute_energy() } } for (int k=0; kminimize->init(); // cannot use PRD with a changing box + // removing this restriction would require saving/restoring box params if (domain->box_change) error->all(FLERR,"Cannot use PRD with a changing box"); diff --git a/src/RIGID/compute_rigid_local.cpp b/src/RIGID/compute_rigid_local.cpp index b746e5e6d5..9b2f7c7da1 100644 --- a/src/RIGID/compute_rigid_local.cpp +++ b/src/RIGID/compute_rigid_local.cpp @@ -33,7 +33,8 @@ enum{ID,MOL,MASS,X,Y,Z,XU,YU,ZU,VX,VY,VZ,FX,FY,FZ,IX,IY,IZ, /* ---------------------------------------------------------------------- */ ComputeRigidLocal::ComputeRigidLocal(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg), rstyle(NULL), idrigid(NULL), fixrigid(NULL) + Compute(lmp, narg, arg), + rstyle(NULL), idrigid(NULL), fixrigid(NULL), vlocal(NULL), alocal(NULL) { if (narg < 5) error->all(FLERR,"Illegal compute rigid/local command"); @@ -88,16 +89,16 @@ ComputeRigidLocal::ComputeRigidLocal(LAMMPS *lmp, int narg, char **arg) : } ncount = nmax = 0; - vector = NULL; - array = NULL; + vlocal = NULL; + alocal = NULL; } /* ---------------------------------------------------------------------- */ ComputeRigidLocal::~ComputeRigidLocal() { - memory->destroy(vector); - memory->destroy(array); + memory->destroy(vlocal); + memory->destroy(alocal); delete [] idrigid; delete [] rstyle; } @@ -169,8 +170,8 @@ int ComputeRigidLocal::compute_rigid(int flag) body = &fixrigid->body[ibody]; if (flag) { - if (nvalues == 1) ptr = &vector[m]; - else ptr = array[m]; + if (nvalues == 1) ptr = &vlocal[m]; + else ptr = alocal[m]; for (n = 0; n < nvalues; n++) { switch (rstyle[n]) { @@ -293,18 +294,18 @@ int ComputeRigidLocal::compute_rigid(int flag) void ComputeRigidLocal::reallocate(int n) { - // grow vector or array + // grow vector_local or array_local while (nmax < n) nmax += DELTA; if (nvalues == 1) { - memory->destroy(vector); - memory->create(vector,nmax,"rigid/local:vector"); - vector_local = vector; + memory->destroy(vlocal); + memory->create(vlocal,nmax,"rigid/local:vector_local"); + vector_local = vlocal; } else { - memory->destroy(array); - memory->create(array,nmax,nvalues,"rigid/local:array"); - array_local = array; + memory->destroy(alocal); + memory->create(alocal,nmax,nvalues,"rigid/local:array_local"); + array_local = alocal; } } diff --git a/src/RIGID/compute_rigid_local.h b/src/RIGID/compute_rigid_local.h index d7d545d1fd..0672f5975c 100644 --- a/src/RIGID/compute_rigid_local.h +++ b/src/RIGID/compute_rigid_local.h @@ -41,6 +41,8 @@ class ComputeRigidLocal : public Compute { class FixRigidSmall *fixrigid; int nmax; + double *vlocal; + double **alocal; int compute_rigid(int); void reallocate(int); diff --git a/src/compute.cpp b/src/compute.cpp index 88ae70c067..d1db42a24e 100644 --- a/src/compute.cpp +++ b/src/compute.cpp @@ -38,10 +38,11 @@ int Compute::instance_total = 0; /* ---------------------------------------------------------------------- */ -Compute::Compute(LAMMPS *lmp, int narg, char **arg) : Pointers(lmp), +Compute::Compute(LAMMPS *lmp, int narg, char **arg) : + Pointers(lmp), id(NULL), style(NULL), vector(NULL), array(NULL), vector_atom(NULL), - array_atom(NULL), vector_local(NULL), array_local(NULL), + array_atom(NULL), vector_local(NULL), array_local(NULL), extlist(NULL), tlist(NULL), vbiasall(NULL) { instance_me = instance_total++; diff --git a/src/compute_angle_local.cpp b/src/compute_angle_local.cpp index c5478780ae..1e249554dc 100644 --- a/src/compute_angle_local.cpp +++ b/src/compute_angle_local.cpp @@ -33,7 +33,8 @@ using namespace MathConst; /* ---------------------------------------------------------------------- */ ComputeAngleLocal::ComputeAngleLocal(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg) + Compute(lmp, narg, arg), + vlocal(NULL), alocal(NULL) { if (narg < 4) error->all(FLERR,"Illegal compute angle/local command"); @@ -55,14 +56,16 @@ ComputeAngleLocal::ComputeAngleLocal(LAMMPS *lmp, int narg, char **arg) : } nmax = 0; + vlocal = NULL; + alocal = NULL; } /* ---------------------------------------------------------------------- */ ComputeAngleLocal::~ComputeAngleLocal() { - memory->destroy(vector); - memory->destroy(array); + memory->destroy(vlocal); + memory->destroy(alocal); } /* ---------------------------------------------------------------------- */ @@ -130,12 +133,12 @@ int ComputeAngleLocal::compute_angles(int flag) if (flag) { if (nvalues == 1) { - if (tflag >= 0) tbuf = vector; - if (eflag >= 0) ebuf = vector; + if (tflag >= 0) tbuf = vlocal; + if (eflag >= 0) ebuf = vlocal; } else { - if (tflag >= 0 && array) tbuf = &array[0][tflag]; + if (tflag >= 0 && alocal) tbuf = &alocal[0][tflag]; else tbuf = NULL; - if (eflag >= 0 && array) ebuf = &array[0][eflag]; + if (eflag >= 0 && alocal) ebuf = &alocal[0][eflag]; else ebuf = NULL; } } @@ -218,18 +221,18 @@ int ComputeAngleLocal::compute_angles(int flag) void ComputeAngleLocal::reallocate(int n) { - // grow vector or array and indices array + // grow vector_local or array_local while (nmax < n) nmax += DELTA; if (nvalues == 1) { - memory->destroy(vector); - memory->create(vector,nmax,"bond/local:vector"); - vector_local = vector; + memory->destroy(vlocal); + memory->create(vlocal,nmax,"angle/local:vector_local"); + vector_local = vlocal; } else { - memory->destroy(array); - memory->create(array,nmax,nvalues,"bond/local:array"); - array_local = array; + memory->destroy(alocal); + memory->create(alocal,nmax,nvalues,"angle/local:array_local"); + array_local = alocal; } } diff --git a/src/compute_angle_local.h b/src/compute_angle_local.h index cf5d2fbe37..6f59d7f097 100644 --- a/src/compute_angle_local.h +++ b/src/compute_angle_local.h @@ -37,6 +37,8 @@ class ComputeAngleLocal : public Compute { int ncount; int nmax; + double *vlocal; + double **alocal; int compute_angles(int); void reallocate(int); diff --git a/src/compute_bond_local.cpp b/src/compute_bond_local.cpp index 424059aa86..e145ad373d 100644 --- a/src/compute_bond_local.cpp +++ b/src/compute_bond_local.cpp @@ -37,7 +37,7 @@ enum{DIST,VELVIB,OMEGA,ENGTRANS,ENGVIB,ENGROT,ENGPOT,FORCE}; ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg), - bstyle(NULL) + bstyle(NULL), vlocal(NULL), alocal(NULL) { if (narg < 4) error->all(FLERR,"Illegal compute bond/local command"); @@ -78,14 +78,16 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) : } nmax = 0; + vlocal = NULL; + alocal = NULL; } /* ---------------------------------------------------------------------- */ ComputeBondLocal::~ComputeBondLocal() { - memory->destroy(vector); - memory->destroy(array); + memory->destroy(vlocal); + memory->destroy(alocal); delete [] bstyle; } @@ -292,8 +294,8 @@ int ComputeBondLocal::compute_bonds(int flag) engrot *= mvv2e; } - if (nvalues == 1) ptr = &vector[m]; - else ptr = array[m]; + if (nvalues == 1) ptr = &vlocal[m]; + else ptr = alocal[m]; for (n = 0; n < nvalues; n++) { switch (bstyle[n]) { @@ -373,18 +375,18 @@ void ComputeBondLocal::unpack_forward_comm(int n, int first, double *buf) void ComputeBondLocal::reallocate(int n) { - // grow vector or array and indices array + // grow vector_local or array_local while (nmax < n) nmax += DELTA; if (nvalues == 1) { - memory->destroy(vector); - memory->create(vector,nmax,"bond/local:vector"); - vector_local = vector; + memory->destroy(vlocal); + memory->create(vlocal,nmax,"bond/local:vector_local"); + vector_local = vlocal; } else { - memory->destroy(array); - memory->create(array,nmax,nvalues,"bond/local:array"); - array_local = array; + memory->destroy(alocal); + memory->create(alocal,nmax,nvalues,"bond/local:array_local"); + array_local = alocal; } } diff --git a/src/compute_bond_local.h b/src/compute_bond_local.h index 64cafc75e0..9716749917 100644 --- a/src/compute_bond_local.h +++ b/src/compute_bond_local.h @@ -41,6 +41,8 @@ class ComputeBondLocal : public Compute { int singleflag,velflag,ghostvelflag,initflag; int nmax; + double *vlocal; + double **alocal; int compute_bonds(int); void reallocate(int); diff --git a/src/compute_dihedral_local.cpp b/src/compute_dihedral_local.cpp index d8fc1b792c..7a21d0f9d7 100644 --- a/src/compute_dihedral_local.cpp +++ b/src/compute_dihedral_local.cpp @@ -34,7 +34,8 @@ using namespace MathConst; /* ---------------------------------------------------------------------- */ ComputeDihedralLocal::ComputeDihedralLocal(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg) + Compute(lmp, narg, arg), + vlocal(NULL), alocal(NULL) { if (narg < 4) error->all(FLERR,"Illegal compute dihedral/local command"); @@ -56,14 +57,16 @@ ComputeDihedralLocal::ComputeDihedralLocal(LAMMPS *lmp, int narg, char **arg) : } nmax = 0; + vlocal = NULL; + alocal = NULL; } /* ---------------------------------------------------------------------- */ ComputeDihedralLocal::~ComputeDihedralLocal() { - memory->destroy(vector); - memory->destroy(array); + memory->destroy(vlocal); + memory->destroy(alocal); } /* ---------------------------------------------------------------------- */ @@ -129,9 +132,9 @@ int ComputeDihedralLocal::compute_dihedrals(int flag) if (flag) { if (nvalues == 1) { - if (pflag >= 0) pbuf = vector; + if (pflag >= 0) pbuf = vlocal; } else { - if (pflag >= 0 && array) pbuf = &array[0][pflag]; + if (pflag >= 0 && alocal) pbuf = &alocal[0][pflag]; else pbuf = NULL; } } @@ -229,18 +232,18 @@ int ComputeDihedralLocal::compute_dihedrals(int flag) void ComputeDihedralLocal::reallocate(int n) { - // grow vector or array and indices array + // grow vector_local or array_local while (nmax < n) nmax += DELTA; if (nvalues == 1) { - memory->destroy(vector); - memory->create(vector,nmax,"bond/local:vector"); - vector_local = vector; + memory->destroy(vlocal); + memory->create(vlocal,nmax,"dihedral/local:vector_local"); + vector_local = vlocal; } else { - memory->destroy(array); - memory->create(array,nmax,nvalues,"bond/local:array"); - array_local = array; + memory->destroy(alocal); + memory->create(alocal,nmax,nvalues,"dihedral/local:array_local"); + array_local = alocal; } } diff --git a/src/compute_dihedral_local.h b/src/compute_dihedral_local.h index cc4c16ca79..c077e52dbd 100644 --- a/src/compute_dihedral_local.h +++ b/src/compute_dihedral_local.h @@ -37,6 +37,8 @@ class ComputeDihedralLocal : public Compute { int ncount; int nmax; + double *vlocal; + double **alocal; int compute_dihedrals(int); void reallocate(int); diff --git a/src/compute_gyration_chunk.cpp b/src/compute_gyration_chunk.cpp index f7e774566b..ddec1315a2 100644 --- a/src/compute_gyration_chunk.cpp +++ b/src/compute_gyration_chunk.cpp @@ -28,8 +28,8 @@ using namespace LAMMPS_NS; ComputeGyrationChunk::ComputeGyrationChunk(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg), - idchunk(NULL), massproc(NULL), masstotal(NULL), com(NULL), comall(NULL), rg(NULL), - rgall(NULL), rgt(NULL), rgtall(NULL) + idchunk(NULL), massproc(NULL), masstotal(NULL), com(NULL), comall(NULL), + rg(NULL), rgall(NULL), rgt(NULL), rgtall(NULL) { if (narg < 4) error->all(FLERR,"Illegal compute gyration/chunk command"); diff --git a/src/compute_improper_local.cpp b/src/compute_improper_local.cpp index 05a0cb4291..e183405956 100644 --- a/src/compute_improper_local.cpp +++ b/src/compute_improper_local.cpp @@ -35,7 +35,8 @@ using namespace MathConst; /* ---------------------------------------------------------------------- */ ComputeImproperLocal::ComputeImproperLocal(LAMMPS *lmp, int narg, char **arg) : - Compute(lmp, narg, arg) + Compute(lmp, narg, arg), + vlocal(NULL), alocal(NULL) { if (narg < 4) error->all(FLERR,"Illegal compute improper/local command"); @@ -57,14 +58,16 @@ ComputeImproperLocal::ComputeImproperLocal(LAMMPS *lmp, int narg, char **arg) : } nmax = 0; + vlocal = NULL; + alocal = NULL; } /* ---------------------------------------------------------------------- */ ComputeImproperLocal::~ComputeImproperLocal() { - memory->destroy(vector); - memory->destroy(array); + memory->destroy(vlocal); + memory->destroy(alocal); } /* ---------------------------------------------------------------------- */ @@ -130,9 +133,9 @@ int ComputeImproperLocal::compute_impropers(int flag) if (flag) { if (nvalues == 1) { - if (cflag >= 0) cbuf = vector; + if (cflag >= 0) cbuf = vlocal; } else { - if (cflag >= 0 && array) cbuf = &array[0][cflag]; + if (cflag >= 0 && alocal) cbuf = &alocal[0][cflag]; else cbuf = NULL; } } @@ -228,18 +231,18 @@ int ComputeImproperLocal::compute_impropers(int flag) void ComputeImproperLocal::reallocate(int n) { - // grow vector or array and indices array + // grow vector_local or array_local while (nmax < n) nmax += DELTA; if (nvalues == 1) { - memory->destroy(vector); - memory->create(vector,nmax,"bond/local:vector"); - vector_local = vector; + memory->destroy(vlocal); + memory->create(vlocal,nmax,"improper/local:vector_local"); + vector_local = vlocal; } else { - memory->destroy(array); - memory->create(array,nmax,nvalues,"bond/local:array"); - array_local = array; + memory->destroy(alocal); + memory->create(alocal,nmax,nvalues,"improper/local:array_local"); + array_local = alocal; } } diff --git a/src/compute_improper_local.h b/src/compute_improper_local.h index 57f6b9dcf9..29530d603d 100644 --- a/src/compute_improper_local.h +++ b/src/compute_improper_local.h @@ -37,6 +37,8 @@ class ComputeImproperLocal : public Compute { int ncount; int nmax; + double *vlocal; + double **alocal; int compute_impropers(int); void reallocate(int); diff --git a/src/compute_pair_local.cpp b/src/compute_pair_local.cpp index 01a317f4a4..4595175503 100644 --- a/src/compute_pair_local.cpp +++ b/src/compute_pair_local.cpp @@ -37,7 +37,7 @@ enum{TYPE,RADIUS}; ComputePairLocal::ComputePairLocal(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg), - pstyle(NULL), pindex(NULL) + pstyle(NULL), pindex(NULL), vlocal(NULL), alocal(NULL) { if (narg < 4) error->all(FLERR,"Illegal compute pair/local command"); @@ -97,14 +97,16 @@ ComputePairLocal::ComputePairLocal(LAMMPS *lmp, int narg, char **arg) : if (pstyle[i] != DIST) singleflag = 1; nmax = 0; + vlocal = NULL; + alocal = NULL; } /* ---------------------------------------------------------------------- */ ComputePairLocal::~ComputePairLocal() { - memory->destroy(vector); - memory->destroy(array); + memory->destroy(vlocal); + memory->destroy(alocal); delete [] pstyle; delete [] pindex; } @@ -254,8 +256,8 @@ int ComputePairLocal::compute_pairs(int flag) eng = pair->single(i,j,itype,jtype,rsq,factor_coul,factor_lj,fpair); else eng = fpair = 0.0; - if (nvalues == 1) ptr = &vector[m]; - else ptr = array[m]; + if (nvalues == 1) ptr = &vlocal[m]; + else ptr = alocal[m]; for (n = 0; n < nvalues; n++) { switch (pstyle[n]) { @@ -295,18 +297,18 @@ int ComputePairLocal::compute_pairs(int flag) void ComputePairLocal::reallocate(int n) { - // grow vector or array and indices array + // grow vector_local or array_local while (nmax < n) nmax += DELTA; if (nvalues == 1) { - memory->destroy(vector); - memory->create(vector,nmax,"pair/local:vector"); - vector_local = vector; + memory->destroy(vlocal); + memory->create(vlocal,nmax,"pair/local:vector_local"); + vector_local = vlocal; } else { - memory->destroy(array); - memory->create(array,nmax,nvalues,"pair/local:array"); - array_local = array; + memory->destroy(alocal); + memory->create(alocal,nmax,nvalues,"pair/local:array_local"); + array_local = alocal; } } diff --git a/src/compute_pair_local.h b/src/compute_pair_local.h index 2e01cba3ac..df8e6a999b 100644 --- a/src/compute_pair_local.h +++ b/src/compute_pair_local.h @@ -41,6 +41,8 @@ class ComputePairLocal : public Compute { int singleflag; int nmax; + double *vlocal; + double **alocal; class NeighList *list; diff --git a/src/compute_property_atom.cpp b/src/compute_property_atom.cpp index 5d454fef84..d158d00816 100644 --- a/src/compute_property_atom.cpp +++ b/src/compute_property_atom.cpp @@ -361,8 +361,8 @@ ComputePropertyAtom::~ComputePropertyAtom() { delete [] pack_choice; delete [] index; - memory->destroy(vector); - memory->destroy(array); + memory->destroy(vector_atom); + memory->destroy(array_atom); } /* ---------------------------------------------------------------------- */ @@ -386,23 +386,21 @@ void ComputePropertyAtom::compute_peratom() if (atom->nmax > nmax) { nmax = atom->nmax; if (nvalues == 1) { - memory->destroy(vector); - memory->create(vector,nmax,"property/atom:vector"); - vector_atom = vector; + memory->destroy(vector_atom); + memory->create(vector_atom,nmax,"property/atom:vector"); } else { - memory->destroy(array); - memory->create(array,nmax,nvalues,"property/atom:array"); - array_atom = array; + memory->destroy(array_atom); + memory->create(array_atom,nmax,nvalues,"property/atom:array"); } } // fill vector or array with per-atom values if (nvalues == 1) { - buf = vector; + buf = vector_atom; (this->*pack_choice[0])(0); } else { - if (nmax) buf = &array[0][0]; + if (nmax) buf = &array_atom[0][0]; else buf = NULL; for (int n = 0; n < nvalues; n++) (this->*pack_choice[n])(n); diff --git a/src/compute_property_local.cpp b/src/compute_property_local.cpp index d39f8825e3..7814b6cb6b 100644 --- a/src/compute_property_local.cpp +++ b/src/compute_property_local.cpp @@ -35,7 +35,7 @@ enum{TYPE,RADIUS}; ComputePropertyLocal::ComputePropertyLocal(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg), - indices(NULL), pack_choice(NULL) + indices(NULL), pack_choice(NULL), vlocal(NULL), alocal(NULL) { if (narg < 4) error->all(FLERR,"Illegal compute property/local command"); @@ -254,6 +254,8 @@ ComputePropertyLocal::ComputePropertyLocal(LAMMPS *lmp, int narg, char **arg) : error->all(FLERR,"Compute property/local requires atom attribute radius"); nmax = 0; + vlocal = NULL; + alocal = NULL; } /* ---------------------------------------------------------------------- */ @@ -261,8 +263,8 @@ ComputePropertyLocal::ComputePropertyLocal(LAMMPS *lmp, int narg, char **arg) : ComputePropertyLocal::~ComputePropertyLocal() { delete [] pack_choice; - memory->destroy(vector); - memory->destroy(array); + memory->destroy(vlocal); + memory->destroy(alocal); memory->destroy(indices); } @@ -335,10 +337,10 @@ void ComputePropertyLocal::compute_local() // fill vector or array with local values if (nvalues == 1) { - buf = vector; + buf = vlocal; (this->*pack_choice[0])(0); } else { - if (array) buf = &array[0][0]; + if (alocal) buf = &alocal[0][0]; for (int n = 0; n < nvalues; n++) (this->*pack_choice[n])(n); } @@ -620,17 +622,18 @@ int ComputePropertyLocal::count_impropers(int flag) void ComputePropertyLocal::reallocate(int n) { - // grow vector or array and indices array + // grow vector_local or array_local, also indices while (nmax < n) nmax += DELTA; + if (nvalues == 1) { - memory->destroy(vector); - memory->create(vector,nmax,"property/local:vector"); - vector_local = vector; + memory->destroy(vlocal); + memory->create(vlocal,nmax,"property/local:vector_local"); + vector_local = vlocal; } else { - memory->destroy(array); - memory->create(array,nmax,nvalues,"property/local:array"); - array_local = array; + memory->destroy(alocal); + memory->create(alocal,nmax,nvalues,"property/local:array_local"); + array_local = alocal; } memory->destroy(indices); diff --git a/src/compute_property_local.h b/src/compute_property_local.h index fb2c5d2d18..fbc1ea3097 100644 --- a/src/compute_property_local.h +++ b/src/compute_property_local.h @@ -37,6 +37,8 @@ class ComputePropertyLocal : public Compute { int nvalues,kindflag,cutstyle; int nmax; + double *vlocal; + double **alocal; double *buf; class NeighList *list; diff --git a/src/compute_slice.cpp b/src/compute_slice.cpp index da077ddfc0..aa2c1067c4 100644 --- a/src/compute_slice.cpp +++ b/src/compute_slice.cpp @@ -148,10 +148,6 @@ ComputeSlice::ComputeSlice(LAMMPS *lmp, int narg, char **arg) : // for vector, set intensive/extensive to mirror input values // for array, set intensive if all input values are intensive, else extensive - vector = NULL; - array = NULL; - extlist = NULL; - if (nvalues == 1) { vector_flag = 1; size_vector = (nstop-nstart) / nskip; diff --git a/src/domain.cpp b/src/domain.cpp index c348167dfb..54183f6f2c 100644 --- a/src/domain.cpp +++ b/src/domain.cpp @@ -1499,22 +1499,78 @@ void Domain::image_flip(int m, int n, int p) /* ---------------------------------------------------------------------- return 1 if this proc owns atom with coords x, else return 0 x is returned remapped into periodic box + if image flag is passed, flag is updated via remap(x,image) + if image = NULL is passed, no update with remap(x) + if shrinkexceed, atom can be outside shrinkwrap boundaries + called from create_atoms() in library.cpp ------------------------------------------------------------------------- */ -int Domain::ownatom(double *x) +int Domain::ownatom(int id, double *x, imageint *image, int shrinkexceed) { double lamda[3]; - double *coord; + double *coord,*blo,*bhi,*slo,*shi; - remap(x); + if (image) remap(x,*image); + else remap(x); + if (triclinic) { x2lamda(x,lamda); coord = lamda; } else coord = x; - - if (coord[0] >= sublo[0] && coord[0] < subhi[0] && - coord[1] >= sublo[1] && coord[1] < subhi[1] && - coord[2] >= sublo[2] && coord[2] < subhi[2]) return 1; + + // box and subbox bounds for orthogonal vs triclinic + + if (triclinic == 0) { + blo = boxlo; + bhi = boxhi; + slo = sublo; + shi = subhi; + } else { + blo = boxlo_lamda; + bhi = boxhi_lamda; + slo = sublo_lamda; + shi = subhi_lamda; + } + + if (coord[0] >= slo[0] && coord[0] < shi[0] && + coord[1] >= slo[1] && coord[1] < shi[1] && + coord[2] >= slo[2] && coord[2] < shi[2]) return 1; + + // check if atom did not return 1 only b/c it was + // outside a shrink-wrapped boundary + + if (shrinkexceed) { + int outside = 0; + if (coord[0] < blo[0] && boundary[0][0] > 1) outside = 1; + if (coord[0] >= bhi[0] && boundary[0][1] > 1) outside = 1; + if (coord[1] < blo[1] && boundary[1][0] > 1) outside = 1; + if (coord[1] >= bhi[1] && boundary[1][1] > 1) outside = 1; + if (coord[2] < blo[2] && boundary[2][0] > 1) outside = 1; + if (coord[2] >= bhi[2] && boundary[2][1] > 1) outside = 1; + if (!outside) return 0; + + // newcoord = coords pushed back to be on shrink-wrapped boundary + // newcoord is a copy, so caller's x[] is not affected + + double newcoord[3]; + if (coord[0] < blo[0] && boundary[0][0] > 1) newcoord[0] = blo[0]; + else if (coord[0] >= bhi[0] && boundary[0][1] > 1) newcoord[0] = bhi[0]; + else newcoord[0] = coord[0]; + if (coord[1] < blo[1] && boundary[1][1] > 1) newcoord[1] = blo[1]; + else if (coord[1] >= bhi[1] && boundary[1][1] > 1) newcoord[1] = bhi[1]; + else newcoord[1] = coord[1]; + if (coord[2] < blo[2] && boundary[2][2] > 1) newcoord[2] = blo[2]; + else if (coord[2] >= bhi[2] && boundary[2][1] > 1) newcoord[2] = bhi[2]; + else newcoord[2] = coord[2]; + + // re-test for newcoord inside my sub-domain + // use <= test for upper-boundary since may have just put atom at boxhi + + if (newcoord[0] >= slo[0] && newcoord[0] <= shi[0] && + newcoord[1] >= slo[1] && newcoord[1] <= shi[1] && + newcoord[2] >= slo[2] && newcoord[2] <= shi[2]) return 1; + } + return 0; } diff --git a/src/domain.h b/src/domain.h index e2425f4585..ad55f051cf 100644 --- a/src/domain.h +++ b/src/domain.h @@ -121,7 +121,7 @@ class Domain : protected Pointers { void unmap(double *, imageint); void unmap(double *, imageint, double *); void image_flip(int, int, int); - int ownatom(double *); + int ownatom(int, double *, imageint *, int); void set_lattice(int, char **); void add_region(int, char **); @@ -149,7 +149,7 @@ class Domain : protected Pointers { // indicates a special neighbor is actually not in a bond, // but is a far-away image that should be treated as an unbonded neighbor // inline since called from neighbor build inner loop - // + inline int minimum_image_check(double dx, double dy, double dz) { if (xperiodic && fabs(dx) > xprd_half) return 1; if (yperiodic && fabs(dy) > yprd_half) return 1; diff --git a/src/fix_langevin.cpp b/src/fix_langevin.cpp index 22f4021b0f..542d2d6d6f 100644 --- a/src/fix_langevin.cpp +++ b/src/fix_langevin.cpp @@ -158,10 +158,9 @@ FixLangevin::FixLangevin(LAMMPS *lmp, int narg, char **arg) : tforce = NULL; maxatom1 = maxatom2 = 0; - // Setup atom-based array for franprev + // setup atom-based array for franprev // register with Atom class - // No need to set peratom_flag - // as this data is for internal use only + // no need to set peratom_flag, b/c data is for internal use only if (gjfflag) { nvalues = 3; diff --git a/src/integrate.cpp b/src/integrate.cpp index db71ef563c..4aeb6e18c9 100644 --- a/src/integrate.cpp +++ b/src/integrate.cpp @@ -59,9 +59,8 @@ void Integrate::init() // in case input script has reset the run or minimize style explicitly // e.g. invalid to have kokkos pair style with non-kokkos verlet // but OK to have kokkos verlet with non kokkos pair style (just warn) - // ditto for USER-CUDA package verlet with their pair, fix, etc // making these checks would require all the pair, fix, etc styles have - // cuda, kokkos, intel flags + // kokkos, intel flags } /* ---------------------------------------------------------------------- diff --git a/src/library.cpp b/src/library.cpp index 90c7cfa9d5..eac8e059ee 100644 --- a/src/library.cpp +++ b/src/library.cpp @@ -129,7 +129,7 @@ void lammps_open(int argc, char **argv, MPI_Comm communicator, void **ptr) } catch(LAMMPSException & e) { fprintf(stderr, "LAMMPS Exception: %s", e.message.c_str()); - *ptr = (void*) NULL; + *ptr = (void *) NULL; } #else LAMMPS *lmp = new LAMMPS(argc,argv,communicator); @@ -309,6 +309,23 @@ void lammps_free(void *ptr) customize by adding a function here and in library.h header file ------------------------------------------------------------------------- */ +/* ---------------------------------------------------------------------- + extract a LAMMPS setting as an integer + only use for settings that require return of an int + customize by adding names +------------------------------------------------------------------------- */ + +int lammps_extract_setting(void *ptr, char *name) +{ + LAMMPS *lmp = (LAMMPS *) ptr; + + if (strcmp(name,"bigint") == 0) return sizeof(bigint); + if (strcmp(name,"tagint") == 0) return sizeof(tagint); + if (strcmp(name,"imageint") == 0) return sizeof(imageint); + + return -1; +} + /* ---------------------------------------------------------------------- extract a pointer to an internal LAMMPS global entity name = desired quantity, e.g. dt or boxyhi or natoms @@ -325,12 +342,16 @@ void *lammps_extract_global(void *ptr, char *name) LAMMPS *lmp = (LAMMPS *) ptr; if (strcmp(name,"dt") == 0) return (void *) &lmp->update->dt; + if (strcmp(name,"boxlo") == 0) return (void *) lmp->domain->boxlo; + if (strcmp(name,"boxhi") == 0) return (void *) lmp->domain->boxhi; if (strcmp(name,"boxxlo") == 0) return (void *) &lmp->domain->boxlo[0]; if (strcmp(name,"boxxhi") == 0) return (void *) &lmp->domain->boxhi[0]; if (strcmp(name,"boxylo") == 0) return (void *) &lmp->domain->boxlo[1]; if (strcmp(name,"boxyhi") == 0) return (void *) &lmp->domain->boxhi[1]; if (strcmp(name,"boxzlo") == 0) return (void *) &lmp->domain->boxlo[2]; if (strcmp(name,"boxzhi") == 0) return (void *) &lmp->domain->boxhi[2]; + if (strcmp(name,"periodicity") == 0) return (void *) lmp->domain->periodicity; + if (strcmp(name,"xy") == 0) return (void *) &lmp->domain->xy; if (strcmp(name,"xz") == 0) return (void *) &lmp->domain->xz; if (strcmp(name,"yz") == 0) return (void *) &lmp->domain->yz; @@ -344,7 +365,12 @@ void *lammps_extract_global(void *ptr, char *name) if (strcmp(name,"nmax") == 0) return (void *) &lmp->atom->nmax; if (strcmp(name,"ntimestep") == 0) return (void *) &lmp->update->ntimestep; - // update atime can be referenced as a pointer + if (strcmp(name,"units") == 0) return (void *) lmp->update->unit_style; + if (strcmp(name,"triclinic") == 0) return (void *) &lmp->domain->triclinic; + + if (strcmp(name,"q_flag") == 0) return (void *) &lmp->atom->q_flag; + + // update->atime can be referenced as a pointer // thermo "timer" data cannot be, since it is computed on request // lammps_get_thermo() can access all thermo keywords by value @@ -354,6 +380,38 @@ void *lammps_extract_global(void *ptr, char *name) return NULL; } +/* ---------------------------------------------------------------------- + extract simulation box parameters + see domain.h for definition of these arguments + domain->init() call needed to set box_change +------------------------------------------------------------------------- */ + +void lammps_extract_box(void *ptr, double *boxlo, double *boxhi, + double *xy, double *yz, double *xz, + int *periodicity, int *box_change) +{ + LAMMPS *lmp = (LAMMPS *) ptr; + Domain *domain = lmp->domain; + domain->init(); + + boxlo[0] = domain->boxlo[0]; + boxlo[1] = domain->boxlo[1]; + boxlo[2] = domain->boxlo[2]; + boxhi[0] = domain->boxhi[0]; + boxhi[1] = domain->boxhi[1]; + boxhi[2] = domain->boxhi[2]; + + *xy = domain->xy; + *yz = domain->yz; + *xz = domain->xz; + + periodicity[0] = domain->periodicity[0]; + periodicity[1] = domain->periodicity[1]; + periodicity[2] = domain->periodicity[2]; + + *box_change = domain->box_change; +} + /* ---------------------------------------------------------------------- extract a pointer to an internal LAMMPS atom-based entity name = desired quantity, e.g. x or mass @@ -586,6 +644,35 @@ void *lammps_extract_variable(void *ptr, char *name, char *group) return NULL; } + +/* ---------------------------------------------------------------------- + reset simulation box parameters + see domain.h for definition of these arguments + assumes domain->set_intiial_box() has been invoked previously +------------------------------------------------------------------------- */ + +void lammps_reset_box(void *ptr, double *boxlo, double *boxhi, + double xy, double yz, double xz) +{ + LAMMPS *lmp = (LAMMPS *) ptr; + Domain *domain = lmp->domain; + + domain->boxlo[0] = boxlo[0]; + domain->boxlo[1] = boxlo[1]; + domain->boxlo[2] = boxlo[2]; + domain->boxhi[0] = boxhi[0]; + domain->boxhi[1] = boxhi[1]; + domain->boxhi[2] = boxhi[2]; + + domain->xy = xy; + domain->yz = yz; + domain->xz = xz; + + domain->set_global_box(); + lmp->comm->set_proc_grid(); + domain->set_local_box(); +} + /* ---------------------------------------------------------------------- set the value of a STRING variable to str return -1 if variable doesn't exist or not a STRING variable @@ -648,7 +735,7 @@ int lammps_get_natoms(void *ptr) name = desired quantity, e.g. x or charge type = 0 for integer values, 1 for double values count = # of per-atom values, e.g. 1 for type or charge, 3 for x or f - return atom-based values in data, ordered by count, then by atom ID + return atom-based values in 1d data, ordered by count, then by atom ID e.g. x[0][0],x[0][1],x[0][2],x[1][0],x[1][1],x[1][2],x[2][0],... data must be pre-allocated by caller to correct length ------------------------------------------------------------------------- */ @@ -744,7 +831,7 @@ void lammps_gather_atoms(void *ptr, char *name, name = desired quantity, e.g. x or charge type = 0 for integer values, 1 for double values count = # of per-atom values, e.g. 1 for type or charge, 3 for x or f - data = atom-based values in data, ordered by count, then by atom ID + data = atom-based values in 1d data, ordered by count, then by atom ID e.g. x[0][0],x[0][1],x[0][2],x[1][0],x[1][1],x[1][2],x[2][0],... ------------------------------------------------------------------------- */ @@ -823,16 +910,29 @@ void lammps_scatter_atoms(void *ptr, char *name, /* ---------------------------------------------------------------------- create N atoms and assign them to procs based on coords - id = atom IDs (optional, NULL if just use 1 to N) + id = atom IDs (optional, NULL will generate 1 to N) type = N-length vector of atom types (required) - x = 3N-length vector of atom coords (required) - v = 3N-length vector of atom velocities (optional, NULL if just 0.0) + x = 3N-length 1d vector of atom coords (required) + v = 3N-length 1d vector of atom velocities (optional, NULL if just 0.0) + image flags can be treated in two ways: + (a) image = vector of current image flags + each atom will be remapped into periodic box by domain->ownatom() + image flag will be incremented accordingly and stored with atom + (b) image = NULL + each atom will be remapped into periodic box by domain->ownatom() + image flag will be set to 0 by atom->avec->create_atom() + shrinkexceed = 1 allows atoms to be outside a shrinkwrapped boundary + passed to ownatom() which will assign them to boundary proc + important if atoms may be (slightly) outside non-periodic dim + e.g. due to restoring a snapshot from a previous run and previous box + id and image must be 32-bit integers x,v = ordered by xyz, then by atom e.g. x[0][0],x[0][1],x[0][2],x[1][0],x[1][1],x[1][2],x[2][0],... ------------------------------------------------------------------------- */ void lammps_create_atoms(void *ptr, int n, tagint *id, int *type, - double *x, double *v) + double *x, double *v, imageint *image, + int shrinkexceed) { LAMMPS *lmp = (LAMMPS *) ptr; @@ -857,14 +957,19 @@ void lammps_create_atoms(void *ptr, int n, tagint *id, int *type, Domain *domain = lmp->domain; int nlocal = atom->nlocal; - int nprev = nlocal; + bigint natoms_prev = atom->natoms; + int nlocal_prev = nlocal; double xdata[3]; for (int i = 0; i < n; i++) { xdata[0] = x[3*i]; xdata[1] = x[3*i+1]; xdata[2] = x[3*i+2]; - if (!domain->ownatom(xdata)) continue; + if (image) { + if (!domain->ownatom(id[i],xdata,&image[i],shrinkexceed)) continue; + } else { + if (!domain->ownatom(id[i],xdata,NULL,shrinkexceed)) continue; + } atom->avec->create_atom(type[i],xdata); if (id) atom->tag[nlocal] = id[i]; @@ -874,6 +979,7 @@ void lammps_create_atoms(void *ptr, int n, tagint *id, int *type, atom->v[nlocal][1] = v[3*i+1]; atom->v[nlocal][2] = v[3*i+2]; } + if (image) atom->image[nlocal] = image[i]; nlocal++; } @@ -885,8 +991,8 @@ void lammps_create_atoms(void *ptr, int n, tagint *id, int *type, // init per-atom fix/compute/variable values for created atoms - atom->data_fix_compute_variable(nprev,nlocal); - + atom->data_fix_compute_variable(nlocal_prev,nlocal); + // if global map exists, reset it // invoke map_init() b/c atom count has grown @@ -894,6 +1000,16 @@ void lammps_create_atoms(void *ptr, int n, tagint *id, int *type, lmp->atom->map_init(); lmp->atom->map_set(); } + + // warn if new natoms is not correct + + if (lmp->atom->natoms != natoms_prev + n) { + char str[128]; + sprintf(str,"Library warning in lammps_create_atoms, " + "invalid total atoms %ld %ld",lmp->atom->natoms,natoms_prev+n); + if (lmp->comm->me == 0) + lmp->error->warning(FLERR,str); + } } END_CAPTURE } diff --git a/src/library.h b/src/library.h index a05b7b4eb8..e4f9accaa6 100644 --- a/src/library.h +++ b/src/library.h @@ -34,25 +34,41 @@ void lammps_commands_list(void *, int, char **); void lammps_commands_string(void *, char *); void lammps_free(void *); +int lammps_extract_setting(void *, char *); void *lammps_extract_global(void *, char *); +void lammps_extract_box(void *, double *, double *, + double *, double *, double *, int *, int *); void *lammps_extract_atom(void *, char *); void *lammps_extract_compute(void *, char *, int, int); void *lammps_extract_fix(void *, char *, int, int, int, int); void *lammps_extract_variable(void *, char *, char *); +void lammps_reset_box(void *, double *, double *, double, double, double); int lammps_set_variable(void *, char *, char *); double lammps_get_thermo(void *, char *); int lammps_get_natoms(void *); void lammps_gather_atoms(void *, char *, int, int, void *); void lammps_scatter_atoms(void *, char *, int, int, void *); -void lammps_create_atoms(void *, int, int *, int *, double *, double *); + +// lammps_create_atoms() takes tagint and imageint as args +// ifdef insures they are compatible with rest of LAMMPS +// caller must match to how LAMMPS library is built + +#ifdef LAMMPS_BIGBIG +void lammps_create_atoms(void *, int, int64_t *, int *, + double *, double *, int64_t *, int); +#else +void lammps_create_atoms(void *, int, int *, int *, + double *, double *, int *, int); +#endif #ifdef LAMMPS_EXCEPTIONS int lammps_has_error(void *); int lammps_get_last_error_message(void *, char *, int); #endif +#undef LAMMPS #ifdef __cplusplus } #endif diff --git a/src/variable.h b/src/variable.h index cdcc607b18..042ce5b3fe 100644 --- a/src/variable.h +++ b/src/variable.h @@ -80,7 +80,7 @@ class Variable : protected Pointers { class Python *python; // ptr to embedded Python interpreter - struct Tree { // parse tree for atom-style or vector-style variables + struct Tree { // parse tree for atom-style or vector-style vars double value; // single scalar double *array; // per-atom or per-type list of doubles int *iarray; // per-atom list of ints