static real ekrot(rvec x[], rvec v[], real mass[], int isize, int index[])
{
- static real **TCM = nullptr, **L;
+ real TCM[5][5], L[5][5];
double tm, m0, lxx, lxy, lxz, lyy, lyz, lzz, ekrot;
rvec a0, ocm;
dvec dx, b0;
dvec xcm, vcm, acm;
int i, j, m, n;
- if (TCM == nullptr)
- {
- snew(TCM, DIM);
- for (i = 0; i < DIM; i++)
- {
- snew(TCM[i], DIM);
- }
- snew(L, DIM);
- for (i = 0; i < DIM; i++)
- {
- snew(L[i], DIM);
- }
- }
-
clear_dvec(xcm);
clear_dvec(vcm);
clear_dvec(acm);
L[YY][ZZ] = -lyz;
L[ZZ][ZZ] = lxx + lyy;
- m_inv_gen(L, DIM, TCM);
+ m_inv_gen(&L[0][0], DIM, &TCM[0][0]);
/* Compute omega (hoeksnelheid) */
clear_rvec(ocm);
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2012,2014,2015, by the GROMACS development team, led by
+ * Copyright (c) 2012,2014,2015,2017, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
gmx_fatal(FARGS, "Error: Too many iterations in routine JACOBI\n");
}
-int m_inv_gen(real **m, int n, real **minv)
+int m_inv_gen(real *m, int n, real *minv)
{
double **md, **v, *eig, tol, s;
int nzero, i, j, k, nrot;
{
for (j = 0; j < n; j++)
{
- md[i][j] = m[i][j];
+ md[i][j] = m[i*n + j];
}
}
{
s += eig[k]*v[i][k]*v[j][k];
}
- minv[i][j] = s;
+ minv[i*n + j] = s;
}
}
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2010,2014, by the GROMACS development team, led by
+ * Copyright (c) 2010,2014,2017, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
* int *irot = number of jacobi rotations
*/
-int m_inv_gen(real **m, int n, real **minv);
-/* Produces minv, a generalized inverse of m.
+int m_inv_gen(real *m, int n, real *minv);
+/* Produces minv, a generalized inverse of m, both stored as linear arrays.
* Inversion is done via diagonalization,
* eigenvalues smaller than 1e-6 times the average diagonal element
* are assumed to be zero.
{
snew(od->otav, od->nr);
}
- snew(od->tmp, od->nex);
- snew(od->TMP, od->nex);
- for (int ex = 0; ex < od->nex; ex++)
- {
- snew(od->TMP[ex], 5);
- for (int i = 0; i < 5; i++)
- {
- snew(od->TMP[ex][i], 5);
- }
- }
+ snew(od->tmpEq, od->nex);
od->nref = 0;
for (int i = 0; i < mtop->natoms; i++)
{
int nref;
real edt, edt_1, invn, pfac, r2, invr, corrfac, wsv2, sw, dev;
- tensor *S, R, TMP;
- rvec5 *rhs;
- real *mref, ***T;
+ OriresMatEq *matEq;
+ real *mref;
double mtot;
rvec *xref, *xtmp, com, r_unrot, r;
t_oriresdata *od;
bTAV = (od->edt != 0);
edt = od->edt;
edt_1 = od->edt_1;
- S = od->S;
- T = od->TMP;
- rhs = od->tmp;
+ matEq = od->tmpEq;
nref = od->nref;
mref = od->mref;
xref = od->xref;
rvec_dec(xtmp[j], com);
}
/* Calculate the rotation matrix to rotate x to the reference orientation */
- calc_fit_R(DIM, nref, mref, xref, xtmp, R);
- copy_mat(R, od->R);
+ calc_fit_R(DIM, nref, mref, xref, xtmp, od->R);
/* Index restraint data in order of appearance in forceatoms */
int restraintIndex = 0;
{
rvec_sub(x[forceatoms[fa+1]], x[forceatoms[fa+2]], r_unrot);
}
- mvmul(R, r_unrot, r);
+ mvmul(od->R, r_unrot, r);
r2 = norm2(r);
invr = gmx::invsqrt(r2);
/* Calculate the prefactor for the D tensor, this includes the factor 3! */
{
for (int i = 0; i < 5; i++)
{
- rhs[ex][i] = 0;
+ matEq[ex].rhs[i] = 0;
for (int j = 0; j <= i; j++)
{
- T[ex][i][j] = 0;
+ matEq[ex].mat[i][j] = 0;
}
}
}
/* Calculate the vector rhs and half the matrix T for the 5 equations */
for (int i = 0; i < 5; i++)
{
- rhs[ex][i] += Dtav[i]*ip[type].orires.obs*weight;
+ matEq[ex].rhs[i] += Dtav[i]*ip[type].orires.obs*weight;
for (int j = 0; j <= i; j++)
{
- T[ex][i][j] += Dtav[i]*Dtav[j]*weight;
+ matEq[ex].mat[i][j] += Dtav[i]*Dtav[j]*weight;
}
}
restraintIndex++;
/* Now we have all the data we can calculate S */
for (int ex = 0; ex < od->nex; ex++)
{
+ OriresMatEq &eq = matEq[ex];
/* Correct corrfac and copy one half of T to the other half */
for (int i = 0; i < 5; i++)
{
- rhs[ex][i] *= corrfac;
- T[ex][i][i] *= gmx::square(corrfac);
+ eq.rhs[i] *= corrfac;
+ eq.mat[i][i] *= gmx::square(corrfac);
for (int j = 0; j < i; j++)
{
- T[ex][i][j] *= gmx::square(corrfac);
- T[ex][j][i] = T[ex][i][j];
+ eq.mat[i][j] *= gmx::square(corrfac);
+ eq.mat[j][i] = eq.mat[i][j];
}
}
- m_inv_gen(T[ex], 5, T[ex]);
+ m_inv_gen(&eq.mat[0][0], 5, &eq.mat[0][0]);
/* Calculate the orientation tensor S for this experiment */
- S[ex][0][0] = 0;
- S[ex][0][1] = 0;
- S[ex][0][2] = 0;
- S[ex][1][1] = 0;
- S[ex][1][2] = 0;
+ matrix &S = od->S[ex];
+ S[0][0] = 0;
+ S[0][1] = 0;
+ S[0][2] = 0;
+ S[1][1] = 0;
+ S[1][2] = 0;
for (int i = 0; i < 5; i++)
{
- S[ex][0][0] += 1.5*T[ex][0][i]*rhs[ex][i];
- S[ex][0][1] += 1.5*T[ex][1][i]*rhs[ex][i];
- S[ex][0][2] += 1.5*T[ex][2][i]*rhs[ex][i];
- S[ex][1][1] += 1.5*T[ex][3][i]*rhs[ex][i];
- S[ex][1][2] += 1.5*T[ex][4][i]*rhs[ex][i];
+ S[0][0] += 1.5*eq.mat[0][i]*eq.rhs[i];
+ S[0][1] += 1.5*eq.mat[1][i]*eq.rhs[i];
+ S[0][2] += 1.5*eq.mat[2][i]*eq.rhs[i];
+ S[1][1] += 1.5*eq.mat[3][i]*eq.rhs[i];
+ S[1][2] += 1.5*eq.mat[4][i]*eq.rhs[i];
}
- S[ex][1][0] = S[ex][0][1];
- S[ex][2][0] = S[ex][0][2];
- S[ex][2][1] = S[ex][1][2];
- S[ex][2][2] = -S[ex][0][0] - S[ex][1][1];
+ S[1][0] = S[0][1];
+ S[2][0] = S[0][2];
+ S[2][1] = S[1][2];
+ S[2][2] = -S[0][0] - S[1][1];
}
- wsv2 = 0;
- sw = 0;
+ const matrix *S = od->S;
+
+ wsv2 = 0;
+ sw = 0;
/* Index restraint data in order of appearance in forceatoms */
restraintIndex = 0;
/* Rotate the S matrices back, so we get the correct grad(tr(S D)) */
for (int ex = 0; ex < od->nex; ex++)
{
- tmmul(R, S[ex], TMP);
- mmul(TMP, R, S[ex]);
+ matrix RS;
+ tmmul(od->R, od->S[ex], RS);
+ mmul(RS, od->R, od->S[ex]);
}
return od->rmsdev;
*
* Copyright (c) 1991-2000, University of Groningen, The Netherlands.
* Copyright (c) 2001-2004, The GROMACS development team.
- * Copyright (c) 2012,2014,2015,2016, by the GROMACS development team, led by
+ * Copyright (c) 2012,2014,2015,2016,2017, by the GROMACS development team, led by
* Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
* and including many others, as listed in the AUTHORS file in the
* top-level source directory and at http://www.gromacs.org.
const t_iatom *forceatomsStart; /* Pointer to the start of the disre forceatoms */
} t_disresdata;
+/* All coefficients for the matrix equation for the orientation tensor */
+struct OriresMatEq
+{
+ real rhs[5]; /* The right hand side of the matrix equation */
+ real mat[5][5]; /* The matrix */
+};
/* Orientation restraining stuff */
typedef struct t_oriresdata {
- real fc; /* Force constant for the restraints */
- real edt; /* Multiplication factor for time averaging */
- real edt_1; /* 1 - edt */
- real exp_min_t_tau; /* Factor for slowly switching on the force */
- int nr; /* The number of orientation restraints */
- int nex; /* The number of experiments */
- int nref; /* The number of atoms for the fit */
- real *mref; /* The masses of the reference atoms */
- rvec *xref; /* The reference coordinates for the fit (nref) */
- rvec *xtmp; /* Temporary array for fitting (nref) */
- matrix R; /* Rotation matrix to rotate to the reference coor. */
- tensor *S; /* Array of order tensors for each experiment (nexp) */
- rvec5 *Dinsl; /* The order matrix D for all restraints (nr x 5) */
- rvec5 *Dins; /* The ensemble averaged D (nr x 5) */
- rvec5 *Dtav; /* The time and ensemble averaged D (nr x 5) */
- real *oinsl; /* The calculated instantaneous orientations */
- real *oins; /* The calculated emsemble averaged orientations */
- real *otav; /* The calculated time and ensemble averaged orient. */
- real rmsdev; /* The weighted (using kfac) RMS deviation */
- rvec5 *tmp; /* An array of temporary 5-vectors (nex); */
- real ***TMP; /* An array of temporary 5x5 matrices (nex); */
- real *eig; /* Eigenvalues/vectors, for output only (nex x 12) */
+ real fc; /* Force constant for the restraints */
+ real edt; /* Multiplication factor for time averaging */
+ real edt_1; /* 1 - edt */
+ real exp_min_t_tau; /* Factor for slowly switching on the force */
+ int nr; /* The number of orientation restraints */
+ int nex; /* The number of experiments */
+ int nref; /* The number of atoms for the fit */
+ real *mref; /* The masses of the reference atoms */
+ rvec *xref; /* The reference coordinates for the fit (nref) */
+ rvec *xtmp; /* Temporary array for fitting (nref) */
+ matrix R; /* Rotation matrix to rotate to the reference coor. */
+ tensor *S; /* Array of order tensors for each experiment (nexp) */
+ rvec5 *Dinsl; /* The order matrix D for all restraints (nr x 5) */
+ rvec5 *Dins; /* The ensemble averaged D (nr x 5) */
+ rvec5 *Dtav; /* The time and ensemble averaged D (nr x 5) */
+ real *oinsl; /* The calculated instantaneous orientations */
+ real *oins; /* The calculated emsemble averaged orientations */
+ real *otav; /* The calculated time and ensemble averaged orient. */
+ real rmsdev; /* The weighted (using kfac) RMS deviation */
+ OriresMatEq *tmpEq; /* An temporary array of matrix + rhs */
+ real *eig; /* Eigenvalues/vectors, for output only (nex x 12) */
/* variables for diagonalization with diagonalize_orires_tensors()*/
double **M;
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff