src/gromacs/nbnxm/benchmark/bench_system.cpp

   1 /*
   2  * This file is part of the GROMACS molecular simulation package.
   3  *
   4  * Copyright (c) 2019,2020,2021, 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 (!gmx::isPowerOfTwo(multiplicationFactor))
  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",
 119            coordinates1000.size(),
 120            factors[XX],
 121            factors[YY],
 122            factors[ZZ]);
 123
 124     coordinates->resize(factors[XX] * factors[YY] * factors[ZZ] * coordinates1000.size());
 125
 126     int       i = 0;
 127     gmx::RVec shift;
 128     for (int x = 0; x < factors[XX]; x++)
 129     {
 130         shift[XX] = x * box1000[XX][XX];
 131         for (int y = 0; y < factors[YY]; y++)
 132         {
 133             shift[YY] = y * box1000[YY][YY];
 134             for (int z = 0; z < factors[ZZ]; z++)
 135             {
 136                 shift[ZZ] = z * box1000[ZZ][ZZ];
 137
 138                 for (const gmx::RVec& coordOrig : coordinates1000)
 139                 {
 140                     (*coordinates)[i] = coordOrig + shift;
 141                     i++;
 142                 }
 143             }
 144         }
 145     }
 146
 147     for (int d1 = 0; d1 < DIM; d1++)
 148     {
 149         for (int d2 = 0; d2 < DIM; d2++)
 150         {
 151             box[d1][d2] = factors[d1] * box1000[d1][d2];
 152         }
 153     }
 154 }
 155
 156 BenchmarkSystem::BenchmarkSystem(const int multiplicationFactor, const std::string& outputFile)
 157 {
 158     numAtomTypes = 2;
 159     nonbondedParameters.resize(numAtomTypes * numAtomTypes * 2, 0);
 160     nonbondedParameters[0] = c6Oxygen;
 161     nonbondedParameters[1] = c12Oxygen;
 162
 163     generateCoordinates(multiplicationFactor, &coordinates, box);
 164     put_atoms_in_box(PbcType::Xyz, box, coordinates);
 165
 166     int numAtoms = coordinates.size();
 167     GMX_RELEASE_ASSERT(numAtoms % numAtomsInMolecule == 0,
 168                        "Coordinates should match whole molecules");
 169
 170     atomTypes.resize(numAtoms);
 171     charges.resize(numAtoms);
 172     atomInfoAllVdw.resize(numAtoms);
 173     atomInfoOxygenVdw.resize(numAtoms);
 174
 175     for (int a = 0; a < numAtoms; a++)
 176     {
 177         if (a % numAtomsInMolecule == 0)
 178         {
 179             // Oxgygen
 180             atomTypes[a] = typeOxygen;
 181             charges[a]   = chargeOxygen;
 182             SET_CGINFO_HAS_VDW(atomInfoAllVdw[a]);
 183             SET_CGINFO_HAS_VDW(atomInfoOxygenVdw[a]);
 184         }
 185         else
 186         {
 187             // Hydrogen
 188             atomTypes[a] = typeHydrogen;
 189             charges[a]   = chargeHydrogen;
 190             SET_CGINFO_HAS_VDW(atomInfoAllVdw[a]);
 191         }
 192         SET_CGINFO_HAS_Q(atomInfoAllVdw[a]);
 193         SET_CGINFO_HAS_Q(atomInfoOxygenVdw[a]);
 194
 195         excls.pushBackListOfSize(numAtomsInMolecule);
 196         gmx::ArrayRef<int> exclusionsForAtom   = excls.back();
 197         const int          firstAtomInMolecule = a - (a % numAtomsInMolecule);
 198         std::iota(exclusionsForAtom.begin(), exclusionsForAtom.end(), firstAtomInMolecule);
 199     }
 200
 201     forceRec.ntype = numAtomTypes;
 202     forceRec.nbfp  = nonbondedParameters;
 203     forceRec.shift_vec.resize(SHIFTS);
 204     calc_shifts(box, forceRec.shift_vec);
 205     if (!outputFile.empty())
 206     {
 207         csv = fopen(outputFile.c_str(), "w+");
 208     }
 209 }
 210
 211 } // namespace gmx