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 <vector>
  49
  50 #include "gromacs/math/vec.h"
  51 #include "gromacs/mdlib/dispersioncorrection.h"
  52 #include "gromacs/mdtypes/forcerec.h"
  53 #include "gromacs/nbnxm/nbnxm.h"
  54 #include "gromacs/pbcutil/ishift.h"
  55 #include "gromacs/pbcutil/pbc.h"
  56 #include "gromacs/utility/fatalerror.h"
  57
  58 #include "bench_coords.h"
  59
  60 namespace gmx
  61 {
  62
  63 namespace
  64 {
  65
  66 // A 3-site water model
  67 //! The number of atoms in a molecule
  68 constexpr int numAtomsInMolecule = 3;
  69 //! The atom type of the oxygen atom
  70 constexpr int typeOxygen = 0;
  71 //! The atom type of the hydrogen atom
  72 constexpr int typeHydrogen = 1;
  73 //! The charge of the oxygen atom
  74 constexpr real chargeOxygen = -0.8476;
  75 //! The charge of the hydrogen atom
  76 constexpr real chargeHydrogen = 0.4238;
  77 //! The LJ C6 parameter of the Oxygen atom
  78 constexpr real c6Oxygen = 0.0026173456;
  79 //! The LJ C12 parameter of the Oxygen atom
  80 constexpr real c12Oxygen = 2.634129e-06;
  81 // Note that the hydrogen has LJ parameters all zero
  82
  83 } // namespace
  84
  85 //! Generates coordinates and a box for the base system scaled by \p multiplicationFactor
  86 //
  87 // The parameter \p multiplicationFactor should be a power of 2.
  88 // A fatal error is generated when this is not the case.
  89 static void generateCoordinates(int multiplicationFactor, std::vector<gmx::RVec>* coordinates, matrix box)
  90 {
  91     if (multiplicationFactor < 1 || (multiplicationFactor & (multiplicationFactor - 1)) != 0)
  92     {
  93         gmx_fatal(FARGS, "The size factor has to be a power of 2");
  94     }
  95
  96     if (multiplicationFactor == 1)
  97     {
  98         *coordinates = coordinates1000;
  99         copy_mat(box1000, box);
 100
 101         return;
 102     }
 103
 104     ivec factors = { 1, 1, 1 };
 105
 106     int dim = 0;
 107     while (multiplicationFactor > 1)
 108     {
 109         factors[dim] *= 2;
 110         multiplicationFactor /= 2;
 111         dim++;
 112         if (dim == DIM)
 113         {
 114             dim = 0;
 115         }
 116     }
 117     printf("Stacking a box of %zu atoms %d x %d x %d times\n", coordinates1000.size(), factors[XX],
 118            factors[YY], factors[ZZ]);
 119
 120     coordinates->resize(factors[XX] * factors[YY] * factors[ZZ] * coordinates1000.size());
 121
 122     int       i = 0;
 123     gmx::RVec shift;
 124     for (int x = 0; x < factors[XX]; x++)
 125     {
 126         shift[XX] = x * box1000[XX][XX];
 127         for (int y = 0; y < factors[YY]; y++)
 128         {
 129             shift[YY] = y * box1000[YY][YY];
 130             for (int z = 0; z < factors[ZZ]; z++)
 131             {
 132                 shift[ZZ] = z * box1000[ZZ][ZZ];
 133
 134                 for (const gmx::RVec& coordOrig : coordinates1000)
 135                 {
 136                     (*coordinates)[i] = coordOrig + shift;
 137                     i++;
 138                 }
 139             }
 140         }
 141     }
 142
 143     for (int d1 = 0; d1 < DIM; d1++)
 144     {
 145         for (int d2 = 0; d2 < DIM; d2++)
 146         {
 147             box[d1][d2] = factors[d1] * box1000[d1][d2];
 148         }
 149     }
 150 }
 151
 152 BenchmarkSystem::BenchmarkSystem(const int multiplicationFactor)
 153 {
 154     numAtomTypes = 2;
 155     nonbondedParameters.resize(numAtomTypes * numAtomTypes * 2, 0);
 156     nonbondedParameters[0] = c6Oxygen;
 157     nonbondedParameters[1] = c12Oxygen;
 158
 159     generateCoordinates(multiplicationFactor, &coordinates, box);
 160     put_atoms_in_box(epbcXYZ, box, coordinates);
 161
 162     int numAtoms = coordinates.size();
 163     GMX_RELEASE_ASSERT(numAtoms % numAtomsInMolecule == 0,
 164                        "Coordinates should match whole molecules");
 165
 166     atomTypes.resize(numAtoms);
 167     charges.resize(numAtoms);
 168     atomInfoAllVdw.resize(numAtoms);
 169     atomInfoOxygenVdw.resize(numAtoms);
 170     snew(excls.index, numAtoms + 1);
 171     snew(excls.a, numAtoms * numAtomsInMolecule);
 172     excls.index[0] = 0;
 173
 174     for (int a = 0; a < numAtoms; a++)
 175     {
 176         if (a % numAtomsInMolecule == 0)
 177         {
 178             // Oxgygen
 179             atomTypes[a] = typeOxygen;
 180             charges[a]   = chargeOxygen;
 181             SET_CGINFO_HAS_VDW(atomInfoAllVdw[a]);
 182             SET_CGINFO_HAS_VDW(atomInfoOxygenVdw[a]);
 183         }
 184         else
 185         {
 186             // Hydrogen
 187             atomTypes[a] = typeHydrogen;
 188             charges[a]   = chargeHydrogen;
 189             SET_CGINFO_HAS_VDW(atomInfoAllVdw[a]);
 190         }
 191         SET_CGINFO_HAS_Q(atomInfoAllVdw[a]);
 192         SET_CGINFO_HAS_Q(atomInfoOxygenVdw[a]);
 193
 194         const int firstAtomInMolecule = a - (a % numAtomsInMolecule);
 195         for (int aj = 0; aj < numAtomsInMolecule; aj++)
 196         {
 197             excls.a[a * numAtomsInMolecule + aj] = firstAtomInMolecule + aj;
 198         }
 199         excls.index[a + 1] = (a + 1) * numAtomsInMolecule;
 200     }
 201
 202     forceRec.ntype = numAtomTypes;
 203     forceRec.nbfp  = nonbondedParameters;
 204     snew(forceRec.shift_vec, SHIFTS);
 205     calc_shifts(box, forceRec.shift_vec);
 206 }
 207
 208 } // namespace gmx