/*
+ * This file is part of the GROMACS molecular simulation package.
*
- * This source code is part of
+ * Copyright (c) 2009,2010,2011,2012,2013,2014, 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.
*
- * G R O M A C S
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
*
- * GROningen MAchine for Chemical Simulations
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
*
- * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
- * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2009, The GROMACS development team,
- * check out http://www.gromacs.org for more information.
-
- * This program is free software; you can redistribute it and/or
- * modify it under the terms of the GNU General Public License
- * as published by the Free Software Foundation; either version 2
- * of the License, or (at your option) any later version.
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
*
- * If you want to redistribute modifications, please consider that
- * scientific software is very special. Version control is crucial -
- * bugs must be traceable. We will be happy to consider code for
- * inclusion in the official distribution, but derived work must not
- * be called official GROMACS. Details are found in the README & COPYING
- * files - if they are missing, get the official version at www.gromacs.org.
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
*
* To help us fund GROMACS development, we humbly ask that you cite
- * the papers on the package - you can find them in the top README file.
- *
- * For more info, check our website at http://www.gromacs.org
+ * the research papers on the package. Check out http://www.gromacs.org.
*/
/*! \internal \file
* \brief
* For best results, the setup should probably be postponed (at least
* partially) to gmx_ana_poscalc_init_eval().
*
- * \author Teemu Murtola <teemu.murtola@cbr.su.se>
+ * \author Teemu Murtola <teemu.murtola@gmail.com>
* \ingroup module_selection
*/
+#include "gmxpre.h"
+
+#include "poscalc.h"
+
#include <string.h>
-#include "gromacs/legacyheaders/smalloc.h"
-#include "gromacs/legacyheaders/typedefs.h"
-#include "gromacs/legacyheaders/pbc.h"
-#include "gromacs/legacyheaders/vec.h"
+#include <algorithm>
-#include "gromacs/selection/centerofmass.h"
+#include "gromacs/fileio/trx.h"
+#include "gromacs/math/vec.h"
#include "gromacs/selection/indexutil.h"
-#include "gromacs/selection/poscalc.h"
#include "gromacs/selection/position.h"
#include "gromacs/utility/exceptions.h"
#include "gromacs/utility/gmxassert.h"
+#include "gromacs/utility/smalloc.h"
+
+#include "centerofmass.h"
namespace gmx
{
* already been calculated in \p sbase.
* The structure pointed by \p sbase is always a static calculation.
*/
- struct gmx_ana_poscalc_t *sbase;
+ gmx_ana_poscalc_t *sbase;
/** Next structure in the linked list of calculations. */
- struct gmx_ana_poscalc_t *next;
+ gmx_ana_poscalc_t *next;
/** Previous structure in the linked list of calculations. */
- struct gmx_ana_poscalc_t *prev;
+ gmx_ana_poscalc_t *prev;
/** Number of references to this structure. */
- int refcount;
+ int refcount;
/** Collection this calculation belongs to. */
gmx::PositionCalculationCollection::Impl *coll;
};
static e_index_t
index_type_for_poscalc(e_poscalc_t type)
{
- switch(type)
+ switch (type)
{
case POS_ATOM: return INDEX_ATOM;
case POS_RES: return INDEX_RES;
{
*flags &= ~POS_COMPLMAX;
*flags |= POS_COMPLWHOLE;
- ptr = post + 6;
+ ptr = post + 6;
}
else if (post[0] == 'p')
{
*flags &= ~POS_COMPLWHOLE;
*flags |= POS_COMPLMAX;
- ptr = post + 5;
+ ptr = post + 5;
}
else if (post[0] == 'd')
{
*flags &= ~(POS_COMPLMAX | POS_COMPLWHOLE);
- ptr = post + 4;
+ ptr = post + 4;
}
if (ptr[0] == 'r')
gmx_ana_poscalc_t *base;
fprintf(fp, " Base: ");
- j = 1;
+ j = 1;
base = impl_->first_;
while (base && base != pc->sbase)
{
PositionCalculationCollection::createCalculationFromEnum(const char *post, int flags)
{
e_poscalc_t type;
- int cflags = flags;
+ int cflags = flags;
typeFromEnum(post, &type, &cflags);
return impl_->createCalculation(type, cflags);
}
+int PositionCalculationCollection::getHighestRequiredAtomIndex() const
+{
+ int result = 0;
+ gmx_ana_poscalc_t *pc = impl_->first_;
+ while (pc)
+ {
+ // Calculations with a base just copy positions from the base, so
+ // those do not need to be considered in the check.
+ if (!pc->sbase)
+ {
+ gmx_ana_index_t g;
+ gmx_ana_index_set(&g, pc->b.nra, pc->b.a, 0);
+ result = std::max(result, gmx_ana_index_get_max_index(&g));
+ }
+ pc = pc->next;
+ }
+ return result;
+}
+
void PositionCalculationCollection::initEvaluation()
{
if (impl_->bInit_)
while (pc)
{
pc->bEval = false;
- pc = pc->next;
+ pc = pc->next;
}
}
/* Set the type to POS_ATOM if the calculation in fact is such. */
if (pc->b.nr == pc->b.nra)
{
- pc->type = POS_ATOM;
+ pc->type = POS_ATOM;
pc->flags &= ~(POS_MASS | POS_COMPLMAX | POS_COMPLWHOLE);
}
/* Set the POS_COMPLWHOLE flag if the calculation in fact always uses
if ((pc->flags & POS_COMPLWHOLE) && !bBase && pc->b.nra > g->isize)
{
gmx_ana_index_copy(&pc->gmax, g, true);
- sfree(pc->gmax.name);
- pc->gmax.name = NULL;
}
else
{
- gmx_ana_index_set(&pc->gmax, pc->b.nra, pc->b.a, NULL, 0);
+ gmx_ana_index_set(&pc->gmax, pc->b.nra, pc->b.a, 0);
}
}
return false;
}
/* Find the overlap between the calculations. */
- gmx_ana_index_set(&g2, pc2->b.nra, pc2->b.a, NULL, 0);
+ gmx_ana_index_set(&g2, pc2->b.nra, pc2->b.a, 0);
gmx_ana_index_intersection(g, g1, &g2);
/* Do not merge if there is no overlap. */
if (g->isize == 0)
int flags;
flags = pc->flags & ~(POS_DYNAMIC | POS_MASKONLY);
- base = pc->coll->createCalculation(pc->type, flags);
+ base = pc->coll->createCalculation(pc->type, flags);
set_poscalc_maxindex(base, &pc->gmax, true);
- snew(base->p, 1);
+ base->p = new gmx_ana_pos_t();
pc->sbase = base;
pc->coll->removeCalculation(base);
int i, bi, bj, bo;
base->flags |= pc->flags & (POS_VELOCITIES | POS_FORCES);
- gmx_ana_index_set(&gp, pc->b.nra, pc->b.a, NULL, 0);
- gmx_ana_index_set(&gb, base->b.nra, base->b.a, NULL, 0);
+ gmx_ana_index_set(&gp, pc->b.nra, pc->b.a, 0);
+ gmx_ana_index_set(&gb, base->b.nra, base->b.a, 0);
isize = gmx_ana_index_difference_size(&gp, &gb);
if (isize > 0)
{
gmx_ana_index_merge(&g, &gb, &g);
/* Merge the blocks */
srenew(base->b.index, base->b.nr + bnr + 1);
- i = g.isize - 1;
- bi = base->b.nr - 1;
- bj = pc->b.nr - 1;
- bo = base->b.nr + bnr - 1;
+ i = g.isize - 1;
+ bi = base->b.nr - 1;
+ bj = pc->b.nr - 1;
+ bo = base->b.nr + bnr - 1;
base->b.index[bo+1] = i + 1;
while (bo >= 0)
{
base->b.a = g.index;
base->b.nalloc_a = g.isize;
/* Refresh the gmax field */
- gmx_ana_index_set(&base->gmax, base->b.nra, base->b.a, NULL, 0);
+ gmx_ana_index_set(&base->gmax, base->b.nra, base->b.a, 0);
}
}
return;
}
- gmx_ana_index_set(&gp, pc->b.nra, pc->b.a, NULL, 0);
+ gmx_ana_index_set(&gp, pc->b.nra, pc->b.a, 0);
gmx_ana_index_clear(&g);
gmx_ana_index_reserve(&g, pc->b.nra);
pbase = pc;
- base = pc->coll->first_;
+ base = pc->coll->first_;
while (base)
{
/* Save the next calculation so that we can safely delete base */
merge_bases(pbase, base);
}
}
- gmx_ana_index_set(&gp, pbase->b.nra, pbase->b.a, NULL, 0);
+ gmx_ana_index_set(&gp, pbase->b.nra, pbase->b.a, 0);
gmx_ana_index_reserve(&g, pc->b.nra);
}
/* Proceed to the next unchecked calculation */
{
gmx_ana_indexmap_set_static(&p->m, &pc->b);
}
- gmx_ana_pos_reserve(p, p->m.nr, 0);
+ gmx_ana_pos_reserve(p, p->m.mapb.nr, 0);
if (pc->flags & POS_VELOCITIES)
{
gmx_ana_pos_reserve_velocities(p);
{
gmx_ana_pos_reserve_forces(p);
}
- gmx_ana_pos_set_nr(p, p->m.nr);
- gmx_ana_pos_set_evalgrp(p, &pc->gmax);
}
/*!
{
gmx_ana_index_deinit(&pc->gmax);
}
- if (pc->p)
- {
- gmx_ana_pos_free(pc->p);
- }
+ delete pc->p;
if (pc->sbase)
{
gmx_ana_poscalc_free(pc->sbase);
bool
gmx_ana_poscalc_requires_top(gmx_ana_poscalc_t *pc)
{
- if ((pc->flags & POS_MASS) || pc->type == POS_RES || pc->type == POS_MOL)
+ if ((pc->flags & POS_MASS) || pc->type == POS_RES || pc->type == POS_MOL
+ || ((pc->flags & POS_FORCES) && pc->type != POS_ATOM))
{
return true;
}
gmx_ana_index_t *g, t_trxframe *fr, t_pbc *pbc)
{
int i, bi, bj;
-
+
if (pc->bEval == true && !(pc->flags & POS_MASKONLY))
{
return;
{
g = &pc->gmax;
}
- gmx_ana_pos_set_evalgrp(p, g);
/* Update the index map */
if (pc->flags & POS_DYNAMIC)
{
gmx_ana_indexmap_update(&p->m, g, false);
- p->nr = p->m.nr;
}
else if (pc->flags & POS_MASKONLY)
{
gmx_ana_indexmap_update(&p->m, g, true);
if (pc->bEval)
+ {
return;
+ }
}
if (!(pc->flags & POS_DYNAMIC))
{
* loop instead of in the beginning. */
if (pc->flags & POS_DYNAMIC)
{
- for (bi = 0; bi < p->nr; ++bi)
+ for (bi = 0; bi < p->count(); ++bi)
{
bj = pc->baseid[p->m.refid[bi]];
copy_rvec(pc->sbase->p->x[bj], p->x[bi]);
}
if (p->v)
{
- for (bi = 0; bi < p->nr; ++bi)
+ for (bi = 0; bi < p->count(); ++bi)
{
bj = pc->baseid[p->m.refid[bi]];
copy_rvec(pc->sbase->p->v[bj], p->v[bi]);
}
if (p->f)
{
- for (bi = 0; bi < p->nr; ++bi)
+ for (bi = 0; bi < p->count(); ++bi)
{
bj = pc->baseid[p->m.refid[bi]];
copy_rvec(pc->sbase->p->f[bj], p->f[bi]);
}
else
{
- for (bi = 0; bi < p->nr; ++bi)
+ for (bi = 0; bi < p->count(); ++bi)
{
bj = pc->baseid[bi];
copy_rvec(pc->sbase->p->x[bj], p->x[bi]);
}
if (p->v)
{
- for (bi = 0; bi < p->nr; ++bi)
+ for (bi = 0; bi < p->count(); ++bi)
{
bj = pc->baseid[bi];
copy_rvec(pc->sbase->p->v[bj], p->v[bi]);
}
if (p->f)
{
- for (bi = 0; bi < p->nr; ++bi)
+ for (bi = 0; bi < p->count(); ++bi)
{
bj = pc->baseid[bi];
copy_rvec(pc->sbase->p->f[bj], p->f[bi]);
}
/* Here, we assume that the topology has been properly initialized,
* and do not check the return values of gmx_calc_comg*(). */
- t_topology *top = pc->coll->top_;
- bool bMass = pc->flags & POS_MASS;
+ t_topology *top = pc->coll->top_;
+ bool bMass = pc->flags & POS_MASS;
switch (pc->type)
{
- case POS_ATOM:
- for (i = 0; i < pc->b.nra; ++i)
- {
- copy_rvec(fr->x[pc->b.a[i]], p->x[i]);
- }
- if (p->v && fr->bV)
- {
+ case POS_ATOM:
for (i = 0; i < pc->b.nra; ++i)
{
- copy_rvec(fr->v[pc->b.a[i]], p->v[i]);
+ copy_rvec(fr->x[pc->b.a[i]], p->x[i]);
}
- }
- if (p->f && fr->bF)
- {
- for (i = 0; i < pc->b.nra; ++i)
+ if (p->v && fr->bV)
{
- copy_rvec(fr->f[pc->b.a[i]], p->f[i]);
+ for (i = 0; i < pc->b.nra; ++i)
+ {
+ copy_rvec(fr->v[pc->b.a[i]], p->v[i]);
+ }
}
- }
- break;
- case POS_ALL:
- gmx_calc_comg(top, fr->x, pc->b.nra, pc->b.a, bMass, p->x[0]);
- if (p->v && fr->bV)
- {
- gmx_calc_comg(top, fr->v, pc->b.nra, pc->b.a, bMass, p->v[0]);
- }
- if (p->f && fr->bF)
- {
- gmx_calc_comg_f(top, fr->f, pc->b.nra, pc->b.a, bMass, p->f[0]);
- }
- break;
- case POS_ALL_PBC:
- gmx_calc_comg_pbc(top, fr->x, pbc, pc->b.nra, pc->b.a, bMass, p->x[0]);
- if (p->v && fr->bV)
- {
- gmx_calc_comg(top, fr->v, pc->b.nra, pc->b.a, bMass, p->v[0]);
- }
- if (p->f && fr->bF)
- {
- gmx_calc_comg_f(top, fr->f, pc->b.nra, pc->b.a, bMass, p->f[0]);
- }
- break;
- default:
- gmx_calc_comg_blocka(top, fr->x, &pc->b, bMass, p->x);
- if (p->v && fr->bV)
- {
- gmx_calc_comg_blocka(top, fr->v, &pc->b, bMass, p->v);
- }
- if (p->f && fr->bF)
- {
- gmx_calc_comg_blocka(top, fr->f, &pc->b, bMass, p->f);
- }
- break;
+ if (p->f && fr->bF)
+ {
+ for (i = 0; i < pc->b.nra; ++i)
+ {
+ copy_rvec(fr->f[pc->b.a[i]], p->f[i]);
+ }
+ }
+ break;
+ case POS_ALL:
+ gmx_calc_comg(top, fr->x, pc->b.nra, pc->b.a, bMass, p->x[0]);
+ if (p->v && fr->bV)
+ {
+ gmx_calc_comg(top, fr->v, pc->b.nra, pc->b.a, bMass, p->v[0]);
+ }
+ if (p->f && fr->bF)
+ {
+ gmx_calc_comg_f(top, fr->f, pc->b.nra, pc->b.a, bMass, p->f[0]);
+ }
+ break;
+ case POS_ALL_PBC:
+ gmx_calc_comg_pbc(top, fr->x, pbc, pc->b.nra, pc->b.a, bMass, p->x[0]);
+ if (p->v && fr->bV)
+ {
+ gmx_calc_comg(top, fr->v, pc->b.nra, pc->b.a, bMass, p->v[0]);
+ }
+ if (p->f && fr->bF)
+ {
+ gmx_calc_comg_f(top, fr->f, pc->b.nra, pc->b.a, bMass, p->f[0]);
+ }
+ break;
+ default:
+ gmx_calc_comg_blocka(top, fr->x, &pc->b, bMass, p->x);
+ if (p->v && fr->bV)
+ {
+ gmx_calc_comg_blocka(top, fr->v, &pc->b, bMass, p->v);
+ }
+ if (p->f && fr->bF)
+ {
+ gmx_calc_comg_f_blocka(top, fr->f, &pc->b, bMass, p->f);
+ }
+ break;
}
}
}