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

Fixing bpm/sphere error in fix move, displace atoms

Browse Source
This commit is contained in:
jtclemm
2023-12-15 13:33:45 -07:00
parent 7a1466671b
commit 95d1a41ee4
3 changed files with 49 additions and 14 deletions
Download Patch File Download Diff File Expand all files Collapse all files

16
src/displace_atoms.cpp
View File

@ -160,7 +160,7 @@ void DisplaceAtoms::command(int narg, char **arg)
int *mask = atom->mask; int *mask = atom->mask;
int nlocal = atom->nlocal; int nlocal = atom->nlocal;
double fraction,dramp; double fraction, dramp;
for (i = 0; i < nlocal; i++) { for (i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) { if (mask[i] & groupbit) {
@ -255,11 +255,12 @@ void DisplaceAtoms::command(int narg, char **arg)
int line_flag = atom->line_flag; int line_flag = atom->line_flag;
int tri_flag = atom->tri_flag; int tri_flag = atom->tri_flag;
int body_flag = atom->body_flag; int body_flag = atom->body_flag;
int quat_atom_flag = atom->quat_flag;
int theta_flag = 0; int theta_flag = 0;
int quat_flag = 0; int quat_flag = 0;
if (line_flag) theta_flag = 1; if (line_flag) theta_flag = 1;
if (ellipsoid_flag || tri_flag || body_flag) quat_flag = 1; if (ellipsoid_flag || tri_flag || body_flag || quat_atom_flag) quat_flag = 1;
// AtomVec pointers to retrieve per-atom storage of extra quantities // AtomVec pointers to retrieve per-atom storage of extra quantities
@ -269,6 +270,7 @@ void DisplaceAtoms::command(int narg, char **arg)
auto avec_body = dynamic_cast<AtomVecBody *>(atom->style_match("body")); auto avec_body = dynamic_cast<AtomVecBody *>(atom->style_match("body"));
double **x = atom->x; double **x = atom->x;
double **quat_atom = atom->quat;
int *ellipsoid = atom->ellipsoid; int *ellipsoid = atom->ellipsoid;
int *line = atom->line; int *line = atom->line;
int *tri = atom->tri; int *tri = atom->tri;
@ -313,7 +315,7 @@ void DisplaceAtoms::command(int narg, char **arg)
// quats for ellipsoids, tris, and bodies // quats for ellipsoids, tris, and bodies
if (quat_flag) { if (quat_flag && !quat_atom_flag) {
quat = nullptr; quat = nullptr;
if (ellipsoid_flag && ellipsoid[i] >= 0) if (ellipsoid_flag && ellipsoid[i] >= 0)
quat = avec_ellipsoid->bonus[ellipsoid[i]].quat; quat = avec_ellipsoid->bonus[ellipsoid[i]].quat;
@ -322,12 +324,18 @@ void DisplaceAtoms::command(int narg, char **arg)
else if (body_flag && body[i] >= 0) else if (body_flag && body[i] >= 0)
quat = avec_body->bonus[body[i]].quat; quat = avec_body->bonus[body[i]].quat;
if (quat) { if (quat) {
MathExtra::quatquat(qrotate,quat,qnew); MathExtra::quatquat(qrotate, quat, qnew);
quat[0] = qnew[0]; quat[0] = qnew[0];
quat[1] = qnew[1]; quat[1] = qnew[1];
quat[2] = qnew[2]; quat[2] = qnew[2];
quat[3] = qnew[3]; quat[3] = qnew[3];
} }
} else if (quat_atom_flag) {
MathExtra::quatquat(qrotate, quat_atom[i], qnew);
quat_atom[i][0] = qnew[0];
quat_atom[i][1] = qnew[1];
quat_atom[i][2] = qnew[2];
quat_atom[i][3] = qnew[3];
} }
} }
} }

45
src/fix_move.cpp
View File

@ -276,10 +276,11 @@ FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) :
line_flag = atom->line_flag; line_flag = atom->line_flag;
tri_flag = atom->tri_flag; tri_flag = atom->tri_flag;
body_flag = atom->body_flag; body_flag = atom->body_flag;
quat_atom_flag = atom->quat_flag;
theta_flag = quat_flag = 0; theta_flag = quat_flag = 0;
if (line_flag) theta_flag = 1; if (line_flag) theta_flag = 1;
if (ellipsoid_flag || tri_flag || body_flag) quat_flag = 1; if (ellipsoid_flag || tri_flag || body_flag || quat_atom_flag) quat_flag = 1;
extra_flag = 0; extra_flag = 0;
if (omega_flag || angmom_flag || theta_flag || quat_flag) extra_flag = 1; if (omega_flag || angmom_flag || theta_flag || quat_flag) extra_flag = 1;
@ -329,7 +330,7 @@ FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) :
} }
} }
if (quat_flag) { if (quat_flag && !quat_atom_flag) {
double *quat; double *quat;
for (int i = 0; i < nlocal; i++) { for (int i = 0; i < nlocal; i++) {
quat = nullptr; quat = nullptr;
@ -349,6 +350,16 @@ FixMove::FixMove(LAMMPS *lmp, int narg, char **arg) :
} else } else
qoriginal[i][0] = qoriginal[i][1] = qoriginal[i][2] = qoriginal[i][3] = 0.0; qoriginal[i][0] = qoriginal[i][1] = qoriginal[i][2] = qoriginal[i][3] = 0.0;
} }
} else if (quat_atom_flag) {
double **quat_atom = atom->quat;
for (int i = 0; i < nlocal; i++) {
if (mask[i] & groupbit) {
qoriginal[i][0] = quat_atom[i][0];
qoriginal[i][1] = quat_atom[i][1];
qoriginal[i][2] = quat_atom[i][2];
qoriginal[i][3] = quat_atom[i][3];
}
}
} }
// nrestart = size of per-atom restart data // nrestart = size of per-atom restart data
@ -521,6 +532,7 @@ void FixMove::initial_integrate(int /*vflag*/)
double *radius = atom->radius; double *radius = atom->radius;
double *rmass = atom->rmass; double *rmass = atom->rmass;
double *mass = atom->mass; double *mass = atom->mass;
double **quat_atom = atom->quat;
int *type = atom->type; int *type = atom->type;
int *ellipsoid = atom->ellipsoid; int *ellipsoid = atom->ellipsoid;
int *line = atom->line; int *line = atom->line;
@ -749,9 +761,9 @@ void FixMove::initial_integrate(int /*vflag*/)
avec_line->bonus[atom->line[i]].theta = theta_new; avec_line->bonus[atom->line[i]].theta = theta_new;
} }
// quats for ellipsoids, tris, and bodies // quats for ellipsoids, tris, bodies, and bpm/sphere
if (quat_flag) { if (quat_flag && !quat_atom_flag) {
quat = nullptr; quat = nullptr;
if (ellipsoid_flag && ellipsoid[i] >= 0) if (ellipsoid_flag && ellipsoid[i] >= 0)
quat = avec_ellipsoid->bonus[ellipsoid[i]].quat; quat = avec_ellipsoid->bonus[ellipsoid[i]].quat;
@ -760,6 +772,8 @@ void FixMove::initial_integrate(int /*vflag*/)
else if (body_flag && body[i] >= 0) else if (body_flag && body[i] >= 0)
quat = avec_body->bonus[body[i]].quat; quat = avec_body->bonus[body[i]].quat;
if (quat) MathExtra::quatquat(qrotate, qoriginal[i], quat); if (quat) MathExtra::quatquat(qrotate, qoriginal[i], quat);
} else if (quat_atom_flag) {
MathExtra::quatquat(qrotate, qoriginal[i], quat_atom[i]);
} }
} }
@ -880,9 +894,9 @@ void FixMove::initial_integrate(int /*vflag*/)
avec_line->bonus[atom->line[i]].theta = theta_new; avec_line->bonus[atom->line[i]].theta = theta_new;
} }
// quats for ellipsoids, tris, and bodies // quats for ellipsoids, tris, bodies, and bpm/sphere
if (quat_flag) { if (quat_flag && !quat_atom_flag) {
quat = nullptr; quat = nullptr;
if (ellipsoid_flag && ellipsoid[i] >= 0) if (ellipsoid_flag && ellipsoid[i] >= 0)
quat = avec_ellipsoid->bonus[ellipsoid[i]].quat; quat = avec_ellipsoid->bonus[ellipsoid[i]].quat;
@ -891,6 +905,8 @@ void FixMove::initial_integrate(int /*vflag*/)
else if (body_flag && body[i] >= 0) else if (body_flag && body[i] >= 0)
quat = avec_body->bonus[body[i]].quat; quat = avec_body->bonus[body[i]].quat;
if (quat) MathExtra::quatquat(qrotate, qoriginal[i], quat); if (quat) MathExtra::quatquat(qrotate, qoriginal[i], quat);
} else if (quat_atom_flag) {
MathExtra::quatquat(qrotate, qoriginal[i], quat_atom[i]);
} }
} }
@ -1263,6 +1279,7 @@ void FixMove::set_arrays(int i)
double *quat; double *quat;
double **x = atom->x; double **x = atom->x;
double **quat_atom = atom->quat;
imageint *image = atom->image; imageint *image = atom->image;
int *ellipsoid = atom->ellipsoid; int *ellipsoid = atom->ellipsoid;
int *line = atom->line; int *line = atom->line;
@ -1341,9 +1358,9 @@ void FixMove::set_arrays(int i)
toriginal[i] = theta - 0.0; // NOTE: edit this line toriginal[i] = theta - 0.0; // NOTE: edit this line
} }
// quats for ellipsoids, tris, and bodies // quats for ellipsoids, tris, bodies, and bpm/sphere
if (quat_flag) { if (quat_flag & !quat_atom_flag) {
quat = nullptr; quat = nullptr;
if (ellipsoid_flag && ellipsoid[i] >= 0) if (ellipsoid_flag && ellipsoid[i] >= 0)
quat = avec_ellipsoid->bonus[ellipsoid[i]].quat; quat = avec_ellipsoid->bonus[ellipsoid[i]].quat;
@ -1354,6 +1371,11 @@ void FixMove::set_arrays(int i)
if (quat) { if (quat) {
// qoriginal = f(quat,-delta); // NOTE: edit this line // qoriginal = f(quat,-delta); // NOTE: edit this line
} }
} else if (quat_atom_flag) {
// qoriginal[0] = quat_atom[i][0]; // NOTE: edit this line
// qoriginal[1] = quat_atom[i][1]; // NOTE: edit this line
// qoriginal[2] = quat_atom[i][2]; // NOTE: edit this line
// qoriginal[3] = quat_atom[i][3]; // NOTE: edit this line
} }
} }
xoriginal[i][0] -= vx * delta; xoriginal[i][0] -= vx * delta;
@ -1400,7 +1422,7 @@ void FixMove::set_arrays(int i)
// quats for ellipsoids, tris, and bodies // quats for ellipsoids, tris, and bodies
if (quat_flag) { if (quat_flag && !quat_atom_flag) {
quat = nullptr; quat = nullptr;
if (ellipsoid_flag && ellipsoid[i] >= 0) if (ellipsoid_flag && ellipsoid[i] >= 0)
quat = avec_ellipsoid->bonus[ellipsoid[i]].quat; quat = avec_ellipsoid->bonus[ellipsoid[i]].quat;
@ -1411,6 +1433,11 @@ void FixMove::set_arrays(int i)
if (quat) { if (quat) {
// qoriginal = f(quat,-delta); // NOTE: edit this line // qoriginal = f(quat,-delta); // NOTE: edit this line
} }
} else if (quat_atom_flag) {
// qoriginal[0] = quat_atom[i][0]; // NOTE: edit this line
// qoriginal[1] = quat_atom[i][1]; // NOTE: edit this line
// qoriginal[2] = quat_atom[i][2]; // NOTE: edit this line
// qoriginal[3] = quat_atom[i][3]; // NOTE: edit this line
} }
} }
} }

2
src/fix_move.h
View File

@ -61,7 +61,7 @@ class FixMove : public Fix {
int xvar, yvar, zvar, vxvar, vyvar, vzvar; int xvar, yvar, zvar, vxvar, vyvar, vzvar;
int xvarstyle, yvarstyle, zvarstyle, vxvarstyle, vyvarstyle, vzvarstyle; int xvarstyle, yvarstyle, zvarstyle, vxvarstyle, vyvarstyle, vzvarstyle;
int extra_flag, omega_flag, angmom_flag; int extra_flag, omega_flag, angmom_flag;
int radius_flag, ellipsoid_flag, line_flag, tri_flag, body_flag; int radius_flag, ellipsoid_flag, line_flag, tri_flag, body_flag, quat_atom_flag;
int theta_flag, quat_flag; int theta_flag, quat_flag;
int nlevels_respa, nrestart; int nlevels_respa, nrestart;
int time_origin; int time_origin;