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 655 : LinkCells::LinkCells( Communicator& cc ) :
29 655 : comm(cc),
30 655 : cutoffwasset(false),
31 655 : nopbc(false),
32 655 : link_cutoff(0.0),
33 655 : ncells(3),
34 1310 : nstride(3) {
35 655 : }
36 :
37 655 : void LinkCells::setCutoff( const double& lcut ) {
38 655 : cutoffwasset=true;
39 655 : link_cutoff=lcut;
40 655 : }
41 :
42 513 : double LinkCells::getCutoff() const {
43 513 : plumed_assert( cutoffwasset );
44 513 : return link_cutoff;
45 : }
46 :
47 12094 : void LinkCells::buildCellLists( const std::vector<Vector>& pos, const std::vector<unsigned>& indices, const Pbc& pbc ) {
48 12094 : plumed_assert( cutoffwasset && pos.size()==indices.size() );
49 :
50 : // Create an orthorhombic box around the atomic positions that encompasses every atomic position if there are no pbc
51 12094 : auto box = pbc.getBox();
52 12094 : 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) {
53 129 : Vector minp, maxp;
54 129 : minp = maxp = pos[0];
55 516 : for(unsigned k=0; k<3; ++k) {
56 20262 : for(unsigned i=1; i<pos.size(); ++i) {
57 19875 : if( pos[i][k]>maxp[k] ) {
58 1875 : maxp[k] = pos[i][k];
59 : }
60 19875 : if( pos[i][k]<minp[k] ) {
61 0 : minp[k] = pos[i][k];
62 : }
63 : }
64 387 : if( link_cutoff<std::sqrt(std::numeric_limits<double>::max()) ) {
65 381 : box[k][k] = link_cutoff*( 1 + std::ceil( (maxp[k] - minp[k])/link_cutoff ) );
66 : } else {
67 6 : box[k][k] = maxp[k] - minp[k] + 1;
68 : }
69 387 : origin[k] = ( minp[k] + maxp[k] ) / 2;
70 : }
71 129 : nopbc=true;
72 : // Setup the pbc object by copying it from action if the Pbc were set
73 : } else {
74 11965 : auto determinant = box.determinant();
75 11965 : nopbc=false;
76 11965 : plumed_assert(determinant > epsilon) <<"Cell lists cannot be built when passing a box with null volume. Volume is "<<determinant;
77 : }
78 12094 : mypbc.setBox( box );
79 :
80 : // Setup the lists
81 12094 : if( pos.size()!=allcells.size() ) {
82 585 : allcells.resize( pos.size() );
83 585 : lcell_lists.resize( pos.size() );
84 : }
85 :
86 : {
87 : // This is the reciprocal lattice
88 : // notice that reciprocal.getRow(0) is a vector that is orthogonal to b and c
89 : // This allows to use linked cells in non orthorhomic boxes
90 12094 : Tensor reciprocal(transpose(mypbc.getInvBox()));
91 12094 : ncells[0] = std::floor( 1.0/ reciprocal.getRow(0).modulo() / link_cutoff );
92 12094 : if( ncells[0]==0 ) {
93 11058 : ncells[0]=1;
94 : }
95 12094 : ncells[1] = std::floor( 1.0/ reciprocal.getRow(1).modulo() / link_cutoff );
96 12094 : if( ncells[1]==0 ) {
97 11058 : ncells[1]=1;
98 : }
99 12094 : ncells[2] = std::floor( 1.0/ reciprocal.getRow(2).modulo() / link_cutoff );
100 12094 : if( ncells[2]==0 ) {
101 11058 : ncells[2]=1;
102 : }
103 : }
104 : // Setup the strides
105 12094 : nstride[0]=1;
106 12094 : nstride[1]=ncells[0];
107 12094 : nstride[2]=ncells[0]*ncells[1];
108 :
109 : // Setup the storage for link cells
110 12094 : unsigned ncellstot=ncells[0]*ncells[1]*ncells[2];
111 12094 : if( lcell_tots.size()!=ncellstot ) {
112 583 : lcell_tots.resize( ncellstot );
113 583 : lcell_starts.resize( ncellstot );
114 : }
115 : // Clear nlcells
116 12094 : lcell_tots.assign( lcell_tots.size(), 0 );
117 : // Clear allcells
118 12094 : allcells.assign( allcells.size(), 0 );
119 :
120 : // Find out what cell everyone is in
121 12094 : unsigned rank=comm.Get_rank(), size=comm.Get_size();
122 1202201 : for(unsigned i=rank; i<pos.size(); i+=size) {
123 1190107 : allcells[i]=findCell( pos[i] );
124 1190107 : lcell_tots[allcells[i]]++;
125 : }
126 : // And gather all this information on every node
127 12094 : comm.Sum( allcells );
128 12094 : comm.Sum( lcell_tots );
129 :
130 : // Now prepare the link cell lists
131 12094 : unsigned tot=0;
132 1550283 : for(unsigned i=0; i<lcell_tots.size(); ++i) {
133 1538189 : lcell_starts[i]=tot;
134 1538189 : tot+=lcell_tots[i];
135 1538189 : lcell_tots[i]=0;
136 : }
137 12094 : plumed_assert( tot==pos.size() ) <<"Total number of atoms found in link cells is "<<tot<<" number of atoms is "<<pos.size();
138 :
139 : // And setup the link cells properly
140 1274607 : for(unsigned j=0; j<pos.size(); ++j) {
141 1262513 : unsigned myind = lcell_starts[ allcells[j] ] + lcell_tots[ allcells[j] ];
142 1262513 : lcell_lists[ myind ] = indices[j];
143 1262513 : lcell_tots[allcells[j]]++;
144 : }
145 12094 : }
146 :
147 : #define LINKC_MIN(n) ((n<2)? 0 : -1)
148 : #define LINKC_MAX(n) ((n<3)? 1 : 2)
149 : #define LINKC_PBC(n,num) ((n<0)? num-1 : n%num )
150 :
151 266161 : void LinkCells::addRequiredCells( const std::array<unsigned,3>& celn, unsigned& ncells_required,
152 : std::vector<unsigned>& cells_required ) const {
153 : unsigned nnew_cells=0;
154 1398271 : for(int nx=LINKC_MIN(ncells[0]); nx<LINKC_MAX(ncells[0]); ++nx) {
155 468785 : int xval = celn[0] + nx;
156 468785 : xval=LINKC_PBC(xval,ncells[0])*nstride[0];
157 3195397 : for(int ny=LINKC_MIN(ncells[1]); ny<LINKC_MAX(ncells[1]); ++ny) {
158 1071047 : int yval = celn[1] + ny;
159 1071047 : yval=LINKC_PBC(yval,ncells[1])*nstride[1];
160 8536012 : for(int nz=LINKC_MIN(ncells[2]); nz<LINKC_MAX(ncells[2]); ++nz) {
161 2858634 : int zval = celn[2] + nz;
162 2858634 : zval=LINKC_PBC(zval,ncells[2])*nstride[2];
163 :
164 2858634 : unsigned mybox=xval+yval+zval;
165 : bool added=false;
166 2858634 : for(unsigned k=0; k<ncells_required; ++k) {
167 0 : if( mybox==cells_required[k] ) {
168 : added=true;
169 : break;
170 : }
171 : }
172 2858634 : if( !added ) {
173 2858634 : cells_required[ncells_required+nnew_cells]=mybox;
174 2858634 : nnew_cells++;
175 : }
176 : }
177 : }
178 : }
179 266161 : ncells_required += nnew_cells;
180 266161 : }
181 :
182 6750 : void LinkCells::retrieveNeighboringAtoms( const Vector& pos, std::vector<unsigned>& cell_list,
183 : unsigned& natomsper, std::vector<unsigned>& atoms ) const {
184 6750 : if( cell_list.size()!=getNumberOfCells() ) {
185 243 : cell_list.resize( getNumberOfCells() );
186 : }
187 6750 : unsigned ncellt=0;
188 6750 : addRequiredCells( findMyCell( pos ), ncellt, cell_list );
189 6750 : retrieveAtomsInCells( ncellt, cell_list, natomsper, atoms );
190 6750 : }
191 :
192 266161 : void LinkCells::retrieveAtomsInCells( const unsigned& ncells_required,
193 : const std::vector<unsigned>& cells_required,
194 : unsigned& natomsper, std::vector<unsigned>& atoms ) const {
195 3124795 : for(unsigned i=0; i<ncells_required; ++i) {
196 2858634 : unsigned mybox=cells_required[i];
197 26260493 : for(unsigned k=0; k<lcell_tots[mybox]; ++k) {
198 23401859 : unsigned myatom = lcell_lists[lcell_starts[mybox]+k];
199 23401859 : if( myatom!=atoms[0] ) { // Ideally would provide an option to not do this
200 23148066 : atoms[natomsper]=myatom;
201 23148066 : natomsper++;
202 : }
203 : }
204 : }
205 266161 : }
206 :
207 1456268 : std::array<unsigned,3> LinkCells::findMyCell( const Vector& pos ) const {
208 : std::array<unsigned,3> celn;
209 1456268 : Vector mypos = pos;
210 1456268 : if( nopbc ) {
211 10137 : mypos = pos - origin;
212 : }
213 1456268 : Vector fpos=mypbc.realToScaled( mypos );
214 5825072 : for(unsigned j=0; j<3; ++j) {
215 4368804 : celn[j] = std::floor( ( Tools::pbc(fpos[j]) + 0.5 ) * ncells[j] );
216 4368804 : plumed_assert( celn[j]>=0 && celn[j]<ncells[j] ) <<"in link cell "<<celn[j]<<" but should be between 0 and "<<ncells[j]<<" link cell cutoff is "<<link_cutoff<<" position is "<<fpos[0]<<" "<<fpos[1]<<" "<<fpos[2]<<" box is "<<mypbc.getBox()(0,0)<<" "<<mypbc.getBox()(1,1)<<" "<<mypbc.getBox()(2,2);
217 : }
218 1456268 : return celn;
219 : }
220 :
221 1190107 : unsigned LinkCells::convertIndicesToIndex( const unsigned& nx, const unsigned& ny, const unsigned& nz ) const {
222 1190107 : return nx*nstride[0] + ny*nstride[1] + nz*nstride[2];
223 : }
224 :
225 1190107 : unsigned LinkCells::findCell( const Vector& pos ) const {
226 1190107 : std::array<unsigned,3> celn( findMyCell(pos ) );
227 1190107 : return convertIndicesToIndex( celn[0], celn[1], celn[2] );
228 : }
229 :
230 11750 : unsigned LinkCells::getMaxInCell() const {
231 11750 : unsigned maxn = lcell_tots[0];
232 1527297 : for(unsigned i=1; i<lcell_tots.size(); ++i) {
233 1515547 : if( lcell_tots[i]>maxn ) {
234 : maxn=lcell_tots[i];
235 : }
236 : }
237 11750 : return maxn;
238 : }
239 :
240 : }
|
