From 36eb11f499d2ebfa896e14bd0caa8792561788ad Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Sat, 20 May 2023 12:59:34 +0300 Subject: [PATCH 01/38] Include method for bond contribution & variables to compute_stress_mop.h --- src/EXTRA-COMPUTE/compute_stress_mop.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop.h b/src/EXTRA-COMPUTE/compute_stress_mop.h index 5357d36371..c9f95c1996 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop.h +++ b/src/EXTRA-COMPUTE/compute_stress_mop.h @@ -38,11 +38,13 @@ class ComputeStressMop : public Compute { private: void compute_pairs(); + void compute_bonds(); int me, nvalues, dir; int *which; double *values_local, *values_global; + double *bond_local, *bond_global; double pos, pos1, dt, nktv2p, ftm2v; double area; class NeighList *list; From 5da65bbd0a158e89b06e959a17cde3af02474d71 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Sat, 20 May 2023 13:06:19 +0300 Subject: [PATCH 02/38] Code for bond contribution to stress/mop --- src/EXTRA-COMPUTE/compute_stress_mop.cpp | 137 ++++++++++++++++++++++- 1 file changed, 132 insertions(+), 5 deletions(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop.cpp b/src/EXTRA-COMPUTE/compute_stress_mop.cpp index 60f2d76e06..7900a62bef 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop.cpp @@ -19,10 +19,13 @@ #include "compute_stress_mop.h" #include "atom.h" +#include "atom_vec.h" +#include "bond.h" #include "domain.h" #include "error.h" #include "force.h" #include "memory.h" +#include "molecule.h" #include "neigh_list.h" #include "neighbor.h" #include "pair.h" @@ -118,12 +121,16 @@ ComputeStressMop::ComputeStressMop(LAMMPS *lmp, int narg, char **arg) : // Initialize some variables values_local = values_global = vector = nullptr; + bond_local = nullptr; + bond_global = nullptr; // this fix produces a global vector memory->create(vector,nvalues,"stress/mop:vector"); memory->create(values_local,nvalues,"stress/mop/spatial:values_local"); memory->create(values_global,nvalues,"stress/mop/spatial:values_global"); + memory->create(bond_local,nvalues,"stress/mop/spatial:bond_local"); + memory->create(bond_global,nvalues,"stress/mop/spatial:bond_global"); size_vector = nvalues; vector_flag = 1; @@ -140,6 +147,8 @@ ComputeStressMop::~ComputeStressMop() memory->destroy(values_local); memory->destroy(values_global); + memory->destroy(bond_local); + memory->destroy(bond_global); memory->destroy(vector); } @@ -185,8 +194,8 @@ void ComputeStressMop::init() //Compute stress/mop only accounts for pair interactions. // issue a warning if any intramolecular potential or Kspace is defined. - if (force->bond!=nullptr) - error->warning(FLERR,"compute stress/mop does not account for bond potentials"); + //if (force->bond!=nullptr) + // error->warning(FLERR,"compute stress/mop does not account for bond potentials"); if (force->angle!=nullptr) error->warning(FLERR,"compute stress/mop does not account for angle potentials"); if (force->dihedral!=nullptr) @@ -224,9 +233,14 @@ void ComputeStressMop::compute_vector() MPI_Allreduce(values_local,values_global,nvalues, MPI_DOUBLE,MPI_SUM,world); - int m; - for (m=0; mx; + tagint *tag = atom->tag; + int *num_bond = atom->num_bond; + tagint **bond_atom = atom->bond_atom; + int **bond_type = atom->bond_type; + int *mask = atom->mask; + + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; + + int nlocal = atom->nlocal; + int newton_bond = force->newton_bond; + int molecular = atom->molecular; + + Bond *bond = force->bond; + + double dx[3] {0}; + double x_bond_1[3] {0}; + double x_bond_2[3] {0}; + double local_contribution[3] {0}; + + // initialization + for (int i {0}; i < nvalues; i++) bond_local[i] = 0.0; + + // loop over all bonded atoms in the current proc + for (atom1 = 0; atom1 < nlocal; atom1++) { + if (!(mask[atom1] & groupbit)) continue; + + if (molecular == 1) + nb = num_bond[atom1]; + else { + if (molindex[atom1] < 0) continue; + imol = molindex[atom1]; + iatom = molatom[atom1]; + nb = onemols[imol]->num_bond[iatom]; + } + + for (i = 0; i < nb; i++) { + if (molecular == 1) { + btype = bond_type[atom1][i]; + atom2 = atom->map(bond_atom[atom1][i]); + } else { + tagprev = tag[atom1] - iatom - 1; + btype = onemols[imol]->bond_type[iatom][i]; + atom2 = atom->map(onemols[imol]->bond_atom[iatom][i] + tagprev); + } + + if (atom2 < 0 || !(mask[atom2] & groupbit)) continue; + if (newton_bond == 0 && tag[atom1] > tag[atom2]) continue; + if (btype <= 0) continue; + + // minimum image of atom1 with respect to the plane of interest + dx[0] = x[atom1][0]; + dx[1] = x[atom1][1]; + dx[2] = x[atom1][2]; + dx[dir] -= pos; + domain->minimum_image(dx[0], dx[1], dx[2]); + x_bond_1[0] = dx[0]; + x_bond_1[1] = dx[1]; + x_bond_1[2] = dx[2]; + x_bond_1[dir] += pos; + + // minimum image of atom2 with respect to atom1 + dx[0] = x[atom2][0] - x_bond_1[0]; + dx[1] = x[atom2][1] - x_bond_1[1]; + dx[2] = x[atom2][2] - x_bond_1[2]; + domain->minimum_image(dx[0], dx[1], dx[2]); + x_bond_2[0] = x_bond_1[0] + dx[0]; + x_bond_2[1] = x_bond_1[1] + dx[1]; + x_bond_2[2] = x_bond_1[2] + dx[2]; + + // check if the bond vector crosses the plane of interest + double tau = (x_bond_1[dir] - pos) / (x_bond_1[dir] - x_bond_2[dir]); + if ((tau <= 1) and (tau >= 0)) + { + //std::cout << "I have found one crossing bond " << tau << std::endl; + dx[0] = x_bond_1[0] - x_bond_2[0]; + dx[1] = x_bond_1[1] - x_bond_2[1]; + dx[2] = x_bond_1[2] - x_bond_2[2]; + rsq = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]; + bond->single(btype, rsq, atom1, atom2, fpair); + + // check the correct contribution with the + or - sign + local_contribution[0] += fpair*dx[0]/area*nktv2p; + local_contribution[1] += fpair*dx[1]/area*nktv2p; + local_contribution[2] += fpair*dx[2]/area*nktv2p; + } + } + } + + // loop over the keywords and if necessary add the bond contribution + int m {0}; + while (m Date: Sat, 20 May 2023 15:12:42 +0300 Subject: [PATCH 03/38] Update compute_stress_mop.cpp --- src/EXTRA-COMPUTE/compute_stress_mop.cpp | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop.cpp b/src/EXTRA-COMPUTE/compute_stress_mop.cpp index 7900a62bef..89ef4e8b2b 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop.cpp @@ -525,16 +525,13 @@ void ComputeStressMop::compute_bonds() // check if the bond vector crosses the plane of interest double tau = (x_bond_1[dir] - pos) / (x_bond_1[dir] - x_bond_2[dir]); - if ((tau <= 1) and (tau >= 0)) - { - //std::cout << "I have found one crossing bond " << tau << std::endl; + if ((tau <= 1) and (tau >= 0)) { dx[0] = x_bond_1[0] - x_bond_2[0]; dx[1] = x_bond_1[1] - x_bond_2[1]; dx[2] = x_bond_1[2] - x_bond_2[2]; rsq = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]; bond->single(btype, rsq, atom1, atom2, fpair); - // check the correct contribution with the + or - sign local_contribution[0] += fpair*dx[0]/area*nktv2p; local_contribution[1] += fpair*dx[1]/area*nktv2p; local_contribution[2] += fpair*dx[2]/area*nktv2p; @@ -552,5 +549,4 @@ void ComputeStressMop::compute_bonds() } m += 3; } - return; } From 1d7a6f813bd623ba1d483da70e13d1c40ef24a43 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Sat, 20 May 2023 15:14:47 +0300 Subject: [PATCH 04/38] Update compute_stress_mop.cpp --- src/EXTRA-COMPUTE/compute_stress_mop.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop.cpp b/src/EXTRA-COMPUTE/compute_stress_mop.cpp index 89ef4e8b2b..8391622242 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop.cpp @@ -525,7 +525,7 @@ void ComputeStressMop::compute_bonds() // check if the bond vector crosses the plane of interest double tau = (x_bond_1[dir] - pos) / (x_bond_1[dir] - x_bond_2[dir]); - if ((tau <= 1) and (tau >= 0)) { + if ((tau <= 1) && (tau >= 0)) { dx[0] = x_bond_1[0] - x_bond_2[0]; dx[1] = x_bond_1[1] - x_bond_2[1]; dx[2] = x_bond_1[2] - x_bond_2[2]; From 0692ed3bd72bfbd9f0524baa8ea45dadefb08407 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Wed, 24 May 2023 11:47:15 +0300 Subject: [PATCH 05/38] @evoyiatzis Include method for angle contribution & variables to compute_stress_mop.h --- src/EXTRA-COMPUTE/compute_stress_mop.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop.h b/src/EXTRA-COMPUTE/compute_stress_mop.h index c9f95c1996..d1074e002a 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop.h +++ b/src/EXTRA-COMPUTE/compute_stress_mop.h @@ -39,12 +39,14 @@ class ComputeStressMop : public Compute { private: void compute_pairs(); void compute_bonds(); + void compute_angles(); int me, nvalues, dir; int *which; double *values_local, *values_global; double *bond_local, *bond_global; + double *angle_local, *angle_global; double pos, pos1, dt, nktv2p, ftm2v; double area; class NeighList *list; From df708a67a5f56f6ac64bd03ac2f6b8f3033f1cb4 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Wed, 24 May 2023 11:55:08 +0300 Subject: [PATCH 06/38] Code for angle contribution to stress/mop --- src/EXTRA-COMPUTE/compute_stress_mop.cpp | 186 ++++++++++++++++++++++- 1 file changed, 184 insertions(+), 2 deletions(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop.cpp b/src/EXTRA-COMPUTE/compute_stress_mop.cpp index 8391622242..139c5b5903 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop.cpp @@ -18,6 +18,7 @@ #include "compute_stress_mop.h" +#include "angle.h" #include "atom.h" #include "atom_vec.h" #include "bond.h" @@ -123,6 +124,8 @@ ComputeStressMop::ComputeStressMop(LAMMPS *lmp, int narg, char **arg) : values_local = values_global = vector = nullptr; bond_local = nullptr; bond_global = nullptr; + angle_local = nullptr; + angle_global = nullptr; // this fix produces a global vector @@ -131,6 +134,8 @@ ComputeStressMop::ComputeStressMop(LAMMPS *lmp, int narg, char **arg) : memory->create(values_global,nvalues,"stress/mop/spatial:values_global"); memory->create(bond_local,nvalues,"stress/mop/spatial:bond_local"); memory->create(bond_global,nvalues,"stress/mop/spatial:bond_global"); + memory->create(angle_local,nvalues,"stress/mop/spatial:angle_local"); + memory->create(angle_global,nvalues,"stress/mop/spatial:angle_global"); size_vector = nvalues; vector_flag = 1; @@ -149,6 +154,8 @@ ComputeStressMop::~ComputeStressMop() memory->destroy(values_global); memory->destroy(bond_local); memory->destroy(bond_global); + memory->destroy(angle_local); + memory->destroy(angle_global); memory->destroy(vector); } @@ -197,7 +204,8 @@ void ComputeStressMop::init() //if (force->bond!=nullptr) // error->warning(FLERR,"compute stress/mop does not account for bond potentials"); if (force->angle!=nullptr) - error->warning(FLERR,"compute stress/mop does not account for angle potentials"); + if (force->angle->born_matrix_enable == 0) + error->warning(FLERR,"compute stress/mop does not account for angle potentials"); if (force->dihedral!=nullptr) error->warning(FLERR,"compute stress/mop does not account for dihedral potentials"); if (force->improper!=nullptr) @@ -239,8 +247,14 @@ void ComputeStressMop::compute_vector() // sum bond contribution over all procs MPI_Allreduce(bond_local,bond_global,nvalues,MPI_DOUBLE,MPI_SUM,world); + //Compute angle contribution on separate procs + compute_angles(); + + // sum angle contribution over all procs + MPI_Allreduce(angle_local,angle_global,nvalues,MPI_DOUBLE,MPI_SUM,world); + for (int m=0; mx; + 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}; + double dx_left[3] {0}; + double dx_right[3] {0}; + double x_angle_left[3] {0}; + double x_angle_middle[3] {0}; + double x_angle_right[3] {0}; + double dcos_theta[3] {0}; + double local_contribution[3] {0}; + + // initialization + for (int i {0}; i < nvalues; i++) angle_local[i] = 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; + + // minimum image of atom1 with respect to the plane of interest + dx[0] = x[atom1][0]; + dx[1] = x[atom1][1]; + dx[2] = x[atom1][2]; + dx[dir] -= pos; + domain->minimum_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) + { + dcos_theta[0] = (dx_right[0]*cos_theta/r2 - dx_left[0]/r1)/r2; + dcos_theta[1] = (dx_right[1]*cos_theta/r2 - dx_left[1]/r1)/r2; + dcos_theta[2] = (dx_right[1]*cos_theta/r2 - dx_left[2]/r1)/r2; + } + + // only left bond crossing the plane + if (!right_cross && left_cross) + { + dcos_theta[0] = -(dx_left[0]*cos_theta/r1 - dx_right[0]/r2)/r1; + dcos_theta[1] = -(dx_left[1]*cos_theta/r1 - dx_right[1]/r2)/r1; + dcos_theta[2] = -(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 + dcos_theta[0] = -(dx_right[0]*cos_theta/r2 - dx_left[0]/r1)/r2; + dcos_theta[1] = -(dx_right[1]*cos_theta/r2 - dx_left[1]/r1)/r2; + dcos_theta[2] = -(dx_right[1]*cos_theta/r2 - dx_left[2]/r1)/r2; + + // due to left bond + dcos_theta[0] += (dx_left[0]*cos_theta/r1 - dx_right[0]/r2)/r1; + dcos_theta[1] += (dx_left[1]*cos_theta/r1 - dx_right[1]/r2)/r1; + dcos_theta[2] += (dx_left[2]*cos_theta/r1 - dx_right[2]/r2)/r1; + } + + // final contribution of the given angle term + local_contribution[0] += duang*dcos_theta[0]/area*nktv2p; + local_contribution[1] += duang*dcos_theta[1]/area*nktv2p; + local_contribution[2] += duang*dcos_theta[2]/area*nktv2p; + } + } + // loop over the keywords and if necessary add the angle contribution + int m {0}; + while (m Date: Wed, 24 May 2023 16:24:33 +0300 Subject: [PATCH 07/38] fixing bug with sign for angle contribution in compute_stress_mop.cpp --- src/EXTRA-COMPUTE/compute_stress_mop.cpp | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop.cpp b/src/EXTRA-COMPUTE/compute_stress_mop.cpp index 139c5b5903..ff7c6d9e1b 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop.cpp @@ -705,14 +705,14 @@ void ComputeStressMop::compute_angles() if (right_cross && left_cross) { // due to right bond - dcos_theta[0] = -(dx_right[0]*cos_theta/r2 - dx_left[0]/r1)/r2; - dcos_theta[1] = -(dx_right[1]*cos_theta/r2 - dx_left[1]/r1)/r2; - dcos_theta[2] = -(dx_right[1]*cos_theta/r2 - dx_left[2]/r1)/r2; + dcos_theta[0] = (dx_right[0]*cos_theta/r2 - dx_left[0]/r1)/r2; + dcos_theta[1] = (dx_right[1]*cos_theta/r2 - dx_left[1]/r1)/r2; + dcos_theta[2] = (dx_right[1]*cos_theta/r2 - dx_left[2]/r1)/r2; // due to left bond - dcos_theta[0] += (dx_left[0]*cos_theta/r1 - dx_right[0]/r2)/r1; - dcos_theta[1] += (dx_left[1]*cos_theta/r1 - dx_right[1]/r2)/r1; - dcos_theta[2] += (dx_left[2]*cos_theta/r1 - dx_right[2]/r2)/r1; + dcos_theta[0] -= (dx_left[0]*cos_theta/r1 - dx_right[0]/r2)/r1; + dcos_theta[1] -= (dx_left[1]*cos_theta/r1 - dx_right[1]/r2)/r1; + dcos_theta[2] -= (dx_left[2]*cos_theta/r1 - dx_right[2]/r2)/r1; } // final contribution of the given angle term From c7c8b065a259979e7fcd3b6ccc777d3068a68d52 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Wed, 24 May 2023 16:49:35 +0300 Subject: [PATCH 08/38] fixing bug with sign issue for bond contribution in compute_stress_mop.cpp --- src/EXTRA-COMPUTE/compute_stress_mop.cpp | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop.cpp b/src/EXTRA-COMPUTE/compute_stress_mop.cpp index ff7c6d9e1b..2454f93b7a 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop.cpp @@ -546,9 +546,10 @@ void ComputeStressMop::compute_bonds() rsq = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]; bond->single(btype, rsq, atom1, atom2, fpair); - local_contribution[0] += fpair*dx[0]/area*nktv2p; - local_contribution[1] += fpair*dx[1]/area*nktv2p; - local_contribution[2] += fpair*dx[2]/area*nktv2p; + double sgn = copysign(1.0, x_bond_1[dir] - pos); + local_contribution[0] += sgn*fpair*dx[0]/area*nktv2p; + local_contribution[1] += sgn*fpair*dx[1]/area*nktv2p; + local_contribution[2] += sgn*fpair*dx[2]/area*nktv2p; } } } From 9ee40cceef79043d8c6c8c41e41093e227f97ff6 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Wed, 24 May 2023 17:01:13 +0300 Subject: [PATCH 09/38] fixing indexing issue and more sign problems for angle contributions --- src/EXTRA-COMPUTE/compute_stress_mop.cpp | 27 +++++++++++++----------- 1 file changed, 15 insertions(+), 12 deletions(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop.cpp b/src/EXTRA-COMPUTE/compute_stress_mop.cpp index 2454f93b7a..63d671eee7 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop.cpp @@ -689,31 +689,34 @@ void ComputeStressMop::compute_angles() // only right bond crossing the plane if (right_cross && !left_cross) { - dcos_theta[0] = (dx_right[0]*cos_theta/r2 - dx_left[0]/r1)/r2; - dcos_theta[1] = (dx_right[1]*cos_theta/r2 - dx_left[1]/r1)/r2; - dcos_theta[2] = (dx_right[1]*cos_theta/r2 - dx_left[2]/r1)/r2; + 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) { - dcos_theta[0] = -(dx_left[0]*cos_theta/r1 - dx_right[0]/r2)/r1; - dcos_theta[1] = -(dx_left[1]*cos_theta/r1 - dx_right[1]/r2)/r1; - dcos_theta[2] = -(dx_left[2]*cos_theta/r1 - dx_right[2]/r2)/r1; + 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 - dcos_theta[0] = (dx_right[0]*cos_theta/r2 - dx_left[0]/r1)/r2; - dcos_theta[1] = (dx_right[1]*cos_theta/r2 - dx_left[1]/r1)/r2; - dcos_theta[2] = (dx_right[1]*cos_theta/r2 - dx_left[2]/r1)/r2; + 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] -= (dx_left[0]*cos_theta/r1 - dx_right[0]/r2)/r1; - dcos_theta[1] -= (dx_left[1]*cos_theta/r1 - dx_right[1]/r2)/r1; - dcos_theta[2] -= (dx_left[2]*cos_theta/r1 - dx_right[2]/r2)/r1; + 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 From 8e6615918b0c7414664e17cea86378921d132069 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Thu, 25 May 2023 14:39:03 +0300 Subject: [PATCH 10/38] avoid crashing when no bonds or no angles exist --- src/EXTRA-COMPUTE/compute_stress_mop.h | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop.h b/src/EXTRA-COMPUTE/compute_stress_mop.h index d1074e002a..791e18d624 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop.h +++ b/src/EXTRA-COMPUTE/compute_stress_mop.h @@ -44,6 +44,8 @@ class ComputeStressMop : public Compute { int me, nvalues, dir; int *which; + int bondflag, angleflag; + double *values_local, *values_global; double *bond_local, *bond_global; double *angle_local, *angle_global; From f26f397e087c183facd9e1cb108fe7f1b777a221 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Thu, 25 May 2023 14:47:20 +0300 Subject: [PATCH 11/38] avoid crashing when there are no bonds or no angles --- src/EXTRA-COMPUTE/compute_stress_mop.cpp | 27 ++++++++++++++++++------ 1 file changed, 20 insertions(+), 7 deletions(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop.cpp b/src/EXTRA-COMPUTE/compute_stress_mop.cpp index 63d671eee7..e7a602e583 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop.cpp @@ -51,6 +51,9 @@ ComputeStressMop::ComputeStressMop(LAMMPS *lmp, int narg, char **arg) : MPI_Comm_rank(world,&me); + bondflag = 0; + angleflag = 0; + // set compute mode and direction of plane(s) for pressure calculation if (strcmp(arg[3],"x")==0) { @@ -201,11 +204,13 @@ void ComputeStressMop::init() //Compute stress/mop only accounts for pair interactions. // issue a warning if any intramolecular potential or Kspace is defined. - //if (force->bond!=nullptr) - // error->warning(FLERR,"compute stress/mop does not account for bond potentials"); + if (force->bond!=nullptr) bondflag = 1; if (force->angle!=nullptr) - if (force->angle->born_matrix_enable == 0) + if (force->angle->born_matrix_enable == 0) { error->warning(FLERR,"compute stress/mop does not account for angle potentials"); + } else { + angleflag = 1; + } if (force->dihedral!=nullptr) error->warning(FLERR,"compute stress/mop does not account for dihedral potentials"); if (force->improper!=nullptr) @@ -241,14 +246,22 @@ void ComputeStressMop::compute_vector() MPI_Allreduce(values_local,values_global,nvalues, MPI_DOUBLE,MPI_SUM,world); - //Compute bond contribution on separate procs - compute_bonds(); + if (bondflag) { + //Compute bond contribution on separate procs + compute_bonds(); + } else { + for (int i=0; i Date: Sat, 27 May 2023 14:33:46 +0300 Subject: [PATCH 12/38] Updating unit test for mop to reflect the contribution from bonds --- unittest/commands/test_compute_global.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/unittest/commands/test_compute_global.cpp b/unittest/commands/test_compute_global.cpp index 6c365c8c2b..66a24f049f 100644 --- a/unittest/commands/test_compute_global.cpp +++ b/unittest/commands/test_compute_global.cpp @@ -225,9 +225,9 @@ TEST_F(ComputeGlobalTest, Geometry) EXPECT_DOUBLE_EQ(mom2[0], -0.022332069630161717); EXPECT_DOUBLE_EQ(mom2[1], -0.056896553865696115); EXPECT_DOUBLE_EQ(mom2[2], 0.069179891052881484); - EXPECT_DOUBLE_EQ(mop1[0], 3522311.3572200728); - EXPECT_DOUBLE_EQ(mop1[1], 2871104.9055934539); - EXPECT_DOUBLE_EQ(mop1[2], -4136077.5224247416); + EXPECT_DOUBLE_EQ(mop1[0], 3536584.0880458541); + EXPECT_DOUBLE_EQ(mop1[1], 2887485.033995091); + EXPECT_DOUBLE_EQ(mop1[2], -4154145.8952306858); EXPECT_DOUBLE_EQ(mop2[0][0], -8.0869239999999998); EXPECT_DOUBLE_EQ(mop2[0][1], 0.0); EXPECT_DOUBLE_EQ(mop2[0][2], 0.0); From b28ee36f000af8e90a70390b72f6d132312de43c Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Sat, 27 May 2023 14:50:31 +0300 Subject: [PATCH 13/38] update documentation for compute stress/mop --- doc/src/compute_stress_mop.rst | 9 +++++++-- 1 file changed, 7 insertions(+), 2 deletions(-) diff --git a/doc/src/compute_stress_mop.rst b/doc/src/compute_stress_mop.rst index 4ad2261bb0..a8a3bc5660 100644 --- a/doc/src/compute_stress_mop.rst +++ b/doc/src/compute_stress_mop.rst @@ -117,8 +117,13 @@ 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. In -particular, it does not work with more than two-body pair interactions, -intra-molecular interactions, and long range (kspace) interactions. +particular, compute *stress/mop/profile* does not work with more than +two-body pair interactions, intra-molecular interactions, and long range +(kspace) 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. Related commands """""""""""""""" From 86743bc0a6fa09ffd630d4d9220c76f60823d18d Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Mon, 29 May 2023 10:59:18 +0300 Subject: [PATCH 14/38] Update compute_stress_mop.cpp The angle was computed using the dot product of the vectors x[atom2] - x[atom1] and x[atom3] - x[atom2]. This is not consistent with the lammps convention where the angle is computed using the dot product between x[atom1]-x[atom2] and x[atom3]-x[atom2]. --- src/EXTRA-COMPUTE/compute_stress_mop.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop.cpp b/src/EXTRA-COMPUTE/compute_stress_mop.cpp index e7a602e583..d512ce4810 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop.cpp @@ -692,7 +692,7 @@ void ComputeStressMop::compute_angles() // 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); + 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; From 3b38145d91b130b72a6770992efa35b282f0cfee Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Mon, 29 May 2023 16:34:44 +0300 Subject: [PATCH 15/38] Update compute_stress_mop.cpp Fixing sign issues because I was considering the theta angle to be formed by vectors x[atom2] - x[atom1] & x[atom3] - x[atom2] instead of x[atom1] - x[atom2] & x[atom3] - x[atom2] as done in lammps --- src/EXTRA-COMPUTE/compute_stress_mop.cpp | 24 ++++++++++++------------ 1 file changed, 12 insertions(+), 12 deletions(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop.cpp b/src/EXTRA-COMPUTE/compute_stress_mop.cpp index d512ce4810..a5f5969317 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop.cpp @@ -703,18 +703,18 @@ void ComputeStressMop::compute_angles() 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; + 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; + 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 @@ -722,14 +722,14 @@ void ComputeStressMop::compute_angles() { // 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; + 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; + 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 From b3e9efcb50ded57acff46cc52305c59a0adb04a8 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Mon, 29 May 2023 17:55:24 +0300 Subject: [PATCH 16/38] Use system periodicity to find an equivalent position of the plane --- src/EXTRA-COMPUTE/compute_stress_mop.cpp | 16 ++++++++++++---- 1 file changed, 12 insertions(+), 4 deletions(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop.cpp b/src/EXTRA-COMPUTE/compute_stress_mop.cpp index a5f5969317..9249e4ce2d 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop.cpp @@ -74,6 +74,18 @@ ComputeStressMop::ComputeStressMop(LAMMPS *lmp, int narg, char **arg) : pos = 0.5*(domain->boxlo[dir]+domain->boxhi[dir]); } else pos = utils::numeric(FLERR,arg[4],false,lmp); + // plane inside the box + if (pos >domain->boxhi[dir] || pos boxlo[dir]) { + error->warning(FLERR,"The specified initial plane lies outside of the simulation box"); + double dx[3] {0}; + dx[dir] = pos - 0.5*(domain->boxhi[dir] + domain->boxlo[dir]); + domain->minimum_image(dx[0], dx[1], dx[2]); + pos = 0.5*(domain->boxhi[dir] + domain->boxlo[dir]) + dx[dir]; + + if (pos >domain->boxhi[dir] || pos boxlo[dir]) + error->all(FLERR, "Plane for compute stress/mop is out of bounds"); + } + if (pos < (domain->boxlo[dir]+domain->prd_half[dir])) { pos1 = pos + domain->prd[dir]; } else { @@ -117,10 +129,6 @@ ComputeStressMop::ComputeStressMop(LAMMPS *lmp, int narg, char **arg) : // orthogonal simulation box if (domain->triclinic != 0) error->all(FLERR, "Compute stress/mop incompatible with triclinic simulation box"); - // plane inside the box - if (pos >domain->boxhi[dir] || pos boxlo[dir]) - error->all(FLERR, "Plane for compute stress/mop is out of bounds"); - // Initialize some variables From ea6ece510edc25cc8f381d7d67e546f62bbd8058 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Thu, 1 Jun 2023 10:22:01 +0300 Subject: [PATCH 17/38] turning warning into errors for unsupported styles --- src/EXTRA-COMPUTE/compute_stress_mop.cpp | 13 ++++++------- 1 file changed, 6 insertions(+), 7 deletions(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop.cpp b/src/EXTRA-COMPUTE/compute_stress_mop.cpp index 9249e4ce2d..0199a49c01 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop.cpp @@ -205,26 +205,25 @@ void ComputeStressMop::init() if (force->pair->single_enable == 0) error->all(FLERR,"Pair style does not support compute stress/mop"); - // Warnings + // Errors if (me==0) { - //Compute stress/mop only accounts for pair interactions. - // issue a warning if any intramolecular potential or Kspace is defined. + // issue an error for unimplemented intramolecular potentials or Kspace. if (force->bond!=nullptr) bondflag = 1; if (force->angle!=nullptr) if (force->angle->born_matrix_enable == 0) { - error->warning(FLERR,"compute stress/mop does not account for angle potentials"); + error->all(FLERR,"compute stress/mop does not account for angle potentials"); } else { angleflag = 1; } if (force->dihedral!=nullptr) - error->warning(FLERR,"compute stress/mop does not account for dihedral potentials"); + error->all(FLERR,"compute stress/mop does not account for dihedral potentials"); if (force->improper!=nullptr) - error->warning(FLERR,"compute stress/mop does not account for improper potentials"); + error->all(FLERR,"compute stress/mop does not account for improper potentials"); if (force->kspace!=nullptr) - error->warning(FLERR,"compute stress/mop does not account for kspace contributions"); + error->all(FLERR,"compute stress/mop does not account for kspace contributions"); } // need an occasional half neighbor list From b01db47b2d83a5d4656ef279c318c717cfefe192 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Thu, 1 Jun 2023 10:33:50 +0300 Subject: [PATCH 18/38] consistency in issuing errors between mop and mop/profile --- src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp index 709f109cb1..ffcc712548 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp @@ -178,23 +178,23 @@ void ComputeStressMopProfile::init() if (force->pair->single_enable == 0) error->all(FLERR,"Pair style does not support compute stress/mop/profile"); - // Warnings + // Errors if (me==0) { //Compute stress/mop/profile only accounts for pair interactions. - // issue a warning if any intramolecular potential or Kspace is defined. + // issue an error if any intramolecular potential or Kspace is defined. if (force->bond!=nullptr) - error->warning(FLERR,"compute stress/mop/profile does not account for bond potentials"); + error->all(FLERR,"compute stress/mop/profile does not account for bond potentials"); if (force->angle!=nullptr) - error->warning(FLERR,"compute stress/mop/profile does not account for angle potentials"); + error->all(FLERR,"compute stress/mop/profile does not account for angle potentials"); if (force->dihedral!=nullptr) - error->warning(FLERR,"compute stress/mop/profile does not account for dihedral potentials"); + error->all(FLERR,"compute stress/mop/profile does not account for dihedral potentials"); if (force->improper!=nullptr) - error->warning(FLERR,"compute stress/mop/profile does not account for improper potentials"); + error->all(FLERR,"compute stress/mop/profile does not account for improper potentials"); if (force->kspace!=nullptr) - error->warning(FLERR,"compute stress/mop/profile does not account for kspace contributions"); + error->all(FLERR,"compute stress/mop/profile does not account for kspace contributions"); } // need an occasional half neighbor list From c1cec4565250c5a1b28c850badae0785b140c60d Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Thu, 1 Jun 2023 15:26:43 +0300 Subject: [PATCH 19/38] add keywords to specify contributions to stress/mop --- src/EXTRA-COMPUTE/compute_stress_mop.cpp | 23 +++++++++++++++++++---- 1 file changed, 19 insertions(+), 4 deletions(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop.cpp b/src/EXTRA-COMPUTE/compute_stress_mop.cpp index 0199a49c01..071bdf31ee 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop.cpp @@ -38,7 +38,7 @@ using namespace LAMMPS_NS; enum{X,Y,Z}; -enum{TOTAL,CONF,KIN}; +enum{TOTAL,CONF,KIN,PAIR,BOND,ANGLE}; #define BIG 1000000000 @@ -115,6 +115,21 @@ ComputeStressMop::ComputeStressMop(LAMMPS *lmp, int narg, char **arg) : which[nvalues] = TOTAL; nvalues++; } + } else if (strcmp(arg[iarg],"pair") == 0) { + for (i=0; i<3; i++) { + which[nvalues] = PAIR; + nvalues++; + } + } else if (strcmp(arg[iarg],"bond") == 0) { + for (i=0; i<3; i++) { + which[nvalues] = BOND; + nvalues++; + } + } else if (strcmp(arg[iarg],"angle") == 0) { + for (i=0; i<3; i++) { + which[nvalues] = ANGLE; + nvalues++; + } } else error->all(FLERR, "Illegal compute stress/mop command"); //break; iarg++; @@ -329,7 +344,7 @@ void ComputeStressMop::compute_pairs() m = 0; while (m Date: Thu, 1 Jun 2023 15:34:41 +0300 Subject: [PATCH 20/38] Update compute_stress_mop.rst to reflect the added keywords --- doc/src/compute_stress_mop.rst | 10 ++++++---- 1 file changed, 6 insertions(+), 4 deletions(-) diff --git a/doc/src/compute_stress_mop.rst b/doc/src/compute_stress_mop.rst index a8a3bc5660..b4c15257af 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* (one of more can be specified) +* keywords = *kin* or *conf* or *total* or *pair* or *bond* or *angle* (one or more can be specified) .. parsed-literal:: @@ -64,9 +64,11 @@ atoms have crossed the plane if their positions at times :math:`t-\Delta t` and :math:`t` are one on either side of the plane, and uses the velocity at time :math:`t-\Delta t/2` given by the velocity Verlet algorithm. -Between one and three keywords can be used to indicate which -contributions to the stress must be computed: kinetic stress (kin), -configurational stress (conf), and/or total stress (total). +Between one and six 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 last three keywords are currently available only for the stress/mop command and not the stress/mop/profile. NOTE 1: The configurational stress is computed considering all pairs of atoms where at least one atom belongs to group group-ID. From 08ffd268bff9734217ea14be6d223810fd983991 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Fri, 2 Jun 2023 09:43:58 +0300 Subject: [PATCH 21/38] remove unused symbolic constant --- src/EXTRA-COMPUTE/compute_stress_mop.cpp | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop.cpp b/src/EXTRA-COMPUTE/compute_stress_mop.cpp index 071bdf31ee..3ae758b875 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop.cpp @@ -40,8 +40,6 @@ using namespace LAMMPS_NS; enum{X,Y,Z}; enum{TOTAL,CONF,KIN,PAIR,BOND,ANGLE}; -#define BIG 1000000000 - /* ---------------------------------------------------------------------- */ ComputeStressMop::ComputeStressMop(LAMMPS *lmp, int narg, char **arg) : From 70507462e9e04e8466e9e7c48a56f065235f871d Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Fri, 2 Jun 2023 19:22:28 +0300 Subject: [PATCH 22/38] Include method for bond contribution & variables to compute_stress_mop_profile.h --- src/EXTRA-COMPUTE/compute_stress_mop_profile.h | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop_profile.h b/src/EXTRA-COMPUTE/compute_stress_mop_profile.h index b58f762c12..7bb87fa4c6 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop_profile.h +++ b/src/EXTRA-COMPUTE/compute_stress_mop_profile.h @@ -38,16 +38,20 @@ class ComputeStressMopProfile : public Compute { private: void compute_pairs(); + void compute_bonds(); void setup_bins(); int me, nvalues, dir; int *which; + int bondflag; + int originflag; double origin, delta, offset, invdelta; int nbins; double **coord, **coordp; double **values_local, **values_global; + double **bond_local, **bond_global; double dt, nktv2p, ftm2v; double area; From e7ae02dd2a53b12f8b6de5eabea4d9750d98b085 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Fri, 2 Jun 2023 19:33:15 +0300 Subject: [PATCH 23/38] Code for bond contribution to stress/mop/profile --- .../compute_stress_mop_profile.cpp | 169 +++++++++++++++++- 1 file changed, 162 insertions(+), 7 deletions(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp index ffcc712548..d15067a266 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp @@ -19,8 +19,11 @@ #include "compute_stress_mop_profile.h" #include "atom.h" +#include "atom_vec.h" +#include "bond.h" #include "update.h" #include "domain.h" +#include "molecule.h" #include "neighbor.h" #include "force.h" #include "pair.h" @@ -35,9 +38,7 @@ using namespace LAMMPS_NS; enum{X,Y,Z}; enum{LOWER,CENTER,UPPER,COORD}; -enum{TOTAL,CONF,KIN}; - -#define BIG 1000000000 +enum{TOTAL,CONF,KIN,PAIR,BOND}; /* ---------------------------------------------------------------------- */ @@ -48,6 +49,8 @@ ComputeStressMopProfile::ComputeStressMopProfile(LAMMPS *lmp, int narg, char **a MPI_Comm_rank(world,&me); + bondflag = 0; + // set compute mode and direction of plane(s) for pressure calculation if (strcmp(arg[3],"x")==0) { @@ -92,6 +95,16 @@ ComputeStressMopProfile::ComputeStressMopProfile(LAMMPS *lmp, int narg, char **a which[nvalues] = TOTAL; nvalues++; } + } else if (strcmp(arg[iarg],"pair") == 0) { + for (i=0; i<3; i++) { + which[nvalues] = PAIR; + nvalues++; + } + } else if (strcmp(arg[iarg],"bond") == 0) { + for (i=0; i<3; i++) { + which[nvalues] = PAIR; + nvalues++; + } } else error->all(FLERR, "Illegal compute stress/mop/profile command"); //break; iarg++; @@ -114,6 +127,8 @@ ComputeStressMopProfile::ComputeStressMopProfile(LAMMPS *lmp, int narg, char **a nbins = 0; coord = coordp = nullptr; values_local = values_global = array = nullptr; + bond_local = nullptr; + bond_global = nullptr; // bin setup @@ -140,6 +155,8 @@ ComputeStressMopProfile::~ComputeStressMopProfile() memory->destroy(coordp); memory->destroy(values_local); memory->destroy(values_global); + memory->destroy(bond_local); + memory->destroy(bond_global); memory->destroy(array); } @@ -185,8 +202,8 @@ void ComputeStressMopProfile::init() //Compute stress/mop/profile only accounts for pair interactions. // issue an error if any intramolecular potential or Kspace is defined. - if (force->bond!=nullptr) - error->all(FLERR,"compute stress/mop/profile does not account for bond potentials"); + if (force->bond!=nullptr) bondflag = 1; + if (force->angle!=nullptr) error->all(FLERR,"compute stress/mop/profile does not account for angle potentials"); if (force->dihedral!=nullptr) @@ -225,6 +242,21 @@ void ComputeStressMopProfile::compute_array() MPI_Allreduce(&values_local[0][0],&values_global[0][0],nbins*nvalues, MPI_DOUBLE,MPI_SUM,world); + if (bondflag) { + //Compute bond contribution on separate procs + compute_bonds(); + } else { + for (int m = 0; m < nbins; m++) { + for (int i = 0; i < nvalues; i++) { + bond_local[m][i] = 0.0; + } + } + } + + // sum bond contribution over all procs + MPI_Allreduce(&bond_local[0][0],&bond_global[0][0],nbins*nvalues, + MPI_DOUBLE,MPI_SUM,world); + int ibin,m,mo; for (ibin=0; ibinx; + tagint *tag = atom->tag; + int *num_bond = atom->num_bond; + tagint **bond_atom = atom->bond_atom; + int **bond_type = atom->bond_type; + int *mask = atom->mask; + + int *molindex = atom->molindex; + int *molatom = atom->molatom; + Molecule **onemols = atom->avec->onemols; + + int nlocal = atom->nlocal; + int newton_bond = force->newton_bond; + int molecular = atom->molecular; + + Bond *bond = force->bond; + + double dx[3] {0}; + double x_bond_1[3] {0}; + double x_bond_2[3] {0}; + double local_contribution[nbins][3] {0}; + + // initialization + for (int m {0}; m < nbins; m++) { + for (int i {0}; i < nvalues; i++) { + bond_local[m][i] = 0.0; + } + } + + // loop over all bonded atoms in the current proc + for (atom1 = 0; atom1 < nlocal; atom1++) { + if (!(mask[atom1] & groupbit)) continue; + + if (molecular == 1) + nb = num_bond[atom1]; + else { + if (molindex[atom1] < 0) continue; + imol = molindex[atom1]; + iatom = molatom[atom1]; + nb = onemols[imol]->num_bond[iatom]; + } + + for (i = 0; i < nb; i++) { + if (molecular == 1) { + btype = bond_type[atom1][i]; + atom2 = atom->map(bond_atom[atom1][i]); + } else { + tagprev = tag[atom1] - iatom - 1; + btype = onemols[imol]->bond_type[iatom][i]; + atom2 = atom->map(onemols[imol]->bond_atom[iatom][i] + tagprev); + } + + if (atom2 < 0 || !(mask[atom2] & groupbit)) continue; + if (newton_bond == 0 && tag[atom1] > tag[atom2]) continue; + if (btype <= 0) continue; + + for (int ibin {0}; ibinminimum_image(dx[0], dx[1], dx[2]); + x_bond_1[0] = dx[0]; + x_bond_1[1] = dx[1]; + x_bond_1[2] = dx[2]; + x_bond_1[dir] += pos; + + // minimum image of atom2 with respect to atom1 + dx[0] = x[atom2][0] - x_bond_1[0]; + dx[1] = x[atom2][1] - x_bond_1[1]; + dx[2] = x[atom2][2] - x_bond_1[2]; + domain->minimum_image(dx[0], dx[1], dx[2]); + x_bond_2[0] = x_bond_1[0] + dx[0]; + x_bond_2[1] = x_bond_1[1] + dx[1]; + x_bond_2[2] = x_bond_1[2] + dx[2]; + + // check if the bond vector crosses the plane of interest + double tau = (x_bond_1[dir] - pos) / (x_bond_1[dir] - x_bond_2[dir]); + if ((tau <= 1) && (tau >= 0)) { + dx[0] = x_bond_1[0] - x_bond_2[0]; + dx[1] = x_bond_1[1] - x_bond_2[1]; + dx[2] = x_bond_1[2] - x_bond_2[2]; + rsq = dx[0] * dx[0] + dx[1] * dx[1] + dx[2] * dx[2]; + bond->single(btype, rsq, atom1, atom2, fpair); + + double sgn = copysign(1.0, x_bond_1[dir] - pos); + local_contribution[ibin][0] += sgn*fpair*dx[0]/area*nktv2p; + local_contribution[ibin][1] += sgn*fpair*dx[1]/area*nktv2p; + local_contribution[ibin][2] += sgn*fpair*dx[2]/area*nktv2p; + } + } + } + } + + // loop over the keywords and if necessary add the bond contribution + int m {0}; + while (mcreate(coordp,nbins,1,"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"); // set bin coordinates for (i = 0; i < nbins; i++) { From c30762ca8bfd35bb8554e2024ea57462298fa2e5 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Fri, 2 Jun 2023 19:44:37 +0300 Subject: [PATCH 24/38] Update documentation for compute stress/mop/profile --- doc/src/compute_stress_mop.rst | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/doc/src/compute_stress_mop.rst b/doc/src/compute_stress_mop.rst index b4c15257af..e1b151ebbd 100644 --- a/doc/src/compute_stress_mop.rst +++ b/doc/src/compute_stress_mop.rst @@ -68,7 +68,7 @@ Between one and six 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 last three keywords are currently available only for the stress/mop command and not the stress/mop/profile. +The angle keyword is currently available only for the stress/mop command and not the stress/mop/profile. NOTE 1: The configurational stress is computed considering all pairs of atoms where at least one atom belongs to group group-ID. @@ -120,8 +120,8 @@ 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. In particular, compute *stress/mop/profile* does not work with more than -two-body pair interactions, intra-molecular interactions, and long range -(kspace) interactions. Similarly, compute *stress/mop* does 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() From 28e3a741a8e84d7512e10c34364ac198cf456266 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Fri, 2 Jun 2023 20:02:52 +0300 Subject: [PATCH 25/38] declare local_contribution as pointer in compute_stress_mop_profile.h --- src/EXTRA-COMPUTE/compute_stress_mop_profile.h | 1 + 1 file changed, 1 insertion(+) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop_profile.h b/src/EXTRA-COMPUTE/compute_stress_mop_profile.h index 7bb87fa4c6..1015d6e3eb 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop_profile.h +++ b/src/EXTRA-COMPUTE/compute_stress_mop_profile.h @@ -52,6 +52,7 @@ class ComputeStressMopProfile : public Compute { double **coord, **coordp; double **values_local, **values_global; double **bond_local, **bond_global; + double **local_contribution; double dt, nktv2p, ftm2v; double area; From 9dc1f45e1e7df4915740479fb29bd1944c99aecf Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Fri, 2 Jun 2023 20:08:59 +0300 Subject: [PATCH 26/38] Create/destroy local_contribution --- src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp index d15067a266..3a60342887 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp @@ -129,6 +129,7 @@ ComputeStressMopProfile::ComputeStressMopProfile(LAMMPS *lmp, int narg, char **a values_local = values_global = array = nullptr; bond_local = nullptr; bond_global = nullptr; + local_contribution = nullptr; // bin setup @@ -157,6 +158,7 @@ ComputeStressMopProfile::~ComputeStressMopProfile() memory->destroy(values_global); memory->destroy(bond_local); memory->destroy(bond_global); + memory->destroy(local_contribution); memory->destroy(array); } @@ -514,13 +516,15 @@ void ComputeStressMopProfile::compute_bonds() double dx[3] {0}; double x_bond_1[3] {0}; double x_bond_2[3] {0}; - double local_contribution[nbins][3] {0}; // initialization for (int m {0}; m < nbins; m++) { for (int i {0}; i < nvalues; i++) { bond_local[m][i] = 0.0; } + local_contribution[m][0] = 0.0; + local_contribution[m][1] = 0.0; + local_contribution[m][2] = 0.0; } // loop over all bonded atoms in the current proc @@ -647,6 +651,7 @@ void ComputeStressMopProfile::setup_bins() 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(local_contribution,nbins,3,"stress/mop/profile:local_contribution"); // set bin coordinates for (i = 0; i < nbins; i++) { From 390888179fa2282ac9cd60f5eb5745260708d1a3 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Fri, 2 Jun 2023 20:24:32 +0300 Subject: [PATCH 27/38] Update compute_stress_mop_profile.cpp --- src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp index 3a60342887..d299633731 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp @@ -102,7 +102,7 @@ ComputeStressMopProfile::ComputeStressMopProfile(LAMMPS *lmp, int narg, char **a } } else if (strcmp(arg[iarg],"bond") == 0) { for (i=0; i<3; i++) { - which[nvalues] = PAIR; + which[nvalues] = BOND; nvalues++; } } else error->all(FLERR, "Illegal compute stress/mop/profile command"); //break; From d85342cd6da16eded485da0e544bc35eddf35b94 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Sat, 3 Jun 2023 16:04:40 +0300 Subject: [PATCH 28/38] Update test_compute_global.cpp --- unittest/commands/test_compute_global.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/unittest/commands/test_compute_global.cpp b/unittest/commands/test_compute_global.cpp index 66a24f049f..04913f51b5 100644 --- a/unittest/commands/test_compute_global.cpp +++ b/unittest/commands/test_compute_global.cpp @@ -169,7 +169,7 @@ TEST_F(ComputeGlobalTest, Geometry) command("compute mom1 all momentum"); command("compute mom2 allwater momentum"); command("compute mop1 all stress/mop x 0.0 total"); - command("compute mop2 all stress/mop/profile z lower 0.5 kin conf"); + command("compute mop2 all stress/mop/profile z lower 0.5 kin pair"); thermo_style += " c_mu1 c_mu2 c_mop1[*] c_mop2[1][1]"; } From a0057d674ff2ef0c7f9e99fb713852f7c1d3ab49 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Sun, 4 Jun 2023 13:40:31 +0300 Subject: [PATCH 29/38] Update compute_stress_mop_profile.cpp --- src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp index d299633731..99f522096c 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp @@ -207,9 +207,11 @@ void ComputeStressMopProfile::init() if (force->bond!=nullptr) bondflag = 1; if (force->angle!=nullptr) - error->all(FLERR,"compute stress/mop/profile does not account for angle potentials"); + 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!=nullptr) - error->all(FLERR,"compute stress/mop/profile does not account for dihedral potentials"); + 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->improper!=nullptr) error->all(FLERR,"compute stress/mop/profile does not account for improper potentials"); if (force->kspace!=nullptr) From d8fad4db1516f2a9a3d5a2ee1e9fda7b6959089c Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Sun, 4 Jun 2023 13:45:00 +0300 Subject: [PATCH 30/38] remove white space from compute_stress_mop.cpp --- src/EXTRA-COMPUTE/compute_stress_mop.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop.cpp b/src/EXTRA-COMPUTE/compute_stress_mop.cpp index 3ae758b875..c9557a7ef3 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop.cpp @@ -83,7 +83,7 @@ ComputeStressMop::ComputeStressMop(LAMMPS *lmp, int narg, char **arg) : if (pos >domain->boxhi[dir] || pos boxlo[dir]) error->all(FLERR, "Plane for compute stress/mop is out of bounds"); } - + if (pos < (domain->boxlo[dir]+domain->prd_half[dir])) { pos1 = pos + domain->prd[dir]; } else { @@ -285,7 +285,7 @@ void ComputeStressMop::compute_vector() // sum angle contribution over all procs MPI_Allreduce(angle_local,angle_global,nvalues,MPI_DOUBLE,MPI_SUM,world); - + for (int m=0; m Date: Sun, 4 Jun 2023 13:45:52 +0300 Subject: [PATCH 31/38] remove whitespace from compute_stress_mop_profile.cpp --- src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp index 99f522096c..4f91f8f7fb 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp @@ -609,7 +609,7 @@ void ComputeStressMopProfile::compute_bonds() } m += 3; } -} +} /* ---------------------------------------------------------------------- setup 1d bins and their extent and coordinates From 3782eeee2b2bca5435ef51f13082e96a1e0151df Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Sun, 4 Jun 2023 13:47:54 +0300 Subject: [PATCH 32/38] remove whitespace from compute_stress_mop.rst --- doc/src/compute_stress_mop.rst | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/doc/src/compute_stress_mop.rst b/doc/src/compute_stress_mop.rst index e1b151ebbd..7ab8c58a76 100644 --- a/doc/src/compute_stress_mop.rst +++ b/doc/src/compute_stress_mop.rst @@ -68,7 +68,7 @@ Between one and six 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. +The angle keyword is currently available only for the stress/mop command and not the stress/mop/profile. NOTE 1: The configurational stress is computed considering all pairs of atoms where at least one atom belongs to group group-ID. @@ -119,12 +119,12 @@ 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. In -particular, compute *stress/mop/profile* does not work with more than +particular, compute *stress/mop/profile* does 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 +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() +single() implemented and all angle interactions with the class method born_matrix() implemented. Related commands From c25999d2087c6b79b71ab1d0f6934a7c4598087a Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Sun, 4 Jun 2023 16:08:33 +0300 Subject: [PATCH 33/38] Update compute_stress_mop_profile.cpp --- src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp index 4f91f8f7fb..2c2260d22c 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp @@ -213,7 +213,8 @@ void ComputeStressMopProfile::init() 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->improper!=nullptr) - error->all(FLERR,"compute stress/mop/profile does not account for improper potentials"); + if ((strcmp(force->improper_style, "zero") != 0) && (strcmp(force->improper_style, "none") != 0)) + error->all(FLERR,"compute stress/mop/profile does not account for improper potentials"); if (force->kspace!=nullptr) error->all(FLERR,"compute stress/mop/profile does not account for kspace contributions"); } From 0cff31060b5342f69d0b8d382fab423f29677fb4 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Sun, 4 Jun 2023 16:10:24 +0300 Subject: [PATCH 34/38] Update compute_stress_mop.cpp --- src/EXTRA-COMPUTE/compute_stress_mop.cpp | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop.cpp b/src/EXTRA-COMPUTE/compute_stress_mop.cpp index c9557a7ef3..50eb3425f8 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop.cpp @@ -232,9 +232,11 @@ void ComputeStressMop::init() angleflag = 1; } if (force->dihedral!=nullptr) - error->all(FLERR,"compute stress/mop does not account for dihedral potentials"); + 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->improper!=nullptr) - error->all(FLERR,"compute stress/mop does not account for improper potentials"); + if ((strcmp(force->improper_style, "zero") != 0) && (strcmp(force->improper_style, "none") != 0)) + error->all(FLERR,"compute stress/mop does not account for improper potentials"); if (force->kspace!=nullptr) error->all(FLERR,"compute stress/mop does not account for kspace contributions"); } From ead5a28d354297f20ea40b090aa19f87968a9ce6 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Mon, 5 Jun 2023 18:53:45 +0300 Subject: [PATCH 35/38] Update compute_stress_mop.cpp --- src/EXTRA-COMPUTE/compute_stress_mop.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop.cpp b/src/EXTRA-COMPUTE/compute_stress_mop.cpp index 50eb3425f8..e5d0e2ed13 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop.cpp @@ -227,7 +227,8 @@ void ComputeStressMop::init() if (force->bond!=nullptr) bondflag = 1; if (force->angle!=nullptr) if (force->angle->born_matrix_enable == 0) { - error->all(FLERR,"compute stress/mop does not account for angle potentials"); + if ((strcmp(force->angle_style, "zero") != 0) && (strcmp(force->angle_style, "none") != 0)) + error->all(FLERR,"compute stress/mop does not account for angle potentials"); } else { angleflag = 1; } @@ -238,7 +239,7 @@ void ComputeStressMop::init() if ((strcmp(force->improper_style, "zero") != 0) && (strcmp(force->improper_style, "none") != 0)) error->all(FLERR,"compute stress/mop does not account for improper potentials"); if (force->kspace!=nullptr) - error->all(FLERR,"compute stress/mop does not account for kspace contributions"); + error->warning(FLERR,"compute stress/mop does not account for kspace contributions"); } // need an occasional half neighbor list From 8eed55b56c8330e46fd3715fca962140a0691b52 Mon Sep 17 00:00:00 2001 From: Evangelos Voyiatzis Date: Mon, 5 Jun 2023 18:54:26 +0300 Subject: [PATCH 36/38] Update compute_stress_mop_profile.cpp --- src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp index 2c2260d22c..807092b37e 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop_profile.cpp @@ -216,7 +216,7 @@ void ComputeStressMopProfile::init() if ((strcmp(force->improper_style, "zero") != 0) && (strcmp(force->improper_style, "none") != 0)) error->all(FLERR,"compute stress/mop/profile does not account for improper potentials"); if (force->kspace!=nullptr) - error->all(FLERR,"compute stress/mop/profile does not account for kspace contributions"); + error->warning(FLERR,"compute stress/mop/profile does not account for kspace contributions"); } // need an occasional half neighbor list From 0f925f7a39ea4fccd20e5bd35cf81dc56dd00296 Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 7 Jun 2023 04:33:08 -0400 Subject: [PATCH 37/38] reformat, add versionadded marker --- doc/src/compute_stress_mop.rst | 123 ++++++++++++++++++--------------- 1 file changed, 69 insertions(+), 54 deletions(-) diff --git a/doc/src/compute_stress_mop.rst b/doc/src/compute_stress_mop.rst index 7ab8c58a76..5662fc76d4 100644 --- a/doc/src/compute_stress_mop.rst +++ b/doc/src/compute_stress_mop.rst @@ -45,92 +45,107 @@ Examples Description """"""""""" -Compute *stress/mop* and compute *stress/mop/profile* define computations that -calculate components of the local stress tensor using the method of -planes :ref:`(Todd) `. Specifically in compute *stress/mop* calculates 3 -components are computed in directions *dir*,\ *x*\ ; *dir*,\ *y*\ ; and -*dir*,\ *z*\ ; where *dir* is the direction normal to the plane, while -in compute *stress/mop/profile* the profile of the stress is computed. +Compute *stress/mop* and compute *stress/mop/profile* define +computations that calculate components of the local stress tensor using +the method of planes :ref:`(Todd) `. Specifically in compute +*stress/mop* calculates 3 components are computed in directions *dir*,\ +*x*\ ; *dir*,\ *y*\ ; and *dir*,\ *z*\ ; where *dir* is the direction +normal to the plane, while in compute *stress/mop/profile* the profile +of the stress is computed. Contrary to methods based on histograms of atomic stress (i.e., using -:doc:`compute stress/atom `), the method of planes is -compatible with mechanical balance in heterogeneous systems and at +:doc:`compute stress/atom `), the method of planes +is compatible with mechanical balance in heterogeneous systems and at interfaces :ref:`(Todd) `. The stress tensor is the sum of a kinetic term and a configurational term, which are given respectively by Eq. (21) and Eq. (16) in -:ref:`(Todd) `. For the kinetic part, the algorithm considers that -atoms have crossed the plane if their positions at times :math:`t-\Delta t` -and :math:`t` are one on either side of the plane, and uses the velocity at -time :math:`t-\Delta t/2` given by the velocity Verlet algorithm. +:ref:`(Todd) `. For the kinetic part, the algorithm considers +that atoms have crossed the plane if their positions at times +:math:`t-\Delta t` and :math:`t` are one on either side of the plane, +and uses the velocity at time :math:`t-\Delta t/2` given by the velocity +Verlet algorithm. -Between one and six 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. +.. versionadded:: TBD -NOTE 1: The configurational stress is computed considering all pairs of atoms where at least one atom belongs to group group-ID. + contributions from bond and angle potentials + +Between one and six 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. + +NOTE 1: The configurational stress is computed considering all pairs of +atoms where at least one atom belongs to group group-ID. NOTE 2: The local stress does not include any Lennard-Jones tail -corrections to the stress added by the :doc:`pair_modify tail yes ` -command, since those are contributions to the global system pressure. +corrections to the stress added by the :doc:`pair_modify tail yes +` command, since those are contributions to the global +system pressure. -NOTE 3: The local stress profile generated by compute *stress/mop/profile* -is similar to that obtained by compute -:doc:`stress/cartesian `. -A key difference is that compute *stress/mop/profile* considers particles -crossing a set of planes, while compute *stress/cartesian* computes averages -for a set of small volumes. More information -on the similarities and differences can be found in -:ref:`(Ikeshoji)`. +NOTE 3: The local stress profile generated by compute +*stress/mop/profile* is similar to that obtained by compute +:doc:`stress/cartesian `. A key difference is +that compute *stress/mop/profile* considers particles crossing a set of +planes, while compute *stress/cartesian* computes averages for a set of +small volumes. More information on the similarities and differences can +be found in :ref:`(Ikeshoji)`. Output info """"""""""" -Compute *stress/mop* calculates a global vector (indices starting at 1), with 3 -values for each declared keyword (in the order the keywords have been -declared). For each keyword, the stress tensor components are ordered as -follows: stress_dir,x, stress_dir,y, and stress_dir,z. +Compute *stress/mop* calculates a global vector (indices starting at 1), +with 3 values for each declared keyword (in the order the keywords have +been declared). For each keyword, the stress tensor components are +ordered as follows: stress_dir,x, stress_dir,y, and stress_dir,z. -Compute *stress/mop/profile* instead calculates a global array, with 1 column -giving the position of the planes where the stress tensor was computed, -and with 3 columns of values for each declared keyword (in the order the -keywords have been declared). For each keyword, the profiles of stress -tensor components are ordered as follows: stress_dir,x; stress_dir,y; -and stress_dir,z. +Compute *stress/mop/profile* instead calculates a global array, with 1 +column giving the position of the planes where the stress tensor was +computed, and with 3 columns of values for each declared keyword (in the +order the keywords have been declared). For each keyword, the profiles +of stress tensor components are ordered as follows: stress_dir,x; +stress_dir,y; and stress_dir,z. The values are in pressure :doc:`units `. -The values produced by this compute can be accessed by various :doc:`output commands `. -For instance, the results can be written to a file using the -:doc:`fix ave/time ` command. Please see the example -in the examples/PACKAGES/mop folder. +The values produced by this compute can be accessed by various +:doc:`output commands `. For instance, the results can be +written to a file using the :doc:`fix ave/time ` +command. Please see the example in the examples/PACKAGES/mop folder. Restrictions """""""""""" -These styles are part of the EXTRA-COMPUTE package. They are only enabled if -LAMMPS is built with that package. See the :doc:`Build package ` -doc page on for more info. +These styles are part of the EXTRA-COMPUTE package. They are only +enabled if LAMMPS is built with that package. See the :doc:`Build +package ` doc page on for more info. The method is only implemented for 3d orthogonal simulation boxes whose 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. In -particular, compute *stress/mop/profile* does 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. +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 +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. Related commands """""""""""""""" -:doc:`compute stress/atom `, :doc:`compute pressure `, :doc:`compute stress/cartesian `, :doc:`compute stress/cylinder `, :doc:`compute stress/spherical ` +:doc:`compute stress/atom `, +:doc:`compute pressure `, +:doc:`compute stress/cartesian `, +:doc:`compute stress/cylinder `, +:doc:`compute stress/spherical ` Default """"""" From 53b1af7720c9246eec9cbadf6ef0cb39bcc30fad Mon Sep 17 00:00:00 2001 From: Axel Kohlmeyer Date: Wed, 7 Jun 2023 05:11:53 -0400 Subject: [PATCH 38/38] LAMMPS programming style/conventions updates --- src/EXTRA-COMPUTE/compute_stress_mop.cpp | 171 ++++++++---------- src/EXTRA-COMPUTE/compute_stress_mop.h | 2 +- .../compute_stress_mop_profile.cpp | 96 +++++----- .../compute_stress_mop_profile.h | 2 +- 4 files changed, 126 insertions(+), 145 deletions(-) diff --git a/src/EXTRA-COMPUTE/compute_stress_mop.cpp b/src/EXTRA-COMPUTE/compute_stress_mop.cpp index e5d0e2ed13..98e3bf7043 100644 --- a/src/EXTRA-COMPUTE/compute_stress_mop.cpp +++ b/src/EXTRA-COMPUTE/compute_stress_mop.cpp @@ -1,4 +1,3 @@ -// clang-format off /* ---------------------------------------------------------------------- LAMMPS - Large-scale Atomic/Molecular Massively Parallel Simulator https://www.lammps.org/, Sandia National Laboratories @@ -14,6 +13,7 @@ /*------------------------------------------------------------------------ Contributing Authors : Romain Vermorel (LFCR), Laurent Joly (ULyon) + Support for bonds and angles added by : Evangelos Voyiatzis (NovaMechanics) --------------------------------------------------------------------------*/ #include "compute_stress_mop.h" @@ -22,6 +22,7 @@ #include "atom.h" #include "atom_vec.h" #include "bond.h" +#include "comm.h" #include "domain.h" #include "error.h" #include "force.h" @@ -37,50 +38,49 @@ using namespace LAMMPS_NS; -enum{X,Y,Z}; -enum{TOTAL,CONF,KIN,PAIR,BOND,ANGLE}; +enum { X, Y, Z }; +enum { TOTAL, CONF, KIN, PAIR, BOND, ANGLE }; +// clang-format off /* ---------------------------------------------------------------------- */ ComputeStressMop::ComputeStressMop(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { - if (narg < 6) error->all(FLERR,"Illegal compute stress/mop command"); - - MPI_Comm_rank(world,&me); + if (narg < 6) utils::missing_cmd_args(FLERR, "compute stress/mop", error); bondflag = 0; angleflag = 0; // set compute mode and direction of plane(s) for pressure calculation - if (strcmp(arg[3],"x")==0) { + if (strcmp(arg[3],"x") == 0) { dir = X; - } else if (strcmp(arg[3],"y")==0) { + } else if (strcmp(arg[3],"y") == 0) { dir = Y; - } else if (strcmp(arg[3],"z")==0) { + } else if (strcmp(arg[3],"z") == 0) { dir = Z; } else error->all(FLERR,"Illegal compute stress/mop command"); // Position of the plane - if (strcmp(arg[4],"lower")==0) { + if (strcmp(arg[4],"lower") == 0) { pos = domain->boxlo[dir]; - } else if (strcmp(arg[4],"upper")==0) { + } else if (strcmp(arg[4],"upper") == 0) { pos = domain->boxhi[dir]; - } else if (strcmp(arg[4],"center")==0) { + } else if (strcmp(arg[4],"center") == 0) { pos = 0.5*(domain->boxlo[dir]+domain->boxhi[dir]); } else pos = utils::numeric(FLERR,arg[4],false,lmp); // plane inside the box - if (pos >domain->boxhi[dir] || pos boxlo[dir]) { - error->warning(FLERR,"The specified initial plane lies outside of the simulation box"); - double dx[3] {0}; + if ((pos > domain->boxhi[dir]) || (pos < domain->boxlo[dir])) { + error->warning(FLERR, "The specified initial plane lies outside of the simulation box"); + double dx[3] = {0.0, 0.0, 0.0}; dx[dir] = pos - 0.5*(domain->boxhi[dir] + domain->boxlo[dir]); domain->minimum_image(dx[0], dx[1], dx[2]); pos = 0.5*(domain->boxhi[dir] + domain->boxlo[dir]) + dx[dir]; - if (pos >domain->boxhi[dir] || pos boxlo[dir]) + if ((pos > domain->boxhi[dir]) || (pos < domain->boxlo[dir])) error->all(FLERR, "Plane for compute stress/mop is out of bounds"); } @@ -133,15 +133,15 @@ ComputeStressMop::ComputeStressMop(LAMMPS *lmp, int narg, char **arg) : iarg++; } - // Error check + // Error checks // 3D only - if (domain->dimension < 3) - error->all(FLERR, "Compute stress/mop incompatible with simulation dimension"); + if (domain->dimension != 3) + error->all(FLERR, "Compute stress/mop requires a 3d system"); // orthogonal simulation box if (domain->triclinic != 0) - error->all(FLERR, "Compute stress/mop incompatible with triclinic simulation box"); + error->all(FLERR, "Compute stress/mop is incompatible with triclinic simulation box"); // Initialize some variables @@ -154,12 +154,12 @@ ComputeStressMop::ComputeStressMop(LAMMPS *lmp, int narg, char **arg) : // this fix produces a global vector memory->create(vector,nvalues,"stress/mop:vector"); - memory->create(values_local,nvalues,"stress/mop/spatial:values_local"); - memory->create(values_global,nvalues,"stress/mop/spatial:values_global"); - memory->create(bond_local,nvalues,"stress/mop/spatial:bond_local"); - memory->create(bond_global,nvalues,"stress/mop/spatial:bond_global"); - memory->create(angle_local,nvalues,"stress/mop/spatial:angle_local"); - memory->create(angle_global,nvalues,"stress/mop/spatial:angle_global"); + memory->create(values_local,nvalues,"stress/mop:values_local"); + memory->create(values_global,nvalues,"stress/mop:values_global"); + memory->create(bond_local,nvalues,"stress/mop:bond_local"); + 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"); size_vector = nvalues; vector_flag = 1; @@ -171,8 +171,7 @@ ComputeStressMop::ComputeStressMop(LAMMPS *lmp, int narg, char **arg) : ComputeStressMop::~ComputeStressMop() { - - delete [] which; + delete[] which; memory->destroy(values_local); memory->destroy(values_global); @@ -209,36 +208,39 @@ void ComputeStressMop::init() // Compute stress/mop requires fixed simulation box if (domain->box_change_size || domain->box_change_shape || domain->deform_flag) - error->all(FLERR, "Compute stress/mop requires a fixed simulation box"); + error->all(FLERR, "Compute stress/mop requires a fixed size simulation box"); // This compute requires a pair style with pair_single method implemented - if (force->pair == nullptr) + if (!force->pair) error->all(FLERR,"No pair style is defined for compute stress/mop"); if (force->pair->single_enable == 0) error->all(FLERR,"Pair style does not support compute stress/mop"); // Errors - if (me==0) { + if (comm->me == 0) { // issue an error for unimplemented intramolecular potentials or Kspace. - if (force->bond!=nullptr) bondflag = 1; - if (force->angle!=nullptr) + if (force->bond) bondflag = 1; + 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 does not account for angle potentials"); } else { angleflag = 1; } - if (force->dihedral!=nullptr) + } + 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->improper!=nullptr) + } + if (force->improper) { if ((strcmp(force->improper_style, "zero") != 0) && (strcmp(force->improper_style, "none") != 0)) error->all(FLERR,"compute stress/mop does not account for improper potentials"); - if (force->kspace!=nullptr) + } + if (force->kspace) error->warning(FLERR,"compute stress/mop does not account for kspace contributions"); } @@ -266,8 +268,7 @@ void ComputeStressMop::compute_vector() compute_pairs(); // sum pressure contributions over all procs - MPI_Allreduce(values_local,values_global,nvalues, - MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(values_local,values_global,nvalues,MPI_DOUBLE,MPI_SUM,world); if (bondflag) { //Compute bond contribution on separate procs @@ -292,7 +293,6 @@ void ComputeStressMop::compute_vector() for (int m=0; mpos) && (xj[dir]pos1) && (xj[dir] pos) && (xj[dir] < pos)) || ((xi[dir] > pos1) && (xj[dir] < pos1))) { pair->single(i,j,itype,jtype,rsq,factor_coul,factor_lj,fpair); - values_local[m] += fpair*(xi[0]-xj[0])/area*nktv2p; values_local[m+1] += fpair*(xi[1]-xj[1])/area*nktv2p; values_local[m+2] += fpair*(xi[2]-xj[2])/area*nktv2p; - } - else if (((xi[dir]pos)) || ((xi[dir]pos1))) { - + } else if (((xi[dir] < pos) && (xj[dir] > pos)) || ((xi[dir] < pos1) && (xj[dir] > pos1))) { pair->single(i,j,itype,jtype,rsq,factor_coul,factor_lj,fpair); - values_local[m] -= fpair*(xi[0]-xj[0])/area*nktv2p; values_local[m+1] -= fpair*(xi[1]-xj[1])/area*nktv2p; values_local[m+2] -= fpair*(xi[2]-xj[2])/area*nktv2p; } - } else { - - if (((xi[dir]>pos) && (xj[dir]pos1) && (xj[dir] pos) && (xj[dir] < pos)) || ((xi[dir] > pos1) && (xj[dir] < pos1))) { pair->single(i,j,itype,jtype,rsq,factor_coul,factor_lj,fpair); - values_local[m] += fpair*(xi[0]-xj[0])/area*nktv2p; values_local[m+1] += fpair*(xi[1]-xj[1])/area*nktv2p; values_local[m+2] += fpair*(xi[2]-xj[2])/area*nktv2p; } - } - } - } - } // Compute kinetic contribution to pressure // counts local particles transfers across the plane - if (which[m] == KIN || which[m] == TOTAL) { + if ((which[m] == KIN) || (which[m] == TOTAL)) { double sgn; for (int i = 0; i < nlocal; i++) { @@ -447,13 +434,13 @@ void ComputeStressMop::compute_pairs() //coordinates at t-dt (based on Velocity-Verlet alg.) if (rmass) { - xj[0] = xi[0]-vi[0]*dt+fi[0]/2/rmass[i]*dt*dt*ftm2v; - xj[1] = xi[1]-vi[1]*dt+fi[1]/2/rmass[i]*dt*dt*ftm2v; - xj[2] = xi[2]-vi[2]*dt+fi[2]/2/rmass[i]*dt*dt*ftm2v; + xj[0] = xi[0]-vi[0]*dt+fi[0]/2.0/rmass[i]*dt*dt*ftm2v; + xj[1] = xi[1]-vi[1]*dt+fi[1]/2.0/rmass[i]*dt*dt*ftm2v; + xj[2] = xi[2]-vi[2]*dt+fi[2]/2.0/rmass[i]*dt*dt*ftm2v; } else { - xj[0] = xi[0]-vi[0]*dt+fi[0]/2/mass[itype]*dt*dt*ftm2v; - xj[1] = xi[1]-vi[1]*dt+fi[1]/2/mass[itype]*dt*dt*ftm2v; - xj[2] = xi[2]-vi[2]*dt+fi[2]/2/mass[itype]*dt*dt*ftm2v; + xj[0] = xi[0]-vi[0]*dt+fi[0]/2.0/mass[itype]*dt*dt*ftm2v; + xj[1] = xi[1]-vi[1]*dt+fi[1]/2.0/mass[itype]*dt*dt*ftm2v; + xj[2] = xi[2]-vi[2]*dt+fi[2]/2.0/mass[itype]*dt*dt*ftm2v; } // because LAMMPS does not put atoms back in the box @@ -463,21 +450,21 @@ void ComputeStressMop::compute_pairs() double pos_temp = pos+copysign(1.0,domain->prd_half[dir]-pos)*domain->prd[dir]; if (fabs(xi[dir]-pos)bond; - double dx[3] {0}; - double x_bond_1[3] {0}; - double x_bond_2[3] {0}; - double local_contribution[3] {0}; + double dx[3] = {0.0, 0.0, 0.0}; + double x_bond_1[3] = {0.0, 0.0, 0.0}; + double x_bond_2[3] = {0.0, 0.0, 0.0}; + double local_contribution[3] = {0.0, 0.0, 0.0}; // initialization - for (int i {0}; i < nvalues; i++) bond_local[i] = 0.0; + for (int i = 0; i < nvalues; i++) bond_local[i] = 0.0; // loop over all bonded atoms in the current proc for (atom1 = 0; atom1 < nlocal; atom1++) { @@ -591,7 +578,7 @@ void ComputeStressMop::compute_bonds() } // loop over the keywords and if necessary add the bond contribution - int m {0}; + int m = 0; while (mangle; double duang, du2ang; - double dx[3] {0}; - double dx_left[3] {0}; - double dx_right[3] {0}; - double x_angle_left[3] {0}; - double x_angle_middle[3] {0}; - double x_angle_right[3] {0}; - double dcos_theta[3] {0}; - double local_contribution[3] {0}; + 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}; + double local_contribution[3] = {0.0, 0.0, 0.0}; // initialization - for (int i {0}; i < nvalues; i++) angle_local[i] = 0.0; + for (int i = 0; i < nvalues; i++) angle_local[i] = 0.0; for (atom2 = 0; atom2 < nlocal; atom2++) { if (!(mask[atom2] & groupbit)) continue; @@ -762,8 +749,8 @@ void ComputeStressMop::compute_angles() } } // loop over the keywords and if necessary add the angle contribution - int m {0}; - while (m #include using namespace LAMMPS_NS; -enum{X,Y,Z}; -enum{LOWER,CENTER,UPPER,COORD}; -enum{TOTAL,CONF,KIN,PAIR,BOND}; +enum { X, Y, Z }; +enum { LOWER, CENTER, UPPER, COORD }; +enum { TOTAL, CONF, KIN, PAIR, BOND }; +// clang-format off /* ---------------------------------------------------------------------- */ ComputeStressMopProfile::ComputeStressMopProfile(LAMMPS *lmp, int narg, char **arg) : Compute(lmp, narg, arg) { - if (narg < 7) error->all(FLERR,"Illegal compute stress/mop/profile command"); - - MPI_Comm_rank(world,&me); + if (narg < 7) utils::missing_cmd_args(FLERR, "compute stress/mop/profile", error); bondflag = 0; // set compute mode and direction of plane(s) for pressure calculation - if (strcmp(arg[3],"x")==0) { + if (strcmp(arg[3],"x") == 0) { dir = X; - } else if (strcmp(arg[3],"y")==0) { + } else if (strcmp(arg[3],"y") == 0) { dir = Y; - } else if (strcmp(arg[3],"z")==0) { + } else if (strcmp(arg[3],"z") == 0) { dir = Z; } else error->all(FLERR,"Illegal compute stress/mop/profile command"); @@ -67,8 +67,7 @@ ComputeStressMopProfile::ComputeStressMopProfile(LAMMPS *lmp, int narg, char **a 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 (originflag == COORD) origin = utils::numeric(FLERR,arg[4],false,lmp); delta = utils::numeric(FLERR,arg[5],false,lmp); invdelta = 1.0/delta; @@ -149,8 +148,7 @@ ComputeStressMopProfile::ComputeStressMopProfile(LAMMPS *lmp, int narg, char **a ComputeStressMopProfile::~ComputeStressMopProfile() { - - delete [] which; + delete[] which; memory->destroy(coord); memory->destroy(coordp); @@ -166,7 +164,6 @@ ComputeStressMopProfile::~ComputeStressMopProfile() void ComputeStressMopProfile::init() { - // conversion constants nktv2p = force->nktv2p; @@ -192,30 +189,30 @@ void ComputeStressMopProfile::init() //This compute requires a pair style with pair_single method implemented - if (force->pair == nullptr) + if (!force->pair) error->all(FLERR,"No pair style is defined for compute stress/mop/profile"); if (force->pair->single_enable == 0) error->all(FLERR,"Pair style does not support compute stress/mop/profile"); // Errors - if (me==0) { + if (comm->me == 0) { //Compute stress/mop/profile only accounts for pair interactions. // issue an error if any intramolecular potential or Kspace is defined. - if (force->bond!=nullptr) bondflag = 1; + if (force->bond) bondflag = 1; - if (force->angle!=nullptr) + 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!=nullptr) + 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->improper!=nullptr) + if (force->improper) if ((strcmp(force->improper_style, "zero") != 0) && (strcmp(force->improper_style, "none") != 0)) error->all(FLERR,"compute stress/mop/profile does not account for improper potentials"); - if (force->kspace!=nullptr) + if (force->kspace) error->warning(FLERR,"compute stress/mop/profile does not account for kspace contributions"); } @@ -244,8 +241,7 @@ void ComputeStressMopProfile::compute_array() compute_pairs(); // sum pressure contributions over all procs - MPI_Allreduce(&values_local[0][0],&values_global[0][0],nbins*nvalues, - MPI_DOUBLE,MPI_SUM,world); + MPI_Allreduce(&values_local[0][0],&values_global[0][0],nbins*nvalues,MPI_DOUBLE,MPI_SUM,world); if (bondflag) { //Compute bond contribution on separate procs @@ -259,16 +255,14 @@ 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); + MPI_Allreduce(&bond_local[0][0],&bond_global[0][0],nbins*nvalues,MPI_DOUBLE,MPI_SUM,world); - int ibin,m,mo; - for (ibin=0; ibinbond; - double dx[3] {0}; - double x_bond_1[3] {0}; - double x_bond_2[3] {0}; + double dx[3] = {0.0, 0.0, 0.0}; + double x_bond_1[3] = {0.0, 0.0, 0.0}; + double x_bond_2[3] = {0.0, 0.0, 0.0}; // initialization - for (int m {0}; m < nbins; m++) { - for (int i {0}; i < nvalues; i++) { + for (int m = 0; m < nbins; m++) { + for (int i = 0; i < nvalues; i++) { bond_local[m][i] = 0.0; } local_contribution[m][0] = 0.0; @@ -557,7 +551,7 @@ void ComputeStressMopProfile::compute_bonds() if (newton_bond == 0 && tag[atom1] > tag[atom2]) continue; if (btype <= 0) continue; - for (int ibin {0}; ibin