*/
#define CALC_SPLINE(order) \
{ \
- int j, k, l; \
- real dr, div; \
- real data[PME_ORDER_MAX]; \
- real ddata[PME_ORDER_MAX]; \
- \
- for (j = 0; (j < DIM); j++) \
+ for (int j = 0; (j < DIM); j++) \
{ \
+ real dr, div; \
+ real data[PME_ORDER_MAX]; \
+ \
dr = xptr[j]; \
\
/* dr is relative offset from lower cell limit */ \
data[1] = dr; \
data[0] = 1 - dr; \
\
- for (k = 3; (k < order); k++) \
+ for (int k = 3; (k < order); k++) \
{ \
div = 1.0/(k - 1.0); \
data[k-1] = div*dr*data[k-2]; \
- for (l = 1; (l < (k-1)); l++) \
+ for (int l = 1; (l < (k-1)); l++) \
{ \
data[k-l-1] = div*((dr+l)*data[k-l-2]+(k-l-dr)* \
data[k-l-1]); \
data[0] = div*(1-dr)*data[0]; \
} \
/* differentiate */ \
- ddata[0] = -data[0]; \
- for (k = 1; (k < order); k++) \
+ dtheta[j][i*order+0] = -data[0]; \
+ for (int k = 1; (k < order); k++) \
{ \
- ddata[k] = data[k-1] - data[k]; \
+ dtheta[j][i*order+k] = data[k-1] - data[k]; \
} \
\
div = 1.0/(order - 1); \
data[order-1] = div*dr*data[order-2]; \
- for (l = 1; (l < (order-1)); l++) \
+ for (int l = 1; (l < (order-1)); l++) \
{ \
data[order-l-1] = div*((dr+l)*data[order-l-2]+ \
(order-l-dr)*data[order-l-1]); \
} \
data[0] = div*(1 - dr)*data[0]; \
\
- for (k = 0; k < order; k++) \
+ for (int k = 0; k < order; k++) \
{ \
theta[j][i*order+k] = data[k]; \
- dtheta[j][i*order+k] = ddata[k]; \
} \
} \
}
if (bDoSplines || coefficient[ii] != 0.0)
{
xptr = fractx[ii];
+ assert(order >= 4 && order <= PME_ORDER_MAX);
switch (order)
{
case 4: CALC_SPLINE(4); break;
for (i = 0; i < il->nr; i += 1+NRAL(ft))
{
const t_iparams *ip;
- real cam[5] = {0}, inv_mass, coeff, m_aj;
- int a1, j, aj;
+ real inv_mass, coeff, m_aj;
+ int a1, aj;
ip = &ffparams->iparams[il->iatoms[i]];
if (ft != F_VSITEN)
{
- for (j = 1; j < NRAL(ft); j++)
+ /* Only vsiten can have more than four
+ constructing atoms, so NRAL(ft) <= 5 */
+ int j;
+ real *cam;
+ const int maxj = NRAL(ft);
+
+ snew(cam, maxj);
+ assert(maxj <= 5);
+ for (j = 1; j < maxj; j++)
{
cam[j] = moltype->atoms.atom[il->iatoms[i+1+j]].m;
if (cam[j] == 0)
il->iatoms[i+1+j]+1);
}
}
- }
- switch (ft)
- {
- case F_VSITE2:
- /* Exact */
- vsite_m[a1] = (cam[1]*cam[2])/(cam[2]*sqr(1-ip->vsite.a) + cam[1]*sqr(ip->vsite.a));
- break;
- case F_VSITE3:
- /* Exact */
- vsite_m[a1] = (cam[1]*cam[2]*cam[3])/(cam[2]*cam[3]*sqr(1-ip->vsite.a-ip->vsite.b) + cam[1]*cam[3]*sqr(ip->vsite.a) + cam[1]*cam[2]*sqr(ip->vsite.b));
- break;
- case F_VSITEN:
- /* Exact */
- inv_mass = 0;
- for (j = 0; j < 3*ffparams->iparams[il->iatoms[i]].vsiten.n; j += 3)
- {
- aj = il->iatoms[i+j+2];
- coeff = ffparams->iparams[il->iatoms[i+j]].vsiten.a;
- if (moltype->atoms.atom[aj].ptype == eptVSite)
- {
- m_aj = vsite_m[aj];
- }
- else
- {
- m_aj = moltype->atoms.atom[aj].m;
- }
- if (m_aj <= 0)
+ switch (ft)
+ {
+ case F_VSITE2:
+ /* Exact */
+ vsite_m[a1] = (cam[1]*cam[2])/(cam[2]*sqr(1-ip->vsite.a) + cam[1]*sqr(ip->vsite.a));
+ break;
+ case F_VSITE3:
+ /* Exact */
+ vsite_m[a1] = (cam[1]*cam[2]*cam[3])/(cam[2]*cam[3]*sqr(1-ip->vsite.a-ip->vsite.b) + cam[1]*cam[3]*sqr(ip->vsite.a) + cam[1]*cam[2]*sqr(ip->vsite.b));
+ break;
+ case F_VSITEN:
+ gmx_incons("Invalid vsite type");
+ break;
+ default:
+ /* Use the mass of the lightest constructing atom.
+ * This is an approximation.
+ * If the distance of the virtual site to the
+ * constructing atom is less than all distances
+ * between constructing atoms, this is a safe
+ * over-estimate of the displacement of the vsite.
+ * This condition holds for all H mass replacement
+ * vsite constructions, except for SP2/3 groups.
+ * In SP3 groups one H will have a F_VSITE3
+ * construction, so even there the total drift
+ * estimate shouldn't be far off.
+ */
+ vsite_m[a1] = cam[1];
+ for (j = 2; j < maxj; j++)
{
- gmx_incons("The mass of a vsiten constructing atom is <= 0");
+ vsite_m[a1] = min(vsite_m[a1], cam[j]);
}
- inv_mass += coeff*coeff/m_aj;
+ (*n_nonlin_vsite)++;
+ break;
+ }
+ sfree(cam);
+ }
+ else
+ {
+ int j;
+
+ /* Exact */
+ inv_mass = 0;
+ for (j = 0; j < 3*ffparams->iparams[il->iatoms[i]].vsiten.n; j += 3)
+ {
+ aj = il->iatoms[i+j+2];
+ coeff = ffparams->iparams[il->iatoms[i+j]].vsiten.a;
+ if (moltype->atoms.atom[aj].ptype == eptVSite)
+ {
+ m_aj = vsite_m[aj];
+ }
+ else
+ {
+ m_aj = moltype->atoms.atom[aj].m;
}
- vsite_m[a1] = 1/inv_mass;
- /* Correct for loop increment of i */
- i += j - 1 - NRAL(ft);
- break;
- default:
- /* Use the mass of the lightest constructing atom.
- * This is an approximation.
- * If the distance of the virtual site to the
- * constructing atom is less than all distances
- * between constructing atoms, this is a safe
- * over-estimate of the displacement of the vsite.
- * This condition holds for all H mass replacement
- * vsite constructions, except for SP2/3 groups.
- * In SP3 groups one H will have a F_VSITE3
- * construction, so even there the total drift
- * estimate shouldn't be far off.
- */
- assert(j >= 1);
- vsite_m[a1] = cam[1];
- for (j = 2; j < NRAL(ft); j++)
+ if (m_aj <= 0)
{
- vsite_m[a1] = min(vsite_m[a1], cam[j]);
+ gmx_incons("The mass of a vsiten constructing atom is <= 0");
}
- (*n_nonlin_vsite)++;
- break;
+ inv_mass += coeff*coeff/m_aj;
+ }
+ vsite_m[a1] = 1/inv_mass;
+ /* Correct for loop increment of i */
+ i += j - 1 - NRAL(ft);
}
if (gmx_debug_at)
{
#include "pull.h"
+#include <assert.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "gromacs/pbcutil/pbc.h"
#include "gromacs/topology/mtop_util.h"
#include "gromacs/utility/futil.h"
+#include "gromacs/utility/gmxassert.h"
#include "gromacs/utility/smalloc.h"
static void pull_print_group_x(FILE *out, ivec dim, const t_pull_group *pgrp)
/* Determine if we need to take PBC into account for calculating
* the COM's of the pull groups.
*/
+ GMX_RELEASE_ASSERT(pull->npbcdim <= DIM, "npbcdim must be <= the number of dimensions");
for (m = 0; m < pull->npbcdim; m++)
{
+ GMX_RELEASE_ASSERT(m <= DIM, "npbcdim must be <= the number of dimensions");
if (pulldim[m] == 1 && pgrp->nat > 1)
{
if (pgrp->pbcatom >= 0)
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff