diff --git a/src/group.cpp b/src/group.cpp index 2c482d31b7..9c369101ad 100644 --- a/src/group.cpp +++ b/src/group.cpp @@ -646,6 +646,31 @@ void Group::vcm(int igroup, double masstotal, double *cm) cm[2] /= masstotal; } +/* ---------------------------------------------------------------------- + compute the total force on group of atoms +------------------------------------------------------------------------- */ + +void Group::fcm(int igroup, double *cm) +{ + int groupbit = bitmask[igroup]; + + double **f = atom->f; + int *mask = atom->mask; + int nlocal = atom->nlocal; + + double flocal[3]; + flocal[0] = flocal[1] = flocal[2] = 0.0; + + for (int i = 0; i < nlocal; i++) + if (mask[i] & groupbit) { + flocal[0] += f[i][0]; + flocal[1] += f[i][1]; + flocal[2] += f[i][2]; + } + + MPI_Allreduce(flocal,cm,3,MPI_DOUBLE,MPI_SUM,world); +} + /* ---------------------------------------------------------------------- compute the radius-of-gyration of group around center-of-mass cm must unwrap atoms to compute Rg correctly diff --git a/src/group.h b/src/group.h index f44e13883a..7d22d0c38c 100644 --- a/src/group.h +++ b/src/group.h @@ -41,6 +41,7 @@ class Group : protected Pointers { void bounds(int, double *); // bounds of atoms in group void xcm(int, double, double *); // center-of-mass coords of group void vcm(int, double, double *); // center-of-mass velocity of group + void fcm(int, double *); // total force on group double gyration(int, double, double *); // radius-of-gyration of group void angmom(int, double *, double *); // angular momentum of group void inertia(int, double *, double [3][3]); // inertia tensor diff --git a/src/style.h b/src/style.h index 8c1e094c08..971d0c7700 100644 --- a/src/style.h +++ b/src/style.h @@ -145,6 +145,7 @@ DumpStyle(xyz,DumpXYZ) #include "fix_nph.h" #include "fix_npt.h" #include "fix_nve.h" +#include "fix_nve_noforce.h" #include "fix_nvt.h" #include "fix_plane_force.h" #include "fix_print.h" @@ -192,6 +193,7 @@ FixStyle(msd,FixMSD) FixStyle(nph,FixNPH) FixStyle(npt,FixNPT) FixStyle(nve,FixNVE) +FixStyle(nve,FixNVENoforce) FixStyle(nvt,FixNVT) FixStyle(orient/fcc,FixOrientFCC) FixStyle(print,FixPrint) diff --git a/src/variable.cpp b/src/variable.cpp index dd687eb6d4..cca6d31046 100644 --- a/src/variable.cpp +++ b/src/variable.cpp @@ -569,7 +569,7 @@ double Variable::evaluate(char *str, Tree *tree) // customize by adding group function to this list and to if statement // mass(group),charge(group),xcm(group,dim),vcm(group,dim), - // bound(group,xmin),gyration(group) + // fcm(group,dim),bound(group,xmin),gyration(group) } else if (strcmp(func,"mass") == 0) { if (arg2) error->all("Cannot evaluate variable"); @@ -610,6 +610,17 @@ double Variable::evaluate(char *str, Tree *tree) else if (strcmp(arg2,"z") == 0) answer = vcm[2]; else error->all("Cannot evaluate variable"); + } else if (strcmp(func,"fcm") == 0) { + if (!arg2) error->all("Cannot evaluate variable"); + int igroup = group->find(arg1); + if (igroup == -1) error->all("Variable group ID does not exist"); + double fcm[3]; + group->fcm(igroup,fcm); + if (strcmp(arg2,"x") == 0) answer = fcm[0]; + else if (strcmp(arg2,"y") == 0) answer = fcm[1]; + else if (strcmp(arg2,"z") == 0) answer = fcm[2]; + else error->all("Cannot evaluate variable"); + } else if (strcmp(func,"bound") == 0) { if (!arg2) error->all("Cannot evaluate variable"); int igroup = group->find(arg1);