Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2014-2019 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 : comm(cc),
30 : cutoffwasset(false),
31 : link_cutoff(0.0),
32 : ncells(3),
33 658 : nstride(3)
34 : {
35 658 : }
36 :
37 352 : void LinkCells::setCutoff( const double& lcut ) {
38 352 : cutoffwasset=true; link_cutoff=lcut;
39 352 : }
40 :
41 8 : double LinkCells::getCutoff() const {
42 8 : plumed_assert( cutoffwasset ); return link_cutoff;
43 : }
44 :
45 1693 : void LinkCells::buildCellLists( const std::vector<Vector>& pos, const std::vector<unsigned>& indices, const Pbc& pbc ) {
46 3386 : plumed_assert( cutoffwasset && pos.size()==indices.size() );
47 :
48 : // Must be able to check that pbcs are not nonsensical in some way?? -- GAT
49 :
50 : // Setup the pbc object by copying it from action
51 1693 : mypbc.setBox( pbc.getBox() );
52 :
53 : // Setup the lists
54 1693 : if( pos.size()!=allcells.size() ) {
55 840 : allcells.resize( pos.size() ); lcell_lists.resize( pos.size() );
56 : }
57 :
58 : {
59 : // This is the reciprocal lattice
60 : // notice that reciprocal.getRow(0) is a vector that is orthogonal to b and c
61 : // This allows to use linked cells in non orthorhomic boxes
62 1693 : Tensor reciprocal(transpose(mypbc.getInvBox()));
63 3386 : ncells[0] = std::floor( 1.0/ reciprocal.getRow(0).modulo() / link_cutoff );
64 1693 : if( ncells[0]==0 ) ncells[0]=1;
65 3386 : ncells[1] = std::floor( 1.0/ reciprocal.getRow(1).modulo() / link_cutoff );
66 1693 : if( ncells[1]==0 ) ncells[1]=1;
67 3386 : ncells[2] = std::floor( 1.0/ reciprocal.getRow(2).modulo() / link_cutoff );
68 1693 : if( ncells[2]==0 ) ncells[2]=1;
69 : }
70 : // Setup the strides
71 6772 : nstride[0]=1; nstride[1]=ncells[0]; nstride[2]=ncells[0]*ncells[1];
72 :
73 : // Setup the storage for link cells
74 1693 : unsigned ncellstot=ncells[0]*ncells[1]*ncells[2];
75 1693 : if( lcell_tots.size()!=ncellstot ) {
76 409 : lcell_tots.resize( ncellstot ); lcell_starts.resize( ncellstot );
77 : }
78 : // Clear nlcells
79 601493 : for(unsigned i=0; i<ncellstot; ++i) lcell_tots[i]=0;
80 : // Clear allcells
81 3386 : allcells.assign( allcells.size(), 0 );
82 :
83 : // Find out what cell everyone is in
84 1693 : unsigned rank=comm.Get_rank(), size=comm.Get_size();
85 1239509 : for(unsigned i=rank; i<pos.size(); i+=size) {
86 412041 : allcells[i]=findCell( pos[i] );
87 824082 : lcell_tots[allcells[i]]++;
88 : }
89 : // And gather all this information on every node
90 1693 : comm.Sum( allcells ); comm.Sum( lcell_tots );
91 :
92 : // Now prepare the link cell lists
93 : unsigned tot=0;
94 1202986 : for(unsigned i=0; i<lcell_tots.size(); ++i) { lcell_starts[i]=tot; tot+=lcell_tots[i]; lcell_tots[i]=0; }
95 3386 : plumed_assert( tot==pos.size() );
96 :
97 : // And setup the link cells properly
98 1375346 : for(unsigned j=0; j<pos.size(); ++j) {
99 1371960 : unsigned myind = lcell_starts[ allcells[j] ] + lcell_tots[ allcells[j] ];
100 914640 : lcell_lists[ myind ] = indices[j];
101 914640 : lcell_tots[allcells[j]]++;
102 : }
103 1693 : }
104 :
105 : #define LINKC_MIN(n) ((n<2)? 0 : -1)
106 : #define LINKC_MAX(n) ((n<3)? 1 : 2)
107 : #define LINKC_PBC(n,num) ((n<0)? num-1 : n%num )
108 :
109 249696 : void LinkCells::addRequiredCells( const std::vector<unsigned>& celn, unsigned& ncells_required,
110 : std::vector<unsigned>& cells_required ) const {
111 : unsigned nnew_cells=0;
112 1301554 : for(int nx=LINKC_MIN(ncells[0]); nx<LINKC_MAX(ncells[0]); ++nx) {
113 401081 : int xval = celn[0] + nx;
114 802162 : xval=LINKC_PBC(xval,ncells[0])*nstride[0];
115 2504948 : for(int ny=LINKC_MIN(ncells[1]); ny<LINKC_MAX(ncells[1]); ++ny) {
116 851393 : int yval = celn[1] + ny;
117 1702786 : yval=LINKC_PBC(yval,ncells[1])*nstride[1];
118 6074728 : for(int nz=LINKC_MIN(ncells[2]); nz<LINKC_MAX(ncells[2]); ++nz) {
119 2185971 : int zval = celn[2] + nz;
120 4371942 : zval=LINKC_PBC(zval,ncells[2])*nstride[2];
121 :
122 2185971 : unsigned mybox=xval+yval+zval; bool added=false;
123 2185971 : for(unsigned k=0; k<ncells_required; ++k) {
124 0 : if( mybox==cells_required[k] ) { added=true; break; }
125 : }
126 4371942 : if( !added ) { cells_required[ncells_required+nnew_cells]=mybox; nnew_cells++; }
127 : }
128 : }
129 : }
130 249696 : ncells_required += nnew_cells;
131 249696 : }
132 :
133 14884 : void LinkCells::retrieveNeighboringAtoms( const Vector& pos, std::vector<unsigned>& cell_list,
134 : unsigned& natomsper, std::vector<unsigned>& atoms ) const {
135 14884 : if( cell_list.size()!=getNumberOfCells() ) cell_list.resize( getNumberOfCells() );
136 29768 : unsigned ncellt=0; addRequiredCells( findMyCell( pos ), ncellt, cell_list );
137 14884 : retrieveAtomsInCells( ncellt, cell_list, natomsper, atoms );
138 14884 : }
139 :
140 249696 : void LinkCells::retrieveAtomsInCells( const unsigned& ncells_required,
141 : const std::vector<unsigned>& cells_required,
142 : unsigned& natomsper, std::vector<unsigned>& atoms ) const {
143 249696 : plumed_assert( natomsper==1 || natomsper==2 ); // This is really a bug. If you are trying to reuse this ask GAT for help
144 4621638 : for(unsigned i=0; i<ncells_required; ++i) {
145 4371942 : unsigned mybox=cells_required[i];
146 1209196941 : for(unsigned k=0; k<lcell_tots[mybox]; ++k) {
147 803216666 : unsigned myatom = lcell_lists[lcell_starts[mybox]+k];
148 401608333 : if( myatom!=atoms[0] ) { // Ideally would provide an option to not do this
149 802906958 : atoms[natomsper]=myatom;
150 401453479 : natomsper++;
151 : }
152 : }
153 : }
154 249696 : }
155 :
156 661737 : std::vector<unsigned> LinkCells::findMyCell( const Vector& pos ) const {
157 661737 : Vector fpos=mypbc.realToScaled( pos ); std::vector<unsigned> celn(3);
158 4632159 : for(unsigned j=0; j<3; ++j) {
159 7940844 : celn[j] = std::floor( ( Tools::pbc(fpos[j]) + 0.5 ) * ncells[j] );
160 3970422 : plumed_assert( celn[j]>=0 && celn[j]<ncells[j] ); // Check that atom is in box
161 : }
162 661737 : return celn;
163 : }
164 :
165 412041 : unsigned LinkCells::convertIndicesToIndex( const unsigned& nx, const unsigned& ny, const unsigned& nz ) const {
166 824082 : return nx*nstride[0] + ny*nstride[1] + nz*nstride[2];
167 : }
168 :
169 412041 : unsigned LinkCells::findCell( const Vector& pos ) const {
170 412041 : std::vector<unsigned> celn( findMyCell(pos ) );
171 824082 : return convertIndicesToIndex( celn[0], celn[1], celn[2] );
172 : }
173 :
174 :
175 4839 : }
|
