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++) {
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];
}
}
}
*/
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;
}
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;
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);
}
}
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;
}
* 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 */
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);
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);
}
-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
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];
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);
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);
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;
}
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) */
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];
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 */
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");
}
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;
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 */
} 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);
}
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 */
#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)
{
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;
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;
}
/* 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))));
}
}
}
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);
}
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);
}
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);
}
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);
}
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);
}
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);
}
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);
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);
}
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;
}
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));
}
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;
}
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];
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;
}
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));
}
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;
free(*bp);
*bp = NULL;
}
+ return NULL;
}
int
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]; */
OSDN Git Service
