ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Merge branch 'develop' into delete-atoms-porosity-group

Browse Source
This commit is contained in:
Axel Kohlmeyer
2021-12-02 15:15:47 -05:00
parent 4a05628938 1a4511bb8d
commit 637c6bf28a
62 changed files with 687 additions and 473 deletions
Download Patch File Download Diff File Expand all files Collapse all files

4
src/DRUDE/compute_temp_drude.cpp
View File

@ -89,9 +89,7 @@ void ComputeTempDrude::dof_compute()
int dim = domain->dimension;
int *drudetype = fix_drude->drudetype;
fix_dof = 0;
for (int i = 0; i < modify->nfix; i++)
fix_dof += modify->fix[i]->dof(igroup);
adjust_dof_fix();
bigint dof_core_loc = 0, dof_drude_loc = 0;
for (int i = 0; i < nlocal; i++) {

1
src/DRUDE/compute_temp_drude.h
View File

@ -35,7 +35,6 @@ class ComputeTempDrude : public Compute {
int modify_param(int, char **);
private:
int fix_dof;
class FixDrude *fix_drude;
char *id_temp;
class Compute *temperature;

1
src/EXTRA-COMPUTE/compute_temp_rotate.h
View File

@ -43,7 +43,6 @@ class ComputeTempRotate : public Compute {
double memory_usage();
private:
int fix_dof;
double tfactor, masstotal;
double **vbiasall; // stored velocity bias for all atoms
int maxbias; // size of vbiasall array

3
src/MACHDYN/atom_vec_smd.cpp
View File

@ -115,6 +115,7 @@ void AtomVecSMD::grow_pointers()
vfrac = atom->vfrac;
rmass = atom->rmass;
x0 = atom->x0;
x = atom->x;
radius = atom->radius;
contact_radius = atom->contact_radius;
molecule = atom->molecule;
@ -129,13 +130,11 @@ void AtomVecSMD::grow_pointers()
/* ----------------------------------------------------------------------
clear extra forces starting at atom N
nbytes = # of bytes to clear for a per-atom vector
NOTE: does f need to be re-cleared?
------------------------------------------------------------------------- */
void AtomVecSMD::force_clear(int n, size_t nbytes)
{
memset(&desph[n],0,nbytes);
memset(&f[n][0],0,3*nbytes);
}
/* ----------------------------------------------------------------------

2
src/MACHDYN/compute_smd_triangle_vertices.cpp
View File

@ -53,7 +53,7 @@ ComputeSMDTriangleVertices::ComputeSMDTriangleVertices(LAMMPS *lmp, int narg, ch
/* ---------------------------------------------------------------------- */
ComputeSMDTriangleVertices::~ComputeSMDTriangleVertices() {
memory->sfree(outputVector);
memory->destroy(outputVector);
}
/* ---------------------------------------------------------------------- */

3
src/MACHDYN/pair_smd_ulsph.cpp
View File

@ -85,7 +85,8 @@ PairULSPH::PairULSPH(LAMMPS *lmp) :
PairULSPH::~PairULSPH() {
if (allocated) {
//printf("... deallocating\n");
memory->destroy(setflag);
memory->destroy(cutsq);
memory->destroy(Q1);
memory->destroy(rho0);
memory->destroy(eos);

17
src/compute_bond_local.cpp
View File

@ -34,7 +34,7 @@ using namespace LAMMPS_NS;
#define DELTA 10000
#define EPSILON 1.0e-12
enum{DIST,VELVIB,OMEGA,ENGTRANS,ENGVIB,ENGROT,ENGPOT,FORCE,FX,FY,FZ,VARIABLE};
enum{DIST,DX,DY,DZ,VELVIB,OMEGA,ENGTRANS,ENGVIB,ENGROT,ENGPOT,FORCE,FX,FY,FZ,VARIABLE};
/* ---------------------------------------------------------------------- */
@ -63,6 +63,9 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) :
int iarg;
for (iarg = 3; iarg < narg; iarg++) {
if (strcmp(arg[iarg],"dist") == 0) bstyle[nvalues++] = DIST;
else if (strcmp(arg[iarg],"dx") == 0) bstyle[nvalues++] = DX;
else if (strcmp(arg[iarg],"dy") == 0) bstyle[nvalues++] = DY;
else if (strcmp(arg[iarg],"dz") == 0) bstyle[nvalues++] = DZ;
else if (strcmp(arg[iarg],"engpot") == 0) bstyle[nvalues++] = ENGPOT;
else if (strcmp(arg[iarg],"force") == 0) bstyle[nvalues++] = FORCE;
else if (strcmp(arg[iarg],"fx") == 0) bstyle[nvalues++] = FX;
@ -384,11 +387,23 @@ int ComputeBondLocal::compute_bonds(int flag)
if (dstr) input->variable->internal_set(dvar,sqrt(rsq));
}
// to make sure dx, dy and dz are always from the lower to the higher id
double directionCorrection = tag[atom1] > tag[atom2] ? -1.0 : 1.0;
for (int n = 0; n < nvalues; n++) {
switch (bstyle[n]) {
case DIST:
ptr[n] = sqrt(rsq);
break;
case DX:
ptr[n] = dx*directionCorrection;
break;
case DY:
ptr[n] = dy*directionCorrection;
break;
case DZ:
ptr[n] = dz*directionCorrection;
break;
case ENGPOT:
ptr[n] = engpot;
break;

16
src/compute_pair_local.cpp
View File

@ -31,7 +31,7 @@ using namespace LAMMPS_NS;
#define DELTA 10000
enum{DIST,ENG,FORCE,FX,FY,FZ,PN};
enum{DIST,ENG,FORCE,FX,FY,FZ,PN,DX,DY,DZ};
enum{TYPE,RADIUS};
/* ---------------------------------------------------------------------- */
@ -56,6 +56,9 @@ ComputePairLocal::ComputePairLocal(LAMMPS *lmp, int narg, char **arg) :
else if (strcmp(arg[iarg],"fx") == 0) pstyle[nvalues++] = FX;
else if (strcmp(arg[iarg],"fy") == 0) pstyle[nvalues++] = FY;
else if (strcmp(arg[iarg],"fz") == 0) pstyle[nvalues++] = FZ;
else if (strcmp(arg[iarg],"dx") == 0) pstyle[nvalues++] = DX;
else if (strcmp(arg[iarg],"dy") == 0) pstyle[nvalues++] = DY;
else if (strcmp(arg[iarg],"dz") == 0) pstyle[nvalues++] = DZ;
else if (arg[iarg][0] == 'p') {
int n = atoi(&arg[iarg][1]);
if (n <= 0) error->all(FLERR,
@ -92,7 +95,7 @@ ComputePairLocal::ComputePairLocal(LAMMPS *lmp, int narg, char **arg) :
singleflag = 0;
for (int i = 0; i < nvalues; i++)
if (pstyle[i] != DIST) singleflag = 1;
if (pstyle[i] != DIST && pstyle[i] != DX && pstyle[i] != DY && pstyle[i] != DZ) singleflag = 1;
if (nvalues == 1) size_local_cols = 0;
else size_local_cols = nvalues;
@ -264,11 +267,20 @@ int ComputePairLocal::compute_pairs(int flag)
if (nvalues == 1) ptr = &vlocal[m];
else ptr = alocal[m];
// to make sure dx, dy and dz are always from the lower to the higher id
double directionCorrection = itag > jtag ? -1.0 : 1.0;
for (n = 0; n < nvalues; n++) {
switch (pstyle[n]) {
case DIST:
ptr[n] = sqrt(rsq);
break;
case DX:
ptr[n] = delx*directionCorrection;
case DY:
ptr[n] = dely*directionCorrection;
case DZ:
ptr[n] = delz*directionCorrection;
case ENG:
ptr[n] = eng;
break;

10
src/fix_deposit.cpp
View File

@ -569,10 +569,12 @@ void FixDeposit::pre_exchange()
// coord is new position of geometric center of mol, not COM
// FixShake::set_molecule stores shake info for molecule
if (rigidflag)
fixrigid->set_molecule(nlocalprev,maxtag_all,imol,coord,vnew,quat);
else if (shakeflag)
fixshake->set_molecule(nlocalprev,maxtag_all,imol,coord,vnew,quat);
if (mode == MOLECULE) {
if (rigidflag)
fixrigid->set_molecule(nlocalprev,maxtag_all,imol,coord,vnew,quat);
else if (shakeflag)
fixshake->set_molecule(nlocalprev,maxtag_all,imol,coord,vnew,quat);
}
success = 1;
break;

24
src/pair.cpp
View File

@ -41,6 +41,7 @@ using namespace LAMMPS_NS;
using namespace MathConst;
enum{NONE,RLINEAR,RSQ,BMP};
static const std::string mixing_rule_names[Pair::SIXTHPOWER+1] = {"geometric", "arithmetic", "sixthpower" };
// allocate space for static class instance variable and initialize it
@ -217,11 +218,9 @@ void Pair::init()
if (tail_flag && domain->nonperiodic && comm->me == 0)
error->warning(FLERR,"Using pair tail corrections with non-periodic system");
if (!compute_flag && tail_flag && comm->me == 0)
error->warning(FLERR,"Using pair tail corrections with "
"pair_modify compute no");
error->warning(FLERR,"Using pair tail corrections with pair_modify compute no");
if (!compute_flag && offset_flag && comm->me == 0)
error->warning(FLERR,"Using pair potential shift with "
"pair_modify compute no");
error->warning(FLERR,"Using pair potential shift with pair_modify compute no");
// for manybody potentials
// check if bonded exclusions could invalidate the neighbor list
@ -259,13 +258,18 @@ void Pair::init()
etail = ptail = 0.0;
mixed_flag = 1;
double cut;
int mixed_count = 0;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
if ((i != j) && setflag[i][j]) mixed_flag = 0;
did_mix = false;
cut = init_one(i,j);
cutsq[i][j] = cutsq[j][i] = cut*cut;
cutforce = MAX(cutforce,cut);
if (i != j) {
if (setflag[i][j]) mixed_flag = 0;
if (did_mix) ++mixed_count;
}
if (tail_flag) {
etail += etail_ij;
ptail += ptail_ij;
@ -275,6 +279,12 @@ void Pair::init()
}
}
}
if (!manybody_flag && (comm->me == 0)) {
const int num_mixed_pairs = atom->ntypes * (atom->ntypes - 1) / 2;
utils::logmesg(lmp," generated {} of {} mixed pair_coeff terms from {} mixing rule\n",
mixed_count, num_mixed_pairs, mixing_rule_names[mix_flag]);
}
}
/* ----------------------------------------------------------------------
@ -681,6 +691,7 @@ void Pair::free_disp_tables()
double Pair::mix_energy(double eps1, double eps2, double sig1, double sig2)
{
did_mix = true;
if (mix_flag == GEOMETRIC)
return sqrt(eps1*eps2);
else if (mix_flag == ARITHMETIC)
@ -688,7 +699,8 @@ double Pair::mix_energy(double eps1, double eps2, double sig1, double sig2)
else if (mix_flag == SIXTHPOWER)
return (2.0 * sqrt(eps1*eps2) *
pow(sig1,3.0) * pow(sig2,3.0) / (pow(sig1,6.0) + pow(sig2,6.0)));
else return 0.0;
else did_mix = false;
return 0.0;
}
/* ----------------------------------------------------------------------

1
src/pair.h
View File

@ -110,6 +110,7 @@ class Pair : protected Pointers {
// public so external driver can check
int compute_flag; // 0 if skip compute()
int mixed_flag; // 1 if all itype != jtype coeffs are from mixing
bool did_mix; // set to true by mix_energy() to indicate that mixing was performed
enum { GEOMETRIC, ARITHMETIC, SIXTHPOWER }; // mixing options

7
src/pair_hybrid.cpp
View File

@ -706,11 +706,10 @@ double PairHybrid::init_one(int i, int j)
for (int k = 0; k < nmap[i][j]; k++) {
map[j][i][k] = map[i][j][k];
double cut = styles[map[i][j][k]]->init_one(i,j);
styles[map[i][j][k]]->cutsq[i][j] =
styles[map[i][j][k]]->cutsq[j][i] = cut*cut;
if (styles[map[i][j][k]]->did_mix) did_mix = true;
styles[map[i][j][k]]->cutsq[i][j] = styles[map[i][j][k]]->cutsq[j][i] = cut*cut;
if (styles[map[i][j][k]]->ghostneigh)
cutghost[i][j] = cutghost[j][i] =
MAX(cutghost[i][j],styles[map[i][j][k]]->cutghost[i][j]);
cutghost[i][j] = cutghost[j][i] = MAX(cutghost[i][j],styles[map[i][j][k]]->cutghost[i][j]);
if (tail_flag) {
etail_ij += styles[map[i][j][k]]->etail_ij;
ptail_ij += styles[map[i][j][k]]->ptail_ij;

24
src/platform.cpp
View File

@ -446,11 +446,11 @@ int platform::putenv(const std::string &vardef)
auto found = vardef.find_first_of('=');
#ifdef _WIN32
// must assign a value to variable with _putenv()
// must assign a value to variable with _putenv_s()
if (found == std::string::npos)
return _putenv(utils::strdup(vardef + "=1"));
return _putenv_s(vardef.c_str(), "1");
else
return _putenv(utils::strdup(vardef));
return _putenv_s(vardef.substr(0, found).c_str(), vardef.substr(found+1).c_str());
#else
if (found == std::string::npos)
return setenv(vardef.c_str(), "", 1);
@ -460,6 +460,24 @@ int platform::putenv(const std::string &vardef)
return -1;
}
/* ----------------------------------------------------------------------
unset environment variable
------------------------------------------------------------------------- */
int platform::unsetenv(const std::string &variable)
{
if (variable.size() == 0) return -1;
#ifdef _WIN32
// emulate POSIX semantics by returning -1 on trying to unset non-existing variable
const char *ptr = getenv(variable.c_str());
if (!ptr) return -1;
// empty _putenv_s() definition deletes variable
return _putenv_s(variable.c_str(),"");
#else
return ::unsetenv(variable.c_str());
#endif
}
/* ----------------------------------------------------------------------
split a "path" environment variable into a list
------------------------------------------------------------------------- */

7
src/platform.h
View File

@ -125,6 +125,13 @@ namespace platform {
int putenv(const std::string &vardef);
/*! Delete variable from the environment
*
* \param variable variable name
* \return -1 if failure otherwise 0 */
int unsetenv(const std::string &variable);
/*! Get list of entries in a path environment variable
*
* This provides a list of strings of the entries in an environment

2
src/verlet.cpp
View File

@ -54,7 +54,7 @@ void Verlet::init()
bool do_time_integrate = false;
for (const auto &fix : modify->get_fix_list())
if (fix->time_integrate) do_time_integrate;
if (fix->time_integrate) do_time_integrate = true;
if (!do_time_integrate && (comm->me == 0))
error->warning(FLERR,"No fixes with time integration, atoms won't move");