Sort all includes in src/gromacs
[alexxy/gromacs.git] / src / gromacs / fileio / tngio.cpp
index 44f68752bfa8d381a0c20aabfe7d442dcb8563d0..1303872d3d8364a0ee2a3405b1efc82a2d0c9ee8 100644 (file)
@@ -1,7 +1,7 @@
 /*
  * This file is part of the GROMACS molecular simulation package.
  *
- * Copyright (c) 2013, by the GROMACS development team, led by
+ * Copyright (c) 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.
  * To help us fund GROMACS development, we humbly ask that you cite
  * the research papers on the package. Check out http://www.gromacs.org.
  */
+#include "gmxpre.h"
+
 #include "tngio.h"
 
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
+#include "config.h"
 
 #ifdef HAVE_UNISTD_H
 #include <unistd.h>
 #endif
 
 #ifdef GMX_USE_TNG
-#include "../../external/tng_io/include/tng_io.h"
+#include "tng/tng_io.h"
 #endif
 
+#include "gromacs/fileio/gmxfio.h"
 #include "gromacs/legacyheaders/copyrite.h"
-#include "gromacs/legacyheaders/gmx_fatal.h"
-#include "gromacs/legacyheaders/main.h"
-#include "gromacs/legacyheaders/physics.h"
-#include "gromacs/utility/gmxassert.h"
-#include "gromacs/utility/programinfo.h"
+#include "gromacs/legacyheaders/types/ifunc.h"
+#include "gromacs/math/units.h"
 #include "gromacs/math/utilities.h"
-#include "gmxfio.h"
+#include "gromacs/topology/topology.h"
+#include "gromacs/utility/basenetwork.h"
+#include "gromacs/utility/common.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/gmxassert.h"
+#include "gromacs/utility/programcontext.h"
 
 static const char *modeToVerb(char mode)
 {
+    const char *p;
     switch (mode)
     {
         case 'r':
-            return "reading";
+            p = "reading";
             break;
         case 'w':
-            return "writing";
+            p = "writing";
             break;
         case 'a':
-            return "appending";
+            p = "appending";
             break;
         default:
             gmx_fatal(FARGS, "Invalid file opening mode %c", mode);
-            return "";
+            p = "";
+            break;
     }
+    return p;
 }
 
 void gmx_tng_open(const char       *filename,
@@ -129,7 +135,8 @@ void gmx_tng_open(const char       *filename,
 #ifdef GMX_DOUBLE
         precisionString = " (double precision)";
 #endif
-        sprintf(programInfo, "%.100s, %.128s%.24s", gmx::ProgramInfo::getInstance().displayName().c_str(),
+        sprintf(programInfo, "%.100s, %.128s%.24s",
+                gmx::getProgramContext().displayName(),
                 GromacsVersion(), precisionString);
         if (mode == 'w')
         {
@@ -142,12 +149,14 @@ void gmx_tng_open(const char       *filename,
 //             tng_last_program_name_set(*tng, programInfo);
 //         }
 
-#ifdef HAVE_UNISTD_H
+#if defined(HAVE_UNISTD_H) && !defined(__MINGW32__)
         char username[256];
-        getlogin_r(username, 256);
-        if (mode == 'w')
+        if (!getlogin_r(username, 256))
         {
-            tng_first_user_name_set(*tng, username);
+            if (mode == 'w')
+            {
+                tng_first_user_name_set(*tng, username);
+            }
         }
 /* TODO: This should be implemented when the above fixme is done (adding data to
  * the header). */
@@ -187,6 +196,9 @@ static void addTngMoleculeFromTopology(tng_trajectory_t     tng,
                                        gmx_int64_t          numMolecules,
                                        tng_molecule_t      *tngMol)
 {
+    tng_chain_t      tngChain = NULL;
+    tng_residue_t    tngRes   = NULL;
+
     if (tng_molecule_add(tng, moleculeName, tngMol) != TNG_SUCCESS)
     {
         gmx_file("Cannot add molecule to TNG molecular system.");
@@ -194,36 +206,44 @@ static void addTngMoleculeFromTopology(tng_trajectory_t     tng,
 
     /* FIXME: The TNG atoms should contain mass and atomB info (for free
      * energy calculations), i.e. in when it's available in TNG (2.0). */
-    for (int atomIt = 0; atomIt < atoms->nr; atomIt++)
+    for (int atomIndex = 0; atomIndex < atoms->nr; atomIndex++)
     {
-        const t_atom *at = &atoms->atom[atomIt];
+        const t_atom *at = &atoms->atom[atomIndex];
         /* FIXME: Currently the TNG API can only add atoms belonging to a
          * residue and chain. Wait for TNG 2.0*/
         if (atoms->nres > 0)
         {
             const t_resinfo *resInfo        = &atoms->resinfo[at->resind];
             char             chainName[2]   = {resInfo->chainid, 0};
-            tng_chain_t      tngChain       = NULL;
-            tng_residue_t    tngRes         = NULL;
             tng_atom_t       tngAtom        = NULL;
+            t_atom          *prevAtom;
 
-            if (tng_molecule_chain_find (tng, *tngMol, chainName,
-                                         (gmx_int64_t)-1, &tngChain) !=
-                TNG_SUCCESS)
+            if (atomIndex > 0)
+            {
+                prevAtom = &atoms->atom[atomIndex - 1];
+            }
+            else
             {
-                tng_molecule_chain_add (tng, *tngMol, chainName,
-                                        &tngChain);
+                prevAtom = 0;
             }
 
-            /* FIXME: When TNG supports both residue index and residue
-             * number the latter should be used. Wait for TNG 2.0*/
-            if (tng_chain_residue_find(tng, tngChain, *resInfo->name,
-                                       at->resind + 1, &tngRes)
-                != TNG_SUCCESS)
+            /* If this is the first atom or if the residue changed add the
+             * residue to the TNG molecular system. */
+            if (!prevAtom || resInfo != &atoms->resinfo[prevAtom->resind])
             {
+                /* If this is the first atom or if the chain changed add
+                 * the chain to the TNG molecular system. */
+                if (!prevAtom || resInfo->chainid !=
+                    atoms->resinfo[prevAtom->resind].chainid)
+                {
+                    tng_molecule_chain_add(tng, *tngMol, chainName,
+                                           &tngChain);
+                }
+                /* FIXME: When TNG supports both residue index and residue
+                 * number the latter should be used. Wait for TNG 2.0*/
                 tng_chain_residue_add(tng, tngChain, *resInfo->name, &tngRes);
             }
-            tng_residue_atom_add(tng, tngRes, *(atoms->atomname[atomIt]), *(atoms->atomtype[atomIt]), &tngAtom);
+            tng_residue_atom_add(tng, tngRes, *(atoms->atomname[atomIndex]), *(atoms->atomtype[atomIndex]), &tngAtom);
         }
     }
     tng_molecule_cnt_set(tng, *tngMol, numMolecules);
@@ -242,18 +262,18 @@ void gmx_tng_add_mtop(tng_trajectory_t  tng,
         return;
     }
 
-    for (int molIt = 0; molIt < mtop->nmolblock; molIt++)
+    for (int molIndex = 0; molIndex < mtop->nmolblock; molIndex++)
     {
         tng_molecule_t       tngMol  = NULL;
         const gmx_moltype_t *molType =
-            &mtop->moltype[mtop->molblock[molIt].type];
+            &mtop->moltype[mtop->molblock[molIndex].type];
 
         /* Add a molecule to the TNG trajectory with the same name as the
          * current molecule. */
         addTngMoleculeFromTopology(tng,
                                    *(molType->name),
                                    &molType->atoms,
-                                   mtop->molblock[molIt].nmol,
+                                   mtop->molblock[molIndex].nmol,
                                    &tngMol);
 
         /* Bonds have to be deduced from interactions (constraints etc). Different
@@ -485,8 +505,44 @@ void gmx_tng_prepare_md_writing(tng_trajectory_t  tng,
 }
 
 #ifdef GMX_USE_TNG
-/* Create a TNG molecule representing the selection groups
- * to write */
+/* Check if all atoms in the molecule system are selected
+ * by a selection group of type specified by gtype. */
+static gmx_bool all_atoms_selected(const gmx_mtop_t *mtop,
+                                   const int         gtype)
+{
+    const gmx_moltype_t     *molType;
+    const t_atoms           *atoms;
+
+    /* Iterate through all atoms in the molecule system and
+     * check if they belong to a selection group of the
+     * requested type. */
+    for (int molIndex = 0, i = 0; molIndex < mtop->nmoltype; molIndex++)
+    {
+        molType = &mtop->moltype[mtop->molblock[molIndex].type];
+
+        atoms = &molType->atoms;
+
+        for (int j = 0; j < mtop->molblock[molIndex].nmol; j++)
+        {
+            for (int atomIndex = 0; atomIndex < atoms->nr; atomIndex++, i++)
+            {
+                if (ggrpnr(&mtop->groups, gtype, i) != 0)
+                {
+                    return FALSE;
+                }
+            }
+        }
+    }
+    return TRUE;
+}
+
+/* Create TNG molecules which will represent each of the selection
+ * groups for writing. But do that only if there is actually a
+ * specified selection group and it is not the whole system.
+ * TODO: Currently the only selection that is taken into account
+ * is egcCompressedX, but other selections should be added when
+ * e.g. writing energies is implemented.
+ */
 static void add_selection_groups(tng_trajectory_t  tng,
                                  const gmx_mtop_t *mtop)
 {
@@ -495,7 +551,7 @@ static void add_selection_groups(tng_trajectory_t  tng,
     const t_atom            *at;
     const t_resinfo         *resInfo;
     const t_ilist           *ilist;
-    int                      nAtoms      = 0, i = 0, j, molIt, atomIt, nameIndex;
+    int                      nameIndex;
     int                      atom_offset = 0;
     tng_molecule_t           mol, iterMol;
     tng_chain_t              chain;
@@ -505,6 +561,25 @@ static void add_selection_groups(tng_trajectory_t  tng,
     gmx_int64_t              nMols;
     char                    *groupName;
 
+    /* TODO: When the TNG molecules block is more flexible TNG selection
+     * groups should not need all atoms specified. It should be possible
+     * just to specify what molecules are selected (and/or which atoms in
+     * the molecule) and how many (if applicable). */
+
+    /* If no atoms are selected we do not need to create a
+     * TNG selection group molecule. */
+    if (mtop->groups.ngrpnr[egcCompressedX] == 0)
+    {
+        return;
+    }
+
+    /* If all atoms are selected we do not have to create a selection
+     * group molecule in the TNG molecule system. */
+    if (all_atoms_selected(mtop, egcCompressedX))
+    {
+        return;
+    }
+
     /* The name of the TNG molecule containing the selection group is the
      * same as the name of the selection group. */
     nameIndex = *mtop->groups.grps[egcCompressedX].nm_ind;
@@ -513,16 +588,16 @@ static void add_selection_groups(tng_trajectory_t  tng,
     tng_molecule_alloc(tng, &mol);
     tng_molecule_name_set(tng, mol, groupName);
     tng_molecule_chain_add(tng, mol, "", &chain);
-    for (molIt = 0; molIt < mtop->nmoltype; molIt++)
+    for (int molIndex = 0, i = 0; molIndex < mtop->nmoltype; molIndex++)
     {
-        molType = &mtop->moltype[mtop->molblock[molIt].type];
+        molType = &mtop->moltype[mtop->molblock[molIndex].type];
 
         atoms = &molType->atoms;
 
-        for (j = 0; j < mtop->molblock[molIt].nmol; j++)
+        for (int j = 0; j < mtop->molblock[molIndex].nmol; j++)
         {
             bool bAtomsAdded = FALSE;
-            for (atomIt = 0; atomIt < atoms->nr; atomIt++, i++)
+            for (int atomIndex = 0; atomIndex < atoms->nr; atomIndex++, i++)
             {
                 char *res_name;
                 int   res_id;
@@ -531,7 +606,7 @@ static void add_selection_groups(tng_trajectory_t  tng,
                 {
                     continue;
                 }
-                at = &atoms->atom[atomIt];
+                at = &atoms->atom[atomIndex];
                 if (atoms->nres > 0)
                 {
                     resInfo = &atoms->resinfo[at->resind];
@@ -552,10 +627,9 @@ static void add_selection_groups(tng_trajectory_t  tng,
                      * original residue IDs - otherwise there might be conflicts. */
                     tng_chain_residue_add(tng, chain, res_name, &res);
                 }
-                tng_residue_atom_w_id_add(tng, res, *(atoms->atomname[atomIt]),
-                                          *(atoms->atomtype[atomIt]),
-                                          atom_offset + atomIt, &atom);
-                nAtoms++;
+                tng_residue_atom_w_id_add(tng, res, *(atoms->atomname[atomIndex]),
+                                          *(atoms->atomtype[atomIndex]),
+                                          atom_offset + atomIndex, &atom);
                 bAtomsAdded = TRUE;
             }
             /* Add bonds. */
@@ -616,24 +690,17 @@ static void add_selection_groups(tng_trajectory_t  tng,
             atom_offset += atoms->nr;
         }
     }
-    if (nAtoms != i)
+    tng_molecule_existing_add(tng, &mol);
+    tng_molecule_cnt_set(tng, mol, 1);
+    tng_num_molecule_types_get(tng, &nMols);
+    for (gmx_int64_t k = 0; k < nMols; k++)
     {
-        tng_molecule_existing_add(tng, &mol);
-        tng_molecule_cnt_set(tng, mol, 1);
-        tng_num_molecule_types_get(tng, &nMols);
-        for (gmx_int64_t k = 0; k < nMols; k++)
+        tng_molecule_of_index_get(tng, k, &iterMol);
+        if (iterMol == mol)
         {
-            tng_molecule_of_index_get(tng, k, &iterMol);
-            if (iterMol == mol)
-            {
-                continue;
-            }
-            tng_molecule_cnt_set(tng, iterMol, 0);
+            continue;
         }
-    }
-    else
-    {
-        tng_molecule_free(tng, &mol);
+        tng_molecule_cnt_set(tng, iterMol, 0);
     }
 }
 #endif
@@ -642,7 +709,7 @@ void gmx_tng_set_compression_precision(tng_trajectory_t tng,
                                        real             prec)
 {
 #ifdef GMX_USE_TNG
-    tng_compression_precision_set(tng, 1.0/prec);
+    tng_compression_precision_set(tng, prec);
 #else
     GMX_UNUSED_VALUE(tng);
     GMX_UNUSED_VALUE(prec);
@@ -692,9 +759,9 @@ void gmx_fwrite_tng(tng_trajectory_t tng,
 #else
     static write_data_func_pointer           write_data           = tng_util_generic_with_time_write;
 #endif
-    double                  elapsedSeconds = elapsedPicoSeconds * PICO;
-    gmx_int64_t             nParticles;
-    char                    compression;
+    double                                   elapsedSeconds = elapsedPicoSeconds * PICO;
+    gmx_int64_t                              nParticles;
+    char                                     compression;
 
 
     if (!tng)