src/gromacs/nbnxm/benchmark/bench_system.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2019, by the GROMACS development team, led by
   5  * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
   6  * and including many others, as listed in the AUTHORS file in the
   7  * top-level source directory and at http://www.gromacs.org.
   8  *
   9  * GROMACS is free software; you can redistribute it and/or
  10  * modify it under the terms of the GNU Lesser General Public License
  11  * as published by the Free Software Foundation; either version 2.1
  12  * of the License, or (at your option) any later version.
  13  *
  14  * GROMACS is distributed in the hope that it will be useful,
  15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
  16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  17  * Lesser General Public License for more details.
  18  *
  19  * You should have received a copy of the GNU Lesser General Public
  20  * License along with GROMACS; if not, see
  21  * http://www.gnu.org/licenses, or write to the Free Software Foundation,
  22  * Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA.
  23  *
  24  * If you want to redistribute modifications to GROMACS, please
  25  * consider that scientific software is very special. Version
  26  * control is crucial - bugs must be traceable. We will be happy to
  27  * consider code for inclusion in the official distribution, but
  28  * derived work must not be called official GROMACS. Details are found
  29  * in the README & COPYING files - if they are missing, get the
  30  * official version at http://www.gromacs.org.
  31  *
  32  * To help us fund GROMACS development, we humbly ask that you cite
  33  * the research papers on the package. Check out http://www.gromacs.org.
  34  */
  35
  36 /*! \internal \file
  37  * \brief
  38  * This file defines functions for setting up a benchmark system
  39  *
  40  * \author Berk Hess <hess@kth.se>
  41  * \ingroup module_nbnxm
  42  */
  43
  44 #include "gmxpre.h"
  45
  46 #include "bench_system.h"
  47
  48 #include <numeric>
  49 #include <vector>
  50
  51 #include "gromacs/math/vec.h"
  52 #include "gromacs/mdlib/dispersioncorrection.h"
  53 #include "gromacs/mdtypes/forcerec.h"
  54 #include "gromacs/nbnxm/nbnxm.h"
  55 #include "gromacs/pbcutil/ishift.h"
  56 #include "gromacs/pbcutil/pbc.h"
  57 #include "gromacs/utility/fatalerror.h"
  58
  59 #include "bench_coords.h"
  60
  61 namespace gmx
  62 {
  63
  64 namespace
  65 {
  66
  67 // A 3-site water model
  68 //! The number of atoms in a molecule
  69 constexpr int numAtomsInMolecule = 3;
  70 //! The atom type of the oxygen atom
  71 constexpr int typeOxygen = 0;
  72 //! The atom type of the hydrogen atom
  73 constexpr int typeHydrogen = 1;
  74 //! The charge of the oxygen atom
  75 constexpr real chargeOxygen = -0.8476;
  76 //! The charge of the hydrogen atom
  77 constexpr real chargeHydrogen = 0.4238;
  78 //! The LJ C6 parameter of the Oxygen atom
  79 constexpr real c6Oxygen = 0.0026173456;
  80 //! The LJ C12 parameter of the Oxygen atom
  81 constexpr real c12Oxygen = 2.634129e-06;
  82 // Note that the hydrogen has LJ parameters all zero
  83
  84 } // namespace
  85
  86 //! Generates coordinates and a box for the base system scaled by \p multiplicationFactor
  87 //
  88 // The parameter \p multiplicationFactor should be a power of 2.
  89 // A fatal error is generated when this is not the case.
  90 static void generateCoordinates(int multiplicationFactor, std::vector<gmx::RVec>* coordinates, matrix box)
  91 {
  92     if (multiplicationFactor < 1 || (multiplicationFactor & (multiplicationFactor - 1)) != 0)
  93     {
  94         gmx_fatal(FARGS, "The size factor has to be a power of 2");
  95     }
  96
  97     if (multiplicationFactor == 1)
  98     {
  99         *coordinates = coordinates1000;
 100         copy_mat(box1000, box);
 101
 102         return;
 103     }
 104
 105     ivec factors = { 1, 1, 1 };
 106
 107     int dim = 0;
 108     while (multiplicationFactor > 1)
 109     {
 110         factors[dim] *= 2;
 111         multiplicationFactor /= 2;
 112         dim++;
 113         if (dim == DIM)
 114         {
 115             dim = 0;
 116         }
 117     }
 118     printf("Stacking a box of %zu atoms %d x %d x %d times\n", coordinates1000.size(), factors[XX],
 119            factors[YY], factors[ZZ]);
 120
 121     coordinates->resize(factors[XX] * factors[YY] * factors[ZZ] * coordinates1000.size());
 122
 123     int       i = 0;
 124     gmx::RVec shift;
 125     for (int x = 0; x < factors[XX]; x++)
 126     {
 127         shift[XX] = x * box1000[XX][XX];
 128         for (int y = 0; y < factors[YY]; y++)
 129         {
 130             shift[YY] = y * box1000[YY][YY];
 131             for (int z = 0; z < factors[ZZ]; z++)
 132             {
 133                 shift[ZZ] = z * box1000[ZZ][ZZ];
 134
 135                 for (const gmx::RVec& coordOrig : coordinates1000)
 136                 {
 137                     (*coordinates)[i] = coordOrig + shift;
 138                     i++;
 139                 }
 140             }
 141         }
 142     }
 143
 144     for (int d1 = 0; d1 < DIM; d1++)
 145     {
 146         for (int d2 = 0; d2 < DIM; d2++)
 147         {
 148             box[d1][d2] = factors[d1] * box1000[d1][d2];
 149         }
 150     }
 151 }
 152
 153 BenchmarkSystem::BenchmarkSystem(const int multiplicationFactor)
 154 {
 155     numAtomTypes = 2;
 156     nonbondedParameters.resize(numAtomTypes * numAtomTypes * 2, 0);
 157     nonbondedParameters[0] = c6Oxygen;
 158     nonbondedParameters[1] = c12Oxygen;
 159
 160     generateCoordinates(multiplicationFactor, &coordinates, box);
 161     put_atoms_in_box(epbcXYZ, box, coordinates);
 162
 163     int numAtoms = coordinates.size();
 164     GMX_RELEASE_ASSERT(numAtoms % numAtomsInMolecule == 0,
 165                        "Coordinates should match whole molecules");
 166
 167     atomTypes.resize(numAtoms);
 168     charges.resize(numAtoms);
 169     atomInfoAllVdw.resize(numAtoms);
 170     atomInfoOxygenVdw.resize(numAtoms);
 171
 172     for (int a = 0; a < numAtoms; a++)
 173     {
 174         if (a % numAtomsInMolecule == 0)
 175         {
 176             // Oxgygen
 177             atomTypes[a] = typeOxygen;
 178             charges[a]   = chargeOxygen;
 179             SET_CGINFO_HAS_VDW(atomInfoAllVdw[a]);
 180             SET_CGINFO_HAS_VDW(atomInfoOxygenVdw[a]);
 181         }
 182         else
 183         {
 184             // Hydrogen
 185             atomTypes[a] = typeHydrogen;
 186             charges[a]   = chargeHydrogen;
 187             SET_CGINFO_HAS_VDW(atomInfoAllVdw[a]);
 188         }
 189         SET_CGINFO_HAS_Q(atomInfoAllVdw[a]);
 190         SET_CGINFO_HAS_Q(atomInfoOxygenVdw[a]);
 191
 192         excls.pushBackListOfSize(numAtomsInMolecule);
 193         gmx::ArrayRef<int> exclusionsForAtom   = excls.back();
 194         const int          firstAtomInMolecule = a - (a % numAtomsInMolecule);
 195         std::iota(exclusionsForAtom.begin(), exclusionsForAtom.end(), firstAtomInMolecule);
 196     }
 197
 198     forceRec.ntype = numAtomTypes;
 199     forceRec.nbfp  = nonbondedParameters;
 200     snew(forceRec.shift_vec, SHIFTS);
 201     calc_shifts(box, forceRec.shift_vec);
 202 }
 203
 204 } // namespace gmx