f2dbaf4c1badf37f056e950aba1802af7d585242
[alexxy/gromacs.git] / src / gromacs / nbnxm / gridset.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  *
38  * \brief
39  * Implements the GridSet class.
40  *
41  * \author Berk Hess <hess@kth.se>
42  * \ingroup module_nbnxm
43  */
44
45 #include "gmxpre.h"
46
47 #include "gridset.h"
48
49 #include "gromacs/mdlib/gmx_omp_nthreads.h"
50 #include "gromacs/mdlib/updategroupscog.h"
51 #include "gromacs/nbnxm/atomdata.h"
52 #include "gromacs/utility/fatalerror.h"
53
54 namespace Nbnxm
55 {
56
57 //! Returns the number of DD zones
58 static int numDDZones(const std::array<bool, DIM> &haveDomDecPerDim)
59 {
60     int  numDDZones = 1;
61     for (auto haveDD : haveDomDecPerDim)
62     {
63         if (haveDD)
64         {
65             numDDZones *= 2;
66         }
67     }
68
69     return numDDZones;
70 }
71
72 GridSet::GridSet(const std::array<bool, DIM> &haveDomDecPerDim,
73                  const PairlistType           pairlistType,
74                  const bool                   haveFep,
75                  const int                    numThreads) :
76     grids_(numDDZones(haveDomDecPerDim), pairlistType),
77     haveFep_(haveFep),
78     numRealAtomsLocal_(0),
79     numRealAtomsTotal_(0),
80     gridWork_(numThreads)
81 {
82     clear_mat(box_);
83 }
84
85 void GridSet::setLocalAtomOrder()
86 {
87     /* Set the atom order for the home cell (index 0) */
88     const Nbnxm::Grid &grid = grids_[0];
89
90     int                atomIndex = 0;
91     for (int cxy = 0; cxy < grid.numColumns(); cxy++)
92     {
93         const int numAtoms  = grid.numAtomsInColumn(cxy);
94         int       cellIndex = grid.firstCellInColumn(cxy)*grid.geometry().numAtomsPerCell;
95         for (int i = 0; i < numAtoms; i++)
96         {
97             atomIndices_[cellIndex] = atomIndex;
98             cells_[atomIndex]       = cellIndex;
99             atomIndex++;
100             cellIndex++;
101         }
102     }
103 }
104
105 // TODO: Move to gridset.cpp
106 void GridSet::putOnGrid(const matrix                    box,
107                         const int                       ddZone,
108                         const rvec                      lowerCorner,
109                         const rvec                      upperCorner,
110                         const gmx::UpdateGroupsCog     *updateGroupsCog,
111                         const int                       atomStart,
112                         const int                       atomEnd,
113                         real                            atomDensity,
114                         const int                      *atinfo,
115                         gmx::ArrayRef<const gmx::RVec>  x,
116                         const int                       numAtomsMoved,
117                         const int                      *move,
118                         nbnxn_atomdata_t               *nbat)
119 {
120     Nbnxm::Grid  &grid = grids_[ddZone];
121
122     int           cellOffset;
123     if (ddZone == 0)
124     {
125         cellOffset = 0;
126     }
127     else
128     {
129         const Nbnxm::Grid &previousGrid = grids_[ddZone - 1];
130         cellOffset = previousGrid.atomIndexEnd()/previousGrid.geometry().numAtomsPerCell;
131     }
132
133     const int n = atomEnd - atomStart;
134
135     real      maxAtomGroupRadius;
136     if (ddZone == 0)
137     {
138         copy_mat(box, box_);
139
140         numRealAtomsLocal_ = atomEnd - numAtomsMoved;
141         /* We assume that nbnxn_put_on_grid is called first
142          * for the local atoms (ddZone=0).
143          */
144         numRealAtomsTotal_ = atomEnd - numAtomsMoved;
145
146         maxAtomGroupRadius = (updateGroupsCog ? updateGroupsCog->maxUpdateGroupRadius() : 0);
147
148         if (debug)
149         {
150             fprintf(debug, "natoms_local = %5d atom_density = %5.1f\n",
151                     numRealAtomsLocal_, atomDensity);
152         }
153     }
154     else
155     {
156         const Nbnxm::Grid::Dimensions &dimsGrid0 = grids_[0].dimensions();
157         atomDensity        = dimsGrid0.atomDensity;
158         maxAtomGroupRadius = dimsGrid0.maxAtomGroupRadius;
159
160         numRealAtomsTotal_ = std::max(numRealAtomsTotal_, atomEnd);
161     }
162
163     /* We always use the home zone (grid[0]) for setting the cell size,
164      * since determining densities for non-local zones is difficult.
165      */
166     grid.setDimensions(ddZone, n - numAtomsMoved,
167                        lowerCorner, upperCorner,
168                        atomDensity,
169                        maxAtomGroupRadius,
170                        haveFep_);
171
172     for (GridWork &work : gridWork_)
173     {
174         work.numAtomsPerColumn.resize(grid.numColumns() + 1);
175     }
176
177     /* Make space for the new cell indices */
178     cells_.resize(atomEnd);
179
180     const int nthread = gmx_omp_nthreads_get(emntPairsearch);
181
182 #pragma omp parallel for num_threads(nthread) schedule(static)
183     for (int thread = 0; thread < nthread; thread++)
184     {
185         try
186         {
187             Grid::calcColumnIndices(grid.dimensions(),
188                                     updateGroupsCog,
189                                     atomStart, atomEnd, x,
190                                     ddZone, move, thread, nthread,
191                                     cells_, gridWork_[thread].numAtomsPerColumn);
192         }
193         GMX_CATCH_ALL_AND_EXIT_WITH_FATAL_ERROR;
194     }
195
196     GridSetData gridSetData = getGridSetData();
197
198     /* Copy the already computed cell indices to the grid and sort, when needed */
199     grid.setCellIndices(ddZone, cellOffset, &gridSetData, gridWork_,
200                         atomStart, atomEnd, atinfo, x, numAtomsMoved, nbat);
201
202     if (ddZone == 0)
203     {
204         nbat->natoms_local = nbat->numAtoms();
205     }
206     if (ddZone == gmx::ssize(grids_) - 1)
207     {
208         /* We are done setting up all grids, we can resize the force buffers */
209         nbat->resizeForceBuffers();
210     }
211 }
212
213 } // namespace Nbnxm