Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2014-2023 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed 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
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #include "LinkCells.h"
23 : #include "Communicator.h"
24 : #include "Tools.h"
25 :
26 : namespace PLMD {
27 :
28 654 : LinkCells::LinkCells( Communicator& cc ) :
29 654 : comm(cc),
30 654 : cutoffwasset(false),
31 654 : link_cutoff(0.0),
32 654 : ncells(3),
33 1308 : nstride(3) {
34 654 : }
35 :
36 654 : void LinkCells::setCutoff( const double& lcut ) {
37 654 : cutoffwasset=true;
38 654 : link_cutoff=lcut;
39 654 : }
40 :
41 513 : double LinkCells::getCutoff() const {
42 513 : plumed_assert( cutoffwasset );
43 513 : return link_cutoff;
44 : }
45 :
46 11969 : void LinkCells::buildCellLists( const std::vector<Vector>& pos, const std::vector<unsigned>& indices, const Pbc& pbc ) {
47 11969 : plumed_assert( cutoffwasset && pos.size()==indices.size() );
48 :
49 : // Must be able to check that pbcs are not nonsensical in some way?? -- GAT
50 11969 : auto box = pbc.getBox();
51 11969 : if(box(0,0)==0.0 && box(0,1)==0.0 && box(0,2)==0.0 && box(1,0)==0.0 && box(1,1)==0.0 && box(1,2)==0.0 && box(2,0)==0.0 && box(2,1)==0 && box(2,2)==0) {
52 : // If the box is not set then we can't use link cells. We thus set the link cell cutoff and box vectors equal to 23 (because it is the best number).
53 : // Setting everything this way ensures that the link cells are a 1x1x1 box. Notice that if it is a one by one by one box then we are hard coded to return
54 : // 0 in findCell. GAT
55 : // Note: any non infinite - non zero number should work correctly here! Apparently Gareth likes numerology. I do as well, so let's keep this. GB
56 4 : box(0,0) = box(1,1) = box(2,2) = link_cutoff = 23;
57 : } else {
58 11965 : auto determinant = box.determinant();
59 11965 : plumed_assert(determinant > epsilon) <<"Cell lists cannot be built when passing a box with null volume. Volume is "<<determinant;
60 : }
61 : // Setup the pbc object by copying it from action
62 11969 : mypbc.setBox( box );
63 :
64 : // Setup the lists
65 11969 : if( pos.size()!=allcells.size() ) {
66 460 : allcells.resize( pos.size() );
67 460 : lcell_lists.resize( pos.size() );
68 : }
69 :
70 : {
71 : // This is the reciprocal lattice
72 : // notice that reciprocal.getRow(0) is a vector that is orthogonal to b and c
73 : // This allows to use linked cells in non orthorhomic boxes
74 11969 : Tensor reciprocal(transpose(mypbc.getInvBox()));
75 11969 : ncells[0] = std::floor( 1.0/ reciprocal.getRow(0).modulo() / link_cutoff );
76 11969 : if( ncells[0]==0 ) {
77 11056 : ncells[0]=1;
78 : }
79 11969 : ncells[1] = std::floor( 1.0/ reciprocal.getRow(1).modulo() / link_cutoff );
80 11969 : if( ncells[1]==0 ) {
81 11056 : ncells[1]=1;
82 : }
83 11969 : ncells[2] = std::floor( 1.0/ reciprocal.getRow(2).modulo() / link_cutoff );
84 11969 : if( ncells[2]==0 ) {
85 11056 : ncells[2]=1;
86 : }
87 : }
88 : // Setup the strides
89 11969 : nstride[0]=1;
90 11969 : nstride[1]=ncells[0];
91 11969 : nstride[2]=ncells[0]*ncells[1];
92 :
93 : // Setup the storage for link cells
94 11969 : unsigned ncellstot=ncells[0]*ncells[1]*ncells[2];
95 11969 : if( lcell_tots.size()!=ncellstot ) {
96 464 : lcell_tots.resize( ncellstot );
97 464 : lcell_starts.resize( ncellstot );
98 : }
99 : // Clear nlcells
100 11969 : lcell_tots.assign( lcell_tots.size(), 0 );
101 : // Clear allcells
102 11969 : allcells.assign( allcells.size(), 0 );
103 :
104 : // Find out what cell everyone is in
105 11969 : unsigned rank=comm.Get_rank(), size=comm.Get_size();
106 1195326 : for(unsigned i=rank; i<pos.size(); i+=size) {
107 1183357 : allcells[i]=findCell( pos[i] );
108 1183357 : lcell_tots[allcells[i]]++;
109 : }
110 : // And gather all this information on every node
111 11969 : comm.Sum( allcells );
112 11969 : comm.Sum( lcell_tots );
113 :
114 : // Now prepare the link cell lists
115 11969 : unsigned tot=0;
116 1542788 : for(unsigned i=0; i<lcell_tots.size(); ++i) {
117 1530819 : lcell_starts[i]=tot;
118 1530819 : tot+=lcell_tots[i];
119 1530819 : lcell_tots[i]=0;
120 : }
121 11969 : plumed_assert( tot==pos.size() ) <<"Total number of atoms found in link cells is "<<tot<<" number of atoms is "<<pos.size();
122 :
123 : // And setup the link cells properly
124 1267732 : for(unsigned j=0; j<pos.size(); ++j) {
125 1255763 : unsigned myind = lcell_starts[ allcells[j] ] + lcell_tots[ allcells[j] ];
126 1255763 : lcell_lists[ myind ] = indices[j];
127 1255763 : lcell_tots[allcells[j]]++;
128 : }
129 11969 : }
130 :
131 : #define LINKC_MIN(n) ((n<2)? 0 : -1)
132 : #define LINKC_MAX(n) ((n<3)? 1 : 2)
133 : #define LINKC_PBC(n,num) ((n<0)? num-1 : n%num )
134 :
135 262786 : void LinkCells::addRequiredCells( const std::array<unsigned,3>& celn, unsigned& ncells_required,
136 : std::vector<unsigned>& cells_required ) const {
137 : unsigned nnew_cells=0;
138 1369741 : for(int nx=LINKC_MIN(ncells[0]); nx<LINKC_MAX(ncells[0]); ++nx) {
139 459029 : int xval = celn[0] + nx;
140 459029 : xval=LINKC_PBC(xval,ncells[0])*nstride[0];
141 3112973 : for(int ny=LINKC_MIN(ncells[1]); ny<LINKC_MAX(ncells[1]); ++ny) {
142 1042855 : int yval = celn[1] + ny;
143 1042855 : yval=LINKC_PBC(yval,ncells[1])*nstride[1];
144 8297964 : for(int nz=LINKC_MIN(ncells[2]); nz<LINKC_MAX(ncells[2]); ++nz) {
145 2777194 : int zval = celn[2] + nz;
146 2777194 : zval=LINKC_PBC(zval,ncells[2])*nstride[2];
147 :
148 2777194 : unsigned mybox=xval+yval+zval;
149 : bool added=false;
150 2777194 : for(unsigned k=0; k<ncells_required; ++k) {
151 0 : if( mybox==cells_required[k] ) {
152 : added=true;
153 : break;
154 : }
155 : }
156 2777194 : if( !added ) {
157 2777194 : cells_required[ncells_required+nnew_cells]=mybox;
158 2777194 : nnew_cells++;
159 : }
160 : }
161 : }
162 : }
163 262786 : ncells_required += nnew_cells;
164 262786 : }
165 :
166 3375 : void LinkCells::retrieveNeighboringAtoms( const Vector& pos, std::vector<unsigned>& cell_list,
167 : unsigned& natomsper, std::vector<unsigned>& atoms ) const {
168 3375 : if( cell_list.size()!=getNumberOfCells() ) {
169 124 : cell_list.resize( getNumberOfCells() );
170 : }
171 3375 : unsigned ncellt=0;
172 3375 : addRequiredCells( findMyCell( pos ), ncellt, cell_list );
173 3375 : retrieveAtomsInCells( ncellt, cell_list, natomsper, atoms );
174 3375 : }
175 :
176 262786 : void LinkCells::retrieveAtomsInCells( const unsigned& ncells_required,
177 : const std::vector<unsigned>& cells_required,
178 : unsigned& natomsper, std::vector<unsigned>& atoms ) const {
179 3039980 : for(unsigned i=0; i<ncells_required; ++i) {
180 2777194 : unsigned mybox=cells_required[i];
181 26085109 : for(unsigned k=0; k<lcell_tots[mybox]; ++k) {
182 23307915 : unsigned myatom = lcell_lists[lcell_starts[mybox]+k];
183 23307915 : if( myatom!=atoms[0] ) { // Ideally would provide an option to not do this
184 23054122 : atoms[natomsper]=myatom;
185 23054122 : natomsper++;
186 : }
187 : }
188 : }
189 262786 : }
190 :
191 1446143 : std::array<unsigned,3> LinkCells::findMyCell( const Vector& pos ) const {
192 : std::array<unsigned,3> celn;
193 1446143 : if( ncells[0]*ncells[1]*ncells[2] == 1 ) {
194 366260 : celn[0]=celn[1]=celn[2]=0;
195 366260 : return celn;
196 : }
197 1079883 : Vector fpos=mypbc.realToScaled( pos );
198 4319532 : for(unsigned j=0; j<3; ++j) {
199 3239649 : celn[j] = std::floor( ( Tools::pbc(fpos[j]) + 0.5 ) * ncells[j] );
200 3239649 : plumed_assert( celn[j]>=0 && celn[j]<ncells[j] ); // Check that atom is in box
201 : }
202 1079883 : return celn;
203 : }
204 :
205 1183357 : unsigned LinkCells::convertIndicesToIndex( const unsigned& nx, const unsigned& ny, const unsigned& nz ) const {
206 1183357 : return nx*nstride[0] + ny*nstride[1] + nz*nstride[2];
207 : }
208 :
209 1183357 : unsigned LinkCells::findCell( const Vector& pos ) const {
210 1183357 : std::array<unsigned,3> celn( findMyCell(pos ) );
211 1183357 : return convertIndicesToIndex( celn[0], celn[1], celn[2] );
212 : }
213 :
214 11750 : unsigned LinkCells::getMaxInCell() const {
215 11750 : unsigned maxn = lcell_tots[0];
216 1527297 : for(unsigned i=1; i<lcell_tots.size(); ++i) {
217 1515547 : if( lcell_tots[i]>maxn ) {
218 : maxn=lcell_tots[i];
219 : }
220 : }
221 11750 : return maxn;
222 : }
223 :
224 : }
|
