From 36e4e3e7462512cc04fbcee7dcfbb345080863af Mon Sep 17 00:00:00 2001 From: Tim Bernhard Date: Thu, 18 Nov 2021 17:22:32 +0100 Subject: [PATCH] Add ddx, dy and dz computes to compute bond/local and property/local --- doc/src/compute_bond_local.rst | 6 +++++- doc/src/compute_pair_local.rst | 14 +++++++++----- src/compute_bond_local.cpp | 17 ++++++++++++++++- src/compute_pair_local.cpp | 13 +++++++++++-- 4 files changed, 41 insertions(+), 9 deletions(-) diff --git a/doc/src/compute_bond_local.rst b/doc/src/compute_bond_local.rst index 8bdde70dd9..cf98eb4ef3 100644 --- a/doc/src/compute_bond_local.rst +++ b/doc/src/compute_bond_local.rst @@ -13,7 +13,7 @@ Syntax * ID, group-ID are documented in :doc:`compute ` command * bond/local = style name of this compute command * one or more values may be appended -* value = *dist* or *engpot* or *force* or *fx* or *fy* or *fz* or *engvib* or *engrot* or *engtrans* or *omega* or *velvib* or *v_name* +* value = *dist* or *dx* or *dy* or *dz* or *engpot* or *force* or *fx* or *fy* or *fz* or *engvib* or *engrot* or *engtrans* or *omega* or *velvib* or *v_name* .. parsed-literal:: @@ -21,6 +21,7 @@ Syntax *engpot* = bond potential energy *force* = bond force + *dx*,\ *dy*,\ *dz* = components of pairwise distance *fx*,\ *fy*,\ *fz* = components of bond force *engvib* = bond kinetic energy of vibration *engrot* = bond kinetic energy of rotation @@ -63,6 +64,9 @@ whether the 2 atoms represent a simple diatomic molecule, or are part of some larger molecule. The value *dist* is the current length of the bond. +The values *dx*, *dy*, and *dz* are the xyz components of the +*distance* between the pair of atoms. This value is always the +distance from the atom of lower to the one of higher id. The value *engpot* is the potential energy for the bond, based on the current separation of the pair of atoms in the bond. diff --git a/doc/src/compute_pair_local.rst b/doc/src/compute_pair_local.rst index f464c7cec6..fec35b1c41 100644 --- a/doc/src/compute_pair_local.rst +++ b/doc/src/compute_pair_local.rst @@ -13,11 +13,12 @@ Syntax * ID, group-ID are documented in :doc:`compute ` command * pair/local = style name of this compute command * one or more values may be appended -* value = *dist* or *eng* or *force* or *fx* or *fy* or *fz* or *pN* +* value = *dist* or *dx* or *dy* or *dz* or *eng* or *force* or *fx* or *fy* or *fz* or *pN* .. parsed-literal:: *dist* = pairwise distance + *dx*,\ *dy*,\ *dz* = components of pairwise distance *eng* = pairwise energy *force* = pairwise force *fx*,\ *fy*,\ *fz* = components of pairwise force @@ -56,6 +57,9 @@ force cutoff distance for that interaction, as defined by the commands. The value *dist* is the distance between the pair of atoms. +The values *dx*, *dy*, and *dz* are the xyz components of the +*distance* between the pair of atoms. This value can be inconsistent +due to changing interactions and neighbor list orders. The value *eng* is the interaction energy for the pair of atoms. @@ -89,10 +93,10 @@ from the second of the two sub-styles. If the referenced *pN* is not computed for the specific pairwise interaction (based on atom types), then the output will be 0.0. -The value *dist* will be in distance :doc:`units `. The value -*eng* will be in energy :doc:`units `. The values *force*, *fx*, -*fy*, and *fz* will be in force :doc:`units `. The values *pN* -will be in whatever units the pair style defines. +The value *dist*, *dx*, *dy* and *dz* will be in distance :doc:`units `. +The value *eng* will be in energy :doc:`units `. +The values *force*, *fx*, *fy*, and *fz* will be in force :doc:`units `. +The values *pN* will be in whatever units the pair style defines. The optional *cutoff* keyword determines how the force cutoff distance for an interaction is determined. For the default setting of *type*, diff --git a/src/compute_bond_local.cpp b/src/compute_bond_local.cpp index ab0fa3fb0a..d721d28247 100644 --- a/src/compute_bond_local.cpp +++ b/src/compute_bond_local.cpp @@ -34,7 +34,7 @@ using namespace LAMMPS_NS; #define DELTA 10000 #define EPSILON 1.0e-12 -enum{DIST,VELVIB,OMEGA,ENGTRANS,ENGVIB,ENGROT,ENGPOT,FORCE,FX,FY,FZ,VARIABLE}; +enum{DIST,DX,DY,DZ,VELVIB,OMEGA,ENGTRANS,ENGVIB,ENGROT,ENGPOT,FORCE,FX,FY,FZ,VARIABLE}; /* ---------------------------------------------------------------------- */ @@ -63,6 +63,9 @@ ComputeBondLocal::ComputeBondLocal(LAMMPS *lmp, int narg, char **arg) : int iarg; for (iarg = 3; iarg < narg; iarg++) { if (strcmp(arg[iarg],"dist") == 0) bstyle[nvalues++] = DIST; + else if (strcmp(arg[iarg],"dx") == 0) bstyle[nvalues++] = DX; + else if (strcmp(arg[iarg],"dy") == 0) bstyle[nvalues++] = DY; + else if (strcmp(arg[iarg],"dz") == 0) bstyle[nvalues++] = DZ; else if (strcmp(arg[iarg],"engpot") == 0) bstyle[nvalues++] = ENGPOT; else if (strcmp(arg[iarg],"force") == 0) bstyle[nvalues++] = FORCE; else if (strcmp(arg[iarg],"fx") == 0) bstyle[nvalues++] = FX; @@ -384,11 +387,23 @@ int ComputeBondLocal::compute_bonds(int flag) if (dstr) input->variable->internal_set(dvar,sqrt(rsq)); } + // to make sure dx, dy and dz are always from the lower to the higher id + double directionCorrection = atom1 > atom2 ? -1.0 : 1.0; + for (int n = 0; n < nvalues; n++) { switch (bstyle[n]) { case DIST: ptr[n] = sqrt(rsq); break; + case DX: + ptr[n] = dx*directionCorrection; + break; + case DY: + ptr[n] = dy*directionCorrection; + break; + case DZ: + ptr[n] = dz*directionCorrection; + break; case ENGPOT: ptr[n] = engpot; break; diff --git a/src/compute_pair_local.cpp b/src/compute_pair_local.cpp index af0f20014c..f783b8435b 100644 --- a/src/compute_pair_local.cpp +++ b/src/compute_pair_local.cpp @@ -31,7 +31,7 @@ using namespace LAMMPS_NS; #define DELTA 10000 -enum{DIST,ENG,FORCE,FX,FY,FZ,PN}; +enum{DIST,ENG,FORCE,FX,FY,FZ,PN,DX,DY,DZ}; enum{TYPE,RADIUS}; /* ---------------------------------------------------------------------- */ @@ -56,6 +56,9 @@ ComputePairLocal::ComputePairLocal(LAMMPS *lmp, int narg, char **arg) : else if (strcmp(arg[iarg],"fx") == 0) pstyle[nvalues++] = FX; else if (strcmp(arg[iarg],"fy") == 0) pstyle[nvalues++] = FY; else if (strcmp(arg[iarg],"fz") == 0) pstyle[nvalues++] = FZ; + else if (strcmp(arg[iarg],"dx") == 0) pstyle[nvalues++] = DX; + else if (strcmp(arg[iarg],"dy") == 0) pstyle[nvalues++] = DY; + else if (strcmp(arg[iarg],"dz") == 0) pstyle[nvalues++] = DZ; else if (arg[iarg][0] == 'p') { int n = atoi(&arg[iarg][1]); if (n <= 0) error->all(FLERR, @@ -92,7 +95,7 @@ ComputePairLocal::ComputePairLocal(LAMMPS *lmp, int narg, char **arg) : singleflag = 0; for (int i = 0; i < nvalues; i++) - if (pstyle[i] != DIST) singleflag = 1; + if (pstyle[i] != DIST && pstyle[i] != DX && pstyle[i] != DY && pstyle[i] != DZ) singleflag = 1; if (nvalues == 1) size_local_cols = 0; else size_local_cols = nvalues; @@ -269,6 +272,12 @@ int ComputePairLocal::compute_pairs(int flag) case DIST: ptr[n] = sqrt(rsq); break; + case DX: + ptr[n] = delx; + case DY: + ptr[n] = dely; + case DZ: + ptr[n] = delz; case ENG: ptr[n] = eng; break;