Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-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 "CoordinationBase.h"
23 : #include <array>
24 : #include <algorithm>
25 : #include "tools/NeighborList.h"
26 : #include "tools/Communicator.h"
27 : #include "tools/OpenMP.h"
28 :
29 : namespace PLMD {
30 : namespace colvar {
31 :
32 266 : void CoordinationBase::registerKeywords( Keywords& keys ) {
33 266 : Colvar::registerKeywords(keys);
34 266 : keys.addFlag("SERIAL",false,"Perform the calculation in serial - for debug purpose");
35 266 : keys.addFlag("PAIR",false,"Pair only 1st element of the 1st group with 1st element in the second, etc");
36 266 : keys.addFlag("NLIST",false,"Use a neighbor list to speed up the calculation");
37 266 : keys.addFlag("NLISTCELLS",false,"Use a neighbor list to speed up the calculation - use the cell list implementation instead of the classical one");
38 266 : keys.add("optional","NL_CUTOFF","The cutoff for the neighbor list");
39 266 : keys.add("optional","NL_STRIDE","The frequency with which we are updating the atoms in the neighbor list");
40 266 : keys.add("atoms","GROUPA","First list of atoms");
41 266 : keys.add("atoms","GROUPB","Second list of atoms (if empty, N*(N-1)/2 pairs in GROUPA are counted)");
42 266 : }
43 :
44 254 : CoordinationBase::CoordinationBase(const ActionOptions&ao):
45 : PLUMED_COLVAR_INIT(ao),
46 254 : pbc(true),
47 254 : serial(false),
48 254 : invalidateList(true),
49 254 : firsttime(true) {
50 :
51 254 : parseFlag("SERIAL",serial);
52 :
53 254 : std::vector<AtomNumber> ga_lista,gb_lista;
54 254 : parseAtomList("GROUPA",ga_lista);
55 254 : parseAtomList("GROUPB",gb_lista);
56 :
57 254 : bool nopbc=!pbc;
58 254 : parseFlag("NOPBC",nopbc);
59 254 : pbc=!nopbc;
60 :
61 : // pair stuff
62 254 : bool dopair=false;
63 254 : parseFlag("PAIR",dopair);
64 :
65 : // neighbor list stuff
66 254 : bool doneigh_classic=false;
67 254 : double nl_cut=0.0;
68 254 : int nl_st=0;
69 254 : parseFlag("NLIST",doneigh_classic);
70 254 : bool doneighcells=false;
71 254 : parseFlag("NLISTCELLS",doneighcells);
72 : //temporary message
73 254 : plumed_assert(!(doneighcells && doneigh_classic)) << "Please activate only one of the two version of the NL";
74 254 : plumed_assert(!(doneighcells && dopair)) << "Pair is not compatible with the CELLS implementation of the NL";
75 254 : bool doneigh=doneighcells||doneigh_classic;
76 : if(doneigh) {
77 45 : parse("NL_CUTOFF",nl_cut);
78 45 : if(nl_cut<=0.0) {
79 0 : error("NL_CUTOFF should be explicitly specified and positive");
80 : }
81 45 : parse("NL_STRIDE",nl_st);
82 45 : if(nl_st<=0) {
83 0 : error("NL_STRIDE should be explicitly specified and positive");
84 : }
85 : }
86 :
87 254 : addValueWithDerivatives();
88 254 : setNotPeriodic();
89 254 : if(gb_lista.size()>0) {
90 239 : if(doneigh) {
91 82 : nl=Tools::make_unique<NeighborList>(ga_lista,gb_lista,serial,dopair,pbc,getPbc(),comm,nl_cut,nl_st,doneighcells);
92 : } else {
93 396 : nl=Tools::make_unique<NeighborList>(ga_lista,gb_lista,serial,dopair,pbc,getPbc(),comm);
94 : }
95 : } else {
96 15 : if(doneigh) {
97 8 : nl=Tools::make_unique<NeighborList>(ga_lista,serial,pbc,getPbc(),comm,nl_cut,nl_st,doneighcells);
98 : } else {
99 22 : nl=Tools::make_unique<NeighborList>(ga_lista,serial,pbc,getPbc(),comm);
100 : }
101 : }
102 :
103 254 : requestAtoms(nl->getFullAtomList());
104 :
105 254 : log.printf(" between two groups of %u and %u atoms\n",static_cast<unsigned>(ga_lista.size()),static_cast<unsigned>(gb_lista.size()));
106 254 : log.printf(" first group:\n");
107 7317 : for(unsigned int i=0; i<ga_lista.size(); ++i) {
108 7063 : if ( (i+1) % 25 == 0 ) {
109 221 : log.printf(" \n");
110 : }
111 7063 : log.printf(" %d", ga_lista[i].serial());
112 : }
113 254 : log.printf(" \n second group:\n");
114 12061 : for(unsigned int i=0; i<gb_lista.size(); ++i) {
115 11807 : if ( (i+1) % 25 == 0 ) {
116 387 : log.printf(" \n");
117 : }
118 11807 : log.printf(" %d", gb_lista[i].serial());
119 : }
120 254 : log.printf(" \n");
121 254 : if(pbc) {
122 250 : log.printf(" using periodic boundary conditions\n");
123 : } else {
124 4 : log.printf(" without periodic boundary conditions\n");
125 : }
126 254 : if(dopair) {
127 6 : log.printf(" with PAIR option\n");
128 : }
129 254 : if(doneigh) {
130 45 : log.printf(" using neighbor lists with\n");
131 45 : log.printf(" update every %d steps and cutoff %f\n",nl_st,nl_cut);
132 : }
133 254 : }
134 :
135 254 : CoordinationBase::~CoordinationBase() {
136 : // destructor required to delete forward declared class
137 254 : }
138 :
139 4616 : void CoordinationBase::prepare() {
140 4616 : std::tie(firsttime,invalidateList) =nl->prepare(this,firsttime, invalidateList).get();
141 4616 : }
142 :
143 : // calculator
144 4391 : void CoordinationBase::calculate() {
145 4391 : if(pbc) {
146 4371 : calculateImpl<true>();
147 : } else {
148 20 : calculateImpl<false>();
149 : }
150 4391 : }
151 :
152 : template<bool useBatch>
153 4391 : void CoordinationBase::calculateImpl() {
154 :
155 4391 : double ncoord=0.;
156 4391 : Tensor virial;
157 4391 : std::vector<Vector> deriv(getNumberOfAtoms());
158 :
159 4391 : if(nl->getStride()>0 && invalidateList) {
160 251 : nl->update(getPositions());
161 : }
162 :
163 : unsigned stride;
164 : unsigned rank;
165 4391 : if(serial) {
166 : stride=1;
167 : rank=0;
168 : } else {
169 4391 : stride=comm.Get_size();
170 4391 : rank=comm.Get_rank();
171 : }
172 :
173 4391 : unsigned nt=OpenMP::getNumThreads();
174 4391 : const unsigned nn=nl->size();
175 4391 : if(nt*stride*10>nn) {
176 : nt=1;
177 : }
178 :
179 4391 : const unsigned elementsPerRank = std::ceil(double(nn)/stride);
180 4391 : const unsigned int start= rank*elementsPerRank;
181 4391 : const unsigned int end = ((start + elementsPerRank)< nn)?(start + elementsPerRank): nn;
182 :
183 4391 : #pragma omp parallel num_threads(nt)
184 : {
185 : // the per-thread buffers are only needed when several threads accumulate (otherwise deriv and virial are used directly)
186 : std::vector<Vector> omp_deriv(nt>1 ? getPositions().size() : 0);
187 : Tensor omp_virial;
188 :
189 : if constexpr (useBatch) {
190 : // Keep scratch storage bounded and leave enough chunks for all threads.
191 : constexpr unsigned batchSize = 1024;
192 : const unsigned work=end>start ? end-start : 0;
193 : const unsigned batchStep=nt>1 ? std::min(batchSize,std::max(1u,work/(4*nt))) : batchSize;
194 : std::array<Vector, batchSize> distances;
195 : std::array<unsigned, batchSize> firstAtoms, secondAtoms;
196 : #pragma omp for reduction(+:ncoord) nowait
197 : for(unsigned begin=start; begin<end; begin+=batchStep) {
198 : const unsigned limit=std::min(batchStep,end-begin);
199 : unsigned count=0;
200 : for(unsigned j=0; j<limit; ++j) {
201 : const auto [i0,i1]=nl->getUpdatedPair(begin+j);
202 : if(getAbsoluteIndex(i0)==getAbsoluteIndex(i1)) {
203 : continue;
204 : }
205 : firstAtoms[count]=i0;
206 : secondAtoms[count]=i1;
207 : distances[count]=delta(getPosition(i0),getPosition(i1));
208 : ++count;
209 : }
210 : if(count) {
211 : getPbc().apply(VectorView(distances[0].data(),count));
212 : }
213 : for(unsigned j=0; j<count; ++j) {
214 : const auto i0=firstAtoms[j],i1=secondAtoms[j];
215 : const auto& distance=distances[j];
216 : double dfunc=0.;
217 : ncoord += pairing(distance.modulo2(),dfunc,i0,i1);
218 : Vector dd(dfunc*distance);
219 : Tensor vv(dd,distance);
220 : if(nt>1) {
221 : omp_deriv[i0]-=dd;
222 : omp_deriv[i1]+=dd;
223 : omp_virial-=vv;
224 : } else {
225 : deriv[i0]-=dd;
226 : deriv[i1]+=dd;
227 : virial-=vv;
228 : }
229 : }
230 : }
231 : } else {
232 : #pragma omp for reduction(+:ncoord) nowait
233 : for(unsigned int i=start; i<end; ++i) {
234 :
235 : Vector distance;
236 : const auto [i0, i1] = nl->getUpdatedPair(i);
237 :
238 : if(getAbsoluteIndex(i0)==getAbsoluteIndex(i1)) {
239 : continue;
240 : }
241 :
242 : distance=delta(getPosition(i0),getPosition(i1));
243 :
244 : double dfunc=0.;
245 : ncoord += pairing(distance.modulo2(), dfunc,i0,i1);
246 :
247 : Vector dd(dfunc*distance);
248 : Tensor vv(dd,distance);
249 : if(nt>1) {
250 : omp_deriv[i0]-=dd;
251 : omp_deriv[i1]+=dd;
252 : omp_virial-=vv;
253 : } else {
254 : deriv[i0]-=dd;
255 : deriv[i1]+=dd;
256 : virial-=vv;
257 : }
258 :
259 : }
260 : }
261 : #pragma omp critical
262 : if(nt>1) {
263 : for(unsigned i=0; i<getPositions().size(); i++) {
264 : deriv[i]+=omp_deriv[i];
265 : }
266 : virial+=omp_virial;
267 : }
268 : }
269 :
270 4391 : if(!serial) {
271 4391 : comm.Sum(ncoord);
272 4391 : if(!deriv.empty()) {
273 4391 : comm.Sum(&deriv[0][0],3*deriv.size());
274 : }
275 4391 : comm.Sum(virial);
276 : }
277 :
278 414197 : for(unsigned i=0; i<deriv.size(); ++i) {
279 409806 : setAtomsDerivatives(i,deriv[i]);
280 : }
281 4391 : setValue (ncoord);
282 4391 : setBoxDerivatives (virial);
283 :
284 4391 : }
285 : }
286 : }
|
