Sort all includes in src/gromacs
[alexxy/gromacs.git] / src / gromacs / selection / poscalc.cpp
index 1b9debee8f6492c9687bb265bb2aa599b61e5356..f193abd42c577b1445fc4d4707a106080528671e 100644 (file)
@@ -1,32 +1,36 @@
 /*
+ * 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
 {
@@ -184,13 +192,13 @@ struct gmx_ana_poscalc_t
      * 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;
 };
@@ -217,7 +225,7 @@ const char * const gmx::PositionCalculationCollection::typeEnumValues[] = {
 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;
@@ -249,18 +257,18 @@ PositionCalculationCollection::typeFromEnum(const char *post,
     {
         *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')
@@ -517,7 +525,7 @@ PositionCalculationCollection::printTree(FILE *fp) const
             gmx_ana_poscalc_t *base;
 
             fprintf(fp, "   Base: ");
-            j = 1;
+            j    = 1;
             base = impl_->first_;
             while (base && base != pc->sbase)
             {
@@ -550,11 +558,30 @@ gmx_ana_poscalc_t *
 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_)
@@ -614,7 +641,7 @@ void PositionCalculationCollection::initFrame()
     while (pc)
     {
         pc->bEval = false;
-        pc = pc->next;
+        pc        = pc->next;
     }
 }
 
@@ -637,7 +664,7 @@ set_poscalc_maxindex(gmx_ana_poscalc_t *pc, gmx_ana_index_t *g, bool bBase)
     /* 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
@@ -654,12 +681,10 @@ set_poscalc_maxindex(gmx_ana_poscalc_t *pc, gmx_ana_index_t *g, bool bBase)
     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);
     }
 }
 
@@ -724,7 +749,7 @@ should_merge(gmx_ana_poscalc_t *pc1, gmx_ana_poscalc_t *pc2,
         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)
@@ -765,10 +790,10 @@ create_simple_base(gmx_ana_poscalc_t *pc)
     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);
@@ -797,8 +822,8 @@ merge_to_base(gmx_ana_poscalc_t *base, gmx_ana_poscalc_t *pc)
     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)
     {
@@ -822,10 +847,10 @@ merge_to_base(gmx_ana_poscalc_t *base, gmx_ana_poscalc_t *pc)
         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)
         {
@@ -853,7 +878,7 @@ merge_to_base(gmx_ana_poscalc_t *base, gmx_ana_poscalc_t *pc)
         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);
     }
 }
 
@@ -908,11 +933,11 @@ setup_base(gmx_ana_poscalc_t *pc)
         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 */
@@ -959,7 +984,7 @@ setup_base(gmx_ana_poscalc_t *pc)
                     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 */
@@ -1041,7 +1066,7 @@ gmx_ana_poscalc_init_pos(gmx_ana_poscalc_t *pc, gmx_ana_pos_t *p)
     {
         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);
@@ -1050,8 +1075,6 @@ gmx_ana_poscalc_init_pos(gmx_ana_poscalc_t *pc, gmx_ana_pos_t *p)
     {
         gmx_ana_pos_reserve_forces(p);
     }
-    gmx_ana_pos_set_nr(p, p->m.nr);
-    gmx_ana_pos_set_evalgrp(p, &pc->gmax);
 }
 
 /*!
@@ -1086,10 +1109,7 @@ gmx_ana_poscalc_free(gmx_ana_poscalc_t *pc)
     {
         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);
@@ -1106,7 +1126,8 @@ gmx_ana_poscalc_free(gmx_ana_poscalc_t *pc)
 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;
     }
@@ -1129,7 +1150,7 @@ gmx_ana_poscalc_update(gmx_ana_poscalc_t *pc, gmx_ana_pos_t *p,
                        gmx_ana_index_t *g, t_trxframe *fr, t_pbc *pbc)
 {
     int  i, bi, bj;
-    
+
     if (pc->bEval == true && !(pc->flags & POS_MASKONLY))
     {
         return;
@@ -1146,19 +1167,19 @@ gmx_ana_poscalc_update(gmx_ana_poscalc_t *pc, gmx_ana_pos_t *p,
     {
         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))
     {
@@ -1172,14 +1193,14 @@ gmx_ana_poscalc_update(gmx_ana_poscalc_t *pc, gmx_ana_pos_t *p,
          * 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]);
@@ -1187,7 +1208,7 @@ gmx_ana_poscalc_update(gmx_ana_poscalc_t *pc, gmx_ana_pos_t *p,
             }
             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]);
@@ -1196,14 +1217,14 @@ gmx_ana_poscalc_update(gmx_ana_poscalc_t *pc, gmx_ana_pos_t *p,
         }
         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]);
@@ -1211,7 +1232,7 @@ gmx_ana_poscalc_update(gmx_ana_poscalc_t *pc, gmx_ana_pos_t *p,
             }
             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]);
@@ -1244,63 +1265,63 @@ gmx_ana_poscalc_update(gmx_ana_poscalc_t *pc, gmx_ana_pos_t *p,
         }
         /* 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;
         }
     }
 }