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 658 : LinkCells::LinkCells( Communicator& cc ) :
29 658 : comm(cc),
30 658 : cutoffwasset(false),
31 658 : link_cutoff(0.0),
32 658 : ncells(3),
33 1316 : nstride(3) {
34 658 : }
35 :
36 658 : void LinkCells::setCutoff( const double& lcut ) {
37 658 : cutoffwasset=true;
38 658 : link_cutoff=lcut;
39 658 : }
40 :
41 518 : double LinkCells::getCutoff() const {
42 518 : plumed_assert( cutoffwasset );
43 518 : return link_cutoff;
44 : }
45 :
46 11949 : void LinkCells::buildCellLists( const std::vector<Vector>& pos, const std::vector<unsigned>& indices, const Pbc& pbc ) {
47 11949 : plumed_assert( cutoffwasset && pos.size()==indices.size() );
48 :
49 : // Must be able to check that pbcs are not nonsensical in some way?? -- GAT
50 11949 : auto box = pbc.getBox();
51 11949 : 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 11945 : auto determinant = box.determinant();
59 11945 : 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 11949 : mypbc.setBox( box );
63 :
64 : // Setup the lists
65 11949 : if( pos.size()!=allcells.size() ) {
66 462 : allcells.resize( pos.size() );
67 462 : 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 11949 : Tensor reciprocal(transpose(mypbc.getInvBox()));
75 11949 : ncells[0] = std::floor( 1.0/ reciprocal.getRow(0).modulo() / link_cutoff );
76 11949 : if( ncells[0]==0 ) {
77 11051 : ncells[0]=1;
78 : }
79 11949 : ncells[1] = std::floor( 1.0/ reciprocal.getRow(1).modulo() / link_cutoff );
80 11949 : if( ncells[1]==0 ) {
81 11051 : ncells[1]=1;
82 : }
83 11949 : ncells[2] = std::floor( 1.0/ reciprocal.getRow(2).modulo() / link_cutoff );
84 11949 : if( ncells[2]==0 ) {
85 11051 : ncells[2]=1;
86 : }
87 : }
88 : // Setup the strides
89 11949 : nstride[0]=1;
90 11949 : nstride[1]=ncells[0];
91 11949 : nstride[2]=ncells[0]*ncells[1];
92 :
93 : // Setup the storage for link cells
94 11949 : unsigned ncellstot=ncells[0]*ncells[1]*ncells[2];
95 11949 : if( lcell_tots.size()!=ncellstot ) {
96 466 : lcell_tots.resize( ncellstot );
97 466 : lcell_starts.resize( ncellstot );
98 : }
99 : // Clear nlcells
100 11949 : lcell_tots.assign( lcell_tots.size(), 0 );
101 : // Clear allcells
102 11949 : allcells.assign( allcells.size(), 0 );
103 :
104 : // Find out what cell everyone is in
105 11949 : unsigned rank=comm.Get_rank(), size=comm.Get_size();
106 1108641 : for(unsigned i=rank; i<pos.size(); i+=size) {
107 1096692 : allcells[i]=findCell( pos[i] );
108 1096692 : lcell_tots[allcells[i]]++;
109 : }
110 : // And gather all this information on every node
111 11949 : comm.Sum( allcells );
112 11949 : comm.Sum( lcell_tots );
113 :
114 : // Now prepare the link cell lists
115 11949 : unsigned tot=0;
116 1485187 : for(unsigned i=0; i<lcell_tots.size(); ++i) {
117 1473238 : lcell_starts[i]=tot;
118 1473238 : tot+=lcell_tots[i];
119 1473238 : lcell_tots[i]=0;
120 : }
121 11949 : 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 1181047 : for(unsigned j=0; j<pos.size(); ++j) {
125 1169098 : unsigned myind = lcell_starts[ allcells[j] ] + lcell_tots[ allcells[j] ];
126 1169098 : lcell_lists[ myind ] = indices[j];
127 1169098 : lcell_tots[allcells[j]]++;
128 : }
129 11949 : }
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 518329 : 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 3670176 : for(int nx=LINKC_MIN(ncells[0]); nx<LINKC_MAX(ncells[0]); ++nx) {
139 1226518 : int xval = celn[0] + nx;
140 1226518 : xval=LINKC_PBC(xval,ncells[0])*nstride[0];
141 10015842 : for(int ny=LINKC_MIN(ncells[1]); ny<LINKC_MAX(ncells[1]); ++ny) {
142 3345166 : int yval = celn[1] + ny;
143 3345166 : yval=LINKC_PBC(yval,ncells[1])*nstride[1];
144 29004071 : for(int nz=LINKC_MIN(ncells[2]); nz<LINKC_MAX(ncells[2]); ++nz) {
145 9681939 : int zval = celn[2] + nz;
146 9681939 : zval=LINKC_PBC(zval,ncells[2])*nstride[2];
147 :
148 9681939 : unsigned mybox=xval+yval+zval;
149 : bool added=false;
150 9681939 : for(unsigned k=0; k<ncells_required; ++k) {
151 0 : if( mybox==cells_required[k] ) {
152 : added=true;
153 : break;
154 : }
155 : }
156 9681939 : if( !added ) {
157 9681939 : cells_required[ncells_required+nnew_cells]=mybox;
158 9681939 : nnew_cells++;
159 : }
160 : }
161 : }
162 : }
163 518329 : ncells_required += nnew_cells;
164 518329 : }
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 518329 : void LinkCells::retrieveAtomsInCells( const unsigned& ncells_required,
177 : const std::vector<unsigned>& cells_required,
178 : unsigned& natomsper, std::vector<unsigned>& atoms ) const {
179 10200268 : for(unsigned i=0; i<ncells_required; ++i) {
180 9681939 : unsigned mybox=cells_required[i];
181 43377970 : for(unsigned k=0; k<lcell_tots[mybox]; ++k) {
182 33696031 : unsigned myatom = lcell_lists[lcell_starts[mybox]+k];
183 33696031 : if( myatom!=atoms[0] ) { // Ideally would provide an option to not do this
184 33186703 : atoms[natomsper]=myatom;
185 33186703 : natomsper++;
186 : }
187 : }
188 : }
189 518329 : }
190 :
191 1615021 : std::array<unsigned,3> LinkCells::findMyCell( const Vector& pos ) const {
192 : std::array<unsigned,3> celn;
193 1615021 : if( ncells[0]*ncells[1]*ncells[2] == 1 ) {
194 364258 : celn[0]=celn[1]=celn[2]=0;
195 364258 : return celn;
196 : }
197 1250763 : Vector fpos=mypbc.realToScaled( pos );
198 5003052 : for(unsigned j=0; j<3; ++j) {
199 3752289 : celn[j] = std::floor( ( Tools::pbc(fpos[j]) + 0.5 ) * ncells[j] );
200 3752289 : plumed_assert( celn[j]>=0 && celn[j]<ncells[j] ); // Check that atom is in box
201 : }
202 1250763 : return celn;
203 : }
204 :
205 1096692 : unsigned LinkCells::convertIndicesToIndex( const unsigned& nx, const unsigned& ny, const unsigned& nz ) const {
206 1096692 : return nx*nstride[0] + ny*nstride[1] + nz*nstride[2];
207 : }
208 :
209 1096692 : unsigned LinkCells::findCell( const Vector& pos ) const {
210 1096692 : std::array<unsigned,3> celn( findMyCell(pos ) );
211 1096692 : return convertIndicesToIndex( celn[0], celn[1], celn[2] );
212 : }
213 :
214 11730 : unsigned LinkCells::getMaxInCell() const {
215 11730 : unsigned maxn = lcell_tots[0];
216 1469716 : for(unsigned i=1; i<lcell_tots.size(); ++i) {
217 1457986 : if( lcell_tots[i]>maxn ) {
218 : maxn=lcell_tots[i];
219 : }
220 : }
221 11730 : return maxn;
222 : }
223 :
224 : }
|
