From: toshinagata1964 <toshinagata1964@a2be9bc6-48de-4e38-9406-05402d4bc13c>
Date: Wed, 7 Dec 2011 11:07:33 +0000 (+0000)
Subject: The structure Mat33 and Transform are now column-first arrangement (no change on... 
X-Git-Tag: v1.0.2~473
X-Git-Url: http://git.osdn.net/view?a=commitdiff_plain;h=0ff736f4965794c17ba6c634fb12e1d6a6db4d95;p=molby%2FMolby.git

The structure Mat33 and Transform are now column-first arrangement (no change on Ruby interface and mbsf file format)

git-svn-id: svn+ssh://svn.sourceforge.jp/svnroot/molby/trunk@161 a2be9bc6-48de-4e38-9406-05402d4bc13c
---

diff --git a/MolLib/MD/MDCore.c b/MolLib/MD/MDCore.c
index b5ba1b0..9691ea6 100644
--- a/MolLib/MD/MDCore.c
+++ b/MolLib/MD/MDCore.c
@@ -165,13 +165,13 @@ md_fit_coordinates(MDArena *arena, Int refno, Double *weights, Transform trans)
 		VecSub(v1, ap[i].r, org1);
 		VecSub(v2, rvf[3 * i], org2);
 		r[0] += w1 * v2.x * v1.x;
-		r[1] += w1 * v2.x * v1.y;
-		r[2] += w1 * v2.x * v1.z;
-		r[3] += w1 * v2.y * v1.x;
+		r[1] += w1 * v2.y * v1.x;
+		r[2] += w1 * v2.z * v1.x;
+		r[3] += w1 * v2.x * v1.y;
 		r[4] += w1 * v2.y * v1.y;
-		r[5] += w1 * v2.y * v1.z;
-		r[6] += w1 * v2.z * v1.x;
-		r[7] += w1 * v2.z * v1.y;
+		r[5] += w1 * v2.z * v1.y;
+		r[6] += w1 * v2.x * v1.z;
+		r[7] += w1 * v2.y * v1.z;
 		r[8] += w1 * v2.z * v1.z;
 /*		for (j = 0; j < 3; j++) {
 			for (k = 0; k < 3; k++) {
@@ -189,7 +189,7 @@ md_fit_coordinates(MDArena *arena, Int refno, Double *weights, Transform trans)
 	for (i = 0; i < 3; i++) {
 		for (j = 0; j < 3; j++) {
 			for (k = 0; k < 3; k++) {
-				q[i*3+j] += r[i*3+k] * r[j*3+k];
+				q[j*3+i] += r[k*3+i] * r[k*3+j];
 			}
 		}
 	}
@@ -232,13 +232,13 @@ md_fit_coordinates(MDArena *arena, Int refno, Double *weights, Transform trans)
 */
 	for (k = 0; k < 3; k++) {
 		u[0] += s[k].x * eigen_vec[k].x;
-		u[1] += s[k].y * eigen_vec[k].x;
-		u[2] += s[k].z * eigen_vec[k].x;
-		u[3] += s[k].x * eigen_vec[k].y;
+		u[1] += s[k].x * eigen_vec[k].y;
+		u[2] += s[k].x * eigen_vec[k].z;
+		u[3] += s[k].y * eigen_vec[k].x;
 		u[4] += s[k].y * eigen_vec[k].y;
-		u[5] += s[k].z * eigen_vec[k].y;
-		u[6] += s[k].x * eigen_vec[k].z;
-		u[7] += s[k].y * eigen_vec[k].z;
+		u[5] += s[k].y * eigen_vec[k].z;
+		u[6] += s[k].z * eigen_vec[k].x;
+		u[7] += s[k].z * eigen_vec[k].y;
 		u[8] += s[k].z * eigen_vec[k].z;
 	}
 	
diff --git a/MolLib/MD/MDPressure.c b/MolLib/MD/MDPressure.c
index 5647bed..98f757c 100644
--- a/MolLib/MD/MDPressure.c
+++ b/MolLib/MD/MDPressure.c
@@ -147,7 +147,7 @@ pressure_control(MDArena *arena)
 		memmove(tf0, tf, sizeof(tf));
 		w0 = (md_rand() - 0.5) * pressure->trial_width;
 		for (i = 0; i < 3; i++) {
-			static char idx[] = {4, 8, 5, 7, 0, 8, 2, 6, 0, 4, 1, 3};
+			static char idx[] = {4, 8, 7, 5, 0, 8, 6, 2, 0, 4, 3, 1};
 			flex = pressure->cell_flexibility[i + 3];
 			if (flex > 0) {
 				w = w0 = (md_rand() + md_rand() - 1.0) * pressure->trial_width * 0.5;
@@ -172,7 +172,7 @@ pressure_control(MDArena *arena)
 				w = w0;
 			} else w = 0;
 			for (j = 0; j < 3; j++) {
-				tf[j * 3 + i] *= 1 + w * fabs(flex);
+				tf[i * 3 + j] *= 1 + w * fabs(flex);
 			}
 		}
 		
diff --git a/MolLib/MD/MDSurface.c b/MolLib/MD/MDSurface.c
index beffc19..eab06eb 100644
--- a/MolLib/MD/MDSurface.c
+++ b/MolLib/MD/MDSurface.c
@@ -788,30 +788,30 @@ s_calc_gradient_vectors(gradient_record *gp, const MDArena *arena, const crossin
 		w4 = 1.0 / (w2 * w2sq);
 		w5 = w3 * w4;
 		gp->mx[0] = csr1->vij.z * w5 * (-wxx * w2 + wyy * wdd);
-		gp->mx[3] = -csr1->vij.x * csr1->vij.y * csr1->vij.z * w5 * (w2 + wdd);
+		gp->mx[1] = -csr1->vij.x * csr1->vij.y * csr1->vij.z * w5 * (w2 + wdd);
 		w6 = csr1->vij.z * csr1->vij.z * w2 * w5;
-		gp->mx[6] = -csr1->vij.x * w6;
-		gp->my[0] = gp->mx[3];
-		gp->my[3] = csr1->vij.z * w5 * (-wyy * w2 + wxx * wdd);
-		gp->my[6] = -csr1->vij.y * w6;
+		gp->mx[2] = -csr1->vij.x * w6;
+		gp->my[0] = gp->mx[1];
+		gp->my[1] = csr1->vij.z * w5 * (-wyy * w2 + wxx * wdd);
+		gp->my[2] = -csr1->vij.y * w6;
 		w6 = w2sq * w3;
 		gp->mz[0] = csr1->vij.x * w6;
-		gp->mz[3] = csr1->vij.y * w6;
-		gp->mz[6] = csr1->vij.z * w6;
-		gp->mx[1] = csr1->vij.x * csr1->vij.y * w4;
+		gp->mz[1] = csr1->vij.y * w6;
+		gp->mz[2] = csr1->vij.z * w6;
+		gp->mx[3] = csr1->vij.x * csr1->vij.y * w4;
 		gp->mx[4] = wyy * w4;
-		gp->my[1] = -wxx * w4;
-		gp->my[4] = -gp->mx[1];
-		gp->mx[7] = gp->my[7] = 0.0;
-		gp->mz[1] = gp->mz[4] = gp->mz[7] = 0.0;
-		gp->mx[2] = (wyy + wzz) * w3;
-		gp->mx[5] = -csr1->vij.x * csr1->vij.y * w3;
+		gp->my[3] = -wxx * w4;
+		gp->my[4] = -gp->mx[3];
+		gp->mx[5] = gp->my[5] = 0.0;
+		gp->mz[3] = gp->mz[4] = gp->mz[5] = 0.0;
+		gp->mx[6] = (wyy + wzz) * w3;
+		gp->mx[7] = -csr1->vij.x * csr1->vij.y * w3;
 		gp->mx[8] = -csr1->vij.x * csr1->vij.z * w3;
-		gp->my[2] = gp->mx[5];
-		gp->my[5] = (wxx + wzz) * w3;
+		gp->my[6] = gp->mx[7];
+		gp->my[7] = (wxx + wzz) * w3;
 		gp->my[8] = -csr1->vij.y * csr1->vij.z * w3;
-		gp->mz[2] = gp->mx[8];
-		gp->mz[5] = gp->my[8];
+		gp->mz[6] = gp->mx[8];
+		gp->mz[7] = gp->my[8];
 		gp->mz[8] = (wxx + wyy) * w3;
 	}
 	
@@ -821,15 +821,15 @@ s_calc_gradient_vectors(gradient_record *gp, const MDArena *arena, const crossin
 	 *  Also for gvky_xj, gvkz_xj  */
 	va1 = csr2->vij;
 	w1 = csr2->gj / csr2->dij;
-	gvkx_xj.x = w1 * (gp->mx[0] * va1.x + gp->mx[3] * va1.y + gp->mx[6] * va1.z);
-	gvkx_xj.y = w1 * (gp->my[0] * va1.x + gp->my[3] * va1.y + gp->my[6] * va1.z);
-	gvkx_xj.z = w1 * (gp->mz[0] * va1.x + gp->mz[3] * va1.y + gp->mz[6] * va1.z);
-	gvky_xj.x = w1 * (gp->mx[1] * va1.x + gp->mx[4] * va1.y + gp->mx[7] * va1.z);
-	gvky_xj.y = w1 * (gp->my[1] * va1.x + gp->my[4] * va1.y);
+	gvkx_xj.x = w1 * (gp->mx[0] * va1.x + gp->mx[1] * va1.y + gp->mx[2] * va1.z);
+	gvkx_xj.y = w1 * (gp->my[0] * va1.x + gp->my[1] * va1.y + gp->my[2] * va1.z);
+	gvkx_xj.z = w1 * (gp->mz[0] * va1.x + gp->mz[1] * va1.y + gp->mz[2] * va1.z);
+	gvky_xj.x = w1 * (gp->mx[3] * va1.x + gp->mx[4] * va1.y + gp->mx[5] * va1.z);
+	gvky_xj.y = w1 * (gp->my[3] * va1.x + gp->my[4] * va1.y);
 	gvky_xj.z = 0.0;
-	gvkz_xj.x = w1 * (gp->mx[2] * va1.x + gp->mx[5] * va1.y + gp->mx[8] * va1.z);
-	gvkz_xj.y = w1 * (gp->my[2] * va1.x + gp->my[5] * va1.y + gp->my[8] * va1.z);
-	gvkz_xj.z = w1 * (gp->mz[2] * va1.x + gp->mz[5] * va1.y + gp->mz[8] * va1.z);
+	gvkz_xj.x = w1 * (gp->mx[6] * va1.x + gp->mx[7] * va1.y + gp->mx[8] * va1.z);
+	gvkz_xj.y = w1 * (gp->my[6] * va1.x + gp->my[7] * va1.y + gp->my[8] * va1.z);
+	gvkz_xj.z = w1 * (gp->mz[6] * va1.x + gp->mz[7] * va1.y + gp->mz[8] * va1.z);
 
 #if 0
 	/*  Intermediate vectors: gradients of gvk.x, gvk.y, gvk.y by xj  */
@@ -1318,37 +1318,37 @@ s_calc_surface_area(MDArena *arena)
 				w3 = -sin(ir1->s) * sinomega_inv;
 				w4 = cos(ir1->s) * sinomega_inv;
 				grad.a_xj.x += w2 * grad.s_xj.x
-					- w3 * (grad.mx[0]*p2.x + grad.mx[3]*p2.y + grad.mx[6]*p2.z)
-					- w4 * (grad.mx[1]*p2.x + grad.mx[4]*p2.y + grad.mx[7]*p2.z);
+					- w3 * (grad.mx[0]*p2.x + grad.mx[1]*p2.y + grad.mx[2]*p2.z)
+					- w4 * (grad.mx[3]*p2.x + grad.mx[4]*p2.y + grad.mx[5]*p2.z);
 				grad.a_xj.y += w2 * grad.s_xj.y
-					- w3 * (grad.my[0]*p2.x + grad.my[3]*p2.y + grad.my[6]*p2.z)
-					- w4 * (grad.my[1]*p2.x + grad.my[4]*p2.y + grad.my[7]*p2.z);
+					- w3 * (grad.my[0]*p2.x + grad.my[1]*p2.y + grad.my[2]*p2.z)
+					- w4 * (grad.my[3]*p2.x + grad.my[4]*p2.y + grad.my[5]*p2.z);
 				grad.a_xj.z += w2 * grad.s_xj.z
-					- w3 * (grad.mz[0]*p2.x + grad.mz[3]*p2.y + grad.mz[6]*p2.z)
-					- w4 * (grad.mz[1]*p2.x + grad.mz[4]*p2.y + grad.mz[7]*p2.z);
+					- w3 * (grad.mz[0]*p2.x + grad.mz[1]*p2.y + grad.mz[2]*p2.z)
+					- w4 * (grad.mz[3]*p2.x + grad.mz[4]*p2.y + grad.mz[5]*p2.z);
 				VecScaleInc(grad.a_xk, grad.s_xk, w2);
 			}
 
 		#if 0
 		{
 			s_sp_debug(arena, "[%d,%d,%d] p1.x = %10.6f, d(p1.x)/dxj = {%10.6f,%10.6f,%10.6f} d(p1.x)/dxk = {%10.6f,%10.6f,%10.6f}", ir1->i+1, ir1->j+1, ir1->k+1, p1.x,
-				(grad.mx[0]*w3+grad.mx[1]*w4)*sinomega - pw1.x/csr1->rj*grad.s_xj.x,
-				(grad.my[0]*w3+grad.my[1]*w4)*sinomega - pw1.x/csr1->rj*grad.s_xj.y,
-				(grad.mz[0]*w3+grad.mz[1]*w4)*sinomega - pw1.x/csr1->rj*grad.s_xj.z,
+				(grad.mx[0]*w3+grad.mx[3]*w4)*sinomega - pw1.x/csr1->rj*grad.s_xj.x,
+				(grad.my[0]*w3+grad.my[3]*w4)*sinomega - pw1.x/csr1->rj*grad.s_xj.y,
+				(grad.mz[0]*w3+grad.mz[3]*w4)*sinomega - pw1.x/csr1->rj*grad.s_xj.z,
 				-pw1.x/csr1->rj*grad.s_xk.x,
 				-pw1.x/csr1->rj*grad.s_xk.y,
 				-pw1.x/csr1->rj*grad.s_xk.z);
 			s_sp_debug(arena, "[%d,%d,%d] p1.y = %10.6f, d(p1.y)/dxj = {%10.6f,%10.6f,%10.6f} d(p1.y)/dxk = {%10.6f,%10.6f,%10.6f}", ir1->i+1, ir1->j+1, ir1->k+1, p1.y,
-				(grad.mx[3]*w3+grad.mx[4]*w4)*sinomega - pw1.y/csr1->rj*grad.s_xj.x,
-				(grad.my[3]*w3+grad.my[4]*w4)*sinomega - pw1.y/csr1->rj*grad.s_xj.y,
-				(grad.mz[3]*w3+grad.mz[4]*w4)*sinomega - pw1.y/csr1->rj*grad.s_xj.z,
+				(grad.mx[1]*w3+grad.mx[4]*w4)*sinomega - pw1.y/csr1->rj*grad.s_xj.x,
+				(grad.my[1]*w3+grad.my[4]*w4)*sinomega - pw1.y/csr1->rj*grad.s_xj.y,
+				(grad.mz[1]*w3+grad.mz[4]*w4)*sinomega - pw1.y/csr1->rj*grad.s_xj.z,
 				-pw1.y/csr1->rj*grad.s_xk.x,
 				-pw1.y/csr1->rj*grad.s_xk.y,
 				-pw1.y/csr1->rj*grad.s_xk.z);
 			s_sp_debug(arena, "[%d,%d,%d] p1.z = %10.6f, d(p1.z)/dxj = {%10.6f,%10.6f,%10.6f} d(p1.z)/dxk = {%10.6f,%10.6f,%10.6f}", ir1->i+1, ir1->j+1, ir1->k+1, p1.y,
-				(grad.mx[6]*w3+grad.mx[7]*w4)*sinomega - pw1.z/csr1->rj*grad.s_xj.x,
-				(grad.my[6]*w3+grad.my[7]*w4)*sinomega - pw1.z/csr1->rj*grad.s_xj.y,
-				(grad.mz[6]*w3+grad.mz[7]*w4)*sinomega - pw1.z/csr1->rj*grad.s_xj.z,
+				(grad.mx[2]*w3+grad.mx[5]*w4)*sinomega - pw1.z/csr1->rj*grad.s_xj.x,
+				(grad.my[2]*w3+grad.my[5]*w4)*sinomega - pw1.z/csr1->rj*grad.s_xj.y,
+				(grad.mz[2]*w3+grad.mz[5]*w4)*sinomega - pw1.z/csr1->rj*grad.s_xj.z,
 				-pw1.z/csr1->rj*grad.s_xk.x,
 				-pw1.z/csr1->rj*grad.s_xk.y,
 				-pw1.z/csr1->rj*grad.s_xk.z);
@@ -1368,14 +1368,14 @@ s_calc_surface_area(MDArena *arena)
 				w3 = -sin(ir2->s) * sinomega_inv;
 				w4 = cos(ir2->s) * sinomega_inv;
 				grad.a_xk.x += w2 * grad.s_xj.x
-					- w3 * (grad.mx[0]*p1.x + grad.mx[3]*p1.y + grad.mx[6]*p1.z)
-					- w4 * (grad.mx[1]*p1.x + grad.mx[4]*p1.y + grad.mx[7]*p1.z);
+					- w3 * (grad.mx[0]*p1.x + grad.mx[1]*p1.y + grad.mx[2]*p1.z)
+					- w4 * (grad.mx[3]*p1.x + grad.mx[4]*p1.y + grad.mx[5]*p1.z);
 				grad.a_xk.y += w2 * grad.s_xj.y
-					- w3 * (grad.my[0]*p1.x + grad.my[3]*p1.y + grad.my[6]*p1.z)
-					- w4 * (grad.my[1]*p1.x + grad.my[4]*p1.y + grad.my[7]*p1.z);
+					- w3 * (grad.my[0]*p1.x + grad.my[1]*p1.y + grad.my[2]*p1.z)
+					- w4 * (grad.my[3]*p1.x + grad.my[4]*p1.y + grad.my[5]*p1.z);
 				grad.a_xk.z += w2 * grad.s_xj.z
-					- w3 * (grad.mz[0]*p1.x + grad.mz[3]*p1.y + grad.mz[6]*p1.z)
-					- w4 * (grad.mz[1]*p1.x + grad.mz[4]*p1.y + grad.mz[7]*p1.z);
+					- w3 * (grad.mz[0]*p1.x + grad.mz[1]*p1.y + grad.mz[2]*p1.z)
+					- w4 * (grad.mz[3]*p1.x + grad.mz[4]*p1.y + grad.mz[5]*p1.z);
 				VecScaleInc(grad.a_xj, grad.s_xk, w2);
 				SA_DEBUG("OMEGA", ir1);
 			}
@@ -1386,23 +1386,23 @@ s_calc_surface_area(MDArena *arena)
 				-pw2.x/csr2->rj*grad.s_xk.x,
 				-pw2.x/csr2->rj*grad.s_xk.y,
 				-pw2.x/csr2->rj*grad.s_xk.z,
-				(grad.mx[0]*w3+grad.mx[1]*w4)*sinomega - pw2.x/csr2->rj*grad.s_xj.x,
-				(grad.my[0]*w3+grad.my[1]*w4)*sinomega - pw2.x/csr2->rj*grad.s_xj.y,
-				(grad.mz[0]*w3+grad.mz[1]*w4)*sinomega - pw2.x/csr2->rj*grad.s_xj.z);
+				(grad.mx[0]*w3+grad.mx[3]*w4)*sinomega - pw2.x/csr2->rj*grad.s_xj.x,
+				(grad.my[0]*w3+grad.my[3]*w4)*sinomega - pw2.x/csr2->rj*grad.s_xj.y,
+				(grad.mz[0]*w3+grad.mz[3]*w4)*sinomega - pw2.x/csr2->rj*grad.s_xj.z);
 			s_sp_debug(arena, "[%d,%d,%d] p2.y = %10.6f, d(p2.y)/dxj = {%10.6f,%10.6f,%10.6f} d(p2.y)/dxk = {%10.6f,%10.6f,%10.6f}", ir1->i+1, ir1->j+1, ir1->k+1, p2.y,
 				-pw2.y/csr2->rj*grad.s_xk.x,
 				-pw2.y/csr2->rj*grad.s_xk.y,
 				-pw2.y/csr2->rj*grad.s_xk.z,
-				(grad.mx[3]*w3+grad.mx[4]*w4)*sinomega - pw2.y/csr2->rj*grad.s_xj.x,
-				(grad.my[3]*w3+grad.my[4]*w4)*sinomega - pw2.y/csr2->rj*grad.s_xj.y,
-				(grad.mz[3]*w3+grad.mz[4]*w4)*sinomega - pw2.y/csr2->rj*grad.s_xj.z);
+				(grad.mx[1]*w3+grad.mx[4]*w4)*sinomega - pw2.y/csr2->rj*grad.s_xj.x,
+				(grad.my[1]*w3+grad.my[4]*w4)*sinomega - pw2.y/csr2->rj*grad.s_xj.y,
+				(grad.mz[1]*w3+grad.mz[4]*w4)*sinomega - pw2.y/csr2->rj*grad.s_xj.z);
 			s_sp_debug(arena, "[%d,%d,%d] p2.z = %10.6f, d(p2.z)/dxj = {%10.6f,%10.6f,%10.6f} d(p2.z)/dxk = {%10.6f,%10.6f,%10.6f}", ir1->i+1, ir1->j+1, ir1->k+1, p2.z,
 				-pw2.z/csr2->rj*grad.s_xk.x,
 				-pw2.z/csr2->rj*grad.s_xk.y,
 				-pw2.z/csr2->rj*grad.s_xk.z,
-				(grad.mx[6]*w3+grad.mx[7]*w4)*sinomega - pw2.z/csr2->rj*grad.s_xj.x,
-				(grad.my[6]*w3+grad.my[7]*w4)*sinomega - pw2.z/csr2->rj*grad.s_xj.y,
-				(grad.mz[6]*w3+grad.mz[7]*w4)*sinomega - pw2.z/csr2->rj*grad.s_xj.z);
+				(grad.mx[2]*w3+grad.mx[5]*w4)*sinomega - pw2.z/csr2->rj*grad.s_xj.x,
+				(grad.my[2]*w3+grad.my[5]*w4)*sinomega - pw2.z/csr2->rj*grad.s_xj.y,
+				(grad.mz[2]*w3+grad.mz[5]*w4)*sinomega - pw2.z/csr2->rj*grad.s_xj.z);
 		}
 		#endif
 
diff --git a/MolLib/MainView.c b/MolLib/MainView.c
index e45be50..1b34c1f 100755
--- a/MolLib/MainView.c
+++ b/MolLib/MainView.c
@@ -1184,12 +1184,12 @@ drawAtom(MainView *mview, int i1, int selected, const Vector *dragOffset, const
 				MatrixMul(pmat2, mview->mol->cell->rtr, ap->aniso->pmat);
 				MatrixMul(pmat2, *((Mat33 *)trp), pmat2);
 				MatrixMul(pmat2, mview->mol->cell->tr, pmat2);
-				MatrixTranspose(pmat2, pmat2);
+				for (i = 0; i < 9; i++)
+					elip[i] = pmat2[i] * mview->probabilityScale;
 			} else {
-				MatrixTranspose(pmat2, ap->aniso->pmat);
+				for (i = 0; i < 9; i++)
+					elip[i] = ap->aniso->pmat[i] * mview->probabilityScale;
 			}
-			for (i = 0; i < 9; i++)
-				elip[i] = pmat2[i] * mview->probabilityScale;
 		/*	for (i = 0; i < 3; i++) {
 				Double w = ap->aniso->val[i];
 				Vector vv = ap->aniso->axis[i];
@@ -1845,14 +1845,14 @@ MainView_drawModel(MainView *mview)
     glGetDoublev(GL_MODELVIEW_MATRIX, mview->modelview_matrix);
     glGetDoublev(GL_PROJECTION_MATRIX, mview->projection_matrix);
 	pp = mview->modelview_matrix;
-	mtr[0] = pp[0]; mtr[1] = pp[4]; mtr[2] = pp[8];
-	mtr[3] = pp[1]; mtr[4] = pp[5]; mtr[5] = pp[9];
-	mtr[6] = pp[2]; mtr[7] = pp[6]; mtr[8] = pp[10];
+	mtr[0] = pp[0]; mtr[1] = pp[1]; mtr[2] = pp[2];
+	mtr[3] = pp[4]; mtr[4] = pp[5]; mtr[5] = pp[6];
+	mtr[6] = pp[8]; mtr[7] = pp[9]; mtr[8] = pp[10];
 	mtr[9] = pp[12]; mtr[10] = pp[13]; mtr[11] = pp[14];
 	TransformInvert(mtr, mtr);
 	mview->camera.x = mtr[9]; mview->camera.y = mtr[10]; mview->camera.z = mtr[11];
-	mview->lookto.x = mtr[2]; mview->lookto.y = mtr[5]; mview->lookto.z = mtr[8];
-	mview->up.x = mtr[1]; mview->up.y = mtr[4]; mview->up.z = mtr[7];
+	mview->lookto.x = mtr[6]; mview->lookto.y = mtr[7]; mview->lookto.z = mtr[8];
+	mview->up.x = mtr[3]; mview->up.y = mtr[4]; mview->up.z = mtr[5];
 
 	/*  Draw labels  */
 	glDisable(GL_LIGHTING);
diff --git a/MolLib/Molecule.c b/MolLib/Molecule.c
index 8c62292..a8a3a94 100755
--- a/MolLib/Molecule.c
+++ b/MolLib/Molecule.c
@@ -868,9 +868,15 @@ MoleculeLoadMbsfFile(Molecule *mp, const char *fname, char *errbuf, int errbufsi
 					snprintf(errbuf, errbufsize, "line %d: bad symmetry_operation format", lineNumber);
 					goto exit;
 				}
-				tr[i * 3] = dbuf[0];
-				tr[i * 3 + 1] = dbuf[1];
-				tr[i * 3 + 2] = dbuf[2];
+				if (i < 3) {
+					tr[i] = dbuf[0];
+					tr[i + 3] = dbuf[1];
+					tr[i + 6] = dbuf[2];
+				} else {
+					tr[9] = dbuf[0];
+					tr[10] = dbuf[1];
+					tr[11] = dbuf[2];
+				}
 				i++;
 				if (i == 4) {
 					AssignArray(&mp->syms, &mp->nsyms, sizeof(Transform), mp->nsyms, tr);
@@ -886,7 +892,7 @@ MoleculeLoadMbsfFile(Molecule *mp, const char *fname, char *errbuf, int errbufsi
 					continue;
 				if (buf[0] == '\n')
 					break;
-				/* a11 a12 a13; a21 a22 a23; a31 a32 a33; t1 t2 t3 */
+				/* ax ay az; bx by bz; cx cy cz; ox oy oz; fx fy fz */
 				if (sscanf(buf, "%lf %lf %lf", &dbuf[0], &dbuf[1], &dbuf[2]) < 3) {
 					snprintf(errbuf, errbufsize, "line %d: bad symmetry_operation format", lineNumber);
 					goto exit;
@@ -1570,13 +1576,13 @@ sMoleculeSymopStringsToTransform(char **symops, Transform tr)
 				}
 				tr[9 + i] = d * sn;
 			} else if (*symop == 'x' || *symop == 'X') {
-				tr[i * 3] = sn;
+				tr[i] = sn;
 				symop++;
 			} else if (*symop == 'y' || *symop == 'Y') {
-				tr[i * 3 + 1] = sn;
+				tr[i + 3] = sn;
 				symop++;
 			} else if (*symop == 'z' || *symop == 'Z') {
-				tr[i * 3 + 2] = sn;
+				tr[i + 6] = sn;
 				symop++;
 			} else return 1;  /*  Bad format  */
 		} /* end while (*symop != 0) */
@@ -1661,13 +1667,13 @@ MoleculeLoadTepFile(Molecule *mp, const char *fname, char *errbuf, int errbufsiz
 			if (cellType == 0) {
 				ReadFormat(buf, "I1F14F3F3F3F15F3F3F3F15F3F3F3", ibuf, fbuf, fbuf+1, fbuf+2, fbuf+3, fbuf+4, fbuf+5, fbuf+6, fbuf+7, fbuf+8, fbuf+9, fbuf+10, fbuf+11);
 				tr[0] = fbuf[1];
-				tr[1] = fbuf[2];
-				tr[2] = fbuf[3];
-				tr[3] = fbuf[5];
+				tr[3] = fbuf[2];
+				tr[6] = fbuf[3];
+				tr[1] = fbuf[5];
 				tr[4] = fbuf[6];
-				tr[5] = fbuf[7];
-				tr[6] = fbuf[9];
-				tr[7] = fbuf[10];
+				tr[7] = fbuf[7];
+				tr[2] = fbuf[9];
+				tr[5] = fbuf[10];
 				tr[8] = fbuf[11];
 				tr[9] = fbuf[0];
 				tr[10] = fbuf[4];
@@ -1893,8 +1899,8 @@ MoleculeLoadShelxFile(Molecule *mp, const char *fname, char *errbuf, int errbufs
 		static Transform tr_c = {1, 0, 0, 0, 1, 0, 0, 0, 1, 0.5, 0.5, 0};
 		static Transform tr_a = {1, 0, 0, 0, 1, 0, 0, 0, 1, 0, 0.5, 0.5};
 		static Transform tr_b = {1, 0, 0, 0, 1, 0, 0, 0, 1, 0.5, 0, 0.5};
-		static Transform tr_r1 = {0, -1, 0, 1, -1, 0, 0, 0, 1, 0, 0, 0};
-		static Transform tr_r2 = {-1, 1, 0, -1, 0, 0, 0, 0, 1, 0, 0, 0};
+		static Transform tr_r1 = {0, 1, 0, -1, -1, 0, 0, 0, 1, 0, 0, 0};
+		static Transform tr_r2 = {-1, -1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 0};
 		case 1:  /* P */
 			break;
 		case 2:  /* I */
@@ -3505,8 +3511,9 @@ MoleculeWriteToMbsfFile(Molecule *mp, const char *fname, char *errbuf, int errbu
 		fprintf(fp, "! a11 a12 a13; a21 a22 a23; a31 a32 a33; t1 t2 t3\n");
 		for (i = 0; i < mp->nsyms; i++) {
 			Transform *tp = mp->syms + i;
+			const unsigned char s_index_order[12] = {0, 3, 6, 1, 4, 7, 2, 5, 8, 9, 10, 11};
 			for (j = 0; j < 12; j++)
-				fprintf(fp, "%11.6f%c", (*tp)[j], (j % 3 == 2 ? '\n' : ' '));
+				fprintf(fp, "%11.6f%c", (*tp)[s_index_order[j]], (j % 3 == 2 ? '\n' : ' '));
 		}
 		fprintf(fp, "\n");
 	}
@@ -8297,9 +8304,9 @@ MoleculeCalculateCellFromAxes(XtalCell *cp, int calc_abc)
 		vp1 = &(cp->axes[n1]);
 		vp2 = &(cp->axes[n2]);
 		vp3 = &(cp->axes[n3]);
-		cp->tr[n1] = vp1->x;
-		cp->tr[n1 + 3] = vp1->y;
-		cp->tr[n1 + 6] = vp1->z;
+		cp->tr[n1*3] = vp1->x;
+		cp->tr[n1*3+1] = vp1->y;
+		cp->tr[n1*3+2] = vp1->z;
 		cp->tr[9] = cp->origin.x;
 		cp->tr[10] = cp->origin.y;
 		cp->tr[11] = cp->origin.z;
@@ -8328,34 +8335,34 @@ MoleculeCalculateCellFromAxes(XtalCell *cp, int calc_abc)
 				if (NormalizeVec(&v2, &v2))
 					return -1;  /*  Non-regular transform  */
 			}
-			cp->tr[n2] = v1.x;
-			cp->tr[n2 + 3] = v1.y;
-			cp->tr[n2 + 6] = v1.z;
-			cp->tr[n3] = v2.x;
-			cp->tr[n3 + 3] = v2.y;
-			cp->tr[n3 + 6] = v2.z;
+			cp->tr[n2*3] = v1.x;
+			cp->tr[n2*3+1] = v1.y;
+			cp->tr[n2*3+2] = v1.z;
+			cp->tr[n3*3] = v2.x;
+			cp->tr[n3*3+1] = v2.y;
+			cp->tr[n3*3+2] = v2.z;
 		} else {
 			VecCross(v1, *vp1, *vp2);
 			if (fabs(VecDot(v1, *vp3)) < 1e-7)
 				return -1;  /*  Non-regular transform  */
-			cp->tr[n2] = vp2->x;
-			cp->tr[n2 + 3] = vp2->y;
-			cp->tr[n2 + 6] = vp2->z;
-			cp->tr[n3] = vp3->x;
-			cp->tr[n3 + 3] = vp3->y;
-			cp->tr[n3 + 6] = vp3->z;
+			cp->tr[n2*3] = vp2->x;
+			cp->tr[n2*3+1] = vp2->y;
+			cp->tr[n2*3+2] = vp2->z;
+			cp->tr[n3*3] = vp3->x;
+			cp->tr[n3*3+1] = vp3->y;
+			cp->tr[n3*3+2] = vp3->z;
 		}
 	}
 	if (TransformInvert(cp->rtr, cp->tr))
 		return -1;  /*  Non-regular transform  */
 
 	/*  Calculate the reciprocal cell parameters  */
-	cp->rcell[0] = sqrt(cp->rtr[0] * cp->rtr[0] + cp->rtr[3] * cp->rtr[3] + cp->rtr[6] * cp->rtr[6]);
-	cp->rcell[1] = sqrt(cp->rtr[1] * cp->rtr[1] + cp->rtr[4] * cp->rtr[4] + cp->rtr[7] * cp->rtr[7]);
-	cp->rcell[2] = sqrt(cp->rtr[2] * cp->rtr[2] + cp->rtr[5] * cp->rtr[5] + cp->rtr[8] * cp->rtr[8]);
-	cp->rcell[3] = acos((cp->rtr[1] * cp->rtr[2] + cp->rtr[4] * cp->rtr[5] + cp->rtr[7] * cp->rtr[8]) / (cp->rcell[1] * cp->rcell[2])) * kRad2Deg;
-	cp->rcell[4] = acos((cp->rtr[2] * cp->rtr[0] + cp->rtr[5] * cp->rtr[3] + cp->rtr[8] * cp->rtr[6]) / (cp->rcell[2] * cp->rcell[0])) * kRad2Deg;
-	cp->rcell[5] = acos((cp->rtr[0] * cp->rtr[1] + cp->rtr[3] * cp->rtr[4] + cp->rtr[6] * cp->rtr[7]) / (cp->rcell[0] * cp->rcell[1])) * kRad2Deg;
+	cp->rcell[0] = sqrt(cp->rtr[0] * cp->rtr[0] + cp->rtr[1] * cp->rtr[1] + cp->rtr[2] * cp->rtr[2]);
+	cp->rcell[1] = sqrt(cp->rtr[3] * cp->rtr[3] + cp->rtr[4] * cp->rtr[4] + cp->rtr[5] * cp->rtr[5]);
+	cp->rcell[2] = sqrt(cp->rtr[6] * cp->rtr[6] + cp->rtr[7] * cp->rtr[7] + cp->rtr[8] * cp->rtr[8]);
+	cp->rcell[3] = acos((cp->rtr[3] * cp->rtr[6] + cp->rtr[4] * cp->rtr[7] + cp->rtr[5] * cp->rtr[8]) / (cp->rcell[1] * cp->rcell[2])) * kRad2Deg;
+	cp->rcell[4] = acos((cp->rtr[6] * cp->rtr[0] + cp->rtr[7] * cp->rtr[1] + cp->rtr[8] * cp->rtr[2]) / (cp->rcell[2] * cp->rcell[0])) * kRad2Deg;
+	cp->rcell[5] = acos((cp->rtr[0] * cp->rtr[3] + cp->rtr[1] * cp->rtr[4] + cp->rtr[2] * cp->rtr[5]) / (cp->rcell[0] * cp->rcell[1])) * kRad2Deg;
 	
 	if (calc_abc) {
 		/*  Calculate a, b, c, alpha, beta, gamma  */
@@ -8563,9 +8570,9 @@ MoleculeSetAniso(Molecule *mp, int n1, int type, Double x11, Double x22, Double
 		} else {
 			val[u] = 1 / sqrt(val[u]);
 		}
-		anp->pmat[u] = axis[u].x * val[u];
-		anp->pmat[u+3] = axis[u].y * val[u];
-		anp->pmat[u+6] = axis[u].z * val[u];
+		anp->pmat[u*3] = axis[u].x * val[u];
+		anp->pmat[u*3+1] = axis[u].y * val[u];
+		anp->pmat[u*3+2] = axis[u].z * val[u];
 	}
 	__MoleculeUnlock(mp);
 }
diff --git a/MolLib/Molecule.h b/MolLib/Molecule.h
index 2330eb4..86c721e 100755
--- a/MolLib/Molecule.h
+++ b/MolLib/Molecule.h
@@ -43,7 +43,7 @@ typedef struct Aniso {
 	Double  bij[6];    /*  b11, b22, b33, b12, b23, b31 (ORTEP type 0) */
 	char has_bsig;     /*  Has sigma values?  */
 	Double  bsig[6];   /*  sigma values  */
-	Mat33  pmat;      /*  A 3x3 matrix whose three row vectors are the principal axes of the ellipsoid. Note: If the B matrix is not positive definite, the axis length corresponding to the negative eigenvalue is replaced with 0.001.  */
+	Mat33  pmat;      /*  A 3x3 matrix whose three column vectors are the principal axes of the ellipsoid. Note: If the B matrix is not positive definite, the axis length corresponding to the negative eigenvalue is replaced with 0.001.  */
 } Aniso;
 
 /*  Symmetry operation  */
diff --git a/MolLib/Ruby_bind/ruby_types.c b/MolLib/Ruby_bind/ruby_types.c
index a9d056f..2f0c644 100644
--- a/MolLib/Ruby_bind/ruby_types.c
+++ b/MolLib/Ruby_bind/ruby_types.c
@@ -521,8 +521,6 @@ s_Vector3D_Inspect(VALUE self)
 
 #pragma mark ====== Transform Class ======
 
-static int s_index_order[16] = {0, 3, 6, 1, 4, 7, 2, 5, 8, 9, 10, 11};
-
 void
 TransformFromValue(VALUE val, Transform *tp)
 {
@@ -539,7 +537,7 @@ TransformFromValue(VALUE val, Transform *tp)
 		VALUE *valp = RARRAY_PTR(val);
 		if (len == 12) {
 			for (i = 0; i < 12; i++)
-				(*tp)[s_index_order[i]] = NUM2DBL(rb_Float(valp[i]));
+				(*tp)[i] = NUM2DBL(rb_Float(valp[i]));
 			return;
 		} else if (len == 4) {
 			VALUE val2, *valp2;
@@ -551,7 +549,7 @@ TransformFromValue(VALUE val, Transform *tp)
 					goto array_format_error;
 				valp2 = RARRAY_PTR(val2);
 				for (j = 0; j < 3; j++)
-					(*tp)[s_index_order[3 * i + j]] = NUM2DBL(rb_Float(valp2[j]));
+					(*tp)[3 * i + j] = NUM2DBL(rb_Float(valp2[j]));
 			}
 			return;
 		}
@@ -568,12 +566,12 @@ TransformFromValue(VALUE val, Transform *tp)
 			/*  Matrix-type object  */
 			for (i = 0; i < 4; i++) {
 				for (j = 0; j < 3; j++)
-					(*tp)[s_index_order[i * 3 + j]] = NUM2DBL(rb_Float(rb_funcall(val, index_mid, 2, INT2FIX(i), INT2FIX(j))));
+					(*tp)[i * 3 + j] = NUM2DBL(rb_Float(rb_funcall(val, index_mid, 2, INT2FIX(i), INT2FIX(j))));
 			}
 		} else {
 			/*  Other "array-like" object  */
 			for (i = 0; i < 12; i++)
-				(*tp)[s_index_order[i]] = NUM2DBL(rb_Float(rb_funcall(val, index_mid, 1, INT2FIX(i))));
+				(*tp)[i] = NUM2DBL(rb_Float(rb_funcall(val, index_mid, 1, INT2FIX(i))));
 		}
 	}
 }
@@ -677,7 +675,7 @@ s_Transform_NewFromColumns(VALUE klass, VALUE val)
 				w[j] = NUM2DBL(rb_Float(valpp[j]));
 		}
 		for (j = 0; j < 3; j++)
-			tr[s_index_order[i * 3 + j]] = w[j];
+			tr[i * 3 + j] = w[j];
 	}
 	return s_Transform_NewFromTransform(&tr);
 }
@@ -718,7 +716,7 @@ s_Transform_NewFromRows(VALUE klass, VALUE val)
 				w[j] = NUM2DBL(rb_Float(valpp[j]));
 		}
 		for (j = 0; j < 4; j++)
-			tr[s_index_order[j * 3 + i]] = w[j];
+			tr[j * 3 + i] = w[j];
 	}
 	return s_Transform_NewFromTransform(&tr);
 }
@@ -741,7 +739,7 @@ s_Transform_ElementAtIndex(VALUE self, VALUE val1, VALUE val2)
 	Data_Get_Struct(self, Transform, tp);
 	if (n1 < 0 || n1 >= 4 || n2 < 0 || n2 >= 3)
 		rb_raise(rb_eMolbyError, "index to Transform out of range");
-	w = (*tp)[s_index_order[n1 * 3 + n2]];
+	w = (*tp)[n1 * 3 + n2];
 	return rb_float_new(w);
 }
 
@@ -764,7 +762,7 @@ s_Transform_SetElementAtIndex(VALUE self, VALUE idx1, VALUE idx2, VALUE val)
 	if (n1 < 0 || n1 >= 4 || n2 < 0 || n2 >= 3)
 		rb_raise(rb_eMolbyError, "index to Transform out of range");
 	w = NUM2DBL(rb_Float(val));
-	(*tp)[s_index_order[n1 * 3 + n2]] = w;
+	(*tp)[n1 * 3 + n2] = w;
 	return rb_float_new(w);
 }
 
@@ -1029,9 +1027,9 @@ s_Transform_Column(VALUE self, VALUE val)
 	Data_Get_Struct(self, Transform, tp1);
 	if (n < 0 || n >= 4)
 		rb_raise(rb_eMolbyError, "row index out of range");
-	v.x = (*tp1)[n];
-	v.y = (*tp1)[n + 3];
-	v.z = (*tp1)[n + 6];
+	v.x = (*tp1)[n * 3];
+	v.y = (*tp1)[n * 3 + 1];
+	v.z = (*tp1)[n * 3 + 2];
 	return ValueFromVector(&v);
 }
 
@@ -1082,7 +1080,7 @@ s_Transform_ToArray(VALUE self)
 	int i;
 	Data_Get_Struct(self, Transform, tp1);
 	for (i = 0; i < 12; i++)
-		val[i] = rb_float_new((*tp1)[s_index_order[i]]);
+		val[i] = rb_float_new((*tp1)[i]);
 	return rb_ary_new4(12, val);
 }
 
@@ -1109,7 +1107,7 @@ s_Transform_Inspect(VALUE self)
 		rb_str_cat(val, "[", 1);
 		for (j = 0; j < 3; j++) {
 			double f;
-			f = (*tp)[s_index_order[i * 3 + j]];
+			f = (*tp)[i * 3 + j];
 			rb_funcall(val, mid, 1, rb_funcall(rb_float_new(f), mid2, 0));
 			if (j < 2)
 				rb_str_cat(val, ",", 1);
diff --git a/MolLib/Types.c b/MolLib/Types.c
index 9e2d920..523605d 100755
--- a/MolLib/Types.c
+++ b/MolLib/Types.c
@@ -90,13 +90,13 @@ MatrixRotation(Mat33 dst, const Vector *axis, Double angle)
 	double rcos = cos((double)angle);
 	double rsin = sin((double)angle);
 	dst[0] =            rcos + axis->x*axis->x*(1-rcos);
-	dst[1] =  axis->z * rsin + axis->y*axis->x*(1-rcos);
-	dst[2] = -axis->y * rsin + axis->z*axis->x*(1-rcos);
-	dst[3] = -axis->z * rsin + axis->x*axis->y*(1-rcos);
+	dst[1] = -axis->z * rsin + axis->x*axis->y*(1-rcos);
+	dst[2] =  axis->y * rsin + axis->x*axis->z*(1-rcos);
+	dst[3] =  axis->z * rsin + axis->y*axis->x*(1-rcos);
 	dst[4] =            rcos + axis->y*axis->y*(1-rcos);
-	dst[5] =  axis->x * rsin + axis->z*axis->y*(1-rcos);
-	dst[6] =  axis->y * rsin + axis->x*axis->z*(1-rcos);
-	dst[7] = -axis->x * rsin + axis->y*axis->z*(1-rcos);
+	dst[5] = -axis->x * rsin + axis->y*axis->z*(1-rcos);
+	dst[6] = -axis->y * rsin + axis->z*axis->x*(1-rcos);
+	dst[7] =  axis->x * rsin + axis->z*axis->y*(1-rcos);
 	dst[8] =            rcos + axis->z*axis->z*(1-rcos);
 }
 
@@ -104,9 +104,9 @@ void
 MatrixVec(Vector *dst, const Mat33 mat, const Vector *src)
 {
 	Vector temp;
-	temp.x = mat[0] * src->x + mat[1] * src->y + mat[2] * src->z;
-	temp.y = mat[3] * src->x + mat[4] * src->y + mat[5] * src->z;
-	temp.z = mat[6] * src->x + mat[7] * src->y + mat[8] * src->z;
+	temp.x = mat[0] * src->x + mat[3] * src->y + mat[6] * src->z;
+	temp.y = mat[1] * src->x + mat[4] * src->y + mat[7] * src->z;
+	temp.z = mat[2] * src->x + mat[5] * src->y + mat[8] * src->z;
 	*dst = temp;
 }
 
@@ -114,15 +114,15 @@ void
 MatrixMul(Mat33 dst, const Mat33 src1, const Mat33 src2)
 {
 	Mat33 temp;
-	temp[0] = src1[0] * src2[0] + src1[1] * src2[3] + src1[2] * src2[6];
-	temp[1] = src1[0] * src2[1] + src1[1] * src2[4] + src1[2] * src2[7];
-	temp[2] = src1[0] * src2[2] + src1[1] * src2[5] + src1[2] * src2[8];
-	temp[3] = src1[3] * src2[0] + src1[4] * src2[3] + src1[5] * src2[6];
-	temp[4] = src1[3] * src2[1] + src1[4] * src2[4] + src1[5] * src2[7];
-	temp[5] = src1[3] * src2[2] + src1[4] * src2[5] + src1[5] * src2[8];
-	temp[6] = src1[6] * src2[0] + src1[7] * src2[3] + src1[8] * src2[6];
-	temp[7] = src1[6] * src2[1] + src1[7] * src2[4] + src1[8] * src2[7];
-	temp[8] = src1[6] * src2[2] + src1[7] * src2[5] + src1[8] * src2[8];
+	temp[0] = src1[0] * src2[0] + src1[3] * src2[1] + src1[6] * src2[2];
+	temp[1] = src1[1] * src2[0] + src1[4] * src2[1] + src1[7] * src2[2];
+	temp[2] = src1[2] * src2[0] + src1[5] * src2[1] + src1[8] * src2[2];
+	temp[3] = src1[0] * src2[3] + src1[3] * src2[4] + src1[6] * src2[5];
+	temp[4] = src1[1] * src2[3] + src1[4] * src2[4] + src1[7] * src2[5];
+	temp[5] = src1[2] * src2[3] + src1[5] * src2[4] + src1[8] * src2[5];
+	temp[6] = src1[0] * src2[6] + src1[3] * src2[7] + src1[6] * src2[8];
+	temp[7] = src1[1] * src2[6] + src1[4] * src2[7] + src1[7] * src2[8];
+	temp[8] = src1[2] * src2[6] + src1[5] * src2[7] + src1[8] * src2[8];
 	memmove(dst, temp, sizeof(Mat33));
 }
 
@@ -184,9 +184,9 @@ MatrixScale(Mat33 dst, const Mat33 src, Double factor)
 void
 MatrixGeneralRotation(Mat33 dst, const Vector *v1, const Vector *v2, const Vector *v3)
 {
-	dst[0] = v1->x; dst[1] = v2->x; dst[2] = v3->x;
-	dst[3] = v1->y; dst[4] = v2->y; dst[5] = v3->y;
-	dst[6] = v1->z; dst[7] = v2->z; dst[8] = v3->z;
+	dst[0] = v1->x; dst[3] = v2->x; dst[6] = v3->x;
+	dst[1] = v1->y; dst[4] = v2->y; dst[7] = v3->y;
+	dst[2] = v1->z; dst[5] = v2->z; dst[8] = v3->z;
 }
 
 
@@ -202,18 +202,18 @@ int
 MatrixSymDiagonalize(Mat33 mat, Double *out_values, Vector *out_vectors)
 {
 	__CLPK_integer n, lda, lwork, info;
-	__CLPK_real a[9], w[3], work[9];
+	__CLPK_doublereal a[9], w[3], work[9];
 	int i, j;
 	for (i = 0; i < 3; i++) {
 		for (j = 0; j < 3; j++) {
-			a[j * 3 + i] = mat[i * 3 + j];
+			a[i * 3 + j] = mat[i * 3 + j];
 		}
 	}
 	n = lda = 3;
 	lwork = 9;
 	/*  For the meanings of the arguments, consult the LAPACK source; 
-		http://www.netlib.org/lapack/single/ssyev.f */
-	ssyev_("V", "U", &n, a, &lda, w, work, &lwork, &info);
+		http://www.netlib.org/lapack/double/dsyev.f */
+	dsyev_("V", "U", &n, a, &lda, w, work, &lwork, &info);
 	if (info == 0) {
 		for (i = 0; i < 3; i++) {
 			out_values[i] = w[i];
@@ -229,9 +229,9 @@ void
 TransformVec(Vector *dst, const Transform tf, const Vector *src)
 {
 	Vector temp;
-	temp.x = tf[0] * src->x + tf[1] * src->y + tf[2] * src->z + tf[9];
-	temp.y = tf[3] * src->x + tf[4] * src->y + tf[5] * src->z + tf[10];
-	temp.z = tf[6] * src->x + tf[7] * src->y + tf[8] * src->z + tf[11];
+	temp.x = tf[0] * src->x + tf[3] * src->y + tf[6] * src->z + tf[9];
+	temp.y = tf[1] * src->x + tf[4] * src->y + tf[7] * src->z + tf[10];
+	temp.z = tf[2] * src->x + tf[5] * src->y + tf[8] * src->z + tf[11];
 	*dst = temp;
 }
 
@@ -239,18 +239,18 @@ void
 TransformMul(Transform dst, const Transform src1, const Transform src2)
 {
 	Transform temp;
-	temp[0] = src1[0] * src2[0] + src1[1] * src2[3] + src1[2] * src2[6];
-	temp[1] = src1[0] * src2[1] + src1[1] * src2[4] + src1[2] * src2[7];
-	temp[2] = src1[0] * src2[2] + src1[1] * src2[5] + src1[2] * src2[8];
-	temp[3] = src1[3] * src2[0] + src1[4] * src2[3] + src1[5] * src2[6];
-	temp[4] = src1[3] * src2[1] + src1[4] * src2[4] + src1[5] * src2[7];
-	temp[5] = src1[3] * src2[2] + src1[4] * src2[5] + src1[5] * src2[8];
-	temp[6] = src1[6] * src2[0] + src1[7] * src2[3] + src1[8] * src2[6];
-	temp[7] = src1[6] * src2[1] + src1[7] * src2[4] + src1[8] * src2[7];
-	temp[8] = src1[6] * src2[2] + src1[7] * src2[5] + src1[8] * src2[8];
-	temp[9] = src1[0] * src2[9] + src1[1] * src2[10] + src1[2] * src2[11] + src1[9];
-	temp[10] = src1[3] * src2[9] + src1[4] * src2[10] + src1[5] * src2[11] + src1[10];
-	temp[11] = src1[6] * src2[9] + src1[7] * src2[10] + src1[8] * src2[11] + src1[11];
+	temp[0] = src1[0] * src2[0] + src1[3] * src2[1] + src1[6] * src2[2];
+	temp[1] = src1[1] * src2[0] + src1[4] * src2[1] + src1[7] * src2[2];
+	temp[2] = src1[2] * src2[0] + src1[5] * src2[1] + src1[8] * src2[2];
+	temp[3] = src1[0] * src2[3] + src1[3] * src2[4] + src1[6] * src2[5];
+	temp[4] = src1[1] * src2[3] + src1[4] * src2[4] + src1[7] * src2[5];
+	temp[5] = src1[2] * src2[3] + src1[5] * src2[4] + src1[8] * src2[5];
+	temp[6] = src1[0] * src2[6] + src1[3] * src2[7] + src1[6] * src2[8];
+	temp[7] = src1[1] * src2[6] + src1[4] * src2[7] + src1[7] * src2[8];
+	temp[8] = src1[2] * src2[6] + src1[5] * src2[7] + src1[8] * src2[8];
+	temp[9] = src1[0] * src2[9] + src1[3] * src2[10] + src1[6] * src2[11] + src1[9];
+	temp[10] = src1[1] * src2[9] + src1[4] * src2[10] + src1[7] * src2[11] + src1[10];
+	temp[11] = src1[2] * src2[9] + src1[5] * src2[10] + src1[8] * src2[11] + src1[11];
 	memmove(dst, temp, sizeof(Transform));
 }
 
@@ -281,9 +281,9 @@ TransformInvert(Transform dst, const Transform src)
 	Transform temp;
 	int n = MatrixInvert(temp, src);
 	if (n == 0) {
-		temp[9] = -temp[0] * src[9] - temp[1] * src[10] - temp[2] * src[11];
-		temp[10] = -temp[3] * src[9] - temp[4] * src[10] - temp[5] * src[11];
-		temp[11] = -temp[6] * src[9] - temp[7] * src[10] - temp[8] * src[11];
+		temp[9] = -temp[0] * src[9] - temp[3] * src[10] - temp[6] * src[11];
+		temp[10] = -temp[1] * src[9] - temp[4] * src[10] - temp[7] * src[11];
+		temp[11] = -temp[2] * src[9] - temp[5] * src[10] - temp[8] * src[11];
 		memmove(dst, temp, sizeof(Transform));
 		return 0;
 	} else return n;
@@ -432,6 +432,7 @@ DeleteArray(void *base, Int *count, int item_size, int idx, int nitems, void *ou
 		free(*bp);
 		*bp = NULL;
 	}
+	return NULL;
 }
 
 int
diff --git a/MolLib/Types.h b/MolLib/Types.h
index 78bfa8e..dd9ac80 100755
--- a/MolLib/Types.h
+++ b/MolLib/Types.h
@@ -64,13 +64,13 @@ typedef unsigned int UInt;
 
 typedef struct Vector { Double x, y, z; } Vector;
 typedef struct Quat { Double x, y, z, w; } Quat;   /*  A quaternion  */
-typedef Double Mat33[9];
+typedef Double Mat33[9];  /*  Columns first!  */
 
-/*  Mat33 (rot) and Vector (translation)  */
-/*  Transform[0..8] are row-first, and Transform[9..11] are the fourth column. So that, the
+/*  Mat33 (rot) and Vector (translation) */
+/*  Transform[0..8] are the same as Mat33, and Transform[9..11] are the fourth column. So that, the
     matrix elements are in the following order;
-    {a11, a12, a13, a21, a22, a23, a31, a32, a33, a14, a24, a34}
-   Sorry for confusion!  */
+    {a11, a21, a31, a12, a22, a32, a13, a23, a33, a14, a24, a34}  */
+/*  (Changed from row-first Mat33 on 2011.12.7.)  */
 typedef Double Transform[12];
 typedef Double *TransformPtr;
 /*typedef Double Matrix[4][4]; */