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

Merge branch 'develop' into Iximiel/develop

Browse Source
This commit is contained in:
Axel Kohlmeyer
2022-04-27 14:15:10 -04:00
parent 47f3f9d44a 691ba89b6f
commit 75d42a87e5
4535 changed files with 343096 additions and 522642 deletions
Download Patch File Download Diff File Expand all files Collapse all files

53
src/.gitignore vendored
View File

@ -1,6 +1,7 @@
/Makefile.package
/Makefile.package.settings
/MAKE/MINE
/ADIOS/Makefile.lammps
/Make.py.last
/lmp_*
@ -253,6 +254,21 @@
/pair_mesont_tpm.cpp
/pair_mesont_tpm.h
/atom_vec_bpm_sphere.cpp
/atom_vec_bpm_sphere.h
/bond_bpm.cpp
/bond_bpm.h
/bond_bpm_rotational.cpp
/bond_bpm_rotational.h
/bond_bpm_spring.cpp
/bond_bpm_spring.h
/compute_nbond_atom.cpp
/compute_nbond_atom.h
/fix_nve_bpm_sphere.cpp
/fix_nve_bpm_sphere.h
/pair_bpm_spring.cpp
/pair_bpm_spring.h
/compute_adf.cpp
/compute_adf.h
/compute_contact_atom.cpp
@ -267,6 +283,8 @@
/fix_drag.h
/fix_numdiff.cpp
/fix_numdiff.h
/fix_numdiff_virial.cpp
/fix_numdiff_virial.h
/fix_spring_rg.cpp
/fix_spring_rg.h
/fix_temp_csld.cpp
@ -393,6 +411,8 @@
/bond_fene.h
/bond_fene_expand.cpp
/bond_fene_expand.h
/bond_fene_nm.cpp
/bond_fene_nm.h
/bond_gaussian.cpp
/bond_gaussian.h
/bond_gromos.cpp
@ -427,10 +447,14 @@
/commgrid.h
/compute_ackland_atom.cpp
/compute_ackland_atom.h
/compute_ave_sphere_atom.cpp
/compute_ave_sphere_atom.h
/compute_basal_atom.cpp
/compute_basal_atom.h
/compute_body_local.cpp
/compute_body_local.h
/compute_born_matrix.cpp
/compute_born_matrix.h
/compute_cnp_atom.cpp
/compute_cnp_atom.h
/compute_damage_atom.cpp
@ -453,6 +477,8 @@
/compute_fabric.h
/compute_fep.cpp
/compute_fep.h
/compute_fep_ta.cpp
/compute_fep_ta.h
/compute_force_tally.cpp
/compute_force_tally.h
/compute_gyration_shape.cpp
@ -485,8 +511,6 @@
/compute_plasticity_atom.h
/compute_pressure_bocs.cpp
/compute_pressure_bocs.h
/compute_pressure_cylinder.cpp
/compute_pressure_cylinder.h
/compute_pressure_grem.cpp
/compute_pressure_grem.h
/compute_ptm_atom.cpp
@ -497,10 +521,16 @@
/compute_smd_triangle_vertices.h
/compute_spec_atom.cpp
/compute_spec_atom.h
/compute_stress_cartesian.cpp
/compute_stress_cartesian.h
/compute_stress_cylinder.cpp
/compute_stress_cylinder.h
/compute_stress_mop.cpp
/compute_stress_mop.h
/compute_stress_mop_profile.cpp
/compute_stress_mop_profile.h
/compute_stress_spherical.cpp
/compute_stress_spherical.h
/compute_stress_tally.cpp
/compute_stress_tally.h
/compute_temp_asphere.cpp
@ -589,6 +619,8 @@
/dump_xyz_mpiio.h
/dump_xyz_zstd.cpp
/dump_xyz_zstd.h
/dump_yaml.cpp
/dump_yaml.h
/dynamical_matrix.cpp
/dynamical_matrix.h
/ewald.cpp
@ -777,8 +809,6 @@
/fix_orient_eco.h
/fix_orient_fcc.cpp
/fix_orient_fcc.h
/fix_pair_tracker.cpp
/fix_pair_tracker.h
/fix_peri_neigh.cpp
/fix_peri_neigh.h
/fix_phonon.cpp
@ -926,6 +956,8 @@
/msm_cg.h
/neb.cpp
/neb.h
/netcdf_units.cpp
/netcdf_units.h
/pair_adp.cpp
/pair_adp.h
/pair_agni.cpp
@ -1061,6 +1093,10 @@
/pair_hdnnp.h
/pair_ilp_graphene_hbn.cpp
/pair_ilp_graphene_hbn.h
/pair_ilp_graphene_hbn_opt.cpp
/pair_ilp_graphene_hbn_opt.h
/pair_ilp_tmd.cpp
/pair_ilp_tmd.h
/pair_kolmogorov_crespi_full.cpp
/pair_kolmogorov_crespi_full.h
/pair_kolmogorov_crespi_z.cpp
@ -1167,6 +1203,8 @@
/pair_nm_cut_coul_cut.h
/pair_nm_cut_coul_long.cpp
/pair_nm_cut_coul_long.h
/pair_nm_cut_split.cpp
/pair_nm_cut_split.h
/pair_oxdna_*.cpp
/pair_oxdna_*.h
/pair_oxdna2_*.cpp
@ -1192,6 +1230,8 @@
/pair_rebo.h
/pair_resquared.cpp
/pair_resquared.h
/pair_saip_metal.cpp
/pair_saip_metal.h
/pair_sdpd_taitwater_isothermal.cpp
/pair_sdpd_taitwater_isothermal.h
/pair_sph_heatconduction.cpp
@ -1208,6 +1248,8 @@
/pair_sph_taitwater_morris.h
/pair_sw.cpp
/pair_sw.h
/pair_sw_mod.cpp
/pair_sw_mod.h
/pair_tersoff.cpp
/pair_tersoff.h
/pair_tersoff_mod.cpp
@ -1481,6 +1523,8 @@
/pair_dpd_fdt.h
/pair_dpd_fdt_energy.cpp
/pair_dpd_fdt_energy.h
/pair_harmonic_cut.cpp
/pair_harmonic_cut.h
/pair_lennard_mdf.cpp
/pair_lennard_mdf.h
/pair_lj_cut_coul_long_cs.cpp
@ -1497,3 +1541,4 @@
/pair_smtbq.h
/pair_vashishta*.cpp
/pair_vashishta*.h

6
src/ADIOS/Install.sh
View File

@ -31,17 +31,17 @@ done
if (test $1 = 1) then
CONFIGSCRIPT=none
if ( test `which adios2-config 2>> /dev/null` ) then
if ( type adios2-config > /dev/null 2>&1 ) then
CONFIGSCRIPT=adios2-config
elif ( ! test -z "$ADIOS2_DIR" ) then
if ( test `which $ADIOS2_DIR/bin/adios2-config` ) then
if ( type $ADIOS2_DIR/bin/adios2-config > /dev/null 2>&1 ) then
CONFIGSCRIPT=$ADIOS2_DIR/bin/adios2-config
else
echo "ERROR: ADIOS2_DIR environment variable is set but" \
"\$ADIOS2_DIR/bin/adios2-config does not exist"
fi
elif ( ! test -z "$ADIOS_DIR" ) then
if ( test `which $ADIOS_DIR/bin/adios2-config` ) then
if ( type $ADIOS_DIR/bin/adios2-config > /dev/null 2>&1 ) then
CONFIGSCRIPT=$ADIOS_DIR/bin/adios2-config
else
echo "ERROR: ADIOS_DIR environment variable is set but" \

39
src/ADIOS/adios_common.h Normal file
View File

@ -0,0 +1,39 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
Contributed by Norbert Podhorszki (Oak Ridge National Laboratory)
------------------------------------------------------------------------- */
#ifndef LMP_ADIOS_COMMON_H
#define LMP_ADIOS_COMMON_H
// common definitions for all ADIOS package classes
static const char default_config[] = "<?xml version=\"1.0\"?>\n"
"<adios-config>\n"
" <io name=\"atom\">\n"
" <engine type=\"BP4\">\n"
" <parameter key=\"substreams\" value=\"1\"/>\n"
" </engine>\n"
" </io>\n"
" <io name=\"custom\">\n"
" <engine type=\"BP4\">\n"
" <parameter key=\"substreams\" value=\"1\"/>\n"
" </engine>\n"
" </io>\n"
" <io name=\"read_dump\">\n"
" <engine type=\"BP4\">\n"
" </engine>\n"
" </io>\n"
"</adios-config>\n";
#endif

500
src/ADIOS/dump_atom_adios.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -24,333 +23,294 @@
#include "memory.h"
#include "universe.h"
#include "update.h"
#include <cstring>
#include "adios2.h"
#include "adios_common.h"
using namespace LAMMPS_NS;
#define MAX_TEXT_HEADER_SIZE 4096
#define DUMP_BUF_CHUNK_SIZE 16384
#define DUMP_BUF_INCREMENT_SIZE 4096
namespace LAMMPS_NS {
class DumpAtomADIOSInternal {
namespace LAMMPS_NS
{
class DumpAtomADIOSInternal
{
public:
DumpAtomADIOSInternal(){};
~DumpAtomADIOSInternal() = default;
public:
DumpAtomADIOSInternal() {};
~DumpAtomADIOSInternal() = default;
// name of adios group, referrable in adios2_config.xml
const std::string ioName = "atom";
adios2::ADIOS *ad = nullptr; // adios object
adios2::IO io; // adios group of variables and attributes in this dump
adios2::Engine fh; // adios file/stream handle object
// one ADIOS output variable we need to change every step
adios2::Variable<double> varAtoms;
// name of adios group, referrable in adios2_config.xml
const std::string ioName = "atom";
adios2::ADIOS *ad = nullptr; // adios object
adios2::IO io; // adios group of variables and attributes in this dump
adios2::Engine fh; // adios file/stream handle object
// one ADIOS output variable we need to change every step
adios2::Variable<double> varAtoms;
};
} // namespace LAMMPS_NS
} // namespace LAMMPS_NS
/* ---------------------------------------------------------------------- */
DumpAtomADIOS::DumpAtomADIOS(LAMMPS *lmp, int narg, char **arg)
: DumpAtom(lmp, narg, arg)
DumpAtomADIOS::DumpAtomADIOS(LAMMPS *lmp, int narg, char **arg) : DumpAtom(lmp, narg, arg)
{
internal = new DumpAtomADIOSInternal();
try {
internal->ad =
new adios2::ADIOS("adios2_config.xml", world, adios2::DebugON);
} catch (std::ios_base::failure &e) {
char str[256];
snprintf(str, sizeof(str), "ADIOS initialization failed with error: %s",
e.what());
error->one(FLERR, str);
}
// create a default adios2_config.xml if it doesn't exist yet.
FILE *cfgfp = fopen("adios2_config.xml", "r");
if (!cfgfp) {
cfgfp = fopen("adios2_config.xml", "w");
if (cfgfp) fputs(default_config, cfgfp);
}
if (cfgfp) fclose(cfgfp);
internal = new DumpAtomADIOSInternal();
try {
internal->ad = new adios2::ADIOS("adios2_config.xml", world, adios2::DebugON);
} catch (std::ios_base::failure &e) {
error->all(FLERR, "ADIOS initialization failed with error: {}", e.what());
}
}
/* ---------------------------------------------------------------------- */
DumpAtomADIOS::~DumpAtomADIOS()
{
if (internal->fh) {
internal->fh.Close();
}
delete internal->ad;
delete internal;
if (internal->fh) internal->fh.Close();
delete internal->ad;
delete internal;
}
/* ---------------------------------------------------------------------- */
void DumpAtomADIOS::openfile()
{
if (multifile) {
// if one file per timestep, replace '*' with current timestep
char *filestar = strdup(filename);
char *filecurrent = new char[strlen(filestar) + 16];
char *ptr = strchr(filestar, '*');
*ptr = '\0';
if (padflag == 0)
snprintf(filecurrent, sizeof(filecurrent), "%s" BIGINT_FORMAT "%s",
filestar, update->ntimestep, ptr + 1);
else {
char bif[8], pad[16];
strcpy(bif, BIGINT_FORMAT);
snprintf(pad, sizeof(pad), "%%s%%0%d%s%%s", padflag, &bif[1]);
snprintf(filecurrent, sizeof(filecurrent), pad, filestar,
update->ntimestep, ptr + 1);
}
internal->fh =
internal->io.Open(filecurrent, adios2::Mode::Write, world);
if (!internal->fh) {
char str[128];
snprintf(str, sizeof(str), "Cannot open dump file %s", filecurrent);
error->one(FLERR, str);
}
free(filestar);
delete[] filecurrent;
} else {
if (!singlefile_opened) {
internal->fh =
internal->io.Open(filename, adios2::Mode::Write, world);
if (!internal->fh) {
char str[128];
snprintf(str, sizeof(str), "Cannot open dump file %s",
filename);
error->one(FLERR, str);
}
singlefile_opened = 1;
}
if (multifile) {
// if one file per timestep, replace '*' with current timestep
auto filecurrent = utils::star_subst(filename, update->ntimestep, padflag);
internal->fh = internal->io.Open(filecurrent, adios2::Mode::Write, world);
if (!internal->fh) error->one(FLERR, "Cannot open dump file {}", filecurrent);
} else {
if (!singlefile_opened) {
internal->fh = internal->io.Open(filename, adios2::Mode::Write, world);
if (!internal->fh) error->one(FLERR, "Cannot open dump file {}", filename);
singlefile_opened = 1;
}
}
}
/* ---------------------------------------------------------------------- */
void DumpAtomADIOS::write()
{
if (domain->triclinic == 0) {
boxxlo = domain->boxlo[0];
boxxhi = domain->boxhi[0];
boxylo = domain->boxlo[1];
boxyhi = domain->boxhi[1];
boxzlo = domain->boxlo[2];
boxzhi = domain->boxhi[2];
} else {
boxxlo = domain->boxlo_bound[0];
boxxhi = domain->boxhi_bound[0];
boxylo = domain->boxlo_bound[1];
boxyhi = domain->boxhi_bound[1];
boxzlo = domain->boxlo_bound[2];
boxzhi = domain->boxhi_bound[2];
boxxy = domain->xy;
boxxz = domain->xz;
boxyz = domain->yz;
if (domain->triclinic == 0) {
boxxlo = domain->boxlo[0];
boxxhi = domain->boxhi[0];
boxylo = domain->boxlo[1];
boxyhi = domain->boxhi[1];
boxzlo = domain->boxlo[2];
boxzhi = domain->boxhi[2];
} else {
boxxlo = domain->boxlo_bound[0];
boxxhi = domain->boxhi_bound[0];
boxylo = domain->boxlo_bound[1];
boxyhi = domain->boxhi_bound[1];
boxzlo = domain->boxlo_bound[2];
boxzhi = domain->boxhi_bound[2];
boxxy = domain->xy;
boxxz = domain->xz;
boxyz = domain->yz;
}
// nme = # of dump lines this proc contributes to dump
nme = count();
// ntotal = total # of atoms in snapshot
// atomOffset = sum of # of atoms up to this proc (exclusive prefix sum)
bigint bnme = nme;
MPI_Allreduce(&bnme, &ntotal, 1, MPI_LMP_BIGINT, MPI_SUM, world);
bigint atomOffset; // sum of all atoms on processes 0..me-1
MPI_Scan(&bnme, &atomOffset, 1, MPI_LMP_BIGINT, MPI_SUM, world);
atomOffset -= nme; // exclusive prefix sum needed
// Now we know the global size and the local subset size and offset
// of the atoms table
auto nAtomsGlobal = static_cast<size_t>(ntotal);
auto startRow = static_cast<size_t>(atomOffset);
auto nAtomsLocal = static_cast<size_t>(nme);
auto nColumns = static_cast<size_t>(size_one);
internal->varAtoms.SetShape({nAtomsGlobal, nColumns});
internal->varAtoms.SetSelection({{startRow, 0}, {nAtomsLocal, nColumns}});
// insure buf is sized for packing
// adios does not limit per-process data size so nme*size_one is not
// constrained to int
// if sorting on IDs also request ID list from pack()
// sort buf as needed
if (nme > maxbuf) {
maxbuf = nme;
memory->destroy(buf);
memory->create(buf, (maxbuf * size_one), "dump:buf");
}
if (sort_flag && sortcol == 0 && nme > maxids) {
maxids = nme;
memory->destroy(ids);
memory->create(ids, maxids, "dump:ids");
}
if (sort_flag && sortcol == 0)
pack(ids);
else
pack(nullptr);
if (sort_flag) sort();
openfile();
internal->fh.BeginStep();
// write info on data as scalars (by me==0)
if (me == 0) {
internal->fh.Put<uint64_t>("ntimestep", update->ntimestep);
internal->fh.Put<int>("nprocs", nprocs);
internal->fh.Put<double>("boxxlo", boxxlo);
internal->fh.Put<double>("boxxhi", boxxhi);
internal->fh.Put<double>("boxylo", boxylo);
internal->fh.Put<double>("boxyhi", boxyhi);
internal->fh.Put<double>("boxzlo", boxzlo);
internal->fh.Put<double>("boxzhi", boxzhi);
if (domain->triclinic) {
internal->fh.Put<double>("boxxy", boxxy);
internal->fh.Put<double>("boxxz", boxxz);
internal->fh.Put<double>("boxyz", boxyz);
}
}
// Everyone needs to write scalar variables that are used as dimensions and
// offsets of arrays
internal->fh.Put<uint64_t>("natoms", ntotal);
internal->fh.Put<int>("ncolumns", size_one);
internal->fh.Put<uint64_t>("nme", bnme);
internal->fh.Put<uint64_t>("offset", atomOffset);
// now write the atoms
internal->fh.Put<double>(internal->varAtoms, buf);
internal->fh.EndStep(); // I/O will happen now...
// nme = # of dump lines this proc contributes to dump
nme = count();
// ntotal = total # of atoms in snapshot
// atomOffset = sum of # of atoms up to this proc (exclusive prefix sum)
bigint bnme = nme;
MPI_Allreduce(&bnme, &ntotal, 1, MPI_LMP_BIGINT, MPI_SUM, world);
bigint atomOffset; // sum of all atoms on processes 0..me-1
MPI_Scan(&bnme, &atomOffset, 1, MPI_LMP_BIGINT, MPI_SUM, world);
atomOffset -= nme; // exclusive prefix sum needed
// Now we know the global size and the local subset size and offset
// of the atoms table
size_t nAtomsGlobal = static_cast<size_t>(ntotal);
size_t startRow = static_cast<size_t>(atomOffset);
size_t nAtomsLocal = static_cast<size_t>(nme);
size_t nColumns = static_cast<size_t>(size_one);
internal->varAtoms.SetShape({nAtomsGlobal, nColumns});
internal->varAtoms.SetSelection({{startRow, 0}, {nAtomsLocal, nColumns}});
// insure buf is sized for packing
// adios does not limit per-process data size so nme*size_one is not
// constrained to int
// if sorting on IDs also request ID list from pack()
// sort buf as needed
if (nme > maxbuf) {
maxbuf = nme;
memory->destroy(buf);
memory->create(buf, (maxbuf * size_one), "dump:buf");
}
if (sort_flag && sortcol == 0 && nme > maxids) {
maxids = nme;
memory->destroy(ids);
memory->create(ids, maxids, "dump:ids");
}
if (sort_flag && sortcol == 0)
pack(ids);
else
pack(nullptr);
if (sort_flag)
sort();
openfile();
internal->fh.BeginStep();
// write info on data as scalars (by me==0)
if (me == 0) {
internal->fh.Put<uint64_t>("ntimestep", update->ntimestep);
internal->fh.Put<int>("nprocs", nprocs);
internal->fh.Put<double>("boxxlo", boxxlo);
internal->fh.Put<double>("boxxhi", boxxhi);
internal->fh.Put<double>("boxylo", boxylo);
internal->fh.Put<double>("boxyhi", boxyhi);
internal->fh.Put<double>("boxzlo", boxzlo);
internal->fh.Put<double>("boxzhi", boxzhi);
if (domain->triclinic) {
internal->fh.Put<double>("boxxy", boxxy);
internal->fh.Put<double>("boxxz", boxxz);
internal->fh.Put<double>("boxyz", boxyz);
}
}
// Everyone needs to write scalar variables that are used as dimensions and
// offsets of arrays
internal->fh.Put<uint64_t>("natoms", ntotal);
internal->fh.Put<int>("ncolumns", size_one);
internal->fh.Put<uint64_t>("nme", bnme);
internal->fh.Put<uint64_t>("offset", atomOffset);
// now write the atoms
internal->fh.Put<double>(internal->varAtoms, buf);
internal->fh.EndStep(); // I/O will happen now...
if (multifile) {
internal->fh.Close();
}
if (multifile) internal->fh.Close();
}
/* ---------------------------------------------------------------------- */
void DumpAtomADIOS::init_style()
{
if (image_flag == 0)
size_one = 5;
else
size_one = 8;
if (image_flag == 0)
size_one = 5;
else
size_one = 8;
// setup boundary string
// setup boundary string
domain->boundary_string(boundstr);
domain->boundary_string(boundstr);
// remove % from filename since ADIOS always writes a global file with
// data/metadata
int len = strlen(filename);
char *ptr = strchr(filename, '%');
if (ptr) {
*ptr = '\0';
char *s = new char[len - 1];
snprintf(s, sizeof(s), "%s%s", filename, ptr + 1);
strncpy(filename, s, len);
// remove % from filename since ADIOS always writes a global file with
// data/metadata.
char *ptr = strchr(filename, '%');
if (ptr) {
while (*ptr) {
ptr[0] = ptr[1];
++ptr;
}
}
// setup column string
// setup column string
std::vector<std::string> columnNames;
std::vector<std::string> columnNames;
if (scale_flag == 0 && image_flag == 0) {
columns = (char *)"id type x y z";
columnNames = {"id", "type", "x", "y", "z"};
} else if (scale_flag == 0 && image_flag == 1) {
columns = (char *)"id type x y z ix iy iz";
columnNames = {"id", "type", "x", "y", "z", "ix", "iy", "iz"};
} else if (scale_flag == 1 && image_flag == 0) {
columns = (char *)"id type xs ys zs";
columnNames = {"id", "type", "xs", "ys", "zs"};
} else if (scale_flag == 1 && image_flag == 1) {
columns = (char *)"id type xs ys zs ix iy iz";
columnNames = {"id", "type", "xs", "ys", "zs", "ix", "iy", "iz"};
}
if (scale_flag == 0 && image_flag == 0) {
columns = (char *) "id type x y z";
columnNames = {"id", "type", "x", "y", "z"};
} else if (scale_flag == 0 && image_flag == 1) {
columns = (char *) "id type x y z ix iy iz";
columnNames = {"id", "type", "x", "y", "z", "ix", "iy", "iz"};
} else if (scale_flag == 1 && image_flag == 0) {
columns = (char *) "id type xs ys zs";
columnNames = {"id", "type", "xs", "ys", "zs"};
} else if (scale_flag == 1 && image_flag == 1) {
columns = (char *) "id type xs ys zs ix iy iz";
columnNames = {"id", "type", "xs", "ys", "zs", "ix", "iy", "iz"};
}
// setup function ptrs
for (int icol = 0; icol < (int) columnNames.size(); ++icol)
if (keyword_user[icol].size()) columnNames[icol] = keyword_user[icol];
if (scale_flag == 1 && image_flag == 0 && domain->triclinic == 0)
pack_choice = &DumpAtomADIOS::pack_scale_noimage;
else if (scale_flag == 1 && image_flag == 1 && domain->triclinic == 0)
pack_choice = &DumpAtomADIOS::pack_scale_image;
else if (scale_flag == 1 && image_flag == 0 && domain->triclinic == 1)
pack_choice = &DumpAtomADIOS::pack_scale_noimage_triclinic;
else if (scale_flag == 1 && image_flag == 1 && domain->triclinic == 1)
pack_choice = &DumpAtomADIOS::pack_scale_image_triclinic;
else if (scale_flag == 0 && image_flag == 0)
pack_choice = &DumpAtomADIOS::pack_noscale_noimage;
else if (scale_flag == 0 && image_flag == 1)
pack_choice = &DumpAtomADIOS::pack_noscale_image;
// setup function ptrs
/* Define the group of variables for the atom style here since it's a fixed
if (scale_flag == 1 && image_flag == 0 && domain->triclinic == 0)
pack_choice = &DumpAtomADIOS::pack_scale_noimage;
else if (scale_flag == 1 && image_flag == 1 && domain->triclinic == 0)
pack_choice = &DumpAtomADIOS::pack_scale_image;
else if (scale_flag == 1 && image_flag == 0 && domain->triclinic == 1)
pack_choice = &DumpAtomADIOS::pack_scale_noimage_triclinic;
else if (scale_flag == 1 && image_flag == 1 && domain->triclinic == 1)
pack_choice = &DumpAtomADIOS::pack_scale_image_triclinic;
else if (scale_flag == 0 && image_flag == 0)
pack_choice = &DumpAtomADIOS::pack_noscale_noimage;
else if (scale_flag == 0 && image_flag == 1)
pack_choice = &DumpAtomADIOS::pack_noscale_image;
/* Define the group of variables for the atom style here since it's a fixed
* set */
internal->io = internal->ad->DeclareIO(internal->ioName);
if (!internal->io.InConfigFile()) {
// if not defined by user, we can change the default settings
// BPFile is the default writer
internal->io.SetEngine("BPFile");
int num_aggregators = multiproc;
if (num_aggregators == 0)
num_aggregators = 1;
char nstreams[128];
snprintf(nstreams, sizeof(nstreams), "%d", num_aggregators);
internal->io.SetParameters({{"substreams", nstreams}});
if (me == 0 && screen)
fprintf(
screen,
"ADIOS method for %s is n-to-m (aggregation with %s writers)\n",
filename, nstreams);
}
internal->io = internal->ad->DeclareIO(internal->ioName);
if (!internal->io.InConfigFile()) {
// if not defined by user, we can change the default settings
// BPFile is the default writer
internal->io.SetEngine("BPFile");
int num_aggregators = multiproc;
if (num_aggregators == 0) num_aggregators = 1;
auto nstreams = std::to_string(num_aggregators);
internal->io.SetParameters({{"substreams", nstreams}});
if (me == 0)
utils::logmesg(lmp, "ADIOS method for {} is n-to-m (aggregation with {} writers)\n", filename,
nstreams);
}
internal->io.DefineVariable<uint64_t>("ntimestep");
internal->io.DefineVariable<uint64_t>("natoms");
internal->io.DefineVariable<uint64_t>("ntimestep");
internal->io.DefineVariable<uint64_t>("natoms");
internal->io.DefineVariable<int>("nprocs");
internal->io.DefineVariable<int>("ncolumns");
internal->io.DefineVariable<int>("nprocs");
internal->io.DefineVariable<int>("ncolumns");
internal->io.DefineVariable<double>("boxxlo");
internal->io.DefineVariable<double>("boxxhi");
internal->io.DefineVariable<double>("boxylo");
internal->io.DefineVariable<double>("boxyhi");
internal->io.DefineVariable<double>("boxzlo");
internal->io.DefineVariable<double>("boxzhi");
internal->io.DefineVariable<double>("boxxlo");
internal->io.DefineVariable<double>("boxxhi");
internal->io.DefineVariable<double>("boxylo");
internal->io.DefineVariable<double>("boxyhi");
internal->io.DefineVariable<double>("boxzlo");
internal->io.DefineVariable<double>("boxzhi");
internal->io.DefineVariable<double>("boxxy");
internal->io.DefineVariable<double>("boxxz");
internal->io.DefineVariable<double>("boxyz");
internal->io.DefineVariable<double>("boxxy");
internal->io.DefineVariable<double>("boxxz");
internal->io.DefineVariable<double>("boxyz");
internal->io.DefineAttribute<int>("triclinic", domain->triclinic);
internal->io.DefineAttribute<int>("scaled", scale_flag);
internal->io.DefineAttribute<int>("image", image_flag);
internal->io.DefineAttribute<int>("triclinic", domain->triclinic);
internal->io.DefineAttribute<int>("scaled", scale_flag);
internal->io.DefineAttribute<int>("image", image_flag);
int *boundaryptr = reinterpret_cast<int *>(domain->boundary);
internal->io.DefineAttribute<int>("boundary", boundaryptr, 6);
int *boundaryptr = reinterpret_cast<int *>(domain->boundary);
internal->io.DefineAttribute<int>("boundary", boundaryptr, 6);
size_t nColumns = static_cast<size_t>(size_one);
internal->io.DefineAttribute<std::string>("columns", columnNames.data(),
nColumns);
internal->io.DefineAttribute<std::string>("columnstr", columns);
internal->io.DefineAttribute<std::string>("boundarystr", boundstr);
internal->io.DefineAttribute<std::string>("LAMMPS/dump_style", "atom");
internal->io.DefineAttribute<std::string>("LAMMPS/version",
lmp->version);
internal->io.DefineAttribute<std::string>("LAMMPS/num_ver",
std::to_string(lmp->num_ver));
auto nColumns = static_cast<size_t>(size_one);
internal->io.DefineAttribute<std::string>("columns", columnNames.data(), nColumns);
internal->io.DefineAttribute<std::string>("columnstr", columns);
internal->io.DefineAttribute<std::string>("boundarystr", boundstr);
internal->io.DefineAttribute<std::string>("LAMMPS/dump_style", "atom");
internal->io.DefineAttribute<std::string>("LAMMPS/version", lmp->version);
internal->io.DefineAttribute<std::string>("LAMMPS/num_ver", std::to_string(lmp->num_ver));
internal->io.DefineVariable<uint64_t>(
"nme", {adios2::LocalValueDim}); // local dimension variable
internal->io.DefineVariable<uint64_t>(
"offset", {adios2::LocalValueDim}); // local dimension variable
// local dimension variables
internal->io.DefineVariable<uint64_t>("nme", {adios2::LocalValueDim});
internal->io.DefineVariable<uint64_t>("offset", {adios2::LocalValueDim});
// atom table size is not known at the moment
// it will be correctly defined at the moment of write
size_t UnknownSizeYet = 1;
internal->varAtoms = internal->io.DefineVariable<double>(
"atoms", {UnknownSizeYet, nColumns}, {UnknownSizeYet, 0},
{UnknownSizeYet, nColumns});
// atom table size is not known at the moment
// it will be correctly defined at the moment of write
size_t UnknownSizeYet = 1;
internal->varAtoms = internal->io.DefineVariable<double>(
"atoms", {UnknownSizeYet, nColumns}, {UnknownSizeYet, 0}, {UnknownSizeYet, nColumns});
}

22
src/ADIOS/dump_atom_adios.h
View File

@ -30,12 +30,12 @@ class DumpAtomADIOS : public DumpAtom {
public:
DumpAtomADIOS(class LAMMPS *, int, char **);
virtual ~DumpAtomADIOS();
~DumpAtomADIOS() override;
protected:
virtual void openfile();
virtual void write();
virtual void init_style();
void openfile() override;
void write() override;
void init_style() override;
private:
DumpAtomADIOSInternal *internal;
@ -44,17 +44,3 @@ class DumpAtomADIOS : public DumpAtom {
#endif
#endif
/* ERROR/WARNING messages:
E: Cannot open dump file %s
The output file for the dump command cannot be opened. Check that the
path and name are correct.
E: Too much per-proc info for dump
Number of local atoms times number of columns must fit in a 32-bit
integer for dump.
*/

602
src/ADIOS/dump_custom_adios.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -17,427 +16,320 @@
------------------------------------------------------------------------- */
#include "dump_custom_adios.h"
#include "atom.h"
#include "compute.h"
#include "domain.h"
#include "error.h"
#include "fix.h"
#include "force.h"
#include "group.h"
#include "input.h"
#include "memory.h"
#include "modify.h"
#include "region.h"
#include "universe.h"
#include "update.h"
#include "variable.h"
#include <cmath>
#include <cstring>
#include "adios2.h"
#include "adios_common.h"
using namespace LAMMPS_NS;
#define MAX_TEXT_HEADER_SIZE 4096
#define DUMP_BUF_CHUNK_SIZE 16384
#define DUMP_BUF_INCREMENT_SIZE 4096
namespace LAMMPS_NS {
class DumpCustomADIOSInternal {
enum {
ID,
MOL,
TYPE,
ELEMENT,
MASS,
X,
Y,
Z,
XS,
YS,
ZS,
XSTRI,
YSTRI,
ZSTRI,
XU,
YU,
ZU,
XUTRI,
YUTRI,
ZUTRI,
XSU,
YSU,
ZSU,
XSUTRI,
YSUTRI,
ZSUTRI,
IX,
IY,
IZ,
VX,
VY,
VZ,
FX,
FY,
FZ,
Q,
MUX,
MUY,
MUZ,
MU,
RADIUS,
DIAMETER,
OMEGAX,
OMEGAY,
OMEGAZ,
ANGMOMX,
ANGMOMY,
ANGMOMZ,
TQX,
TQY,
TQZ,
SPIN,
ERADIUS,
ERVEL,
ERFORCE,
COMPUTE,
FIX,
VARIABLE
public:
DumpCustomADIOSInternal(){};
~DumpCustomADIOSInternal() = default;
// name of adios group, referrable in adios2_config.xml
const std::string ioName = "custom";
adios2::ADIOS *ad = nullptr; // adios object
adios2::IO io; // adios group of variables and attributes in this dump
adios2::Engine fh; // adios file/stream handle object
// one ADIOS output variable we need to change every step
adios2::Variable<double> varAtoms;
// list of column names for the atom table
// (individual list of 'columns' string)
std::vector<std::string> columnNames;
};
enum { LT, LE, GT, GE, EQ, NEQ };
enum { INT, DOUBLE, STRING, BIGINT }; // same as in DumpCustom
namespace LAMMPS_NS
{
class DumpCustomADIOSInternal
{
public:
DumpCustomADIOSInternal() {};
~DumpCustomADIOSInternal() = default;
// name of adios group, referrable in adios2_config.xml
const std::string ioName = "custom";
adios2::ADIOS *ad = nullptr; // adios object
adios2::IO io; // adios group of variables and attributes in this dump
adios2::Engine fh; // adios file/stream handle object
// one ADIOS output variable we need to change every step
adios2::Variable<double> varAtoms;
// list of column names for the atom table
// (individual list of 'columns' string)
std::vector<std::string> columnNames;
};
} // namespace LAMMPS_NS
} // namespace LAMMPS_NS
/* ---------------------------------------------------------------------- */
DumpCustomADIOS::DumpCustomADIOS(LAMMPS *lmp, int narg, char **arg)
: DumpCustom(lmp, narg, arg)
DumpCustomADIOS::DumpCustomADIOS(LAMMPS *lmp, int narg, char **arg) : DumpCustom(lmp, narg, arg)
{
internal = new DumpCustomADIOSInternal();
try {
internal->ad =
new adios2::ADIOS("adios2_config.xml", world, adios2::DebugON);
} catch (std::ios_base::failure &e) {
char str[256];
snprintf(str, sizeof(str), "ADIOS initialization failed with error: %s",
e.what());
error->one(FLERR, str);
}
// create a default adios2_config.xml if it doesn't exist yet.
FILE *cfgfp = fopen("adios2_config.xml", "r");
if (!cfgfp) {
cfgfp = fopen("adios2_config.xml", "w");
if (cfgfp) fputs(default_config, cfgfp);
}
if (cfgfp) fclose(cfgfp);
// if (screen) fprintf(screen, "DumpCustomADIOS constructor: nvariable=%d
// id_variable=%p, variables=%p, nfield=%d, earg=%p\n", nvariable,
// id_variable, variable, nfield, earg);
internal->columnNames.reserve(nfield);
for (int i = 0; i < nfield; ++i) {
internal->columnNames.push_back(earg[i]);
// if (screen) fprintf(screen, "earg[%d] = '%s'\n", i, earg[i]);
}
internal = new DumpCustomADIOSInternal();
try {
internal->ad = new adios2::ADIOS("adios2_config.xml", world, adios2::DebugON);
} catch (std::ios_base::failure &e) {
error->all(FLERR, "ADIOS initialization failed with error: {}", e.what());
}
internal->columnNames.reserve(nfield);
for (int i = 0; i < nfield; ++i) { internal->columnNames.push_back(earg[i]); }
}
/* ---------------------------------------------------------------------- */
DumpCustomADIOS::~DumpCustomADIOS()
{
internal->columnNames.clear();
if (internal->fh) {
internal->fh.Close();
}
delete internal->ad;
delete internal;
internal->columnNames.clear();
if (internal->fh) { internal->fh.Close(); }
delete internal->ad;
delete internal;
}
/* ---------------------------------------------------------------------- */
void DumpCustomADIOS::openfile()
{
if (multifile) {
// if one file per timestep, replace '*' with current timestep
char *filestar = strdup(filename);
char *filecurrent = new char[strlen(filestar) + 16];
char *ptr = strchr(filestar, '*');
*ptr = '\0';
if (padflag == 0)
sprintf(filecurrent, "%s" BIGINT_FORMAT "%s", filestar,
update->ntimestep, ptr + 1);
else {
char bif[8], pad[16];
strcpy(bif, BIGINT_FORMAT);
sprintf(pad, "%%s%%0%d%s%%s", padflag, &bif[1]);
sprintf(filecurrent, pad, filestar, update->ntimestep, ptr + 1);
}
internal->fh =
internal->io.Open(filecurrent, adios2::Mode::Write, world);
if (!internal->fh) {
char str[128];
sprintf(str, "Cannot open dump file %s", filecurrent);
error->one(FLERR, str);
}
free(filestar);
delete[] filecurrent;
} else {
if (!singlefile_opened) {
internal->fh =
internal->io.Open(filename, adios2::Mode::Write, world);
if (!internal->fh) {
char str[128];
sprintf(str, "Cannot open dump file %s", filename);
error->one(FLERR, str);
}
singlefile_opened = 1;
}
if (multifile) {
// if one file per timestep, replace '*' with current timestep
auto filecurrent = utils::star_subst(filename, update->ntimestep, padflag);
internal->fh = internal->io.Open(filecurrent, adios2::Mode::Write, world);
if (!internal->fh) error->one(FLERR, "Cannot open dump file {}", filecurrent);
} else {
if (!singlefile_opened) {
internal->fh = internal->io.Open(filename, adios2::Mode::Write, world);
if (!internal->fh) error->one(FLERR, "Cannot open dump file {}", filename);
singlefile_opened = 1;
}
}
}
/* ---------------------------------------------------------------------- */
void DumpCustomADIOS::write()
{
if (domain->triclinic == 0) {
boxxlo = domain->boxlo[0];
boxxhi = domain->boxhi[0];
boxylo = domain->boxlo[1];
boxyhi = domain->boxhi[1];
boxzlo = domain->boxlo[2];
boxzhi = domain->boxhi[2];
} else {
boxxlo = domain->boxlo_bound[0];
boxxhi = domain->boxhi_bound[0];
boxylo = domain->boxlo_bound[1];
boxyhi = domain->boxhi_bound[1];
boxzlo = domain->boxlo_bound[2];
boxzhi = domain->boxhi_bound[2];
boxxy = domain->xy;
boxxz = domain->xz;
boxyz = domain->yz;
if (domain->triclinic == 0) {
boxxlo = domain->boxlo[0];
boxxhi = domain->boxhi[0];
boxylo = domain->boxlo[1];
boxyhi = domain->boxhi[1];
boxzlo = domain->boxlo[2];
boxzhi = domain->boxhi[2];
} else {
boxxlo = domain->boxlo_bound[0];
boxxhi = domain->boxhi_bound[0];
boxylo = domain->boxlo_bound[1];
boxyhi = domain->boxhi_bound[1];
boxzlo = domain->boxlo_bound[2];
boxzhi = domain->boxhi_bound[2];
boxxy = domain->xy;
boxxz = domain->xz;
boxyz = domain->yz;
}
// nme = # of dump lines this proc contributes to dump
nme = count();
// ntotal = total # of atoms in snapshot
// atomOffset = sum of # of atoms up to this proc (exclusive prefix sum)
bigint bnme = nme;
MPI_Allreduce(&bnme, &ntotal, 1, MPI_LMP_BIGINT, MPI_SUM, world);
bigint atomOffset; // sum of all atoms on processes 0..me-1
MPI_Scan(&bnme, &atomOffset, 1, MPI_LMP_BIGINT, MPI_SUM, world);
atomOffset -= nme; // exclusive prefix sum needed
// Now we know the global size and the local subset size and offset
// of the atoms table
auto nAtomsGlobal = static_cast<size_t>(ntotal);
auto startRow = static_cast<size_t>(atomOffset);
auto nAtomsLocal = static_cast<size_t>(nme);
auto nColumns = static_cast<size_t>(size_one);
internal->varAtoms.SetShape({nAtomsGlobal, nColumns});
internal->varAtoms.SetSelection({{startRow, 0}, {nAtomsLocal, nColumns}});
// insure filewriter proc can receive everyone's info
// limit nmax*size_one to int since used as arg in MPI_Rsend() below
// pack my data into buf
// if sorting on IDs also request ID list from pack()
// sort buf as needed
if (nme > maxbuf) {
if ((bigint) nme * size_one > MAXSMALLINT) error->all(FLERR, "Too much per-proc info for dump");
maxbuf = nme;
memory->destroy(buf);
memory->create(buf, (maxbuf * size_one), "dump:buf");
}
if (sort_flag && sortcol == 0 && nme > maxids) {
maxids = nme;
memory->destroy(ids);
memory->create(ids, maxids, "dump:ids");
}
if (sort_flag && sortcol == 0)
pack(ids);
else
pack(nullptr);
if (sort_flag) sort();
openfile();
internal->fh.BeginStep();
// write info on data as scalars (by me==0)
if (me == 0) {
internal->fh.Put<uint64_t>("ntimestep", update->ntimestep);
internal->fh.Put<int>("nprocs", nprocs);
internal->fh.Put<double>("boxxlo", boxxlo);
internal->fh.Put<double>("boxxhi", boxxhi);
internal->fh.Put<double>("boxylo", boxylo);
internal->fh.Put<double>("boxyhi", boxyhi);
internal->fh.Put<double>("boxzlo", boxzlo);
internal->fh.Put<double>("boxzhi", boxzhi);
if (domain->triclinic) {
internal->fh.Put<double>("boxxy", boxxy);
internal->fh.Put<double>("boxxz", boxxz);
internal->fh.Put<double>("boxyz", boxyz);
}
}
// Everyone needs to write scalar variables that are used as dimensions and
// offsets of arrays
internal->fh.Put<uint64_t>("natoms", ntotal);
internal->fh.Put<int>("ncolumns", size_one);
internal->fh.Put<uint64_t>("nme", bnme);
internal->fh.Put<uint64_t>("offset", atomOffset);
// now write the atoms
internal->fh.Put<double>("atoms", buf);
internal->fh.EndStep(); // I/O will happen now...
// nme = # of dump lines this proc contributes to dump
nme = count();
// ntotal = total # of atoms in snapshot
// atomOffset = sum of # of atoms up to this proc (exclusive prefix sum)
bigint bnme = nme;
MPI_Allreduce(&bnme, &ntotal, 1, MPI_LMP_BIGINT, MPI_SUM, world);
bigint atomOffset; // sum of all atoms on processes 0..me-1
MPI_Scan(&bnme, &atomOffset, 1, MPI_LMP_BIGINT, MPI_SUM, world);
atomOffset -= nme; // exclusive prefix sum needed
// Now we know the global size and the local subset size and offset
// of the atoms table
size_t nAtomsGlobal = static_cast<size_t>(ntotal);
size_t startRow = static_cast<size_t>(atomOffset);
size_t nAtomsLocal = static_cast<size_t>(nme);
size_t nColumns = static_cast<size_t>(size_one);
internal->varAtoms.SetShape({nAtomsGlobal, nColumns});
internal->varAtoms.SetSelection({{startRow, 0}, {nAtomsLocal, nColumns}});
// insure filewriter proc can receive everyone's info
// limit nmax*size_one to int since used as arg in MPI_Rsend() below
// pack my data into buf
// if sorting on IDs also request ID list from pack()
// sort buf as needed
if (nme > maxbuf) {
if ((bigint)nme * size_one > MAXSMALLINT)
error->all(FLERR, "Too much per-proc info for dump");
maxbuf = nme;
memory->destroy(buf);
memory->create(buf, (maxbuf * size_one), "dump:buf");
}
if (sort_flag && sortcol == 0 && nme > maxids) {
maxids = nme;
memory->destroy(ids);
memory->create(ids, maxids, "dump:ids");
}
if (sort_flag && sortcol == 0)
pack(ids);
else
pack(nullptr);
if (sort_flag)
sort();
openfile();
internal->fh.BeginStep();
// write info on data as scalars (by me==0)
if (me == 0) {
internal->fh.Put<uint64_t>("ntimestep", update->ntimestep);
internal->fh.Put<int>("nprocs", nprocs);
internal->fh.Put<double>("boxxlo", boxxlo);
internal->fh.Put<double>("boxxhi", boxxhi);
internal->fh.Put<double>("boxylo", boxylo);
internal->fh.Put<double>("boxyhi", boxyhi);
internal->fh.Put<double>("boxzlo", boxzlo);
internal->fh.Put<double>("boxzhi", boxzhi);
if (domain->triclinic) {
internal->fh.Put<double>("boxxy", boxxy);
internal->fh.Put<double>("boxxz", boxxz);
internal->fh.Put<double>("boxyz", boxyz);
}
}
// Everyone needs to write scalar variables that are used as dimensions and
// offsets of arrays
internal->fh.Put<uint64_t>("natoms", ntotal);
internal->fh.Put<int>("ncolumns", size_one);
internal->fh.Put<uint64_t>("nme", bnme);
internal->fh.Put<uint64_t>("offset", atomOffset);
// now write the atoms
internal->fh.Put<double>("atoms", buf);
internal->fh.EndStep(); // I/O will happen now...
if (multifile) {
internal->fh.Close();
}
if (multifile) { internal->fh.Close(); }
}
/* ---------------------------------------------------------------------- */
void DumpCustomADIOS::init_style()
{
// assemble column string from defaults and user values
// setup boundary string
delete[] columns;
std::string combined;
int icol = 0;
for (auto item : utils::split_words(columns_default)) {
if (combined.size()) combined += " ";
if (keyword_user[icol].size())
combined += keyword_user[icol];
else
combined += item;
++icol;
}
columns = utils::strdup(combined);
domain->boundary_string(boundstr);
// setup boundary string
// remove % from filename since ADIOS always writes a global file with
// data/metadata
int len = strlen(filename);
char *ptr = strchr(filename, '%');
if (ptr) {
*ptr = '\0';
char *s = new char[len - 1];
sprintf(s, "%s%s", filename, ptr + 1);
strncpy(filename, s, len);
domain->boundary_string(boundstr);
// remove % from filename since ADIOS always writes a global file with
// data/metadata
char *ptr = strchr(filename, '%');
if (ptr) {
while (*ptr) {
ptr[0] = ptr[1];
++ptr;
}
}
/* The next four loops are copied from dump_custom_mpiio, but nothing is
* done with them.
* It is unclear why we need them here.
* For metadata, variable[] will be written out as an ADIOS attribute if
* nvariable>0
*/
// find current ptr for each compute,fix,variable
// check that fix frequency is acceptable
int icompute;
for (int i = 0; i < ncompute; i++) {
icompute = modify->find_compute(id_compute[i]);
if (icompute < 0)
error->all(FLERR, "Could not find dump custom compute ID");
compute[i] = modify->compute[icompute];
}
/* The next four loops are copied from dump_custom_mpiio, but nothing is
* done with them.
* It is unclear why we need them here.
* For metadata, variable[] will be written out as an ADIOS attribute if
* nvariable>0
*/
// find current ptr for each compute,fix,variable
// check that fix frequency is acceptable
for (int i = 0; i < ncompute; i++) {
compute[i] = modify->get_compute_by_id(id_compute[i]);
if (!compute[i])
error->all(FLERR, "Could not find dump custom/adios compute ID {}", id_compute[i]);
}
int ifix;
for (int i = 0; i < nfix; i++) {
ifix = modify->find_fix(id_fix[i]);
if (ifix < 0)
error->all(FLERR, "Could not find dump custom fix ID");
fix[i] = modify->fix[ifix];
if (nevery % modify->fix[ifix]->peratom_freq)
error->all(FLERR,
"Dump custom and fix not computed at compatible times");
}
for (int i = 0; i < nfix; i++) {
fix[i] = modify->get_fix_by_id(id_fix[i]);
if (!fix[i]) error->all(FLERR, "Could not find dump custom/adios fix ID {}", id_fix[i]);
if (nevery % fix[i]->peratom_freq)
error->all(FLERR, "dump custom/adios and fix {} with ID {} not computed at compatible times",
fix[i]->style, id_fix[i]);
}
int ivariable;
for (int i = 0; i < nvariable; i++) {
ivariable = input->variable->find(id_variable[i]);
if (ivariable < 0)
error->all(FLERR, "Could not find dump custom variable name");
variable[i] = ivariable;
}
int ivariable;
for (int i = 0; i < nvariable; i++) {
ivariable = input->variable->find(id_variable[i]);
if (ivariable < 0) error->all(FLERR, "Could not find dump custom/adios variable name");
variable[i] = ivariable;
}
// set index and check validity of region
if (iregion >= 0) {
iregion = domain->find_region(idregion);
if (iregion == -1)
error->all(FLERR, "Region ID for dump custom does not exist");
}
// set index and check validity of region
if (idregion && !domain->get_region_by_id(idregion))
error->all(FLERR, "Region {} for dump custom/adios does not exist", idregion);
/* Define the group of variables for the atom style here since it's a fixed
* set */
internal->io = internal->ad->DeclareIO(internal->ioName);
if (!internal->io.InConfigFile()) {
// if not defined by user, we can change the default settings
// BPFile is the default writer
internal->io.SetEngine("BPFile");
int num_aggregators = multiproc;
if (num_aggregators == 0)
num_aggregators = 1;
char nstreams[128];
sprintf(nstreams, "%d", num_aggregators);
internal->io.SetParameters({{"substreams", nstreams}});
if (me == 0 && screen)
fprintf(
screen,
"ADIOS method for %s is n-to-m (aggregation with %s writers)\n",
filename, nstreams);
}
/* Define the group of variables for the atom style here since it's a fixed
* set */
internal->io = internal->ad->DeclareIO(internal->ioName);
if (!internal->io.InConfigFile()) {
// if not defined by user, we can change the default settings
// BPFile is the default writer
internal->io.SetEngine("BPFile");
int num_aggregators = multiproc;
if (num_aggregators == 0) num_aggregators = 1;
auto nstreams = std::to_string(num_aggregators);
internal->io.SetParameters({{"substreams", nstreams}});
if (me == 0)
utils::logmesg(lmp, "ADIOS method for {} is n-to-m (aggregation with {} writers)\n", filename,
nstreams);
}
internal->io.DefineVariable<uint64_t>("ntimestep");
internal->io.DefineVariable<uint64_t>("natoms");
internal->io.DefineVariable<uint64_t>("ntimestep");
internal->io.DefineVariable<uint64_t>("natoms");
internal->io.DefineVariable<int>("nprocs");
internal->io.DefineVariable<int>("ncolumns");
internal->io.DefineVariable<int>("nprocs");
internal->io.DefineVariable<int>("ncolumns");
internal->io.DefineVariable<double>("boxxlo");
internal->io.DefineVariable<double>("boxxhi");
internal->io.DefineVariable<double>("boxylo");
internal->io.DefineVariable<double>("boxyhi");
internal->io.DefineVariable<double>("boxzlo");
internal->io.DefineVariable<double>("boxzhi");
internal->io.DefineVariable<double>("boxxlo");
internal->io.DefineVariable<double>("boxxhi");
internal->io.DefineVariable<double>("boxylo");
internal->io.DefineVariable<double>("boxyhi");
internal->io.DefineVariable<double>("boxzlo");
internal->io.DefineVariable<double>("boxzhi");
internal->io.DefineVariable<double>("boxxy");
internal->io.DefineVariable<double>("boxxz");
internal->io.DefineVariable<double>("boxyz");
internal->io.DefineVariable<double>("boxxy");
internal->io.DefineVariable<double>("boxxz");
internal->io.DefineVariable<double>("boxyz");
internal->io.DefineAttribute<int>("triclinic", domain->triclinic);
internal->io.DefineAttribute<int>("triclinic", domain->triclinic);
int *boundaryptr = reinterpret_cast<int *>(domain->boundary);
internal->io.DefineAttribute<int>("boundary", boundaryptr, 6);
int *boundaryptr = reinterpret_cast<int *>(domain->boundary);
internal->io.DefineAttribute<int>("boundary", boundaryptr, 6);
size_t nColumns = static_cast<size_t>(size_one);
internal->io.DefineAttribute<std::string>(
"columns", internal->columnNames.data(), nColumns);
internal->io.DefineAttribute<std::string>("columnstr", columns);
internal->io.DefineAttribute<std::string>("boundarystr", boundstr);
internal->io.DefineAttribute<std::string>("LAMMPS/dump_style", "custom");
internal->io.DefineAttribute<std::string>("LAMMPS/version",
lmp->version);
internal->io.DefineAttribute<std::string>("LAMMPS/num_ver",
std::to_string(lmp->num_ver));
auto nColumns = static_cast<size_t>(size_one);
internal->io.DefineAttribute<std::string>("columns", internal->columnNames.data(), nColumns);
internal->io.DefineAttribute<std::string>("columnstr", columns);
internal->io.DefineAttribute<std::string>("boundarystr", boundstr);
internal->io.DefineAttribute<std::string>("LAMMPS/dump_style", "custom");
internal->io.DefineAttribute<std::string>("LAMMPS/version", lmp->version);
internal->io.DefineAttribute<std::string>("LAMMPS/num_ver", std::to_string(lmp->num_ver));
internal->io.DefineVariable<uint64_t>(
"nme", {adios2::LocalValueDim}); // local dimension variable
internal->io.DefineVariable<uint64_t>(
"offset", {adios2::LocalValueDim}); // local dimension variable
internal->io.DefineVariable<uint64_t>("nme",
{adios2::LocalValueDim}); // local dimension variable
internal->io.DefineVariable<uint64_t>("offset",
{adios2::LocalValueDim}); // local dimension variable
// atom table size is not known at the moment
// it will be correctly defined at the moment of write
size_t UnknownSizeYet = 1;
internal->varAtoms = internal->io.DefineVariable<double>(
"atoms", {UnknownSizeYet, nColumns}, {UnknownSizeYet, 0},
{UnknownSizeYet, nColumns});
// atom table size is not known at the moment
// it will be correctly defined at the moment of write
size_t UnknownSizeYet = 1;
internal->varAtoms = internal->io.DefineVariable<double>(
"atoms", {UnknownSizeYet, nColumns}, {UnknownSizeYet, 0}, {UnknownSizeYet, nColumns});
}

48
src/ADIOS/dump_custom_adios.h
View File

@ -29,12 +29,12 @@ class DumpCustomADIOSInternal;
class DumpCustomADIOS : public DumpCustom {
public:
DumpCustomADIOS(class LAMMPS *, int, char **);
virtual ~DumpCustomADIOS();
~DumpCustomADIOS() override;
protected:
virtual void openfile();
virtual void write();
virtual void init_style();
void openfile() override;
void write() override;
void init_style() override;
private:
DumpCustomADIOSInternal *internal;
@ -43,43 +43,3 @@ class DumpCustomADIOS : public DumpCustom {
#endif
#endif
/* ERROR/WARNING messages:
E: Cannot open dump file %s
The output file for the dump command cannot be opened. Check that the
path and name are correct.
E: Too much per-proc info for dump
Number of local atoms times number of columns must fit in a 32-bit
integer for dump.
E: Dump_modify format string is too short
There are more fields to be dumped in a line of output than your
format string specifies.
E: Could not find dump custom compute ID
Self-explanatory.
E: Could not find dump custom fix ID
Self-explanatory.
E: Dump custom and fix not computed at compatible times
The fix must produce per-atom quantities on timesteps that dump custom
needs them.
E: Could not find dump custom variable name
Self-explanatory.
E: Region ID for dump custom does not exist
Self-explanatory.
*/

660
src/ADIOS/reader_adios.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -17,20 +16,20 @@
------------------------------------------------------------------------- */
#include "reader_adios.h"
#include "comm.h"
#include "error.h"
#include "memory.h"
#include <cmath>
#include <cstdlib>
#include <cstring>
#include "adios2.h"
#include "math_const.h"
#include "adios_common.h"
using namespace LAMMPS_NS;
using namespace MathConst;
// also in read_dump.cpp
enum { ID, TYPE, X, Y, Z, VX, VY, VZ, Q, IX, IY, IZ, FX, FY, FZ };
@ -38,78 +37,69 @@ enum { UNSET, NOSCALE_NOWRAP, NOSCALE_WRAP, SCALE_NOWRAP, SCALE_WRAP };
#define SMALL 1.0e-6
// true if the difference between two floats is "small".
// cannot use fabsf() since it is not fully portable.
static bool is_smalldiff(const float &val1, const float &val2)
{
return (fabs(static_cast<double>(val1 - val2)) < SMALL);
}
namespace LAMMPS_NS {
class ReadADIOSInternal {
namespace LAMMPS_NS
{
class ReadADIOSInternal
{
public:
ReadADIOSInternal(){};
~ReadADIOSInternal() = default;
public:
ReadADIOSInternal() {};
~ReadADIOSInternal() = default;
// name of adios group, referrable in adios2_config.xml
const std::string ioName = "read_dump";
adios2::ADIOS *ad = nullptr; // adios object
adios2::IO io; // adios group of variables and attributes in this dump
adios2::Engine fh; // adios file/stream handle object
// ADIOS input variables we need to change every step
adios2::Variable<uint64_t> varNtimestep;
adios2::Variable<uint64_t> varNatoms;
adios2::Variable<double> varAtoms;
// list of column names for the atom table
// (individual list of 'columns' string)
std::vector<std::string> columnNames;
float timeout = 0.0;
// name of adios group, referrable in adios2_config.xml
const std::string ioName = "read_dump";
adios2::ADIOS *ad = nullptr; // adios object
adios2::IO io; // adios group of variables and attributes in this dump
adios2::Engine fh; // adios file/stream handle object
// ADIOS input variables we need to change every step
adios2::Variable<uint64_t> varNtimestep;
adios2::Variable<uint64_t> varNatoms;
adios2::Variable<double> varAtoms;
// list of column names for the atom table
// (individual list of 'columns' string)
std::vector<std::string> columnNames;
float timeout = 0.0;
};
} // namespace LAMMPS_NS
} // namespace LAMMPS_NS
/* ---------------------------------------------------------------------- */
ReaderADIOS::ReaderADIOS(LAMMPS *lmp) : Reader(lmp)
{
fieldindex = nullptr;
nAtoms = 0;
nAtomsTotal = 0;
atomOffset = 0;
nstep = 0;
nid = 0;
me = comm->me;
fieldindex = nullptr;
nAtoms = 0;
nAtomsTotal = 0;
atomOffset = 0;
nstep = 0;
nid = 0;
me = comm->me;
internal = new ReadADIOSInternal();
try {
internal->ad =
new adios2::ADIOS("adios2_config.xml", world, adios2::DebugON);
} catch (std::ios_base::failure &e) {
char str[256];
snprintf(str, sizeof(str), "ADIOS initialization failed with error: %s",
e.what());
error->one(FLERR, str);
}
// create a default adios2_config.xml if it doesn't exist yet.
FILE *cfgfp = fopen("adios2_config.xml", "r");
if (!cfgfp) {
cfgfp = fopen("adios2_config.xml", "w");
if (cfgfp) fputs(default_config, cfgfp);
}
if (cfgfp) fclose(cfgfp);
/* Define the group holding all variables and attributes */
internal->io = internal->ad->DeclareIO(internal->ioName);
internal = new ReadADIOSInternal();
try {
internal->ad = new adios2::ADIOS("adios2_config.xml", world, adios2::DebugON);
} catch (std::ios_base::failure &e) {
error->one(FLERR, "ADIOS initialization failed with error: {}", e.what());
}
/* Define the group holding all variables and attributes */
internal->io = internal->ad->DeclareIO(internal->ioName);
}
/* ---------------------------------------------------------------------- */
ReaderADIOS::~ReaderADIOS()
{
if (me == 0) {
memory->destroy(fieldindex);
}
internal->columnNames.clear();
if (internal->fh) {
internal->fh.Close();
}
delete internal->ad;
delete internal;
if (me == 0) memory->destroy(fieldindex);
internal->columnNames.clear();
if (internal->fh) internal->fh.Close();
delete internal->ad;
delete internal;
}
/* ----------------------------------------------------------------------
@ -118,23 +108,19 @@ ReaderADIOS::~ReaderADIOS()
------------------------------------------------------------------------- */
void ReaderADIOS::settings(int narg, char **arg)
{
int idx = 0;
while (idx < narg) {
if (!strcmp(arg[idx], "timeout")) {
if (idx + 1 < narg) {
internal->timeout = std::stof(arg[idx + 1]);
internal->io.SetParameter("OpenTimeoutSecs", arg[idx + 1]);
++idx;
} else {
char str[128];
snprintf(str, sizeof(str),
"Missing value for 'timeout' option for ADIOS "
"read_dump command");
error->one(FLERR, str);
}
}
int idx = 0;
while (idx < narg) {
if (!strcmp(arg[idx], "timeout")) {
if (idx + 1 < narg) {
internal->timeout = std::stof(arg[idx + 1]);
internal->io.SetParameter("OpenTimeoutSecs", arg[idx + 1]);
++idx;
} else {
error->one(FLERR, "Missing value for 'timeout' option for ADIOS read_dump command");
}
}
++idx;
}
}
/* ----------------------------------------------------------------------
@ -142,26 +128,17 @@ void ReaderADIOS::settings(int narg, char **arg)
Every process must call this Collective operation
------------------------------------------------------------------------- */
void ReaderADIOS::open_file(const char *file)
void ReaderADIOS::open_file(const std::string &file)
{
int rv;
char str[1024];
// close open file, if needed.
if (internal->fh) internal->fh.Close();
// close open file, if needed.
if (internal->fh)
internal->fh.Close();
try {
internal->fh = internal->io.Open(file, adios2::Mode::Read, world);
} catch (std::ios_base::failure &e) {
char str[256];
snprintf(str, sizeof(str), "%s", e.what());
error->one(FLERR, str);
}
if (!internal->fh) {
snprintf(str, sizeof(str), "Cannot open file %s using ADIOS", file);
error->one(FLERR, str);
}
try {
internal->fh = internal->io.Open(file, adios2::Mode::Read, world);
} catch (std::ios_base::failure &e) {
error->one(FLERR, "Error opening file {}: {}", file, e.what());
}
if (!internal->fh) error->one(FLERR, "Cannot open file {} using ADIOS", file);
}
/* ----------------------------------------------------------------------
@ -171,10 +148,8 @@ void ReaderADIOS::open_file(const char *file)
void ReaderADIOS::close_file()
{
// close open file, if needed.
if (internal->fh) {
internal->fh.Close();
}
// close open file, if needed.
if (internal->fh) { internal->fh.Close(); }
}
/* ----------------------------------------------------------------------
@ -185,34 +160,24 @@ void ReaderADIOS::close_file()
int ReaderADIOS::read_time(bigint &ntimestep)
{
char str[1024];
adios2::StepStatus status = internal->fh.BeginStep(adios2::StepMode::Read, internal->timeout);
adios2::StepStatus status =
internal->fh.BeginStep(adios2::StepMode::Read, internal->timeout);
switch (status) {
switch (status) {
case adios2::StepStatus::EndOfStream:
case adios2::StepStatus::NotReady:
case adios2::StepStatus::OtherError:
return 1;
return 1;
default:
break;
}
break;
}
internal->varNtimestep =
internal->io.InquireVariable<uint64_t>("ntimestep");
internal->varNtimestep = internal->io.InquireVariable<uint64_t>("ntimestep");
if (!internal->varNtimestep) {
snprintf(str, sizeof(str),
"Did not find 'ntimestep' variable in ADIOS file %s",
internal->fh.Name().c_str());
error->one(FLERR, str);
}
if (!internal->varNtimestep)
error->one(FLERR, "Did not find 'ntimestep' variable in ADIOS file {}", internal->fh.Name());
ntimestep = static_cast<bigint>(internal->varNtimestep.Max());
// std::cerr << " **** ReaderADIOS::read_time found step " << ntimestep
// << " **** " << std::endl;
return 0;
ntimestep = static_cast<bigint>(internal->varNtimestep.Max());
return 0;
}
/* ----------------------------------------------------------------------
@ -220,7 +185,10 @@ int ReaderADIOS::read_time(bigint &ntimestep)
Called by all processors.
------------------------------------------------------------------------- */
void ReaderADIOS::skip() { internal->fh.EndStep(); }
void ReaderADIOS::skip()
{
internal->fh.EndStep();
}
/* ----------------------------------------------------------------------
read remaining header info:
@ -236,234 +204,205 @@ void ReaderADIOS::skip() { internal->fh.EndStep(); }
only called by proc 0
------------------------------------------------------------------------- */
bigint ReaderADIOS::read_header(double box[3][3], int &boxinfo, int &triclinic,
int fieldinfo, int nfield, int *fieldtype,
char **fieldlabel, int scaleflag, int wrapflag,
int &fieldflag, int &xflag, int &yflag,
int &zflag)
bigint ReaderADIOS::read_header(double box[3][3], int &boxinfo, int &triclinic, int fieldinfo,
int nfield, int *fieldtype, char **fieldlabel, int scaleflag,
int wrapflag, int &fieldflag, int &xflag, int &yflag, int &zflag)
{
char str[1024];
nid = 0;
nid = 0;
// signal that we have no box info at all so far.
// signal that we have no box info at all so far.
internal->varNatoms = internal->io.InquireVariable<uint64_t>("natoms");
if (!internal->varNatoms) {
snprintf(str, sizeof(str),
"Did not find 'natoms' variable in ADIOS file %s",
internal->fh.Name().c_str());
error->one(FLERR, str);
internal->varNatoms = internal->io.InquireVariable<uint64_t>("natoms");
if (!internal->varNatoms)
error->one(FLERR, "Did not find 'natoms' variable in ADIOS file {}", internal->fh.Name());
/* nAtoms */
nAtomsTotal = internal->varNatoms.Max();
uint64_t rem = nAtomsTotal % comm->nprocs;
nAtoms = nAtomsTotal / comm->nprocs;
atomOffset = comm->me * nAtoms;
if (comm->me < (int) rem) {
++nAtoms;
atomOffset += comm->me;
} else {
atomOffset += rem;
}
/* triclinic */
adios2::Attribute<int32_t> attTriclinic = internal->io.InquireAttribute<int32_t>("triclinic");
if (!attTriclinic)
error->one(FLERR, "Did not find 'triclinic' attribute in ADIOS file {}", internal->fh.Name());
triclinic = attTriclinic.Data()[0];
/* read Box */
adios2::Variable<double> varBoxxlo = internal->io.InquireVariable<double>("boxxlo");
adios2::Variable<double> varBoxxhi = internal->io.InquireVariable<double>("boxxhi");
adios2::Variable<double> varBoxylo = internal->io.InquireVariable<double>("boxylo");
adios2::Variable<double> varBoxyhi = internal->io.InquireVariable<double>("boxyhi");
adios2::Variable<double> varBoxzlo = internal->io.InquireVariable<double>("boxzlo");
adios2::Variable<double> varBoxzhi = internal->io.InquireVariable<double>("boxzhi");
box[0][0] = varBoxxlo.Max();
box[0][1] = varBoxxhi.Max();
box[0][2] = 0.0;
box[1][0] = varBoxylo.Max();
box[1][1] = varBoxyhi.Max();
box[1][2] = 0.0;
box[2][0] = varBoxzlo.Max();
box[2][1] = varBoxzhi.Max();
box[2][2] = 0.0;
if (triclinic) {
adios2::Variable<double> varBoxxy = internal->io.InquireVariable<double>("boxxy");
adios2::Variable<double> varBoxxz = internal->io.InquireVariable<double>("boxxz");
adios2::Variable<double> varBoxyz = internal->io.InquireVariable<double>("boxyz");
box[0][2] = varBoxxy.Max();
box[1][2] = varBoxxz.Max();
box[2][2] = varBoxyz.Max();
}
boxinfo = 1;
// if no field info requested, just return
if (!fieldinfo) return nAtoms;
memory->create(fieldindex, nfield, "read_dump:fieldindex");
/* Columns */
adios2::Attribute<std::string> attColumns = internal->io.InquireAttribute<std::string>("columns");
std::vector<std::string> labelVector = attColumns.Data();
int nwords = labelVector.size();
std::map<std::string, int> labels;
for (int i = 0; i < nwords; ++i) { labels.emplace(labelVector[i], i); }
int s_index, u_index, su_index;
xflag = UNSET;
yflag = UNSET;
zflag = UNSET;
// copy fieldtype list for supported fields
for (int i = 0; i < nfield; i++) {
if (fieldlabel[i]) {
fieldindex[i] = find_label(fieldlabel[i], labels);
if (fieldtype[i] == X)
xflag = 2 * scaleflag + wrapflag + 1;
else if (fieldtype[i] == Y)
yflag = 2 * scaleflag + wrapflag + 1;
else if (fieldtype[i] == Z)
zflag = 2 * scaleflag + wrapflag + 1;
}
/* nAtoms */
nAtomsTotal = internal->varNatoms.Max();
uint64_t rem = nAtomsTotal % comm->nprocs;
nAtoms = nAtomsTotal / comm->nprocs;
atomOffset = comm->me * nAtoms;
if (comm->me < rem) {
++nAtoms;
atomOffset += comm->me;
} else {
atomOffset += rem;
}
else if (fieldtype[i] == ID)
fieldindex[i] = find_label("id", labels);
else if (fieldtype[i] == TYPE)
fieldindex[i] = find_label("type", labels);
/* triclinic */
adios2::Attribute<int32_t> attTriclinic =
internal->io.InquireAttribute<int32_t>("triclinic");
if (!attTriclinic) {
snprintf(str, sizeof(str),
"Did not find 'triclinic' attribute in ADIOS file %s",
internal->fh.Name().c_str());
error->one(FLERR, str);
}
triclinic = attTriclinic.Data()[0];
/* read Box */
adios2::Variable<double> varBoxxlo =
internal->io.InquireVariable<double>("boxxlo");
adios2::Variable<double> varBoxxhi =
internal->io.InquireVariable<double>("boxxhi");
adios2::Variable<double> varBoxylo =
internal->io.InquireVariable<double>("boxylo");
adios2::Variable<double> varBoxyhi =
internal->io.InquireVariable<double>("boxyhi");
adios2::Variable<double> varBoxzlo =
internal->io.InquireVariable<double>("boxzlo");
adios2::Variable<double> varBoxzhi =
internal->io.InquireVariable<double>("boxzhi");
box[0][0] = varBoxxlo.Max();
box[0][1] = varBoxxhi.Max();
box[0][2] = 0.0;
box[1][0] = varBoxylo.Max();
box[1][1] = varBoxyhi.Max();
box[1][2] = 0.0;
box[2][0] = varBoxzlo.Max();
box[2][1] = varBoxzhi.Max();
box[2][2] = 0.0;
if (triclinic) {
adios2::Variable<double> varBoxxy =
internal->io.InquireVariable<double>("boxxy");
adios2::Variable<double> varBoxxz =
internal->io.InquireVariable<double>("boxxz");
adios2::Variable<double> varBoxyz =
internal->io.InquireVariable<double>("boxyz");
box[0][2] = varBoxxy.Max();
box[1][2] = varBoxxz.Max();
box[2][2] = varBoxyz.Max();
}
boxinfo = 1;
// if no field info requested, just return
if (!fieldinfo)
return nAtoms;
memory->create(fieldindex, nfield, "read_dump:fieldindex");
/* Columns */
adios2::Attribute<std::string> attColumns =
internal->io.InquireAttribute<std::string>("columns");
std::vector<std::string> labelVector = attColumns.Data();
int nwords = labelVector.size();
std::map<std::string, int> labels;
for (int i = 0; i < nwords; ++i) {
labels.emplace(labelVector[i], i);
}
int s_index, u_index, su_index;
xflag = UNSET;
yflag = UNSET;
zflag = UNSET;
// copy fieldtype list for supported fields
for (int i = 0; i < nfield; i++) {
if (fieldlabel[i]) {
fieldindex[i] = find_label(fieldlabel[i], labels);
if (fieldtype[i] == X)
xflag = 2 * scaleflag + wrapflag + 1;
else if (fieldtype[i] == Y)
yflag = 2 * scaleflag + wrapflag + 1;
else if (fieldtype[i] == Z)
zflag = 2 * scaleflag + wrapflag + 1;
else if (fieldtype[i] == X) {
fieldindex[i] = find_label("x", labels);
xflag = NOSCALE_WRAP;
if (fieldindex[i] < 0) {
fieldindex[i] = nwords;
s_index = find_label("xs", labels);
u_index = find_label("xu", labels);
su_index = find_label("xsu", labels);
if (s_index >= 0 && s_index < fieldindex[i]) {
fieldindex[i] = s_index;
xflag = SCALE_WRAP;
}
if (u_index >= 0 && u_index < fieldindex[i]) {
fieldindex[i] = u_index;
xflag = NOSCALE_NOWRAP;
}
if (su_index >= 0 && su_index < fieldindex[i]) {
fieldindex[i] = su_index;
xflag = SCALE_NOWRAP;
}
}
if (fieldindex[i] == nwords) fieldindex[i] = -1;
else if (fieldtype[i] == ID)
fieldindex[i] = find_label("id", labels);
else if (fieldtype[i] == TYPE)
fieldindex[i] = find_label("type", labels);
} else if (fieldtype[i] == Y) {
fieldindex[i] = find_label("y", labels);
yflag = NOSCALE_WRAP;
if (fieldindex[i] < 0) {
fieldindex[i] = nwords;
s_index = find_label("ys", labels);
u_index = find_label("yu", labels);
su_index = find_label("ysu", labels);
if (s_index >= 0 && s_index < fieldindex[i]) {
fieldindex[i] = s_index;
yflag = SCALE_WRAP;
}
if (u_index >= 0 && u_index < fieldindex[i]) {
fieldindex[i] = u_index;
yflag = NOSCALE_NOWRAP;
}
if (su_index >= 0 && su_index < fieldindex[i]) {
fieldindex[i] = su_index;
yflag = SCALE_NOWRAP;
}
}
if (fieldindex[i] == nwords) fieldindex[i] = -1;
else if (fieldtype[i] == X) {
fieldindex[i] = find_label("x", labels);
xflag = NOSCALE_WRAP;
if (fieldindex[i] < 0) {
fieldindex[i] = nwords;
s_index = find_label("xs", labels);
u_index = find_label("xu", labels);
su_index = find_label("xsu", labels);
if (s_index >= 0 && s_index < fieldindex[i]) {
fieldindex[i] = s_index;
xflag = SCALE_WRAP;
}
if (u_index >= 0 && u_index < fieldindex[i]) {
fieldindex[i] = u_index;
xflag = NOSCALE_NOWRAP;
}
if (su_index >= 0 && su_index < fieldindex[i]) {
fieldindex[i] = su_index;
xflag = SCALE_NOWRAP;
}
}
if (fieldindex[i] == nwords)
fieldindex[i] = -1;
} else if (fieldtype[i] == Z) {
fieldindex[i] = find_label("z", labels);
zflag = NOSCALE_WRAP;
if (fieldindex[i] < 0) {
fieldindex[i] = nwords;
s_index = find_label("zs", labels);
u_index = find_label("zu", labels);
su_index = find_label("zsu", labels);
if (s_index >= 0 && s_index < fieldindex[i]) {
fieldindex[i] = s_index;
zflag = SCALE_WRAP;
}
if (u_index >= 0 && u_index < fieldindex[i]) {
fieldindex[i] = u_index;
zflag = NOSCALE_NOWRAP;
}
if (su_index >= 0 && su_index < fieldindex[i]) {
fieldindex[i] = su_index;
zflag = SCALE_NOWRAP;
}
}
if (fieldindex[i] == nwords) fieldindex[i] = -1;
} else if (fieldtype[i] == Y) {
fieldindex[i] = find_label("y", labels);
yflag = NOSCALE_WRAP;
if (fieldindex[i] < 0) {
fieldindex[i] = nwords;
s_index = find_label("ys", labels);
u_index = find_label("yu", labels);
su_index = find_label("ysu", labels);
if (s_index >= 0 && s_index < fieldindex[i]) {
fieldindex[i] = s_index;
yflag = SCALE_WRAP;
}
if (u_index >= 0 && u_index < fieldindex[i]) {
fieldindex[i] = u_index;
yflag = NOSCALE_NOWRAP;
}
if (su_index >= 0 && su_index < fieldindex[i]) {
fieldindex[i] = su_index;
yflag = SCALE_NOWRAP;
}
}
if (fieldindex[i] == nwords)
fieldindex[i] = -1;
} else if (fieldtype[i] == VX)
fieldindex[i] = find_label("vx", labels);
else if (fieldtype[i] == VY)
fieldindex[i] = find_label("vy", labels);
else if (fieldtype[i] == VZ)
fieldindex[i] = find_label("vz", labels);
} else if (fieldtype[i] == Z) {
fieldindex[i] = find_label("z", labels);
zflag = NOSCALE_WRAP;
if (fieldindex[i] < 0) {
fieldindex[i] = nwords;
s_index = find_label("zs", labels);
u_index = find_label("zu", labels);
su_index = find_label("zsu", labels);
if (s_index >= 0 && s_index < fieldindex[i]) {
fieldindex[i] = s_index;
zflag = SCALE_WRAP;
}
if (u_index >= 0 && u_index < fieldindex[i]) {
fieldindex[i] = u_index;
zflag = NOSCALE_NOWRAP;
}
if (su_index >= 0 && su_index < fieldindex[i]) {
fieldindex[i] = su_index;
zflag = SCALE_NOWRAP;
}
}
if (fieldindex[i] == nwords)
fieldindex[i] = -1;
else if (fieldtype[i] == FX)
fieldindex[i] = find_label("fx", labels);
else if (fieldtype[i] == FY)
fieldindex[i] = find_label("fy", labels);
else if (fieldtype[i] == FZ)
fieldindex[i] = find_label("fz", labels);
} else if (fieldtype[i] == VX)
fieldindex[i] = find_label("vx", labels);
else if (fieldtype[i] == VY)
fieldindex[i] = find_label("vy", labels);
else if (fieldtype[i] == VZ)
fieldindex[i] = find_label("vz", labels);
else if (fieldtype[i] == Q)
fieldindex[i] = find_label("q", labels);
else if (fieldtype[i] == FX)
fieldindex[i] = find_label("fx", labels);
else if (fieldtype[i] == FY)
fieldindex[i] = find_label("fy", labels);
else if (fieldtype[i] == FZ)
fieldindex[i] = find_label("fz", labels);
else if (fieldtype[i] == IX)
fieldindex[i] = find_label("ix", labels);
else if (fieldtype[i] == IY)
fieldindex[i] = find_label("iy", labels);
else if (fieldtype[i] == IZ)
fieldindex[i] = find_label("iz", labels);
}
else if (fieldtype[i] == Q)
fieldindex[i] = find_label("q", labels);
// set fieldflag = -1 if any unfound fields
else if (fieldtype[i] == IX)
fieldindex[i] = find_label("ix", labels);
else if (fieldtype[i] == IY)
fieldindex[i] = find_label("iy", labels);
else if (fieldtype[i] == IZ)
fieldindex[i] = find_label("iz", labels);
}
fieldflag = 0;
for (int i = 0; i < nfield; i++)
if (fieldindex[i] < 0) fieldflag = -1;
// set fieldflag = -1 if any unfound fields
fieldflag = 0;
for (int i = 0; i < nfield; i++)
if (fieldindex[i] < 0)
fieldflag = -1;
return nAtoms;
return nAtoms;
}
/* ----------------------------------------------------------------------
@ -475,41 +414,33 @@ bigint ReaderADIOS::read_header(double box[3][3], int &boxinfo, int &triclinic,
void ReaderADIOS::read_atoms(int n, int nfield, double **fields)
{
char str[1024];
/* Read Atoms */
/* This is the firsts (and last) read of array data, so we
* call EndStep() here instead of PerformGets()
*/
/* Read Atoms */
/* This is the firsts (and last) read of array data, so we
* call EndStep() here instead of PerformGets()
*/
adios2::Variable<double> varAtoms = internal->io.InquireVariable<double>("atoms");
adios2::Variable<double> varAtoms =
internal->io.InquireVariable<double>("atoms");
if ((uint64_t) n != nAtoms)
error->one(FLERR,
"ReaderADIOS::read_atoms() expects 'n={}' equal to the number of "
"atoms (={}) for process {} in ADIOS file {}.",
n, nAtoms, comm->me, internal->fh.Name());
if (n != nAtoms) {
snprintf(
str, sizeof(str),
"ReaderADIOS::read_atoms() expects 'n=%d' equal to the number of "
"atoms (=%" PRIu64 ") for process %d in ADIOS file %s.",
n, nAtoms, comm->me, internal->fh.Name().c_str());
error->one(FLERR, str);
}
size_t ncols = varAtoms.Count()[1];
varAtoms.SetSelection({{atomOffset, 0}, {nAtoms, ncols}});
size_t ncols = varAtoms.Count()[1];
varAtoms.SetSelection({{atomOffset, 0}, {nAtoms, ncols}});
std::vector<double> table;
std::vector<double> table;
internal->fh.Get<double>(varAtoms, table);
// EndStep or PerformGets required to make the read happen
internal->fh.EndStep();
internal->fh.Get<double>(varAtoms, table);
// EndStep or PerformGets required to make the read happen
internal->fh.EndStep();
size_t idx;
for (int i = 0; i < nAtoms; i++) {
idx = i * ncols;
for (int m = 0; m < nfield; m++) {
fields[i][m] = table[idx + fieldindex[m]];
}
}
size_t idx;
for (uint64_t i = 0; i < nAtoms; i++) {
idx = i * ncols;
for (int m = 0; m < nfield; m++) { fields[i][m] = table[idx + fieldindex[m]]; }
}
}
/* ----------------------------------------------------------------------
@ -517,12 +448,9 @@ void ReaderADIOS::read_atoms(int n, int nfield, double **fields)
return index of match or -1 if no match
------------------------------------------------------------------------- */
int ReaderADIOS::find_label(const std::string &label,
const std::map<std::string, int> &labels)
int ReaderADIOS::find_label(const std::string &label, const std::map<std::string, int> &labels)
{
std::map<std::string, int>::const_iterator it = labels.find(label);
if (it != labels.end()) {
return it->second;
}
return -1;
auto it = labels.find(label);
if (it != labels.end()) { return it->second; }
return -1;
}

18
src/ADIOS/reader_adios.h
View File

@ -34,18 +34,18 @@ class ReadADIOSInternal;
class ReaderADIOS : public Reader {
public:
ReaderADIOS(class LAMMPS *);
virtual ~ReaderADIOS();
~ReaderADIOS() override;
virtual void settings(int, char **);
void settings(int, char **) override;
virtual int read_time(bigint &);
virtual void skip();
virtual bigint read_header(double[3][3], int &, int &, int, int, int *, char **, int, int, int &,
int &, int &, int &);
virtual void read_atoms(int, int, double **);
int read_time(bigint &) override;
void skip() override;
bigint read_header(double[3][3], int &, int &, int, int, int *, char **, int, int, int &, int &,
int &, int &) override;
void read_atoms(int, int, double **) override;
virtual void open_file(const char *);
virtual void close_file();
void open_file(const std::string &) override;
void close_file() override;
private:
int *fieldindex; // mapping of input fields to dump

6
src/ASPHERE/compute_erotate_asphere.cpp
View File

@ -43,9 +43,9 @@ void ComputeERotateAsphere::init()
{
// error check
avec_ellipsoid = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
avec_line = (AtomVecLine *) atom->style_match("line");
avec_tri = (AtomVecTri *) atom->style_match("tri");
avec_ellipsoid = dynamic_cast<AtomVecEllipsoid *>( atom->style_match("ellipsoid"));
avec_line = dynamic_cast<AtomVecLine *>( atom->style_match("line"));
avec_tri = dynamic_cast<AtomVecTri *>( atom->style_match("tri"));
if (!avec_ellipsoid && !avec_line && !avec_tri)
error->all(FLERR,"Compute erotate/asphere requires "
"atom style ellipsoid or line or tri");

22
src/ASPHERE/compute_erotate_asphere.h
View File

@ -27,8 +27,8 @@ namespace LAMMPS_NS {
class ComputeERotateAsphere : public Compute {
public:
ComputeERotateAsphere(class LAMMPS *, int, char **);
void init();
double compute_scalar();
void init() override;
double compute_scalar() override;
private:
double pfactor;
@ -41,21 +41,3 @@ class ComputeERotateAsphere : public Compute {
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Compute erotate/asphere requires atom style ellipsoid or line or tri
Self-explanatory.
E: Compute erotate/asphere requires extended particles
This compute cannot be used with point particles.
*/

2
src/ASPHERE/compute_temp_asphere.cpp
View File

@ -94,7 +94,7 @@ void ComputeTempAsphere::init()
{
// error check
avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
avec = dynamic_cast<AtomVecEllipsoid *>( atom->style_match("ellipsoid"));
if (!avec)
error->all(FLERR,"Compute temp/asphere requires atom style ellipsoid");

58
src/ASPHERE/compute_temp_asphere.h
View File

@ -27,16 +27,16 @@ namespace LAMMPS_NS {
class ComputeTempAsphere : public Compute {
public:
ComputeTempAsphere(class LAMMPS *, int, char **);
~ComputeTempAsphere();
void init();
void setup();
double compute_scalar();
void compute_vector();
~ComputeTempAsphere() override;
void init() override;
void setup() override;
double compute_scalar() override;
void compute_vector() override;
void remove_bias(int, double *);
void restore_bias(int, double *);
void remove_bias_thr(int, double *, double *);
void restore_bias_thr(int, double *, double *);
void remove_bias(int, double *) override;
void restore_bias(int, double *) override;
void remove_bias_thr(int, double *, double *) override;
void restore_bias_thr(int, double *, double *) override;
private:
int mode;
@ -52,43 +52,3 @@ class ComputeTempAsphere : public Compute {
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Compute temp/asphere requires atom style ellipsoid
Self-explanatory.
E: Compute temp/asphere requires extended particles
This compute cannot be used with point particles.
E: Could not find compute ID for temperature bias
Self-explanatory.
E: Bias compute does not calculate temperature
The specified compute must compute temperature.
E: Bias compute does not calculate a velocity bias
The specified compute must compute a bias for temperature.
E: Bias compute group does not match compute group
The specified compute must operate on the same group as the parent
compute.
E: Temperature compute degrees of freedom < 0
This should not happen if you are calculating the temperature
on a valid set of atoms.
*/

2
src/ASPHERE/fix_nh_asphere.cpp
View File

@ -36,7 +36,7 @@ FixNHAsphere::FixNHAsphere(LAMMPS *lmp, int narg, char **arg) :
void FixNHAsphere::init()
{
avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
avec = dynamic_cast<AtomVecEllipsoid *>( atom->style_match("ellipsoid"));
if (!avec)
error->all(FLERR,
"Compute nvt/nph/npt asphere requires atom style ellipsoid");

22
src/ASPHERE/fix_nh_asphere.h
View File

@ -21,31 +21,17 @@ namespace LAMMPS_NS {
class FixNHAsphere : public FixNH {
public:
FixNHAsphere(class LAMMPS *, int, char **);
virtual ~FixNHAsphere() {}
void init();
void init() override;
protected:
double dtq;
class AtomVecEllipsoid *avec;
void nve_v();
void nve_x();
void nh_v_temp();
void nve_v() override;
void nve_x() override;
void nh_v_temp() override;
};
} // namespace LAMMPS_NS
#endif
/* ERROR/WARNING messages:
E: Compute nvt/nph/npt asphere requires atom style ellipsoid
Self-explanatory.
E: Fix nvt/nph/npt asphere requires extended particles
The shape setting for a particle in the fix group has shape = 0.0,
which means it is a point particle.
*/

13
src/ASPHERE/fix_nph_asphere.h
View File

@ -27,22 +27,9 @@ namespace LAMMPS_NS {
class FixNPHAsphere : public FixNHAsphere {
public:
FixNPHAsphere(class LAMMPS *, int, char **);
~FixNPHAsphere() {}
};
} // namespace LAMMPS_NS
#endif
#endif
/* ERROR/WARNING messages:
E: Temperature control can not be used with fix nph/asphere
Self-explanatory.
E: Pressure control must be used with fix nph/asphere
Self-explanatory.
*/

13
src/ASPHERE/fix_npt_asphere.h
View File

@ -27,22 +27,9 @@ namespace LAMMPS_NS {
class FixNPTAsphere : public FixNHAsphere {
public:
FixNPTAsphere(class LAMMPS *, int, char **);
~FixNPTAsphere() {}
};
} // namespace LAMMPS_NS
#endif
#endif
/* ERROR/WARNING messages:
E: Temperature control must be used with fix npt/asphere
Self-explanatory.
E: Pressure control must be used with fix npt/asphere
Self-explanatory.
*/

2
src/ASPHERE/fix_nve_asphere.cpp
View File

@ -37,7 +37,7 @@ FixNVEAsphere::FixNVEAsphere(LAMMPS *lmp, int narg, char **arg) :
void FixNVEAsphere::init()
{
avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
avec = dynamic_cast<AtomVecEllipsoid *>( atom->style_match("ellipsoid"));
if (!avec)
error->all(FLERR,"Compute nve/asphere requires atom style ellipsoid");

18
src/ASPHERE/fix_nve_asphere.h
View File

@ -27,9 +27,9 @@ namespace LAMMPS_NS {
class FixNVEAsphere : public FixNVE {
public:
FixNVEAsphere(class LAMMPS *, int, char **);
void init();
void initial_integrate(int);
void final_integrate();
void init() override;
void initial_integrate(int) override;
void final_integrate() override;
private:
double dtq;
@ -39,15 +39,3 @@ class FixNVEAsphere : public FixNVE {
} // namespace LAMMPS_NS
#endif
#endif
/* ERROR/WARNING messages:
E: Compute nve/asphere requires atom style ellipsoid
Self-explanatory.
E: Fix nve/asphere requires extended particles
This fix can only be used for particles with a shape setting.
*/

2
src/ASPHERE/fix_nve_asphere_noforce.cpp
View File

@ -37,7 +37,7 @@ void FixNVEAsphereNoforce::init()
{
// error check
avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
avec = dynamic_cast<AtomVecEllipsoid *>( atom->style_match("ellipsoid"));
if (!atom->ellipsoid_flag)
error->all(FLERR,"Fix nve/asphere/noforce requires atom style ellipsoid");

22
src/ASPHERE/fix_nve_asphere_noforce.h
View File

@ -27,8 +27,8 @@ namespace LAMMPS_NS {
class FixNVEAsphereNoforce : public FixNVENoforce {
public:
FixNVEAsphereNoforce(class LAMMPS *, int, char **);
void initial_integrate(int);
void init();
void initial_integrate(int) override;
void init() override;
private:
double dtq;
@ -39,21 +39,3 @@ class FixNVEAsphereNoforce : public FixNVENoforce {
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Fix nve/asphere/noforce requires atom style ellipsoid
Self-explanatory.
E: Fix nve/asphere/noforce requires extended particles
One of the particles is not an ellipsoid.
*/

2
src/ASPHERE/fix_nve_line.cpp
View File

@ -56,7 +56,7 @@ void FixNVELine::init()
{
// error checks
avec = (AtomVecLine *) atom->style_match("line");
avec = dynamic_cast<AtomVecLine *>( atom->style_match("line"));
if (!avec) error->all(FLERR,"Fix nve/line requires atom style line");
if (domain->dimension != 2)

31
src/ASPHERE/fix_nve_line.h
View File

@ -27,11 +27,10 @@ namespace LAMMPS_NS {
class FixNVELine : public FixNVE {
public:
FixNVELine(class LAMMPS *, int, char **);
~FixNVELine() {}
int setmask();
void init();
void initial_integrate(int);
void final_integrate();
int setmask() override;
void init() override;
void initial_integrate(int) override;
void final_integrate() override;
private:
double MINUSPI, TWOPI;
@ -42,25 +41,3 @@ class FixNVELine : public FixNVE {
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Fix nve/line requires atom style line
Self-explanatory.
E: Fix nve/line can only be used for 2d simulations
Self-explanatory.
E: Fix nve/line requires line particles
Self-explanatory.
*/

2
src/ASPHERE/fix_nve_tri.cpp
View File

@ -50,7 +50,7 @@ void FixNVETri::init()
{
// error checks
avec = (AtomVecTri *) atom->style_match("tri");
avec = dynamic_cast<AtomVecTri *>( atom->style_match("tri"));
if (!avec) error->all(FLERR,"Fix nve/tri requires atom style tri");
if (domain->dimension != 3)

31
src/ASPHERE/fix_nve_tri.h
View File

@ -27,11 +27,10 @@ namespace LAMMPS_NS {
class FixNVETri : public FixNVE {
public:
FixNVETri(class LAMMPS *, int, char **);
~FixNVETri() {}
int setmask();
void init();
void initial_integrate(int);
void final_integrate();
int setmask() override;
void init() override;
void initial_integrate(int) override;
void final_integrate() override;
private:
double dtq;
@ -42,25 +41,3 @@ class FixNVETri : public FixNVE {
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Fix nve/tri requires atom style tri
Self-explanatory.
E: Fix nve/tri can only be used for 3d simulations
Self-explanatory.
E: Fix nve/tri requires tri particles
Self-explanatory.
*/

13
src/ASPHERE/fix_nvt_asphere.h
View File

@ -27,22 +27,9 @@ namespace LAMMPS_NS {
class FixNVTAsphere : public FixNHAsphere {
public:
FixNVTAsphere(class LAMMPS *, int, char **);
~FixNVTAsphere() {}
};
} // namespace LAMMPS_NS
#endif
#endif
/* ERROR/WARNING messages:
E: Temperature control must be used with fix nvt/asphere
Self-explanatory.
E: Pressure control can not be used with fix nvt/asphere
Self-explanatory.
*/

18
src/ASPHERE/pair_gayberne.cpp
View File

@ -18,18 +18,18 @@
#include "pair_gayberne.h"
#include <cmath>
#include "math_extra.h"
#include "atom.h"
#include "atom_vec_ellipsoid.h"
#include "comm.h"
#include "force.h"
#include "neighbor.h"
#include "neigh_list.h"
#include "citeme.h"
#include "memory.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "math_extra.h"
#include "memory.h"
#include "neigh_list.h"
#include "neighbor.h"
#include <cmath>
using namespace LAMMPS_NS;
@ -345,10 +345,10 @@ void PairGayBerne::coeff(int narg, char **arg)
void PairGayBerne::init_style()
{
avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
avec = dynamic_cast<AtomVecEllipsoid *>( atom->style_match("ellipsoid"));
if (!avec) error->all(FLERR,"Pair gayberne requires atom style ellipsoid");
neighbor->request(this,instance_me);
neighbor->add_request(this,NeighConst::REQ_DEFAULT);
// per-type shape precalculations
// require that atom shapes are identical within each type

55
src/ASPHERE/pair_gayberne.h
View File

@ -27,18 +27,18 @@ namespace LAMMPS_NS {
class PairGayBerne : public Pair {
public:
PairGayBerne(LAMMPS *lmp);
virtual ~PairGayBerne();
virtual void compute(int, int);
virtual void settings(int, char **);
void coeff(int, char **);
virtual void init_style();
double init_one(int, int);
void write_restart(FILE *);
void read_restart(FILE *);
void write_restart_settings(FILE *);
void read_restart_settings(FILE *);
void write_data(FILE *);
void write_data_all(FILE *);
~PairGayBerne() override;
void compute(int, int) override;
void settings(int, char **) override;
void coeff(int, char **) override;
void init_style() override;
double init_one(int, int) override;
void write_restart(FILE *) override;
void read_restart(FILE *) override;
void write_restart_settings(FILE *) override;
void read_restart_settings(FILE *) override;
void write_data(FILE *) override;
void write_data_all(FILE *) override;
protected:
enum { SPHERE_SPHERE, SPHERE_ELLIPSE, ELLIPSE_SPHERE, ELLIPSE_ELLIPSE };
@ -72,34 +72,3 @@ class PairGayBerne : public Pair {
} // namespace LAMMPS_NS
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
E: Pair gayberne requires atom style ellipsoid
Self-explanatory.
E: Pair gayberne requires atoms with same type have same shape
Self-explanatory.
E: Pair gayberne epsilon a,b,c coeffs are not all set
Each atom type involved in pair_style gayberne must
have these 3 coefficients set at least once.
E: Bad matrix inversion in mldivide3
This error should not occur unless the matrix is badly formed.
*/

4
src/ASPHERE/pair_line_lj.cpp
View File

@ -402,10 +402,10 @@ void PairLineLJ::coeff(int narg, char **arg)
void PairLineLJ::init_style()
{
avec = (AtomVecLine *) atom->style_match("line");
avec = dynamic_cast<AtomVecLine *>( atom->style_match("line"));
if (!avec) error->all(FLERR,"Pair line/lj requires atom style line");
neighbor->request(this,instance_me);
neighbor->add_request(this,NeighConst::REQ_DEFAULT);
}
/* ----------------------------------------------------------------------

35
src/ASPHERE/pair_line_lj.h
View File

@ -27,12 +27,12 @@ namespace LAMMPS_NS {
class PairLineLJ : public Pair {
public:
PairLineLJ(class LAMMPS *);
virtual ~PairLineLJ();
virtual void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
virtual void init_style();
double init_one(int, int);
~PairLineLJ() override;
void compute(int, int) override;
void settings(int, char **) override;
void coeff(int, char **) override;
void init_style() override;
double init_one(int, int) override;
protected:
double cut_global;
@ -62,26 +62,3 @@ class PairLineLJ : public Pair {
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
E: Pair line/lj requires atom style line
Self-explanatory.
E: All pair coeffs are not set
All pair coefficients must be set in the data file or by the
pair_coeff command before running a simulation.
*/

4
src/ASPHERE/pair_resquared.cpp
View File

@ -317,10 +317,10 @@ void PairRESquared::coeff(int narg, char **arg)
void PairRESquared::init_style()
{
avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
avec = dynamic_cast<AtomVecEllipsoid *>(atom->style_match("ellipsoid"));
if (!avec) error->all(FLERR, "Pair resquared requires atom style ellipsoid");
neighbor->request(this, instance_me);
neighbor->add_request(this, NeighConst::REQ_DEFAULT);
// per-type shape precalculations
// require that atom shapes are identical within each type

54
src/ASPHERE/pair_resquared.h
View File

@ -27,16 +27,16 @@ namespace LAMMPS_NS {
class PairRESquared : public Pair {
public:
PairRESquared(LAMMPS *lmp);
virtual ~PairRESquared();
virtual void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
virtual void init_style();
double init_one(int, int);
void write_restart(FILE *);
void read_restart(FILE *);
void write_restart_settings(FILE *);
void read_restart_settings(FILE *);
~PairRESquared() override;
void compute(int, int) override;
void settings(int, char **) override;
void coeff(int, char **) override;
void init_style() override;
double init_one(int, int) override;
void write_restart(FILE *) override;
void read_restart(FILE *) override;
void write_restart_settings(FILE *) override;
void read_restart_settings(FILE *) override;
protected:
enum { SPHERE_SPHERE, SPHERE_ELLIPSE, ELLIPSE_SPHERE, ELLIPSE_ELLIPSE };
@ -92,37 +92,3 @@ class PairRESquared : public Pair {
} // namespace LAMMPS_NS
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
E: Pair resquared requires atom style ellipsoid
Self-explanatory.
E: Pair resquared requires atoms with same type have same shape
Self-explanatory.
E: Pair resquared epsilon a,b,c coeffs are not all set
Self-explanatory.
E: Pair resquared epsilon and sigma coeffs are not all set
Self-explanatory.
E: Bad matrix inversion in mldivide3
This error should not occur unless the matrix is badly formed.
*/

4
src/ASPHERE/pair_tri_lj.cpp
View File

@ -467,10 +467,10 @@ void PairTriLJ::coeff(int narg, char **arg)
void PairTriLJ::init_style()
{
avec = (AtomVecTri *) atom->style_match("tri");
avec = dynamic_cast<AtomVecTri *>( atom->style_match("tri"));
if (!avec) error->all(FLERR,"Pair tri/lj requires atom style tri");
neighbor->request(this,instance_me);
neighbor->add_request(this,NeighConst::REQ_DEFAULT);
}
/* ----------------------------------------------------------------------

30
src/ASPHERE/pair_tri_lj.h
View File

@ -27,12 +27,12 @@ namespace LAMMPS_NS {
class PairTriLJ : public Pair {
public:
PairTriLJ(class LAMMPS *);
virtual ~PairTriLJ();
virtual void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
virtual void init_style();
double init_one(int, int);
~PairTriLJ() override;
void compute(int, int) override;
void settings(int, char **) override;
void coeff(int, char **) override;
void init_style() override;
double init_one(int, int) override;
protected:
double cut_global;
@ -60,21 +60,3 @@ class PairTriLJ : public Pair {
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
E: Pair tri/lj requires atom style tri
Self-explanatory.
*/

11
src/ATC/fix_atc.cpp
View File

@ -18,7 +18,6 @@
#include "comm.h"
#include "error.h"
#include "group.h"
#include "neigh_request.h"
#include "neighbor.h"
#include "ATC_Method.h"
@ -558,11 +557,7 @@ int FixATC::modify_param(int narg, char** arg)
void FixATC::init()
{
// Guarantee construction of full neighborlist
int irequest = neighbor->request(this,instance_me);
neighbor->requests[irequest]->pair = 0;
neighbor->requests[irequest]->fix = 1;
neighbor->requests[irequest]->half = 0;
neighbor->requests[irequest]->full = 1;
neighbor->add_request(this, NeighConst::REQ_FULL);
// create computes, if necessary
atc_->init_computes();
@ -575,7 +570,7 @@ void FixATC::min_setup(int vflag)
void FixATC::setup(int /* vflag */)
{
comm->forward_comm_fix(this);
comm->forward_comm(this);
try {
atc_->initialize();
@ -814,7 +809,7 @@ void FixATC::pre_neighbor()
{
try {
atc_->pre_neighbor();
comm->forward_comm_fix(this);
comm->forward_comm(this);
}
catch (ATC::ATC_Error& atcError) {
ATC::LammpsInterface::instance()->print_msg(atcError.error_description());

78
src/ATC/fix_atc.h
View File

@ -37,91 +37,91 @@ class FixATC : public Fix {
public:
/** constructor & destructor */
FixATC(class LAMMPS *, int, char **);
~FixATC();
~FixATC() override;
/** initialization functions */
void init();
void init_list(int id, NeighList *ptr);
void setup(int vflag);
void min_setup(int vflag);
void init() override;
void init_list(int id, NeighList *ptr) override;
void setup(int vflag) override;
void min_setup(int vflag) override;
/** setmask: tell LAMMPS which fix methods to call */
int setmask();
int setmask() override;
/** initial_integrate */
void initial_integrate(int vflag);
void initial_integrate(int vflag) override;
/** after first integrate phase */
void post_integrate();
void post_integrate() override;
/** final_integrate */
void final_integrate();
void final_integrate() override;
/** end of step for run or minimize */
void end_of_step();
void end_of_step() override;
/** pre_exchange is used to modify fix-specific data
and is called before domain->pbc() and comm->exchange(). */
void setup_pre_exchange();
void pre_exchange();
void min_pre_exchange();
void setup_pre_exchange() override;
void pre_exchange() override;
void min_pre_exchange() override;
double memory_usage();
void grow_arrays(int);
void copy_arrays(int, int, int);
double memory_usage() override;
void grow_arrays(int) override;
void copy_arrays(int, int, int) override;
/** pack_exchange called from atom_vec->pack_exchange()
and packs fix-specific data for a given real (local)
atom being sent to another processor. */
int pack_exchange(int, double *);
int pack_exchange(int, double *) override;
/** unpack_exchange called from atom_vec->unpack_exchange()
and unpacks fix-specific data for a given real (local)
atom received from another processor. */
int unpack_exchange(int, double *);
int unpack_exchange(int, double *) override;
/** pack_comm called from comm->forward_comm_fix and
/** pack_comm called from comm->forward_comm and
packs fix-specific data for a given ghost atom
from exchange with another proc */
int pack_forward_comm(int, int *, double *, int, int *);
int pack_forward_comm(int, int *, double *, int, int *) override;
/** unpack_comm called from comm->forward_comm_fix and
/** unpack_comm called from comm->forward_comm and
unpacks fix-specific data for a given ghost atom
from exchange with another proc */
void unpack_forward_comm(int, int, double *);
void unpack_forward_comm(int, int, double *) override;
/** pre_neighbor is used to modify fix-specific data
and is called before neighbor list is built in
neighbor->build(). */
void pre_neighbor();
void setup_pre_neighbor();
void pre_neighbor() override;
void setup_pre_neighbor() override;
/** pre/post_force is used to modify fix-specific data
and is before/after the various force computations. */
void pre_force(int vflag);
void post_force(int vflag);
void pre_force(int vflag) override;
void post_force(int vflag) override;
/** post_run is called after a run completes */
void post_run();
void post_run() override;
/** min_pre_force is called before forces are calculated in minimize */
void min_pre_force(int vflag);
void min_pre_force(int vflag) override;
/** min_post_force is called after forces are calculated in minimize */
void min_post_force(int vflag);
void min_post_force(int vflag) override;
/** modify atc parameters (parser) */
int modify_param(int narg, char **arg);
int modify_param(int narg, char **arg) override;
/** calls ATC_Method to handle restarting/checkpointing */
/** these four methods are for writing per-atom quantities */
int pack_restart(int, double *);
void unpack_restart(int, int);
int size_restart(int);
int maxsize_restart();
int pack_restart(int, double *) override;
void unpack_restart(int, int) override;
int size_restart(int) override;
int maxsize_restart() override;
/** these two methods are for writing all other quantities */
void write_restart(FILE *);
void restart(char *);
void write_restart(FILE *) override;
void restart(char *) override;
/** accessor function for ATC_Method class pointer */
const ATC::ATC_Method *atc() { return atc_; }
@ -130,9 +130,9 @@ class FixATC : public Fix {
LAMMPS *lammps_;
/** functions for "thermo" output */
virtual double compute_scalar();
virtual double compute_vector(int n);
virtual double compute_array(int irow, int icol);
double compute_scalar() override;
double compute_vector(int n) override;
double compute_array(int irow, int icol) override;
double dtv, dtf;
ATC::ATC_Method *atc_;
};

78
src/AWPMD/atom_vec_wavepacket.cpp
View File

@ -1,4 +1,3 @@
// clang-format off
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
@ -20,8 +19,6 @@
#include "atom.h"
#include <cstring>
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
@ -32,33 +29,29 @@ AtomVecWavepacket::AtomVecWavepacket(LAMMPS *lmp) : AtomVec(lmp)
molecular = Atom::ATOMIC;
forceclearflag = 1;
atom->wavepacket_flag = 1;
atom->electron_flag = 1; // compatible with eff
atom->q_flag = atom->spin_flag = atom->eradius_flag =
atom->ervel_flag = atom->erforce_flag = 1;
atom->cs_flag = atom->csforce_flag =
atom->vforce_flag = atom->ervelforce_flag = atom->etag_flag = 1;
atom->wavepacket_flag = atom->q_flag = atom->spin_flag = atom->eradius_flag = 1;
atom->ervel_flag = atom->erforce_flag = atom->cs_flag = atom->csforce_flag = 1;
atom->vforce_flag = atom->ervelforce_flag = atom->etag_flag = 1;
// strings with peratom variables to include in each AtomVec method
// strings cannot contain fields in corresponding AtomVec default strings
// order of fields in a string does not matter
// except: fields_data_atom & fields_data_vel must match data file
fields_grow = (char *)
"q spin eradius ervel erforce cs csforce "
"vforce ervelforce etag";
fields_copy = (char *) "q spin eradius ervel cs etag";
fields_comm = (char *) "eradius";
fields_comm_vel = (char *) "eradius ervel cs";
fields_reverse = (char *) "erforce ervelforce vforce csforce";
fields_border = (char *) "q spin eradius etag";
fields_border_vel = (char *) "q spin eradius etag ervel cs";
fields_exchange = (char *) "q spin eradius ervel etag cs";
fields_restart = (char *) "q spin eradius ervel etag cs";
fields_create = (char *) "q spin eradius ervel etag cs";
fields_data_atom = (char *) "id type q spin eradius etag cs x";
fields_data_vel = (char *) "id v ervel";
fields_grow = {"q", "spin", "eradius", "ervel", "erforce",
"cs", "csforce", "vforce", "ervelforce", "etag"};
fields_copy = {"q", "spin", "eradius", "ervel", "cs", "etag"};
fields_comm = {"eradius"};
fields_comm_vel = {"eradius", "ervel", "cs"};
fields_reverse = {"erforce", "ervelforce", "vforce", "csforce"};
fields_border = {"q", "spin", "eradius", "etag"};
fields_border_vel = {"q", "spin", "eradius", "etag", "ervel", "cs"};
fields_exchange = {"q", "spin", "eradius", "ervel", "etag", "cs"};
fields_restart = {"q", "spin", "eradius", "ervel", "etag", "cs"};
fields_create = {"q", "spin", "eradius", "ervel", "etag", "cs"};
fields_data_atom = {"id", "type", "q", "spin", "eradius", "etag", "cs", "x"};
fields_data_vel = {"id", "v", "ervel"};
setup_fields();
}
@ -84,7 +77,7 @@ void AtomVecWavepacket::grow_pointers()
void AtomVecWavepacket::force_clear(int n, size_t nbytes)
{
memset(&erforce[n],0,nbytes);
memset(&erforce[n], 0, nbytes);
}
/* ----------------------------------------------------------------------
@ -112,12 +105,12 @@ void AtomVecWavepacket::data_atom_post(int ilocal)
return -1 if name is unknown to this atom style
------------------------------------------------------------------------- */
int AtomVecWavepacket::property_atom(char *name)
int AtomVecWavepacket::property_atom(const std::string &name)
{
if (strcmp(name,"spin") == 0) return 0;
if (strcmp(name,"eradius") == 0) return 1;
if (strcmp(name,"ervel") == 0) return 2;
if (strcmp(name,"erforce") == 0) return 3;
if (name == "spin") return 0;
if (name == "eradius") return 1;
if (name == "ervel") return 2;
if (name == "erforce") return 3;
return -1;
}
@ -126,34 +119,41 @@ int AtomVecWavepacket::property_atom(char *name)
index maps to data specific to this atom style
------------------------------------------------------------------------- */
void AtomVecWavepacket::pack_property_atom(int index, double *buf,
int nvalues, int groupbit)
void AtomVecWavepacket::pack_property_atom(int index, double *buf, int nvalues, int groupbit)
{
int nlocal = atom->nlocal;
int n = 0;
if (index == 0) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = spin[i];
else buf[n] = 0.0;
if (mask[i] & groupbit)
buf[n] = spin[i];
else
buf[n] = 0.0;
n += nvalues;
}
} else if (index == 1) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = eradius[i];
else buf[n] = 0.0;
if (mask[i] & groupbit)
buf[n] = eradius[i];
else
buf[n] = 0.0;
n += nvalues;
}
} else if (index == 2) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = ervel[i];
else buf[n] = 0.0;
if (mask[i] & groupbit)
buf[n] = ervel[i];
else
buf[n] = 0.0;
n += nvalues;
}
} else if (index == 3) {
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) buf[n] = erforce[i];
else buf[n] = 0.0;
if (mask[i] & groupbit)
buf[n] = erforce[i];
else
buf[n] = 0.0;
n += nvalues;
}
}

12
src/AWPMD/atom_vec_wavepacket.h
View File

@ -28,12 +28,12 @@ class AtomVecWavepacket : public AtomVec {
public:
AtomVecWavepacket(class LAMMPS *);
void grow_pointers();
void force_clear(int, size_t);
void create_atom_post(int);
void data_atom_post(int);
int property_atom(char *);
void pack_property_atom(int, double *, int, int);
void grow_pointers() override;
void force_clear(int, size_t) override;
void create_atom_post(int) override;
void data_atom_post(int) override;
int property_atom(const std::string &) override;
void pack_property_atom(int, double *, int, int) override;
private:
int *spin;

4
src/AWPMD/fix_nve_awpmd.cpp
View File

@ -61,9 +61,9 @@ void FixNVEAwpmd::init()
dtf = 0.5 * update->dt * force->ftm2v;
if (utils::strmatch(update->integrate_style,"^respa"))
step_respa = ((Respa *) update->integrate)->step;
step_respa = (dynamic_cast<Respa *>( update->integrate))->step;
awpmd_pair=(PairAWPMDCut *)force->pair;
awpmd_pair=dynamic_cast<PairAWPMDCut *>(force->pair);
awpmd_pair->wpmd->norm_needed=1;
}

14
src/AWPMD/fix_nve_awpmd.h
View File

@ -31,13 +31,13 @@ namespace LAMMPS_NS {
class FixNVEAwpmd : public Fix {
public:
FixNVEAwpmd(class LAMMPS *, int, char **);
int setmask();
virtual void init();
virtual void initial_integrate(int);
virtual void final_integrate();
void initial_integrate_respa(int, int, int);
void final_integrate_respa(int, int);
void reset_dt();
int setmask() override;
void init() override;
void initial_integrate(int) override;
void final_integrate() override;
void initial_integrate_respa(int, int, int) override;
void final_integrate_respa(int, int) override;
void reset_dt() override;
protected:
double dtv, dtf;

41
src/AWPMD/pair_awpmd_cut.cpp
View File

@ -26,7 +26,6 @@
#include "memory.h"
#include "min.h"
#include "neigh_list.h"
#include "neigh_request.h"
#include "neighbor.h"
#include "update.h"
@ -96,10 +95,7 @@ struct cmp_x{
else if (d>tol)
return false;
d=xx[left.second][2]-xx[right.second][2];
if (d<-tol)
return true;
else
return false;
return d < -tol;
}
else
return left.first<right.first;
@ -513,45 +509,18 @@ void PairAWPMDCut::init_style()
if (!atom->q_flag || !atom->spin_flag ||
!atom->eradius_flag || !atom->erforce_flag ) // TO DO: adjust this to match approximation used
error->all(FLERR,"Pair awpmd/cut requires atom attributes "
"q, spin, eradius, erforce");
/*
if (vflag_atom) { // can't compute virial per atom
//warning->
error->all(FLERR,"Pair style awpmd can't compute per atom virials");
}*/
error->all(FLERR,"Pair awpmd/cut requires atom attributes q, spin, eradius, erforce");
// add hook to minimizer for eradius and erforce
if (update->whichflag == 2)
int ignore = update->minimize->request(this,1,0.01);
// make sure to use the appropriate timestep when using real units
/*if (update->whichflag == 1) {
if (force->qqr2e == 332.06371 && update->dt == 1.0)
error->all(FLERR,"You must lower the default real units timestep for pEFF ");
}*/
// need a half neigh list and optionally a granular history neigh list
//int irequest = neighbor->request(this,instance_me);
//if (atom->tag_enable == 0)
// error->all(FLERR,"Pair style reax requires atom IDs");
//if (force->newton_pair == 0)
//error->all(FLERR,"Pair style awpmd requires newton pair on");
//if (strcmp(update->unit_style,"real") != 0 && comm->me == 0)
//error->warning(FLERR,"Not using real units with pair reax");
int irequest = neighbor->request(this,instance_me);
neighbor->requests[irequest]->newton = 2;
neighbor->add_request(this, NeighConst::REQ_NEWTON_OFF);
if (force->e_mass==0. || force->hhmrr2e==0. || force->mvh2r==0.)
error->all(FLERR,"Pair style awpmd requires e_mass and conversions hhmrr2e, mvh2r to be properly set for unit system");
error->all(FLERR,"Pair style awpmd requires e_mass and conversions "
"hhmrr2e, mvh2r to be properly set for unit system");
wpmd->me=force->e_mass;
wpmd->h2_me=force->hhmrr2e/force->e_mass;

28
src/AWPMD/pair_awpmd_cut.h
View File

@ -34,22 +34,22 @@ class PairAWPMDCut : public Pair {
public:
PairAWPMDCut(class LAMMPS *);
virtual ~PairAWPMDCut();
virtual void compute(int, int);
virtual void settings(int, char **);
void coeff(int, char **);
void init_style();
~PairAWPMDCut() override;
void compute(int, int) override;
void settings(int, char **) override;
void coeff(int, char **) override;
void init_style() override;
void min_pointers(double **, double **);
double init_one(int, int);
void write_restart(FILE *);
void read_restart(FILE *);
virtual void write_restart_settings(FILE *);
virtual void read_restart_settings(FILE *);
double init_one(int, int) override;
void write_restart(FILE *) override;
void read_restart(FILE *) override;
void write_restart_settings(FILE *) override;
void read_restart_settings(FILE *) override;
void min_xf_pointers(int, double **, double **);
void min_xf_get(int);
void min_x_set(int);
double memory_usage();
void min_xf_pointers(int, double **, double **) override;
void min_xf_get(int) override;
void min_x_set(int) override;
double memory_usage() override;
private:
int flexible_pressure_flag;

26
src/BOCS/compute_pressure_bocs.cpp
View File

@ -233,23 +233,14 @@ double ComputePressureBocs::get_cg_p_corr(double ** grid, int basis_type,
double vCG)
{
int i = find_index(grid[0],vCG);
double correction, deltax = vCG - grid[0][i];
double deltax = vCG - grid[0][i];
if (basis_type == BASIS_LINEAR_SPLINE)
{
correction = grid[1][i] + (deltax) *
( grid[1][i+1] - grid[1][i] ) / ( grid[0][i+1] - grid[0][i] );
}
return grid[1][i] + (deltax) * ( grid[1][i+1] - grid[1][i] ) / ( grid[0][i+1] - grid[0][i] );
else if (basis_type == BASIS_CUBIC_SPLINE)
{
correction = grid[1][i] + (grid[2][i] * deltax) +
(grid[3][i] * pow(deltax,2)) + (grid[4][i] * pow(deltax,3));
}
else
{
error->all(FLERR,"bad spline type passed to get_cg_p_corr()\n");
}
return correction;
return grid[1][i] + (grid[2][i] * deltax) + (grid[3][i] * pow(deltax,2)) + (grid[4][i] * pow(deltax,3));
else error->all(FLERR,"bad spline type passed to get_cg_p_corr()\n");
return 0.0;
}
/* ----------------------------------------------------------------------
@ -261,11 +252,8 @@ void ComputePressureBocs::send_cg_info(int basis_type, int sent_N_basis,
double *sent_phi_coeff, int sent_N_mol,
double sent_vavg)
{
if (basis_type == BASIS_ANALYTIC) { p_basis_type = BASIS_ANALYTIC; }
else
{
error->all(FLERR,"Incorrect basis type passed to ComputePressureBocs\n");
}
if (basis_type == BASIS_ANALYTIC) p_basis_type = BASIS_ANALYTIC;
else error->all(FLERR,"Incorrect basis type passed to ComputePressureBocs\n");
p_match_flag = 1;

49
src/BOCS/compute_pressure_bocs.h
View File

@ -38,11 +38,11 @@ const int NUM_CUBIC_SPLINE_COLUMNS = 5; // cubic spline columns passed to co
class ComputePressureBocs : public Compute {
public:
ComputePressureBocs(class LAMMPS *, int, char **);
virtual ~ComputePressureBocs();
virtual void init();
virtual double compute_scalar();
virtual void compute_vector();
void reset_extra_compute_fix(const char *);
~ComputePressureBocs() override;
void init() override;
double compute_scalar() override;
void compute_vector() override;
void reset_extra_compute_fix(const char *) override;
double compute_cg_scalar();
double get_cg_p_corr(int, double *, int, double, double);
@ -80,42 +80,3 @@ class ComputePressureBocs : public Compute {
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Compute pressure must use group all
Virial contributions computed by potentials (pair, bond, etc) are
computed on all atoms.
E: Could not find compute pressure temperature ID
The compute ID for calculating temperature does not exist.
E: Compute pressure temperature ID does not compute temperature
The compute ID assigned to a pressure computation must compute
temperature.
E: Compute pressure requires temperature ID to include kinetic energy
The keflag cannot be used unless a temperature compute is provided.
E: Virial was not tallied on needed timestep
You are using a thermo keyword that requires potentials to
have tallied the virial, but they didn't on this timestep. See the
variable doc page for ideas on how to make this work.
E: Must use 'kspace_modify pressure/scalar no' for tensor components with kspace_style msm
Otherwise MSM will compute only a scalar pressure. See the kspace_modify
command for details on this setting.
*/

16
src/BOCS/fix_bocs.cpp
View File

@ -489,7 +489,7 @@ void FixBocs::init()
{
for (int i = 0; i < modify->nfix; i++)
if (strcmp(modify->fix[i]->style,"deform") == 0) {
int *dimflag = ((FixDeform *) modify->fix[i])->dimflag;
int *dimflag = (dynamic_cast<FixDeform *>( modify->fix[i]))->dimflag;
if ((p_flag[0] && dimflag[0]) || (p_flag[1] && dimflag[1]) ||
(p_flag[2] && dimflag[2]) || (p_flag[3] && dimflag[3]) ||
(p_flag[4] && dimflag[4]) || (p_flag[5] && dimflag[5]))
@ -523,12 +523,12 @@ void FixBocs::init()
{
if (p_basis_type == BASIS_ANALYTIC)
{
((ComputePressureBocs *)pressure)->send_cg_info(p_basis_type,
(dynamic_cast<ComputePressureBocs *>(pressure))->send_cg_info(p_basis_type,
N_p_match, p_match_coeffs, N_mol, vavg);
}
else if (p_basis_type == BASIS_LINEAR_SPLINE || p_basis_type == BASIS_CUBIC_SPLINE)
{
((ComputePressureBocs *)pressure)->send_cg_info(p_basis_type,
(dynamic_cast<ComputePressureBocs *>(pressure))->send_cg_info(p_basis_type,
splines, spline_length);
}
}
@ -589,8 +589,8 @@ void FixBocs::init()
else kspace_flag = 0;
if (utils::strmatch(update->integrate_style,"^respa")) {
nlevels_respa = ((Respa *) update->integrate)->nlevels;
step_respa = ((Respa *) update->integrate)->step;
nlevels_respa = (dynamic_cast<Respa *>( update->integrate))->nlevels;
step_respa = (dynamic_cast<Respa *>( update->integrate))->step;
dto = 0.5*step_respa[0];
}
@ -1452,7 +1452,7 @@ int FixBocs::pack_restart_data(double *list)
void FixBocs::restart(char *buf)
{
int n = 0;
double *list = (double *) buf;
auto list = (double *) buf;
int flag = static_cast<int> (list[n++]);
if (flag) {
int m = static_cast<int> (list[n++]);
@ -1551,12 +1551,12 @@ int FixBocs::modify_param(int narg, char **arg)
{
if (p_basis_type == BASIS_ANALYTIC)
{
((ComputePressureBocs *)pressure)->send_cg_info(p_basis_type, N_p_match,
(dynamic_cast<ComputePressureBocs *>(pressure))->send_cg_info(p_basis_type, N_p_match,
p_match_coeffs, N_mol, vavg);
}
else if (p_basis_type == BASIS_LINEAR_SPLINE || p_basis_type == BASIS_CUBIC_SPLINE )
{
((ComputePressureBocs *)pressure)->send_cg_info(p_basis_type, splines, spline_length );
(dynamic_cast<ComputePressureBocs *>(pressure))->send_cg_info(p_basis_type, splines, spline_length );
}
}

171
src/BOCS/fix_bocs.h
View File

@ -30,26 +30,26 @@ namespace LAMMPS_NS {
class FixBocs : public Fix {
public:
FixBocs(class LAMMPS *, int, char **); // MRD NJD
virtual ~FixBocs(); // MRD NJD
int setmask();
virtual void init();
virtual void setup(int);
virtual void initial_integrate(int);
virtual void final_integrate();
void initial_integrate_respa(int, int, int);
void pre_force_respa(int, int, int);
void final_integrate_respa(int, int);
virtual void pre_exchange();
double compute_scalar();
virtual double compute_vector(int);
void write_restart(FILE *);
~FixBocs() override; // MRD NJD
int setmask() override;
void init() override;
void setup(int) override;
void initial_integrate(int) override;
void final_integrate() override;
void initial_integrate_respa(int, int, int) override;
void pre_force_respa(int, int, int) override;
void final_integrate_respa(int, int) override;
void pre_exchange() override;
double compute_scalar() override;
double compute_vector(int) override;
void write_restart(FILE *) override;
virtual int pack_restart_data(double *); // pack restart data
virtual void restart(char *);
int modify_param(int, char **);
void reset_target(double);
void reset_dt();
virtual void *extract(const char *, int &);
double memory_usage();
void restart(char *) override;
int modify_param(int, char **) override;
void reset_target(double) override;
void reset_dt() override;
void *extract(const char *, int &) override;
double memory_usage() override;
protected:
int dimension, which;
@ -169,136 +169,3 @@ class FixBocs : public Fix {
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: CG basis type XXX is not recognized
See second line of message for supported basis types.
E: Target temperature for fix bocs cannot be 0.0
Self-explanatory.
E: Invalid fix bocs command for a 2d simulation
Cannot control z dimension in a 2d model.
E: Fix bocs dilate group ID does not exist
Self-explanatory.
E: Invalid fix bocs command pressure settings
If multiple dimensions are coupled, those dimensions must be
specified.
E: Cannot use fix bocs on a non-periodic dimension
When specifying a diagonal pressure component, the dimension must be
periodic.
E: Cannot use fix bocs on a 2nd non-periodic dimension
When specifying an off-diagonal pressure component, the 2nd of the two
dimensions must be periodic. E.g. if the xy component is specified,
then the y dimension must be periodic.
E: Cannot use fix bocs with yz scaling when z is non-periodic dimension
The 2nd dimension in the barostatted tilt factor must be periodic.
E: Cannot use fix bocs with xz scaling when z is non-periodic dimension
The 2nd dimension in the barostatted tilt factor must be periodic.
E: Cannot use fix bocs with xy scaling when y is non-periodic dimension
The 2nd dimension in the barostatted tilt factor must be periodic.
E: Cannot use fix bocs with both yz dynamics and yz scaling
Self-explanatory.
E: Cannot use fix bocs with both xz dynamics and xz scaling
Self-explanatory.
E: Cannot use fix bocs with both xy dynamics and xy scaling
Self-explanatory.
E: Can not specify Pxy/Pxz/Pyz in fix bocs with non-triclinic box
Only triclinic boxes can be used with off-diagonal pressure components.
See the region prism command for details.
E: Invalid fix bocs pressure settings
Settings for coupled dimensions must be the same.
E: Using update dipole flag requires atom style sphere
Self-explanatory.
E: Using update dipole flag requires atom attribute mu
Self-explanatory.
E: The dlm flag must be used with update dipole
Self-explanatory.
E: Fix bocs damping parameters must be > 0.0
Self-explanatory.
E: Cannot use fix npt and fix deform on same component of stress tensor
This would be changing the same box dimension twice.
E: Temperature ID for fix bocs does not exist
Self-explanatory.
E: Pressure ID for fix bocs does not exist
Self-explanatory.
E: Fix bocs has tilted box too far in one step - periodic cell is too far from equilibrium state
Self-explanatory. The change in the box tilt is too extreme
on a short timescale.
E: Could not find fix_modify temperature ID
The compute ID for computing temperature does not exist.
E: Fix_modify temperature ID does not compute temperature
The compute ID assigned to the fix must compute temperature.
W: Temperature for fix modify is not for group all
The temperature compute is being used with a pressure calculation
which does operate on group all, so this may be inconsistent.
E: Pressure ID for fix modify does not exist
Self-explanatory.
E: Could not find fix_modify pressure ID
The compute ID for computing pressure does not exist.
E: Fix_modify pressure ID does not compute pressure
The compute ID assigned to the fix must compute pressure.
*/

46
src/BODY/body_nparticle.h
View File

@ -28,21 +28,21 @@ namespace LAMMPS_NS {
class BodyNparticle : public Body {
public:
BodyNparticle(class LAMMPS *, int, char **);
~BodyNparticle();
~BodyNparticle() override;
int nsub(struct AtomVecBody::Bonus *);
double *coords(struct AtomVecBody::Bonus *);
int pack_border_body(struct AtomVecBody::Bonus *, double *);
int unpack_border_body(struct AtomVecBody::Bonus *, double *);
void data_body(int, int, int, int *, double *);
int pack_data_body(tagint, int, double *);
int write_data_body(FILE *, double *);
double radius_body(int, int, int *, double *);
int pack_border_body(struct AtomVecBody::Bonus *, double *) override;
int unpack_border_body(struct AtomVecBody::Bonus *, double *) override;
void data_body(int, int, int, int *, double *) override;
int pack_data_body(tagint, int, double *) override;
int write_data_body(FILE *, double *) override;
double radius_body(int, int, int *, double *) override;
int noutrow(int);
int noutcol();
void output(int, int, double *);
int image(int, double, double, int *&, double **&);
int noutrow(int) override;
int noutcol() override;
void output(int, int, double *) override;
int image(int, double, double, int *&, double **&) override;
private:
int *imflag;
@ -53,27 +53,3 @@ class BodyNparticle : public Body {
#endif
#endif
/* ERROR/WARNING messages:
E: Invalid body nparticle command
Arguments in atom-style command are not correct.
E: Incorrect # of integer values in Bodies section of data file
See doc page for body style.
E: Incorrect integer value in Bodies section of data file
See doc page for body style.
E: Incorrect # of floating-point values in Bodies section of data file
See doc page for body style.
E: Insufficient Jacobi rotations for body nparticle
Eigensolve for rigid body was not sufficiently accurate.
*/

51
src/BODY/body_rounded_polygon.h
View File

@ -28,7 +28,7 @@ namespace LAMMPS_NS {
class BodyRoundedPolygon : public Body {
public:
BodyRoundedPolygon(class LAMMPS *, int, char **);
~BodyRoundedPolygon();
~BodyRoundedPolygon() override;
int nsub(struct AtomVecBody::Bonus *);
double *coords(struct AtomVecBody::Bonus *);
int nedges(struct AtomVecBody::Bonus *);
@ -36,17 +36,17 @@ class BodyRoundedPolygon : public Body {
double enclosing_radius(struct AtomVecBody::Bonus *);
double rounded_radius(struct AtomVecBody::Bonus *);
int pack_border_body(struct AtomVecBody::Bonus *, double *);
int unpack_border_body(struct AtomVecBody::Bonus *, double *);
void data_body(int, int, int, int *, double *);
int pack_data_body(tagint, int, double *);
int write_data_body(FILE *, double *);
double radius_body(int, int, int *, double *);
int pack_border_body(struct AtomVecBody::Bonus *, double *) override;
int unpack_border_body(struct AtomVecBody::Bonus *, double *) override;
void data_body(int, int, int, int *, double *) override;
int pack_data_body(tagint, int, double *) override;
int write_data_body(FILE *, double *) override;
double radius_body(int, int, int *, double *) override;
int noutrow(int);
int noutcol();
void output(int, int, double *);
int image(int, double, double, int *&, double **&);
int noutrow(int) override;
int noutcol() override;
void output(int, int, double *) override;
int image(int, double, double, int *&, double **&) override;
private:
int *imflag;
@ -57,32 +57,3 @@ class BodyRoundedPolygon : public Body {
#endif
#endif
/* ERROR/WARNING messages:
E: Invalid body rounded/polygon command
Arguments in atom-style command are not correct.
E: Invalid format in Bodies section of data file
The specified number of integer or floating point values does not
appear.
E: Incorrect # of integer values in Bodies section of data file
See doc page for body style.
E: Incorrect integer value in Bodies section of data file
See doc page for body style.
E: Incorrect # of floating-point values in Bodies section of data file
See doc page for body style.
E: Insufficient Jacobi rotations for body nparticle
Eigensolve for rigid body was not sufficiently accurate.
*/

51
src/BODY/body_rounded_polyhedron.h
View File

@ -28,7 +28,7 @@ namespace LAMMPS_NS {
class BodyRoundedPolyhedron : public Body {
public:
BodyRoundedPolyhedron(class LAMMPS *, int, char **);
~BodyRoundedPolyhedron();
~BodyRoundedPolyhedron() override;
int nsub(struct AtomVecBody::Bonus *);
double *coords(struct AtomVecBody::Bonus *);
int nedges(struct AtomVecBody::Bonus *);
@ -38,17 +38,17 @@ class BodyRoundedPolyhedron : public Body {
double enclosing_radius(struct AtomVecBody::Bonus *);
double rounded_radius(struct AtomVecBody::Bonus *);
int pack_border_body(struct AtomVecBody::Bonus *, double *);
int unpack_border_body(struct AtomVecBody::Bonus *, double *);
void data_body(int, int, int, int *, double *);
int pack_data_body(tagint, int, double *);
int write_data_body(FILE *, double *);
double radius_body(int, int, int *, double *);
int pack_border_body(struct AtomVecBody::Bonus *, double *) override;
int unpack_border_body(struct AtomVecBody::Bonus *, double *) override;
void data_body(int, int, int, int *, double *) override;
int pack_data_body(tagint, int, double *) override;
int write_data_body(FILE *, double *) override;
double radius_body(int, int, int *, double *) override;
int noutrow(int);
int noutcol();
void output(int, int, double *);
int image(int, double, double, int *&, double **&);
int noutrow(int) override;
int noutcol() override;
void output(int, int, double *) override;
int image(int, double, double, int *&, double **&) override;
private:
int *imflag;
@ -59,32 +59,3 @@ class BodyRoundedPolyhedron : public Body {
#endif
#endif
/* ERROR/WARNING messages:
E: Invalid body rounded/polyhedron command
Arguments in atom-style command are not correct.
E: Invalid format in Bodies section of data file
The specified number of integer or floating point values does not
appear.
E: Incorrect # of integer values in Bodies section of data file
See doc page for body style.
E: Incorrect integer value in Bodies section of data file
See doc page for body style.
E: Incorrect # of floating-point values in Bodies section of data file
See doc page for body style.
E: Insufficient Jacobi rotations for body rounded/polyhedron
Eigensolve for rigid body was not sufficiently accurate.
*/

4
src/BODY/compute_body_local.cpp
View File

@ -53,7 +53,7 @@ ComputeBodyLocal::ComputeBodyLocal(LAMMPS *lmp, int narg, char **arg) :
}
}
avec = (AtomVecBody *) atom->style_match("body");
avec = dynamic_cast<AtomVecBody *>( atom->style_match("body"));
if (!avec) error->all(FLERR,"Compute body/local requires atom style body");
bptr = avec->bptr;
@ -153,7 +153,7 @@ int ComputeBodyLocal::compute_body(int flag)
// perform computation and fill output vector/array
int m,n,ibonus;
double *values = new double[bptr->noutcol()];
auto values = new double[bptr->noutcol()];
double **x = atom->x;
tagint *tag = atom->tag;

30
src/BODY/compute_body_local.h
View File

@ -27,10 +27,10 @@ namespace LAMMPS_NS {
class ComputeBodyLocal : public Compute {
public:
ComputeBodyLocal(class LAMMPS *, int, char **);
~ComputeBodyLocal();
void init();
void compute_local();
double memory_usage();
~ComputeBodyLocal() override;
void init() override;
void compute_local() override;
double memory_usage() override;
private:
int nvalues;
@ -49,25 +49,3 @@ class ComputeBodyLocal : public Compute {
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Compute body/local requires atom style body
Self-explanatory.
E: Invalid index in compute body/local command
Self-explanatory.
E: Invalid index for non-body particles in compute body/local command
Only indices 1,2,3 can be used for non-body particles.
*/

2
src/BODY/compute_temp_body.cpp
View File

@ -87,7 +87,7 @@ void ComputeTempBody::init()
{
// error check
avec = (AtomVecBody *) atom->style_match("body");
avec = dynamic_cast<AtomVecBody *>( atom->style_match("body"));
if (!avec)
error->all(FLERR,"Compute temp/body requires atom style body");

54
src/BODY/compute_temp_body.h
View File

@ -27,14 +27,14 @@ namespace LAMMPS_NS {
class ComputeTempBody : public Compute {
public:
ComputeTempBody(class LAMMPS *, int, char **);
~ComputeTempBody();
void init();
void setup();
double compute_scalar();
void compute_vector();
~ComputeTempBody() override;
void init() override;
void setup() override;
double compute_scalar() override;
void compute_vector() override;
void remove_bias(int, double *);
void restore_bias(int, double *);
void remove_bias(int, double *) override;
void restore_bias(int, double *) override;
private:
int mode;
@ -50,43 +50,3 @@ class ComputeTempBody : public Compute {
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Compute temp/body requires atom style body
Self-explanatory.
E: Compute temp/body requires bodies
This compute can only be applied to body particles.
E: Could not find compute ID for temperature bias
Self-explanatory.
E: Bias compute does not calculate temperature
The specified compute must compute temperature.
E: Bias compute does not calculate a velocity bias
The specified compute must compute a bias for temperature.
E: Bias compute group does not match compute group
The specified compute must operate on the same group as the parent
compute.
E: Temperature compute degrees of freedom < 0
This should not happen if you are calculating the temperature
on a valid set of atoms.
*/

2
src/BODY/fix_nh_body.cpp
View File

@ -37,7 +37,7 @@ FixNHBody::FixNHBody(LAMMPS *lmp, int narg, char **arg) :
void FixNHBody::init()
{
avec = (AtomVecBody *) atom->style_match("body");
avec = dynamic_cast<AtomVecBody *>( atom->style_match("body"));
if (!avec)
error->all(FLERR,
"Compute nvt/nph/npt body requires atom style body");

21
src/BODY/fix_nh_body.h
View File

@ -21,30 +21,17 @@ namespace LAMMPS_NS {
class FixNHBody : public FixNH {
public:
FixNHBody(class LAMMPS *, int, char **);
virtual ~FixNHBody() {}
void init();
void init() override;
protected:
double dtq;
class AtomVecBody *avec;
void nve_v();
void nve_x();
void nh_v_temp();
void nve_v() override;
void nve_x() override;
void nh_v_temp() override;
};
} // namespace LAMMPS_NS
#endif
/* ERROR/WARNING messages:
E: Compute nvt/nph/npt body requires atom style body
Self-explanatory.
E: Fix nvt/nph/npt body requires bodies
Self-explanatory.
*/

13
src/BODY/fix_nph_body.h
View File

@ -27,22 +27,9 @@ namespace LAMMPS_NS {
class FixNPHBody : public FixNHBody {
public:
FixNPHBody(class LAMMPS *, int, char **);
~FixNPHBody() {}
};
} // namespace LAMMPS_NS
#endif
#endif
/* ERROR/WARNING messages:
E: Temperature control can not be used with fix nph/body
Self-explanatory.
E: Pressure control must be used with fix nph/body
Self-explanatory.
*/

13
src/BODY/fix_npt_body.h
View File

@ -27,22 +27,9 @@ namespace LAMMPS_NS {
class FixNPTBody : public FixNHBody {
public:
FixNPTBody(class LAMMPS *, int, char **);
~FixNPTBody() {}
};
} // namespace LAMMPS_NS
#endif
#endif
/* ERROR/WARNING messages:
E: Temperature control must be used with fix npt/body
Self-explanatory.
E: Pressure control must be used with fix npt/body
Self-explanatory.
*/

2
src/BODY/fix_nve_body.cpp
View File

@ -30,7 +30,7 @@ FixNVEBody::FixNVEBody(LAMMPS *lmp, int narg, char **arg) :
void FixNVEBody::init()
{
avec = (AtomVecBody *) atom->style_match("body");
avec = dynamic_cast<AtomVecBody *>( atom->style_match("body"));
if (!avec) error->all(FLERR,"Fix nve/body requires atom style body");
// check that all particles are bodies

18
src/BODY/fix_nve_body.h
View File

@ -27,9 +27,9 @@ namespace LAMMPS_NS {
class FixNVEBody : public FixNVE {
public:
FixNVEBody(class LAMMPS *, int, char **);
void init();
void initial_integrate(int);
void final_integrate();
void init() override;
void initial_integrate(int) override;
void final_integrate() override;
private:
double dtq;
@ -39,15 +39,3 @@ class FixNVEBody : public FixNVE {
} // namespace LAMMPS_NS
#endif
#endif
/* ERROR/WARNING messages:
E: Fix nve/body requires atom style body
Self-explanatory.
E: Fix nve/body requires bodies
This fix can only be used for particles that are bodies.
*/

13
src/BODY/fix_nvt_body.h
View File

@ -27,22 +27,9 @@ namespace LAMMPS_NS {
class FixNVTBody : public FixNHBody {
public:
FixNVTBody(class LAMMPS *, int, char **);
~FixNVTBody() {}
};
} // namespace LAMMPS_NS
#endif
#endif
/* ERROR/WARNING messages:
E: Temperature control must be used with fix nvt/body
Self-explanatory.
E: Pressure control can not be used with fix nvt/body
Self-explanatory.
*/

4
src/BODY/fix_wall_body_polygon.cpp
View File

@ -179,13 +179,13 @@ void FixWallBodyPolygon::init()
{
dt = update->dt;
avec = (AtomVecBody *) atom->style_match("body");
avec = dynamic_cast<AtomVecBody *>( atom->style_match("body"));
if (!avec)
error->all(FLERR,"Pair body/rounded/polygon requires atom style body");
if (strcmp(avec->bptr->style,"rounded/polygon") != 0)
error->all(FLERR,"Pair body/rounded/polygon requires "
"body style rounded/polygon");
bptr = (BodyRoundedPolygon *) avec->bptr;
bptr = dynamic_cast<BodyRoundedPolygon *>( avec->bptr);
// set pairstyle from body/polygonular pair style

43
src/BODY/fix_wall_body_polygon.h
View File

@ -27,12 +27,12 @@ namespace LAMMPS_NS {
class FixWallBodyPolygon : public Fix {
public:
FixWallBodyPolygon(class LAMMPS *, int, char **);
virtual ~FixWallBodyPolygon();
int setmask();
void init();
void setup(int);
virtual void post_force(int);
void reset_dt();
~FixWallBodyPolygon() override;
int setmask() override;
void init() override;
void setup(int) override;
void post_force(int) override;
void reset_dt() override;
struct Contact {
int ibody, jbody; // body (i.e. atom) indices (not tags)
@ -93,34 +93,3 @@ class FixWallBodyPolygon : public Fix {
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Fix wall/body/polygon requires atom style body rounded/polygon
Self-explanatory.
E: Cannot use wall in periodic dimension
Self-explanatory.
E: Cannot wiggle and shear fix wall/body/polygon
Cannot specify both options at the same time.
E: Invalid wiggle direction for fix wall/body/polygon
Self-explanatory.
E: Fix wall/body/polygon is incompatible with Pair style
Must use a body pair style to define the parameters needed for
this fix.
*/

4
src/BODY/fix_wall_body_polyhedron.cpp
View File

@ -184,13 +184,13 @@ void FixWallBodyPolyhedron::init()
{
dt = update->dt;
avec = (AtomVecBody *) atom->style_match("body");
avec = dynamic_cast<AtomVecBody *>( atom->style_match("body"));
if (!avec)
error->all(FLERR,"Pair body/rounded/polyhedron requires atom style body");
if (strcmp(avec->bptr->style,"rounded/polyhedron") != 0)
error->all(FLERR,"Pair body/rounded/polyhedron requires "
"body style rounded/polyhedron");
bptr = (BodyRoundedPolyhedron *) avec->bptr;
bptr = dynamic_cast<BodyRoundedPolyhedron *>( avec->bptr);
// set pairstyle from body/polyhedronular pair style

43
src/BODY/fix_wall_body_polyhedron.h
View File

@ -27,12 +27,12 @@ namespace LAMMPS_NS {
class FixWallBodyPolyhedron : public Fix {
public:
FixWallBodyPolyhedron(class LAMMPS *, int, char **);
virtual ~FixWallBodyPolyhedron();
int setmask();
void init();
void setup(int);
virtual void post_force(int);
void reset_dt();
~FixWallBodyPolyhedron() override;
int setmask() override;
void init() override;
void setup(int) override;
void post_force(int) override;
void reset_dt() override;
struct Contact {
int ibody, jbody; // body (i.e. atom) indices (not tags)
@ -104,34 +104,3 @@ class FixWallBodyPolyhedron : public Fix {
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Fix wall/body/polyhedron requires atom style body rounded/polyhedron
Self-explanatory.
E: Cannot use wall in periodic dimension
Self-explanatory.
E: Cannot wiggle and shear fix wall/body/polygon
Cannot specify both options at the same time.
E: Invalid wiggle direction for fix wall/body/polygon
Self-explanatory.
E: Fix wall/body/polygon is incompatible with Pair style
Must use a body pair style to define the parameters needed for
this fix.
*/

6
src/BODY/pair_body_nparticle.cpp
View File

@ -416,13 +416,13 @@ void PairBodyNparticle::coeff(int narg, char **arg)
void PairBodyNparticle::init_style()
{
avec = (AtomVecBody *) atom->style_match("body");
avec = dynamic_cast<AtomVecBody *>( atom->style_match("body"));
if (!avec) error->all(FLERR,"Pair body/nparticle requires atom style body");
if (strcmp(avec->bptr->style,"nparticle") != 0)
error->all(FLERR,"Pair body/nparticle requires body style nparticle");
bptr = (BodyNparticle *) avec->bptr;
bptr = dynamic_cast<BodyNparticle *>( avec->bptr);
neighbor->request(this,instance_me);
neighbor->add_request(this);
}
/* ----------------------------------------------------------------------

34
src/BODY/pair_body_nparticle.h
View File

@ -27,12 +27,12 @@ namespace LAMMPS_NS {
class PairBodyNparticle : public Pair {
public:
PairBodyNparticle(class LAMMPS *);
~PairBodyNparticle();
void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
void init_style();
double init_one(int, int);
~PairBodyNparticle() override;
void compute(int, int) override;
void settings(int, char **) override;
void coeff(int, char **) override;
void init_style() override;
double init_one(int, int) override;
protected:
double cut_global;
@ -58,25 +58,3 @@ class PairBodyNparticle : public Pair {
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
E: Pair body requires atom style body
Self-explanatory.
E: Pair body requires body style nparticle
This pair style is specific to the nparticle body style.
*/

6
src/BODY/pair_body_rounded_polygon.cpp
View File

@ -411,13 +411,13 @@ void PairBodyRoundedPolygon::coeff(int narg, char **arg)
void PairBodyRoundedPolygon::init_style()
{
avec = (AtomVecBody *) atom->style_match("body");
avec = dynamic_cast<AtomVecBody *>( atom->style_match("body"));
if (!avec)
error->all(FLERR,"Pair body/rounded/polygon requires atom style body");
if (strcmp(avec->bptr->style,"rounded/polygon") != 0)
error->all(FLERR,"Pair body/rounded/polygon requires "
"body style rounded/polygon");
bptr = (BodyRoundedPolygon *) avec->bptr;
bptr = dynamic_cast<BodyRoundedPolygon *>( avec->bptr);
if (force->newton_pair == 0)
error->all(FLERR,"Pair style body/rounded/polygon requires "
@ -427,7 +427,7 @@ void PairBodyRoundedPolygon::init_style()
error->all(FLERR,"Pair body/rounded/polygon requires "
"ghost atoms store velocity");
neighbor->request(this);
neighbor->add_request(this);
// find the maximum enclosing radius for each atom type

34
src/BODY/pair_body_rounded_polygon.h
View File

@ -27,12 +27,12 @@ namespace LAMMPS_NS {
class PairBodyRoundedPolygon : public Pair {
public:
PairBodyRoundedPolygon(class LAMMPS *);
~PairBodyRoundedPolygon();
void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
void init_style();
double init_one(int, int);
~PairBodyRoundedPolygon() override;
void compute(int, int) override;
void settings(int, char **) override;
void coeff(int, char **) override;
void init_style() override;
double init_one(int, int) override;
struct Contact {
int ibody, jbody; // body (i.e. atom) indices (not tags)
@ -108,25 +108,3 @@ class PairBodyRoundedPolygon : public Pair {
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
E: Pair body/rounded/polygon requires atom style body rounded/polygon
Self-explanatory.
E: Pair body requires body style rounded/polygon
This pair style is specific to the rounded/polygon body style.
*/

6
src/BODY/pair_body_rounded_polyhedron.cpp
View File

@ -390,13 +390,13 @@ void PairBodyRoundedPolyhedron::coeff(int narg, char **arg)
void PairBodyRoundedPolyhedron::init_style()
{
avec = (AtomVecBody *) atom->style_match("body");
avec = dynamic_cast<AtomVecBody *>( atom->style_match("body"));
if (!avec) error->all(FLERR,"Pair body/rounded/polyhedron requires "
"atom style body");
if (strcmp(avec->bptr->style,"rounded/polyhedron") != 0)
error->all(FLERR,"Pair body/rounded/polyhedron requires "
"body style rounded/polyhedron");
bptr = (BodyRoundedPolyhedron *) avec->bptr;
bptr = dynamic_cast<BodyRoundedPolyhedron *>( avec->bptr);
if (force->newton_pair == 0)
error->all(FLERR,"Pair style body/rounded/polyhedron requires "
@ -406,7 +406,7 @@ void PairBodyRoundedPolyhedron::init_style()
error->all(FLERR,"Pair body/rounded/polyhedron requires "
"ghost atoms store velocity");
neighbor->request(this);
neighbor->add_request(this);
// find the maximum enclosing radius for each atom type

34
src/BODY/pair_body_rounded_polyhedron.h
View File

@ -27,12 +27,12 @@ namespace LAMMPS_NS {
class PairBodyRoundedPolyhedron : public Pair {
public:
PairBodyRoundedPolyhedron(class LAMMPS *);
~PairBodyRoundedPolyhedron();
void compute(int, int);
void settings(int, char **);
void coeff(int, char **);
void init_style();
double init_one(int, int);
~PairBodyRoundedPolyhedron() override;
void compute(int, int) override;
void settings(int, char **) override;
void coeff(int, char **) override;
void init_style() override;
double init_one(int, int) override;
virtual void kernel_force(double R, int itype, int jtype, double &energy, double &fpair);
@ -162,25 +162,3 @@ class PairBodyRoundedPolyhedron : public Pair {
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal ... command
Self-explanatory. Check the input script syntax and compare to the
documentation for the command. You can use -echo screen as a
command-line option when running LAMMPS to see the offending line.
E: Incorrect args for pair coefficients
Self-explanatory. Check the input script or data file.
E: Pair body/rounded/polyhedron requires atom style body rounded/polyhedron
Self-explanatory.
E: Pair body requires body style rounded/polyhedron
This pair style is specific to the rounded/polyhedron body style.
*/

248
src/BPM/atom_vec_bpm_sphere.cpp Normal file
View File

@ -0,0 +1,248 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "atom_vec_bpm_sphere.h"
#include "atom.h"
#include "comm.h"
#include "error.h"
#include "fix.h"
#include "fix_adapt.h"
#include "math_const.h"
#include "modify.h"
#include "utils.h"
#include <cstring>
using namespace LAMMPS_NS;
using MathConst::MY_PI;
/* ---------------------------------------------------------------------- */
AtomVecBPMSphere::AtomVecBPMSphere(LAMMPS *_lmp) : AtomVec(_lmp)
{
mass_type = PER_ATOM;
molecular = Atom::MOLECULAR;
bonds_allow = 1;
atom->molecule_flag = 1;
atom->sphere_flag = 1;
atom->radius_flag = atom->rmass_flag = atom->omega_flag = atom->torque_flag = atom->quat_flag = 1;
// strings with peratom variables to include in each AtomVec method
// strings cannot contain fields in corresponding AtomVec default strings
// order of fields in a string does not matter
// except: fields_data_atom & fields_data_vel must match data file
fields_grow = {"molecule", "num_bond", "bond_type", "bond_atom", "nspecial", "special",
"radius", "rmass", "omega", "torque", "quat"};
fields_copy = {"molecule", "num_bond", "bond_type", "bond_atom", "nspecial",
"special", "radius", "rmass", "omega", "quat"};
fields_comm_vel = {"omega", "quat"};
fields_reverse = {"torque"};
fields_border = {"molecule", "radius", "rmass"};
fields_border_vel = {"molecule", "radius", "rmass", "omega", "quat"};
fields_exchange = {"molecule", "num_bond", "bond_type", "bond_atom", "nspecial",
"special", "radius", "rmass", "omega", "quat"};
fields_restart = {"molecule", "num_bond", "bond_type", "bond_atom",
"radius", "rmass", "omega", "quat"};
fields_create = {"molecule", "num_bond", "nspecial", "radius", "rmass", "omega", "quat"};
fields_data_atom = {"id", "molecule", "type", "radius", "rmass", "x"};
fields_data_vel = {"id", "v", "omega"};
bond_per_atom = 0;
bond_negative = NULL;
}
/* ----------------------------------------------------------------------
process sub-style args
optional arg = 0/1 for static/dynamic particle radii
------------------------------------------------------------------------- */
void AtomVecBPMSphere::process_args(int narg, char **arg)
{
if (narg != 0 && narg != 1) error->all(FLERR, "Illegal atom_style bpm/sphere command");
radvary = 0;
if (narg == 1) {
radvary = utils::numeric(FLERR, arg[0], true, lmp);
if (radvary < 0 || radvary > 1) error->all(FLERR, "Illegal atom_style bpm/sphere command");
}
// dynamic particle radius and mass must be communicated every step
if (radvary) {
fields_comm = {"radius", "rmass"};
fields_comm_vel = {"radius", "rmass", "omega"};
}
// delay setting up of fields until now
setup_fields();
}
/* ---------------------------------------------------------------------- */
void AtomVecBPMSphere::init()
{
AtomVec::init();
// check if optional radvary setting should have been set to 1
for (auto ifix : modify->get_fix_by_style("^adapt")) {
if (radvary == 0) {
if ((strcmp(ifix->style, "adapt") == 0) && (dynamic_cast<FixAdapt *>(ifix)->diamflag))
error->all(FLERR, "Fix adapt changes atom radii but atom_style bpm/sphere is not dynamic");
// cannot properly check for fix adapt/fep since its header is optional
if ((strcmp(ifix->style, "adapt/fep") == 0) && (comm->me == 0))
error->warning(
FLERR, "Fix adapt/fep may change atom radii but atom_style bpm/sphere is not dynamic");
}
}
}
/* ----------------------------------------------------------------------
set local copies of all grow ptrs used by this class, except defaults
needed in replicate when 2 atom classes exist and it calls pack_restart()
------------------------------------------------------------------------- */
void AtomVecBPMSphere::grow_pointers()
{
radius = atom->radius;
rmass = atom->rmass;
omega = atom->omega;
quat = atom->quat;
num_bond = atom->num_bond;
bond_type = atom->bond_type;
nspecial = atom->nspecial;
}
/* ----------------------------------------------------------------------
initialize non-zero atom quantities
------------------------------------------------------------------------- */
void AtomVecBPMSphere::create_atom_post(int ilocal)
{
radius[ilocal] = 0.5;
rmass[ilocal] = 4.0 * MY_PI / 3.0 * 0.5 * 0.5 * 0.5;
quat[ilocal][0] = 1.0;
quat[ilocal][1] = 0.0;
quat[ilocal][2] = 0.0;
quat[ilocal][3] = 0.0;
}
/* ----------------------------------------------------------------------
modify values for AtomVec::pack_restart() to pack
------------------------------------------------------------------------- */
void AtomVecBPMSphere::pack_restart_pre(int ilocal)
{
// insure bond_negative vector is needed length
if (bond_per_atom < atom->bond_per_atom) {
delete[] bond_negative;
bond_per_atom = atom->bond_per_atom;
bond_negative = new int[bond_per_atom];
}
// flip any negative types to positive and flag which ones
any_bond_negative = 0;
for (int m = 0; m < num_bond[ilocal]; m++) {
if (bond_type[ilocal][m] < 0) {
bond_negative[m] = 1;
bond_type[ilocal][m] = -bond_type[ilocal][m];
any_bond_negative = 1;
} else
bond_negative[m] = 0;
}
}
/* ----------------------------------------------------------------------
unmodify values packed by AtomVec::pack_restart()
------------------------------------------------------------------------- */
void AtomVecBPMSphere::pack_restart_post(int ilocal)
{
// restore the flagged types to their negative values
if (any_bond_negative) {
for (int m = 0; m < num_bond[ilocal]; m++)
if (bond_negative[m]) bond_type[ilocal][m] = -bond_type[ilocal][m];
}
}
/* ----------------------------------------------------------------------
initialize other atom quantities after AtomVec::unpack_restart()
------------------------------------------------------------------------- */
void AtomVecBPMSphere::unpack_restart_init(int ilocal)
{
nspecial[ilocal][0] = 0;
nspecial[ilocal][1] = 0;
nspecial[ilocal][2] = 0;
}
/* ----------------------------------------------------------------------
modify what AtomVec::data_atom() just unpacked
or initialize other atom quantities
------------------------------------------------------------------------- */
void AtomVecBPMSphere::data_atom_post(int ilocal)
{
radius_one = 0.5 * atom->radius[ilocal];
radius[ilocal] = radius_one;
if (radius_one > 0.0) rmass[ilocal] *= 4.0 * MY_PI / 3.0 * radius_one * radius_one * radius_one;
if (rmass[ilocal] <= 0.0) error->one(FLERR, "Invalid density in Atoms section of data file");
omega[ilocal][0] = 0.0;
omega[ilocal][1] = 0.0;
omega[ilocal][2] = 0.0;
quat[ilocal][0] = 1.0;
quat[ilocal][1] = 0.0;
quat[ilocal][2] = 0.0;
quat[ilocal][3] = 0.0;
num_bond[ilocal] = 0;
nspecial[ilocal][0] = 0;
nspecial[ilocal][1] = 0;
nspecial[ilocal][2] = 0;
}
/* ----------------------------------------------------------------------
modify values for AtomVec::pack_data() to pack
------------------------------------------------------------------------- */
void AtomVecBPMSphere::pack_data_pre(int ilocal)
{
radius_one = radius[ilocal];
rmass_one = rmass[ilocal];
radius[ilocal] *= 2.0;
if (radius_one != 0.0)
rmass[ilocal] = rmass_one / (4.0 * MY_PI / 3.0 * radius_one * radius_one * radius_one);
}
/* ----------------------------------------------------------------------
unmodify values packed by AtomVec::pack_data()
------------------------------------------------------------------------- */
void AtomVecBPMSphere::pack_data_post(int ilocal)
{
radius[ilocal] = radius_one;
rmass[ilocal] = rmass_one;
}

61
src/BPM/atom_vec_bpm_sphere.h Normal file
View File

@ -0,0 +1,61 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef ATOM_CLASS
// clang-format off
AtomStyle(bpm/sphere,AtomVecBPMSphere);
// clang-format on
#else
#ifndef LMP_ATOM_VEC_BPM_SPHERE_H
#define LMP_ATOM_VEC_BPM_SPHERE_H
#include "atom_vec.h"
namespace LAMMPS_NS {
class AtomVecBPMSphere : public AtomVec {
public:
AtomVecBPMSphere(class LAMMPS *);
void process_args(int, char **) override;
void init() override;
void grow_pointers() override;
void create_atom_post(int) override;
void pack_restart_pre(int) override;
void pack_restart_post(int) override;
void unpack_restart_init(int) override;
void data_atom_post(int) override;
void pack_data_pre(int) override;
void pack_data_post(int) override;
private:
int *num_bond;
int **bond_type;
int **nspecial;
double *radius, *rmass;
double **omega, **quat;
int any_bond_negative;
int bond_per_atom;
int *bond_negative;
int radvary;
double radius_one, rmass_one;
};
} // namespace LAMMPS_NS
#endif
#endif

415
src/BPM/bond_bpm.cpp Normal file
View File

@ -0,0 +1,415 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "bond_bpm.h"
#include "atom.h"
#include "domain.h"
#include "error.h"
#include "fix_bond_history.h"
#include "fix_store_local.h"
#include "fix_update_special_bonds.h"
#include "force.h"
#include "memory.h"
#include "modify.h"
#include "update.h"
#include <vector>
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
BondBPM::BondBPM(LAMMPS *_lmp) :
Bond(_lmp), id_fix_dummy(nullptr), id_fix_dummy2(nullptr), id_fix_update(nullptr),
id_fix_bond_history(nullptr), id_fix_store_local(nullptr), id_fix_prop_atom(nullptr),
fix_store_local(nullptr), fix_bond_history(nullptr), fix_update_special_bonds(nullptr),
pack_choice(nullptr), output_data(nullptr)
{
overlay_flag = 0;
prop_atom_flag = 0;
nvalues = 0;
r0_max_estimate = 0.0;
max_stretch = 1.0;
// create dummy fix as placeholder for FixUpdateSpecialBonds & BondHistory
// this is so final order of Modify:fix will conform to input script
// BondHistory technically only needs this if updateflag = 1
id_fix_dummy = utils::strdup("BPM_DUMMY");
modify->add_fix(fmt::format("{} all DUMMY ", id_fix_dummy));
id_fix_dummy2 = utils::strdup("BPM_DUMMY2");
modify->add_fix(fmt::format("{} all DUMMY ", id_fix_dummy2));
}
/* ---------------------------------------------------------------------- */
BondBPM::~BondBPM()
{
delete[] pack_choice;
if (id_fix_dummy) modify->delete_fix(id_fix_dummy);
if (id_fix_dummy2) modify->delete_fix(id_fix_dummy2);
if (id_fix_update) modify->delete_fix(id_fix_update);
if (id_fix_bond_history) modify->delete_fix(id_fix_bond_history);
if (id_fix_store_local) modify->delete_fix(id_fix_store_local);
if (id_fix_prop_atom) modify->delete_fix(id_fix_prop_atom);
delete[] id_fix_dummy;
delete[] id_fix_dummy2;
delete[] id_fix_update;
delete[] id_fix_bond_history;
delete[] id_fix_store_local;
delete[] id_fix_prop_atom;
memory->destroy(output_data);
}
/* ---------------------------------------------------------------------- */
void BondBPM::init_style()
{
if (id_fix_store_local) {
auto ifix = modify->get_fix_by_id(id_fix_store_local);
if (!ifix) error->all(FLERR, "Cannot find fix store/local");
if (strcmp(ifix->style, "STORE_LOCAL") != 0)
error->all(FLERR, "Incorrect fix style matched, not store/local");
fix_store_local = dynamic_cast<FixStoreLocal *>(ifix);
fix_store_local->nvalues = nvalues;
}
if (overlay_flag) {
if (force->special_lj[1] != 1.0)
error->all(FLERR,
"With overlay/pair, BPM bond styles require special_bonds weight of 1.0 for "
"first neighbors");
if (id_fix_update) {
modify->delete_fix(id_fix_update);
delete[] id_fix_update;
id_fix_update = nullptr;
}
} else {
// Require atoms know about all of their bonds and if they break
if (force->newton_bond)
error->all(FLERR, "Without overlay/pair, BPM bond styles require Newton bond off");
// special lj must be 0 1 1 to censor pair forces between bonded particles
// special coulomb must be 1 1 1 to ensure all pairs are included in the
// neighbor list and 1-3 and 1-4 special bond lists are skipped
if (force->special_lj[1] != 0.0 || force->special_lj[2] != 1.0 || force->special_lj[3] != 1.0)
error->all(FLERR,
"Without overlay/pair, BPM bond sytles requires special LJ weights = 0,1,1");
if (force->special_coul[1] != 1.0 || force->special_coul[2] != 1.0 ||
force->special_coul[3] != 1.0)
error->all(FLERR,
"Without overlay/pair, BPM bond sytles requires special Coulomb weights = 1,1,1");
if (id_fix_dummy) {
id_fix_update = utils::strdup("BPM_UPDATE_SPECIAL_BONDS");
fix_update_special_bonds = dynamic_cast<FixUpdateSpecialBonds *>(modify->replace_fix(
id_fix_dummy, fmt::format("{} all UPDATE_SPECIAL_BONDS", id_fix_update), 1));
delete[] id_fix_dummy;
id_fix_dummy = nullptr;
}
}
if (force->angle || force->dihedral || force->improper)
error->all(FLERR, "Bond style bpm cannot be used with 3,4-body interactions");
if (atom->molecular == 2)
error->all(FLERR, "Bond style bpm cannot be used with atom style template");
// special 1-3 and 1-4 weights must be 1 to prevent building 1-3 and 1-4 special bond lists
if (force->special_lj[2] != 1.0 || force->special_lj[3] != 1.0 || force->special_coul[2] != 1.0 ||
force->special_coul[3] != 1.0)
error->all(FLERR, "Bond style bpm requires 1-3 and 1-4 special weights of 1.0");
}
/* ----------------------------------------------------------------------
global settings
All args before store/local command are saved for potential args
for specific bond BPM substyles
All args after optional store/local command are variables stored
in the compute store/local
------------------------------------------------------------------------- */
void BondBPM::settings(int narg, char **arg)
{
leftover_iarg.clear();
int iarg = 0;
while (iarg < narg) {
if (strcmp(arg[iarg], "store/local") == 0) {
nvalues = 0;
id_fix_store_local = utils::strdup(arg[iarg + 1]);
store_local_freq = utils::inumeric(FLERR, arg[iarg + 2], false, lmp);
pack_choice = new FnPtrPack[narg - iarg - 1];
iarg += 3;
while (iarg < narg) {
if (strcmp(arg[iarg], "id1") == 0) {
pack_choice[nvalues++] = &BondBPM::pack_id1;
} else if (strcmp(arg[iarg], "id2") == 0) {
pack_choice[nvalues++] = &BondBPM::pack_id2;
} else if (strcmp(arg[iarg], "time") == 0) {
pack_choice[nvalues++] = &BondBPM::pack_time;
} else if (strcmp(arg[iarg], "x") == 0) {
pack_choice[nvalues++] = &BondBPM::pack_x;
} else if (strcmp(arg[iarg], "y") == 0) {
pack_choice[nvalues++] = &BondBPM::pack_y;
} else if (strcmp(arg[iarg], "z") == 0) {
pack_choice[nvalues++] = &BondBPM::pack_z;
} else if (strcmp(arg[iarg], "x/ref") == 0) {
pack_choice[nvalues++] = &BondBPM::pack_x_ref;
prop_atom_flag = 1;
} else if (strcmp(arg[iarg], "y/ref") == 0) {
pack_choice[nvalues++] = &BondBPM::pack_y_ref;
prop_atom_flag = 1;
} else if (strcmp(arg[iarg], "z/ref") == 0) {
pack_choice[nvalues++] = &BondBPM::pack_z_ref;
prop_atom_flag = 1;
} else {
break;
}
iarg++;
}
} else if (strcmp(arg[iarg], "overlay/pair") == 0) {
overlay_flag = 1;
iarg++;
} else {
leftover_iarg.push_back(iarg);
iarg++;
}
}
if (id_fix_store_local) {
if (nvalues == 0)
error->all(FLERR, "Storing local data must include at least one value to output");
memory->create(output_data, nvalues, "bond/bpm:output_data");
auto ifix = modify->get_fix_by_id(id_fix_store_local);
if (!ifix)
ifix = modify->add_fix(
fmt::format("{} all STORE_LOCAL {} {}", id_fix_store_local, store_local_freq, nvalues));
fix_store_local = dynamic_cast<FixStoreLocal *>(ifix);
// Use property/atom to save reference positions as it can transfer to ghost atoms
// This won't work for instances where bonds are added (e.g. fix pour) but in those cases
// a reference state isn't well defined
if (prop_atom_flag == 1) {
id_fix_prop_atom = utils::strdup("BPM_property_atom");
char *x_ref_id = utils::strdup("BPM_X_REF");
char *y_ref_id = utils::strdup("BPM_Y_REF");
char *z_ref_id = utils::strdup("BPM_Z_REF");
ifix = modify->get_fix_by_id(id_fix_prop_atom);
if (!ifix)
ifix = modify->add_fix(fmt::format("{} all property/atom {} {} {} ghost yes",
id_fix_prop_atom, x_ref_id, y_ref_id, z_ref_id));
int type_flag;
int col_flag;
index_x_ref = atom->find_custom(x_ref_id, type_flag, col_flag);
index_y_ref = atom->find_custom(y_ref_id, type_flag, col_flag);
index_z_ref = atom->find_custom(z_ref_id, type_flag, col_flag);
delete[] x_ref_id;
delete[] y_ref_id;
delete[] z_ref_id;
if (ifix->restart_reset) {
ifix->restart_reset = 0;
} else {
double *x_ref = atom->dvector[index_x_ref];
double *y_ref = atom->dvector[index_y_ref];
double *z_ref = atom->dvector[index_z_ref];
double **x = atom->x;
for (int i = 0; i < atom->nlocal; i++) {
x_ref[i] = x[i][0];
y_ref[i] = x[i][1];
z_ref[i] = x[i][2];
}
}
}
}
}
/* ----------------------------------------------------------------------
used to check bond communiction cutoff - not perfect, estimates based on local-local only
------------------------------------------------------------------------- */
double BondBPM::equilibrium_distance(int /*i*/)
{
// Ghost atoms may not yet be communicated, this may only be an estimate
if (r0_max_estimate == 0) {
int type, j;
double delx, dely, delz, r;
double **x = atom->x;
for (int i = 0; i < atom->nlocal; i++) {
for (int m = 0; m < atom->num_bond[i]; m++) {
type = atom->bond_type[i][m];
if (type == 0) continue;
j = atom->map(atom->bond_atom[i][m]);
if (j == -1) continue;
delx = x[i][0] - x[j][0];
dely = x[i][1] - x[j][1];
delz = x[i][2] - x[j][2];
domain->minimum_image(delx, dely, delz);
r = sqrt(delx * delx + dely * dely + delz * delz);
if (r > r0_max_estimate) r0_max_estimate = r;
}
}
double temp;
MPI_Allreduce(&r0_max_estimate, &temp, 1, MPI_DOUBLE, MPI_MAX, world);
r0_max_estimate = temp;
}
// Divide out heuristic prefactor added in comm class
return max_stretch * r0_max_estimate / 1.5;
}
/* ---------------------------------------------------------------------- */
void BondBPM::process_broken(int i, int j)
{
if (fix_store_local) {
for (int n = 0; n < nvalues; n++) (this->*pack_choice[n])(n, i, j);
fix_store_local->add_data(output_data, i, j);
}
if (fix_update_special_bonds) fix_update_special_bonds->add_broken_bond(i, j);
// Manually search and remove from atom arrays
// need to remove in case special bonds arrays rebuilt
int m, n;
int nlocal = atom->nlocal;
tagint *tag = atom->tag;
tagint **bond_atom = atom->bond_atom;
int **bond_type = atom->bond_type;
int *num_bond = atom->num_bond;
if (i < nlocal) {
for (m = 0; m < num_bond[i]; m++) {
if (bond_atom[i][m] == tag[j]) {
bond_type[i][m] = 0;
n = num_bond[i];
bond_type[i][m] = bond_type[i][n - 1];
bond_atom[i][m] = bond_atom[i][n - 1];
fix_bond_history->shift_history(i, m, n - 1);
fix_bond_history->delete_history(i, n - 1);
num_bond[i]--;
break;
}
}
}
if (j < nlocal) {
for (m = 0; m < num_bond[j]; m++) {
if (bond_atom[j][m] == tag[i]) {
bond_type[j][m] = 0;
n = num_bond[j];
bond_type[j][m] = bond_type[j][n - 1];
bond_atom[j][m] = bond_atom[j][n - 1];
fix_bond_history->shift_history(j, m, n - 1);
fix_bond_history->delete_history(j, n - 1);
num_bond[j]--;
break;
}
}
}
}
/* ----------------------------------------------------------------------
one method for every keyword bond bpm can output
the atom property is packed into array or vector
------------------------------------------------------------------------- */
void BondBPM::pack_id1(int n, int i, int /*j*/)
{
tagint *tag = atom->tag;
output_data[n] = tag[i];
}
/* ---------------------------------------------------------------------- */
void BondBPM::pack_id2(int n, int /*i*/, int j)
{
tagint *tag = atom->tag;
output_data[n] = tag[j];
}
/* ---------------------------------------------------------------------- */
void BondBPM::pack_time(int n, int /*i*/, int /*j*/)
{
bigint time = update->ntimestep;
output_data[n] = time;
}
/* ---------------------------------------------------------------------- */
void BondBPM::pack_x(int n, int i, int j)
{
double **x = atom->x;
output_data[n] = (x[i][0] + x[j][0]) * 0.5;
}
/* ---------------------------------------------------------------------- */
void BondBPM::pack_y(int n, int i, int j)
{
double **x = atom->x;
output_data[n] = (x[i][1] + x[j][1]) * 0.5;
}
/* ---------------------------------------------------------------------- */
void BondBPM::pack_z(int n, int i, int j)
{
double **x = atom->x;
output_data[n] = (x[i][2] + x[j][2]) * 0.5;
}
/* ---------------------------------------------------------------------- */
void BondBPM::pack_x_ref(int n, int i, int j)
{
double *x = atom->dvector[index_x_ref];
output_data[n] = (x[i] + x[j]) * 0.5;
}
/* ---------------------------------------------------------------------- */
void BondBPM::pack_y_ref(int n, int i, int j)
{
double *y = atom->dvector[index_y_ref];
output_data[n] = (y[i] + y[j]) * 0.5;
}
/* ---------------------------------------------------------------------- */
void BondBPM::pack_z_ref(int n, int i, int j)
{
double *z = atom->dvector[index_z_ref];
output_data[n] = (z[i] + z[j]) * 0.5;
}

72
src/BPM/bond_bpm.h Normal file
View File

@ -0,0 +1,72 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifndef LMP_BOND_BPM_H
#define LMP_BOND_BPM_H
#include "bond.h"
#include <vector>
namespace LAMMPS_NS {
class BondBPM : public Bond {
public:
BondBPM(class LAMMPS *);
~BondBPM() override;
void compute(int, int) override = 0;
void coeff(int, char **) override = 0;
void init_style() override;
void settings(int, char **) override;
double equilibrium_distance(int) override;
void write_restart(FILE *) override{};
void read_restart(FILE *) override{};
void write_data(FILE *) override{};
double single(int, double, int, int, double &) override = 0;
protected:
double r0_max_estimate;
double max_stretch;
int store_local_freq;
std::vector<int> leftover_iarg;
char *id_fix_dummy, *id_fix_dummy2;
char *id_fix_update, *id_fix_bond_history;
char *id_fix_store_local, *id_fix_prop_atom;
class FixStoreLocal *fix_store_local;
class FixBondHistory *fix_bond_history;
class FixUpdateSpecialBonds *fix_update_special_bonds;
void process_broken(int, int);
typedef void (BondBPM::*FnPtrPack)(int, int, int);
FnPtrPack *pack_choice; // ptrs to pack functions
double *output_data;
int prop_atom_flag, nvalues, overlay_flag;
int index_x_ref, index_y_ref, index_z_ref;
void pack_id1(int, int, int);
void pack_id2(int, int, int);
void pack_time(int, int, int);
void pack_x(int, int, int);
void pack_y(int, int, int);
void pack_z(int, int, int);
void pack_x_ref(int, int, int);
void pack_y_ref(int, int, int);
void pack_z_ref(int, int, int);
};
} // namespace LAMMPS_NS
#endif

796
src/BPM/bond_bpm_rotational.cpp Normal file
View File

@ -0,0 +1,796 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "bond_bpm_rotational.h"
#include "atom.h"
#include "comm.h"
#include "domain.h"
#include "error.h"
#include "fix_bond_history.h"
#include "force.h"
#include "math_const.h"
#include "math_extra.h"
#include "memory.h"
#include "modify.h"
#include "neighbor.h"
#include <cmath>
#define EPSILON 1e-10
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
BondBPMRotational::BondBPMRotational(LAMMPS *_lmp) : BondBPM(_lmp)
{
Kr = nullptr;
Ks = nullptr;
Kt = nullptr;
Kb = nullptr;
Fcr = nullptr;
Fcs = nullptr;
Tct = nullptr;
Tcb = nullptr;
gnorm = nullptr;
gslide = nullptr;
groll = nullptr;
gtwist = nullptr;
partial_flag = 1;
smooth_flag = 1;
}
/* ---------------------------------------------------------------------- */
BondBPMRotational::~BondBPMRotational()
{
if (allocated) {
memory->destroy(setflag);
memory->destroy(Kr);
memory->destroy(Ks);
memory->destroy(Kt);
memory->destroy(Kb);
memory->destroy(Fcr);
memory->destroy(Fcs);
memory->destroy(Tct);
memory->destroy(Tcb);
memory->destroy(gnorm);
memory->destroy(gslide);
memory->destroy(groll);
memory->destroy(gtwist);
}
}
/* ---------------------------------------------------------------------- */
double BondBPMRotational::acos_limit(double c)
{
if (c > 1.0) c = 1.0;
if (c < -1.0) c = -1.0;
return acos(c);
}
/* ----------------------------------------------------------------------
Store data for a single bond - if bond added after LAMMPS init (e.g. pour)
------------------------------------------------------------------------- */
double BondBPMRotational::store_bond(int n, int i, int j)
{
double delx, dely, delz, r, rinv;
double **x = atom->x;
tagint *tag = atom->tag;
double **bondstore = fix_bond_history->bondstore;
if (tag[i] < tag[j]) {
delx = x[i][0] - x[j][0];
dely = x[i][1] - x[j][1];
delz = x[i][2] - x[j][2];
} else {
delx = x[j][0] - x[i][0];
dely = x[j][1] - x[i][1];
delz = x[j][2] - x[i][2];
}
r = sqrt(delx * delx + dely * dely + delz * delz);
rinv = 1.0 / r;
bondstore[n][0] = r;
bondstore[n][1] = delx * rinv;
bondstore[n][2] = dely * rinv;
bondstore[n][3] = delz * rinv;
if (i < atom->nlocal) {
for (int m = 0; m < atom->num_bond[i]; m++) {
if (atom->bond_atom[i][m] == tag[j]) {
fix_bond_history->update_atom_value(i, m, 0, r);
fix_bond_history->update_atom_value(i, m, 1, delx * rinv);
fix_bond_history->update_atom_value(i, m, 2, dely * rinv);
fix_bond_history->update_atom_value(i, m, 3, delz * rinv);
}
}
}
if (j < atom->nlocal) {
for (int m = 0; m < atom->num_bond[j]; m++) {
if (atom->bond_atom[j][m] == tag[i]) {
fix_bond_history->update_atom_value(j, m, 0, r);
fix_bond_history->update_atom_value(j, m, 1, delx * rinv);
fix_bond_history->update_atom_value(j, m, 2, dely * rinv);
fix_bond_history->update_atom_value(j, m, 3, delz * rinv);
}
}
}
return r;
}
/* ----------------------------------------------------------------------
Store data for all bonds called once
------------------------------------------------------------------------- */
void BondBPMRotational::store_data()
{
int i, j, m, type;
double delx, dely, delz, r, rinv;
double **x = atom->x;
int **bond_type = atom->bond_type;
tagint *tag = atom->tag;
for (i = 0; i < atom->nlocal; i++) {
for (m = 0; m < atom->num_bond[i]; m++) {
type = bond_type[i][m];
//Skip if bond was turned off
if (type < 0) continue;
// map to find index n for tag
j = atom->map(atom->bond_atom[i][m]);
if (j == -1) error->one(FLERR, "Atom missing in BPM bond");
// Save orientation as pointing towards small tag
if (tag[i] < tag[j]) {
delx = x[i][0] - x[j][0];
dely = x[i][1] - x[j][1];
delz = x[i][2] - x[j][2];
} else {
delx = x[j][0] - x[i][0];
dely = x[j][1] - x[i][1];
delz = x[j][2] - x[i][2];
}
// Get closest image in case bonded with ghost
domain->minimum_image(delx, dely, delz);
r = sqrt(delx * delx + dely * dely + delz * delz);
rinv = 1.0 / r;
fix_bond_history->update_atom_value(i, m, 0, r);
fix_bond_history->update_atom_value(i, m, 1, delx * rinv);
fix_bond_history->update_atom_value(i, m, 2, dely * rinv);
fix_bond_history->update_atom_value(i, m, 3, delz * rinv);
}
}
fix_bond_history->post_neighbor();
}
/* ----------------------------------------------------------------------
Calculate forces using formulation in:
1) Y. Wang Acta Geotechnica 2009
2) P. Mora & Y. Wang Advances in Geomcomputing 2009
---------------------------------------------------------------------- */
double BondBPMRotational::elastic_forces(int i1, int i2, int type, double &Fr, double r_mag,
double r0_mag, double r_mag_inv, double * /*rhat*/,
double *r, double *r0, double *force1on2,
double *torque1on2, double *torque2on1)
{
double breaking, temp, r0_dot_rb, c, gamma;
double psi, theta, cos_phi, sin_phi;
double mag_in_plane, mag_out_plane;
double Fs_mag, Tt_mag, Tb_mag;
double q1[4], q2[4];
double q2inv[4], mq[4], mqinv[4], qp21[4], q21[4], qtmp[4];
double rb[3], rb_x_r0[3], s[3], t[3];
double Fs[3], Fsp[3], F_rot[3], Ftmp[3];
double Ts[3], Tb[3], Tt[3], Tbp[3], Ttp[3], Tsp[3], T_rot[3], Ttmp[3];
double **quat = atom->quat;
q1[0] = quat[i1][0];
q1[1] = quat[i1][1];
q1[2] = quat[i1][2];
q1[3] = quat[i1][3];
q2[0] = quat[i2][0];
q2[1] = quat[i2][1];
q2[2] = quat[i2][2];
q2[3] = quat[i2][3];
// Calculate normal forces, rb = bond vector in particle 1's frame
MathExtra::qconjugate(q2, q2inv);
MathExtra::quatrotvec(q2inv, r, rb);
Fr = Kr[type] * (r_mag - r0_mag);
MathExtra::scale3(Fr * r_mag_inv, rb, F_rot);
// Calculate forces due to tangential displacements (no rotation)
r0_dot_rb = MathExtra::dot3(r0, rb);
c = r0_dot_rb * r_mag_inv / r0_mag;
gamma = acos_limit(c);
MathExtra::cross3(rb, r0, rb_x_r0);
MathExtra::cross3(rb, rb_x_r0, s);
MathExtra::norm3(s);
MathExtra::scale3(Ks[type] * r_mag * gamma, s, Fs);
// Calculate torque due to tangential displacements
MathExtra::cross3(r0, rb, t);
MathExtra::norm3(t);
MathExtra::scale3(0.5 * r_mag * Ks[type] * r_mag * gamma, t, Ts);
// Relative rotation force/torque
// Use representation of X'Y'Z' rotations from Wang, Mora 2009
temp = r_mag + rb[2];
if (temp < 0.0) temp = 0.0;
mq[0] = sqrt(2) * 0.5 * sqrt(temp * r_mag_inv);
temp = sqrt(rb[0] * rb[0] + rb[1] * rb[1]);
if (temp != 0.0) {
mq[1] = -sqrt(2) * 0.5 / temp;
temp = r_mag - rb[2];
if (temp < 0.0) temp = 0.0;
mq[1] *= sqrt(temp * r_mag_inv);
mq[2] = -mq[1];
mq[1] *= rb[1];
mq[2] *= rb[0];
} else {
// If aligned along z axis, x,y terms zero (r_mag-rb[2] = 0)
mq[1] = 0.0;
mq[2] = 0.0;
}
mq[3] = 0.0;
// qp21 = opposite of r^\circ_21 in Wang
// q21 = opposite of r_21 in Wang
MathExtra::quatquat(q2inv, q1, qp21);
MathExtra::qconjugate(mq, mqinv);
MathExtra::quatquat(mqinv, qp21, qtmp);
MathExtra::quatquat(qtmp, mq, q21);
temp = sqrt(q21[0] * q21[0] + q21[3] * q21[3]);
if (temp != 0.0) {
c = q21[0] / temp;
psi = 2.0 * acos_limit(c);
} else {
c = 0.0;
psi = 0.0;
}
// Map negative rotations
if (q21[3] < 0.0) // sin = q21[3]/temp
psi = -psi;
if (q21[3] == 0.0) psi = 0.0;
c = q21[0] * q21[0] - q21[1] * q21[1] - q21[2] * q21[2] + q21[3] * q21[3];
theta = acos_limit(c);
// Separately calculte magnitude of quaternion in x-y and out of x-y planes
// to avoid dividing by zero
mag_out_plane = (q21[0] * q21[0] + q21[3] * q21[3]);
mag_in_plane = (q21[1] * q21[1] + q21[2] * q21[2]);
if (mag_in_plane == 0.0) {
// No rotation => no bending/shear torque or extra shear force
// achieve by setting cos/sin = 0
cos_phi = 0.0;
sin_phi = 0.0;
} else if (mag_out_plane == 0.0) {
// Calculate angle in plane
cos_phi = q21[2] / sqrt(mag_in_plane);
sin_phi = -q21[1] / sqrt(mag_in_plane);
} else {
// Default equations in Mora, Wang 2009
cos_phi = q21[1] * q21[3] + q21[0] * q21[2];
sin_phi = q21[2] * q21[3] - q21[0] * q21[1];
cos_phi /= sqrt(mag_out_plane * mag_in_plane);
sin_phi /= sqrt(mag_out_plane * mag_in_plane);
}
Tbp[0] = -Kb[type] * theta * sin_phi;
Tbp[1] = Kb[type] * theta * cos_phi;
Tbp[2] = 0.0;
Ttp[0] = 0.0;
Ttp[1] = 0.0;
Ttp[2] = Kt[type] * psi;
Fsp[0] = -0.5 * Ks[type] * r_mag * theta * cos_phi;
Fsp[1] = -0.5 * Ks[type] * r_mag * theta * sin_phi;
Fsp[2] = 0.0;
Tsp[0] = 0.25 * Ks[type] * r_mag * r_mag * theta * sin_phi;
Tsp[1] = -0.25 * Ks[type] * r_mag * r_mag * theta * cos_phi;
Tsp[2] = 0.0;
// Rotate forces/torques back to 1st particle's frame
MathExtra::quatrotvec(mq, Fsp, Ftmp);
MathExtra::quatrotvec(mq, Tsp, Ttmp);
for (int m = 0; m < 3; m++) {
Fs[m] += Ftmp[m];
Ts[m] += Ttmp[m];
}
MathExtra::quatrotvec(mq, Tbp, Tb);
MathExtra::quatrotvec(mq, Ttp, Tt);
// Sum forces and calculate magnitudes
F_rot[0] += Fs[0];
F_rot[1] += Fs[1];
F_rot[2] += Fs[2];
MathExtra::quatrotvec(q2, F_rot, force1on2);
T_rot[0] = Ts[0] + Tt[0] + Tb[0];
T_rot[1] = Ts[1] + Tt[1] + Tb[1];
T_rot[2] = Ts[2] + Tt[2] + Tb[2];
MathExtra::quatrotvec(q2, T_rot, torque1on2);
T_rot[0] = Ts[0] - Tt[0] - Tb[0];
T_rot[1] = Ts[1] - Tt[1] - Tb[1];
T_rot[2] = Ts[2] - Tt[2] - Tb[2];
MathExtra::quatrotvec(q2, T_rot, torque2on1);
Fs_mag = MathExtra::len3(Fs);
Tt_mag = MathExtra::len3(Tt);
Tb_mag = MathExtra::len3(Tb);
breaking = Fr / Fcr[type] + Fs_mag / Fcs[type] + Tb_mag / Tcb[type] + Tt_mag / Tct[type];
if (breaking < 0.0) breaking = 0.0;
return breaking;
}
/* ----------------------------------------------------------------------
Calculate damping using formulation in
Y. Wang, F. Alonso-Marroquin, W. Guo 2015
Note: n points towards 1 vs pointing towards 2
---------------------------------------------------------------------- */
void BondBPMRotational::damping_forces(int i1, int i2, int type, double &Fr, double *rhat,
double *r, double *force1on2, double *torque1on2,
double *torque2on1)
{
double v1dotr, v2dotr, w1dotr, w2dotr;
double s1[3], s2[3], tdamp[3], tmp[3];
double vn1[3], vn2[3], vt1[3], vt2[3], vroll[3];
double wxn1[3], wxn2[3], wn1[3], wn2[3];
double **v = atom->v;
double **omega = atom->omega;
// Damp normal velocity difference
v1dotr = MathExtra::dot3(v[i1], rhat);
v2dotr = MathExtra::dot3(v[i2], rhat);
MathExtra::scale3(v1dotr, rhat, vn1);
MathExtra::scale3(v2dotr, rhat, vn2);
MathExtra::sub3(vn1, vn2, tmp);
MathExtra::scale3(gnorm[type], tmp);
Fr = MathExtra::lensq3(tmp);
MathExtra::add3(force1on2, tmp, force1on2);
// Damp tangential objective velocities
MathExtra::sub3(v[i1], vn1, vt1);
MathExtra::sub3(v[i2], vn2, vt2);
MathExtra::sub3(vt2, vt1, tmp);
MathExtra::scale3(0.5, tmp);
MathExtra::cross3(omega[i1], r, s1);
MathExtra::scale3(-0.5, s1);
MathExtra::sub3(s1, tmp, s1); // Eq 12
MathExtra::cross3(omega[i2], r, s2);
MathExtra::scale3(0.5, s2);
MathExtra::add3(s2, tmp, s2); // Eq 13
MathExtra::sub3(s1, s2, tmp);
MathExtra::scale3(gslide[type], tmp);
MathExtra::add3(force1on2, tmp, force1on2);
// Apply corresponding torque
MathExtra::cross3(r, tmp, tdamp);
MathExtra::scale3(0.5, tdamp);
MathExtra::add3(torque1on2, tdamp, torque1on2);
MathExtra::add3(torque2on1, tdamp, torque2on1);
// Damp rolling
MathExtra::cross3(omega[i1], rhat, wxn1);
MathExtra::cross3(omega[i2], rhat, wxn2);
MathExtra::sub3(wxn1, wxn2, vroll); // Eq. 31
MathExtra::cross3(r, vroll, tdamp);
MathExtra::scale3(0.5 * groll[type], tdamp);
MathExtra::add3(torque1on2, tdamp, torque1on2);
MathExtra::scale3(-1.0, tdamp);
MathExtra::add3(torque2on1, tdamp, torque2on1);
// Damp twist
w1dotr = MathExtra::dot3(omega[i1], rhat);
w2dotr = MathExtra::dot3(omega[i2], rhat);
MathExtra::scale3(w1dotr, rhat, wn1);
MathExtra::scale3(w2dotr, rhat, wn2);
MathExtra::sub3(wn1, wn2, tdamp); // Eq. 38
MathExtra::scale3(0.5 * gtwist[type], tdamp);
MathExtra::add3(torque1on2, tdamp, torque1on2);
MathExtra::scale3(-1.0, tdamp);
MathExtra::add3(torque2on1, tdamp, torque2on1);
}
/* ---------------------------------------------------------------------- */
void BondBPMRotational::compute(int eflag, int vflag)
{
if (!fix_bond_history->stored_flag) {
fix_bond_history->stored_flag = true;
store_data();
}
int i1, i2, itmp, n, type;
double r[3], r0[3], rhat[3];
double rsq, r0_mag, r_mag, r_mag_inv;
double Fr, breaking, smooth;
double force1on2[3], torque1on2[3], torque2on1[3];
ev_init(eflag, vflag);
double **x = atom->x;
double **f = atom->f;
double **torque = atom->torque;
tagint *tag = atom->tag;
int **bondlist = neighbor->bondlist;
int nbondlist = neighbor->nbondlist;
int nlocal = atom->nlocal;
int newton_bond = force->newton_bond;
double **bondstore = fix_bond_history->bondstore;
for (n = 0; n < nbondlist; n++) {
// skip bond if already broken
if (bondlist[n][2] <= 0) continue;
i1 = bondlist[n][0];
i2 = bondlist[n][1];
type = bondlist[n][2];
r0_mag = bondstore[n][0];
// Ensure pair is always ordered such that r0 points in
// a consistent direction and to ensure numerical operations
// are identical to minimize the possibility that a bond straddling
// an mpi grid (newton off) doesn't break on one proc but not the other
if (tag[i2] < tag[i1]) {
itmp = i1;
i1 = i2;
i2 = itmp;
}
// If bond hasn't been set - should be initialized to zero
if (r0_mag < EPSILON || std::isnan(r0_mag)) r0_mag = store_bond(n, i1, i2);
r0[0] = bondstore[n][1];
r0[1] = bondstore[n][2];
r0[2] = bondstore[n][3];
MathExtra::scale3(r0_mag, r0);
// Note this is the reverse of Mora & Wang
MathExtra::sub3(x[i1], x[i2], r);
rsq = MathExtra::lensq3(r);
r_mag = sqrt(rsq);
r_mag_inv = 1.0 / r_mag;
MathExtra::scale3(r_mag_inv, r, rhat);
// ------------------------------------------------------//
// Calculate forces, check if bond breaks
// ------------------------------------------------------//
breaking = elastic_forces(i1, i2, type, Fr, r_mag, r0_mag, r_mag_inv, rhat, r, r0, force1on2,
torque1on2, torque2on1);
if (breaking >= 1.0) {
bondlist[n][2] = 0;
process_broken(i1, i2);
continue;
}
damping_forces(i1, i2, type, Fr, rhat, r, force1on2, torque1on2, torque2on1);
if (smooth_flag) {
smooth = breaking * breaking;
smooth = 1.0 - smooth * smooth;
} else {
smooth = 1.0;
}
// ------------------------------------------------------//
// Apply forces and torques to particles
// ------------------------------------------------------//
if (newton_bond || i1 < nlocal) {
f[i1][0] -= force1on2[0] * smooth;
f[i1][1] -= force1on2[1] * smooth;
f[i1][2] -= force1on2[2] * smooth;
torque[i1][0] += torque2on1[0] * smooth;
torque[i1][1] += torque2on1[1] * smooth;
torque[i1][2] += torque2on1[2] * smooth;
}
if (newton_bond || i2 < nlocal) {
f[i2][0] += force1on2[0] * smooth;
f[i2][1] += force1on2[1] * smooth;
f[i2][2] += force1on2[2] * smooth;
torque[i2][0] += torque1on2[0] * smooth;
torque[i2][1] += torque1on2[1] * smooth;
torque[i2][2] += torque1on2[2] * smooth;
}
if (evflag) ev_tally(i1, i2, nlocal, newton_bond, 0.0, Fr * smooth, r[0], r[1], r[2]);
}
}
/* ---------------------------------------------------------------------- */
void BondBPMRotational::allocate()
{
allocated = 1;
const int np1 = atom->nbondtypes + 1;
memory->create(Kr, np1, "bond:Kr");
memory->create(Ks, np1, "bond:Ks");
memory->create(Kt, np1, "bond:Kt");
memory->create(Kb, np1, "bond:Kb");
memory->create(Fcr, np1, "bond:Fcr");
memory->create(Fcs, np1, "bond:Fcs");
memory->create(Tct, np1, "bond:Tct");
memory->create(Tcb, np1, "bond:Tcb");
memory->create(gnorm, np1, "bond:gnorm");
memory->create(gslide, np1, "bond:gslide");
memory->create(groll, np1, "bond:groll");
memory->create(gtwist, np1, "bond:gtwist");
memory->create(setflag, np1, "bond:setflag");
for (int i = 1; i < np1; i++) setflag[i] = 0;
}
/* ----------------------------------------------------------------------
set coeffs for one or more types
------------------------------------------------------------------------- */
void BondBPMRotational::coeff(int narg, char **arg)
{
if (narg != 13) error->all(FLERR, "Incorrect args for bond coefficients");
if (!allocated) allocate();
int ilo, ihi;
utils::bounds(FLERR, arg[0], 1, atom->nbondtypes, ilo, ihi, error);
double Kr_one = utils::numeric(FLERR, arg[1], false, lmp);
double Ks_one = utils::numeric(FLERR, arg[2], false, lmp);
double Kt_one = utils::numeric(FLERR, arg[3], false, lmp);
double Kb_one = utils::numeric(FLERR, arg[4], false, lmp);
double Fcr_one = utils::numeric(FLERR, arg[5], false, lmp);
double Fcs_one = utils::numeric(FLERR, arg[6], false, lmp);
double Tct_one = utils::numeric(FLERR, arg[7], false, lmp);
double Tcb_one = utils::numeric(FLERR, arg[8], false, lmp);
double gnorm_one = utils::numeric(FLERR, arg[9], false, lmp);
double gslide_one = utils::numeric(FLERR, arg[10], false, lmp);
double groll_one = utils::numeric(FLERR, arg[11], false, lmp);
double gtwist_one = utils::numeric(FLERR, arg[12], false, lmp);
int count = 0;
for (int i = ilo; i <= ihi; i++) {
Kr[i] = Kr_one;
Ks[i] = Ks_one;
Kt[i] = Kt_one;
Kb[i] = Kb_one;
Fcr[i] = Fcr_one;
Fcs[i] = Fcs_one;
Tct[i] = Tct_one;
Tcb[i] = Tcb_one;
gnorm[i] = gnorm_one;
gslide[i] = gslide_one;
groll[i] = groll_one;
gtwist[i] = gtwist_one;
setflag[i] = 1;
count++;
if (Fcr[i] / Kr[i] > max_stretch) max_stretch = Fcr[i] / Kr[i];
}
if (count == 0) error->all(FLERR, "Incorrect args for bond coefficients");
}
/* ----------------------------------------------------------------------
check for correct settings and create fix
------------------------------------------------------------------------- */
void BondBPMRotational::init_style()
{
BondBPM::init_style();
if (!atom->quat_flag || !atom->sphere_flag)
error->all(FLERR, "Bond bpm/rotational requires atom style bpm/sphere");
if (comm->ghost_velocity == 0)
error->all(FLERR, "Bond bpm/rotational requires ghost atoms store velocity");
if (domain->dimension == 2)
error->warning(FLERR, "Bond style bpm/rotational not intended for 2d use");
if (!id_fix_bond_history) {
id_fix_bond_history = utils::strdup("HISTORY_BPM_ROTATIONAL");
fix_bond_history = dynamic_cast<FixBondHistory *>(modify->replace_fix(
id_fix_dummy2, fmt::format("{} all BOND_HISTORY 0 4", id_fix_bond_history), 1));
delete[] id_fix_dummy2;
id_fix_dummy2 = nullptr;
}
}
/* ---------------------------------------------------------------------- */
void BondBPMRotational::settings(int narg, char **arg)
{
BondBPM::settings(narg, arg);
int iarg;
for (std::size_t i = 0; i < leftover_iarg.size(); i++) {
iarg = leftover_iarg[i];
if (strcmp(arg[iarg], "smooth") == 0) {
if (iarg + 1 > narg) error->all(FLERR, "Illegal bond bpm command");
smooth_flag = utils::logical(FLERR, arg[iarg + 1], false, lmp);
i += 1;
} else {
error->all(FLERR, "Illegal bond_style command");
}
}
}
/* ----------------------------------------------------------------------
proc 0 writes out coeffs to restart file
------------------------------------------------------------------------- */
void BondBPMRotational::write_restart(FILE *fp)
{
fwrite(&Kr[1], sizeof(double), atom->nbondtypes, fp);
fwrite(&Ks[1], sizeof(double), atom->nbondtypes, fp);
fwrite(&Kt[1], sizeof(double), atom->nbondtypes, fp);
fwrite(&Kb[1], sizeof(double), atom->nbondtypes, fp);
fwrite(&Fcr[1], sizeof(double), atom->nbondtypes, fp);
fwrite(&Fcs[1], sizeof(double), atom->nbondtypes, fp);
fwrite(&Tct[1], sizeof(double), atom->nbondtypes, fp);
fwrite(&Tcb[1], sizeof(double), atom->nbondtypes, fp);
fwrite(&gnorm[1], sizeof(double), atom->nbondtypes, fp);
fwrite(&gslide[1], sizeof(double), atom->nbondtypes, fp);
fwrite(&groll[1], sizeof(double), atom->nbondtypes, fp);
fwrite(&gtwist[1], sizeof(double), atom->nbondtypes, fp);
}
/* ----------------------------------------------------------------------
proc 0 reads coeffs from restart file, bcasts them
------------------------------------------------------------------------- */
void BondBPMRotational::read_restart(FILE *fp)
{
allocate();
if (comm->me == 0) {
utils::sfread(FLERR, &Kr[1], sizeof(double), atom->nbondtypes, fp, nullptr, error);
utils::sfread(FLERR, &Ks[1], sizeof(double), atom->nbondtypes, fp, nullptr, error);
utils::sfread(FLERR, &Kt[1], sizeof(double), atom->nbondtypes, fp, nullptr, error);
utils::sfread(FLERR, &Kb[1], sizeof(double), atom->nbondtypes, fp, nullptr, error);
utils::sfread(FLERR, &Fcr[1], sizeof(double), atom->nbondtypes, fp, nullptr, error);
utils::sfread(FLERR, &Fcs[1], sizeof(double), atom->nbondtypes, fp, nullptr, error);
utils::sfread(FLERR, &Tct[1], sizeof(double), atom->nbondtypes, fp, nullptr, error);
utils::sfread(FLERR, &Tcb[1], sizeof(double), atom->nbondtypes, fp, nullptr, error);
utils::sfread(FLERR, &gnorm[1], sizeof(double), atom->nbondtypes, fp, nullptr, error);
utils::sfread(FLERR, &gslide[1], sizeof(double), atom->nbondtypes, fp, nullptr, error);
utils::sfread(FLERR, &groll[1], sizeof(double), atom->nbondtypes, fp, nullptr, error);
utils::sfread(FLERR, &gtwist[1], sizeof(double), atom->nbondtypes, fp, nullptr, error);
}
MPI_Bcast(&Kr[1], atom->nbondtypes, MPI_DOUBLE, 0, world);
MPI_Bcast(&Ks[1], atom->nbondtypes, MPI_DOUBLE, 0, world);
MPI_Bcast(&Kt[1], atom->nbondtypes, MPI_DOUBLE, 0, world);
MPI_Bcast(&Kb[1], atom->nbondtypes, MPI_DOUBLE, 0, world);
MPI_Bcast(&Fcr[1], atom->nbondtypes, MPI_DOUBLE, 0, world);
MPI_Bcast(&Fcs[1], atom->nbondtypes, MPI_DOUBLE, 0, world);
MPI_Bcast(&Tct[1], atom->nbondtypes, MPI_DOUBLE, 0, world);
MPI_Bcast(&Tcb[1], atom->nbondtypes, MPI_DOUBLE, 0, world);
MPI_Bcast(&gnorm[1], atom->nbondtypes, MPI_DOUBLE, 0, world);
MPI_Bcast(&gslide[1], atom->nbondtypes, MPI_DOUBLE, 0, world);
MPI_Bcast(&groll[1], atom->nbondtypes, MPI_DOUBLE, 0, world);
MPI_Bcast(&gtwist[1], atom->nbondtypes, MPI_DOUBLE, 0, world);
for (int i = 1; i <= atom->nbondtypes; i++) setflag[i] = 1;
}
/* ----------------------------------------------------------------------
proc 0 writes to data file
------------------------------------------------------------------------- */
void BondBPMRotational::write_data(FILE *fp)
{
for (int i = 1; i <= atom->nbondtypes; i++)
fprintf(fp, "%d %g %g %g %g %g %g %g %g %g %g %g %g\n", i, Kr[i], Ks[i], Kt[i], Kb[i], Fcr[i],
Fcs[i], Tct[i], Tcb[i], gnorm[i], gslide[i], groll[i], gtwist[i]);
}
/* ---------------------------------------------------------------------- */
double BondBPMRotational::single(int type, double rsq, int i, int j, double &fforce)
{
// Not yet enabled
if (type <= 0) return 0.0;
int itmp;
if (atom->tag[j] < atom->tag[i]) {
itmp = i;
i = j;
j = itmp;
}
double r0_mag, r_mag, r_mag_inv;
double r0[3], r[3], rhat[3];
for (int n = 0; n < atom->num_bond[i]; n++) {
if (atom->bond_atom[i][n] == atom->tag[j]) {
r0_mag = fix_bond_history->get_atom_value(i, n, 0);
r0[0] = fix_bond_history->get_atom_value(i, n, 1);
r0[1] = fix_bond_history->get_atom_value(i, n, 2);
r0[2] = fix_bond_history->get_atom_value(i, n, 3);
}
}
double **x = atom->x;
MathExtra::scale3(r0_mag, r0);
MathExtra::sub3(x[i], x[j], r);
r_mag = sqrt(rsq);
r_mag_inv = 1.0 / r_mag;
MathExtra::scale3(r_mag_inv, r, rhat);
double breaking, smooth, Fr;
double force1on2[3], torque1on2[3], torque2on1[3];
breaking = elastic_forces(i, j, type, Fr, r_mag, r0_mag, r_mag_inv, rhat, r, r0, force1on2,
torque1on2, torque2on1);
fforce = Fr;
damping_forces(i, j, type, Fr, rhat, r, force1on2, torque1on2, torque2on1);
fforce += Fr;
if (smooth_flag) {
smooth = breaking * breaking;
smooth = 1.0 - smooth * smooth;
fforce *= smooth;
}
return 0.0;
}

59
src/BPM/bond_bpm_rotational.h Normal file
View File

@ -0,0 +1,59 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef BOND_CLASS
// clang-format off
BondStyle(bpm/rotational,BondBPMRotational);
// clang-format on
#else
#ifndef LMP_BOND_BPM_ROTATIONAL_H
#define LMP_BOND_BPM_ROTATIONAL_H
#include "bond_bpm.h"
namespace LAMMPS_NS {
class BondBPMRotational : public BondBPM {
public:
BondBPMRotational(class LAMMPS *);
~BondBPMRotational() override;
void compute(int, int) override;
void coeff(int, char **) override;
void init_style() override;
void settings(int, char **) override;
void write_restart(FILE *) override;
void read_restart(FILE *) override;
void write_data(FILE *) override;
double single(int, double, int, int, double &) override;
protected:
double *Kr, *Ks, *Kt, *Kb, *gnorm, *gslide, *groll, *gtwist;
double *Fcr, *Fcs, *Tct, *Tcb;
int smooth_flag;
double acos_limit(double);
double elastic_forces(int, int, int, double &, double, double, double, double *, double *,
double *, double *, double *, double *);
void damping_forces(int, int, int, double &, double *, double *, double *, double *, double *);
void allocate();
void store_data();
double store_bond(int, int, int);
};
} // namespace LAMMPS_NS
#endif
#endif

381
src/BPM/bond_bpm_spring.cpp Normal file
View File

@ -0,0 +1,381 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "bond_bpm_spring.h"
#include "atom.h"
#include "comm.h"
#include "domain.h"
#include "error.h"
#include "fix_bond_history.h"
#include "force.h"
#include "memory.h"
#include "modify.h"
#include "neighbor.h"
#define EPSILON 1e-10
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
BondBPMSpring::BondBPMSpring(LAMMPS *_lmp) :
BondBPM(_lmp), k(nullptr), ecrit(nullptr), gamma(nullptr)
{
partial_flag = 1;
smooth_flag = 1;
}
/* ---------------------------------------------------------------------- */
BondBPMSpring::~BondBPMSpring()
{
if (allocated) {
memory->destroy(setflag);
memory->destroy(k);
memory->destroy(ecrit);
memory->destroy(gamma);
}
}
/* ----------------------------------------------------------------------
Store data for a single bond - if bond added after LAMMPS init (e.g. pour)
------------------------------------------------------------------------- */
double BondBPMSpring::store_bond(int n, int i, int j)
{
double delx, dely, delz, r;
double **x = atom->x;
double **bondstore = fix_bond_history->bondstore;
tagint *tag = atom->tag;
delx = x[i][0] - x[j][0];
dely = x[i][1] - x[j][1];
delz = x[i][2] - x[j][2];
r = sqrt(delx * delx + dely * dely + delz * delz);
bondstore[n][0] = r;
if (i < atom->nlocal) {
for (int m = 0; m < atom->num_bond[i]; m++) {
if (atom->bond_atom[i][m] == tag[j]) { fix_bond_history->update_atom_value(i, m, 0, r); }
}
}
if (j < atom->nlocal) {
for (int m = 0; m < atom->num_bond[j]; m++) {
if (atom->bond_atom[j][m] == tag[i]) { fix_bond_history->update_atom_value(j, m, 0, r); }
}
}
return r;
}
/* ----------------------------------------------------------------------
Store data for all bonds called once
------------------------------------------------------------------------- */
void BondBPMSpring::store_data()
{
int i, j, m, type;
double delx, dely, delz, r;
double **x = atom->x;
int **bond_type = atom->bond_type;
for (i = 0; i < atom->nlocal; i++) {
for (m = 0; m < atom->num_bond[i]; m++) {
type = bond_type[i][m];
//Skip if bond was turned off
if (type < 0) continue;
// map to find index n
j = atom->map(atom->bond_atom[i][m]);
if (j == -1) error->one(FLERR, "Atom missing in BPM bond");
delx = x[i][0] - x[j][0];
dely = x[i][1] - x[j][1];
delz = x[i][2] - x[j][2];
// Get closest image in case bonded with ghost
domain->minimum_image(delx, dely, delz);
r = sqrt(delx * delx + dely * dely + delz * delz);
fix_bond_history->update_atom_value(i, m, 0, r);
}
}
fix_bond_history->post_neighbor();
}
/* ---------------------------------------------------------------------- */
void BondBPMSpring::compute(int eflag, int vflag)
{
if (!fix_bond_history->stored_flag) {
fix_bond_history->stored_flag = true;
store_data();
}
int i1, i2, itmp, n, type;
double delx, dely, delz, delvx, delvy, delvz;
double e, rsq, r, r0, rinv, smooth, fbond, dot;
ev_init(eflag, vflag);
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
tagint *tag = atom->tag;
int **bondlist = neighbor->bondlist;
int nbondlist = neighbor->nbondlist;
int nlocal = atom->nlocal;
int newton_bond = force->newton_bond;
double **bondstore = fix_bond_history->bondstore;
for (n = 0; n < nbondlist; n++) {
// skip bond if already broken
if (bondlist[n][2] <= 0) continue;
i1 = bondlist[n][0];
i2 = bondlist[n][1];
type = bondlist[n][2];
r0 = bondstore[n][0];
// Ensure pair is always ordered to ensure numerical operations
// are identical to minimize the possibility that a bond straddling
// an mpi grid (newton off) doesn't break on one proc but not the other
if (tag[i2] < tag[i1]) {
itmp = i1;
i1 = i2;
i2 = itmp;
}
// If bond hasn't been set - should be initialized to zero
if (r0 < EPSILON || std::isnan(r0)) r0 = store_bond(n, i1, i2);
delx = x[i1][0] - x[i2][0];
dely = x[i1][1] - x[i2][1];
delz = x[i1][2] - x[i2][2];
rsq = delx * delx + dely * dely + delz * delz;
r = sqrt(rsq);
e = (r - r0) / r0;
if (fabs(e) > ecrit[type]) {
bondlist[n][2] = 0;
process_broken(i1, i2);
continue;
}
rinv = 1.0 / r;
fbond = k[type] * (r0 - r);
delvx = v[i1][0] - v[i2][0];
delvy = v[i1][1] - v[i2][1];
delvz = v[i1][2] - v[i2][2];
dot = delx * delvx + dely * delvy + delz * delvz;
fbond -= gamma[type] * dot * rinv;
fbond *= rinv;
if (smooth_flag) {
smooth = (r - r0) / (r0 * ecrit[type]);
smooth *= smooth;
smooth *= smooth;
smooth *= smooth;
smooth = 1 - smooth;
fbond *= smooth;
}
if (newton_bond || i1 < nlocal) {
f[i1][0] += delx * fbond;
f[i1][1] += dely * fbond;
f[i1][2] += delz * fbond;
}
if (newton_bond || i2 < nlocal) {
f[i2][0] -= delx * fbond;
f[i2][1] -= dely * fbond;
f[i2][2] -= delz * fbond;
}
if (evflag) ev_tally(i1, i2, nlocal, newton_bond, 0.0, fbond, delx, dely, delz);
}
}
/* ---------------------------------------------------------------------- */
void BondBPMSpring::allocate()
{
allocated = 1;
const int np1 = atom->nbondtypes + 1;
memory->create(k, np1, "bond:k");
memory->create(ecrit, np1, "bond:ecrit");
memory->create(gamma, np1, "bond:gamma");
memory->create(setflag, np1, "bond:setflag");
for (int i = 1; i < np1; i++) setflag[i] = 0;
}
/* ----------------------------------------------------------------------
set coeffs for one or more types
------------------------------------------------------------------------- */
void BondBPMSpring::coeff(int narg, char **arg)
{
if (narg != 4) error->all(FLERR, "Incorrect args for bond coefficients");
if (!allocated) allocate();
int ilo, ihi;
utils::bounds(FLERR, arg[0], 1, atom->nbondtypes, ilo, ihi, error);
double k_one = utils::numeric(FLERR, arg[1], false, lmp);
double ecrit_one = utils::numeric(FLERR, arg[2], false, lmp);
double gamma_one = utils::numeric(FLERR, arg[3], false, lmp);
int count = 0;
for (int i = ilo; i <= ihi; i++) {
k[i] = k_one;
ecrit[i] = ecrit_one;
gamma[i] = gamma_one;
setflag[i] = 1;
count++;
if (1.0 + ecrit[i] > max_stretch) max_stretch = 1.0 + ecrit[i];
}
if (count == 0) error->all(FLERR, "Incorrect args for bond coefficients");
}
/* ----------------------------------------------------------------------
check for correct settings and create fix
------------------------------------------------------------------------- */
void BondBPMSpring::init_style()
{
BondBPM::init_style();
if (comm->ghost_velocity == 0)
error->all(FLERR, "Bond bpm/spring requires ghost atoms store velocity");
if (!id_fix_bond_history) {
id_fix_bond_history = utils::strdup("HISTORY_BPM_SPRING");
fix_bond_history = dynamic_cast<FixBondHistory *>(modify->replace_fix(
id_fix_dummy2, fmt::format("{} all BOND_HISTORY 0 1", id_fix_bond_history), 1));
delete[] id_fix_dummy2;
id_fix_dummy2 = nullptr;
}
}
/* ---------------------------------------------------------------------- */
void BondBPMSpring::settings(int narg, char **arg)
{
BondBPM::settings(narg, arg);
int iarg;
for (std::size_t i = 0; i < leftover_iarg.size(); i++) {
iarg = leftover_iarg[i];
if (strcmp(arg[iarg], "smooth") == 0) {
if (iarg + 1 > narg) error->all(FLERR, "Illegal bond bpm command");
smooth_flag = utils::logical(FLERR, arg[iarg + 1], false, lmp);
i += 1;
} else {
error->all(FLERR, "Illegal bond_style command");
}
}
}
/* ----------------------------------------------------------------------
proc 0 writes out coeffs to restart file
------------------------------------------------------------------------- */
void BondBPMSpring::write_restart(FILE *fp)
{
fwrite(&k[1], sizeof(double), atom->nbondtypes, fp);
fwrite(&ecrit[1], sizeof(double), atom->nbondtypes, fp);
fwrite(&gamma[1], sizeof(double), atom->nbondtypes, fp);
}
/* ----------------------------------------------------------------------
proc 0 reads coeffs from restart file, bcasts them
------------------------------------------------------------------------- */
void BondBPMSpring::read_restart(FILE *fp)
{
allocate();
if (comm->me == 0) {
utils::sfread(FLERR, &k[1], sizeof(double), atom->nbondtypes, fp, nullptr, error);
utils::sfread(FLERR, &ecrit[1], sizeof(double), atom->nbondtypes, fp, nullptr, error);
utils::sfread(FLERR, &gamma[1], sizeof(double), atom->nbondtypes, fp, nullptr, error);
}
MPI_Bcast(&k[1], atom->nbondtypes, MPI_DOUBLE, 0, world);
MPI_Bcast(&ecrit[1], atom->nbondtypes, MPI_DOUBLE, 0, world);
MPI_Bcast(&gamma[1], atom->nbondtypes, MPI_DOUBLE, 0, world);
for (int i = 1; i <= atom->nbondtypes; i++) setflag[i] = 1;
}
/* ----------------------------------------------------------------------
proc 0 writes to data file
------------------------------------------------------------------------- */
void BondBPMSpring::write_data(FILE *fp)
{
for (int i = 1; i <= atom->nbondtypes; i++)
fprintf(fp, "%d %g %g %g\n", i, k[i], ecrit[i], gamma[i]);
}
/* ---------------------------------------------------------------------- */
double BondBPMSpring::single(int type, double rsq, int i, int j, double &fforce)
{
if (type <= 0) return 0.0;
double r0;
for (int n = 0; n < atom->num_bond[i]; n++) {
if (atom->bond_atom[i][n] == atom->tag[j]) r0 = fix_bond_history->get_atom_value(i, n, 0);
}
double r = sqrt(rsq);
double rinv = 1.0 / r;
fforce = k[type] * (r0 - r);
double **x = atom->x;
double **v = atom->v;
double delx = x[i][0] - x[j][0];
double dely = x[i][1] - x[j][1];
double delz = x[i][2] - x[j][2];
double delvx = v[i][0] - v[j][0];
double delvy = v[i][1] - v[j][1];
double delvz = v[i][2] - v[j][2];
double dot = delx * delvx + dely * delvy + delz * delvz;
fforce -= gamma[type] * dot * rinv;
fforce *= rinv;
if (smooth_flag) {
double smooth = (r - r0) / (r0 * ecrit[type]);
smooth *= smooth;
smooth *= smooth;
smooth *= smooth;
smooth = 1 - smooth;
fforce *= smooth;
}
return 0.0;
}

52
src/BPM/bond_bpm_spring.h Normal file
View File

@ -0,0 +1,52 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef BOND_CLASS
// clang-format off
BondStyle(bpm/spring,BondBPMSpring);
// clang-format on
#else
#ifndef LMP_BOND_BPM_SPRING_H
#define LMP_BOND_BPM_SPRING_H
#include "bond_bpm.h"
namespace LAMMPS_NS {
class BondBPMSpring : public BondBPM {
public:
BondBPMSpring(class LAMMPS *);
~BondBPMSpring() override;
void compute(int, int) override;
void coeff(int, char **) override;
void init_style() override;
void settings(int, char **) override;
void write_restart(FILE *) override;
void read_restart(FILE *) override;
void write_data(FILE *) override;
double single(int, double, int, int, double &) override;
protected:
double *k, *ecrit, *gamma;
int smooth_flag;
void allocate();
void store_data();
double store_bond(int, int, int);
};
} // namespace LAMMPS_NS
#endif
#endif

145
src/BPM/compute_nbond_atom.cpp Normal file
View File

@ -0,0 +1,145 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "compute_nbond_atom.h"
#include "atom.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "memory.h"
#include "update.h"
#include <string.h>
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
ComputeNBondAtom::ComputeNBondAtom(LAMMPS *_lmp, int narg, char **arg) :
Compute(_lmp, narg, arg), nbond(nullptr)
{
if (narg < 3) error->all(FLERR, "Illegal compute nbond/atom command");
peratom_flag = 1;
size_peratom_cols = 0;
peatomflag = 1;
timeflag = 1;
comm_reverse = 1;
nmax = 0;
}
/* ---------------------------------------------------------------------- */
ComputeNBondAtom::~ComputeNBondAtom()
{
memory->destroy(nbond);
}
/* ---------------------------------------------------------------------- */
void ComputeNBondAtom::compute_peratom()
{
invoked_peratom = update->ntimestep;
if (update->eflag_atom != invoked_peratom)
error->all(FLERR, "Per-atom nbond was not tallied on needed timestep");
// grow local nbond array if necessary
// needs to be atom->nmax in length
if (atom->nmax > nmax) {
memory->destroy(nbond);
nmax = atom->nmax;
memory->create(nbond, nmax, "nbond/atom:nbond");
vector_atom = nbond;
}
// npair includes ghosts if either newton flag is set
// b/c some bonds/dihedrals call pair::ev_tally with pairwise info
// nbond includes ghosts if newton_bond is set
// ntotal includes ghosts if either newton flag is set
// KSpace includes ghosts if tip4pflag is set
int nlocal = atom->nlocal;
tagint **bond_atom = atom->bond_atom;
int **bond_type = atom->bond_type;
int ntotal = nlocal;
if (force->newton) ntotal += atom->nghost;
// set local nbond array
int i, j, k;
int *num_bond = atom->num_bond;
int newton_bond = force->newton_bond;
for (i = 0; i < ntotal; i++) nbond[i] = 0;
for (i = 0; i < nlocal; i++) {
for (j = 0; j < num_bond[i]; j++) {
if (bond_type[i][j] <= 0) continue;
k = atom->map(bond_atom[i][j]);
if (k < 0) continue;
nbond[i] += 1;
if (newton_bond) nbond[k] += 1;
}
}
// communicate ghost nbond between neighbor procs
if (force->newton) comm->reverse_comm(this);
// zero nbond of atoms not in group
// only do this after comm since ghost contributions must be included
int *mask = atom->mask;
for (i = 0; i < nlocal; i++)
if (!(mask[i] & groupbit)) nbond[i] = 0.0;
}
/* ---------------------------------------------------------------------- */
int ComputeNBondAtom::pack_reverse_comm(int n, int first, double *buf)
{
int i, m, last;
m = 0;
last = first + n;
for (i = first; i < last; i++) buf[m++] = nbond[i];
return m;
}
/* ---------------------------------------------------------------------- */
void ComputeNBondAtom::unpack_reverse_comm(int n, int *list, double *buf)
{
int i, j, m;
m = 0;
for (i = 0; i < n; i++) {
j = list[i];
nbond[j] += buf[m++];
}
}
/* ----------------------------------------------------------------------
memory usage of local atom-based array
------------------------------------------------------------------------- */
double ComputeNBondAtom::memory_usage()
{
double bytes = nmax * sizeof(double);
return bytes;
}

45
src/BPM/compute_nbond_atom.h Normal file
View File

@ -0,0 +1,45 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef COMPUTE_CLASS
// clang-format off
ComputeStyle(nbond/atom,ComputeNBondAtom);
// clang-format on
#else
#ifndef LMP_COMPUTE_NBOND_ATOM_H
#define LMP_COMPUTE_NBOND_ATOM_H
#include "compute.h"
namespace LAMMPS_NS {
class ComputeNBondAtom : public Compute {
public:
ComputeNBondAtom(class LAMMPS *, int, char **);
~ComputeNBondAtom() override;
void init() override {}
void compute_peratom() override;
int pack_reverse_comm(int, int, double *) override;
void unpack_reverse_comm(int, int *, double *) override;
double memory_usage() override;
private:
int nmax;
double *nbond;
};
} // namespace LAMMPS_NS
#endif
#endif

157
src/BPM/fix_nve_bpm_sphere.cpp Normal file
View File

@ -0,0 +1,157 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "fix_nve_bpm_sphere.h"
#include "atom.h"
#include "atom_vec.h"
#include "domain.h"
#include "error.h"
#include "force.h"
#include "math_extra.h"
#include <cmath>
#include <cstring>
using namespace LAMMPS_NS;
using namespace FixConst;
using namespace MathExtra;
/* ---------------------------------------------------------------------- */
FixNVEBPMSphere::FixNVEBPMSphere(LAMMPS *_lmp, int narg, char **arg) : FixNVE(_lmp, narg, arg)
{
if (narg < 3) error->all(FLERR, "Illegal fix nve/bpm/sphere command");
time_integrate = 1;
// process extra keywords
// inertia = moment of inertia prefactor for sphere or disc
inertia = 0.4;
int iarg = 3;
while (iarg < narg) {
if (strcmp(arg[iarg], "disc") == 0) {
inertia = 0.5;
if (domain->dimension != 2)
error->all(FLERR, "Fix nve/bpm/sphere disc requires 2d simulation");
iarg++;
} else
error->all(FLERR, "Illegal fix nve/bpm/sphere command");
}
inv_inertia = 1.0 / inertia;
// error checks
if (!atom->quat_flag || !atom->sphere_flag)
error->all(FLERR, "Fix nve/bpm/sphere requires atom style bpm/sphere");
}
/* ---------------------------------------------------------------------- */
void FixNVEBPMSphere::init()
{
FixNVE::init();
// check that all particles are finite-size spheres
// no point particles allowed
double *radius = atom->radius;
int *mask = atom->mask;
int nlocal = atom->nlocal;
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit)
if (radius[i] == 0.0) error->one(FLERR, "Fix nve/bpm/sphere requires extended particles");
}
/* ---------------------------------------------------------------------- */
void FixNVEBPMSphere::initial_integrate(int /*vflag*/)
{
double dtq, dtfm, dtirotate, particle_inertia;
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
double **omega = atom->omega;
double **torque = atom->torque;
double **quat = atom->quat;
double *radius = atom->radius;
double *rmass = atom->rmass;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
// set timestep here since dt may have changed or come via rRESPA
dtq = 0.5 * dtv;
// update v,x,omega,quat for all particles
// d_omega/dt = torque / inertia
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
dtfm = dtf / rmass[i];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
x[i][0] += dtv * v[i][0];
x[i][1] += dtv * v[i][1];
x[i][2] += dtv * v[i][2];
particle_inertia = inertia * (radius[i] * radius[i] * rmass[i]);
dtirotate = dtf / particle_inertia;
omega[i][0] += dtirotate * torque[i][0];
omega[i][1] += dtirotate * torque[i][1];
omega[i][2] += dtirotate * torque[i][2];
MathExtra::richardson_sphere(quat[i], omega[i], dtq);
}
}
}
/* ---------------------------------------------------------------------- */
void FixNVEBPMSphere::final_integrate()
{
double dtfm, dtirotate, particle_inertia;
double **v = atom->v;
double **f = atom->f;
double **omega = atom->omega;
double **torque = atom->torque;
double *rmass = atom->rmass;
double *radius = atom->radius;
int *mask = atom->mask;
int nlocal = atom->nlocal;
if (igroup == atom->firstgroup) nlocal = atom->nfirst;
// update v,omega for all particles
// d_omega/dt = torque / inertia
for (int i = 0; i < nlocal; i++)
if (mask[i] & groupbit) {
dtfm = dtf / rmass[i];
v[i][0] += dtfm * f[i][0];
v[i][1] += dtfm * f[i][1];
v[i][2] += dtfm * f[i][2];
particle_inertia = inertia * (radius[i] * radius[i] * rmass[i]);
dtirotate = dtf / particle_inertia;
omega[i][0] += dtirotate * torque[i][0];
omega[i][1] += dtirotate * torque[i][1];
omega[i][2] += dtirotate * torque[i][2];
}
}

44
src/BPM/fix_nve_bpm_sphere.h Normal file
View File

@ -0,0 +1,44 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef FIX_CLASS
// clang-format off
FixStyle(nve/bpm/sphere,FixNVEBPMSphere);
// clang-format on
#else
#ifndef LMP_FIX_NVE_BPM_SPHERE_H
#define LMP_FIX_NVE_BPM_SPHERE_H
#include "fix_nve.h"
namespace LAMMPS_NS {
class FixNVEBPMSphere : public FixNVE {
public:
FixNVEBPMSphere(class LAMMPS *, int, char **);
void init() override;
void initial_integrate(int) override;
void final_integrate() override;
protected:
double inertia, inv_inertia;
int extra;
int dlm;
};
} // namespace LAMMPS_NS
#endif
#endif

326
src/BPM/pair_bpm_spring.cpp Normal file
View File

@ -0,0 +1,326 @@
/* ----------------------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#include "pair_bpm_spring.h"
#include "atom.h"
#include "comm.h"
#include "error.h"
#include "force.h"
#include "memory.h"
#include "neigh_list.h"
#include "neighbor.h"
using namespace LAMMPS_NS;
/* ---------------------------------------------------------------------- */
PairBPMSpring::PairBPMSpring(LAMMPS *_lmp) : Pair(_lmp)
{
writedata = 1;
}
/* ---------------------------------------------------------------------- */
PairBPMSpring::~PairBPMSpring()
{
if (allocated) {
memory->destroy(setflag);
memory->destroy(cutsq);
memory->destroy(k);
memory->destroy(cut);
memory->destroy(gamma);
}
}
/* ---------------------------------------------------------------------- */
void PairBPMSpring::compute(int eflag, int vflag)
{
int i, j, ii, jj, inum, jnum, itype, jtype;
double xtmp, ytmp, ztmp, delx, dely, delz, evdwl, fpair;
double r, rsq, rinv, factor_lj;
int *ilist, *jlist, *numneigh, **firstneigh;
double vxtmp, vytmp, vztmp, delvx, delvy, delvz, dot, smooth;
evdwl = 0.0;
if (eflag || vflag)
ev_setup(eflag, vflag);
else
evflag = vflag_fdotr = 0;
double **x = atom->x;
double **v = atom->v;
double **f = atom->f;
int *type = atom->type;
int nlocal = atom->nlocal;
int newton_pair = force->newton_pair;
double *special_lj = force->special_lj;
inum = list->inum;
ilist = list->ilist;
numneigh = list->numneigh;
firstneigh = list->firstneigh;
// loop over neighbors of my atoms
for (ii = 0; ii < inum; ii++) {
i = ilist[ii];
xtmp = x[i][0];
ytmp = x[i][1];
ztmp = x[i][2];
vxtmp = v[i][0];
vytmp = v[i][1];
vztmp = v[i][2];
itype = type[i];
jlist = firstneigh[i];
jnum = numneigh[i];
for (jj = 0; jj < jnum; jj++) {
j = jlist[jj];
factor_lj = special_lj[sbmask(j)];
if (factor_lj == 0) continue;
j &= NEIGHMASK;
delx = xtmp - x[j][0];
dely = ytmp - x[j][1];
delz = ztmp - x[j][2];
rsq = delx * delx + dely * dely + delz * delz;
jtype = type[j];
if (rsq < cutsq[itype][jtype]) {
r = sqrt(rsq);
rinv = 1.0 / r;
fpair = k[itype][jtype] * (cut[itype][jtype] - r);
smooth = rsq / cutsq[itype][jtype];
smooth *= smooth;
smooth *= smooth;
smooth = 1.0 - smooth;
delvx = vxtmp - v[j][0];
delvy = vytmp - v[j][1];
delvz = vztmp - v[j][2];
dot = delx * delvx + dely * delvy + delz * delvz;
fpair -= gamma[itype][jtype] * dot * smooth * rinv;
fpair *= factor_lj * rinv;
if (eflag) evdwl = 0.0;
f[i][0] += delx * fpair;
f[i][1] += dely * fpair;
f[i][2] += delz * fpair;
if (newton_pair || j < nlocal) {
f[j][0] -= delx * fpair;
f[j][1] -= dely * fpair;
f[j][2] -= delz * fpair;
}
if (evflag) ev_tally(i, j, nlocal, newton_pair, evdwl, 0.0, fpair, delx, dely, delz);
}
}
}
if (vflag_fdotr) virial_fdotr_compute();
}
/* ----------------------------------------------------------------------
allocate all arrays
------------------------------------------------------------------------- */
void PairBPMSpring::allocate()
{
allocated = 1;
const int np1 = atom->ntypes + 1;
memory->create(setflag, np1, np1, "pair:setflag");
for (int i = 1; i < np1; i++)
for (int j = i; j < np1; j++) setflag[i][j] = 0;
memory->create(cutsq, np1, np1, "pair:cutsq");
memory->create(k, np1, np1, "pair:k");
memory->create(cut, np1, np1, "pair:cut");
memory->create(gamma, np1, np1, "pair:gamma");
}
/* ----------------------------------------------------------------------
global settings
------------------------------------------------------------------------- */
void PairBPMSpring::settings(int narg, char ** /*arg*/)
{
if (narg != 0) error->all(FLERR, "Illegal pair_style command");
}
/* ----------------------------------------------------------------------
set coeffs for one or more type pairs
------------------------------------------------------------------------- */
void PairBPMSpring::coeff(int narg, char **arg)
{
if (narg != 5) error->all(FLERR, "Incorrect args for pair coefficients");
if (!allocated) allocate();
int ilo, ihi, jlo, jhi;
utils::bounds(FLERR, arg[0], 1, atom->ntypes, ilo, ihi, error);
utils::bounds(FLERR, arg[1], 1, atom->ntypes, jlo, jhi, error);
double k_one = utils::numeric(FLERR, arg[2], false, lmp);
double cut_one = utils::numeric(FLERR, arg[3], false, lmp);
double gamma_one = utils::numeric(FLERR, arg[4], false, lmp);
if (cut_one <= 0.0) error->all(FLERR, "Incorrect args for pair coefficients");
int count = 0;
for (int i = ilo; i <= ihi; i++) {
for (int j = MAX(jlo, i); j <= jhi; j++) {
k[i][j] = k_one;
cut[i][j] = cut_one;
gamma[i][j] = gamma_one;
setflag[i][j] = 1;
count++;
}
}
if (count == 0) error->all(FLERR, "Incorrect args for pair coefficients");
}
/* ----------------------------------------------------------------------
init for one type pair i,j and corresponding j,i
------------------------------------------------------------------------- */
double PairBPMSpring::init_one(int i, int j)
{
if (setflag[i][j] == 0) {
cut[i][j] = mix_distance(cut[i][i], cut[j][j]);
k[i][j] = mix_energy(k[i][i], k[j][j], cut[i][i], cut[j][j]);
gamma[i][j] = mix_energy(gamma[i][i], gamma[j][j], cut[i][i], cut[j][j]);
}
cut[j][i] = cut[i][j];
k[j][i] = k[i][j];
gamma[j][i] = gamma[i][j];
return cut[i][j];
}
/* ----------------------------------------------------------------------
proc 0 writes to restart file
------------------------------------------------------------------------- */
void PairBPMSpring::write_restart(FILE *fp)
{
write_restart_settings(fp);
int i, j;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
fwrite(&setflag[i][j], sizeof(int), 1, fp);
if (setflag[i][j]) {
fwrite(&k[i][j], sizeof(double), 1, fp);
fwrite(&cut[i][j], sizeof(double), 1, fp);
fwrite(&gamma[i][j], sizeof(double), 1, fp);
}
}
}
/* ----------------------------------------------------------------------
proc 0 reads from restart file, bcasts
------------------------------------------------------------------------- */
void PairBPMSpring::read_restart(FILE *fp)
{
read_restart_settings(fp);
allocate();
int i, j;
int me = comm->me;
for (i = 1; i <= atom->ntypes; i++)
for (j = i; j <= atom->ntypes; j++) {
if (me == 0) utils::sfread(FLERR, &setflag[i][j], sizeof(int), 1, fp, nullptr, error);
MPI_Bcast(&setflag[i][j], 1, MPI_INT, 0, world);
if (setflag[i][j]) {
if (me == 0) {
utils::sfread(FLERR, &k[i][j], sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR, &cut[i][j], sizeof(double), 1, fp, nullptr, error);
utils::sfread(FLERR, &gamma[i][j], sizeof(double), 1, fp, nullptr, error);
}
MPI_Bcast(&k[i][j], 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&cut[i][j], 1, MPI_DOUBLE, 0, world);
MPI_Bcast(&gamma[i][j], 1, MPI_DOUBLE, 0, world);
}
}
}
/* ----------------------------------------------------------------------
proc 0 writes to data file
------------------------------------------------------------------------- */
void PairBPMSpring::write_data(FILE *fp)
{
for (int i = 1; i <= atom->ntypes; i++)
fprintf(fp, "%d %g %g %g\n", i, k[i][i], cut[i][i], gamma[i][i]);
}
/* ----------------------------------------------------------------------
proc 0 writes all pairs to data file
------------------------------------------------------------------------- */
void PairBPMSpring::write_data_all(FILE *fp)
{
for (int i = 1; i <= atom->ntypes; i++)
for (int j = i; j <= atom->ntypes; j++)
fprintf(fp, "%d %d %g %g %g\n", i, j, k[i][j], cut[i][j], gamma[i][j]);
}
/* ---------------------------------------------------------------------- */
double PairBPMSpring::single(int i, int j, int itype, int jtype, double rsq, double /*factor_coul*/,
double factor_lj, double &fforce)
{
double fpair, r, rinv;
double delx, dely, delz, delvx, delvy, delvz, dot, smooth;
if (rsq > cutsq[itype][jtype]) return 0.0;
double **x = atom->x;
double **v = atom->v;
r = sqrt(rsq);
rinv = 1.0 / r;
fpair = k[itype][jtype] * (cut[itype][jtype] - r);
smooth = rsq / cutsq[itype][jtype];
smooth *= smooth;
smooth = 1.0 - smooth;
delx = x[i][0] - x[j][0];
dely = x[i][1] - x[j][1];
delz = x[i][2] - x[j][2];
delvx = v[i][0] - v[j][0];
delvy = v[i][1] - v[j][1];
delvz = v[i][2] - v[j][2];
dot = delx * delvx + dely * delvy + delz * delvz;
fpair -= gamma[itype][jtype] * dot * rinv * smooth;
fpair *= factor_lj;
fforce = fpair;
return 0.0;
}

50
src/BPM/pair_bpm_spring.h Normal file
View File

@ -0,0 +1,50 @@
/* -*- c++ -*- ----------------------------------------------------------
LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator
https://www.lammps.org/, Sandia National Laboratories
Steve Plimpton, sjplimp@sandia.gov
Copyright (2003) Sandia Corporation. Under the terms of Contract
DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
certain rights in this software. This software is distributed under
the GNU General Public License.
See the README file in the top-level LAMMPS directory.
------------------------------------------------------------------------- */
#ifdef PAIR_CLASS
// clang-format off
PairStyle(bpm/spring,PairBPMSpring);
// clang-format on
#else
#ifndef LMP_PAIR_BPM_SPRING_H
#define LMP_PAIR_BPM_SPRING_H
#include "pair.h"
namespace LAMMPS_NS {
class PairBPMSpring : public Pair {
public:
PairBPMSpring(class LAMMPS *);
~PairBPMSpring() override;
void compute(int, int) override;
void settings(int, char **) override;
void coeff(int, char **) override;
double init_one(int, int) override;
void write_restart(FILE *) override;
void read_restart(FILE *) override;
void write_data(FILE *) override;
void write_data_all(FILE *) override;
double single(int, int, int, int, double, double, double, double &) override;
protected:
double **k, **cut, **gamma;
void allocate();
};
} // namespace LAMMPS_NS
#endif
#endif

5
src/BROWNIAN/fix_brownian.cpp
View File

@ -33,7 +33,8 @@ using namespace FixConst;
FixBrownian::FixBrownian(LAMMPS *lmp, int narg, char **arg) : FixBrownianBase(lmp, narg, arg)
{
if (dipole_flag || gamma_t_eigen_flag || gamma_r_eigen_flag || gamma_r_flag) {
if (dipole_flag || gamma_t_eigen_flag || gamma_r_eigen_flag || gamma_r_flag || rot_temp_flag ||
planar_rot_flag) {
error->all(FLERR, "Illegal fix brownian command.");
}
if (!gamma_t_flag) { error->all(FLERR, "Illegal fix brownian command."); }
@ -45,7 +46,7 @@ void FixBrownian::init()
{
FixBrownianBase::init();
g1 /= gamma_t;
g2 /= sqrt(gamma_t);
g2 *= sqrt(temp / gamma_t);
}
/* ---------------------------------------------------------------------- */

26
src/BROWNIAN/fix_brownian.h
View File

@ -27,9 +27,9 @@ namespace LAMMPS_NS {
class FixBrownian : public FixBrownianBase {
public:
FixBrownian(class LAMMPS *, int, char **);
virtual ~FixBrownian(){};
void init();
void initial_integrate(int);
void init() override;
void initial_integrate(int) override;
private:
template <int Tp_UNIFORM, int Tp_GAUSS, int Tp_2D> void initial_integrate_templated();
@ -38,23 +38,3 @@ class FixBrownian : public FixBrownianBase {
} // namespace LAMMPS_NS
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal fix brownian command.
Wrong number/type of input arguments.
E: Fix brownian viscous drag coefficient must be > 0.
Self-explanatory.
E: Fix brownian diffusion coefficient must be > 0.
Self-explanatory.
E: Fix brownian seed must be > 0.
Self-explanatory.
*/

82
src/BROWNIAN/fix_brownian_asphere.cpp
View File

@ -21,6 +21,7 @@
#include "atom.h"
#include "atom_vec_ellipsoid.h"
#include "comm.h"
#include "domain.h"
#include "error.h"
#include "math_extra.h"
@ -42,13 +43,20 @@ FixBrownianAsphere::FixBrownianAsphere(LAMMPS *lmp, int narg, char **arg) :
if (dipole_flag && !atom->mu_flag)
error->all(FLERR, "Fix brownian/asphere dipole requires atom attribute mu");
if (!atom->ellipsoid_flag)
error->all(FLERR, "Fix brownian/asphere requires atom style ellipsoid");
if (planar_rot_flag && (comm->me == 0)) {
error->warning(FLERR, "Ignoring first two entries of gamma_r_eigen since rotation is planar.");
}
}
/* ---------------------------------------------------------------------- */
void FixBrownianAsphere::init()
{
avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
avec = dynamic_cast<AtomVecEllipsoid *>(atom->style_match("ellipsoid"));
if (!avec) error->all(FLERR, "Compute brownian/asphere requires atom style ellipsoid");
// check that all particles are finite-size ellipsoids
@ -85,6 +93,9 @@ void FixBrownianAsphere::init()
}
FixBrownianBase::init();
g4 = g2 * sqrt(rot_temp);
g2 *= sqrt(temp);
}
/* ---------------------------------------------------------------------- */
@ -94,37 +105,55 @@ void FixBrownianAsphere::initial_integrate(int /*vflag */)
if (domain->dimension == 2) {
if (dipole_flag) {
if (!noise_flag) {
initial_integrate_templated<0, 0, 1, 1>();
initial_integrate_templated<0, 0, 1, 1, 0>();
} else if (gaussian_noise_flag) {
initial_integrate_templated<0, 1, 1, 1>();
initial_integrate_templated<0, 1, 1, 1, 0>();
} else {
initial_integrate_templated<1, 0, 1, 1>();
initial_integrate_templated<1, 0, 1, 1, 0>();
}
} else {
if (!noise_flag) {
initial_integrate_templated<0, 0, 0, 1>();
initial_integrate_templated<0, 0, 0, 1, 0>();
} else if (gaussian_noise_flag) {
initial_integrate_templated<0, 1, 0, 1>();
initial_integrate_templated<0, 1, 0, 1, 0>();
} else {
initial_integrate_templated<1, 0, 0, 1>();
initial_integrate_templated<1, 0, 0, 1, 0>();
}
}
} else if (planar_rot_flag) {
if (dipole_flag) {
if (!noise_flag) {
initial_integrate_templated<0, 0, 1, 0, 1>();
} else if (gaussian_noise_flag) {
initial_integrate_templated<0, 1, 1, 0, 1>();
} else {
initial_integrate_templated<1, 0, 1, 0, 1>();
}
} else {
if (!noise_flag) {
initial_integrate_templated<0, 0, 0, 0, 1>();
} else if (gaussian_noise_flag) {
initial_integrate_templated<0, 1, 0, 0, 1>();
} else {
initial_integrate_templated<1, 0, 0, 0, 1>();
}
}
} else {
if (dipole_flag) {
if (!noise_flag) {
initial_integrate_templated<0, 0, 1, 0>();
initial_integrate_templated<0, 0, 1, 0, 0>();
} else if (gaussian_noise_flag) {
initial_integrate_templated<0, 1, 1, 0>();
initial_integrate_templated<0, 1, 1, 0, 0>();
} else {
initial_integrate_templated<1, 0, 1, 0>();
initial_integrate_templated<1, 0, 1, 0, 0>();
}
} else {
if (!noise_flag) {
initial_integrate_templated<0, 0, 0, 0>();
initial_integrate_templated<0, 0, 0, 0, 0>();
} else if (gaussian_noise_flag) {
initial_integrate_templated<0, 1, 0, 0>();
initial_integrate_templated<0, 1, 0, 0, 0>();
} else {
initial_integrate_templated<1, 0, 0, 0>();
initial_integrate_templated<1, 0, 0, 0, 0>();
}
}
}
@ -133,7 +162,7 @@ void FixBrownianAsphere::initial_integrate(int /*vflag */)
/* ---------------------------------------------------------------------- */
template <int Tp_UNIFORM, int Tp_GAUSS, int Tp_DIPOLE, int Tp_2D>
template <int Tp_UNIFORM, int Tp_GAUSS, int Tp_DIPOLE, int Tp_2D, int Tp_2Drot>
void FixBrownianAsphere::initial_integrate_templated()
{
double **x = atom->x;
@ -172,21 +201,30 @@ void FixBrownianAsphere::initial_integrate_templated()
if (Tp_2D) {
tmp[0] = tmp[1] = 0.0;
if (Tp_UNIFORM) {
tmp[2] = g1 * tmp[2] * gamma_r_inv[2] + gamma_r_invsqrt[2] * (rng->uniform() - 0.5) * g2;
tmp[2] = g1 * tmp[2] * gamma_r_inv[2] + gamma_r_invsqrt[2] * (rng->uniform() - 0.5) * g4;
} else if (Tp_GAUSS) {
tmp[2] = g1 * tmp[2] * gamma_r_inv[2] + gamma_r_invsqrt[2] * rng->gaussian() * g2;
tmp[2] = g1 * tmp[2] * gamma_r_inv[2] + gamma_r_invsqrt[2] * rng->gaussian() * g4;
} else {
tmp[2] = g1 * tmp[2] * gamma_r_inv[2];
}
} else if (Tp_2Drot) {
tmp[0] = tmp[1] = 0.0;
if (Tp_UNIFORM) {
tmp[2] = g1 * tmp[2] * gamma_r_inv[2] + gamma_r_invsqrt[2] * (rng->uniform() - 0.5) * g4;
} else if (Tp_GAUSS) {
tmp[2] = g1 * tmp[2] * gamma_r_inv[2] + gamma_r_invsqrt[2] * rng->gaussian() * g4;
} else {
tmp[2] = g1 * tmp[2] * gamma_r_inv[2];
}
} else {
if (Tp_UNIFORM) {
tmp[0] = g1 * tmp[0] * gamma_r_inv[0] + gamma_r_invsqrt[0] * (rng->uniform() - 0.5) * g2;
tmp[1] = g1 * tmp[1] * gamma_r_inv[1] + gamma_r_invsqrt[1] * (rng->uniform() - 0.5) * g2;
tmp[2] = g1 * tmp[2] * gamma_r_inv[2] + gamma_r_invsqrt[2] * (rng->uniform() - 0.5) * g2;
tmp[0] = g1 * tmp[0] * gamma_r_inv[0] + gamma_r_invsqrt[0] * (rng->uniform() - 0.5) * g4;
tmp[1] = g1 * tmp[1] * gamma_r_inv[1] + gamma_r_invsqrt[1] * (rng->uniform() - 0.5) * g4;
tmp[2] = g1 * tmp[2] * gamma_r_inv[2] + gamma_r_invsqrt[2] * (rng->uniform() - 0.5) * g4;
} else if (Tp_GAUSS) {
tmp[0] = g1 * tmp[0] * gamma_r_inv[0] + gamma_r_invsqrt[0] * rng->gaussian() * g2;
tmp[1] = g1 * tmp[1] * gamma_r_inv[1] + gamma_r_invsqrt[1] * rng->gaussian() * g2;
tmp[2] = g1 * tmp[2] * gamma_r_inv[2] + gamma_r_invsqrt[2] * rng->gaussian() * g2;
tmp[0] = g1 * tmp[0] * gamma_r_inv[0] + gamma_r_invsqrt[0] * rng->gaussian() * g4;
tmp[1] = g1 * tmp[1] * gamma_r_inv[1] + gamma_r_invsqrt[1] * rng->gaussian() * g4;
tmp[2] = g1 * tmp[2] * gamma_r_inv[2] + gamma_r_invsqrt[2] * rng->gaussian() * g4;
} else {
tmp[0] = g1 * tmp[0] * gamma_r_inv[0];
tmp[1] = g1 * tmp[1] * gamma_r_inv[1];

47
src/BROWNIAN/fix_brownian_asphere.h
View File

@ -27,55 +27,18 @@ namespace LAMMPS_NS {
class FixBrownianAsphere : public FixBrownianBase {
public:
FixBrownianAsphere(class LAMMPS *, int, char **);
virtual ~FixBrownianAsphere(){};
void initial_integrate(int);
void init();
void initial_integrate(int) override;
void init() override;
protected:
class AtomVecEllipsoid *avec;
private:
template <int Tp_UNIFORM, int Tp_GAUSS, int Tp_DIPOLE, int Tp_2D>
template <int Tp_UNIFORM, int Tp_GAUSS, int Tp_DIPOLE, int Tp_2D, int Tp_2Drot>
void initial_integrate_templated();
double g4;
};
} // namespace LAMMPS_NS
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal fix brownian/asphere command.
Wrong number/type of input arguments.
E: Compute brownian/asphere requires atom style sphere
Self-explanatory.
E: Compute brownian/asphere requires atom style ellipsoid
Self-explanatory.
E: Compute brownian/asphere dipole requires atom attribute mu
Self-explanatory.
E: Fix brownian/asphere translational viscous drag coefficient must be > 0.
Self-explanatory.
E: Fix brownian/asphere rotational viscous drag coefficient must be > 0.
Self-explanatory.
E: Fix brownian/asphere translational diffusion coefficient must be > 0.
Self-explanatory.
E: Fix brownian/asphere rotational diffusion coefficient must be > 0.
Self-explanatory.
E: Fix brownian/asphere seed must be > 0.
*/

23
src/BROWNIAN/fix_brownian_base.cpp
View File

@ -43,6 +43,8 @@ FixBrownianBase::FixBrownianBase(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, n
gamma_t_flag = gamma_r_flag = 0;
gamma_t_eigen_flag = gamma_r_eigen_flag = 0;
dipole_flag = 0;
rot_temp_flag = 0;
planar_rot_flag = 0;
g2 = 0.0;
if (narg < 5) error->all(FLERR, "Illegal fix brownian command.");
@ -154,10 +156,26 @@ FixBrownianBase::FixBrownianBase(LAMMPS *lmp, int narg, char **arg) : Fix(lmp, n
if (gamma_r <= 0) error->all(FLERR, "Fix brownian gamma_r must be > 0.");
iarg = iarg + 2;
} else if (strcmp(arg[iarg], "rotation_temp") == 0) {
if (narg == iarg + 1) { error->all(FLERR, "Illegal fix brownian command."); }
rot_temp_flag = 1;
rot_temp = utils::numeric(FLERR, arg[iarg + 1], false, lmp);
if (rot_temp <= 0) error->all(FLERR, "Fix brownian rotation_temp must be > 0.");
iarg = iarg + 2;
} else if (strcmp(arg[iarg], "planar_rotation") == 0) {
planar_rot_flag = 1;
if (domain->dimension == 2)
error->all(FLERR, "The planar_rotation keyword is not allowed for 2D simulations");
iarg = iarg + 1;
} else {
error->all(FLERR, "Illegal fix brownian command.");
}
}
if (!rot_temp_flag) rot_temp = temp;
// initialize Marsaglia RNG with processor-unique seed
rng = new RanMars(lmp, seed + comm->me);
@ -196,14 +214,13 @@ void FixBrownianBase::init()
{
dt = update->dt;
sqrtdt = sqrt(dt);
g1 = force->ftm2v;
if (noise_flag == 0) {
g2 = 0.0;
} else if (gaussian_noise_flag == 1) {
g2 = sqrt(2 * force->boltz * temp / dt / force->mvv2e);
g2 = sqrt(2 * force->boltz / dt / force->mvv2e);
} else {
g2 = sqrt(24 * force->boltz * temp / dt / force->mvv2e);
g2 = sqrt(24 * force->boltz / dt / force->mvv2e);
}
}

47
src/BROWNIAN/fix_brownian_base.h
View File

@ -21,10 +21,10 @@ namespace LAMMPS_NS {
class FixBrownianBase : public Fix {
public:
FixBrownianBase(class LAMMPS *, int, char **);
virtual ~FixBrownianBase();
void init();
int setmask();
void reset_dt();
~FixBrownianBase() override;
void init() override;
int setmask() override;
void reset_dt() override;
protected:
int seed; // RNG seed
@ -33,6 +33,8 @@ class FixBrownianBase : public Fix {
int gamma_r_flag; // 0/1 if isotropic rotational damping is unset/set
int gamma_t_eigen_flag; // 0/1 if anisotropic translational damping is unset/set
int gamma_r_eigen_flag; // 0/1 if anisotropic rotational damping is unset/set
int rot_temp_flag; // 0/1 if rotational temperature is unset/set
int planar_rot_flag; // 0/1 if rotation is constrained to 2D (xy) plane
double gamma_t, gamma_r; // translational and rotational (isotropic) damping params
double *gamma_t_inv; // anisotropic damping parameter eigenvalues
@ -46,43 +48,12 @@ class FixBrownianBase : public Fix {
int noise_flag; // 0/1 for noise off/on
int gaussian_noise_flag; // 0/1 for uniform/gaussian noise
double temp; // temperature
double g1, g2; // prefactors in time stepping
double temp; // temperature
double rot_temp; // temperature
double g1, g2; // prefactors in time stepping
class RanMars *rng;
};
} // namespace LAMMPS_NS
#endif
/* ERROR/WARNING messages:
E: Illegal fix brownian command.
Wrong number/type of input arguments.
E: Fix brownian gamma_t_eigen values must be > 0.
Self-explanatory.
E: Fix brownian gamma_r_eigen values must be > 0.
Self-explanatory.
E: Fix brownian seed must be > 0.
Self-explanatory.
E: Fix brownian temp must be > 0.
Self-explanatory.
E: Fix brownian gamma_t must be > 0.
Self-explanatory.
E: Fix brownian gamma_r must be > 0.
Self-explanatory.
*/

45
src/BROWNIAN/fix_brownian_sphere.cpp
View File

@ -41,6 +41,7 @@ FixBrownianSphere::FixBrownianSphere(LAMMPS *lmp, int narg, char **arg) :
if (!gamma_t_flag || !gamma_r_flag) { error->all(FLERR, "Illegal fix brownian command."); }
if (!atom->mu_flag) error->all(FLERR, "Fix brownian/sphere requires atom attribute mu");
if (!atom->sphere_flag) error->all(FLERR, "Fix brownian/sphere requires atom style sphere");
}
/* ---------------------------------------------------------------------- */
@ -50,9 +51,9 @@ void FixBrownianSphere::init()
FixBrownianBase::init();
g3 = g1 / gamma_r;
g4 = g2 / sqrt(gamma_r);
g4 = g2 * sqrt(rot_temp / gamma_r);
g1 /= gamma_t;
g2 /= sqrt(gamma_t);
g2 *= sqrt(temp / gamma_t);
}
/* ---------------------------------------------------------------------- */
@ -61,19 +62,27 @@ void FixBrownianSphere::initial_integrate(int /*vflag */)
{
if (domain->dimension == 2) {
if (!noise_flag) {
initial_integrate_templated<0, 0, 1>();
initial_integrate_templated<0, 0, 1, 0>();
} else if (gaussian_noise_flag) {
initial_integrate_templated<0, 1, 1>();
initial_integrate_templated<0, 1, 1, 0>();
} else {
initial_integrate_templated<1, 0, 1>();
initial_integrate_templated<1, 0, 1, 0>();
}
} else if (planar_rot_flag) {
if (!noise_flag) {
initial_integrate_templated<0, 0, 0, 1>();
} else if (gaussian_noise_flag) {
initial_integrate_templated<0, 1, 0, 1>();
} else {
initial_integrate_templated<1, 0, 0, 1>();
}
} else {
if (!noise_flag) {
initial_integrate_templated<0, 0, 0>();
initial_integrate_templated<0, 0, 0, 0>();
} else if (gaussian_noise_flag) {
initial_integrate_templated<0, 1, 0>();
initial_integrate_templated<0, 1, 0, 0>();
} else {
initial_integrate_templated<1, 0, 0>();
initial_integrate_templated<1, 0, 0, 0>();
}
}
return;
@ -81,7 +90,7 @@ void FixBrownianSphere::initial_integrate(int /*vflag */)
/* ---------------------------------------------------------------------- */
template <int Tp_UNIFORM, int Tp_GAUSS, int Tp_2D>
template <int Tp_UNIFORM, int Tp_GAUSS, int Tp_2D, int Tp_2Drot>
void FixBrownianSphere::initial_integrate_templated()
{
double **x = atom->x;
@ -116,6 +125,24 @@ void FixBrownianSphere::initial_integrate_templated()
dy = dt * g1 * f[i][1];
wz = 0;
}
} else if (Tp_2Drot) {
wx = wy = 0;
if (Tp_UNIFORM) {
dx = dt * (g1 * f[i][0] + g2 * (rng->uniform() - 0.5));
dy = dt * (g1 * f[i][1] + g2 * (rng->uniform() - 0.5));
dz = dt * (g1 * f[i][2] + g2 * (rng->uniform() - 0.5));
wz = (rng->uniform() - 0.5) * g4;
} else if (Tp_GAUSS) {
dx = dt * (g1 * f[i][0] + g2 * rng->gaussian());
dy = dt * (g1 * f[i][1] + g2 * rng->gaussian());
dz = dt * (g1 * f[i][2] + g2 * rng->gaussian());
wz = rng->gaussian() * g4;
} else {
dx = dt * g1 * f[i][0];
dy = dt * g1 * f[i][1];
dz = dt * g1 * f[i][2];
wz = 0;
}
} else {
if (Tp_UNIFORM) {
dx = dt * (g1 * f[i][0] + g2 * (rng->uniform() - 0.5));

43
src/BROWNIAN/fix_brownian_sphere.h
View File

@ -27,48 +27,15 @@ namespace LAMMPS_NS {
class FixBrownianSphere : public FixBrownianBase {
public:
FixBrownianSphere(class LAMMPS *, int, char **);
virtual ~FixBrownianSphere(){};
void init();
void initial_integrate(int);
void init() override;
void initial_integrate(int) override;
private:
template <int Tp_UNIFORM, int Tp_GAUSS, int Tp_2D> void initial_integrate_templated();
template <int Tp_UNIFORM, int Tp_GAUSS, int Tp_2D, int Tp_2Drot>
void initial_integrate_templated();
double g3, g4;
};
} // namespace LAMMPS_NS
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal fix brownian/sphere command.
Wrong number/type of input arguments.
E: Compute brownian/sphere requires atom style sphere
Self-explanatory.
E: Compute brownian/sphere requires atom attribute mu
Self-explanatory.
E: Fix brownian/sphere translational viscous drag coefficient must be > 0.
Self-explanatory.
E: Fix brownian/sphere rotational viscous drag coefficient must be > 0.
Self-explanatory.
E: Fix brownian/sphere translational diffusion coefficient must be > 0.
Self-explanatory.
E: Fix brownian/sphere rotational diffusion coefficient must be > 0.
Self-explanatory.
E: Fix brownian/sphere seed must be > 0.
*/

2
src/BROWNIAN/fix_propel_self.cpp
View File

@ -100,7 +100,7 @@ void FixPropelSelf::init()
error->all(FLERR, "Fix propel/self requires atom attribute mu with option dipole");
if (mode == QUAT) {
avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
avec = dynamic_cast<AtomVecEllipsoid *>(atom->style_match("ellipsoid"));
if (!avec) error->all(FLERR, "Fix propel/self requires atom style ellipsoid with option quat");
// check that all particles are finite-size ellipsoids

30
src/BROWNIAN/fix_propel_self.h
View File

@ -26,11 +26,11 @@ namespace LAMMPS_NS {
class FixPropelSelf : public Fix {
public:
FixPropelSelf(class LAMMPS *, int, char **);
virtual ~FixPropelSelf(){};
void init();
void post_force(int);
void setup(int);
int setmask();
void init() override;
void post_force(int) override;
void setup(int) override;
int setmask() override;
private:
double magnitude;
@ -46,23 +46,3 @@ class FixPropelSelf : public Fix {
} // namespace LAMMPS_NS
#endif
#endif
/* ERROR/WARNING messages:
E: Illegal fix propel/self command.
Wrong number/type of input arguments.
E: Fix propel/self requires atom attribute mu with option dipole.
Self-explanatory.
E: Fix propel/self requires atom style ellipsoid with option quat.
Self-explanatory.
Fix propel/self requires extended particles with option quat.
Self-explanatory.
*/

19
src/CG-DNA/atom_vec_oxdna.cpp
View File

@ -33,18 +33,13 @@ AtomVecOxdna::AtomVecOxdna(LAMMPS *lmp) : AtomVec(lmp)
// order of fields in a string does not matter
// except: fields_data_atom & fields_data_vel must match data file
fields_grow = (char *) "id5p";
fields_copy = (char *) "id5p";
fields_comm = (char *) "";
fields_comm_vel = (char *) "";
fields_reverse = (char *) "";
fields_border = (char *) "id5p";
fields_border_vel = (char *) "";
fields_exchange = (char *) "id5p";
fields_restart = (char *) "id5p";
fields_create = (char *) "";
fields_data_atom = (char *) "id type x";
fields_data_vel = (char *) "id v";
fields_grow = {"id5p"};
fields_copy = {"id5p"};
fields_border = {"id5p"};
fields_exchange = {"id5p"};
fields_restart = {"id5p"};
fields_data_atom = {"id", "type", "x"};
fields_data_vel = {"id", "v"};
setup_fields();

11
src/CG-DNA/atom_vec_oxdna.h
View File

@ -27,11 +27,10 @@ namespace LAMMPS_NS {
class AtomVecOxdna : public AtomVec {
public:
AtomVecOxdna(class LAMMPS *);
~AtomVecOxdna() = default;
void grow_pointers();
void data_atom_post(int);
void data_bonds_post(int, int, tagint, tagint, tagint);
void grow_pointers() override;
void data_atom_post(int) override;
void data_bonds_post(int, int, tagint, tagint, tagint) override;
private:
tagint *id5p;
@ -41,7 +40,3 @@ class AtomVecOxdna : public AtomVec {
#endif
#endif
/* ERROR/WARNING messages:
*/

30
src/CG-DNA/bond_oxdna2_fene.h
View File

@ -27,38 +27,10 @@ namespace LAMMPS_NS {
class BondOxdna2Fene : public BondOxdnaFene {
public:
BondOxdna2Fene(class LAMMPS *lmp) : BondOxdnaFene(lmp) {}
virtual ~BondOxdna2Fene() {}
virtual void compute_interaction_sites(double *, double *, double *, double *) const;
void compute_interaction_sites(double *, double *, double *, double *) const override;
};
} // namespace LAMMPS_NS
#endif
#endif
/* ERROR/WARNING messages:
W: FENE bond too long: %ld %d %d %g
A FENE bond has stretched dangerously far. It's interaction strength
will be truncated to attempt to prevent the bond from blowing up.
E: Bad FENE bond
Two atoms in a FENE bond have become so far apart that the bond cannot
be computed.
E: Incorrect args for bond coefficients
Self-explanatory. Check the input script or data file.
W: Use special bonds = 0,1,1 with bond style oxdna
Most FENE models need this setting for the special_bonds command.
W: FENE bond too long: %ld %g
A FENE bond has stretched dangerously far. It's interaction strength
will be truncated to attempt to prevent the bond from blowing up.
*/

40
src/CG-DNA/bond_oxdna_fene.cpp
View File

@ -24,8 +24,8 @@
#include "neighbor.h"
#include "update.h"
#include "atom_vec_ellipsoid.h"
#include "math_extra.h"
#include "pair.h"
#include <cmath>
@ -152,18 +152,14 @@ void BondOxdnaFene::compute(int eflag, int vflag)
double r, rr0, rr0sq;
// vectors COM-backbone site in lab frame
double ra_cs[3], rb_cs[3];
double *qa, ax[3], ay[3], az[3];
double *qb, bx[3], by[3], bz[3];
// Cartesian unit vectors in lab frame
double ax[3], ay[3], az[3];
double bx[3], by[3], bz[3];
double **x = atom->x;
double **f = atom->f;
double **torque = atom->torque;
AtomVecEllipsoid *avec = (AtomVecEllipsoid *) atom->style_match("ellipsoid");
AtomVecEllipsoid::Bonus *bonus = avec->bonus;
int *ellipsoid = atom->ellipsoid;
int **bondlist = neighbor->bondlist;
int nbondlist = neighbor->nbondlist;
int nlocal = atom->nlocal;
@ -172,6 +168,12 @@ void BondOxdnaFene::compute(int eflag, int vflag)
ebond = 0.0;
ev_init(eflag, vflag);
// n(x/y/z)_xtrct = extracted local unit vectors in lab frame from oxdna_excv
int dim;
nx_xtrct = (double **) force->pair->extract("nx", dim);
ny_xtrct = (double **) force->pair->extract("ny", dim);
nz_xtrct = (double **) force->pair->extract("nz", dim);
// loop over FENE bonds
for (in = 0; in < nbondlist; in++) {
@ -180,10 +182,24 @@ void BondOxdnaFene::compute(int eflag, int vflag)
b = bondlist[in][0];
type = bondlist[in][2];
qa = bonus[ellipsoid[a]].quat;
MathExtra::q_to_exyz(qa, ax, ay, az);
qb = bonus[ellipsoid[b]].quat;
MathExtra::q_to_exyz(qb, bx, by, bz);
ax[0] = nx_xtrct[a][0];
ax[1] = nx_xtrct[a][1];
ax[2] = nx_xtrct[a][2];
ay[0] = ny_xtrct[a][0];
ay[1] = ny_xtrct[a][1];
ay[2] = ny_xtrct[a][2];
az[0] = nz_xtrct[a][0];
az[1] = nz_xtrct[a][1];
az[2] = nz_xtrct[a][2];
bx[0] = nx_xtrct[b][0];
bx[1] = nx_xtrct[b][1];
bx[2] = nx_xtrct[b][2];
by[0] = ny_xtrct[b][0];
by[1] = ny_xtrct[b][1];
by[2] = ny_xtrct[b][2];
bz[0] = nz_xtrct[b][0];
bz[1] = nz_xtrct[b][1];
bz[2] = nz_xtrct[b][2];
// vector COM-backbone site a and b
compute_interaction_sites(ax, ay, az, ra_cs);

48
src/CG-DNA/bond_oxdna_fene.h
View File

@ -27,19 +27,20 @@ namespace LAMMPS_NS {
class BondOxdnaFene : public Bond {
public:
BondOxdnaFene(class LAMMPS *lmp) : Bond(lmp) {}
virtual ~BondOxdnaFene();
~BondOxdnaFene() override;
virtual void compute_interaction_sites(double *, double *, double *, double *) const;
virtual void compute(int, int);
void coeff(int, char **);
void init_style();
double equilibrium_distance(int);
void write_restart(FILE *);
void read_restart(FILE *);
void write_data(FILE *);
double single(int, double, int, int, double &);
void compute(int, int) override;
void coeff(int, char **) override;
void init_style() override;
double equilibrium_distance(int) override;
void write_restart(FILE *) override;
void read_restart(FILE *) override;
void write_data(FILE *) override;
double single(int, double, int, int, double &) override;
protected:
double *k, *Delta, *r0; // FENE
double *k, *Delta, *r0; // FENE
double **nx_xtrct, **ny_xtrct, **nz_xtrct; // per-atom arrays for local unit vectors
void allocate();
void ev_tally_xyz(int, int, int, int, double, double, double, double, double, double, double);
@ -49,30 +50,3 @@ class BondOxdnaFene : public Bond {
#endif
#endif
/* ERROR/WARNING messages:
W: FENE bond too long: %ld %d %d %g
A FENE bond has stretched dangerously far. It's interaction strength
will be truncated to attempt to prevent the bond from blowing up.
E: Bad FENE bond
Two atoms in a FENE bond have become so far apart that the bond cannot
be computed.
E: Incorrect args for bond coefficients
Self-explanatory. Check the input script or data file.
W: Use special bonds = 0,1,1 with bond style oxdna
Most FENE models need this setting for the special_bonds command.
W: FENE bond too long: %ld %g
A FENE bond has stretched dangerously far. It's interaction strength
will be truncated to attempt to prevent the bond from blowing up.
*/

30
src/CG-DNA/bond_oxrna2_fene.h
View File

@ -27,39 +27,11 @@ namespace LAMMPS_NS {
class BondOxrna2Fene : public BondOxdnaFene {
public:
BondOxrna2Fene(class LAMMPS *lmp) : BondOxdnaFene(lmp) {}
virtual ~BondOxrna2Fene() {}
virtual void compute_interaction_sites(double *, double *, double *, double *) const;
void compute_interaction_sites(double *, double *, double *, double *) const override;
};
} // namespace LAMMPS_NS
#endif
#endif
/* ERROR/WARNING messages:
W: FENE bond too long: %ld %d %d %g
A FENE bond has stretched dangerously far. It's interaction strength
will be truncated to attempt to prevent the bond from blowing up.
E: Bad FENE bond
Two atoms in a FENE bond have become so far apart that the bond cannot
be computed.
E: Incorrect args for bond coefficients
Self-explanatory. Check the input script or data file.
W: Use special bonds = 0,1,1 with bond style oxrna
Most FENE models need this setting for the special_bonds command.
W: FENE bond too long: %ld %g
A FENE bond has stretched dangerously far. It's interaction strength
will be truncated to attempt to prevent the bond from blowing up.
*/

Some files were not shown because too many files have changed in this diff Show More