jianmu
/
lammps
forked from lijiext/lammps
1
0
Fork 0
Code Issues Pull Requests Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@6077 f3b2605a-c512-4ea7-a41b-209d697bcdaa

This commit is contained in:
sjplimp 2011-05-05 16:46:57 +00:00
parent b1d3687d40
commit bce09af7a3
4 changed files with 27 additions and 55 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

6
src/ASPHERE/pair_resquared.cpp
View File

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

20
src/fix_rigid_nve.cpp
View File

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

20
src/fix_rigid_nvt.cpp
View File

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

36
src/math_extra.h
View File

@ -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
------------------------------------------------------------------------- */