ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Add ddx, dy and dz computes to compute bond/local and property/local

Browse Source
This commit is contained in:
Tim Bernhard
2021-11-18 17:22:32 +01:00
parent 515ef7bece
commit 36e4e3e746
4 changed files with 41 additions and 9 deletions
Download Patch File Download Diff File Expand all files Collapse all files

6
doc/src/compute_bond_local.rst
View File

@ -13,7 +13,7 @@ Syntax
* ID, group-ID are documented in :doc:`compute <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.

14
doc/src/compute_pair_local.rst
View File

@ -13,11 +13,12 @@ Syntax
* ID, group-ID are documented in :doc:`compute <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 <units>`. The value
*eng* will be in energy :doc:`units <units>`. The values *force*, *fx*,
*fy*, and *fz* will be in force :doc:`units <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 <units>`.
The value *eng* will be in energy :doc:`units <units>`.
The values *force*, *fx*, *fy*, and *fz* will be in force :doc:`units <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*,

17
src/compute_bond_local.cpp
View File

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

13
src/compute_pair_local.cpp
View File

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