lijiext
/
lammps
mirror of https://github.com/lammps/lammps.git
1
0
Fork 1
Code Issues Packages Projects Releases Wiki Activity

Merge branch 'develop' into mliappy_unified

This commit is contained in:
Axel Kohlmeyer 2022-09-24 03:07:31 -04:00
parent 49a3be17cb e230a3bb5b
commit c581fca39a
No known key found for this signature in database
GPG Key ID: D9B44E93BF0C375A
18 changed files with 212 additions and 284 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

4
cmake/Modules/Packages/MDI.cmake
View File

@ -8,8 +8,8 @@ option(DOWNLOAD_MDI "Download and compile the MDI library instead of using an al
if(DOWNLOAD_MDI)
message(STATUS "MDI download requested - we will build our own")
set(MDI_URL "https://github.com/MolSSI-MDI/MDI_Library/archive/v1.4.11.tar.gz" CACHE STRING "URL for MDI tarball")
set(MDI_MD5 "3791fe5081405c14aac07d4687f1cc58" CACHE STRING "MD5 checksum for MDI tarball")
set(MDI_URL "https://github.com/MolSSI-MDI/MDI_Library/archive/v1.4.12.tar.gz" CACHE STRING "URL for MDI tarball")
set(MDI_MD5 "7a222353ae8e03961d5365e6cd48baee" CACHE STRING "MD5 checksum for MDI tarball")
mark_as_advanced(MDI_URL)
mark_as_advanced(MDI_MD5)
enable_language(C)

26
doc/src/pair_gayberne.rst
View File

@ -44,14 +44,15 @@ ellipsoidal and spherical particle via the formulas
U_r = & 4 \epsilon ( \varrho^{12} - \varrho^6) \\
\varrho = & \frac{\sigma}{ h_{12} + \gamma \sigma}
where A1 and A2 are the transformation matrices from the simulation box
frame to the body frame and :math:`r_{12}` is the center to center
vector between the particles. :math:`U_r` controls the shifted distance
dependent interaction based on the distance of closest approach of the
two particles (:math:`h_{12}`) and the user-specified shift parameter
gamma. When both particles are spherical, the formula reduces to the
usual Lennard-Jones interaction (see details below for when Gay-Berne
treats a particle as "spherical").
where :math:`\mathbf{A}_1` and :math:`\mathbf{A}_2` are the
transformation matrices from the simulation box frame to the body frame
and :math:`r_{12}` is the center to center vector between the particles.
:math:`U_r` controls the shifted distance dependent interaction based on
the distance of closest approach of the two particles (:math:`h_{12}`)
and the user-specified shift parameter :math:`\gamma`. When both
particles are spherical, the formula reduces to the usual Lennard-Jones
interaction (see details below for when Gay-Berne treats a particle as
"spherical").
For large uniform molecules it has been shown that the energy
parameters are approximately representable in terms of local contact
@ -74,8 +75,9 @@ listed below and in `this supplementary document <PDF/pair_gayberne_extra.pdf>`_
Use of this pair style requires the NVE, NVT, or NPT fixes with the
*asphere* extension (e.g. :doc:`fix nve/asphere <fix_nve_asphere>`) in
order to integrate particle rotation. Additionally, :doc:`atom_style ellipsoid <atom_style>` should be used since it defines the
rotational state and the size and shape of each ellipsoidal particle.
order to integrate particle rotation. Additionally, :doc:`atom_style
ellipsoid <atom_style>` should be used since it defines the rotational
state and the size and shape of each ellipsoidal particle.
The following coefficients must be defined for each pair of atoms
types via the :doc:`pair_coeff <pair_coeff>` command as in the examples
@ -129,7 +131,7 @@ that type. e.g. in a "pair_coeff I J" command.
.. note::
If the :math:`\epsilon` a = b = c for an atom type, and if the shape
If the :math:`\epsilon_{a}` = :math:`\epsilon_{b}` = :math:`\epsilon_{c}` for an atom type, and if the shape
of the particle itself is spherical, meaning its 3 shape parameters
are all the same, then the particle is treated as an LJ sphere by the
Gay-Berne potential. This is significant because if two LJ spheres
@ -137,7 +139,7 @@ that type. e.g. in a "pair_coeff I J" command.
their interaction energy/force using the specified epsilon and sigma
as the standard LJ parameters. This is much cheaper to compute than
the full Gay-Berne formula. To treat the particle as a LJ sphere
with sigma = D, you should normally set :math:`\epsilon` a = b = c =
with sigma = D, you should normally set :math:`\epsilon_{a}` = :math:`\epsilon_{b}` = :math:`\epsilon_{c}` =
1.0, set the pair_coeff :math:`\sigma = D`, and also set the 3 shape
parameters for the particle to D. The one exception is that if the 3
shape parameters are set to 0.0, which is a valid way in LAMMPS to

3
lib/mdi/Install.py
View File

@ -32,12 +32,13 @@ make lib-mdi args="-m mpi" # build MDI lib with same settings as in the mpi Make
# settings
version = "1.4.11"
version = "1.4.12"
url = "https://github.com/MolSSI-MDI/MDI_Library/archive/v%s.tar.gz" % version
# known checksums for different MDI versions. used to validate the download.
checksums = { \
'1.4.11' : '3791fe5081405c14aac07d4687f1cc58', \
'1.4.12' : '7a222353ae8e03961d5365e6cd48baee', \
}
# print error message or help

5
src/COLVARS/group_ndx.cpp
View File

@ -47,7 +47,7 @@ static int cmptagint(const void *p1, const void *p2)
void Group2Ndx::command(int narg, char **arg)
{
FILE *fp;
FILE *fp = nullptr;
if (narg < 1) error->all(FLERR,"Illegal group2ndx command");
@ -57,8 +57,7 @@ void Group2Ndx::command(int narg, char **arg)
if (comm->me == 0) {
fp = fopen(arg[0], "w");
if (fp == nullptr)
error->one(FLERR,"Cannot open index file for writing: {}",
utils::getsyserror());
error->one(FLERR,"Cannot open index file for writing: {}", utils::getsyserror());
utils::logmesg(lmp,"Writing groups to index file {}:\n",arg[0]);
}

1
src/EFF/fix_nvt_sllod_eff.cpp
View File

@ -51,6 +51,7 @@ FixNVTSllodEff::FixNVTSllodEff(LAMMPS *lmp, int narg, char **arg) :
modify->add_compute(fmt::format("{} {} tmp/deform/eff",
id_temp,group->names[igroup]));
tcomputeflag = 1;
nondeformbias = 0;
}
/* ---------------------------------------------------------------------- */

2
src/KSPACE/ewald_disp.cpp
View File

@ -398,7 +398,7 @@ void EwaldDisp::reallocate()
h[0] = unit[0]*ix;
h[1] = unit[5]*ix+unit[1]*iy;
h[2] = unit[4]*ix+unit[3]*iy+unit[2]*iz;
if ((*(flag++) = h[0]*h[0]+h[1]*h[1]+h[2]*h[2]<=gsqmx)) ++nkvec;
if ((*(flag++) = (h[0]*h[0]+h[1]*h[1]+h[2]*h[2]<=gsqmx))) ++nkvec;
}
if (nkvec>nkvec_max) {

88
src/MC/fix_charge_regulation.cpp
View File

@ -140,6 +140,8 @@ FixChargeRegulation::FixChargeRegulation(LAMMPS *lmp, int narg, char **arg) :
nsalt_successes = 0;
}
/* ---------------------------------------------------------------------- */
FixChargeRegulation::~FixChargeRegulation() {
memory->destroy(ptype_ID);
@ -153,14 +155,23 @@ FixChargeRegulation::~FixChargeRegulation() {
int igroupall = group->find("all");
neighbor->exclusion_group_group_delete(exclusion_group, igroupall);
}
if (groupstrings) {
for (int i = 0; i < ngroups; ++i) delete[] groupstrings[i];
memory->destroy(groupstrings);
}
}
/* ---------------------------------------------------------------------- */
int FixChargeRegulation::setmask() {
int mask = 0;
mask |= PRE_EXCHANGE;
return mask;
}
/* ---------------------------------------------------------------------- */
void FixChargeRegulation::init() {
triclinic = domain->triclinic;
@ -246,6 +257,8 @@ void FixChargeRegulation::init() {
}
}
/* ---------------------------------------------------------------------- */
void FixChargeRegulation::pre_exchange() {
if (next_reneighbor != update->ntimestep) return;
@ -372,16 +385,15 @@ void FixChargeRegulation::pre_exchange() {
next_reneighbor = update->ntimestep + nevery;
}
/* ---------------------------------------------------------------------- */
void FixChargeRegulation::forward_acid() {
double energy_before = energy_stored;
double factor;
double dummyp[3];
double pos[3];
pos[0] = 0;
pos[1] = 0;
pos[2] = 0; // acid/base particle position
double pos_all[3];
double dummyp[3] = {0.0, 0.0, 0.0};
double pos[3] = {0.0, 0.0, 0.0}; // acid/base particle position
double pos_all[3] = {0.0, 0.0, 0.0};
int m1 = -1, m2 = -1;
m1 = get_random_particle(acid_type, 0, 0, dummyp);
@ -432,17 +444,16 @@ void FixChargeRegulation::forward_acid() {
}
}
/* ---------------------------------------------------------------------- */
void FixChargeRegulation::backward_acid() {
double energy_before = energy_stored;
double factor;
int mask_tmp;
double dummyp[3];
double pos[3];
pos[0] = 0;
pos[1] = 0;
pos[2] = 0; // acid/base particle position
double pos_all[3];
double dummyp[3] = {0.0, 0.0, 0.0};
double pos[3] = {0.0, 0.0, 0.0}; // acid/base particle position
double pos_all[3] = {0.0, 0.0, 0.0};
int m1 = -1, m2 = -1;
m1 = get_random_particle(acid_type, -1, 0, dummyp);
@ -510,16 +521,15 @@ void FixChargeRegulation::backward_acid() {
}
}
/* ---------------------------------------------------------------------- */
void FixChargeRegulation::forward_base() {
double energy_before = energy_stored;
double factor;
double dummyp[3];
double pos[3];
pos[0] = 0;
pos[1] = 0;
pos[2] = 0; // acid/base particle position
double pos_all[3];
double dummyp[3] = {0.0, 0.0, 0.0};
double pos[3] = {0.0, 0.0, 0.0}; // acid/base particle position
double pos_all[3] = {0.0, 0.0, 0.0};
int m1 = -1, m2 = -1;
m1 = get_random_particle(base_type, 0, 0, dummyp);
@ -570,17 +580,16 @@ void FixChargeRegulation::forward_base() {
}
}
/* ---------------------------------------------------------------------- */
void FixChargeRegulation::backward_base() {
double energy_before = energy_stored;
double factor;
double dummyp[3];
double dummyp[3] = {0.0, 0.0, 0.0};
double pos[3] = {0.0, 0.0, 0.0}; // acid/base particle position
double pos_all[3] = {0.0, 0.0, 0.0};
int mask_tmp;
double pos[3];
pos[0] = 0;
pos[1] = 0;
pos[2] = 0; // acid/base particle position
double pos_all[3];
int m1 = -1, m2 = -1;
m1 = get_random_particle(base_type, 1, 0, dummyp);
@ -647,11 +656,13 @@ void FixChargeRegulation::backward_base() {
}
}
/* ---------------------------------------------------------------------- */
void FixChargeRegulation::forward_ions() {
double energy_before = energy_stored;
double factor;
double dummyp[3];
double dummyp[3] = {0.0, 0.0, 0.0};
int m1 = -1, m2 = -1;
factor = volume_rx * volume_rx * c10pI_plus * c10pI_minus /
((1 + ncation) * (1 + nanion));
@ -682,13 +693,14 @@ void FixChargeRegulation::forward_ions() {
}
}
/* ---------------------------------------------------------------------- */
void FixChargeRegulation::backward_ions() {
double energy_before = energy_stored;
double factor;
int mask1_tmp = 0, mask2_tmp = 0;
double dummyp[3];
double dummyp[3] = {0.0, 0.0, 0.0};
int m1 = -1, m2 = -1;
m1 = get_random_particle(cation_type, +1, 0, dummyp);
@ -764,11 +776,13 @@ void FixChargeRegulation::backward_ions() {
}
}
/* ---------------------------------------------------------------------- */
void FixChargeRegulation::forward_ions_multival() {
double energy_before = energy_stored;
double factor = 1;
double dummyp[3];
double dummyp[3] = {0.0, 0.0, 0.0};
// particle ID array for all ions to be inserted
auto mm = std::unique_ptr<int[]>(new int[salt_charge_ratio + 1]);
@ -822,11 +836,13 @@ void FixChargeRegulation::forward_ions_multival() {
}
}
/* ---------------------------------------------------------------------- */
void FixChargeRegulation::backward_ions_multival() {
double energy_before = energy_stored;
double factor = 1;
double dummyp[3]; // dummy particle
double dummyp[3] = {0.0, 0.0, 0.0}; // dummy particle
// particle ID array for all deleted ions
auto mm = std::unique_ptr<int[]>(new int[salt_charge_ratio + 1]);
// charge array for all deleted ions
@ -943,6 +959,8 @@ void FixChargeRegulation::backward_ions_multival() {
}
}
/* ---------------------------------------------------------------------- */
int FixChargeRegulation::insert_particle(int ptype, double charge, double rd, double *target) {
// insert a particle of type (ptype) with charge (charge) within distance (rd) of (target)
@ -1013,6 +1031,8 @@ int FixChargeRegulation::insert_particle(int ptype, double charge, double rd, do
return m;
}
/* ---------------------------------------------------------------------- */
int FixChargeRegulation::get_random_particle(int ptype, double charge, double rd, double *target) {
// returns a randomly chosen particle of type (ptype) with charge (charge)
@ -1077,6 +1097,8 @@ int FixChargeRegulation::get_random_particle(int ptype, double charge, double rd
return -1;
}
/* ---------------------------------------------------------------------- */
double FixChargeRegulation::energy_full() {
if (triclinic) domain->x2lamda(atom->nlocal);
domain->pbc();
@ -1136,6 +1158,8 @@ double FixChargeRegulation::energy_full() {
return total_energy;
}
/* ---------------------------------------------------------------------- */
int FixChargeRegulation::particle_number_xrd(int ptype, double charge, double rd, double *target) {
int count = 0;
@ -1165,6 +1189,8 @@ int FixChargeRegulation::particle_number_xrd(int ptype, double charge, double rd
return count_sum;
}
/* ---------------------------------------------------------------------- */
int FixChargeRegulation::particle_number(int ptype, double charge) {
int count = 0;
@ -1177,6 +1203,8 @@ int FixChargeRegulation::particle_number(int ptype, double charge) {
return count_sum;
}
/* ---------------------------------------------------------------------- */
double FixChargeRegulation::compute_vector(int n) {
if (n == 0) {
return nacid_attempts + nbase_attempts + nsalt_attempts;
@ -1253,6 +1281,8 @@ void FixChargeRegulation::restart(char *buf)
error->all(FLERR,"Must not reset timestep when restarting fix gcmc");
}
/* ---------------------------------------------------------------------- */
void FixChargeRegulation::setThermoTemperaturePointer() {
int ifix = -1;
ifix = modify->find_fix(idftemp);
@ -1266,6 +1296,8 @@ void FixChargeRegulation::setThermoTemperaturePointer() {
}
/* ---------------------------------------------------------------------- */
void FixChargeRegulation::assign_tags() {
// Assign tags to ions with zero tags
if (atom->tag_enable) {

302
src/MC/fix_widom.cpp
View File

@ -48,20 +48,18 @@
using namespace LAMMPS_NS;
using namespace FixConst;
using namespace MathConst;
using MathConst::MY_2PI;
#define MAXENERGYTEST 1.0e50
enum{EXCHATOM,EXCHMOL}; // exchmode
enum { EXCHATOM, EXCHMOL }; // exchmode
/* ---------------------------------------------------------------------- */
FixWidom::FixWidom(LAMMPS *lmp, int narg, char **arg) :
Fix(lmp, narg, arg),
region(nullptr), idregion(nullptr), full_flag(false), local_gas_list(nullptr),
molcoords(nullptr),molq(nullptr), molimage(nullptr), random_equal(nullptr)
Fix(lmp, narg, arg), region(nullptr), idregion(nullptr), full_flag(false), molcoords(nullptr),
molq(nullptr), molimage(nullptr), random_equal(nullptr), c_pe(nullptr)
{
if (narg < 8) error->all(FLERR,"Illegal fix widom command");
if (narg < 8) utils::missing_cmd_args(FLERR, "fix widom", error);
if (atom->molecular == Atom::TEMPLATE)
error->all(FLERR,"Fix widom does not (yet) work with atom_style template");
@ -83,11 +81,11 @@ FixWidom::FixWidom(LAMMPS *lmp, int narg, char **arg) :
seed = utils::inumeric(FLERR,arg[6],false,lmp);
insertion_temperature = utils::numeric(FLERR,arg[7],false,lmp);
if (nevery <= 0) error->all(FLERR,"Illegal fix widom command");
if (ninsertions < 0) error->all(FLERR,"Illegal fix widom command");
if (seed <= 0) error->all(FLERR,"Illegal fix widom command");
if (nevery <= 0) error->all(FLERR,"Invalid fix widom every argument: {}", nevery);
if (ninsertions < 0) error->all(FLERR,"Invalid fix widom insertions argument: {}", ninsertions);
if (seed <= 0) error->all(FLERR,"Invalid fix widom seed argument: {}", seed);
if (insertion_temperature < 0.0)
error->all(FLERR,"Illegal fix widom command");
error->all(FLERR,"Invalid fix widom temperature argument: {}", insertion_temperature);
// read options from end of input line
@ -99,13 +97,12 @@ FixWidom::FixWidom(LAMMPS *lmp, int narg, char **arg) :
// error checks on region and its extent being inside simulation box
region_xlo = region_xhi = region_ylo = region_yhi =
region_zlo = region_zhi = 0.0;
region_xlo = region_xhi = region_ylo = region_yhi = region_zlo = region_zhi = 0.0;
if (region) {
if (region->bboxflag == 0)
error->all(FLERR,"Fix widom region does not support a bounding box");
error->all(FLERR,"Fix widom region {} does not support a bounding box", region->id);
if (region->dynamic_check())
error->all(FLERR,"Fix widom region cannot be dynamic");
error->all(FLERR,"Fix widom region {} cannot be dynamic", region->id);
region_xlo = region->extent_xlo;
region_xhi = region->extent_xhi;
@ -117,7 +114,7 @@ FixWidom::FixWidom(LAMMPS *lmp, int narg, char **arg) :
if (region_xlo < domain->boxlo[0] || region_xhi > domain->boxhi[0] ||
region_ylo < domain->boxlo[1] || region_yhi > domain->boxhi[1] ||
region_zlo < domain->boxlo[2] || region_zhi > domain->boxhi[2])
error->all(FLERR,"Fix widom region extends outside simulation box");
error->all(FLERR,"Fix widom region {} extends outside simulation box", region->id);
// estimate region volume using MC trials
@ -132,8 +129,8 @@ FixWidom::FixWidom(LAMMPS *lmp, int narg, char **arg) :
inside++;
}
double max_region_volume = (region_xhi - region_xlo) *
(region_yhi - region_ylo) * (region_zhi - region_zlo);
double max_region_volume = (region_xhi - region_xlo) * (region_yhi - region_ylo)
* (region_zhi - region_zlo);
region_volume = max_region_volume * static_cast<double>(inside) / static_cast<double>(attempts);
}
@ -141,28 +138,25 @@ FixWidom::FixWidom(LAMMPS *lmp, int narg, char **arg) :
// error check and further setup for exchmode = EXCHMOL
if (exchmode == EXCHMOL) {
if (onemols[imol]->xflag == 0)
error->all(FLERR,"Fix widom molecule must have coordinates");
if (onemols[imol]->typeflag == 0)
error->all(FLERR,"Fix widom molecule must have atom types");
if (onemol->xflag == 0)
error->all(FLERR,"Fix widom molecule {} must have coordinates", onemol->id);
if (onemol->typeflag == 0)
error->all(FLERR,"Fix widom molecule {} must have atom types", onemol->id);
if (nwidom_type != 0)
error->all(FLERR,"Atom type must be zero in fix widom mol command");
if (onemols[imol]->qflag == 1 && atom->q == nullptr)
error->all(FLERR,"Fix widom molecule has charges, but atom style does not");
if (onemol->qflag == 1 && atom->q == nullptr)
error->all(FLERR,"Fix widom molecule {} has charges, but atom style does not", onemol->id);
if (atom->molecular == Atom::TEMPLATE && onemols != atom->avec->onemols)
error->all(FLERR,"Fix widom molecule template ID must be same "
"as atom_style template ID");
onemols[imol]->check_attributes();
onemol->check_attributes();
}
if (charge_flag && atom->q == nullptr)
error->all(FLERR,"Fix Widom atom has charge, but atom style does not");
error->all(FLERR,"Fix widom atom has charge, but atom style does not");
// setup of array of coordinates for molecule insertion
if (exchmode == EXCHATOM) natoms_per_molecule = 1;
else natoms_per_molecule = onemols[imol]->natoms;
else natoms_per_molecule = onemol->natoms;
nmaxmolatoms = natoms_per_molecule;
grow_molecule_arrays(nmaxmolatoms);
@ -173,7 +167,7 @@ FixWidom::FixWidom(LAMMPS *lmp, int narg, char **arg) :
// zero out counters
widom_nmax = 0;
local_gas_list = nullptr;
ave_widom_chemical_potential = 0.0;
}
/* ----------------------------------------------------------------------
@ -201,17 +195,16 @@ void FixWidom::options(int narg, char **arg)
int iarg = 0;
while (iarg < narg) {
if (strcmp(arg[iarg],"mol") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix widom command");
imol = atom->find_molecule(arg[iarg+1]);
if (imol == -1)
error->all(FLERR,"Molecule template ID for fix widom does not exist");
if (atom->molecules[imol]->nset > 1 && comm->me == 0)
error->warning(FLERR,"Molecule template for "
"fix widom has multiple molecules");
if (strcmp(arg[iarg],"mol") == 0) {
if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "fix widom mol", error);
auto onemols = atom->get_molecule_by_id(arg[iarg+1]);
if (onemols.size() == 0)
error->all(FLERR,"Molecule template ID {} for fix widom does not exist", arg[iarg+1]);
if (onemols.size() > 1 && comm->me == 0)
error->warning(FLERR,"Molecule template {} for fix widom has multiple molecules; "
"will use only the first molecule", arg[iarg+1]);
exchmode = EXCHMOL;
onemols = atom->molecules;
nmol = onemols[imol]->nset;
onemol = onemols[0];
iarg += 2;
} else if (strcmp(arg[iarg],"region") == 0) {
if (iarg+2 > narg) error->all(FLERR,"Illegal fix widom command");
@ -243,7 +236,6 @@ FixWidom::~FixWidom()
delete[] idregion;
delete random_equal;
memory->destroy(local_gas_list);
memory->destroy(molcoords);
memory->destroy(molq);
memory->destroy(molimage);
@ -292,7 +284,7 @@ void FixWidom::init()
}
}
if (full_flag) c_pe = modify->compute[modify->find_compute("thermo_pe")];
if (full_flag) c_pe = modify->get_compute_by_id("thermo_pe");
if (exchmode == EXCHATOM) {
if (nwidom_type <= 0 || nwidom_type > atom->ntypes)
@ -360,22 +352,21 @@ void FixWidom::init()
if (exchmode == EXCHMOL) {
onemols[imol]->compute_mass();
onemols[imol]->compute_com();
gas_mass = onemols[imol]->masstotal;
for (int i = 0; i < onemols[imol]->natoms; i++) {
onemols[imol]->x[i][0] -= onemols[imol]->com[0];
onemols[imol]->x[i][1] -= onemols[imol]->com[1];
onemols[imol]->x[i][2] -= onemols[imol]->com[2];
onemol->compute_mass();
onemol->compute_com();
gas_mass = onemol->masstotal;
for (int i = 0; i < onemol->natoms; i++) {
onemol->x[i][0] -= onemol->com[0];
onemol->x[i][1] -= onemol->com[1];
onemol->x[i][2] -= onemol->com[2];
}
onemols[imol]->com[0] = 0;
onemols[imol]->com[1] = 0;
onemols[imol]->com[2] = 0;
onemol->com[0] = 0;
onemol->com[1] = 0;
onemol->com[2] = 0;
} else gas_mass = atom->mass[nwidom_type];
if (gas_mass <= 0.0)
error->all(FLERR,"Illegal fix widom gas mass <= 0");
if (gas_mass <= 0.0) error->all(FLERR,"Illegal fix widom gas mass <= 0");
// check that no deletable atoms are in atom->firstgroup
// deleting such an atom would not leave firstgroup atoms first
@ -443,7 +434,6 @@ void FixWidom::pre_exchange()
atom->nghost = 0;
comm->borders();
if (triclinic) domain->lamda2x(atom->nlocal+atom->nghost);
update_gas_atoms_list();
if (full_flag) {
energy_stored = energy_full();
@ -525,12 +515,12 @@ void FixWidom::attempt_atomic_insertion()
if (!domain->inside(coord))
error->one(FLERR,"Fix widom put atom outside box");
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]) proc_flag = 1;
coord[1] >= sublo[1] && coord[1] < subhi[1] &&
coord[2] >= sublo[2] && coord[2] < subhi[2]) proc_flag = 1;
} else {
if (lamda[0] >= sublo[0] && lamda[0] < subhi[0] &&
lamda[1] >= sublo[1] && lamda[1] < subhi[1] &&
lamda[2] >= sublo[2] && lamda[2] < subhi[2]) proc_flag = 1;
lamda[1] >= sublo[1] && lamda[1] < subhi[1] &&
lamda[2] >= sublo[2] && lamda[2] < subhi[2]) proc_flag = 1;
}
if (proc_flag) {
@ -545,9 +535,7 @@ void FixWidom::attempt_atomic_insertion()
double incr_chem_pot = (inst_chem_pot - ave_widom_chemical_potential);
ave_widom_chemical_potential += incr_chem_pot / (imove + 1);
}
}
}
/* ----------------------------------------------------------------------
@ -571,13 +559,13 @@ void FixWidom::attempt_molecule_insertion()
com_coord[2] = region_zlo + random_equal->uniform() *
(region_zhi-region_zlo);
while (region->match(com_coord[0],com_coord[1],
com_coord[2]) == 0) {
com_coord[2]) == 0) {
com_coord[0] = region_xlo + random_equal->uniform() *
(region_xhi-region_xlo);
(region_xhi-region_xlo);
com_coord[1] = region_ylo + random_equal->uniform() *
(region_yhi-region_ylo);
(region_yhi-region_ylo);
com_coord[2] = region_zlo + random_equal->uniform() *
(region_zhi-region_zlo);
(region_zhi-region_zlo);
region_attempt++;
if (region_attempt >= max_region_attempts) return;
}
@ -622,7 +610,7 @@ void FixWidom::attempt_molecule_insertion()
auto procflag = new bool[natoms_per_molecule];
for (int i = 0; i < natoms_per_molecule; i++) {
MathExtra::matvec(rotmat,onemols[imol]->x[i],molcoords[i]);
MathExtra::matvec(rotmat,onemol->x[i],molcoords[i]);
molcoords[i][0] += com_coord[0];
molcoords[i][1] += com_coord[1];
molcoords[i][2] += com_coord[2];
@ -641,40 +629,37 @@ void FixWidom::attempt_molecule_insertion()
procflag[i] = false;
if (triclinic == 0) {
if (xtmp[0] >= sublo[0] && xtmp[0] < subhi[0] &&
xtmp[1] >= sublo[1] && xtmp[1] < subhi[1] &&
xtmp[2] >= sublo[2] && xtmp[2] < subhi[2]) procflag[i] = true;
xtmp[1] >= sublo[1] && xtmp[1] < subhi[1] &&
xtmp[2] >= sublo[2] && xtmp[2] < subhi[2]) procflag[i] = true;
} else {
domain->x2lamda(xtmp,lamda);
if (lamda[0] >= sublo[0] && lamda[0] < subhi[0] &&
lamda[1] >= sublo[1] && lamda[1] < subhi[1] &&
lamda[2] >= sublo[2] && lamda[2] < subhi[2]) procflag[i] = true;
lamda[1] >= sublo[1] && lamda[1] < subhi[1] &&
lamda[2] >= sublo[2] && lamda[2] < subhi[2]) procflag[i] = true;
}
if (procflag[i]) {
int ii = -1;
if (onemols[imol]->qflag == 1) {
ii = atom->nlocal + atom->nghost;
if (ii >= atom->nmax) atom->avec->grow(0);
atom->q[ii] = onemols[imol]->q[i];
if (onemol->qflag == 1) {
ii = atom->nlocal + atom->nghost;
if (ii >= atom->nmax) atom->avec->grow(0);
atom->q[ii] = onemol->q[i];
}
insertion_energy += energy(ii,onemols[imol]->type[i],-1,xtmp);
insertion_energy += energy(ii,onemol->type[i],-1,xtmp);
}
}
double insertion_energy_sum = 0.0;
MPI_Allreduce(&insertion_energy,&insertion_energy_sum,1,
MPI_DOUBLE,MPI_SUM,world);
MPI_DOUBLE,MPI_SUM,world);
// the insertion_energy_sum is the variable with the energy of inserting one molecule
double inst_chem_pot = exp(-insertion_energy_sum*beta);
double incr_chem_pot = (inst_chem_pot - ave_widom_chemical_potential);
ave_widom_chemical_potential += incr_chem_pot / (imove + 1);
delete[] procflag;
}
}
/* ----------------------------------------------------------------------
@ -727,12 +712,12 @@ void FixWidom::attempt_atomic_insertion_full()
if (!domain->inside(coord))
error->one(FLERR,"Fix widom put atom outside box");
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]) proc_flag = 1;
coord[1] >= sublo[1] && coord[1] < subhi[1] &&
coord[2] >= sublo[2] && coord[2] < subhi[2]) proc_flag = 1;
} else {
if (lamda[0] >= sublo[0] && lamda[0] < subhi[0] &&
lamda[1] >= sublo[1] && lamda[1] < subhi[1] &&
lamda[2] >= sublo[2] && lamda[2] < subhi[2]) proc_flag = 1;
lamda[1] >= sublo[1] && lamda[1] < subhi[1] &&
lamda[2] >= sublo[2] && lamda[2] < subhi[2]) proc_flag = 1;
}
if (proc_flag) {
@ -770,11 +755,7 @@ void FixWidom::attempt_atomic_insertion_full()
if (proc_flag) atom->nlocal--;
if (force->kspace) force->kspace->qsum_qsq();
if (force->pair->tail_flag) force->pair->reinit();
update_gas_atoms_list();
}
}
/* ----------------------------------------------------------------------
@ -803,20 +784,13 @@ void FixWidom::attempt_molecule_insertion_full()
double com_coord[3];
if (region) {
int region_attempt = 0;
com_coord[0] = region_xlo + random_equal->uniform() *
(region_xhi-region_xlo);
com_coord[1] = region_ylo + random_equal->uniform() *
(region_yhi-region_ylo);
com_coord[2] = region_zlo + random_equal->uniform() *
(region_zhi-region_zlo);
while (region->match(com_coord[0],com_coord[1],
com_coord[2]) == 0) {
com_coord[0] = region_xlo + random_equal->uniform() *
(region_xhi-region_xlo);
com_coord[1] = region_ylo + random_equal->uniform() *
(region_yhi-region_ylo);
com_coord[2] = region_zlo + random_equal->uniform() *
(region_zhi-region_zlo);
com_coord[0] = region_xlo + random_equal->uniform() * (region_xhi-region_xlo);
com_coord[1] = region_ylo + random_equal->uniform() * (region_yhi-region_ylo);
com_coord[2] = region_zlo + random_equal->uniform() * (region_zhi-region_zlo);
while (region->match(com_coord[0],com_coord[1], com_coord[2]) == 0) {
com_coord[0] = region_xlo + random_equal->uniform() * (region_xhi-region_xlo);
com_coord[1] = region_ylo + random_equal->uniform() * (region_yhi-region_ylo);
com_coord[2] = region_zlo + random_equal->uniform() * (region_zhi-region_zlo);
region_attempt++;
if (region_attempt >= max_region_attempts) return;
}
@ -861,7 +835,7 @@ void FixWidom::attempt_molecule_insertion_full()
for (int i = 0; i < natoms_per_molecule; i++) {
double xtmp[3];
MathExtra::matvec(rotmat,onemols[imol]->x[i],xtmp);
MathExtra::matvec(rotmat,onemol->x[i],xtmp);
xtmp[0] += com_coord[0];
xtmp[1] += com_coord[1];
xtmp[2] += com_coord[2];
@ -876,40 +850,39 @@ void FixWidom::attempt_molecule_insertion_full()
int proc_flag = 0;
if (triclinic == 0) {
if (xtmp[0] >= sublo[0] && xtmp[0] < subhi[0] &&
xtmp[1] >= sublo[1] && xtmp[1] < subhi[1] &&
xtmp[2] >= sublo[2] && xtmp[2] < subhi[2]) proc_flag = 1;
xtmp[1] >= sublo[1] && xtmp[1] < subhi[1] &&
xtmp[2] >= sublo[2] && xtmp[2] < subhi[2]) proc_flag = 1;
} else {
domain->x2lamda(xtmp,lamda);
if (lamda[0] >= sublo[0] && lamda[0] < subhi[0] &&
lamda[1] >= sublo[1] && lamda[1] < subhi[1] &&
lamda[2] >= sublo[2] && lamda[2] < subhi[2]) proc_flag = 1;
lamda[1] >= sublo[1] && lamda[1] < subhi[1] &&
lamda[2] >= sublo[2] && lamda[2] < subhi[2]) proc_flag = 1;
}
if (proc_flag) {
atom->avec->create_atom(onemols[imol]->type[i],xtmp);
atom->avec->create_atom(onemol->type[i],xtmp);
int m = atom->nlocal - 1;
atom->image[m] = imagetmp;
atom->molecule[m] = insertion_molecule;
if (maxtag_all+i+1 >= MAXTAGINT)
error->all(FLERR,"Fix widom ran out of available atom IDs");
error->all(FLERR,"Fix widom ran out of available atom IDs");
atom->tag[m] = maxtag_all + i + 1;
atom->v[m][0] = 0;
atom->v[m][1] = 0;
atom->v[m][2] = 0;
atom->add_molecule_atom(onemols[imol],i,m,maxtag_all);
atom->add_molecule_atom(onemol,i,m,maxtag_all);
modify->create_attribute(m);
}
}
atom->natoms += natoms_per_molecule;
if (atom->natoms < 0)
error->all(FLERR,"Too many total atoms");
atom->nbonds += onemols[imol]->nbonds;
atom->nangles += onemols[imol]->nangles;
atom->ndihedrals += onemols[imol]->ndihedrals;
atom->nimpropers += onemols[imol]->nimpropers;
if (atom->natoms < 0) error->all(FLERR,"Too many total atoms");
atom->nbonds += onemol->nbonds;
atom->nangles += onemol->nangles;
atom->ndihedrals += onemol->ndihedrals;
atom->nimpropers += onemol->nimpropers;
if (atom->map_style != Atom::MAP_NONE) atom->map_init();
atom->nghost = 0;
if (triclinic) domain->x2lamda(atom->nlocal);
@ -924,10 +897,10 @@ void FixWidom::attempt_molecule_insertion_full()
double incr_chem_pot = (inst_chem_pot - ave_widom_chemical_potential);
ave_widom_chemical_potential += incr_chem_pot / (imove + 1);
atom->nbonds -= onemols[imol]->nbonds;
atom->nangles -= onemols[imol]->nangles;
atom->ndihedrals -= onemols[imol]->ndihedrals;
atom->nimpropers -= onemols[imol]->nimpropers;
atom->nbonds -= onemol->nbonds;
atom->nangles -= onemol->nangles;
atom->ndihedrals -= onemol->ndihedrals;
atom->nimpropers -= onemol->nimpropers;
atom->natoms -= natoms_per_molecule;
int i = 0;
@ -939,9 +912,6 @@ void FixWidom::attempt_molecule_insertion_full()
}
if (force->kspace) force->kspace->qsum_qsq();
if (force->pair->tail_flag) force->pair->reinit();
update_gas_atoms_list();
}
}
@ -1061,92 +1031,6 @@ double FixWidom::energy_full()
return total_energy;
}
/* ----------------------------------------------------------------------
update the list of gas atoms
------------------------------------------------------------------------- */
void FixWidom::update_gas_atoms_list()
{
int nlocal = atom->nlocal;
int *mask = atom->mask;
tagint *molecule = atom->molecule;
double **x = atom->x;
if (atom->nmax > widom_nmax) {
memory->sfree(local_gas_list);
widom_nmax = atom->nmax;
local_gas_list = (int *) memory->smalloc(widom_nmax*sizeof(int),
"Widom:local_gas_list");
}
ngas_local = 0;
if (region) {
if (exchmode == EXCHMOL) {
tagint maxmol = 0;
for (int i = 0; i < nlocal; i++) maxmol = MAX(maxmol,molecule[i]);
tagint maxmol_all;
MPI_Allreduce(&maxmol,&maxmol_all,1,MPI_LMP_TAGINT,MPI_MAX,world);
auto comx = new double[maxmol_all];
auto comy = new double[maxmol_all];
auto comz = new double[maxmol_all];
for (int imolecule = 0; imolecule < maxmol_all; imolecule++) {
for (int i = 0; i < nlocal; i++) {
if (molecule[i] == imolecule) {
mask[i] |= molecule_group_bit;
} else {
mask[i] &= molecule_group_inversebit;
}
}
double com[3];
com[0] = com[1] = com[2] = 0.0;
group->xcm(molecule_group,gas_mass,com);
// remap unwrapped com into periodic box
domain->remap(com);
comx[imolecule] = com[0];
comy[imolecule] = com[1];
comz[imolecule] = com[2];
}
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
if (region->match(comx[molecule[i]],
comy[molecule[i]],comz[molecule[i]]) == 1) {
local_gas_list[ngas_local] = i;
ngas_local++;
}
}
}
delete[] comx;
delete[] comy;
delete[] comz;
} else {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
if (region->match(x[i][0],x[i][1],x[i][2]) == 1) {
local_gas_list[ngas_local] = i;
ngas_local++;
}
}
}
}
} else {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
local_gas_list[ngas_local] = i;
ngas_local++;
}
}
}
MPI_Allreduce(&ngas_local,&ngas,1,MPI_INT,MPI_SUM,world);
}
/* ----------------------------------------------------------------------
return acceptance ratios
------------------------------------------------------------------------- */

10
src/MC/fix_widom.h
View File

@ -40,7 +40,6 @@ class FixWidom : public Fix {
double energy(int, int, tagint, double *);
double molecule_energy(tagint);
double energy_full();
void update_gas_atoms_list();
double compute_vector(int) override;
double memory_usage() override;
void write_restart(FILE *) override;
@ -53,10 +52,7 @@ class FixWidom : public Fix {
int exclusion_group, exclusion_group_bit;
int nwidom_type, nevery, seed;
int ninsertions;
int ngas; // # of gas atoms on all procs
int ngas_local; // # of gas atoms on this proc
int exchmode; // exchange ATOM or MOLECULE
int movemode; // move ATOM or MOLECULE
class Region *region; // widom region
char *idregion; // widom region id
bool charge_flag; // true if user specified atomic charge
@ -71,14 +67,13 @@ class FixWidom : public Fix {
int max_region_attempts;
double gas_mass;
double insertion_temperature;
double beta, sigma, volume;
double beta, volume;
double charge;
double xlo, xhi, ylo, yhi, zlo, zhi;
double region_xlo, region_xhi, region_ylo, region_yhi, region_zlo, region_zhi;
double region_volume;
double energy_stored; // full energy of old/current configuration
double *sublo, *subhi;
int *local_gas_list;
double **cutsq;
double **molcoords;
double *molq;
@ -93,8 +88,7 @@ class FixWidom : public Fix {
class Atom *model_atom;
class Molecule **onemols;
int imol, nmol;
class Molecule *onemol;
int triclinic; // 0 = orthog box, 1 = triclinic
class Compute *c_pe;

7
src/OPENMP/fix_nvt_sllod_omp.cpp
View File

@ -56,9 +56,9 @@ FixNVTSllodOMP::FixNVTSllodOMP(LAMMPS *lmp, int narg, char **arg) :
// id = fix-ID + temp
id_temp = utils::strdup(std::string(id) + "_temp");
modify->add_compute(fmt::format("{} {} temp/deform",
id_temp,group->names[igroup]));
modify->add_compute(fmt::format("{} {} temp/deform",id_temp,group->names[igroup]));
tcomputeflag = 1;
nondeformbias = 0;
}
/* ---------------------------------------------------------------------- */
@ -79,8 +79,7 @@ void FixNVTSllodOMP::init()
for (i = 0; i < modify->nfix; i++)
if (utils::strmatch(modify->fix[i]->style,"^deform")) {
if ((dynamic_cast<FixDeform *>(modify->fix[i]))->remapflag != Domain::V_REMAP)
error->all(FLERR,"Using fix nvt/sllod/omp with inconsistent fix "
"deform remap option");
error->all(FLERR,"Using fix nvt/sllod/omp with inconsistent fix deform remap option");
break;
}
if (i == modify->nfix)

8
src/atom.cpp
View File

@ -779,9 +779,11 @@ std::string Atom::get_style()
std::string retval = atom_style;
if (retval == "hybrid") {
auto avec_hybrid = dynamic_cast<AtomVecHybrid *>(avec);
for (int i = 0; i < avec_hybrid->nstyles; i++) {
retval += ' ';
retval += avec_hybrid->keywords[i];
if (avec_hybrid) {
for (int i = 0; i < avec_hybrid->nstyles; i++) {
retval += ' ';
retval += avec_hybrid->keywords[i];
}
}
}
return retval;

7
src/fix_nvt_sllod.cpp
View File

@ -51,9 +51,9 @@ FixNVTSllod::FixNVTSllod(LAMMPS *lmp, int narg, char **arg) :
// id = fix-ID + temp
id_temp = utils::strdup(std::string(id) + "_temp");
modify->add_compute(fmt::format("{} {} temp/deform",
id_temp,group->names[igroup]));
modify->add_compute(fmt::format("{} {} temp/deform",id_temp,group->names[igroup]));
tcomputeflag = 1;
nondeformbias = 0;
}
/* ---------------------------------------------------------------------- */
@ -74,8 +74,7 @@ void FixNVTSllod::init()
for (i = 0; i < modify->nfix; i++)
if (strncmp(modify->fix[i]->style,"deform",6) == 0) {
if ((dynamic_cast<FixDeform *>(modify->fix[i]))->remapflag != Domain::V_REMAP)
error->all(FLERR,"Using fix nvt/sllod with inconsistent fix deform "
"remap option");
error->all(FLERR,"Using fix nvt/sllod with inconsistent fix deform remap option");
break;
}
if (i == modify->nfix)

19
src/fix_pair.cpp
View File

@ -120,7 +120,8 @@ FixPair::FixPair(LAMMPS *lmp, int narg, char **arg) :
atom->add_callback(Atom::GROW);
// zero the vector/array since dump may access it on timestep 0
// zero the vector/array since a variable may access it before first run
// zero the vector/array since a variable may access it before first ru
// initialize lasttime so step 0 will trigger/extract
int nlocal = atom->nlocal;
@ -132,6 +133,8 @@ FixPair::FixPair(LAMMPS *lmp, int narg, char **arg) :
for (int m = 0; m < ncols; m++)
array[i][m] = 0.0;
}
lasttime = -1;
}
/* ---------------------------------------------------------------------- */
@ -188,6 +191,13 @@ void FixPair::setup(int vflag)
/* ---------------------------------------------------------------------- */
void FixPair::min_setup(int vflag)
{
setup(vflag);
}
/* ---------------------------------------------------------------------- */
void FixPair::setup_pre_force(int vflag)
{
pre_force(vflag);
@ -195,11 +205,13 @@ void FixPair::setup_pre_force(int vflag)
/* ----------------------------------------------------------------------
trigger pair style computation on steps which are multiples of Nevery
lasttime prevents mulitiple triggers by min linesearch on same iteration
------------------------------------------------------------------------- */
void FixPair::pre_force(int /*vflag*/)
{
if (update->ntimestep % nevery) return;
if (update->ntimestep == lasttime) return;
// set pair style triggers
@ -215,12 +227,15 @@ void FixPair::min_pre_force(int vflag)
}
/* ----------------------------------------------------------------------
extract results from pair style
extract results from pair style on steps which are multiples of Nevery
lasttime prevents mulitiple extracts by min linesearch on same iteration
------------------------------------------------------------------------- */
void FixPair::post_force(int /*vflag*/)
{
if (update->ntimestep % nevery) return;
if (update->ntimestep == lasttime) return;
lasttime = update->ntimestep;
// extract pair style fields one by one
// store their values in this fix

2
src/fix_pair.h
View File

@ -31,6 +31,7 @@ class FixPair : public Fix {
int setmask() override;
void init() override;
void setup(int) override;
void min_setup(int) override;
void setup_pre_force(int) override;
void pre_force(int) override;
void min_pre_force(int) override;
@ -46,6 +47,7 @@ class FixPair : public Fix {
private:
int nevery,nfield,ncols;
bigint lasttime;
char *pairname;
char **fieldname,**triggername;
int *trigger;

3
src/neigh_request.cpp
View File

@ -99,9 +99,8 @@ NeighRequest::NeighRequest(LAMMPS *_lmp) : Pointers(_lmp)
/* ---------------------------------------------------------------------- */
NeighRequest::NeighRequest(LAMMPS *_lmp, int idx, void *ptr, int num) : NeighRequest(_lmp)
NeighRequest::NeighRequest(LAMMPS *_lmp, void *ptr, int num) : NeighRequest(_lmp)
{
index = idx;
requestor = ptr;
requestor_instance = num;
}

3
src/neigh_request.h
View File

@ -29,7 +29,6 @@ class NeighRequest : protected Pointers {
friend class FixIntel;
protected:
int index; // index of which neigh request this is
void *requestor; // class that made request
int requestor_instance; // instance of that class (only Fix, Compute, Pair)
int id; // ID of request as stored by requestor
@ -122,7 +121,7 @@ class NeighRequest : protected Pointers {
// methods
public:
NeighRequest(class LAMMPS *);
NeighRequest(class LAMMPS *, int, void *, int);
NeighRequest(class LAMMPS *, void *, int);
NeighRequest(NeighRequest *);
~NeighRequest() override;

4
src/neighbor.cpp
View File

@ -2163,13 +2163,13 @@ int Neighbor::request(void *requestor, int instance)
memory->srealloc(requests,maxrequest*sizeof(NeighRequest *), "neighbor:requests");
}
requests[nrequest] = new NeighRequest(lmp, nrequest, requestor, instance);
requests[nrequest] = new NeighRequest(lmp, requestor, instance);
nrequest++;
return nrequest-1;
}
/* ----------------------------------------------------------------------
called by other classes to request a pairwise neighbor list
add_request() is called by other classes to request a pairwise neighbor list
------------------------------------------------------------------------- */
NeighRequest *Neighbor::add_request(Pair *requestor, int flags)

2
unittest/fortran/CMakeLists.txt
View File

@ -33,7 +33,7 @@ if(CMAKE_Fortran_COMPILER)
add_library(flammps STATIC ${LAMMPS_FORTRAN_MODULE})
add_executable(test_fortran_create wrap_create.cpp test_fortran_create.f90)
target_link_libraries(test_fortran_create PRIVATE lammps MPI::MPI_Fortran GTest::GTestMain)
target_link_libraries(test_fortran_create PRIVATE flammps lammps MPI::MPI_Fortran GTest::GTestMain)
target_include_directories(test_fortran_create PRIVATE "${LAMMPS_SOURCE_DIR}/../fortran")
add_test(NAME FortranOpen COMMAND test_fortran_create)