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/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/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;