added functionity to lib interface

This commit is contained in:
Steve Plimpton 2016-12-13 16:22:17 -07:00
parent fb3f597f41
commit ae5764beac
28 changed files with 386 additions and 160 deletions

View File

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

View File

@ -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; k<EWALD_NFUNCS; ++k) energy += c[k]*sum[k]-energy_self[k];
if (slabflag) compute_slabcorr();
}
/* ---------------------------------------------------------------------- */
@ -1480,10 +1480,7 @@ double EwaldDisp::f(double x, double Rc, bigint natoms, double vol, double b2)
double a = Rc*x;
double f = 0.0;
if (function[1] || function[2]) { // LJ
f = (4.0*MY_PI*b2*powint(x,4)/vol/sqrt((double)natoms)*erfc(a) *
(6.0*powint(a,-5) + 6.0*powint(a,-3) + 3.0/a + a) - accuracy);
} else { // dipole
if (function[3]) { // dipole
double rg2 = a*a;
double rg4 = rg2*rg2;
double rg6 = rg4*rg2;
@ -1492,7 +1489,10 @@ double EwaldDisp::f(double x, double Rc, bigint natoms, double vol, double b2)
f = (b2/(sqrt(vol*powint(x,4)*powint(Rc,9)*natoms)) *
sqrt(13.0/6.0*Cc*Cc + 2.0/15.0*Dc*Dc - 13.0/15.0*Cc*Dc) *
exp(-rg2)) - accuracy;
}
} else if (function[1] || function[2]) { // LJ
f = (4.0*MY_PI*b2*powint(x,4)/vol/sqrt((double)natoms)*erfc(a) *
(6.0*powint(a,-5) + 6.0*powint(a,-3) + 3.0/a + a) - accuracy);
}
return f;
}

View File

@ -241,6 +241,7 @@ void PRD::command(int narg, char **arg)
update->minimize->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");

View File

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

View File

@ -41,6 +41,8 @@ class ComputeRigidLocal : public Compute {
class FixRigidSmall *fixrigid;
int nmax;
double *vlocal;
double **alocal;
int compute_rigid(int);
void reallocate(int);

View File

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

View File

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

View File

@ -37,6 +37,8 @@ class ComputeAngleLocal : public Compute {
int ncount;
int nmax;
double *vlocal;
double **alocal;
int compute_angles(int);
void reallocate(int);

View File

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

View File

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

View File

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

View File

@ -37,6 +37,8 @@ class ComputeDihedralLocal : public Compute {
int ncount;
int nmax;
double *vlocal;
double **alocal;
int compute_dihedrals(int);
void reallocate(int);

View File

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

View File

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

View File

@ -37,6 +37,8 @@ class ComputeImproperLocal : public Compute {
int ncount;
int nmax;
double *vlocal;
double **alocal;
int compute_impropers(int);
void reallocate(int);

View File

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

View File

@ -41,6 +41,8 @@ class ComputePairLocal : public Compute {
int singleflag;
int nmax;
double *vlocal;
double **alocal;
class NeighList *list;

View File

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

View File

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

View File

@ -37,6 +37,8 @@ class ComputePropertyLocal : public Compute {
int nvalues,kindflag,cutstyle;
int nmax;
double *vlocal;
double **alocal;
double *buf;
class NeighList *list;

View File

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

View File

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

View File

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

View File

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

View File

@ -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
}
/* ----------------------------------------------------------------------

View File

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

View File

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

View File

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