diff --git a/doc/src/Howto_bpm.rst b/doc/src/Howto_bpm.rst index 5929dbc962..893486c989 100644 --- a/doc/src/Howto_bpm.rst +++ b/doc/src/Howto_bpm.rst @@ -3,7 +3,7 @@ Bonded particle models The BPM package implements bonded particle models which can be used to simulate mesoscale solids. Solids are constructed as a collection of - + particles which each represent a coarse-grained region of space much larger than the atomistic scale. Particles within a solid region are then connected by a network of bonds to provide solid elasticity. @@ -24,7 +24,7 @@ Bonds can be created using a :doc:`read data ` or :doc:`create bond ` command. Alternatively, a :doc:`molecule ` template with bonds can be used with :doc:`fix deposit ` or :doc:`fix pour ` to -create solid grains. +create solid grains. In this implementation, bonds store their reference state when they are first computed in the setup of the first simulation run. Data is then @@ -47,7 +47,7 @@ this, LAMMPS requires :doc:`newton ` bond off such that all processors containing an atom know when a bond breaks. Additionally, one must do either (A) or (B). -(A) +(A) Use the following special bond settings @@ -64,7 +64,7 @@ charges in BPM models, setting a nonzero coul weight for 1-2 bonds ensures all bonded neighbors are still included in the neighbor list in case bonds break between neighbor list builds. -(B) +(B) Alternatively, one can simply overlay pair interactions such that all bonded particles also feel pair interactions. This can be accomplished @@ -109,13 +109,14 @@ velocity damping as its sister bond style. ---------- -While LAMMPS has many untilites to create and delete bonds, only a -select few are compatible with BPM bond styles. They include: +While LAMMPS has many untilites to create and delete bonds, the +following are currently compatible with BPM bond styles: * :doc:`create_bonds ` * :doc:`delete_bonds ` * :doc:`fix bond/create ` * :doc:`fix bond/break ` +* :doc:`fix bond/swap ` Note :doc:`bond_create ` requires certain special_bonds settings. To subtract pair interactions, one will need to switch between different diff --git a/doc/src/Howto_broken_bonds.rst b/doc/src/Howto_broken_bonds.rst index f7a6df6580..bf82088dec 100755 --- a/doc/src/Howto_broken_bonds.rst +++ b/doc/src/Howto_broken_bonds.rst @@ -31,7 +31,9 @@ page. Bonds can also be broken by fixes which change bond topology, including :doc:`fix bond/break ` and -:doc:`fix bond/react `. +:doc:`fix bond/react `. These fixes will automatically +trigger a rebuild of the neighbor list and update special bond data structures +when bonds are broken. Note that when bonds are dumped to a file via the :doc:`dump local ` command, bonds with type 0 are not included. The :doc:`delete_bonds ` command can also be used to query the diff --git a/doc/src/bond_bpm_rotational.rst b/doc/src/bond_bpm_rotational.rst index ad2c4ad265..fd6b448f24 100644 --- a/doc/src/bond_bpm_rotational.rst +++ b/doc/src/bond_bpm_rotational.rst @@ -150,13 +150,18 @@ status of broken bonds or permanently delete them, e.g.: ---------- -Restart +Restart and other info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" This bond style writes the reference state of each bond to :doc:`binary restart files `. Loading a restart file will properly resume bonds. +The single() function of these pair styles returns 0.0 for the energy +of a pairwise interaction, since energy is not conserved in these +dissipative potentials. It also returns only the normal component of +the pairwise interaction force. + Restrictions """""""""""" diff --git a/doc/src/bond_bpm_spring.rst b/doc/src/bond_bpm_spring.rst index c0624d7d35..4bd73d6af0 100644 --- a/doc/src/bond_bpm_spring.rst +++ b/doc/src/bond_bpm_spring.rst @@ -113,13 +113,17 @@ query the status of broken bonds or permanently delete them, e.g.: ---------- -Restart +Restart and other info """"""""""""""""""""""""""""""""""""""""""""""""""""""""""" This bond style writes the reference state of each bond to :doc:`binary restart files `. Loading a restart file will properly resume bonds. +The single() function of these pair styles returns 0.0 for the energy +of a pairwise interaction, since energy is not conserved in these +dissipative potentials. + Restrictions """""""""""" diff --git a/examples/bpm/impact/in.bpm.impact.rotational b/examples/bpm/impact/in.bpm.impact.rotational index 3c4667ca1b..34443fca6f 100644 --- a/examples/bpm/impact/in.bpm.impact.rotational +++ b/examples/bpm/impact/in.bpm.impact.rotational @@ -1,51 +1,52 @@ -units lj -dimension 3 -boundary f f f -atom_style sphere/bpm -special_bonds lj 0.0 1.0 1.0 coul 0.0 1.0 1.0 -newton on off -comm_modify vel yes cutoff 2.6 -lattice fcc 1.0 -region box block -25 15 -22 22 -22 22 -create_box 1 box bond/types 2 extra/bond/per/atom 20 extra/special/per/atom 50 - -region disk cylinder x 0.0 0.0 20.0 -0.5 0.5 -create_atoms 1 region disk -group plate region disk - -region ball sphere 8.0 0.0 0.0 6.0 -create_atoms 1 region ball -group projectile region ball - -displace_atoms all random 0.1 0.1 0.1 134598738 - -neighbor 1.0 bin -pair_style gran/hooke/history 1.0 NULL 0.5 NULL 0.1 1 -bond_style bpm/rotational store/local 2 time id1 id2 -pair_coeff 1 1 -bond_coeff 1 1.0 0.2 0.02 0.02 0.05 0.01 0.01 0.01 0.1 0.02 0.002 0.002 -bond_coeff 2 1.0 0.2 0.02 0.02 0.20 0.04 0.04 0.04 0.1 0.02 0.002 0.002 - -fix 1 all nve/sphere/bpm -fix 2 all store/local 100 3 - -create_bonds many plate plate 1 0.0 1.5 -create_bonds many projectile projectile 2 0.0 1.5 - -neighbor 0.3 bin -special_bonds lj 0.0 1.0 1.0 coul 1.0 1.0 1.0 - -velocity projectile set -0.05 0.0 0.0 -compute nbond all nbond/atom -compute tbond all reduce sum c_nbond - -timestep 0.05 -thermo_style custom step ke pe pxx pyy pzz c_tbond -thermo 100 -thermo_modify lost ignore lost/bond ignore -#dump 1 all custom 100 atomDump id radius x y z c_nbond - -dump 2 all local 100 brokenDump f_2[1] f_2[2] f_2[3] -dump_modify 2 header no - -run 7500 +units lj +dimension 3 +boundary f f f +atom_style sphere/bpm +special_bonds lj 0.0 1.0 1.0 coul 0.0 1.0 1.0 +newton on off +comm_modify vel yes cutoff 2.6 +lattice fcc 1.0 +region box block -25 15 -22 22 -22 22 +create_box 1 box bond/types 2 extra/bond/per/atom 20 extra/special/per/atom 50 + +region disk cylinder x 0.0 0.0 20.0 -0.5 0.5 +create_atoms 1 region disk +group plate region disk + +region ball sphere 8.0 0.0 0.0 6.0 +create_atoms 1 region ball +group projectile region ball + +displace_atoms all random 0.1 0.1 0.1 134598738 + +neighbor 1.0 bin +pair_style gran/hooke/history 1.0 NULL 0.5 NULL 0.1 1 +pair_coeff 1 1 + +fix 1 all nve/sphere/bpm +fix 2 all store/local 100 3 + +create_bonds many plate plate 1 0.0 1.5 +create_bonds many projectile projectile 2 0.0 1.5 + +neighbor 0.3 bin +special_bonds lj 0.0 1.0 1.0 coul 1.0 1.0 1.0 + +bond_style bpm/rotational store/local 2 time id1 id2 +bond_coeff 1 1.0 0.2 0.02 0.02 0.05 0.01 0.01 0.01 0.1 0.02 0.002 0.002 +bond_coeff 2 1.0 0.2 0.02 0.02 0.20 0.04 0.04 0.04 0.1 0.02 0.002 0.002 + +velocity projectile set -0.05 0.0 0.0 +compute nbond all nbond/atom +compute tbond all reduce sum c_nbond + +timestep 0.05 +thermo_style custom step ke pe pxx pyy pzz c_tbond +thermo 100 +thermo_modify lost ignore lost/bond ignore +#dump 1 all custom 100 atomDump id radius x y z c_nbond + +dump 2 all local 100 brokenDump f_2[1] f_2[2] f_2[3] +dump_modify 2 header no + +run 7500 diff --git a/examples/bpm/impact/in.bpm.impact.spring b/examples/bpm/impact/in.bpm.impact.spring index ef5506218c..60822cde01 100644 --- a/examples/bpm/impact/in.bpm.impact.spring +++ b/examples/bpm/impact/in.bpm.impact.spring @@ -1,53 +1,54 @@ -units lj -dimension 3 -boundary f f f -atom_style bond -special_bonds lj 0.0 1.0 1.0 coul 0.0 1.0 1.0 -newton on off -comm_modify vel yes cutoff 2.6 -lattice fcc 1.0 -region box block -25 15 -22 22 -22 22 -create_box 1 box bond/types 2 extra/bond/per/atom 20 extra/special/per/atom 50 - -region disk cylinder x 0.0 0.0 20.0 -0.5 0.5 -create_atoms 1 region disk -group plate region disk - -region ball sphere 8.0 0.0 0.0 6.0 -create_atoms 1 region ball -group projectile region ball - -displace_atoms all random 0.1 0.1 0.1 134598738 - -mass 1 1.0 - -neighbor 1.0 bin -pair_style bpm/spring -bond_style bpm/spring store/local 2 time id1 id2 -pair_coeff 1 1 1.0 1.0 1.0 -bond_coeff 1 1.0 0.04 1.0 -bond_coeff 2 1.0 0.20 1.0 - -fix 1 all nve -fix 2 all store/local 100 3 - -create_bonds many plate plate 1 0.0 1.5 -create_bonds many projectile projectile 2 0.0 1.5 - -neighbor 0.3 bin -special_bonds lj 0.0 1.0 1.0 coul 1.0 1.0 1.0 - -velocity projectile set -0.05 0.0 0.0 -compute nbond all nbond/atom -compute tbond all reduce sum c_nbond - -timestep 0.1 -thermo_style custom step ke pe pxx pyy pzz c_tbond -thermo 100 -thermo_modify lost ignore lost/bond ignore -#dump 1 all custom 100 atomDump id x y z c_nbond - -dump 2 all local 100 brokenDump f_2[1] f_2[2] f_2[3] -dump_modify 2 header no - -run 7500 +units lj +dimension 3 +boundary f f f +atom_style bond +special_bonds lj 0.0 1.0 1.0 coul 0.0 1.0 1.0 +newton on off +comm_modify vel yes cutoff 2.6 +lattice fcc 1.0 +region box block -25 15 -22 22 -22 22 +create_box 1 box bond/types 2 extra/bond/per/atom 20 extra/special/per/atom 50 + +region disk cylinder x 0.0 0.0 20.0 -0.5 0.5 +create_atoms 1 region disk +group plate region disk + +region ball sphere 8.0 0.0 0.0 6.0 +create_atoms 1 region ball +group projectile region ball + +displace_atoms all random 0.1 0.1 0.1 134598738 + +mass 1 1.0 + +neighbor 1.0 bin +pair_style bpm/spring +pair_coeff 1 1 1.0 1.0 1.0 + +fix 1 all nve +fix 2 all store/local 100 3 + +create_bonds many plate plate 1 0.0 1.5 +create_bonds many projectile projectile 2 0.0 1.5 + +neighbor 0.3 bin +special_bonds lj 0.0 1.0 1.0 coul 1.0 1.0 1.0 + +bond_style bpm/spring store/local 2 time id1 id2 +bond_coeff 1 1.0 0.04 1.0 +bond_coeff 2 1.0 0.20 1.0 + +velocity projectile set -0.05 0.0 0.0 +compute nbond all nbond/atom +compute tbond all reduce sum c_nbond + +timestep 0.1 +thermo_style custom step ke pe pxx pyy pzz c_tbond +thermo 100 +thermo_modify lost ignore lost/bond ignore +#dump 1 all custom 100 atomDump id x y z c_nbond + +dump 2 all local 100 brokenDump f_2[1] f_2[2] f_2[3] +dump_modify 2 header no + +run 7500 diff --git a/examples/bpm/pour/in.bpm.pour b/examples/bpm/pour/in.bpm.pour index 98dca9d3a7..c6bb288d30 100644 --- a/examples/bpm/pour/in.bpm.pour +++ b/examples/bpm/pour/in.bpm.pour @@ -1,35 +1,35 @@ -units lj -dimension 3 -boundary m m m -atom_style sphere/bpm -special_bonds lj 0.0 1.0 1.0 coul 1.0 1.0 1.0 -newton on off -comm_modify vel yes cutoff 3.3 -region box block -15 15 -15 15 0 60.0 -create_box 1 box bond/types 1 extra/bond/per/atom 15 extra/special/per/atom 50 - -molecule my_mol "rect.mol" -region wall_cyl cylinder z 0.0 0.0 10.0 EDGE EDGE side in -region dropzone cylinder z 0.0 0.0 10.0 40.0 50.0 side in - -pair_style gran/hertz/history 1.0 NULL 0.5 NULL 0.1 1 -bond_style bpm/rotational -pair_coeff 1 1 -bond_coeff 1 1.0 0.2 0.01 0.01 2.0 0.4 0.02 0.02 0.2 0.04 0.002 0.002 - -compute nbond all nbond/atom -compute tbond all reduce sum c_nbond -compute_modify thermo_temp dynamic/dof yes - -fix 1 all wall/gran hertz/history 1.0 NULL 0.5 NULL 0.1 1 zplane 0.0 NULL -fix 2 all wall/gran/region hertz/history 1.0 NULL 0.5 NULL 0.1 1 region wall_cyl -fix 3 all gravity 1e-4 vector 0 0 -1 -fix 4 all deposit 40 0 1500 712511343 mol my_mol region dropzone near 2.0 vz -0.05 -0.05 -fix 5 all nve/sphere/bpm - -timestep 0.05 -thermo_style custom step ke pe pxx pyy pzz c_tbond -thermo 100 -#dump 1 all custom 500 atomDump id radius x y z c_nbond mol - -run 100000 +units lj +dimension 3 +boundary m m m +atom_style sphere/bpm +special_bonds lj 0.0 1.0 1.0 coul 1.0 1.0 1.0 +newton on off +comm_modify vel yes cutoff 3.3 +region box block -15 15 -15 15 0 60.0 +create_box 1 box bond/types 1 extra/bond/per/atom 15 extra/special/per/atom 50 + +molecule my_mol "rect.mol" +region wall_cyl cylinder z 0.0 0.0 10.0 EDGE EDGE side in +region dropzone cylinder z 0.0 0.0 10.0 40.0 50.0 side in + +pair_style gran/hertz/history 1.0 NULL 0.5 NULL 0.1 1 +bond_style bpm/rotational +pair_coeff 1 1 +bond_coeff 1 1.0 0.2 0.01 0.01 2.0 0.4 0.02 0.02 0.2 0.04 0.002 0.002 + +compute nbond all nbond/atom +compute tbond all reduce sum c_nbond +compute_modify thermo_temp dynamic/dof yes + +fix 1 all wall/gran hertz/history 1.0 NULL 0.5 NULL 0.1 1 zplane 0.0 NULL +fix 2 all wall/gran/region hertz/history 1.0 NULL 0.5 NULL 0.1 1 region wall_cyl +fix 3 all gravity 1e-4 vector 0 0 -1 +fix 4 all deposit 40 0 1500 712511343 mol my_mol region dropzone near 2.0 vz -0.05 -0.05 +fix 5 all nve/sphere/bpm + +timestep 0.05 +thermo_style custom step ke pe pxx pyy pzz c_tbond +thermo 100 +#dump 1 all custom 500 atomDump id radius x y z c_nbond mol + +run 100000 diff --git a/src/.gitignore b/src/.gitignore index 19c5295576..0c6c893234 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -265,8 +265,6 @@ /compute_nbond_atom.h /fix_nve_sphere_bpm.cpp /fix_nve_sphere_bpm.h -/fix_update_special_bonds.cpp -/fix_update_special_bonds.h /pair_bpm_spring.cpp /pair_bpm_spring.h diff --git a/src/BPM/bond_bpm.cpp b/src/BPM/bond_bpm.cpp index 00a5ebc6e4..7ff883a140 100644 --- a/src/BPM/bond_bpm.cpp +++ b/src/BPM/bond_bpm.cpp @@ -25,6 +25,9 @@ #include "modify.h" #include "update.h" +#include +#include + using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ @@ -49,7 +52,7 @@ BondBPM::BondBPM(LAMMPS *lmp) : Bond(lmp) // create dummy fix as placeholder for FixUpdateSpecialBonds // this is so final order of Modify:fix will conform to input script - id_fix_dummy = utils::strdup("BPM_DUMMY"); + id_fix_dummy = utils::strdup("BPM_DUMMY" + std::to_string(instance_me)); modify->add_fix(fmt::format("{} all DUMMY ", id_fix_dummy)); } @@ -58,14 +61,16 @@ BondBPM::BondBPM(LAMMPS *lmp) : Bond(lmp) BondBPM::~BondBPM() { delete [] pack_choice; - delete [] id_fix_store_local; - delete [] id_fix_prop_atom; if (id_fix_dummy) modify->delete_fix(id_fix_dummy); if (id_fix_update) modify->delete_fix(id_fix_update); + 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_update; + delete [] id_fix_store_local; + delete [] id_fix_prop_atom; memory->destroy(output_data); } @@ -110,7 +115,7 @@ void BondBPM::init_style() 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"); + id_fix_update = utils::strdup("BPM_update_special_bonds" + std::to_string(instance_me)); fix_update_special_bonds = (FixUpdateSpecialBonds *) modify->replace_fix(id_fix_dummy, fmt::format("{} all UPDATE_SPECIAL_BONDS", id_fix_update),1); delete [] id_fix_dummy; @@ -141,11 +146,14 @@ void BondBPM::init_style() void BondBPM::settings(int narg, char **arg) { + leftover_args.clear(); + + int local_freq; int iarg = 0; while (iarg < narg) { if (strcmp(arg[iarg], "store/local") == 0) { - id_fix_store_local = utils::strdup(arg[iarg+1]); nvalues = 0; + local_freq = utils::inumeric(FLERR, arg[iarg+1], false, lmp); pack_choice = new FnPtrPack[narg - iarg - 1]; iarg += 2; while (iarg < narg) { @@ -179,31 +187,44 @@ void BondBPM::settings(int narg, char **arg) overlay_flag = 1; iarg ++; } else { - error->all(FLERR, "Illegal pair_style command"); + leftover_args.push_back(iarg); + iarg ++; } } - if (id_fix_store_local) { + if (nvalues != 0 && !id_fix_store_local) { + //Todo, assign ID and create fix id_fix_store_local = utils::strdup(arg[iarg+1]); + if (nvalues == 0) error->all(FLERR, "Bond style bpm/rotational must include at least one value to output"); memory->create(output_data, nvalues, "bond/bpm:output_data"); // Use store property 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"); + id_fix_prop_atom = utils::strdup("BPM_property_atom" + std::to_string(instance_me)); int ifix = modify->find_fix(id_fix_prop_atom); + + char *x_ref_id = utils::strdup("BPM_X_REF" + std::to_string(instance_me)); + char *y_ref_id = utils::strdup("BPM_Y_REF" + std::to_string(instance_me)); + char *z_ref_id = utils::strdup("BPM_Z_REF" + std::to_string(instance_me)); if (ifix < 0) { - modify->add_fix(fmt::format("{} all property/atom " - "d_BPM_X_REF d_BPM_Y_REF d_BPM_Z_REF ghost yes", id_fix_prop_atom)); + modify->add_fix(fmt::format("{} all property/atom {} {} {} ghost yes", + id_fix_prop_atom, x_ref_id, y_ref_id, z_ref_id)); ifix = modify->find_fix(id_fix_prop_atom); } int type_flag; int col_flag; - index_x_ref = atom->find_custom("BPM_X_REF", type_flag, col_flag); - index_y_ref = atom->find_custom("BPM_Y_REF", type_flag, col_flag); - index_z_ref = atom->find_custom("BPM_Z_REF", type_flag, 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 (modify->fix[ifix]->restart_reset) { modify->fix[ifix]->restart_reset = 0; diff --git a/src/BPM/bond_bpm.h b/src/BPM/bond_bpm.h index 4836a08a64..3f71aefa45 100644 --- a/src/BPM/bond_bpm.h +++ b/src/BPM/bond_bpm.h @@ -16,6 +16,8 @@ #include "bond.h" +#include + namespace LAMMPS_NS { class BondBPM : public Bond { @@ -36,6 +38,8 @@ class BondBPM : public Bond { double r0_max_estimate; double max_stretch; + std::vector leftover_args; + char *id_fix_dummy, *id_fix_update; char *id_fix_store_local, *id_fix_prop_atom; class FixStoreLocal *fix_store_local; diff --git a/src/BPM/bond_bpm_rotational.cpp b/src/BPM/bond_bpm_rotational.cpp index 0f62d88966..80fff550ee 100644 --- a/src/BPM/bond_bpm_rotational.cpp +++ b/src/BPM/bond_bpm_rotational.cpp @@ -176,6 +176,270 @@ void BondBPMRotational::store_data() 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) +{ + int m; + 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 = 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 (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::scale3(-0.5,s2); + + 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); // 0.5*r points from particle 2 to midpoint + 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) @@ -186,32 +450,18 @@ void BondBPMRotational::compute(int eflag, int vflag) store_data(); } - int i1,i2,itmp,m,n,type,itype,jtype; - double q1[4], q2[4], r[3], r0[3]; - double rsq, r0_mag, r_mag, r_mag_inv, Fr, Fs_mag; - double Tt_mag, Tb_mag, breaking, smooth; + int i1,i2,itmp,n,type; + double r[3], r0[3], rhat[3]; + double delx, dely, delz, rsq, r0_mag, r_mag, r_mag_inv; + double Fr, breaking, smooth; double force1on2[3], torque1on2[3], torque2on1[3]; - double rhat[3], wn1[3], wn2[3], wxn1[3], wxn2[3], vroll[3]; - double w1dotr, w2dotr, v1dotr, v2dotr; - double vn1[3], vn2[3], vt1[3], vt2[3], tmp[3], s1[3], s2[3], tdamp[3]; - double tor1, tor2, tor3, fs1, fs2, fs3; - - double q2inv[4], rb[3], rb_x_r0[3], s[3], t[3], Fs[3]; - double q21[4], qp21[4], Tbp[3], Ttp[3]; - double Tsp[3], Fsp[3], Tt[3], Tb[3], Ts[3], F_rot[3], T_rot[3]; - double mq[4], mqinv[4], Ttmp[3], Ftmp[3], qtmp[4]; - double r0_dot_rb, gamma, c, psi, theta, sin_phi, cos_phi, temp; - double mag_in_plane, mag_out_plane; ev_init(eflag,vflag); double **x = atom->x; - double **v = atom->v; - double **omega = atom->omega; double **f = atom->f; double **torque = atom->torque; double *radius = atom->radius; - double **quat = atom->quat; tagint *tag = atom->tag; int **bondlist = neighbor->bondlist; int nbondlist = neighbor->nbondlist; @@ -249,16 +499,6 @@ void BondBPMRotational::compute(int eflag, int vflag) r0[2] = bondstore[n][3]; MathExtra::scale3(r0_mag, r0); - 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]; - // Note this is the reverse of Mora & Wang MathExtra::sub3(x[i1], x[i2], r); @@ -268,160 +508,11 @@ void BondBPMRotational::compute(int eflag, int vflag) MathExtra::scale3(r_mag_inv, r, rhat); // ------------------------------------------------------// - // Calculate forces using formulation in: - // 1) Y. Wang Acta Geotechnica 2009 - // 2) P. Mora & Y. Wang Advances in Geomcomputing 2009 + // Calculate forces, check if bond breaks // ------------------------------------------------------// - // 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 = 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 (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); - - // ------------------------------------------------------// - // Check if bond breaks - // ------------------------------------------------------// - - breaking = Fr/Fcr[type] + Fs_mag/Fcs[type] + Tb_mag/Tcb[type] + Tt_mag/Tct[type]; - if (breaking < 0.0) breaking = 0.0; + 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; @@ -429,76 +520,11 @@ void BondBPMRotational::compute(int eflag, int vflag) continue; } + damping_forces(i1, i2, type, Fr, rhat, r, force1on2, torque1on2, torque2on1); + smooth = breaking*breaking; smooth = 1.0 - smooth*smooth; - // ------------------------------------------------------// - // Calculate damping using formulation in - // Y. Wang, F. Alonso-Marroquin, W. Guo 2015 - // ------------------------------------------------------// - // Note: n points towards 1 vs pointing towards 2 - - // 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); - 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::scale3(-0.5,s2); - - 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); // 0.5*r points from particle 2 to midpoint - 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); - // ------------------------------------------------------// // Apply forces and torques to particles // ------------------------------------------------------// @@ -617,7 +643,18 @@ void BondBPMRotational::init_style() if (!fix_bond_history) fix_bond_history = (FixBondHistory *) modify->add_fix( - "BOND_HISTORY_BPM_ROTATIONAL all BOND_HISTORY 0 4"); + "HISTORY_BPM_ROTATIONAL" + std::to_string(instance_me) + " all BOND_HISTORY 0 4"); +} + +/* ---------------------------------------------------------------------- */ + +void BondBPMRotational::settings(int narg, char **arg) +{ + BondBPM::settings(narg, arg); + + for (int iarg : leftover_args) { + error->all(FLERR, "Illegal bond_style command"); + } } /* ---------------------------------------------------------------------- @@ -698,10 +735,42 @@ double BondBPMRotational::single(int type, double rsq, int i, int j, // Not yet enabled 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); - // } - //} + 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; + + smooth = breaking*breaking; + smooth = 1.0 - smooth*smooth; + fforce *= smooth; + return 0.0; } diff --git a/src/BPM/bond_bpm_rotational.h b/src/BPM/bond_bpm_rotational.h index de3ea29240..150c7444bd 100644 --- a/src/BPM/bond_bpm_rotational.h +++ b/src/BPM/bond_bpm_rotational.h @@ -31,6 +31,7 @@ class BondBPMRotational : public BondBPM { virtual void compute(int, int); void coeff(int, char **); void init_style(); + void settings(int, char **); void write_restart(FILE *); void read_restart(FILE *); void write_data(FILE *); @@ -41,6 +42,11 @@ class BondBPMRotational : public BondBPM { double *Fcr, *Fcs, *Tct, *Tcb; 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); diff --git a/src/BPM/bond_bpm_spring.cpp b/src/BPM/bond_bpm_spring.cpp index be4b7d6fc3..021d552297 100644 --- a/src/BPM/bond_bpm_spring.cpp +++ b/src/BPM/bond_bpm_spring.cpp @@ -136,7 +136,7 @@ void BondBPMSpring::compute(int eflag, int vflag) int i1,i2,itmp,m,n,type,itype,jtype; double delx, dely, delz, delvx, delvy, delvz; - double e, rsq, r, r0, rinv, smooth, fbond, ebond, dot; + double e, rsq, r, r0, rinv, smooth, fbond, dot; ev_init(eflag,vflag); @@ -190,7 +190,6 @@ void BondBPMSpring::compute(int eflag, int vflag) rinv = 1.0/r; fbond = k[type]*(r0-r); - if (eflag) ebond = -0.5*fbond*(r0-r); delvx = v[i1][0] - v[i2][0]; delvy = v[i1][1] - v[i2][1]; @@ -218,7 +217,7 @@ void BondBPMSpring::compute(int eflag, int vflag) f[i2][2] -= delz*fbond; } - if (evflag) ev_tally(i1,i2,nlocal,newton_bond,ebond,fbond,delx,dely,delz); + if (evflag) ev_tally(i1,i2,nlocal,newton_bond,0.0,fbond,delx,dely,delz); } } @@ -280,7 +279,18 @@ void BondBPMSpring::init_style() if (!fix_bond_history) fix_bond_history = (FixBondHistory *) modify->add_fix( - "BOND_HISTORY_BPM_SPRING all BOND_HISTORY 0 1"); + "HISTORY_BPM_SPRING_" + std::to_string(instance_me) + " all BOND_HISTORY 0 1"); +} + +/* ---------------------------------------------------------------------- */ + +void BondBPMSpring::settings(int narg, char **arg) +{ + BondBPM::settings(narg, arg); + + for (int iarg : leftover_args) { + error->all(FLERR, "Illegal bond_style command"); + } } /* ---------------------------------------------------------------------- @@ -329,13 +339,37 @@ void BondBPMSpring::write_data(FILE *fp) double BondBPMSpring::single(int type, double rsq, int i, int j, double &fforce) { - // Not yet enabled 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 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; + double e = (r - r0)/r0; + 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; + + double smooth = (r-r0)/(r0*ecrit[type]); + smooth *= smooth; + smooth *= smooth; + smooth *= smooth; + smooth = 1 - smooth; + + fforce *= rinv*smooth; + return 0.0; } diff --git a/src/BPM/bond_bpm_spring.h b/src/BPM/bond_bpm_spring.h index 2b28723fc1..a81f49e81a 100644 --- a/src/BPM/bond_bpm_spring.h +++ b/src/BPM/bond_bpm_spring.h @@ -31,6 +31,7 @@ class BondBPMSpring : public BondBPM { virtual void compute(int, int); void coeff(int, char **); void init_style(); + void settings(int, char **); void write_restart(FILE *); void read_restart(FILE *); void write_data(FILE *); diff --git a/src/BPM/pair_bpm_spring.cpp b/src/BPM/pair_bpm_spring.cpp index abaf868313..f87ddd51e0 100644 --- a/src/BPM/pair_bpm_spring.cpp +++ b/src/BPM/pair_bpm_spring.cpp @@ -112,7 +112,6 @@ void PairBPMSpring::compute(int eflag, int vflag) rinv = 1.0/r; fpair = k[itype][jtype]*(cut[itype][jtype]-r); - if (eflag) evdwl = -0.5 * fpair * (cut[itype][jtype]-r) * factor_lj; smooth = rsq/cutsq[itype][jtype]; smooth *= smooth; @@ -125,6 +124,7 @@ void PairBPMSpring::compute(int eflag, int vflag) 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; @@ -307,7 +307,7 @@ double PairBPMSpring::single(int i, int j, int itype, int jtype, double rsq, double factor_coul, double factor_lj, double &fforce) { - double fpair,philj,r,rinv; + double fpair,r,rinv; double delx, dely, delz, delvx, delvy, delvz, dot, smooth; if(rsq > cutsq[itype][jtype]) return 0.0; @@ -319,7 +319,6 @@ double PairBPMSpring::single(int i, int j, int itype, int jtype, double rsq, rinv = 1.0/r; fpair = k[itype][jtype]*(cut[itype][jtype]-r); - philj = -0.5*k[itype][jtype]*(cut[itype][jtype]-r)*(cut[itype][jtype]-r); smooth = rsq/cutsq[itype][jtype]; smooth *= smooth; @@ -336,5 +335,5 @@ double PairBPMSpring::single(int i, int j, int itype, int jtype, double rsq, fpair *= factor_lj; fforce = fpair; - return philj; + return 0.0; } diff --git a/src/bond.cpp b/src/bond.cpp index 48972aed72..f3c5d3e341 100644 --- a/src/bond.cpp +++ b/src/bond.cpp @@ -28,6 +28,10 @@ using namespace LAMMPS_NS; enum{NONE,LINEAR,SPLINE}; +// allocate space for static class instance variable and initialize it + +int Bond::instance_total = 0; + /* ----------------------------------------------------------------------- set bond contribution to Vdwl energy to 0.0 a particular bond style can override this @@ -35,6 +39,8 @@ enum{NONE,LINEAR,SPLINE}; Bond::Bond(LAMMPS *lmp) : Pointers(lmp) { + instance_me = instance_total++; + energy = 0.0; virial[0] = virial[1] = virial[2] = virial[3] = virial[4] = virial[5] = 0.0; writedata = 1; diff --git a/src/bond.h b/src/bond.h index ea213da476..e0550ead95 100644 --- a/src/bond.h +++ b/src/bond.h @@ -23,6 +23,8 @@ class Bond : protected Pointers { friend class FixOMP; public: + static int instance_total; // # of Bond classes ever instantiated + int allocated; int *setflag; int partial_flag; // 1 if bond type can be set to 0 and deleted @@ -68,6 +70,7 @@ class Bond : protected Pointers { void write_file(int, char **); protected: + int instance_me; // which Bond class instantiation I am int suffix_flag; // suffix compatibility flag int evflag; diff --git a/src/fix_bond_history.cpp b/src/fix_bond_history.cpp index 6a187ca56b..dfadd4201d 100644 --- a/src/fix_bond_history.cpp +++ b/src/fix_bond_history.cpp @@ -308,4 +308,4 @@ void FixBondHistory::shift_bond(int i, int m, int k) int n = atom->num_bond[i]; for (int idata = 0; idata < ndata; idata ++) stored[i][m*ndata+idata] = stored[i][k*ndata+idata]; -} \ No newline at end of file +} diff --git a/src/fix_bond_history.h b/src/fix_bond_history.h index 58c40c2f36..970b1429d6 100644 --- a/src/fix_bond_history.h +++ b/src/fix_bond_history.h @@ -42,7 +42,7 @@ class FixBondHistory : public Fix { void update_atom_value(int, int, int, double); double get_atom_value(int, int, int); void delete_bond(int,int); - void shift_bond(int,int,int); + void shift_bond(int,int,int); double **bondstore; int stored_flag; diff --git a/src/fix_property_atom.cpp b/src/fix_property_atom.cpp index a971faa066..bb44300357 100644 --- a/src/fix_property_atom.cpp +++ b/src/fix_property_atom.cpp @@ -38,6 +38,7 @@ FixPropertyAtom::FixPropertyAtom(LAMMPS *lmp, int narg, char **arg) : restart_peratom = 1; wd_section = 1; + create_attribute = 1; int iarg = 3; nvalue = narg-iarg; @@ -554,24 +555,44 @@ void FixPropertyAtom::copy_arrays(int i, int j, int /*delflag*/) atom->q[j] = atom->q[i]; else if (styles[nv] == RMASS) atom->rmass[j] = atom->rmass[i]; - else if (styles[nv] == IVEC) { + else if (styles[nv] == IVEC) atom->ivector[index[nv]][j] = atom->ivector[index[nv]][i]; - atom->ivector[index[nv]][i] = 0; - } else if (styles[nv] == DVEC) { + else if (styles[nv] == DVEC) atom->dvector[index[nv]][j] = atom->dvector[index[nv]][i]; - atom->dvector[index[nv]][i] = 0.0; - } else if (styles[nv] == IARRAY) { + else if (styles[nv] == IARRAY) { ncol = cols[nv]; - for (k = 0; k < ncol; k++) { + for (k = 0; k < ncol; k++) atom->iarray[index[nv]][j][k] = atom->iarray[index[nv]][i][k]; - atom->iarray[index[nv]][i][k] = 0; - } } else if (styles[nv] == DARRAY) { ncol = cols[nv]; - for (k = 0; k < ncol; k++) { + for (k = 0; k < ncol; k++) atom->darray[index[nv]][j][k] = atom->darray[index[nv]][i][k]; + } + } +} + + +/* ---------------------------------------------------------------------- + initialize one atom's storage values, called when atom is created +------------------------------------------------------------------------- */ + +void FixPropertyAtom::set_arrays(int i) +{ + int k,ncol; + + for (int nv = 0; nv < nvalue; nv++) { + if (styles[nv] == IVEC) + atom->ivector[index[nv]][i] = 0; + else if (styles[nv] == DVEC) + atom->dvector[index[nv]][i] = 0.0; + else if (styles[nv] == IARRAY) { + ncol = cols[nv]; + for (k = 0; k < ncol; k++) + atom->iarray[index[nv]][i][k] = 0; + } else if (styles[nv] == DARRAY) { + ncol = cols[nv]; + for (k = 0; k < ncol; k++) atom->darray[index[nv]][i][k] = 0.0; - } } } } @@ -706,24 +727,16 @@ int FixPropertyAtom::pack_exchange(int i, double *buf) if (styles[nv] == MOLECULE) buf[m++] = ubuf(atom->molecule[i]).d; else if (styles[nv] == CHARGE) buf[m++] = atom->q[i]; else if (styles[nv] == RMASS) buf[m++] = atom->rmass[i]; - else if (styles[nv] == IVEC) { - buf[m++] = ubuf(atom->ivector[index[nv]][i]).d; - atom->ivector[index[nv]][i] = 0; - } else if (styles[nv] == DVEC) { - buf[m++] = atom->dvector[index[nv]][i]; - atom->dvector[index[nv]][i] = 0.0; - } else if (styles[nv] == IARRAY) { + else if (styles[nv] == IVEC) buf[m++] = ubuf(atom->ivector[index[nv]][i]).d; + else if (styles[nv] == DVEC) buf[m++] = atom->dvector[index[nv]][i]; + else if (styles[nv] == IARRAY) { ncol = cols[nv]; - for (k = 0; k < ncol; k++) { + for (k = 0; k < ncol; k++) buf[m++] = ubuf(atom->iarray[index[nv]][i][k]).d; - atom->iarray[index[nv]][i][k] = 0; - } } else if (styles[nv] == DARRAY) { ncol = cols[nv]; - for (k = 0; k < ncol; k++) { + for (k = 0; k < ncol; k++) buf[m++] = atom->darray[index[nv]][i][k]; - atom->darray[index[nv]][i][k] = 0.0; - } } } diff --git a/src/fix_property_atom.h b/src/fix_property_atom.h index 8580865e59..48472b8969 100644 --- a/src/fix_property_atom.h +++ b/src/fix_property_atom.h @@ -40,6 +40,7 @@ class FixPropertyAtom : public Fix { virtual void grow_arrays(int); void copy_arrays(int, int, int); + void set_arrays(int); int pack_border(int, int *, double *); int unpack_border(int, int, double *); int pack_exchange(int, double *); diff --git a/src/fix_store_local.h b/src/fix_store_local.h index 0be752ff52..591ead9486 100644 --- a/src/fix_store_local.h +++ b/src/fix_store_local.h @@ -13,7 +13,7 @@ #ifdef FIX_CLASS // clang-format off -FixStyle(store/local,FixStoreLocal); +FixStyle(STORE_LOCAL,FixStoreLocal); // clang-format on #else diff --git a/src/BPM/fix_update_special_bonds.cpp b/src/fix_update_special_bonds.cpp similarity index 64% rename from src/BPM/fix_update_special_bonds.cpp rename to src/fix_update_special_bonds.cpp index 19b72e06a7..b479619cf5 100644 --- a/src/BPM/fix_update_special_bonds.cpp +++ b/src/fix_update_special_bonds.cpp @@ -24,7 +24,7 @@ #include "pair.h" #include -#include +#include #include using namespace LAMMPS_NS; @@ -36,7 +36,6 @@ FixUpdateSpecialBonds::FixUpdateSpecialBonds(LAMMPS *lmp, int narg, char **arg) Fix(lmp, narg, arg) { if (narg != 3) error->all(FLERR,"Illegal fix update/special/bonds command"); - comm_forward = 1+atom->maxspecial; } /* ---------------------------------------------------------------------- */ @@ -76,6 +75,7 @@ void FixUpdateSpecialBonds::setup(int /*vflag*/) force->special_coul[3] != 1.0) error->all(FLERR,"Fix update/special/bonds requires special Coulomb weights = 1,1,1"); + new_broken_pairs.clear(); broken_pairs.clear(); } @@ -86,7 +86,7 @@ void FixUpdateSpecialBonds::setup(int /*vflag*/) void FixUpdateSpecialBonds::pre_exchange() { int i, j, key, m, n1, n3; - tagint min_tag, max_tag; + tagint tagi, tagj; int nlocal = atom->nlocal; tagint *tag = atom->tag; @@ -94,19 +94,19 @@ void FixUpdateSpecialBonds::pre_exchange() int **nspecial = atom->nspecial; tagint **special = atom->special; - for (auto const &key : broken_pairs) { - min_tag = key.first; - max_tag = key.second; + for (auto const &it : broken_pairs) { + tagi = it.first; + tagj = it.second; - i = atom->map(min_tag); - j = atom->map(max_tag); + i = atom->map(tagi); + j = atom->map(tagj); // remove i from special bond list for atom j and vice versa if (i < nlocal) { slist = special[i]; n1 = nspecial[i][0]; for (m = 0; m < n1; m++) - if (slist[m] == max_tag) break; + if (slist[m] == tagj) break; n3 = nspecial[i][2]; for (; m < n3-1; m++) slist[m] = slist[m+1]; nspecial[i][0]--; @@ -118,7 +118,7 @@ void FixUpdateSpecialBonds::pre_exchange() slist = special[j]; n1 = nspecial[j][0]; for (m = 0; m < n1; m++) - if (slist[m] == min_tag) break; + if (slist[m] == tagi) break; n3 = nspecial[j][2]; for (; m < n3-1; m++) slist[m] = slist[m+1]; nspecial[j][0]--; @@ -127,9 +127,6 @@ void FixUpdateSpecialBonds::pre_exchange() } } - // Forward updated special bond list - comm->forward_comm_fix(this); - broken_pairs.clear(); } @@ -139,82 +136,51 @@ void FixUpdateSpecialBonds::pre_exchange() void FixUpdateSpecialBonds::pre_force(int /*vflag*/) { - int i,j,n,m,ii,jj,inum,jnum; - int *ilist,*jlist,*numneigh,**firstneigh; - tagint min_tag, max_tag; - std::pair key; + int i1,i2,j,jj,jnum; + int *jlist,*numneigh,**firstneigh; + tagint tag1, tag2; - int **bond_type = atom->bond_type; - int *num_bond = atom->num_bond; - tagint **bond_atom = atom->bond_atom; int nlocal = atom->nlocal; tagint *tag = atom->tag; - NeighList *list = force->pair->list; - inum = list->inum; - ilist = list->ilist; + NeighList *list = force->pair->list; // may need to be generalized to work with pair hybrid* numneigh = list->numneigh; firstneigh = list->firstneigh; - for (ii = 0; ii < inum; ii++) { - i = ilist[ii]; - jlist = firstneigh[i]; - jnum = numneigh[i]; + // In theory could communicate a list of broken bonds to neighboring processors here + // to remove restriction that users use Newton bond off - for (jj = 0; jj < jnum; jj++) { - j = jlist[jj]; - j &= NEIGHMASK; - - min_tag = tag[i]; - max_tag = tag[j]; - if (max_tag < min_tag) { - min_tag = tag[j]; - max_tag = tag[i]; + for (auto const &it : new_broken_pairs) { + tag1 = it.first; + tag2 = it.second; + i1 = atom->map(tag1); + i2 = atom->map(tag2); + + // Loop through atoms of owned atoms i j + if (i1 < nlocal) { + jlist = firstneigh[i1]; + jnum = numneigh[i1]; + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= SPECIALMASK; // Clear special bond bits + if (tag[j] == tag2) + jlist[jj] = j; } - key = std::make_pair(min_tag, max_tag); + } - if (broken_pairs.find(key) != broken_pairs.end()) - jlist[jj] = j; // Clear special bond bits + if (i2 < nlocal) { + jlist = firstneigh[i2]; + jnum = numneigh[i2]; + for (jj = 0; jj < jnum; jj++) { + j = jlist[jj]; + j &= SPECIALMASK; // Clear special bond bits + if (tag[j] == tag1) + jlist[jj] = j; + } } } -} -/* ---------------------------------------------------------------------- */ - -int FixUpdateSpecialBonds::pack_forward_comm(int n, int *list, double *buf, - int /*pbc_flag*/, int * /*pbc*/) -{ - int i,j,k,m,ns; - int **nspecial = atom->nspecial; - tagint **special = atom->special; - - m = 0; - for (i = 0; i < n; i++) { - j = list[i]; - ns = nspecial[j][0]; - buf[m++] = ubuf(ns).d; - for (k = 0; k < ns; k++) - buf[m++] = ubuf(special[j][k]).d; - } - return m; -} - -/* ---------------------------------------------------------------------- */ - -void FixUpdateSpecialBonds::unpack_forward_comm(int n, int first, double *buf) -{ - int i,j,m,ns,last; - int **nspecial = atom->nspecial; - tagint **special = atom->special; - - m = 0; - last = first + n; - for (i = first; i < last; i++) { - ns = (int) ubuf(buf[m++]).i; - nspecial[i][0] = ns; - for (j = 0; j < ns; j++) - special[i][j] = (tagint) ubuf(buf[m++]).i; - } + new_broken_pairs.clear(); } /* ---------------------------------------------------------------------- */ @@ -222,14 +188,8 @@ void FixUpdateSpecialBonds::unpack_forward_comm(int n, int first, double *buf) void FixUpdateSpecialBonds::add_broken_bond(int i, int j) { tagint *tag = atom->tag; + std::pair tag_pair = std::make_pair(tag[i],tag[j]); - tagint min_tag = tag[i]; - tagint max_tag = tag[j]; - if (max_tag < min_tag) { - min_tag = tag[j]; - max_tag = tag[i]; - } - std::pair key = std::make_pair(min_tag, max_tag); - - broken_pairs.insert(key); + new_broken_pairs.push_back(tag_pair); + broken_pairs.push_back(tag_pair); } diff --git a/src/BPM/fix_update_special_bonds.h b/src/fix_update_special_bonds.h similarity index 85% rename from src/BPM/fix_update_special_bonds.h rename to src/fix_update_special_bonds.h index 600b226b71..dc6ca1cf9b 100644 --- a/src/BPM/fix_update_special_bonds.h +++ b/src/fix_update_special_bonds.h @@ -22,7 +22,7 @@ FixStyle(UPDATE_SPECIAL_BONDS,FixUpdateSpecialBonds) #include "fix.h" -#include +#include #include namespace LAMMPS_NS { @@ -34,16 +34,14 @@ class FixUpdateSpecialBonds : public Fix { int setmask(); void setup(int); void pre_exchange(); - void pre_force(int); - int pack_forward_comm(int, int *, double *, int, int *); - void unpack_forward_comm(int, int, double *); + void pre_force(int); void add_broken_bond(int,int); protected: - std::set > broken_pairs; - inline int sbmask(int j) const { - return j >> SBBITS & 3; - } + // Create two arrays to store bonds broken this timestep (new) + // and since the last neighbor list build + std::vector > new_broken_pairs; + std::vector > broken_pairs; }; } // namespace LAMMPS_NS diff --git a/src/lmptype.h b/src/lmptype.h index 905323ac2f..e02cd74734 100644 --- a/src/lmptype.h +++ b/src/lmptype.h @@ -55,14 +55,15 @@ namespace LAMMPS_NS { -// reserve 2 hi bits in molecular system neigh list for special bonds flag -// reserve 3rd last bit in neigh list for fix neigh/history flag +// reserve 2 highest bits in molecular system neigh list for special bonds flag +// reserve 3rd highest bit in neigh list for fix neigh/history flag // max local + ghost atoms per processor = 2^29 - 1 #define SBBITS 30 #define HISTBITS 29 #define NEIGHMASK 0x1FFFFFFF #define HISTMASK 0xDFFFFFFF +#define SPECIALMASK 0x3FFFFFFF // default to 32-bit smallint and other ints, 64-bit bigint