diff --git a/doc/src/Commands_compute.rst b/doc/src/Commands_compute.rst index ab4d8d64a2..0352ad5374 100644 --- a/doc/src/Commands_compute.rst +++ b/doc/src/Commands_compute.rst @@ -118,6 +118,7 @@ KOKKOS, o = OPENMP, t = OPT. * :doc:`ptm/atom ` * :doc:`rattlers/atom ` * :doc:`rdf ` + * :doc:`reaxff/atom (k) ` * :doc:`reduce ` * :doc:`reduce/chunk ` * :doc:`reduce/region ` diff --git a/doc/src/compute.rst b/doc/src/compute.rst index afe57d9048..7b620deed7 100644 --- a/doc/src/compute.rst +++ b/doc/src/compute.rst @@ -282,6 +282,7 @@ The individual style names on the :doc:`Commands compute ` pag * :doc:`ptm/atom ` - determines the local lattice structure based on the Polyhedral Template Matching method * :doc:`rattlers/atom ` - identify under-coordinated rattler atoms * :doc:`rdf ` - radial distribution function :math:`g(r)` histogram of group of atoms +* :doc:`reaxff/atom ` - extract ReaxFF bond information * :doc:`reduce ` - combine per-atom quantities into a single global value * :doc:`reduce/chunk ` - reduce per-atom quantities within each chunk * :doc:`reduce/region ` - same as compute reduce, within a region diff --git a/doc/src/compute_reaxff_atom.rst b/doc/src/compute_reaxff_atom.rst new file mode 100644 index 0000000000..997ad02e9f --- /dev/null +++ b/doc/src/compute_reaxff_atom.rst @@ -0,0 +1,97 @@ +.. index:: compute reaxff/atom +.. index:: compute reaxff/atom/kk + +compute reaxff/atom command +=========================== + +Accelerator Variants: *reaxff/atom/kk* + +Syntax +"""""" + +.. code-block:: LAMMPS + + compute ID group-ID reaxff/atom attribute args ... keyword value ... + +* ID, group-ID are documented in :doc:`compute ` command +* reaxff/atom = name of this compute command +* attribute = *pair* + + .. parsed-literal:: + + *pair* args = nsub + nsub = *n*-instance of a sub-style, if a pair style is used multiple times in a hybrid style + +* keyword = *bonds* + + .. parsed-literal:: + + *bonds* value = *no* or *yes* + *no* = ignore list of local bonds + *yes* = include list of local bonds + +Examples +"""""""" + +.. code-block:: LAMMPS + + compute 1 all reaxff/atom bonds yes + +Description +""""""""""" + +.. versionadded:: TBD + +Define a computation that extracts bond information computed by the ReaxFF +potential specified by :doc:`pair_style reaxff `. + +By default, it produces per-atom data that includes the following columns: + +* abo = atom bond order (sum of all bonds) +* nlp = number of lone pairs +* nb = number of bonds + +Bonds will only be included if its atoms are in the group. + +In addition, if ``bonds`` is set to ``yes``, the compute will also produce a +local array of all bonds on the current processor whose atoms are in the group. +The columns of each entry of this local array are: + +* id_i = atom i id of bond +* id_j = atom j id of bond +* bo = bond order of bond + +Output info +""""""""""" + +This compute calculates a per-atom array and local array depending on the +number of keywords. The number of rows in the local array is the number of +bonds as described above. Both per-atom and local array have 3 columns. + +The arrays can be accessed by any command that uses local and per-atom values +from a compute as input. See the :doc:`Howto output ` page for +an overview of LAMMPS output options. + +---------- + +.. include:: accel_styles.rst + +---------- + +Restrictions +"""""""""""" + +The compute reaxff/atom command requires that the :doc:`pair_style reaxff +` is invoked. This fix is part of the REAXFF package. It is only +enabled if LAMMPS was built with that package. See the :doc:`Build package +` page for more info. + +Related commands +"""""""""""""""" + +:doc:`pair_style reaxff ` + +Default +""""""" + +The option defaults are *bonds* = *no*. diff --git a/doc/src/compute_stress_mop.rst b/doc/src/compute_stress_mop.rst index 21c2963545..74d4c618e7 100644 --- a/doc/src/compute_stress_mop.rst +++ b/doc/src/compute_stress_mop.rst @@ -18,7 +18,7 @@ Syntax * style = *stress/mop* or *stress/mop/profile* * dir = *x* or *y* or *z* is the direction normal to the plane * args = argument specific to the compute style -* keywords = *kin* or *conf* or *total* or *pair* or *bond* or *angle* (one or more can be specified) +* keywords = *kin* or *conf* or *total* or *pair* or *bond* or *angle* or *dihedral* (one or more can be specified) .. parsed-literal:: @@ -68,15 +68,13 @@ Verlet algorithm. .. versionadded:: 15Jun2023 - contributions from bond and angle potentials + contributions from bond, angle and dihedral potentials -Between one and six keywords can be used to indicate which contributions +Between one and seven keywords can be used to indicate which contributions to the stress must be computed: total stress (total), kinetic stress (kin), configurational stress (conf), stress due to bond stretching -(bond), stress due to angle bending (angle) and/or due to pairwise -non-bonded interactions (pair). The angle keyword is currently -available only for the *stress/mop* command and **not** the -*stress/mop/profile* command. +(bond), stress due to angle bending (angle), stress due to dihedral terms (dihedral) +and/or due to pairwise non-bonded interactions (pair). NOTE 1: The configurational stress is computed considering all pairs of atoms where at least one atom belongs to group group-ID. @@ -134,14 +132,9 @@ size does not change in time, and axis-aligned planes. The method only works with two-body pair interactions, because it requires the class method ``Pair::single()`` to be implemented, which is not possible for manybody potentials. In particular, compute -*stress/mop/profile* does not work with more than two-body pair +*stress/mop/profile* and *stress/mop* do not work with more than two-body pair interactions, long range (kspace) interactions and -angle/dihedral/improper intramolecular interactions. Similarly, compute -*stress/mop* does not work with more than two-body pair interactions, -long range (kspace) interactions and dihedral/improper intramolecular -interactions but works with all bond interactions with the class method -single() implemented and all angle interactions with the class method -born_matrix() implemented. +improper intramolecular interactions. Related commands """""""""""""""" diff --git a/doc/src/pair_reaxff.rst b/doc/src/pair_reaxff.rst index d28e15b0a2..03d53d1ff4 100644 --- a/doc/src/pair_reaxff.rst +++ b/doc/src/pair_reaxff.rst @@ -373,7 +373,8 @@ Related commands :doc:`pair_coeff `, :doc:`fix qeq/reaxff `, :doc:`fix acks2/reaxff `, :doc:`fix reaxff/bonds `, -:doc:`fix reaxff/species ` +:doc:`fix reaxff/species `, +:doc:`compute reaxff/atom ` Default """"""" diff --git a/examples/reaxff/in.reaxff.tatb b/examples/reaxff/in.reaxff.tatb index 6cf7828cf1..967ed0a1d6 100644 --- a/examples/reaxff/in.reaxff.tatb +++ b/examples/reaxff/in.reaxff.tatb @@ -31,8 +31,15 @@ neigh_modify delay 0 every 5 check no fix 1 all nve fix 2 all qeq/reaxff 1 0.0 10.0 1.0e-6 reaxff fix 4 all reaxff/bonds 5 bonds.reaxff +compute bonds all reaxff/atom bonds yes variable nqeq equal f_2 +# dumps out the local bond information +dump 1 all local 5 bonds_local.reaxff c_bonds[1] c_bonds[2] c_bonds[3] + +# dumps out the peratom bond information +dump 2 all custom 5 bonds_atom.reaxff id type q c_bonds[*] + thermo 5 thermo_style custom step temp epair etotal press & v_eb v_ea v_elp v_emol v_ev v_epen v_ecoa & diff --git a/src/.gitignore b/src/.gitignore index 3b7902303d..112a1486f7 100644 --- a/src/.gitignore +++ b/src/.gitignore @@ -633,6 +633,8 @@ /compute_ptm_atom.h /compute_rattlers_atom.cpp /compute_rattlers_atom.h +/compute_reaxff_atom.cpp +/compute_reaxff_atom.h /compute_rigid_local.cpp /compute_rigid_local.h /compute_slcsa_atom.cpp diff --git a/src/EXTRA-COMPUTE/compute_stress_mop.cpp b/src/EXTRA-COMPUTE/compute_stress_mop.cpp index fc9de602a7..6c35b4ba07 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop.cpp @@ -23,6 +23,7 @@ #include "atom_vec.h" #include "bond.h" #include "comm.h" +#include "dihedral.h" #include "domain.h" #include "error.h" #include "force.h" @@ -38,8 +39,10 @@ using namespace LAMMPS_NS; +#define SMALL 0.001 + enum { X, Y, Z }; -enum { TOTAL, CONF, KIN, PAIR, BOND, ANGLE }; +enum { TOTAL, CONF, KIN, PAIR, BOND, ANGLE, DIHEDRAL }; /* ---------------------------------------------------------------------- */ @@ -49,6 +52,7 @@ ComputeStressMop::ComputeStressMop(LAMMPS *lmp, int narg, char **arg) : Compute( bondflag = 0; angleflag = 0; + dihedralflag = 0; // set compute mode and direction of plane(s) for pressure calculation @@ -129,6 +133,11 @@ ComputeStressMop::ComputeStressMop(LAMMPS *lmp, int narg, char **arg) : Compute( which[nvalues] = ANGLE; nvalues++; } + } else if (strcmp(arg[iarg],"dihedral") == 0) { + for (i=0; i<3; i++) { + which[nvalues] = DIHEDRAL; + nvalues++; + } } else error->all(FLERR, "Illegal compute stress/mop command"); //break; @@ -152,6 +161,8 @@ ComputeStressMop::ComputeStressMop(LAMMPS *lmp, int narg, char **arg) : Compute( bond_global = nullptr; angle_local = nullptr; angle_global = nullptr; + dihedral_local = nullptr; + dihedral_global = nullptr; // this fix produces a global vector @@ -162,6 +173,8 @@ ComputeStressMop::ComputeStressMop(LAMMPS *lmp, int narg, char **arg) : Compute( memory->create(bond_global, nvalues, "stress/mop:bond_global"); memory->create(angle_local, nvalues, "stress/mop:angle_local"); memory->create(angle_global, nvalues, "stress/mop:angle_global"); + memory->create(dihedral_local,nvalues,"stress/mop:dihedral_local"); + memory->create(dihedral_global,nvalues,"stress/mop:dihedral_global"); size_vector = nvalues; vector_flag = 1; @@ -180,6 +193,8 @@ ComputeStressMop::~ComputeStressMop() memory->destroy(bond_global); memory->destroy(angle_local); memory->destroy(angle_global); + memory->destroy(dihedral_local); + memory->destroy(dihedral_global); memory->destroy(vector); } @@ -233,9 +248,13 @@ void ComputeStressMop::init() } } if (force->dihedral) { - if ((strcmp(force->dihedral_style, "zero") != 0) && - (strcmp(force->dihedral_style, "none") != 0)) - error->all(FLERR, "compute stress/mop does not account for dihedral potentials"); + if (force->dihedral->born_matrix_enable == 0) { + if ((strcmp(force->dihedral_style, "zero") != 0) && + (strcmp(force->dihedral_style, "none") != 0)) + error->all(FLERR, "compute stress/mop does not account for dihedral potentials"); + } else { + dihedralflag = 1; + } } if (force->improper) { if ((strcmp(force->improper_style, "zero") != 0) && @@ -297,8 +316,18 @@ void ComputeStressMop::compute_vector() MPI_Allreduce(angle_local, angle_global, nvalues, MPI_DOUBLE, MPI_SUM, world); + if (dihedralflag) { + //Compute dihedral contribution on separate procs + compute_dihedrals(); + } else { + for (int i=0; ix[i][1]; xi[2] = atom->x[i][2]; - // velocities at t + // minimum image of xi with respect to the plane + xi[dir] -= pos; + domain->minimum_image(xi[0], xi[1], xi[2]); + xi[dir] += pos; + + //velocities at t vi[0] = atom->v[i][0]; vi[1] = atom->v[i][1]; @@ -454,10 +488,8 @@ void ComputeStressMop::compute_pairs() // at each timestep, must check atoms going through the // image of the plane that is closest to the box - double pos_temp = pos + copysign(1.0, domain->prd_half[dir] - pos) * domain->prd[dir]; - if (fabs(xi[dir] - pos) < fabs(xi[dir] - pos_temp)) pos_temp = pos; - - if (((xi[dir] - pos_temp) * (xj[dir] - pos_temp)) < 0) { + double tau = (xi[dir] - pos) / (xi[dir] - xj[dir]); + if ((tau <= 1) && (tau >= 0)) { // sgn = copysign(1.0,vi[dir]-vcm[dir]); @@ -786,3 +818,308 @@ void ComputeStressMop::compute_angles() m += 3; } } + +/*------------------------------------------------------------------------ + compute dihedral contribution to pressure of local proc + -------------------------------------------------------------------------*/ + +void ComputeStressMop::compute_dihedrals() +{ + int i, nd, atom1, atom2, atom3, atom4, imol, iatom; + tagint tagprev; + double vb1x, vb1y, vb1z, vb2x, vb2y, vb2z, vb3x, vb3y, vb3z; + double vb2xm, vb2ym, vb2zm; + double sb1, sb2, sb3, rb1, rb3, c0, b1mag2, b1mag, b2mag2; + double b2mag, b3mag2, b3mag, c2mag, ctmp, r12c1, c1mag, r12c2; + double s1, s2, s12, sc1, sc2, a11, a22, a33, a12, a13, a23; + double df[3], f1[3], f2[3], f3[3], f4[3]; + double c, sx2, sy2, sz2, sin2; + + double **x = atom->x; + tagint *tag = atom->tag; + int *num_dihedral = atom->num_dihedral; + tagint **dihedral_atom1 = atom->dihedral_atom1; + tagint **dihedral_atom2 = atom->dihedral_atom2; + tagint **dihedral_atom3 = atom->dihedral_atom3; + tagint **dihedral_atom4 = atom->dihedral_atom4; + int *mask = atom->mask; + + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; + + int nlocal = atom->nlocal; + int molecular = atom->molecular; + + // loop over all atoms and their dihedrals + + Dihedral *dihedral = force->dihedral; + + double dudih, du2dih; + + double diffx[3] = {0.0, 0.0, 0.0}; + double x_atom_1[3] = {0.0, 0.0, 0.0}; + double x_atom_2[3] = {0.0, 0.0, 0.0}; + double x_atom_3[3] = {0.0, 0.0, 0.0}; + double x_atom_4[3] = {0.0, 0.0, 0.0}; + + // initialization + for (int i = 0; i < nvalues; i++) { + dihedral_local[i] = 0.0; + } + double local_contribution[3] = {0.0, 0.0, 0.0}; + + for (atom2 = 0; atom2 < nlocal; atom2++) { + if (!(mask[atom2] & groupbit)) continue; + + if (molecular == Atom::MOLECULAR) + nd = num_dihedral[atom2]; + else { + if (molindex[atom2] < 0) continue; + imol = molindex[atom2]; + iatom = molatom[atom2]; + nd = onemols[imol]->num_dihedral[iatom]; + } + + for (i = 0; i < nd; i++) { + if (molecular == 1) { + if (tag[atom2] != dihedral_atom2[atom2][i]) continue; + atom1 = atom->map(dihedral_atom1[atom2][i]); + atom3 = atom->map(dihedral_atom3[atom2][i]); + atom4 = atom->map(dihedral_atom4[atom2][i]); + } else { + if (tag[atom2] != onemols[imol]->dihedral_atom2[atom2][i]) continue; + tagprev = tag[atom2] - iatom - 1; + atom1 = atom->map(onemols[imol]->dihedral_atom1[atom2][i] + tagprev); + atom3 = atom->map(onemols[imol]->dihedral_atom3[atom2][i] + tagprev); + atom4 = atom->map(onemols[imol]->dihedral_atom4[atom2][i] + tagprev); + } + + if (atom1 < 0 || !(mask[atom1] & groupbit)) continue; + if (atom3 < 0 || !(mask[atom3] & groupbit)) continue; + if (atom4 < 0 || !(mask[atom4] & groupbit)) continue; + + // minimum image of atom1 with respect to the plane of interest + x_atom_1[0] = x[atom1][0]; + x_atom_1[1] = x[atom1][1]; + x_atom_1[2] = x[atom1][2]; + x_atom_1[dir] -= pos; + domain->minimum_image(x_atom_1[0], x_atom_1[1], x_atom_1[2]); + x_atom_1[dir] += pos; + + // minimum image of atom2 with respect to atom1 + diffx[0] = x[atom2][0] - x_atom_1[0]; + diffx[1] = x[atom2][1] - x_atom_1[1]; + diffx[2] = x[atom2][2] - x_atom_1[2]; + domain->minimum_image(diffx[0], diffx[1], diffx[2]); + x_atom_2[0] = x_atom_1[0] + diffx[0]; + x_atom_2[1] = x_atom_1[1] + diffx[1]; + x_atom_2[2] = x_atom_1[2] + diffx[2]; + + // minimum image of atom3 with respect to atom2 + diffx[0] = x[atom3][0] - x_atom_2[0]; + diffx[1] = x[atom3][1] - x_atom_2[1]; + diffx[2] = x[atom3][2] - x_atom_2[2]; + domain->minimum_image(diffx[0], diffx[1], diffx[2]); + x_atom_3[0] = x_atom_2[0] + diffx[0]; + x_atom_3[1] = x_atom_2[1] + diffx[1]; + x_atom_3[2] = x_atom_2[2] + diffx[2]; + + // minimum image of atom3 with respect to atom2 + diffx[0] = x[atom4][0] - x_atom_3[0]; + diffx[1] = x[atom4][1] - x_atom_3[1]; + diffx[2] = x[atom4][2] - x_atom_3[2]; + domain->minimum_image(diffx[0], diffx[1], diffx[2]); + x_atom_4[0] = x_atom_3[0] + diffx[0]; + x_atom_4[1] = x_atom_3[1] + diffx[1]; + x_atom_4[2] = x_atom_3[2] + diffx[2]; + + // check if any bond vector crosses the plane of interest + double tau_right = (x_atom_2[dir] - pos) / (x_atom_2[dir] - x_atom_1[dir]); + double tau_middle = (x_atom_3[dir] - pos) / (x_atom_3[dir] - x_atom_2[dir]); + double tau_left = (x_atom_4[dir] - pos) / (x_atom_4[dir] - x_atom_3[dir]); + bool right_cross = ((tau_right >=0) && (tau_right <= 1)); + bool middle_cross = ((tau_middle >= 0) && (tau_middle <= 1)); + bool left_cross = ((tau_left >=0) && (tau_left <= 1)); + + // no bonds crossing the plane + if (!right_cross && !middle_cross && !left_cross) continue; + + dihedral->born_matrix(i, atom1, atom2, atom3, atom4, dudih, du2dih); + + // first bond + vb1x = x_atom_1[0] - x_atom_2[0]; + vb1y = x_atom_1[1] - x_atom_2[1]; + vb1z = x_atom_1[2] - x_atom_2[2]; + + // second bond + vb2x = x_atom_3[0] - x_atom_2[0]; + vb2y = x_atom_3[1] - x_atom_2[1]; + vb2z = x_atom_3[2] - x_atom_2[2]; + + vb2xm = -vb2x; + vb2ym = -vb2y; + vb2zm = -vb2z; + + // third bond + vb3x = x_atom_4[0] - x_atom_3[0]; + vb3y = x_atom_4[1] - x_atom_3[1]; + vb3z = x_atom_4[2] - x_atom_3[2]; + + // c0 calculation + sb1 = 1.0 / (vb1x*vb1x + vb1y*vb1y + vb1z*vb1z); + sb2 = 1.0 / (vb2x*vb2x + vb2y*vb2y + vb2z*vb2z); + sb3 = 1.0 / (vb3x*vb3x + vb3y*vb3y + vb3z*vb3z); + + rb1 = sqrt(sb1); + rb3 = sqrt(sb3); + + c0 = (vb1x*vb3x + vb1y*vb3y + vb1z*vb3z) * rb1*rb3; + // 1st and 2nd angle + b1mag2 = vb1x*vb1x + vb1y*vb1y + vb1z*vb1z; + b1mag = sqrt(b1mag2); + b2mag2 = vb2x*vb2x + vb2y*vb2y + vb2z*vb2z; + b2mag = sqrt(b2mag2); + b3mag2 = vb3x*vb3x + vb3y*vb3y + vb3z*vb3z; + b3mag = sqrt(b3mag2); + + ctmp = vb1x*vb2x + vb1y*vb2y + vb1z*vb2z; + r12c1 = 1.0 / (b1mag*b2mag); + c1mag = ctmp * r12c1; + + ctmp = vb2xm*vb3x + vb2ym*vb3y + vb2zm*vb3z; + r12c2 = 1.0 / (b2mag*b3mag); + c2mag = ctmp * r12c2; + + // cos and sin of 2 angles and final c + sin2 = MAX(1.0 - c1mag*c1mag,0.0); + sc1 = sqrt(sin2); + if (sc1 < SMALL) sc1 = SMALL; + sc1 = 1.0/sc1; + + sin2 = MAX(1.0 - c2mag*c2mag,0.0); + sc2 = sqrt(sin2); + if (sc2 < SMALL) sc2 = SMALL; + sc2 = 1.0/sc2; + + s1 = sc1 * sc1; + s2 = sc2 * sc2; + s12 = sc1 * sc2; + c = (c0 + c1mag*c2mag) * s12; + + // error check + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + + // forces on each particle + double a = dudih; + c = c * a; + s12 = s12 * a; + a11 = c*sb1*s1; + a22 = -sb2 * (2.0*c0*s12 - c*(s1+s2)); + a33 = c*sb3*s2; + a12 = -r12c1 * (c1mag*c*s1 + c2mag*s12); + a13 = -rb1*rb3*s12; + a23 = r12c2 * (c2mag*c*s2 + c1mag*s12); + + sx2 = a12*vb1x + a22*vb2x + a23*vb3x; + sy2 = a12*vb1y + a22*vb2y + a23*vb3y; + sz2 = a12*vb1z + a22*vb2z + a23*vb3z; + + f1[0] = a11*vb1x + a12*vb2x + a13*vb3x; + f1[1] = a11*vb1y + a12*vb2y + a13*vb3y; + f1[2] = a11*vb1z + a12*vb2z + a13*vb3z; + + f2[0] = -sx2 - f1[0]; + f2[1] = -sy2 - f1[1]; + f2[2] = -sz2 - f1[2]; + + f4[0] = a13*vb1x + a23*vb2x + a33*vb3x; + f4[1] = a13*vb1y + a23*vb2y + a33*vb3y; + f4[2] = a13*vb1z + a23*vb2z + a33*vb3z; + + f3[0] = sx2 - f4[0]; + f3[1] = sy2 - f4[1]; + f3[2] = sz2 - f4[2]; + + // only right bond crossing the plane + if (right_cross && !middle_cross && !left_cross) + { + double sgn = copysign(1.0, x_atom_1[dir] - pos); + df[0] = sgn * f1[0]; + df[1] = sgn * f1[1]; + df[2] = sgn * f1[2]; + } + + // only middle bond crossing the plane + if (!right_cross && middle_cross && !left_cross) + { + double sgn = copysign(1.0, x_atom_2[dir] - pos); + df[0] = sgn * (f2[0] + f1[0]); + df[1] = sgn * (f2[1] + f1[1]); + df[2] = sgn * (f2[2] + f1[2]); + } + + // only left bond crossing the plane + if (!right_cross && !middle_cross && left_cross) + { + double sgn = copysign(1.0, x_atom_4[dir] - pos); + df[0] = sgn * f4[0]; + df[1] = sgn * f4[1]; + df[2] = sgn * f4[2]; + } + + // only right & middle bonds crossing the plane + if (right_cross && middle_cross && !left_cross) + { + double sgn = copysign(1.0, x_atom_2[dir] - pos); + df[0] = sgn * f2[0]; + df[1] = sgn * f2[1]; + df[2] = sgn * f2[2]; + } + + // only right & left bonds crossing the plane + if (right_cross && !middle_cross && left_cross) + { + double sgn = copysign(1.0, x_atom_1[dir] - pos); + df[0] = sgn * (f1[0] + f4[0]); + df[1] = sgn * (f1[1] + f4[1]); + df[2] = sgn * (f1[2] + f4[2]); + } + + // only middle & left bonds crossing the plane + if (!right_cross && middle_cross && left_cross) + { + double sgn = copysign(1.0, x_atom_3[dir] - pos); + df[0] = sgn * f3[0]; + df[1] = sgn * f3[1]; + df[2] = sgn * f3[2]; + } + + // all three bonds crossing the plane + if (right_cross && middle_cross && left_cross) + { + double sgn = copysign(1.0, x_atom_1[dir] - pos); + df[0] = sgn * (f1[0] + f3[0]); + df[1] = sgn * (f1[1] + f3[1]); + df[2] = sgn * (f1[2] + f3[2]); + } + + local_contribution[0] += df[0]/area*nktv2p; + local_contribution[1] += df[1]/area*nktv2p; + local_contribution[2] += df[2]/area*nktv2p; + } + } + + // loop over the keywords and if necessary add the dihedral contribution + int m = 0; + while (m < nvalues) { + if ((which[m] == CONF) || (which[m] == TOTAL) || (which[m] == DIHEDRAL)) { + dihedral_local[m] = local_contribution[0]; + dihedral_local[m+1] = local_contribution[1]; + dihedral_local[m+2] = local_contribution[2]; + } + m += 3; + } + +} diff --git a/src/EXTRA-COMPUTE/compute_stress_mop.h b/src/EXTRA-COMPUTE/compute_stress_mop.h index 86140dc278..0a0ea8b55a 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop.h +++ b/src/EXTRA-COMPUTE/compute_stress_mop.h @@ -40,15 +40,17 @@ class ComputeStressMop : public Compute { void compute_pairs(); void compute_bonds(); void compute_angles(); + void compute_dihedrals(); int nvalues, dir; int *which; - int bondflag, angleflag; + int bondflag, angleflag, dihedralflag; double *values_local, *values_global; double *bond_local, *bond_global; double *angle_local, *angle_global; + double *dihedral_local, *dihedral_global; double pos, pos1, dt, nktv2p, ftm2v; double area; class NeighList *list; diff --git a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp index cc201fdbaa..41b5f64a67 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp @@ -13,15 +13,17 @@ /*------------------------------------------------------------------------ Contributing Authors : Romain Vermorel (LFCR), Laurent Joly (ULyon) - Support for bonds added by : Evangelos Voyiatzis (NovaMechanics) + Support for bonds, angles and dihedrals added by : Evangelos Voyiatzis (NovaMechanics) --------------------------------------------------------------------------*/ #include "compute_stress_mop_profile.h" +#include "angle.h" #include "atom.h" #include "atom_vec.h" #include "bond.h" #include "comm.h" +#include "dihedral.h" #include "domain.h" #include "error.h" #include "force.h" @@ -37,9 +39,10 @@ using namespace LAMMPS_NS; +#define SMALL 0.001 + enum { X, Y, Z }; -enum { LOWER, CENTER, UPPER, COORD }; -enum { TOTAL, CONF, KIN, PAIR, BOND }; +enum { TOTAL, CONF, KIN, PAIR, BOND, ANGLE, DIHEDRAL }; /* ---------------------------------------------------------------------- */ @@ -49,6 +52,8 @@ ComputeStressMopProfile::ComputeStressMopProfile(LAMMPS *lmp, int narg, char **a if (narg < 7) utils::missing_cmd_args(FLERR, "compute stress/mop/profile", error); bondflag = 0; + angleflag = 0; + dihedralflag = 0; // set compute mode and direction of plane(s) for pressure calculation @@ -63,15 +68,15 @@ ComputeStressMopProfile::ComputeStressMopProfile(LAMMPS *lmp, int narg, char **a // bin parameters - if (strcmp(arg[4], "lower") == 0) - originflag = LOWER; - else if (strcmp(arg[4], "center") == 0) - originflag = CENTER; - else if (strcmp(arg[4], "upper") == 0) - originflag = UPPER; - else - originflag = COORD; - if (originflag == COORD) origin = utils::numeric(FLERR, arg[4], false, lmp); + if (strcmp(arg[4], "lower") == 0) { + origin = domain->boxlo[dir]; + } else if (strcmp(arg[4], "center") == 0) { + origin = 0.5 * (domain->boxlo[dir] + domain->boxhi[dir]); + } else if (strcmp(arg[4], "upper") == 0) { + origin = domain->boxhi[dir]; + } else { + origin = utils::numeric(FLERR, arg[4], false, lmp); + } delta = utils::numeric(FLERR, arg[5], false, lmp); invdelta = 1.0 / delta; @@ -108,6 +113,16 @@ ComputeStressMopProfile::ComputeStressMopProfile(LAMMPS *lmp, int narg, char **a which[nvalues] = BOND; nvalues++; } + } else if (strcmp(arg[iarg], "angle") == 0) { + for (i = 0; i < 3; i++) { + which[nvalues] = ANGLE; + nvalues++; + } + } else if (strcmp(arg[iarg],"dihedral") == 0) { + for (i=0; i<3; i++) { + which[nvalues] = DIHEDRAL; + nvalues++; + } } else error->all(FLERR, "Illegal compute stress/mop/profile command"); //break; @@ -133,6 +148,10 @@ ComputeStressMopProfile::ComputeStressMopProfile(LAMMPS *lmp, int narg, char **a values_local = values_global = array = nullptr; bond_local = nullptr; bond_global = nullptr; + angle_local = nullptr; + angle_global = nullptr; + dihedral_local = nullptr; + dihedral_global = nullptr; local_contribution = nullptr; // bin setup @@ -161,6 +180,10 @@ ComputeStressMopProfile::~ComputeStressMopProfile() memory->destroy(values_global); memory->destroy(bond_local); memory->destroy(bond_global); + memory->destroy(angle_local); + memory->destroy(angle_global); + memory->destroy(dihedral_local); + memory->destroy(dihedral_global); memory->destroy(local_contribution); memory->destroy(array); } @@ -208,13 +231,25 @@ void ComputeStressMopProfile::init() if (force->bond) bondflag = 1; - if (force->angle) - if ((strcmp(force->angle_style, "zero") != 0) && (strcmp(force->angle_style, "none") != 0)) - error->all(FLERR, "compute stress/mop/profile does not account for angle potentials"); - if (force->dihedral) - if ((strcmp(force->dihedral_style, "zero") != 0) && - (strcmp(force->dihedral_style, "none") != 0)) - error->all(FLERR, "compute stress/mop/profile does not account for dihedral potentials"); + if (force->angle) { + if (force->angle->born_matrix_enable == 0) { + if ((strcmp(force->angle_style, "zero") != 0) && (strcmp(force->angle_style, "none") != 0)) + error->all(FLERR,"compute stress/mop/profile does not account for angle potentials"); + } else { + angleflag = 1; + } + } + + if (force->dihedral) { + if (force->dihedral->born_matrix_enable == 0) { + if ((strcmp(force->dihedral_style, "zero") != 0) && + (strcmp(force->dihedral_style, "none") != 0)) + error->all(FLERR, "compute stress/mop/profile does not account for dihedral potentials"); + } else { + dihedralflag = 1; + } + } + if (force->improper) if ((strcmp(force->improper_style, "zero") != 0) && (strcmp(force->improper_style, "none") != 0)) @@ -263,16 +298,43 @@ void ComputeStressMopProfile::compute_array() } // sum bond contribution over all procs - MPI_Allreduce(&bond_local[0][0], &bond_global[0][0], nbins * nvalues, MPI_DOUBLE, MPI_SUM, world); + if (angleflag) { + //Compute angle contribution on separate procs + compute_angles(); + } else { + for (int m = 0; m < nbins; m++) { + for (int i = 0; i < nvalues; i++) { + angle_local[m][i] = 0.0; + } + } + } + + // sum angle contribution over all procs + MPI_Allreduce(&angle_local[0][0],&angle_global[0][0],nbins*nvalues,MPI_DOUBLE,MPI_SUM,world); + + if (dihedralflag) { + //Compute dihedral contribution on separate procs + compute_dihedrals(); + } else { + for (int m = 0; m < nbins; m++) { + for (int i = 0; i < nvalues; i++) { + dihedral_local[m][i] = 0.0; + } + } + } + + // sum dihedral contribution over all procs + MPI_Allreduce(&dihedral_local[0][0],&dihedral_global[0][0],nbins*nvalues,MPI_DOUBLE,MPI_SUM,world); + for (int ibin = 0; ibin < nbins; ibin++) { - array[ibin][0] = coord[ibin][0]; + array[ibin][0] = coord[ibin]; int mo = 1; int m = 0; while (m < nvalues) { - array[ibin][m + mo] = values_global[ibin][m] + bond_global[ibin][m]; + array[ibin][m + mo] = values_global[ibin][m] + bond_global[ibin][m] + angle_global[ibin][m] + dihedral_global[ibin][m]; m++; } } @@ -366,8 +428,8 @@ void ComputeStressMopProfile::compute_pairs() if (newton_pair || j < nlocal) { for (ibin = 0; ibin < nbins; ibin++) { - pos = coord[ibin][0]; - pos1 = coordp[ibin][0]; + pos = coord[ibin]; + pos1 = coordp[ibin]; // check if ij pair is across plane, add contribution to pressure @@ -392,8 +454,8 @@ void ComputeStressMopProfile::compute_pairs() } else { for (ibin = 0; ibin < nbins; ibin++) { - pos = coord[ibin][0]; - pos1 = coordp[ibin][0]; + pos = coord[ibin]; + pos1 = coordp[ibin]; //check if ij pair is across plane, add contribution to pressure @@ -454,15 +516,29 @@ void ComputeStressMopProfile::compute_pairs() xj[2] = xi[2] - vi[2] * dt + fi[2] * iterm * dt; for (ibin = 0; ibin < nbins; ibin++) { - pos = coord[ibin][0]; - pos1 = coordp[ibin][0]; + pos = coord[ibin]; + pos1 = coordp[ibin]; - if (((xi[dir] - pos) * (xj[dir] - pos) * (xi[dir] - pos1) * (xj[dir] - pos1) < 0)) { + // minimum image of xi with respect to the plane + xi[dir] -= pos; + domain->minimum_image(xi[0], xi[1], xi[2]); + xi[dir] += pos; + + // minimum image of xj with respect to xi + xj[0] -= xi[0]; + xj[1] -= xi[1]; + xj[2] -= xi[2]; + domain->minimum_image(xi[0], xi[1], xi[2]); + xj[0] += xi[0]; + xj[1] += xi[1]; + xj[2] += xi[2]; + + double tau = (xi[dir] - pos) / (xi[dir] - xj[dir]); + if ((tau <= 1) && (tau >= 0)) { sgn = copysign(1.0, vi[dir]); - // approximate crossing velocity by v(t-dt/2) (based on Velocity-Verlet alg.) - + //approximate crossing velocity by v(t-dt/2) (based on Velocity-Verlet alg.) double vcross[3]; vcross[0] = vi[0] - fi[0] * iterm; vcross[1] = vi[1] - fi[1] * iterm; @@ -549,7 +625,7 @@ void ComputeStressMopProfile::compute_bonds() if (btype <= 0) continue; for (int ibin = 0; ibin < nbins; ibin++) { - double pos = coord[ibin][0]; + double pos = coord[ibin]; // minimum image of atom1 with respect to the plane of interest @@ -607,6 +683,506 @@ void ComputeStressMopProfile::compute_bonds() } } +/*------------------------------------------------------------------------ + compute angle contribution to pressure of local proc + -------------------------------------------------------------------------*/ + +void ComputeStressMopProfile::compute_angles() +{ + int na, atom1, atom2, atom3, imol, iatom, atype; + tagint tagprev; + double r1, r2, cos_theta; + + double **x = atom->x; + tagint *tag = atom->tag; + int *num_angle = atom->num_angle; + tagint **angle_atom1 = atom->angle_atom1; + tagint **angle_atom2 = atom->angle_atom2; + tagint **angle_atom3 = atom->angle_atom3; + int **angle_type = atom->angle_type; + int *mask = atom->mask; + + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; + + int nlocal = atom->nlocal; + int molecular = atom->molecular; + + // loop over all atoms and their angles + Angle *angle = force->angle; + + double duang, du2ang; + double dx[3] = {0.0, 0.0, 0.0}; + double dx_left[3] = {0.0, 0.0, 0.0}; + double dx_right[3] = {0.0, 0.0, 0.0}; + double x_angle_left[3] = {0.0, 0.0, 0.0}; + double x_angle_middle[3] = {0.0, 0.0, 0.0}; + double x_angle_right[3] = {0.0, 0.0, 0.0}; + double dcos_theta[3] = {0.0, 0.0, 0.0}; + + // initialization + for (int m = 0; m < nbins; m++) { + for (int i = 0; i < nvalues; i++) { + angle_local[m][i] = 0.0; + } + local_contribution[m][0] = 0.0; + local_contribution[m][1] = 0.0; + local_contribution[m][2] = 0.0; + } + + + for (atom2 = 0; atom2 < nlocal; atom2++) { + if (!(mask[atom2] & groupbit)) continue; + + if (molecular == 1) + na = num_angle[atom2]; + else { + if (molindex[atom2] < 0) continue; + imol = molindex[atom2]; + iatom = molatom[atom2]; + na = onemols[imol]->num_angle[iatom]; + } + + for (int i = 0; i < na; i++) { + if (molecular == 1) { + if (tag[atom2] != angle_atom2[atom2][i]) continue; + atype = angle_type[atom2][i]; + atom1 = atom->map(angle_atom1[atom2][i]); + atom3 = atom->map(angle_atom3[atom2][i]); + } else { + if (tag[atom2] != onemols[imol]->angle_atom2[atom2][i]) continue; + atype = onemols[imol]->angle_type[atom2][i]; + tagprev = tag[atom2] - iatom - 1; + atom1 = atom->map(onemols[imol]->angle_atom1[atom2][i] + tagprev); + atom3 = atom->map(onemols[imol]->angle_atom3[atom2][i] + tagprev); + } + + if (atom1 < 0 || !(mask[atom1] & groupbit)) continue; + if (atom3 < 0 || !(mask[atom3] & groupbit)) continue; + if (atype <= 0) continue; + + for (int ibin = 0; ibinminimum_image(dx[0], dx[1], dx[2]); + x_angle_left[0] = dx[0]; + x_angle_left[1] = dx[1]; + x_angle_left[2] = dx[2]; + x_angle_left[dir] += pos; + + // minimum image of atom2 with respect to atom1 + dx_left[0] = x[atom2][0] - x_angle_left[0]; + dx_left[1] = x[atom2][1] - x_angle_left[1]; + dx_left[2] = x[atom2][2] - x_angle_left[2]; + domain->minimum_image(dx_left[0], dx_left[1], dx_left[2]); + x_angle_middle[0] = x_angle_left[0] + dx_left[0]; + x_angle_middle[1] = x_angle_left[1] + dx_left[1]; + x_angle_middle[2] = x_angle_left[2] + dx_left[2]; + + // minimum image of atom3 with respect to atom2 + dx_right[0] = x[atom3][0] - x_angle_middle[0]; + dx_right[1] = x[atom3][1] - x_angle_middle[1]; + dx_right[2] = x[atom3][2] - x_angle_middle[2]; + domain->minimum_image(dx_right[0], dx_right[1], dx_right[2]); + x_angle_right[0] = x_angle_middle[0] + dx_right[0]; + x_angle_right[1] = x_angle_middle[1] + dx_right[1]; + x_angle_right[2] = x_angle_middle[2] + dx_right[2]; + + // check if any bond vector crosses the plane of interest + double tau_right = (x_angle_right[dir] - pos) / (x_angle_right[dir] - x_angle_middle[dir]); + double tau_left = (x_angle_middle[dir] - pos) / (x_angle_middle[dir] - x_angle_left[dir]); + bool right_cross = ((tau_right >=0) && (tau_right <= 1)); + bool left_cross = ((tau_left >=0) && (tau_left <= 1)); + + // no bonds crossing the plane + if (!right_cross && !left_cross) continue; + + // compute the cos(theta) of the angle + r1 = sqrt(dx_left[0]*dx_left[0] + dx_left[1]*dx_left[1] + dx_left[2]*dx_left[2]); + r2 = sqrt(dx_right[0]*dx_right[0] + dx_right[1]*dx_right[1] + dx_right[2]*dx_right[2]); + cos_theta = -(dx_right[0]*dx_left[0] + dx_right[1]*dx_left[1] + dx_right[2]*dx_left[2])/(r1*r2); + + if (cos_theta > 1.0) cos_theta = 1.0; + if (cos_theta < -1.0) cos_theta = -1.0; + + // The method returns derivative with regards to cos(theta) + angle->born_matrix(atype, atom1, atom2, atom3, duang, du2ang); + // only right bond crossing the plane + if (right_cross && !left_cross) + { + double sgn = copysign(1.0, x_angle_right[dir] - pos); + dcos_theta[0] = sgn*(dx_right[0]*cos_theta/r2 + dx_left[0]/r1)/r2; + dcos_theta[1] = sgn*(dx_right[1]*cos_theta/r2 + dx_left[1]/r1)/r2; + dcos_theta[2] = sgn*(dx_right[2]*cos_theta/r2 + dx_left[2]/r1)/r2; + } + + // only left bond crossing the plane + if (!right_cross && left_cross) + { + double sgn = copysign(1.0, x_angle_left[dir] - pos); + dcos_theta[0] = -sgn*(dx_left[0]*cos_theta/r1 + dx_right[0]/r2)/r1; + dcos_theta[1] = -sgn*(dx_left[1]*cos_theta/r1 + dx_right[1]/r2)/r1; + dcos_theta[2] = -sgn*(dx_left[2]*cos_theta/r1 + dx_right[2]/r2)/r1; + } + + // both bonds crossing the plane + if (right_cross && left_cross) + { + // due to right bond + double sgn = copysign(1.0, x_angle_middle[dir] - pos); + dcos_theta[0] = -sgn*(dx_right[0]*cos_theta/r2 + dx_left[0]/r1)/r2; + dcos_theta[1] = -sgn*(dx_right[1]*cos_theta/r2 + dx_left[1]/r1)/r2; + dcos_theta[2] = -sgn*(dx_right[2]*cos_theta/r2 + dx_left[2]/r1)/r2; + + // due to left bond + dcos_theta[0] += sgn*(dx_left[0]*cos_theta/r1 + dx_right[0]/r2)/r1; + dcos_theta[1] += sgn*(dx_left[1]*cos_theta/r1 + dx_right[1]/r2)/r1; + dcos_theta[2] += sgn*(dx_left[2]*cos_theta/r1 + dx_right[2]/r2)/r1; + } + + // final contribution of the given angle term + local_contribution[ibin][0] += duang*dcos_theta[0]/area*nktv2p; + local_contribution[ibin][1] += duang*dcos_theta[1]/area*nktv2p; + local_contribution[ibin][2] += duang*dcos_theta[2]/area*nktv2p; + } + } + } + + // loop over the keywords and if necessary add the angle contribution + int m = 0; + while (m < nvalues) { + if (which[m] == CONF || which[m] == TOTAL || which[m] == ANGLE) { + for (int ibin = 0; ibin < nbins; ibin++) { + angle_local[ibin][m] = local_contribution[ibin][0]; + angle_local[ibin][m+1] = local_contribution[ibin][1]; + angle_local[ibin][m+2] = local_contribution[ibin][2]; + } + } + m += 3; + } +} + +/*------------------------------------------------------------------------ + compute dihedral contribution to pressure of local proc + -------------------------------------------------------------------------*/ + +void ComputeStressMopProfile::compute_dihedrals() +{ + int i, nd, atom1, atom2, atom3, atom4, imol, iatom; + tagint tagprev; + double vb1x, vb1y, vb1z, vb2x, vb2y, vb2z, vb3x, vb3y, vb3z; + double vb2xm, vb2ym, vb2zm; + double sb1, sb2, sb3, rb1, rb3, c0, b1mag2, b1mag, b2mag2; + double b2mag, b3mag2, b3mag, c2mag, ctmp, r12c1, c1mag, r12c2; + double s1, s2, s12, sc1, sc2, a11, a22, a33, a12, a13, a23; + double df[3], f1[3], f2[3], f3[3], f4[3]; + double c, sx2, sy2, sz2, sin2; + + double **x = atom->x; + tagint *tag = atom->tag; + int *num_dihedral = atom->num_dihedral; + tagint **dihedral_atom1 = atom->dihedral_atom1; + tagint **dihedral_atom2 = atom->dihedral_atom2; + tagint **dihedral_atom3 = atom->dihedral_atom3; + tagint **dihedral_atom4 = atom->dihedral_atom4; + int *mask = atom->mask; + + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; + + int nlocal = atom->nlocal; + int molecular = atom->molecular; + + // loop over all atoms and their dihedrals + + Dihedral *dihedral = force->dihedral; + + double dudih, du2dih; + + double diffx[3] = {0.0, 0.0, 0.0}; + double x_atom_1[3] = {0.0, 0.0, 0.0}; + double x_atom_2[3] = {0.0, 0.0, 0.0}; + double x_atom_3[3] = {0.0, 0.0, 0.0}; + double x_atom_4[3] = {0.0, 0.0, 0.0}; + + // initialization + for (int m = 0; m < nbins; m++) { + for (int i = 0; i < nvalues; i++) { + dihedral_local[m][i] = 0.0; + } + local_contribution[m][0] = 0.0; + local_contribution[m][1] = 0.0; + local_contribution[m][2] = 0.0; + } + + for (atom2 = 0; atom2 < nlocal; atom2++) { + if (!(mask[atom2] & groupbit)) continue; + + if (molecular == Atom::MOLECULAR) + nd = num_dihedral[atom2]; + else { + if (molindex[atom2] < 0) continue; + imol = molindex[atom2]; + iatom = molatom[atom2]; + nd = onemols[imol]->num_dihedral[iatom]; + } + + for (i = 0; i < nd; i++) { + if (molecular == 1) { + if (tag[atom2] != dihedral_atom2[atom2][i]) continue; + atom1 = atom->map(dihedral_atom1[atom2][i]); + atom3 = atom->map(dihedral_atom3[atom2][i]); + atom4 = atom->map(dihedral_atom4[atom2][i]); + } else { + if (tag[atom2] != onemols[imol]->dihedral_atom2[atom2][i]) continue; + tagprev = tag[atom2] - iatom - 1; + atom1 = atom->map(onemols[imol]->dihedral_atom1[atom2][i] + tagprev); + atom3 = atom->map(onemols[imol]->dihedral_atom3[atom2][i] + tagprev); + atom4 = atom->map(onemols[imol]->dihedral_atom4[atom2][i] + tagprev); + } + + if (atom1 < 0 || !(mask[atom1] & groupbit)) continue; + if (atom3 < 0 || !(mask[atom3] & groupbit)) continue; + if (atom4 < 0 || !(mask[atom4] & groupbit)) continue; + + for (int ibin = 0; ibinminimum_image(x_atom_1[0], x_atom_1[1], x_atom_1[2]); + x_atom_1[dir] += pos; + + // minimum image of atom2 with respect to atom1 + diffx[0] = x[atom2][0] - x_atom_1[0]; + diffx[1] = x[atom2][1] - x_atom_1[1]; + diffx[2] = x[atom2][2] - x_atom_1[2]; + domain->minimum_image(diffx[0], diffx[1], diffx[2]); + x_atom_2[0] = x_atom_1[0] + diffx[0]; + x_atom_2[1] = x_atom_1[1] + diffx[1]; + x_atom_2[2] = x_atom_1[2] + diffx[2]; + + // minimum image of atom3 with respect to atom2 + diffx[0] = x[atom3][0] - x_atom_2[0]; + diffx[1] = x[atom3][1] - x_atom_2[1]; + diffx[2] = x[atom3][2] - x_atom_2[2]; + domain->minimum_image(diffx[0], diffx[1], diffx[2]); + x_atom_3[0] = x_atom_2[0] + diffx[0]; + x_atom_3[1] = x_atom_2[1] + diffx[1]; + x_atom_3[2] = x_atom_2[2] + diffx[2]; + + // minimum image of atom3 with respect to atom2 + diffx[0] = x[atom4][0] - x_atom_3[0]; + diffx[1] = x[atom4][1] - x_atom_3[1]; + diffx[2] = x[atom4][2] - x_atom_3[2]; + domain->minimum_image(diffx[0], diffx[1], diffx[2]); + x_atom_4[0] = x_atom_3[0] + diffx[0]; + x_atom_4[1] = x_atom_3[1] + diffx[1]; + x_atom_4[2] = x_atom_3[2] + diffx[2]; + + // check if any bond vector crosses the plane of interest + double tau_right = (x_atom_2[dir] - pos) / (x_atom_2[dir] - x_atom_1[dir]); + double tau_middle = (x_atom_3[dir] - pos) / (x_atom_3[dir] - x_atom_2[dir]); + double tau_left = (x_atom_4[dir] - pos) / (x_atom_4[dir] - x_atom_3[dir]); + bool right_cross = ((tau_right >=0) && (tau_right <= 1)); + bool middle_cross = ((tau_middle >= 0) && (tau_middle <= 1)); + bool left_cross = ((tau_left >=0) && (tau_left <= 1)); + + // no bonds crossing the plane + if (!right_cross && !middle_cross && !left_cross) continue; + + dihedral->born_matrix(i, atom1, atom2, atom3, atom4, dudih, du2dih); + + // first bond + vb1x = x_atom_1[0] - x_atom_2[0]; + vb1y = x_atom_1[1] - x_atom_2[1]; + vb1z = x_atom_1[2] - x_atom_2[2]; + + // second bond + vb2x = x_atom_3[0] - x_atom_2[0]; + vb2y = x_atom_3[1] - x_atom_2[1]; + vb2z = x_atom_3[2] - x_atom_2[2]; + + vb2xm = -vb2x; + vb2ym = -vb2y; + vb2zm = -vb2z; + + // third bond + vb3x = x_atom_4[0] - x_atom_3[0]; + vb3y = x_atom_4[1] - x_atom_3[1]; + vb3z = x_atom_4[2] - x_atom_3[2]; + + // c0 calculation + sb1 = 1.0 / (vb1x*vb1x + vb1y*vb1y + vb1z*vb1z); + sb2 = 1.0 / (vb2x*vb2x + vb2y*vb2y + vb2z*vb2z); + sb3 = 1.0 / (vb3x*vb3x + vb3y*vb3y + vb3z*vb3z); + + rb1 = sqrt(sb1); + rb3 = sqrt(sb3); + + c0 = (vb1x*vb3x + vb1y*vb3y + vb1z*vb3z) * rb1*rb3; + // 1st and 2nd angle + b1mag2 = vb1x*vb1x + vb1y*vb1y + vb1z*vb1z; + b1mag = sqrt(b1mag2); + b2mag2 = vb2x*vb2x + vb2y*vb2y + vb2z*vb2z; + b2mag = sqrt(b2mag2); + b3mag2 = vb3x*vb3x + vb3y*vb3y + vb3z*vb3z; + b3mag = sqrt(b3mag2); + + ctmp = vb1x*vb2x + vb1y*vb2y + vb1z*vb2z; + r12c1 = 1.0 / (b1mag*b2mag); + c1mag = ctmp * r12c1; + + ctmp = vb2xm*vb3x + vb2ym*vb3y + vb2zm*vb3z; + r12c2 = 1.0 / (b2mag*b3mag); + c2mag = ctmp * r12c2; + + // cos and sin of 2 angles and final c + sin2 = MAX(1.0 - c1mag*c1mag,0.0); + sc1 = sqrt(sin2); + if (sc1 < SMALL) sc1 = SMALL; + sc1 = 1.0/sc1; + + sin2 = MAX(1.0 - c2mag*c2mag,0.0); + sc2 = sqrt(sin2); + if (sc2 < SMALL) sc2 = SMALL; + sc2 = 1.0/sc2; + + s1 = sc1 * sc1; + s2 = sc2 * sc2; + s12 = sc1 * sc2; + c = (c0 + c1mag*c2mag) * s12; + + // error check + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + + // forces on each particle + double a = dudih; + c = c * a; + s12 = s12 * a; + a11 = c*sb1*s1; + a22 = -sb2 * (2.0*c0*s12 - c*(s1+s2)); + a33 = c*sb3*s2; + a12 = -r12c1 * (c1mag*c*s1 + c2mag*s12); + a13 = -rb1*rb3*s12; + a23 = r12c2 * (c2mag*c*s2 + c1mag*s12); + + sx2 = a12*vb1x + a22*vb2x + a23*vb3x; + sy2 = a12*vb1y + a22*vb2y + a23*vb3y; + sz2 = a12*vb1z + a22*vb2z + a23*vb3z; + + f1[0] = a11*vb1x + a12*vb2x + a13*vb3x; + f1[1] = a11*vb1y + a12*vb2y + a13*vb3y; + f1[2] = a11*vb1z + a12*vb2z + a13*vb3z; + + f2[0] = -sx2 - f1[0]; + f2[1] = -sy2 - f1[1]; + f2[2] = -sz2 - f1[2]; + + f4[0] = a13*vb1x + a23*vb2x + a33*vb3x; + f4[1] = a13*vb1y + a23*vb2y + a33*vb3y; + f4[2] = a13*vb1z + a23*vb2z + a33*vb3z; + + f3[0] = sx2 - f4[0]; + f3[1] = sy2 - f4[1]; + f3[2] = sz2 - f4[2]; + + // only right bond crossing the plane + if (right_cross && !middle_cross && !left_cross) + { + double sgn = copysign(1.0, x_atom_1[dir] - pos); + df[0] = sgn * f1[0]; + df[1] = sgn * f1[1]; + df[2] = sgn * f1[2]; + } + + // only middle bond crossing the plane + if (!right_cross && middle_cross && !left_cross) + { + double sgn = copysign(1.0, x_atom_2[dir] - pos); + df[0] = sgn * (f2[0] + f1[0]); + df[1] = sgn * (f2[1] + f1[1]); + df[2] = sgn * (f2[2] + f1[2]); + } + + // only left bond crossing the plane + if (!right_cross && !middle_cross && left_cross) + { + double sgn = copysign(1.0, x_atom_4[dir] - pos); + df[0] = sgn * f4[0]; + df[1] = sgn * f4[1]; + df[2] = sgn * f4[2]; + } + + // only right & middle bonds crossing the plane + if (right_cross && middle_cross && !left_cross) + { + double sgn = copysign(1.0, x_atom_2[dir] - pos); + df[0] = sgn * f2[0]; + df[1] = sgn * f2[1]; + df[2] = sgn * f2[2]; + } + + // only right & left bonds crossing the plane + if (right_cross && !middle_cross && left_cross) + { + double sgn = copysign(1.0, x_atom_1[dir] - pos); + df[0] = sgn * (f1[0] + f4[0]); + df[1] = sgn * (f1[1] + f4[1]); + df[2] = sgn * (f1[2] + f4[2]); + } + + // only middle & left bonds crossing the plane + if (!right_cross && middle_cross && left_cross) + { + double sgn = copysign(1.0, x_atom_3[dir] - pos); + df[0] = sgn * f3[0]; + df[1] = sgn * f3[1]; + df[2] = sgn * f3[2]; + } + + // all three bonds crossing the plane + if (right_cross && middle_cross && left_cross) + { + double sgn = copysign(1.0, x_atom_1[dir] - pos); + df[0] = sgn * (f1[0] + f3[0]); + df[1] = sgn * (f1[1] + f3[1]); + df[2] = sgn * (f1[2] + f3[2]); + } + + local_contribution[ibin][0] += df[0]/area*nktv2p; + local_contribution[ibin][1] += df[1]/area*nktv2p; + local_contribution[ibin][2] += df[2]/area*nktv2p; + } + } + } + + // loop over the keywords and if necessary add the dihedral contribution + int m = 0; + while (m < nvalues) { + if ((which[m] == CONF) || (which[m] == TOTAL) || (which[m] == DIHEDRAL)) { + for (int ibin = 0; ibin < nbins; ibin++) { + dihedral_local[ibin][m] = local_contribution[ibin][0]; + dihedral_local[ibin][m+1] = local_contribution[ibin][1]; + dihedral_local[ibin][m+2] = local_contribution[ibin][2]; + } + } + m += 3; + } + +} + /* ---------------------------------------------------------------------- setup 1d bins and their extent and coordinates called at init() @@ -621,47 +1197,39 @@ void ComputeStressMopProfile::setup_bins() boxlo = domain->boxlo; boxhi = domain->boxhi; - if (originflag == LOWER) - origin = boxlo[dir]; - else if (originflag == UPPER) - origin = boxhi[dir]; - else if (originflag == CENTER) - origin = 0.5 * (boxlo[dir] + boxhi[dir]); + if ((origin > domain->boxhi[dir]) || (origin < domain->boxlo[dir])) + error->all(FLERR, "Origin of bins for compute stress/mop/profile is out of bounds"); - if (origin < boxlo[dir]) { - error->all(FLERR, "Origin of bins for compute stress/mop/profile is out of bounds"); - } else { - n = static_cast((origin - boxlo[dir]) * invdelta); - lo = origin - n * delta; - } - if (origin < boxhi[dir]) { - n = static_cast((boxhi[dir] - origin) * invdelta); - hi = origin + n * delta; - } else { - error->all(FLERR, "Origin of bins for compute stress/mop/profile is out of bounds"); - } + n = static_cast ((origin - boxlo[dir]) * invdelta); + lo = origin - n*delta; + + n = static_cast ((boxhi[dir] - origin) * invdelta); + hi = origin + n*delta; offset = lo; nbins = static_cast((hi - lo) * invdelta + 1.5); - // allocate bin arrays - - memory->create(coord, nbins, 1, "stress/mop/profile:coord"); - memory->create(coordp, nbins, 1, "stress/mop/profile:coordp"); + //allocate bin arrays + memory->create(coord, nbins, "stress/mop/profile:coord"); + memory->create(coordp, nbins, "stress/mop/profile:coordp"); memory->create(values_local, nbins, nvalues, "stress/mop/profile:values_local"); memory->create(values_global, nbins, nvalues, "stress/mop/profile:values_global"); memory->create(bond_local, nbins, nvalues, "stress/mop/profile:bond_local"); memory->create(bond_global, nbins, nvalues, "stress/mop/profile:bond_global"); + memory->create(angle_local, nbins, nvalues, "stress/mop/profile:angle_local"); + memory->create(angle_global, nbins, nvalues, "stress/mop/profile:angle_global"); + memory->create(dihedral_local,nbins,nvalues,"stress/mop/profile:dihedral_local"); + memory->create(dihedral_global,nbins,nvalues,"stress/mop/profile:dihedral_global"); memory->create(local_contribution, nbins, 3, "stress/mop/profile:local_contribution"); // set bin coordinates for (i = 0; i < nbins; i++) { - coord[i][0] = offset + i * delta; - if (coord[i][0] < (domain->boxlo[dir] + domain->prd_half[dir])) { - coordp[i][0] = coord[i][0] + domain->prd[dir]; + coord[i] = offset + i * delta; + if (coord[i] < (domain->boxlo[dir] + domain->prd_half[dir])) { + coordp[i] = coord[i] + domain->prd[dir]; } else { - coordp[i][0] = coord[i][0] - domain->prd[dir]; + coordp[i] = coord[i] - domain->prd[dir]; } } } diff --git a/src/EXTRA-COMPUTE/compute_stress_mop_profile.h b/src/EXTRA-COMPUTE/compute_stress_mop_profile.h index 2b0ffef0f8..b9b97617c0 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop_profile.h +++ b/src/EXTRA-COMPUTE/compute_stress_mop_profile.h @@ -39,19 +39,22 @@ class ComputeStressMopProfile : public Compute { private: void compute_pairs(); void compute_bonds(); + void compute_angles(); + void compute_dihedrals(); void setup_bins(); int nvalues, dir; int *which; - int bondflag; + int bondflag, angleflag, dihedralflag; - int originflag; double origin, delta, offset, invdelta; int nbins; - double **coord, **coordp; + double *coord, *coordp; double **values_local, **values_global; double **bond_local, **bond_global; + double **angle_local, **angle_global; + double **dihedral_local, **dihedral_global; double **local_contribution; double dt, nktv2p, ftm2v; diff --git a/src/EXTRA-MOLECULE/angle_cosine_periodic.cpp b/src/EXTRA-MOLECULE/angle_cosine_periodic.cpp index 15d0575f6d..34a8e9d8e5 100644 --- a/src/EXTRA-MOLECULE/angle_cosine_periodic.cpp +++ b/src/EXTRA-MOLECULE/angle_cosine_periodic.cpp @@ -38,7 +38,10 @@ using namespace MathSpecial; /* ---------------------------------------------------------------------- */ -AngleCosinePeriodic::AngleCosinePeriodic(LAMMPS *lmp) : Angle(lmp) {} +AngleCosinePeriodic::AngleCosinePeriodic(LAMMPS *lmp) : Angle(lmp) +{ + born_matrix_enable = 1; +} /* ---------------------------------------------------------------------- */ @@ -298,3 +301,38 @@ double AngleCosinePeriodic::single(int type, int i1, int i2, int i3) c = cos(acos(c)*multiplicity[type]); return 2.0*k[type]*(1.0-b[type]*powsign(multiplicity[type])*c); } + +/* ---------------------------------------------------------------------- */ + +void AngleCosinePeriodic::born_matrix(int type, int i1, int i2, int i3, double &du, double &du2) +{ + double **x = atom->x; + + double delx1 = x[i1][0] - x[i2][0]; + double dely1 = x[i1][1] - x[i2][1]; + double delz1 = x[i1][2] - x[i2][2]; + domain->minimum_image(delx1,dely1,delz1); + double r1 = sqrt(delx1*delx1 + dely1*dely1 + delz1*delz1); + + double delx2 = x[i3][0] - x[i2][0]; + double dely2 = x[i3][1] - x[i2][1]; + double delz2 = x[i3][2] - x[i2][2]; + domain->minimum_image(delx2,dely2,delz2); + double r2 = sqrt(delx2*delx2 + dely2*dely2 + delz2*delz2); + + double c = delx1*delx2 + dely1*dely2 + delz1*delz2; + c /= r1*r2; + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + double theta = acos(c); + + double s = sqrt(1.0 - c*c); + if (s < SMALL) s = SMALL; + s = 1.0/s; + + double m_angle = multiplicity[type] * theta; + double prefactor = -2.0 * k[type] * b[type] * powsign(multiplicity[type]) * multiplicity[type]; + + du = prefactor * sin(m_angle) / s; + du2 = prefactor * (c * sin(m_angle) - s * cos(m_angle) * multiplicity[type]) / (s * s * s); +} diff --git a/src/EXTRA-MOLECULE/angle_cosine_periodic.h b/src/EXTRA-MOLECULE/angle_cosine_periodic.h index 4e584b4543..f04ed04784 100644 --- a/src/EXTRA-MOLECULE/angle_cosine_periodic.h +++ b/src/EXTRA-MOLECULE/angle_cosine_periodic.h @@ -35,6 +35,7 @@ class AngleCosinePeriodic : public Angle { void read_restart(FILE *) override; void write_data(FILE *) override; double single(int, int, int, int) override; + void born_matrix(int type, int i1, int i2, int i3, double &du, double &du2) override; protected: double *k; diff --git a/src/EXTRA-MOLECULE/angle_fourier.cpp b/src/EXTRA-MOLECULE/angle_fourier.cpp index 549da0c196..c7eb3d4fe4 100644 --- a/src/EXTRA-MOLECULE/angle_fourier.cpp +++ b/src/EXTRA-MOLECULE/angle_fourier.cpp @@ -39,6 +39,7 @@ using namespace MathConst; AngleFourier::AngleFourier(LAMMPS *lmp) : Angle(lmp) { + born_matrix_enable = 1; k = nullptr; C0 = nullptr; C1 = nullptr; diff --git a/src/EXTRA-MOLECULE/angle_quartic.cpp b/src/EXTRA-MOLECULE/angle_quartic.cpp index f28e209a77..eaccdbe608 100644 --- a/src/EXTRA-MOLECULE/angle_quartic.cpp +++ b/src/EXTRA-MOLECULE/angle_quartic.cpp @@ -37,7 +37,10 @@ using namespace MathConst; /* ---------------------------------------------------------------------- */ -AngleQuartic::AngleQuartic(LAMMPS *lmp) : Angle(lmp) {} +AngleQuartic::AngleQuartic(LAMMPS *lmp) : Angle(lmp) +{ + born_matrix_enable = 1; +} /* ---------------------------------------------------------------------- */ @@ -286,3 +289,39 @@ double AngleQuartic::single(int type, int i1, int i2, int i3) double dtheta4 = dtheta3 * dtheta; return k2[type] * dtheta2 + k3[type] * dtheta3 + k4[type] * dtheta4; } + +/* ---------------------------------------------------------------------- */ + +void AngleQuartic::born_matrix(int type, int i1, int i2, int i3, double &du, double &du2) +{ + double **x = atom->x; + + double delx1 = x[i1][0] - x[i2][0]; + double dely1 = x[i1][1] - x[i2][1]; + double delz1 = x[i1][2] - x[i2][2]; + domain->minimum_image(delx1,dely1,delz1); + double r1 = sqrt(delx1*delx1 + dely1*dely1 + delz1*delz1); + + double delx2 = x[i3][0] - x[i2][0]; + double dely2 = x[i3][1] - x[i2][1]; + double delz2 = x[i3][2] - x[i2][2]; + domain->minimum_image(delx2,dely2,delz2); + double r2 = sqrt(delx2*delx2 + dely2*dely2 + delz2*delz2); + + double c = delx1*delx2 + dely1*dely2 + delz1*delz2; + c /= r1*r2; + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + double theta = acos(c); + + double s = sqrt(1.0 - c*c); + if (s < SMALL) s = SMALL; + + double dtheta = theta - theta0[type]; + double dtheta2 = dtheta * dtheta; + double dtheta3 = dtheta2 * dtheta; + + du = -(2.0 * k2[type] * dtheta + 3.0 * k3[type] * dtheta2 + 4.0 * k4[type] * dtheta3) / s; + du2 = (2.0 * k2[type] + 6.0 * k3[type] * dtheta + 12.0 * k4[type] * dtheta2) / (s*s) - + (2.0 * k2[type] * dtheta + 3.0 * k3[type] * dtheta2 + 4.0 * k4[type] * dtheta3) * c / (s*s*s); +} diff --git a/src/EXTRA-MOLECULE/angle_quartic.h b/src/EXTRA-MOLECULE/angle_quartic.h index 3f0396f27b..7de51b24d1 100644 --- a/src/EXTRA-MOLECULE/angle_quartic.h +++ b/src/EXTRA-MOLECULE/angle_quartic.h @@ -35,6 +35,7 @@ class AngleQuartic : public Angle { void read_restart(FILE *) override; void write_data(FILE *) override; double single(int, int, int, int) override; + void born_matrix(int type, int i1, int i2, int i3, double &du, double &du2) override; protected: double *k2, *k3, *k4, *theta0; diff --git a/src/EXTRA-MOLECULE/bond_gaussian.cpp b/src/EXTRA-MOLECULE/bond_gaussian.cpp index baca0b6e1a..9a8546e278 100644 --- a/src/EXTRA-MOLECULE/bond_gaussian.cpp +++ b/src/EXTRA-MOLECULE/bond_gaussian.cpp @@ -35,6 +35,7 @@ BondGaussian::BondGaussian(LAMMPS *lmp) : Bond(lmp), nterms(nullptr), bond_temperature(nullptr), alpha(nullptr), width(nullptr), r0(nullptr) { + born_matrix_enable = 1; } /* ---------------------------------------------------------------------- */ @@ -294,3 +295,45 @@ double BondGaussian::single(int type, double rsq, int /*i*/, int /*j*/, double & return -(force->boltz * bond_temperature[type]) * log(sum_g_i); } + +/* ---------------------------------------------------------------------- */ + +void BondGaussian::born_matrix(int type, double rsq, int /*i*/, int /*j*/, double &du, double &du2) +{ + double r = sqrt(rsq); + + // first derivative of energy with respect to distance + double sum_g_i = 0.0; + double sum_numerator = 0.0; + for (int i = 0; i < nterms[type]; i++) { + double dr = r - r0[type][i]; + double prefactor = (alpha[type][i] / (width[type][i] * sqrt(MY_PI2))); + double exponent = -2 * dr * dr / (width[type][i] * width[type][i]); + double g_i = prefactor * exp(exponent); + sum_g_i += g_i; + sum_numerator += g_i * dr / (width[type][i] * width[type][i]); + } + + if (sum_g_i < SMALL) sum_g_i = SMALL; + du = 4.0 * (force->boltz * bond_temperature[type]) * (sum_numerator / sum_g_i); + + // second derivative of energy with respect to distance + sum_g_i = 0.0; + double sum_dg_i = 0.0; + double sum_d2g_i = 0.0; + for (int i = 0; i < nterms[type]; i++) { + double dr = r - r0[type][i]; + double prefactor = (alpha[type][i] / (width[type][i] * sqrt(MY_PI2))); + double exponent = -2 * dr * dr / (width[type][i] * width[type][i]); + double g_i = prefactor * exp(exponent); + sum_g_i += g_i; + sum_dg_i -= 4.0 * g_i * dr / pow(width[type][i], 2); + sum_d2g_i += 4.0 * g_i * (4.0 * pow(r0[type][i], 2) - 8.0 * r0[type][i] * r - pow(width[type][i], 2) + 4.0 * r * r) / pow(width[type][i], 4) ; + } + + if (sum_g_i < SMALL) sum_g_i = SMALL; + double numerator = sum_d2g_i*sum_g_i - sum_dg_i*sum_dg_i; + double denominator = sum_g_i * sum_g_i; + + du2 = - (force->boltz * bond_temperature[type]) * numerator / denominator; +} diff --git a/src/EXTRA-MOLECULE/bond_gaussian.h b/src/EXTRA-MOLECULE/bond_gaussian.h index 7af6f1f4d9..e466df47d4 100644 --- a/src/EXTRA-MOLECULE/bond_gaussian.h +++ b/src/EXTRA-MOLECULE/bond_gaussian.h @@ -35,6 +35,7 @@ class BondGaussian : public Bond { void read_restart(FILE *) override; void write_data(FILE *) override; double single(int, double, int, int, double &) override; + void born_matrix(int, double, int, int, double &, double &) override; protected: int *nterms; diff --git a/src/EXTRA-MOLECULE/bond_harmonic_shift_cut.cpp b/src/EXTRA-MOLECULE/bond_harmonic_shift_cut.cpp index fedcb95ee8..ebcfdb0258 100644 --- a/src/EXTRA-MOLECULE/bond_harmonic_shift_cut.cpp +++ b/src/EXTRA-MOLECULE/bond_harmonic_shift_cut.cpp @@ -31,7 +31,10 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -BondHarmonicShiftCut::BondHarmonicShiftCut(LAMMPS *lmp) : Bond(lmp) {} +BondHarmonicShiftCut::BondHarmonicShiftCut(LAMMPS *lmp) : Bond(lmp) +{ + born_matrix_enable = 1; +} /* ---------------------------------------------------------------------- */ @@ -219,3 +222,19 @@ double BondHarmonicShiftCut::single(int type, double rsq, int /*i*/, int /*j*/, fforce = -2.0*k[type]*dr/r; return k[type]*(dr*dr - dr2*dr2); } + +/* ---------------------------------------------------------------------- */ + +void BondHarmonicShiftCut::born_matrix(int type, double rsq, int /*i*/, int /*j*/, double &du, double &du2) +{ + du = 0.0; + du2 = 0.0; + + double r = sqrt(rsq); + if (r>r1[type]) return; + + double dr = r - r0[type]; + + du2 = 2 * k[type]; + if (r > 0.0) du = du2 * dr; +} diff --git a/src/EXTRA-MOLECULE/bond_harmonic_shift_cut.h b/src/EXTRA-MOLECULE/bond_harmonic_shift_cut.h index 752ac010d9..09d6ab5330 100644 --- a/src/EXTRA-MOLECULE/bond_harmonic_shift_cut.h +++ b/src/EXTRA-MOLECULE/bond_harmonic_shift_cut.h @@ -35,6 +35,7 @@ class BondHarmonicShiftCut : public Bond { void read_restart(FILE *) override; void write_data(FILE *) override; double single(int, double, int, int, double &) override; + void born_matrix(int, double, int, int, double &, double &) override; protected: double *k, *r0, *r1; diff --git a/src/EXTRA-MOLECULE/dihedral_helix.cpp b/src/EXTRA-MOLECULE/dihedral_helix.cpp index 059bef74a4..1d99de6ba9 100644 --- a/src/EXTRA-MOLECULE/dihedral_helix.cpp +++ b/src/EXTRA-MOLECULE/dihedral_helix.cpp @@ -41,6 +41,7 @@ using namespace MathConst; DihedralHelix::DihedralHelix(LAMMPS *lmp) : Dihedral(lmp) { writedata = 1; + born_matrix_enable = 1; } /* ---------------------------------------------------------------------- */ @@ -324,3 +325,108 @@ void DihedralHelix::write_data(FILE *fp) for (int i = 1; i <= atom->ndihedraltypes; i++) fprintf(fp,"%d %g %g %g\n",i,aphi[i],bphi[i],cphi[i]); } + +/* ----------------------------------------------------------------------*/ + +void DihedralHelix::born_matrix(int nd, int i1, int i2, int i3, int i4, + double &du, double &du2) +{ + double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm; + double sb1,sb3,rb1,rb3,c0,b1mag2,b1mag,b2mag2; + double b2mag,b3mag2,b3mag,ctmp,r12c1,c1mag,r12c2; + double c2mag,sc1,sc2,s12,c; + double cx,cy,cz,cmag,dx,phi,si,siinv,sin2; + + int **dihedrallist = neighbor->dihedrallist; + double **x = atom->x; + + int type = dihedrallist[nd][4]; + + // 1st bond + + vb1x = x[i1][0] - x[i2][0]; + vb1y = x[i1][1] - x[i2][1]; + vb1z = x[i1][2] - x[i2][2]; + + // 2nd bond + + vb2x = x[i3][0] - x[i2][0]; + vb2y = x[i3][1] - x[i2][1]; + vb2z = x[i3][2] - x[i2][2]; + + vb2xm = -vb2x; + vb2ym = -vb2y; + vb2zm = -vb2z; + + // 3rd bond + + vb3x = x[i4][0] - x[i3][0]; + vb3y = x[i4][1] - x[i3][1]; + vb3z = x[i4][2] - x[i3][2]; + + // c0 calculation + + sb1 = 1.0 / (vb1x*vb1x + vb1y*vb1y + vb1z*vb1z); + sb3 = 1.0 / (vb3x*vb3x + vb3y*vb3y + vb3z*vb3z); + + rb1 = sqrt(sb1); + rb3 = sqrt(sb3); + + c0 = (vb1x*vb3x + vb1y*vb3y + vb1z*vb3z) * rb1*rb3; + + // 1st and 2nd angle + + b1mag2 = vb1x*vb1x + vb1y*vb1y + vb1z*vb1z; + b1mag = sqrt(b1mag2); + b2mag2 = vb2x*vb2x + vb2y*vb2y + vb2z*vb2z; + b2mag = sqrt(b2mag2); + b3mag2 = vb3x*vb3x + vb3y*vb3y + vb3z*vb3z; + b3mag = sqrt(b3mag2); + + ctmp = vb1x*vb2x + vb1y*vb2y + vb1z*vb2z; + r12c1 = 1.0 / (b1mag*b2mag); + c1mag = ctmp * r12c1; + + ctmp = vb2xm*vb3x + vb2ym*vb3y + vb2zm*vb3z; + r12c2 = 1.0 / (b2mag*b3mag); + c2mag = ctmp * r12c2; + + // cos and sin of 2 angles and final c + + sin2 = MAX(1.0 - c1mag*c1mag,0.0); + sc1 = sqrt(sin2); + if (sc1 < SMALL) sc1 = SMALL; + sc1 = 1.0/sc1; + + sin2 = MAX(1.0 - c2mag*c2mag,0.0); + sc2 = sqrt(sin2); + if (sc2 < SMALL) sc2 = SMALL; + sc2 = 1.0/sc2; + + s12 = sc1 * sc2; + c = (c0 + c1mag*c2mag) * s12; + + cx = vb1y*vb2z - vb1z*vb2y; + cy = vb1z*vb2x - vb1x*vb2z; + cz = vb1x*vb2y - vb1y*vb2x; + cmag = sqrt(cx*cx + cy*cy + cz*cz); + dx = (cx*vb3x + cy*vb3y + cz*vb3z)/cmag/b3mag; + + // error check + + if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) problem(FLERR, i1, i2, i3, i4); + + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + + phi = acos(c); + if (dx > 0.0) phi *= -1.0; + si = sin(phi); + if (fabs(si) < SMALLER) si = SMALLER; + siinv = 1.0/si; + + du = -aphi[type] + 3.0*bphi[type]*sin(3.0*phi)*siinv + + cphi[type]*sin(phi + MY_PI4)*siinv; + du2 = -(9.0*bphi[type]*cos(3.0*phi) + cphi[type]*cos(phi + MY_PI4))*siinv*siinv + + (3.0*bphi[type]*sin(3.0*phi) + cphi[type]*sin(phi + MY_PI4))*c*siinv*siinv*siinv; +} diff --git a/src/EXTRA-MOLECULE/dihedral_helix.h b/src/EXTRA-MOLECULE/dihedral_helix.h index 436895c5c3..172a8c3469 100644 --- a/src/EXTRA-MOLECULE/dihedral_helix.h +++ b/src/EXTRA-MOLECULE/dihedral_helix.h @@ -33,6 +33,7 @@ class DihedralHelix : public Dihedral { void write_restart(FILE *) override; void read_restart(FILE *) override; void write_data(FILE *) override; + void born_matrix(int, int, int, int, int, double &, double &) override; protected: double *aphi, *bphi, *cphi; diff --git a/src/EXTRA-MOLECULE/dihedral_quadratic.cpp b/src/EXTRA-MOLECULE/dihedral_quadratic.cpp index cbe9e3e3a2..f576e6efdd 100644 --- a/src/EXTRA-MOLECULE/dihedral_quadratic.cpp +++ b/src/EXTRA-MOLECULE/dihedral_quadratic.cpp @@ -41,6 +41,7 @@ using namespace MathConst; DihedralQuadratic::DihedralQuadratic(LAMMPS *lmp) : Dihedral(lmp) { writedata = 1; + born_matrix_enable = 1; } /* ---------------------------------------------------------------------- */ @@ -327,3 +328,112 @@ void DihedralQuadratic::write_data(FILE *fp) for (int i = 1; i <= atom->ndihedraltypes; i++) fprintf(fp,"%d %g %g \n",i,k[i],phi0[i]*180.0/MY_PI); } + +/* ----------------------------------------------------------------------*/ + +void DihedralQuadratic::born_matrix(int nd, int i1, int i2, int i3, int i4, + double &du, double &du2) +{ + double vb1x,vb1y,vb1z,vb2x,vb2y,vb2z,vb3x,vb3y,vb3z,vb2xm,vb2ym,vb2zm; + double sb1,sb3,rb1,rb3,c0,b1mag2,b1mag,b2mag2; + double b2mag,b3mag2,b3mag,ctmp,r12c1,c1mag,r12c2; + double c2mag,sc1,sc2,s12,c; + double s1,s2,cx,cy,cz,cmag,dx,phi,si,siinv,sin2; + + int **dihedrallist = neighbor->dihedrallist; + double **x = atom->x; + + int type = dihedrallist[nd][4]; + + // 1st bond + + vb1x = x[i1][0] - x[i2][0]; + vb1y = x[i1][1] - x[i2][1]; + vb1z = x[i1][2] - x[i2][2]; + + // 2nd bond + + vb2x = x[i3][0] - x[i2][0]; + vb2y = x[i3][1] - x[i2][1]; + vb2z = x[i3][2] - x[i2][2]; + + vb2xm = -vb2x; + vb2ym = -vb2y; + vb2zm = -vb2z; + + // 3rd bond + vb3x = x[i4][0] - x[i3][0]; + vb3y = x[i4][1] - x[i3][1]; + vb3z = x[i4][2] - x[i3][2]; + + // c0 calculation + + sb1 = 1.0 / (vb1x*vb1x + vb1y*vb1y + vb1z*vb1z); + sb3 = 1.0 / (vb3x*vb3x + vb3y*vb3y + vb3z*vb3z); + + rb1 = sqrt(sb1); + rb3 = sqrt(sb3); + + c0 = (vb1x*vb3x + vb1y*vb3y + vb1z*vb3z) * rb1*rb3; + + // 1st and 2nd angle + + b1mag2 = vb1x*vb1x + vb1y*vb1y + vb1z*vb1z; + b1mag = sqrt(b1mag2); + b2mag2 = vb2x*vb2x + vb2y*vb2y + vb2z*vb2z; + b2mag = sqrt(b2mag2); + b3mag2 = vb3x*vb3x + vb3y*vb3y + vb3z*vb3z; + b3mag = sqrt(b3mag2); + + ctmp = vb1x*vb2x + vb1y*vb2y + vb1z*vb2z; + r12c1 = 1.0 / (b1mag*b2mag); + c1mag = ctmp * r12c1; + + ctmp = vb2xm*vb3x + vb2ym*vb3y + vb2zm*vb3z; + r12c2 = 1.0 / (b2mag*b3mag); + c2mag = ctmp * r12c2; + + // cos and sin of 2 angles and final c + + sin2 = MAX(1.0 - c1mag*c1mag,0.0); + sc1 = sqrt(sin2); + if (sc1 < SMALL) sc1 = SMALL; + sc1 = 1.0/sc1; + + sin2 = MAX(1.0 - c2mag*c2mag,0.0); + sc2 = sqrt(sin2); + if (sc2 < SMALL) sc2 = SMALL; + sc2 = 1.0/sc2; + + s1 = sc1 * sc1; + s2 = sc2 * sc2; + s12 = sc1 * sc2; + c = (c0 + c1mag*c2mag) * s12; + + cx = vb1y*vb2z - vb1z*vb2y; + cy = vb1z*vb2x - vb1x*vb2z; + cz = vb1x*vb2y - vb1y*vb2x; + cmag = sqrt(cx*cx + cy*cy + cz*cz); + dx = (cx*vb3x + cy*vb3y + cz*vb3z)/cmag/b3mag; + + // error check + + if (c > 1.0 + TOLERANCE || c < (-1.0 - TOLERANCE)) + problem(FLERR, i1, i2, i3, i4); + + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + + phi = acos(c); + if (dx > 0.0) phi *= -1.0; + si = sin(phi); + if (fabs(si) < SMALLER) si = SMALLER; + siinv = 1.0/si; + + double dphi = phi-phi0[type]; + if (dphi > MY_PI) dphi -= 2*MY_PI; + else if (dphi < -MY_PI) dphi += 2*MY_PI; + + du = - 2.0 * k[type] * dphi * siinv; + du2 = 2.0 * k[type] * siinv * siinv * ( 1.0 - dphi * c * siinv) ; +} diff --git a/src/EXTRA-MOLECULE/dihedral_quadratic.h b/src/EXTRA-MOLECULE/dihedral_quadratic.h index 90d8c3be6e..89f6fa3b25 100644 --- a/src/EXTRA-MOLECULE/dihedral_quadratic.h +++ b/src/EXTRA-MOLECULE/dihedral_quadratic.h @@ -33,6 +33,7 @@ class DihedralQuadratic : public Dihedral { void write_restart(FILE *) override; void read_restart(FILE *) override; void write_data(FILE *) override; + void born_matrix(int, int, int, int, int, double &, double &) override; protected: double *k, *phi0; diff --git a/src/KOKKOS/Install.sh b/src/KOKKOS/Install.sh index 4be74334c9..af80420d7a 100755 --- a/src/KOKKOS/Install.sh +++ b/src/KOKKOS/Install.sh @@ -165,6 +165,8 @@ action fix_qeq_reaxff_kokkos.cpp fix_qeq_reaxff.cpp action fix_qeq_reaxff_kokkos.h fix_qeq_reaxff.h action fix_reaxff_bonds_kokkos.cpp fix_reaxff_bonds.cpp action fix_reaxff_bonds_kokkos.h fix_reaxff_bonds.h +action compute_reaxff_atom_kokkos.cpp compute_reaxff_atom.cpp +action compute_reaxff_atom_kokkos.h compute_reaxff_atom.h action fix_reaxff_species_kokkos.cpp fix_reaxff_species.cpp action fix_reaxff_species_kokkos.h fix_reaxff_species.h action fix_rx_kokkos.cpp fix_rx.cpp diff --git a/src/KOKKOS/compute_reaxff_atom_kokkos.cpp b/src/KOKKOS/compute_reaxff_atom_kokkos.cpp new file mode 100644 index 0000000000..8dbcb9441e --- /dev/null +++ b/src/KOKKOS/compute_reaxff_atom_kokkos.cpp @@ -0,0 +1,195 @@ +// clang-format off +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + 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. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Richard Berger (LANL) +------------------------------------------------------------------------- */ + +#include "compute_reaxff_atom_kokkos.h" +#include "atom.h" +#include "molecule.h" +#include "update.h" +#include "force.h" +#include "memory.h" +#include "error.h" +#include "neigh_list.h" + +#include "memory_kokkos.h" +#include "pair_reaxff_kokkos.h" +#include "reaxff_api.h" + +using namespace LAMMPS_NS; +using namespace ReaxFF; + +/* ---------------------------------------------------------------------- */ + +template +ComputeReaxFFAtomKokkos::ComputeReaxFFAtomKokkos(LAMMPS *lmp, int narg, char **arg) : + ComputeReaxFFAtom(lmp, narg, arg), + nbuf(-1), buf(nullptr) +{ + kokkosable = 1; +} + +/* ---------------------------------------------------------------------- */ + +template +ComputeReaxFFAtomKokkos::~ComputeReaxFFAtomKokkos() +{ + memoryKK->destroy_kokkos(k_buf, buf); +} + +/* ---------------------------------------------------------------------- */ + +template +void ComputeReaxFFAtomKokkos::init() +{ + ComputeReaxFFAtom::init(); + + if (!reaxff || !reaxff->kokkosable) { + error->all(FLERR,"Cannot use compute reaxff/atom/kk without " + "pair_style reaxff/kk"); + } +} + +/* ---------------------------------------------------------------------- */ +template +void ComputeReaxFFAtomKokkos::compute_bonds() +{ + if (atom->nlocal > nlocal) { + memory->destroy(array_atom); + nlocal = atom->nlocal; + memory->create(array_atom, nlocal, 3, "reaxff/atom:array_atom"); + } + + // retrieve bond information from kokkos pair style. the data potentially + // lives on device. it is copied into buf on the host in a condensed format + // compute_local and compute_atom then expand the data from this buffer into + // appropiate arrays for consumption by others (e.g. dump local, dump custom + // or library interface) + + int maxnumbonds = 0; + if (reaxff->execution_space == Device) + device_pair()->FindBond(maxnumbonds, groupbit); + else + host_pair()->FindBond(maxnumbonds, groupbit); + + nbuf = ((store_bonds ? maxnumbonds*2 : 0) + 3)*nlocal; + + if (!buf || k_buf.extent(0) < nbuf) { + memoryKK->destroy_kokkos(k_buf, buf); + memoryKK->create_kokkos(k_buf, buf, nbuf, "reaxff/atom:buf"); + } + + // Pass information to buffer, will sync to host + + int nbuf_local; + if (reaxff->execution_space == Device) + device_pair()->PackReducedBondBuffer(k_buf, nbuf_local, store_bonds); + else + host_pair()->PackReducedBondBuffer(k_buf, nbuf_local, store_bonds); + + // Extract number of bonds from buffer + + nbonds = 0; + int j = 0; + for (int i = 0; i < nlocal; i++) { + int numbonds = static_cast(buf[j+2]); + nbonds += numbonds; + j += (store_bonds ? 2*numbonds : 0) + 3; + } +} + +/* ---------------------------------------------------------------------- */ + +template +void ComputeReaxFFAtomKokkos::compute_local() +{ + invoked_local = update->ntimestep; + + if (invoked_bonds < update->ntimestep) + compute_bonds(); + + if (nbonds > prev_nbonds) { + // grow array_local + memory->destroy(array_local); + memory->create(array_local, nbonds, 3, "reaxff/atom:array_local"); + prev_nbonds = nbonds; + } + + size_local_rows = nbonds; + + // extract local bond information from buffer + + int b = 0; + int j = 0; + auto tag = atom->tag; + + for (int i = 0; i < nlocal; ++i) { + const int numbonds = static_cast(buf[j+2]); + const int neigh_offset = j + 3; + const int bo_offset = neigh_offset + numbonds; + for (int k = 0; k < numbonds; k++) { + auto bond = array_local[b++]; + bond[0] = tag[i]; + bond[1] = static_cast (buf[neigh_offset+k]); + bond[2] = buf[bo_offset+k]; + } + j += 2*numbonds + 3; + } +} + +/* ---------------------------------------------------------------------- */ + +template +void ComputeReaxFFAtomKokkos::compute_peratom() +{ + invoked_peratom = update->ntimestep; + + if (invoked_bonds < update->ntimestep) + compute_bonds(); + + // extract peratom bond information from buffer + + int j = 0; + for (int i = 0; i < nlocal; ++i) { + auto ptr = array_atom[i]; + int numbonds = static_cast(buf[j+2]); + ptr[0] = buf[j]; // sbo + ptr[1] = buf[j+1]; // nlp + ptr[2] = numbonds; + j += (store_bonds ? 2*numbonds : 0) + 3; + } +} + +/* ---------------------------------------------------------------------- + memory usage of local data +------------------------------------------------------------------------- */ + +template +double ComputeReaxFFAtomKokkos::memory_usage() +{ + double bytes = (double)(nlocal*3) * sizeof(double); + if (store_bonds) + bytes += (double)(nbonds*3) * sizeof(double); + bytes += (double)(nbuf > 0 ? nbuf * sizeof(double) : 0); + return bytes; +} + +namespace LAMMPS_NS { +template class ComputeReaxFFAtomKokkos; +#ifdef LMP_KOKKOS_GPU +template class ComputeReaxFFAtomKokkos; +#endif +} diff --git a/src/KOKKOS/compute_reaxff_atom_kokkos.h b/src/KOKKOS/compute_reaxff_atom_kokkos.h new file mode 100644 index 0000000000..7037c7e308 --- /dev/null +++ b/src/KOKKOS/compute_reaxff_atom_kokkos.h @@ -0,0 +1,66 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + 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. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Richard Berger (LANL) +------------------------------------------------------------------------- */ + +#ifdef COMPUTE_CLASS +// clang-format off +ComputeStyle(reaxff/atom/kk,ComputeReaxFFAtomKokkos); +ComputeStyle(reaxff/atom/kk/device,ComputeReaxFFAtomKokkos); +ComputeStyle(reaxff/atom/kk/host,ComputeReaxFFAtomKokkos); +// clang-format on +#else + +#ifndef LMP_COMPUTE_REAXFF_BONDS_KOKKOS_H +#define LMP_COMPUTE_REAXFF_BONDS_KOKKOS_H + +#include "compute_reaxff_atom.h" +#include "pair_reaxff_kokkos.h" +#include "kokkos_type.h" + +namespace LAMMPS_NS { + +template +class ComputeReaxFFAtomKokkos : public ComputeReaxFFAtom { + public: + using device_type = DeviceType; + using AT = ArrayTypes; + + ComputeReaxFFAtomKokkos(class LAMMPS *, int, char **); + ~ComputeReaxFFAtomKokkos() override; + void init() override; + void compute_local() override; + void compute_peratom() override; + void compute_bonds() override; + double memory_usage() override; + + private: + int nbuf; + double *buf; + typename AT::tdual_float_1d k_buf; + + auto device_pair() { + return static_cast*>(reaxff); + } + + auto host_pair() { + return static_cast*>(reaxff); + } +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/KOKKOS/pair_reaxff_kokkos.cpp b/src/KOKKOS/pair_reaxff_kokkos.cpp index c7d54b80cd..11a40970c2 100644 --- a/src/KOKKOS/pair_reaxff_kokkos.cpp +++ b/src/KOKKOS/pair_reaxff_kokkos.cpp @@ -4162,22 +4162,23 @@ double PairReaxFFKokkos::memory_usage() /* ---------------------------------------------------------------------- */ template -void PairReaxFFKokkos::FindBond(int &numbonds) +void PairReaxFFKokkos::FindBond(int &numbonds, int groupbit) { copymode = 1; Kokkos::parallel_for(Kokkos::RangePolicy(0,nmax),*this); bo_cut_bond = api->control->bg_cut; - atomKK->sync(execution_space,TAG_MASK); + atomKK->sync(execution_space,TAG_MASK|MASK_MASK); tag = atomKK->k_tag.view(); + mask = atomKK->k_mask.view(); const int inum = list->inum; NeighListKokkos* k_list = static_cast*>(list); d_ilist = k_list->d_ilist; numbonds = 0; - PairReaxKokkosFindBondFunctor find_bond_functor(this); + PairReaxKokkosFindBondFunctor find_bond_functor(this, groupbit); Kokkos::parallel_reduce(inum,find_bond_functor,numbonds); copymode = 0; } @@ -4194,24 +4195,28 @@ void PairReaxFFKokkos::operator()(TagPairReaxFindBondZero, const int template KOKKOS_INLINE_FUNCTION -void PairReaxFFKokkos::calculate_find_bond_item(int ii, int &numbonds) const +void PairReaxFFKokkos::calculate_find_bond_item(int ii, int &numbonds, int groupbit) const { const int i = d_ilist[ii]; int nj = 0; - const int j_start = d_bo_first[i]; - const int j_end = j_start + d_bo_num[i]; - for (int jj = j_start; jj < j_end; jj++) { - int j = d_bo_list[jj]; - j &= NEIGHMASK; - const tagint jtag = tag[j]; - const int j_index = jj - j_start; - double bo_tmp = d_BO(i,j_index); + if (mask[i] & groupbit) { + const int j_start = d_bo_first[i]; + const int j_end = j_start + d_bo_num[i]; + for (int jj = j_start; jj < j_end; jj++) { + int j = d_bo_list[jj]; + j &= NEIGHMASK; + if (mask[j] & groupbit) { + const tagint jtag = tag[j]; + const int j_index = jj - j_start; + double bo_tmp = d_BO(i,j_index); - if (bo_tmp > bo_cut_bond) { - d_neighid(i,nj) = jtag; - d_abo(i,nj) = bo_tmp; - nj++; + if (bo_tmp > bo_cut_bond) { + d_neighid(i,nj) = jtag; + d_abo(i,nj) = bo_tmp; + nj++; + } + } } } d_numneigh_bonds[i] = nj; @@ -4247,6 +4252,36 @@ void PairReaxFFKokkos::PackBondBuffer(DAT::tdual_ffloat_1d k_buf, in nbuf_local = k_nbuf_local.h_view(); } +/* ---------------------------------------------------------------------- */ + +template +void PairReaxFFKokkos::PackReducedBondBuffer(DAT::tdual_ffloat_1d k_buf, int &nbuf_local, bool store_bonds) +{ + d_buf = k_buf.view(); + k_params_sing.template sync(); + + copymode = 1; + nlocal = atomKK->nlocal; + if (store_bonds) { + PairReaxKokkosPackReducedBondBufferFunctor pack_bond_buffer_functor(this); + Kokkos::parallel_scan(nlocal,pack_bond_buffer_functor); + } else { + PairReaxKokkosPackReducedBondBufferFunctor pack_bond_buffer_functor(this); + Kokkos::parallel_scan(nlocal,pack_bond_buffer_functor); + } + + copymode = 0; + + k_buf.modify(); + k_nbuf_local.modify(); + + k_buf.sync(); + k_nbuf_local.sync(); + nbuf_local = k_nbuf_local.h_view(); +} + +/* ---------------------------------------------------------------------- */ + template KOKKOS_INLINE_FUNCTION void PairReaxFFKokkos::pack_bond_buffer_item(int i, int &j, const bool &final) const @@ -4288,6 +4323,42 @@ void PairReaxFFKokkos::pack_bond_buffer_item(int i, int &j, const bo k_nbuf_local.view()() = j - 1; } +template +template +KOKKOS_INLINE_FUNCTION +void PairReaxFFKokkos::pack_reduced_bond_buffer_item(int i, int &j, const bool &final) const +{ + const int numbonds = d_numneigh_bonds[i]; + if (final) { + d_buf[j] = d_total_bo[i]; + d_buf[j+1] = paramssing(type[i]).nlp_opt - d_Delta_lp[i]; + d_buf[j+2] = numbonds; + } + + j += 3; + + if constexpr(STORE_BONDS) { + if (final) { + for (int k = 0; k < numbonds; ++k) { + d_buf[j+k] = d_neighid(i,k); + } + } + + j += numbonds; + + if (final) { + for (int k = 0; k < numbonds; k++) { + d_buf[j+k] = d_abo(i,k); + } + } + + j += numbonds; + } + + if (final && i == nlocal-1) + k_nbuf_local.view()() = j - 1; +} + /* ---------------------------------------------------------------------- */ template diff --git a/src/KOKKOS/pair_reaxff_kokkos.h b/src/KOKKOS/pair_reaxff_kokkos.h index 421d704d03..fba7c03ec4 100644 --- a/src/KOKKOS/pair_reaxff_kokkos.h +++ b/src/KOKKOS/pair_reaxff_kokkos.h @@ -130,8 +130,9 @@ class PairReaxFFKokkos : public PairReaxFF { void compute(int, int); void init_style(); double memory_usage(); - void FindBond(int &); + void FindBond(int &, int groupbit = 1); void PackBondBuffer(DAT::tdual_ffloat_1d, int &); + void PackReducedBondBuffer(DAT::tdual_ffloat_1d, int &, bool); void FindBondSpecies(); template @@ -284,11 +285,15 @@ class PairReaxFFKokkos : public PairReaxFF { void operator()(TagPairReaxFindBondZero, const int&) const; KOKKOS_INLINE_FUNCTION - void calculate_find_bond_item(int, int&) const; + void calculate_find_bond_item(int, int&, int) const; KOKKOS_INLINE_FUNCTION void pack_bond_buffer_item(int, int&, const bool&) const; + template + KOKKOS_INLINE_FUNCTION + void pack_reduced_bond_buffer_item(int, int&, const bool&) const; + KOKKOS_INLINE_FUNCTION void operator()(TagPairReaxFindBondSpeciesZero, const int&) const; @@ -409,6 +414,7 @@ class PairReaxFFKokkos : public PairReaxFF { typename AT::t_f_array f; typename AT::t_int_1d_randomread type; typename AT::t_tagint_1d_randomread tag; + typename AT::t_int_1d_randomread mask; typename AT::t_float_1d_randomread q; typename AT::t_tagint_1d_randomread molecule; @@ -518,8 +524,9 @@ template struct PairReaxKokkosFindBondFunctor { typedef DeviceType device_type; typedef int value_type; + int groupbit; PairReaxFFKokkos c; - PairReaxKokkosFindBondFunctor(PairReaxFFKokkos* c_ptr):c(*c_ptr) {}; + PairReaxKokkosFindBondFunctor(PairReaxFFKokkos* c_ptr, int groupbit):c(*c_ptr),groupbit(groupbit) {}; KOKKOS_INLINE_FUNCTION void join(int &dst, @@ -529,7 +536,7 @@ struct PairReaxKokkosFindBondFunctor { KOKKOS_INLINE_FUNCTION void operator()(const int ii, int &numbonds) const { - c.calculate_find_bond_item(ii,numbonds); + c.calculate_find_bond_item(ii,numbonds,groupbit); } }; @@ -546,6 +553,19 @@ struct PairReaxKokkosPackBondBufferFunctor { } }; +template +struct PairReaxKokkosPackReducedBondBufferFunctor { + typedef DeviceType device_type; + typedef int value_type; + PairReaxFFKokkos c; + PairReaxKokkosPackReducedBondBufferFunctor(PairReaxFFKokkos* c_ptr):c(*c_ptr) {}; + + KOKKOS_INLINE_FUNCTION + void operator()(const int ii, int &j, const bool &final) const { + c.template pack_reduced_bond_buffer_item(ii,j,final); + } +}; + } #endif diff --git a/src/REAXFF/compute_reaxff_atom.cpp b/src/REAXFF/compute_reaxff_atom.cpp new file mode 100644 index 0000000000..1b54f3b82a --- /dev/null +++ b/src/REAXFF/compute_reaxff_atom.cpp @@ -0,0 +1,254 @@ +// clang-format off +/* ---------------------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + 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. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Richard Berger (LANL) +------------------------------------------------------------------------- */ + +#include "compute_reaxff_atom.h" +#include "atom.h" +#include "molecule.h" +#include "update.h" +#include "force.h" +#include "memory.h" +#include "error.h" +#include "neigh_list.h" + +#include "pair_reaxff.h" +#include "reaxff_api.h" + +using namespace LAMMPS_NS; +using namespace ReaxFF; + +/* ---------------------------------------------------------------------- */ + +ComputeReaxFFAtom::ComputeReaxFFAtom(LAMMPS *lmp, int narg, char **arg) : + Compute(lmp, narg, arg), + abo(nullptr), neighid(nullptr), bondcount(nullptr), reaxff(nullptr) +{ + if (atom->tag_consecutive() == 0) + error->all(FLERR,"Atom IDs must be consecutive for compute reaxff/atom"); + + peratom_flag = 1; + + // initialize output + + nlocal = -1; + prev_nbonds = -1; + + size_peratom_cols = 3; + + size_local_rows = 0; + size_local_cols = 3; + + invoked_bonds = -1; + + store_bonds = false; + nsub = 0; + + int iarg = 3; + while (iarg narg) utils::missing_cmd_args(FLERR, "compute reaxff/atom pair", error); + ++iarg; + + if (isdigit(arg[iarg][0])) { + nsub = utils::inumeric(FLERR, arg[iarg], false, lmp); + ++iarg; + if (nsub > 0) continue; + } + error->all(FLERR, "Illegal compute reaxff/atom command"); + } else if (strcmp(arg[iarg], "bonds") == 0) { + if (iarg+2 > narg) utils::missing_cmd_args(FLERR, "compute reaxff/atom bonds", error); + store_bonds = utils::logical(FLERR, arg[iarg+1], false, lmp); + iarg += 2; + } else error->all(FLERR,"Illegal compute reaxff/atom command"); + } + + local_flag = store_bonds; +} + +/* ---------------------------------------------------------------------- */ + +ComputeReaxFFAtom::~ComputeReaxFFAtom() +{ + memory->destroy(array_local); + memory->destroy(array_atom); + memory->destroy(abo); + memory->destroy(neighid); + memory->destroy(bondcount); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeReaxFFAtom::init() +{ + if (lmp->suffix_enable) { + if (lmp->suffix) + reaxff = dynamic_cast(force->pair_match(fmt::format("^reax../{}", lmp->suffix), 0, nsub)); + if (!reaxff && lmp->suffix2) + reaxff = dynamic_cast(force->pair_match(fmt::format("^reax../{}", lmp->suffix2), 0, nsub)); + } + + if (!reaxff) reaxff = dynamic_cast(force->pair_match("^reax..", 0, nsub)); + + if (!reaxff) error->all(FLERR,"Cannot use compute reaxff/atom without " + "pair_style reaxff or reaxff/omp"); + + if (reaxff->kokkosable && !kokkosable) + error->all(FLERR,"Cannot use compute reaxff/atom with pair_style reaxff/kk. Use reaxff/atom/kk."); +} + +/* ---------------------------------------------------------------------- */ + +int ComputeReaxFFAtom::FindBond() +{ + int *ilist, i, ii, inum; + int j, pj, nj; + tagint jtag; + double bo_tmp,bo_cut; + + inum = reaxff->list->inum; + ilist = reaxff->list->ilist; + bond_data *bo_ij; + bo_cut = reaxff->api->control->bg_cut; + + tagint *tag = atom->tag; + int * mask = atom->mask; + int numbonds = 0; + + for (ii = 0; ii < inum; ii++) { + i = ilist[ii]; + if (mask[i] & groupbit) { + nj = 0; + + for (pj = Start_Index(i, reaxff->api->lists); pj < End_Index(i, reaxff->api->lists); ++pj) { + bo_ij = &(reaxff->api->lists->select.bond_list[pj]); + j = bo_ij->nbr; + if (mask[j] & groupbit) { + jtag = tag[j]; + bo_tmp = bo_ij->bo_data.BO; + + if (bo_tmp > bo_cut) { + if (store_bonds) { + neighid[i][nj] = jtag; + abo[i][nj] = bo_tmp; + } + nj++; + } + } + } + bondcount[i] = nj; + numbonds += nj; + } + } + return numbonds; +} + +/* ---------------------------------------------------------------------- */ + +void ComputeReaxFFAtom::compute_bonds() +{ + invoked_bonds = update->ntimestep; + + if (atom->nlocal > nlocal) { + memory->destroy(abo); + memory->destroy(neighid); + memory->destroy(bondcount); + memory->destroy(array_atom); + nlocal = atom->nlocal; + if (store_bonds) { + memory->create(abo, nlocal, MAXREAXBOND, "reaxff/atom:abo"); + memory->create(neighid, nlocal, MAXREAXBOND, "reaxff/atom:neighid"); + } + memory->create(bondcount, nlocal, "reaxff/atom:bondcount"); + memory->create(array_atom, nlocal, 3, "reaxff/atom:array_atom"); + } + + for (int i = 0; i < nlocal; i++) { + bondcount[i] = 0; + for (int j = 0; store_bonds && j < MAXREAXBOND; j++) { + neighid[i][j] = 0; + abo[i][j] = 0.0; + } + } + + nbonds = FindBond(); +} + +/* ---------------------------------------------------------------------- */ + +void ComputeReaxFFAtom::compute_local() +{ + invoked_local = update->ntimestep; + + if (invoked_bonds < update->ntimestep) + compute_bonds(); + + if (nbonds > prev_nbonds) { + // grow array_local + memory->destroy(array_local); + memory->create(array_local, nbonds, 3, "reaxff/atom:array_local"); + prev_nbonds = nbonds; + } + + size_local_rows = nbonds; + auto tag = atom->tag; + + int b = 0; + + for (int i = 0; i < nlocal; ++i) { + const int numbonds = bondcount[i]; + + for (int k = 0; k < numbonds; k++) { + auto bond = array_local[b++]; + bond[0] = tag[i]; + bond[1] = neighid[i][k]; + bond[2] = abo[i][k]; + } + } +} + +/* ---------------------------------------------------------------------- */ + +void ComputeReaxFFAtom::compute_peratom() +{ + invoked_peratom = update->ntimestep; + + if (invoked_bonds < update->ntimestep) { + compute_bonds(); + } + + for (int i = 0; i < nlocal; ++i) { + auto ptr = array_atom[i]; + ptr[0] = reaxff->api->workspace->total_bond_order[i]; + ptr[1] = reaxff->api->workspace->nlp[i]; + ptr[2] = bondcount[i]; + } +} + +/* ---------------------------------------------------------------------- + memory usage of local data +------------------------------------------------------------------------- */ + +double ComputeReaxFFAtom::memory_usage() +{ + double bytes = (double)(nlocal*3) * sizeof(double); + bytes += (double)(nlocal) * sizeof(int); + if (store_bonds) { + bytes += (double)(2*nlocal*MAXREAXBOND) * sizeof(double); + bytes += (double)(nbonds*3) * sizeof(double); + } + return bytes; +} diff --git a/src/REAXFF/compute_reaxff_atom.h b/src/REAXFF/compute_reaxff_atom.h new file mode 100644 index 0000000000..1f9aaec1ae --- /dev/null +++ b/src/REAXFF/compute_reaxff_atom.h @@ -0,0 +1,61 @@ +/* -*- c++ -*- ---------------------------------------------------------- + LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator + https://www.lammps.org/, Sandia National Laboratories + LAMMPS development team: developers@lammps.org + + 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. +------------------------------------------------------------------------- */ + +/* ---------------------------------------------------------------------- + Contributing author: Richard Berger (LANL) +------------------------------------------------------------------------- */ + +#ifdef COMPUTE_CLASS +// clang-format off +ComputeStyle(reaxff/atom,ComputeReaxFFAtom); +// clang-format on +#else + +#ifndef LMP_COMPUTE_REAXFF_ATOM_H +#define LMP_COMPUTE_REAXFF_ATOM_H + +#include "compute.h" + +namespace LAMMPS_NS { + +class ComputeReaxFFAtom : public Compute { + public: + ComputeReaxFFAtom(class LAMMPS *, int, char **); + ~ComputeReaxFFAtom() override; + void init() override; + void compute_local() override; + void compute_peratom() override; + virtual void compute_bonds(); + double memory_usage() override; + + protected: + bigint invoked_bonds; // last timestep on which compute_bonds() was invoked + int nlocal; + int nbonds; + int prev_nbonds; + int nsub; + bool store_bonds; + + tagint **neighid; + double **abo; + int *bondcount; + class PairReaxFF *reaxff; + + private: + int FindBond(); +}; + +} // namespace LAMMPS_NS + +#endif +#endif diff --git a/src/YAFF/angle_mm3.cpp b/src/YAFF/angle_mm3.cpp index c75a0d8308..af199f6fe9 100644 --- a/src/YAFF/angle_mm3.cpp +++ b/src/YAFF/angle_mm3.cpp @@ -36,7 +36,10 @@ using namespace MathConst; /* ---------------------------------------------------------------------- */ -AngleMM3::AngleMM3(LAMMPS *lmp) : Angle(lmp) {} +AngleMM3::AngleMM3(LAMMPS *lmp) : Angle(lmp) +{ + born_matrix_enable = 1; +} /* ---------------------------------------------------------------------- */ @@ -284,3 +287,43 @@ double AngleMM3::single(int type, int i1, int i2, int i3) return energy; } + +/* ---------------------------------------------------------------------- */ + +void AngleMM3::born_matrix(int type, int i1, int i2, int i3, double &du, double &du2) +{ + double **x = atom->x; + + double delx1 = x[i1][0] - x[i2][0]; + double dely1 = x[i1][1] - x[i2][1]; + double delz1 = x[i1][2] - x[i2][2]; + domain->minimum_image(delx1,dely1,delz1); + double r1 = sqrt(delx1*delx1 + dely1*dely1 + delz1*delz1); + + double delx2 = x[i3][0] - x[i2][0]; + double dely2 = x[i3][1] - x[i2][1]; + double delz2 = x[i3][2] - x[i2][2]; + domain->minimum_image(delx2,dely2,delz2); + double r2 = sqrt(delx2*delx2 + dely2*dely2 + delz2*delz2); + + double c = delx1*delx2 + dely1*dely2 + delz1*delz2; + c /= r1*r2; + if (c > 1.0) c = 1.0; + if (c < -1.0) c = -1.0; + double theta = acos(c); + + double s = sqrt(1.0 - c*c); + if (s < SMALL) s = SMALL; + s = 1.0/s; + + double dtheta = theta - theta0[type]; + double dtheta2 = dtheta*dtheta; + double dtheta3 = dtheta2*dtheta; + double dtheta4 = dtheta3*dtheta; + double dtheta5 = dtheta4*dtheta; + double df = 2.0 * dtheta - 2.406423 * dtheta2 + 0.735348 * dtheta3 - 0.65832 * dtheta4 + 1.42254 * dtheta5; + double d2f = 2.0 - 4.812846 * dtheta + 2.206044 * dtheta2 - 2.63328 * dtheta3 + 7.1127 * dtheta4; + + du = -k2[type] * df / s; + du2 = k2[type] * (d2f - df * c / s) / (s * s) ; +} diff --git a/src/YAFF/angle_mm3.h b/src/YAFF/angle_mm3.h index 95009a9cf6..22f5bd746c 100644 --- a/src/YAFF/angle_mm3.h +++ b/src/YAFF/angle_mm3.h @@ -35,6 +35,7 @@ class AngleMM3 : public Angle { void read_restart(FILE *) override; void write_data(FILE *) override; double single(int, int, int, int) override; + void born_matrix(int type, int i1, int i2, int i3, double &du, double &du2) override; protected: double *theta0, *k2; diff --git a/src/YAFF/bond_mm3.cpp b/src/YAFF/bond_mm3.cpp index a5ef6fb8bc..31ce2dad3e 100644 --- a/src/YAFF/bond_mm3.cpp +++ b/src/YAFF/bond_mm3.cpp @@ -31,7 +31,10 @@ using namespace LAMMPS_NS; /* ---------------------------------------------------------------------- */ -BondMM3::BondMM3(LAMMPS *lmp) : Bond(lmp) {} +BondMM3::BondMM3(LAMMPS *lmp) : Bond(lmp) +{ + born_matrix_enable = 1; +} /* ---------------------------------------------------------------------- */ @@ -219,3 +222,19 @@ double BondMM3::single(int type, double rsq, else fforce = 0.0; return k2[type]*dr2*(1.0+K3*dr+K4*dr2); } + +/* ---------------------------------------------------------------------- */ + +void BondMM3::born_matrix(int type, double rsq, int /*i*/, int /*j*/, double &du, double &du2) +{ + double r = sqrt(rsq); + double dr = r - r0[type]; + double dr2 = dr * dr; + double dr3 = dr2 * dr; + + double K3 = -2.55 * k2[type] /force->angstrom; + double K4 = 7.0 * k2[type] * 2.55 * 2.55 / (12.0 * force->angstrom * force->angstrom); + + du = 2.0 * k2[type] * dr + 3.0 * K3 * dr2 + 4.0 * K4 * dr3; + du2 = 2.0 * k2[type] + 6.0 * K3 * dr + 12.0 * K4 * dr2; +} diff --git a/src/YAFF/bond_mm3.h b/src/YAFF/bond_mm3.h index 302c4052d0..ea89ac826d 100644 --- a/src/YAFF/bond_mm3.h +++ b/src/YAFF/bond_mm3.h @@ -35,6 +35,7 @@ class BondMM3 : public Bond { void read_restart(FILE *) override; void write_data(FILE *) override; double single(int, double, int, int, double &) override; + void born_matrix(int, double, int, int, double &, double &) override; protected: double *r0, *k2;