From 99ed989715adcb1b2ea0602cc42b57581dc57db5 Mon Sep 17 00:00:00 2001
From: sjplimp <sjplimp@f3b2605a-c512-4ea7-a41b-209d697bcdaa>
Date: Tue, 19 Apr 2011 19:33:05 +0000
Subject: [PATCH] git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@5977
 f3b2605a-c512-4ea7-a41b-209d697bcdaa

---
 src/USER-EWALDN/ewald_n.cpp | 34 +++++++++---------
 src/fix_rigid.cpp           | 70 +++++++++++++++----------------------
 src/math_extra.cpp          |  6 ++--
 3 files changed, 48 insertions(+), 62 deletions(-)

diff --git a/src/USER-EWALDN/ewald_n.cpp b/src/USER-EWALDN/ewald_n.cpp
index affca321c3..d45b90f05d 100644
--- a/src/USER-EWALDN/ewald_n.cpp
+++ b/src/USER-EWALDN/ewald_n.cpp
@@ -353,10 +353,10 @@ void EwaldN::init_coeff_sums()				// sums based on atoms
     }
   }
   if (function[3]) {					// dipole
-    int *type = atom->type, *ntype = type+atom->nlocal;
-    double *dipole = atom->dipole;
-    for (int *i=type; i<ntype; ++i)
-      sum_local[9].x2 += dipole[i[0]]*dipole[i[0]];
+    double **mu = atom->mu;
+    int nlocal = atom->nlocal;
+    for (int i = 0; i < nlocal; i++)
+      sum_local[9].x2 += mu[i][3]*mu[i][3];
   }
   MPI_Allreduce(sum_local, sum, 2*EWALD_MAX_NSUMS, MPI_DOUBLE, MPI_SUM, world);
 }
@@ -413,7 +413,7 @@ void EwaldN::compute_ek()
   complex *cek, zxyz, zxy, cx;
   vector mui;
   double *x = atom->x[0], *xn = x+3*atom->nlocal, *q = atom->q, qi, bi, ci[7];
-  double *dipole = atom->dipole, *mu = atom->mu ? atom->mu[0] : NULL;
+  double *mu = atom->mu ? atom->mu[0] : NULL;
   int i, kx, ky, n = nkvec*nsums, *type = atom->type, tri = domain->triclinic;
   int func[EWALD_NFUNCS];
 
@@ -444,7 +444,8 @@ void EwaldN::compute_ek()
     if (func[2]) memcpy(ci, B+7*type[0], 7*sizeof(double));
     if (func[3]) { 
       memcpy(mui, mu, sizeof(vector)); mu += 3;
-      vec_scalar_mult(mui, dipole[*type]);
+      vec_scalar_mult(mui, mu[3]);
+      mu += 4;
       h = hvec;
     }
     for (k=kvec; k<nk; ++k) {				// compute rho(k)
@@ -482,7 +483,7 @@ void EwaldN::compute_force()
   vector sum[EWALD_MAX_NSUMS], mui;
   complex *cek, zc, zx, zxy;
   double *f = atom->f[0], *fn = f+3*atom->nlocal, *q = atom->q, *t = NULL;
-  double *dipole = atom->dipole, *mu = atom->mu ? atom->mu[0] : NULL;
+  double *mu = atom->mu ? atom->mu[0] : NULL;
   double *ke, c[EWALD_NFUNCS] = {
     8.0*M_PI*qqrd2e*scale/volume, 2.0*M_PI*sqrt(M_PI)/(12.0*volume),
     2.0*M_PI*sqrt(M_PI)/(192.0*volume), 8.0*M_PI*mumurd2e/volume};
@@ -500,8 +501,9 @@ void EwaldN::compute_force()
     cek = cek_global;
     memset(sum, 0, EWALD_MAX_NSUMS*sizeof(vector));
     if (func[3]) { 
-      register double di = dipole[*type]*c[3];
+      register double di = mu[3] * c[3];
       mui[0] = di*(mu++)[0]; mui[1] = di*(mu++)[1]; mui[2] = di*(mu++)[2];
+      mu++;
     }
     for (nh = (h = hvec)+nkvec; h<nh; ++h, ++k) {
       if (ky!=k->y) {					// based on order in
@@ -556,20 +558,19 @@ void EwaldN::compute_force()
   }
 }
 
-
 void EwaldN::compute_surface()
 {
   if (!function[3]) return;
 
   vector sum_local = VECTOR_NULL, sum_total;
-  double *mu = atom->mu ? atom->mu[0] : NULL, *dipole = atom->dipole;
-  int *type = atom->type, *ntype = type+atom->nlocal;
+  double **mu = atom->mu;
+  int nlocal = atom->nlocal;
 
-  for (int *i=type; i<ntype; ++i) {
-    register double di = dipole[i[0]];
-    sum_local[0] += di*(mu++)[0];
-    sum_local[1] += di*(mu++)[1];
-    sum_local[2] += di*(mu++)[2];
+  for (int i=0; i < nlocal; i++) {
+    register double di = mu[i][3];
+    sum_local[0] += di*mu[i][0];
+    sum_local[1] += di*mu[i][1];
+    sum_local[2] += di*mu[i][2];
   }
   MPI_Allreduce(sum_local, sum_total, 3, MPI_DOUBLE, MPI_SUM, world);
 
@@ -579,7 +580,6 @@ void EwaldN::compute_surface()
   energy_self[3] -= virial_self[3];
 }
 
-
 void EwaldN::compute_energy(int eflag)
 {
   energy = 0.0;
diff --git a/src/fix_rigid.cpp b/src/fix_rigid.cpp
index 15cb67dd84..882bb25e8a 100644
--- a/src/fix_rigid.cpp
+++ b/src/fix_rigid.cpp
@@ -590,7 +590,8 @@ void FixRigid::init()
 
   // compute 6 moments of inertia of each body
   // dx,dy,dz = coords relative to center-of-mass
-  
+  // symmetric 3x3 inertia tensor stored in Voigt notation as 6-vector
+
   double dx,dy,dz,rad;
 
   for (ibody = 0; ibody < nbody; ibody++)
@@ -624,16 +625,15 @@ void FixRigid::init()
     sum[ibody][0] += massone * (dy*dy + dz*dz);
     sum[ibody][1] += massone * (dx*dx + dz*dz);
     sum[ibody][2] += massone * (dx*dx + dy*dy);
-    sum[ibody][3] -= massone * dx*dy;
-    sum[ibody][4] -= massone * dy*dz;
-    sum[ibody][5] -= massone * dx*dz;
+    sum[ibody][3] -= massone * dy*dz;
+    sum[ibody][4] -= massone * dx*dz;
+    sum[ibody][5] -= massone * dx*dy;
   }
 
   // extended particles may contribute extra terms to moments of inertia
 
   if (extended) {
-    double ex[3],ey[3],ez[3],idiag[3];
-    double p[3][3],ptrans[3][3],ispace[3][3],itemp[3][3];
+    double ivec[6];
     double *shape,*quatatom;
 
     for (i = 0; i < nlocal; i++) {
@@ -652,19 +652,13 @@ void FixRigid::init()
       if (eflags[i] & INERTIA_ELLIPSOID) {
 	shape = ebonus[ellipsoid[i]].shape;
 	quatatom = ebonus[ellipsoid[i]].quat;
-	MathExtra::quat_to_mat(quatatom,p);
-	MathExtra::quat_to_mat_trans(quatatom,ptrans);
-	idiag[0] = 0.2*massone * (shape[1]*shape[1]+shape[2]*shape[2]);
-	idiag[1] = 0.2*massone * (shape[0]*shape[0]+shape[2]*shape[2]);
-	idiag[2] = 0.2*massone * (shape[0]*shape[0]+shape[1]*shape[1]);
-	MathExtra::diag_times3(idiag,ptrans,itemp);
-	MathExtra::times3(p,itemp,ispace);
-	sum[ibody][0] += ispace[0][0];
-	sum[ibody][1] += ispace[1][1];
-	sum[ibody][2] += ispace[2][2];
-	sum[ibody][3] += ispace[0][1];
-	sum[ibody][4] += ispace[1][2];
-	sum[ibody][5] += ispace[0][2];
+	MathExtra::inertia_ellipsoid(shape,quatatom,massone,ivec);
+	sum[ibody][0] += ivec[0];
+	sum[ibody][1] += ivec[1];
+	sum[ibody][2] += ivec[2];
+	sum[ibody][3] += ivec[3];
+	sum[ibody][4] += ivec[4];
+	sum[ibody][5] += ivec[5];
       }
     }
   }
@@ -674,18 +668,17 @@ void FixRigid::init()
   // inertia = 3 eigenvalues = principal moments of inertia
   // ex_space,ey_space,ez_space = 3 eigenvectors = principal axes of rigid body
   
-  double tensor[3][3],evectors[3][3];
-
   int ierror;
   double ez0,ez1,ez2;
+  double tensor[3][3],evectors[3][3];
 
   for (ibody = 0; ibody < nbody; ibody++) {
     tensor[0][0] = all[ibody][0];
     tensor[1][1] = all[ibody][1];
     tensor[2][2] = all[ibody][2];
-    tensor[0][1] = tensor[1][0] = all[ibody][3];
-    tensor[1][2] = tensor[2][1] = all[ibody][4];
-    tensor[0][2] = tensor[2][0] = all[ibody][5];
+    tensor[1][2] = tensor[2][1] = all[ibody][3];
+    tensor[0][2] = tensor[2][0] = all[ibody][4];
+    tensor[0][1] = tensor[1][0] = all[ibody][5];
 
     ierror = MathExtra::jacobi(tensor,inertia[ibody],evectors);
     if (ierror) error->all("Insufficient Jacobi rotations for rigid body");
@@ -818,14 +811,13 @@ void FixRigid::init()
       (displace[i][0]*displace[i][0] + displace[i][2]*displace[i][2]);
     sum[ibody][2] += massone * 
       (displace[i][0]*displace[i][0] + displace[i][1]*displace[i][1]);
-    sum[ibody][3] -= massone * displace[i][0]*displace[i][1];
-    sum[ibody][4] -= massone * displace[i][1]*displace[i][2];
-    sum[ibody][5] -= massone * displace[i][0]*displace[i][2];
+    sum[ibody][3] -= massone * displace[i][1]*displace[i][2];
+    sum[ibody][4] -= massone * displace[i][0]*displace[i][2];
+    sum[ibody][5] -= massone * displace[i][0]*displace[i][1];
   }
 
   if (extended) {
-    double ex[3],ey[3],ez[3],idiag[3];
-    double p[3][3],ptrans[3][3],ispace[3][3],itemp[3][3];
+    double ivec[6];
     double *shape;
 
     for (i = 0; i < nlocal; i++) {
@@ -843,19 +835,13 @@ void FixRigid::init()
       }
       if (eflags[i] & INERTIA_ELLIPSOID) {
 	shape = ebonus[ellipsoid[i]].shape;
-	MathExtra::quat_to_mat(qorient[i],p);
-	MathExtra::quat_to_mat_trans(qorient[i],ptrans);
-	idiag[0] = 0.2*massone * (shape[1]*shape[1]+shape[2]*shape[2]);
-	idiag[1] = 0.2*massone * (shape[0]*shape[0]+shape[2]*shape[2]);
-	idiag[2] = 0.2*massone * (shape[0]*shape[0]+shape[1]*shape[1]);
-	MathExtra::diag_times3(idiag,ptrans,itemp);
-	MathExtra::times3(p,itemp,ispace);
-	sum[ibody][0] += ispace[0][0];
-	sum[ibody][1] += ispace[1][1];
-	sum[ibody][2] += ispace[2][2];
-	sum[ibody][3] += ispace[0][1];
-	sum[ibody][4] += ispace[1][2];
-	sum[ibody][5] += ispace[0][2];
+	MathExtra::inertia_ellipsoid(shape,qorient[i],massone,ivec);
+	sum[ibody][0] += ivec[0];
+	sum[ibody][1] += ivec[1];
+	sum[ibody][2] += ivec[2];
+	sum[ibody][3] += ivec[3];
+	sum[ibody][4] += ivec[4];
+	sum[ibody][5] += ivec[5];
       }
     }
   }
diff --git a/src/math_extra.cpp b/src/math_extra.cpp
index 2f25a1f6c3..5bbf3fdca1 100644
--- a/src/math_extra.cpp
+++ b/src/math_extra.cpp
@@ -421,9 +421,9 @@ void inertia_triangle(double *v0, double *v1, double *v2,
 
   // NOTE: use mass
 
-  inertia[0] = inv24*a * (sum-vtsv[0][0]);
-  inertia[1] = inv24*a * (sum-vtsv[1][1]);
-  inertia[2] = inv24*a * (sum-vtsv[2][2]);
+  inertia[0] = inv24*a*(sum-vtsv[0][0]);
+  inertia[1] = inv24*a*(sum-vtsv[1][1]);
+  inertia[2] = inv24*a*(sum-vtsv[2][2]);
   inertia[3] = -inv24*a*vtsv[1][2];
   inertia[4] = -inv24*a*vtsv[0][2];
   inertia[5] = -inv24*a*vtsv[0][1];