From bce09af7a36669fccdef9265eaac9968973a9067 Mon Sep 17 00:00:00 2001 From: sjplimp Date: Thu, 5 May 2011 16:46:57 +0000 Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@6077 f3b2605a-c512-4ea7-a41b-209d697bcdaa --- src/ASPHERE/pair_resquared.cpp | 6 +++--- src/fix_rigid_nve.cpp | 20 +++++++++---------- src/fix_rigid_nvt.cpp | 20 +++++++++---------- src/math_extra.h | 36 ++++------------------------------ 4 files changed, 27 insertions(+), 55 deletions(-) diff --git a/src/ASPHERE/pair_resquared.cpp b/src/ASPHERE/pair_resquared.cpp index 8241e1a6c7..6831f7bd20 100755 --- a/src/ASPHERE/pair_resquared.cpp +++ b/src/ASPHERE/pair_resquared.cpp @@ -503,8 +503,8 @@ void PairRESquared::precompute_i(const int i,RE2Vars &ws) int *ellipsoid = atom->ellipsoid; AtomVecEllipsoid::Bonus *bonus = avec->bonus; MathExtra::quat_to_mat_trans(bonus[ellipsoid[i]].quat,ws.A); - MathExtra::transpose_times_diag3(ws.A,well[atom->type[i]],ws.aTe); - MathExtra::transpose_times_diag3(ws.A,shape2[atom->type[i]],aTs); + MathExtra::transpose_diag3(ws.A,well[atom->type[i]],ws.aTe); + MathExtra::transpose_diag3(ws.A,shape2[atom->type[i]],aTs); MathExtra::diag_times3(shape2[atom->type[i]],ws.A,ws.sa); MathExtra::times3(aTs,ws.A,ws.gamma); MathExtra::rotation_generator_x(ws.A,ws.lA[0]); @@ -885,7 +885,7 @@ double PairRESquared::resquared_lj(const int i, const int j, scorrect[0] = scorrect[0] * scorrect[0] / 2.0; scorrect[1] = scorrect[1] * scorrect[1] / 2.0; scorrect[2] = scorrect[2] * scorrect[2] / 2.0; - MathExtra::transpose_times_diag3(wi.A,scorrect,aTs); + MathExtra::transpose_diag3(wi.A,scorrect,aTs); MathExtra::times3(aTs,wi.A,gamma); for (int ii=0; ii<3; ii++) MathExtra::times3(aTs,wi.lA[ii],lAtwo[ii]); diff --git a/src/fix_rigid_nve.cpp b/src/fix_rigid_nve.cpp index abdb258a75..3b539f2533 100644 --- a/src/fix_rigid_nve.cpp +++ b/src/fix_rigid_nve.cpp @@ -58,8 +58,8 @@ void FixRigidNVE::setup(int vflag) double mbody[3]; for (int ibody = 0; ibody < nbody; ibody++) { - MathExtra::matvec_rows(ex_space[ibody],ey_space[ibody],ez_space[ibody], - angmom[ibody],mbody); + MathExtra::transpose_matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody], + angmom[ibody],mbody); MathExtra::quatvec(quat[ibody],mbody,conjqm[ibody]); conjqm[ibody][0] *= 2.0; conjqm[ibody][1] *= 2.0; @@ -99,8 +99,8 @@ void FixRigidNVE::initial_integrate(int vflag) torque[ibody][1] *= tflag[ibody][1]; torque[ibody][2] *= tflag[ibody][2]; - MathExtra::matvec_rows(ex_space[ibody],ey_space[ibody],ez_space[ibody], - torque[ibody],tbody); + MathExtra::transpose_matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody], + torque[ibody],tbody); MathExtra::quatvec(quat[ibody],tbody,fquat); conjqm[ibody][0] += dtf2 * fquat[0]; @@ -123,8 +123,8 @@ void FixRigidNVE::initial_integrate(int vflag) MathExtra::q_to_exyz(quat[ibody],ex_space[ibody],ey_space[ibody], ez_space[ibody]); MathExtra::invquatvec(quat[ibody],conjqm[ibody],mbody); - MathExtra::matvec_cols(ex_space[ibody],ey_space[ibody],ez_space[ibody], - mbody,angmom[ibody]); + MathExtra::matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody], + mbody,angmom[ibody]); angmom[ibody][0] *= 0.5; angmom[ibody][1] *= 0.5; @@ -252,8 +252,8 @@ void FixRigidNVE::final_integrate() torque[ibody][1] *= tflag[ibody][1]; torque[ibody][2] *= tflag[ibody][2]; - MathExtra::matvec_rows(ex_space[ibody],ey_space[ibody],ez_space[ibody], - torque[ibody],tbody); + MathExtra::transpose_matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody], + torque[ibody],tbody); MathExtra::quatvec(quat[ibody],tbody,fquat); conjqm[ibody][0] += dtf2 * fquat[0]; @@ -262,8 +262,8 @@ void FixRigidNVE::final_integrate() conjqm[ibody][3] += dtf2 * fquat[3]; MathExtra::invquatvec(quat[ibody],conjqm[ibody],mbody); - MathExtra::matvec_cols(ex_space[ibody],ey_space[ibody],ez_space[ibody], - mbody,angmom[ibody]); + MathExtra::matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody], + mbody,angmom[ibody]); angmom[ibody][0] *= 0.5; angmom[ibody][1] *= 0.5; diff --git a/src/fix_rigid_nvt.cpp b/src/fix_rigid_nvt.cpp index 816ea19a37..4ad18e1e7a 100644 --- a/src/fix_rigid_nvt.cpp +++ b/src/fix_rigid_nvt.cpp @@ -164,8 +164,8 @@ void FixRigidNVT::setup(int vflag) double mbody[3]; for (int ibody = 0; ibody < nbody; ibody++) { - MathExtra::matvec_rows(ex_space[ibody],ey_space[ibody],ez_space[ibody], - angmom[ibody],mbody); + MathExtra::transpose_matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody], + angmom[ibody],mbody); MathExtra::quatvec(quat[ibody],mbody,conjqm[ibody]); conjqm[ibody][0] *= 2.0; conjqm[ibody][1] *= 2.0; @@ -225,8 +225,8 @@ void FixRigidNVT::initial_integrate(int vflag) // step 1.3 - apply torque (body coords) to quaternion momentum - MathExtra::matvec_rows(ex_space[ibody],ey_space[ibody],ez_space[ibody], - torque[ibody],tbody); + MathExtra::transpose_matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody], + torque[ibody],tbody); MathExtra::quatvec(quat[ibody],tbody,fquat); conjqm[ibody][0] += dtf2 * fquat[0]; @@ -253,8 +253,8 @@ void FixRigidNVT::initial_integrate(int vflag) MathExtra::q_to_exyz(quat[ibody],ex_space[ibody],ey_space[ibody], ez_space[ibody]); MathExtra::invquatvec(quat[ibody],conjqm[ibody],mbody); - MathExtra::matvec_cols(ex_space[ibody],ey_space[ibody],ez_space[ibody], - mbody,angmom[ibody]); + MathExtra::matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody], + mbody,angmom[ibody]); angmom[ibody][0] *= 0.5; angmom[ibody][1] *= 0.5; @@ -398,8 +398,8 @@ void FixRigidNVT::final_integrate() // convert torque to the body frame - MathExtra::matvec_rows(ex_space[ibody],ey_space[ibody],ez_space[ibody], - torque[ibody],tbody); + MathExtra::transpose_matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody], + torque[ibody],tbody); // compute "force" for quaternion @@ -416,8 +416,8 @@ void FixRigidNVT::final_integrate() // then convert to the space-fixed frame MathExtra::invquatvec(quat[ibody],conjqm[ibody],mbody); - MathExtra::matvec_cols(ex_space[ibody],ey_space[ibody],ez_space[ibody], - mbody,angmom[ibody]); + MathExtra::matvec(ex_space[ibody],ey_space[ibody],ez_space[ibody], + mbody,angmom[ibody]); angmom[ibody][0] *= 0.5; angmom[ibody][1] *= 0.5; diff --git a/src/math_extra.h b/src/math_extra.h index 1e05c5d728..a5fdc81ceb 100755 --- a/src/math_extra.h +++ b/src/math_extra.h @@ -63,8 +63,8 @@ namespace MathExtra { inline void transpose_matvec(const double *ex, const double *ey, const double *ez, const double *v, double *ans); - inline void transpose_times_diag3(const double mat[3][3], const double*vec, - double ans[3][3]); + inline void transpose_diag3(const double mat[3][3], const double*vec, + double ans[3][3]); inline void vecmat(const double *v, const double m[3][3], double *ans); inline void scalar_times3(const double f, double m[3][3]); @@ -91,10 +91,6 @@ namespace MathExtra { inline void invquatvec(double *a, double *b, double *c); inline void axisangle_to_quat(const double *v, const double angle, double *quat); - inline void matvec_rows(double *x, double *y, double *z, - double *b, double *c); - inline void matvec_cols(double *x, double *y, double *z, - double *b, double *c); void angmom_to_omega(double *m, double *ex, double *ey, double *ez, double *idiag, double *w); @@ -418,8 +414,8 @@ void MathExtra::transpose_matvec(const double *ex, const double *ey, transposed matrix times diagonal matrix ------------------------------------------------------------------------- */ -void MathExtra::transpose_times_diag3(const double m[3][3], - const double *d, double ans[3][3]) +void MathExtra::transpose_diag3(const double m[3][3], const double *d, + double ans[3][3]) { ans[0][0] = m[0][0]*d[0]; ans[0][1] = m[1][0]*d[1]; @@ -562,30 +558,6 @@ void MathExtra::axisangle_to_quat(const double *v, const double angle, quat[3] = v[2]*sina; } -/* ---------------------------------------------------------------------- - matvec_rows: c = Ab, where rows of A are x, y, z -------------------------------------------------------------------------- */ - -void MathExtra::matvec_rows(double *x, double *y, double *z, - double *b, double *c) -{ - c[0] = x[0]*b[0] + x[1]*b[1] + x[2]*b[2]; - c[1] = y[0]*b[0] + y[1]*b[1] + y[2]*b[2]; - c[2] = z[0]*b[0] + z[1]*b[1] + z[2]*b[2]; -} - -/* ---------------------------------------------------------------------- - matvec_cols: c = Ab, where columns of A are x, y, z -------------------------------------------------------------------------- */ - -void MathExtra::matvec_cols(double *x, double *y, double *z, - double *b, double *c) -{ - c[0] = x[0]*b[0] + y[0]*b[1] + z[0]*b[2]; - c[1] = x[1]*b[0] + y[1]*b[1] + z[1]*b[2]; - c[2] = x[2]*b[0] + y[2]*b[1] + z[2]*b[2]; -} - /* ---------------------------------------------------------------------- Apply principal rotation generator about x to rotation matrix m ------------------------------------------------------------------------- */