Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2026 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 "Colvar.h"
23 : #include "core/ActionRegister.h"
24 : #include "tools/Communicator.h"
25 : #include "tools/NeighborList.h"
26 : #include "tools/OpenMP.h"
27 :
28 : #include <algorithm>
29 : #include <cctype>
30 : #include <cmath>
31 : #include <iomanip>
32 : #include <limits>
33 : #include <memory>
34 : #include <set>
35 : #include <sstream>
36 : #include <string>
37 : #include <type_traits>
38 : #include <utility>
39 : #include <vector>
40 :
41 : namespace PLMD {
42 : namespace colvar {
43 :
44 : namespace {
45 :
46 : template<typename KeywordsType>
47 45 : auto setScalarDescriptionImpl(KeywordsType& keys,
48 : const std::string& description, int)
49 : -> decltype(std::declval<KeywordsType&>().setValueDescription(
50 : std::declval<const std::string&>(),
51 : std::declval<const std::string&>())) {
52 45 : keys.setValueDescription("scalar",description);
53 45 : }
54 :
55 : template<typename KeywordsType>
56 : auto setScalarDescriptionImpl(KeywordsType& keys,
57 : const std::string& description, long)
58 : -> decltype(std::declval<KeywordsType&>().setValueDescription(
59 : std::declval<const std::string&>())) {
60 : keys.setValueDescription(description);
61 : }
62 :
63 : template<typename KeywordsType>
64 : void setScalarDescriptionImpl(KeywordsType&, const std::string&, ...) {}
65 :
66 : template<typename NeighborListType>
67 2 : std::unique_ptr<NeighborListType> makeCandidateNeighborListImpl(
68 : const std::vector<AtomNumber>& centers,
69 : const std::vector<AtomNumber>& assigned,
70 : const bool serial, const bool pbc, const Pbc& pbcObject,
71 : Communicator& communicator, const double cutoff, const unsigned stride,
72 : const bool requestCellList, bool& usingCellList, std::true_type) {
73 2 : usingCellList=requestCellList;
74 : return Tools::make_unique<NeighborListType>(
75 2 : centers,assigned,serial,false,pbc,pbcObject,communicator,
76 2 : cutoff,stride,requestCellList);
77 : }
78 :
79 : template<typename NeighborListType>
80 : std::unique_ptr<NeighborListType> makeCandidateNeighborListImpl(
81 : const std::vector<AtomNumber>& centers,
82 : const std::vector<AtomNumber>& assigned,
83 : const bool serial, const bool pbc, const Pbc& pbcObject,
84 : Communicator& communicator, const double cutoff, const unsigned stride,
85 : const bool requestCellList, bool& usingCellList, std::false_type) {
86 : (void) requestCellList;
87 : usingCellList=false;
88 : return Tools::make_unique<NeighborListType>(
89 : centers,assigned,serial,false,pbc,pbcObject,communicator,
90 : cutoff,stride);
91 : }
92 :
93 : template<typename NeighborListType>
94 2 : std::unique_ptr<NeighborListType> makeCandidateNeighborList(
95 : const std::vector<AtomNumber>& centers,
96 : const std::vector<AtomNumber>& assigned,
97 : const bool serial, const bool pbc, const Pbc& pbcObject,
98 : Communicator& communicator, const double cutoff, const unsigned stride,
99 : const bool requestCellList, bool& usingCellList) {
100 : typedef typename std::is_constructible<
101 : NeighborListType,
102 : const std::vector<AtomNumber>&, const std::vector<AtomNumber>&,
103 : bool, bool, bool, const Pbc&, Communicator&, double, unsigned, bool
104 : >::type CellListConstructor;
105 : return makeCandidateNeighborListImpl<NeighborListType>(
106 : centers,assigned,serial,pbc,pbcObject,communicator,cutoff,stride,
107 2 : requestCellList,usingCellList,CellListConstructor());
108 : }
109 :
110 : }
111 :
112 : class SoftVoronoiBase : public Colvar {
113 : protected:
114 : enum InvalidPairKind {
115 : InvalidPairNone=0,
116 : InvalidPairMapping,
117 : InvalidCenterCoordinate,
118 : InvalidAssignedCoordinate,
119 : InvalidPbcBoxNonFinite,
120 : InvalidPbcBoxDegenerate,
121 : InvalidDisplacement,
122 : InvalidLength,
123 : InvalidSeparation,
124 : InvalidScore
125 : };
126 :
127 : struct InvalidPairRecord {
128 : InvalidPairKind kind;
129 : unsigned pairIndex;
130 : unsigned firstNeighborIndex;
131 : unsigned secondNeighborIndex;
132 : unsigned centerAtomSerial;
133 : unsigned assignedAtomSerial;
134 : unsigned rank;
135 : unsigned thread;
136 : bool pbc;
137 : Vector centerPosition;
138 : Vector assignedPosition;
139 : Vector displacement;
140 : double length;
141 : Tensor box;
142 :
143 130 : InvalidPairRecord():
144 130 : kind(InvalidPairNone),
145 130 : pairIndex(std::numeric_limits<unsigned>::max()),
146 130 : firstNeighborIndex(std::numeric_limits<unsigned>::max()),
147 130 : secondNeighborIndex(std::numeric_limits<unsigned>::max()),
148 130 : centerAtomSerial(0),
149 130 : assignedAtomSerial(0),
150 130 : rank(0),
151 130 : thread(0),
152 130 : pbc(false),
153 : centerPosition(std::numeric_limits<double>::quiet_NaN(),
154 : std::numeric_limits<double>::quiet_NaN(),
155 : std::numeric_limits<double>::quiet_NaN()),
156 : assignedPosition(std::numeric_limits<double>::quiet_NaN(),
157 : std::numeric_limits<double>::quiet_NaN(),
158 : std::numeric_limits<double>::quiet_NaN()),
159 : displacement(std::numeric_limits<double>::quiet_NaN(),
160 : std::numeric_limits<double>::quiet_NaN(),
161 : std::numeric_limits<double>::quiet_NaN()),
162 130 : length(std::numeric_limits<double>::quiet_NaN()) {}
163 : };
164 :
165 : struct PairData {
166 : unsigned center;
167 : unsigned assigned;
168 : Vector distance;
169 : double length;
170 : double score;
171 : double weight;
172 : };
173 :
174 49 : struct Assignment {
175 : std::vector<PairData> pairs;
176 : std::vector<double> occupancy;
177 : };
178 :
179 : bool pbc_;
180 : bool serial_;
181 : bool invalidateList_;
182 : bool firstTime_;
183 : bool filterCandidatePairs_;
184 : bool useCellList_;
185 : bool haveNeighborReference_;
186 : double kappa_;
187 : double neighborCutoff_;
188 : double neighborSkin_;
189 : std::unique_ptr<NeighborList> neighborList_;
190 : std::vector<AtomNumber> centers_;
191 : std::vector<AtomNumber> assigned_;
192 : std::vector<double> reference_;
193 : std::vector<Vector> neighborReferencePositions_;
194 : Tensor neighborReferenceBox_;
195 :
196 : static void registerCommonKeywords(Keywords&);
197 : static void setScalarDescription(Keywords&, const std::string&);
198 : static bool finiteVector(const Vector&);
199 : static bool finiteTensor(const Tensor&);
200 : static const char* invalidPairKindName(InvalidPairKind);
201 : std::string invalidPairMessage(const InvalidPairRecord&) const;
202 : void broadcastOrCheck(std::vector<double>&, unsigned, const std::string&) const;
203 : std::string normalizedSign(std::string) const;
204 : std::vector<unsigned> mapSelection(const std::vector<AtomNumber>&,
205 : const std::string&, bool) const;
206 : bool neighborListNeedsUpdate();
207 : void saveNeighborListReference();
208 : Assignment calculateAssignment();
209 : std::vector<double> defects(const Assignment&) const;
210 : void addAssignmentDerivatives(const Assignment&, const std::vector<double>&,
211 : std::vector<Vector>&, Tensor&);
212 : void finishSetup(const std::string&);
213 : void finalize(double, std::vector<Vector>&, Tensor&);
214 : bool ownsDirectDerivatives() const;
215 :
216 : public:
217 : explicit SoftVoronoiBase(const ActionOptions&);
218 : void prepare() override;
219 : };
220 :
221 45 : void SoftVoronoiBase::registerCommonKeywords(Keywords& keys) {
222 45 : keys.add("atoms","CENTERS","Atoms that receive the smooth assignment");
223 45 : keys.add("atoms","ASSIGNED","Atoms that are distributed over CENTERS");
224 45 : keys.add("compulsory","KAPPA","Positive soft-assignment sharpness in inverse PLUMED length units");
225 45 : keys.add("compulsory","REFERENCE",
226 : "One intended occupancy broadcast to all CENTERS, or one value per CENTER in CENTERS order");
227 45 : keys.addFlag("SERIAL",false,"Perform the calculation redundantly on each rank for debugging");
228 45 : keys.addFlag("NLIST",false,"Use an approximate neighbor-list truncation of the assignment candidates");
229 45 : keys.add("optional","NL_CUTOFF","Candidate cutoff in PLUMED length units; every ASSIGNED atom must retain at least one CENTER");
230 45 : keys.add("optional","NL_STRIDE","Number of steps between neighbor-list updates");
231 45 : keys.add("optional","NL_SKIN","Nonnegative Verlet buffer in PLUMED length units; requires NL_STRIDE greater than one");
232 45 : }
233 :
234 45 : void SoftVoronoiBase::setScalarDescription(
235 : Keywords& keys, const std::string& description) {
236 45 : setScalarDescriptionImpl(keys,description,0);
237 45 : }
238 :
239 1925 : bool SoftVoronoiBase::finiteVector(const Vector& value) {
240 7691 : for(unsigned i=0; i<3; ++i) {
241 5769 : if(!std::isfinite(value[i])) {
242 : return false;
243 : }
244 : }
245 : return true;
246 : }
247 :
248 14 : bool SoftVoronoiBase::finiteTensor(const Tensor& value) {
249 56 : for(unsigned i=0; i<3; ++i) {
250 168 : for(unsigned j=0; j<3; ++j) {
251 126 : if(!std::isfinite(value(i,j))) {
252 : return false;
253 : }
254 : }
255 : }
256 : return true;
257 : }
258 :
259 14 : const char* SoftVoronoiBase::invalidPairKindName(
260 : const InvalidPairKind kind) {
261 14 : switch(kind) {
262 : case InvalidPairMapping:
263 : return "neighbor-list-mapping";
264 1 : case InvalidCenterCoordinate:
265 1 : return "center-coordinate-non-finite";
266 1 : case InvalidAssignedCoordinate:
267 1 : return "assigned-coordinate-non-finite";
268 0 : case InvalidPbcBoxNonFinite:
269 0 : return "pbc-box-non-finite";
270 1 : case InvalidPbcBoxDegenerate:
271 1 : return "pbc-box-degenerate";
272 1 : case InvalidDisplacement:
273 1 : return "minimum-image-displacement-non-finite";
274 1 : case InvalidLength:
275 1 : return "distance-non-finite";
276 8 : case InvalidSeparation:
277 8 : return "distance-zero-or-near-zero";
278 1 : case InvalidScore:
279 1 : return "kappa-times-distance-non-finite";
280 0 : case InvalidPairNone:
281 0 : return "none";
282 : }
283 0 : return "unknown";
284 : }
285 :
286 14 : std::string SoftVoronoiBase::invalidPairMessage(
287 : const InvalidPairRecord& record) const {
288 14 : std::ostringstream message;
289 : message << std::setprecision(17);
290 14 : if(record.kind==InvalidPairMapping) {
291 0 : message << "internal neighbor-list index does not map to CENTERS and ASSIGNED; ";
292 14 : } else if(record.kind==InvalidScore) {
293 1 : message << "KAPPA times a CENTER-ASSIGNED distance is too large; ";
294 : } else {
295 13 : message << "a CENTER-ASSIGNED distance is zero or non-finite; ";
296 : }
297 : message << "reactive Voronoi pair validation failed: category="
298 14 : << invalidPairKindName(record.kind)
299 14 : << " step=" << getStep()
300 14 : << " action=" << getLabel()
301 28 : << " pair_index=" << record.pairIndex << "(zero-based)"
302 14 : << " center_atom=";
303 14 : if(record.centerAtomSerial>0) {
304 : message << record.centerAtomSerial;
305 : } else {
306 0 : message << "unavailable";
307 : }
308 14 : message << " assigned_atom=";
309 14 : if(record.assignedAtomSerial>0) {
310 : message << record.assignedAtomSerial;
311 : } else {
312 0 : message << "unavailable";
313 : }
314 14 : message << " neighbor_indices=(" << record.firstNeighborIndex
315 14 : << "," << record.secondNeighborIndex << ")"
316 14 : << " rank=" << record.rank
317 14 : << " thread=" << record.thread
318 : << " coordinate_units=PLUMED-internal"
319 14 : << " center_position=(" << record.centerPosition[0]
320 14 : << "," << record.centerPosition[1]
321 14 : << "," << record.centerPosition[2] << ")"
322 14 : << " assigned_position=(" << record.assignedPosition[0]
323 14 : << "," << record.assignedPosition[1]
324 14 : << "," << record.assignedPosition[2] << ")"
325 14 : << " minimum_image_displacement=(" << record.displacement[0]
326 14 : << "," << record.displacement[1]
327 14 : << "," << record.displacement[2] << ")"
328 14 : << " length=" << record.length
329 14 : << " threshold=" << std::numeric_limits<double>::epsilon()
330 14 : << " pbc=" << (record.pbc ? "on" : "off")
331 26 : << " box_determinant=" << record.box.determinant()
332 14 : << " box=[";
333 56 : for(unsigned i=0; i<3; ++i) {
334 42 : if(i>0) {
335 28 : message << ";";
336 : }
337 168 : for(unsigned j=0; j<3; ++j) {
338 126 : if(j>0) {
339 84 : message << ",";
340 : }
341 126 : message << record.box(i,j);
342 : }
343 : }
344 14 : message << "]";
345 14 : return message.str();
346 14 : }
347 :
348 40 : void SoftVoronoiBase::broadcastOrCheck(std::vector<double>& keywordValues,
349 : const unsigned size,
350 : const std::string& keyword) const {
351 40 : if(keywordValues.size()==1 && size>1) {
352 17 : const double value=keywordValues[0];
353 17 : keywordValues.assign(size,value);
354 : }
355 40 : if(keywordValues.size()!=size) {
356 0 : error(keyword+" must contain one value or exactly "+
357 0 : std::to_string(size)+" values");
358 : }
359 178 : for(unsigned i=0; i<keywordValues.size(); ++i) {
360 138 : if(!std::isfinite(keywordValues[i])) {
361 0 : error(keyword+" contains a non-finite value at position "+
362 0 : std::to_string(i+1));
363 : }
364 : }
365 40 : }
366 :
367 31 : std::string SoftVoronoiBase::normalizedSign(std::string sign) const {
368 31 : std::transform(sign.begin(),sign.end(),sign.begin(),
369 : [](const char value) {
370 118 : return static_cast<char>(std::toupper(static_cast<unsigned char>(value)));
371 : });
372 31 : if(sign!="ALL" && sign!="POSITIVE" && sign!="NEGATIVE") {
373 0 : error("SIGN must be ALL, POSITIVE, or NEGATIVE");
374 : }
375 31 : return sign;
376 : }
377 :
378 33 : SoftVoronoiBase::SoftVoronoiBase(const ActionOptions& ao):
379 : PLUMED_COLVAR_INIT(ao),
380 33 : pbc_(true),
381 33 : serial_(false),
382 33 : invalidateList_(true),
383 33 : firstTime_(true),
384 33 : filterCandidatePairs_(false),
385 33 : useCellList_(false),
386 33 : haveNeighborReference_(false),
387 33 : kappa_(0.0),
388 33 : neighborCutoff_(0.0),
389 33 : neighborSkin_(0.0) {
390 :
391 33 : parseAtomList("CENTERS",centers_);
392 66 : parseAtomList("ASSIGNED",assigned_);
393 33 : if(centers_.empty()) {
394 0 : error("CENTERS must contain at least one atom");
395 : }
396 33 : if(assigned_.empty()) {
397 0 : error("ASSIGNED must contain at least one atom");
398 : }
399 :
400 : std::set<unsigned> centerIds;
401 158 : for(unsigned i=0; i<centers_.size(); ++i) {
402 125 : if(!centerIds.insert(centers_[i].index()).second) {
403 0 : error("CENTERS contains duplicate atoms");
404 : }
405 : }
406 : std::set<unsigned> assignedIds;
407 155 : for(unsigned i=0; i<assigned_.size(); ++i) {
408 122 : if(!assignedIds.insert(assigned_[i].index()).second) {
409 0 : error("ASSIGNED contains duplicate atoms");
410 : }
411 244 : if(centerIds.count(assigned_[i].index())>0) {
412 0 : error("CENTERS and ASSIGNED must be disjoint");
413 : }
414 : }
415 :
416 33 : parse("KAPPA",kappa_);
417 33 : if(!std::isfinite(kappa_) || kappa_<=0.0) {
418 0 : error("KAPPA must be finite and positive");
419 : }
420 33 : parseVector("REFERENCE",reference_);
421 33 : broadcastOrCheck(reference_,centers_.size(),"REFERENCE");
422 :
423 33 : parseFlag("SERIAL",serial_);
424 33 : bool noPbc=false;
425 33 : parseFlag("NOPBC",noPbc);
426 33 : pbc_=!noPbc;
427 :
428 33 : bool useNeighborList=false;
429 33 : double neighborCutoff=0.0;
430 33 : double neighborSkin=0.0;
431 33 : int neighborStride=0;
432 33 : parseFlag("NLIST",useNeighborList);
433 33 : if(useNeighborList) {
434 2 : parse("NL_CUTOFF",neighborCutoff);
435 2 : parse("NL_STRIDE",neighborStride);
436 2 : parse("NL_SKIN",neighborSkin);
437 2 : if(!std::isfinite(neighborCutoff) || neighborCutoff<=0.0) {
438 0 : error("NL_CUTOFF must be finite and positive");
439 : }
440 2 : if(neighborStride<=0) {
441 0 : error("NL_STRIDE must be positive");
442 : }
443 2 : if(!std::isfinite(neighborSkin) || neighborSkin<0.0) {
444 0 : error("NL_SKIN must be finite and nonnegative");
445 : }
446 2 : if(neighborSkin>0.0 && neighborStride==1) {
447 0 : error("NL_SKIN requires NL_STRIDE greater than one");
448 : }
449 : }
450 :
451 33 : if(useNeighborList) {
452 2 : neighborCutoff_=neighborCutoff;
453 2 : neighborSkin_=neighborSkin;
454 2 : filterCandidatePairs_=neighborSkin_>0.0 || neighborStride>1;
455 2 : const double candidateCutoff=neighborCutoff_+neighborSkin_;
456 2 : if(!std::isfinite(candidateCutoff)) {
457 0 : error("NL_CUTOFF plus NL_SKIN must be finite");
458 : }
459 4 : const unsigned neighborThreads=serial_ ? 1 :
460 4 : std::max(1U,OpenMP::getNumThreads());
461 : const unsigned long long fullPairCount=
462 2 : static_cast<unsigned long long>(centers_.size())*assigned_.size();
463 2 : const bool requestCellList=pbc_ &&
464 2 : fullPairCount>=32768ULL*neighborThreads;
465 2 : neighborList_=makeCandidateNeighborList<NeighborList>(
466 2 : centers_,assigned_,serial_,pbc_,getPbc(),comm,
467 : candidateCutoff,neighborStride,requestCellList,
468 2 : useCellList_);
469 2 : filterCandidatePairs_=filterCandidatePairs_ || useCellList_;
470 : } else {
471 62 : neighborList_=Tools::make_unique<NeighborList>(
472 62 : centers_,assigned_,serial_,false,pbc_,getPbc(),comm);
473 : }
474 33 : }
475 :
476 34 : std::vector<unsigned> SoftVoronoiBase::mapSelection(
477 : const std::vector<AtomNumber>& atoms, const std::string& keyword,
478 : const bool defaultAll) const {
479 34 : std::vector<unsigned> mapped;
480 34 : if(atoms.empty()) {
481 28 : if(!defaultAll) {
482 0 : error(keyword+" must contain at least one atom");
483 : }
484 28 : mapped.resize(centers_.size());
485 143 : for(unsigned i=0; i<centers_.size(); ++i) {
486 115 : mapped[i]=i;
487 : }
488 : return mapped;
489 : }
490 :
491 : std::set<unsigned> ids;
492 14 : for(unsigned i=0; i<atoms.size(); ++i) {
493 8 : if(!ids.insert(atoms[i].index()).second) {
494 0 : error(keyword+" contains duplicate atoms");
495 : }
496 : const std::vector<AtomNumber>::const_iterator found=
497 8 : std::find(centers_.begin(),centers_.end(),atoms[i]);
498 8 : if(found==centers_.end()) {
499 0 : error("every atom in "+keyword+" must also be present in CENTERS");
500 : }
501 8 : mapped.push_back(static_cast<unsigned>(found-centers_.begin()));
502 : }
503 : return mapped;
504 0 : }
505 :
506 33 : void SoftVoronoiBase::finishSetup(const std::string& displayName) {
507 33 : checkRead();
508 33 : addValueWithDerivatives();
509 33 : setNotPeriodic();
510 : // PLUMED 2.8 needs the full atom list to preserve absolute-index mapping.
511 33 : requestAtoms(neighborList_->getFullAtomList());
512 33 : log.printf(" %s assigns %u atoms over %u centers with kappa %g\n",
513 : displayName.c_str(),static_cast<unsigned>(assigned_.size()),
514 : static_cast<unsigned>(centers_.size()),kappa_);
515 33 : if(neighborList_->getStride()>0) {
516 2 : log.printf(" neighbor cutoff %g skin %g stride %u using %s\n",
517 : neighborCutoff_,neighborSkin_,neighborList_->getStride(),
518 2 : useCellList_ ? "link cells" : "a pair scan");
519 : }
520 33 : }
521 :
522 49 : void SoftVoronoiBase::prepare() {
523 49 : if(neighborList_->getStride()>0) {
524 4 : if(firstTime_ || getStep()%neighborList_->getStride()==0) {
525 3 : requestAtoms(neighborList_->getFullAtomList());
526 3 : invalidateList_=true;
527 3 : firstTime_=false;
528 : } else {
529 1 : requestAtoms(neighborList_->getFullAtomList());
530 1 : invalidateList_=false;
531 1 : if(getExchangeStep()) {
532 0 : error("neighbor lists must be updated on exchange steps; choose an NL_STRIDE that divides the exchange stride");
533 : }
534 : }
535 4 : if(getExchangeStep()) {
536 0 : firstTime_=true;
537 : }
538 : }
539 49 : }
540 :
541 1 : bool SoftVoronoiBase::neighborListNeedsUpdate() {
542 1 : if(neighborSkin_<=0.0 || !haveNeighborReference_) {
543 : return false;
544 : }
545 : const std::vector<Vector>& positions=getPositions();
546 1 : if(positions.size()!=neighborReferencePositions_.size()) {
547 : return true;
548 : }
549 1 : if(pbc_) {
550 1 : const Tensor box=getPbc().getBox();
551 4 : for(unsigned i=0; i<3; ++i) {
552 12 : for(unsigned j=0; j<3; ++j) {
553 9 : if(box(i,j)!=neighborReferenceBox_(i,j)) {
554 0 : return true;
555 : }
556 : }
557 : }
558 : }
559 1 : const double thresholdSquared=0.25*neighborSkin_*neighborSkin_;
560 3 : for(unsigned i=0; i<positions.size(); ++i) {
561 3 : const Vector displacement=pbc_ ?
562 : pbcDistance(neighborReferencePositions_[i],
563 : positions[i]) :
564 : delta(neighborReferencePositions_[i],
565 : positions[i]);
566 3 : if(displacement.modulo2()>thresholdSquared) {
567 : return true;
568 : }
569 : }
570 : return false;
571 : }
572 :
573 4 : void SoftVoronoiBase::saveNeighborListReference() {
574 4 : if(neighborSkin_<=0.0) {
575 : return;
576 : }
577 2 : neighborReferencePositions_=getPositions();
578 2 : if(pbc_) {
579 2 : neighborReferenceBox_=getPbc().getBox();
580 : }
581 2 : haveNeighborReference_=true;
582 : }
583 :
584 49 : SoftVoronoiBase::Assignment SoftVoronoiBase::calculateAssignment() {
585 : // Refresh the buffered candidate list on schedule or after a safe-skin exit.
586 49 : if(neighborList_->getStride()>0) {
587 4 : const bool updateList=invalidateList_ || neighborListNeedsUpdate();
588 : if(updateList) {
589 4 : neighborList_->update(getPositions());
590 4 : saveNeighborListReference();
591 : }
592 : }
593 :
594 : // Split candidate pairs into contiguous MPI-owned ranges.
595 49 : const unsigned numberOfPairs=neighborList_->size();
596 49 : const unsigned stride=serial_ ? 1 : comm.Get_size();
597 49 : const unsigned rank=serial_ ? 0 : comm.Get_rank();
598 49 : const unsigned communicatorRank=comm.Get_rank();
599 49 : const unsigned pairsPerRank=(numberOfPairs+stride-1)/stride;
600 49 : const unsigned start=std::min(rank*pairsPerRank,numberOfPairs);
601 49 : const unsigned end=std::min(start+pairsPerRank,numberOfPairs);
602 49 : const unsigned localPairCount=end-start;
603 49 : const double negativeInfinity=-std::numeric_limits<double>::infinity();
604 :
605 49 : Assignment result;
606 : unsigned numberOfThreads=1;
607 : #ifdef _OPENMP
608 49 : numberOfThreads=OpenMP::getGoodNumThreads(
609 : static_cast<PairData*>(nullptr),localPairCount);
610 : #endif
611 : std::vector<double> maximumScore(
612 63 : numberOfThreads*assigned_.size(),negativeInfinity);
613 49 : std::vector<InvalidPairRecord> threadFailures(numberOfThreads);
614 49 : const Tensor currentBox=getPbc().getBox();
615 49 : const double currentBoxDeterminant=currentBox.determinant();
616 49 : InvalidPairKind boxFailure=InvalidPairNone;
617 49 : if(pbc_) {
618 14 : if(!finiteTensor(currentBox) ||
619 : !std::isfinite(currentBoxDeterminant)) {
620 0 : boxFailure=InvalidPbcBoxNonFinite;
621 14 : } else if(!getPbc().isSet() ||
622 : std::fabs(currentBoxDeterminant)<=1.0e-14) {
623 1 : boxFailure=InvalidPbcBoxDegenerate;
624 : }
625 : }
626 :
627 : // Evaluate retained pairs and collect per-thread shifted-softmax maxima.
628 643 : const auto evaluatePair=[&](const unsigned pairIndex,
629 : const unsigned thread, PairData& pair) {
630 : const std::pair<unsigned,unsigned> pairIndexes=
631 643 : neighborList_->getClosePair(pairIndex);
632 643 : bool mappingValid=false;
633 643 : bool positionsAvailable=false;
634 643 : bool displacementAvailable=false;
635 643 : bool lengthAvailable=false;
636 643 : Vector centerPosition;
637 643 : Vector assignedPosition;
638 11 : const auto recordFailure=[&](const InvalidPairKind kind) {
639 11 : InvalidPairRecord failure;
640 11 : failure.kind=kind;
641 11 : failure.pairIndex=pairIndex;
642 11 : failure.firstNeighborIndex=pairIndexes.first;
643 11 : failure.secondNeighborIndex=pairIndexes.second;
644 11 : failure.rank=communicatorRank;
645 11 : failure.thread=thread;
646 11 : failure.pbc=pbc_;
647 11 : failure.box=currentBox;
648 11 : if(mappingValid) {
649 11 : failure.centerAtomSerial=centers_[pair.center].serial();
650 11 : failure.assignedAtomSerial=assigned_[pair.assigned].serial();
651 : }
652 11 : if(positionsAvailable) {
653 11 : failure.centerPosition=centerPosition;
654 11 : failure.assignedPosition=assignedPosition;
655 : }
656 11 : if(displacementAvailable) {
657 8 : failure.displacement=pair.distance;
658 : }
659 11 : if(lengthAvailable) {
660 7 : failure.length=pair.length;
661 : }
662 11 : if(failure.pairIndex<threadFailures[thread].pairIndex) {
663 11 : threadFailures[thread]=failure;
664 : }
665 654 : };
666 643 : if(pairIndexes.first>=centers_.size() ||
667 643 : pairIndexes.second<centers_.size()) {
668 0 : recordFailure(InvalidPairMapping);
669 11 : return;
670 : }
671 643 : pair.center=pairIndexes.first;
672 643 : pair.assigned=pairIndexes.second-centers_.size();
673 643 : if(pair.assigned>=assigned_.size()) {
674 0 : recordFailure(InvalidPairMapping);
675 0 : return;
676 : }
677 643 : mappingValid=true;
678 643 : centerPosition=getPosition(pairIndexes.first);
679 643 : assignedPosition=getPosition(pairIndexes.second);
680 643 : positionsAvailable=true;
681 643 : if(!finiteVector(centerPosition)) {
682 1 : recordFailure(InvalidCenterCoordinate);
683 1 : return;
684 : }
685 642 : if(!finiteVector(assignedPosition)) {
686 1 : recordFailure(InvalidAssignedCoordinate);
687 1 : return;
688 : }
689 641 : if(boxFailure!=InvalidPairNone) {
690 1 : recordFailure(boxFailure);
691 1 : return;
692 : }
693 640 : pair.distance=pbc_ ?
694 : pbcDistance(centerPosition,assignedPosition) :
695 : delta(centerPosition,assignedPosition);
696 640 : displacementAvailable=true;
697 640 : if(!finiteVector(pair.distance)) {
698 1 : recordFailure(InvalidDisplacement);
699 1 : return;
700 : }
701 639 : pair.length=pair.distance.modulo();
702 639 : lengthAvailable=true;
703 639 : if(!std::isfinite(pair.length)) {
704 1 : recordFailure(InvalidLength);
705 1 : return;
706 : }
707 638 : if(pair.length<=std::numeric_limits<double>::epsilon()) {
708 5 : recordFailure(InvalidSeparation);
709 5 : return;
710 : }
711 633 : if(filterCandidatePairs_ && pair.length>neighborCutoff_) {
712 0 : pair.length=-1.0;
713 0 : return;
714 : }
715 633 : pair.score=-kappa_*pair.length;
716 633 : if(!std::isfinite(pair.score)) {
717 1 : recordFailure(InvalidScore);
718 1 : return;
719 : }
720 632 : pair.weight=0.0;
721 632 : const unsigned maximumIndex=thread*assigned_.size()+pair.assigned;
722 632 : maximumScore[maximumIndex]=
723 1234 : std::max(maximumScore[maximumIndex],pair.score);
724 49 : };
725 :
726 49 : if(numberOfThreads==1) {
727 44 : result.pairs.reserve(localPairCount);
728 175 : for(unsigned pairIndex=start; pairIndex<end; ++pairIndex) {
729 131 : PairData pair= {};
730 131 : evaluatePair(pairIndex,0,pair);
731 131 : result.pairs.push_back(pair);
732 : }
733 : } else {
734 5 : result.pairs.resize(localPairCount);
735 5 : #pragma omp parallel for num_threads(numberOfThreads)
736 : for(unsigned pairIndex=start; pairIndex<end; ++pairIndex) {
737 : evaluatePair(pairIndex,OpenMP::getThreadNum(),
738 : result.pairs[pairIndex-start]);
739 : }
740 : }
741 :
742 : // Reduce thread-local failures and maxima before the MPI maximum reduction.
743 49 : InvalidPairRecord localFailure;
744 105 : for(unsigned thread=0; thread<numberOfThreads; ++thread) {
745 56 : if(threadFailures[thread].pairIndex<localFailure.pairIndex) {
746 11 : localFailure=threadFailures[thread];
747 : }
748 56 : if(thread>0) {
749 119 : for(unsigned j=0; j<assigned_.size(); ++j) {
750 112 : maximumScore[j]=
751 224 : std::max(maximumScore[j],maximumScore[thread*assigned_.size()+j]);
752 : }
753 : }
754 : }
755 :
756 : // Reuse the existing maximum collective to select the globally first bad
757 : // pair. A zero key means success; larger keys correspond to smaller pair
758 : // indices, so MPI_MAX selects the same pair independently of rank layout.
759 63 : std::vector<double> rankMaximum(assigned_.size()+1,negativeInfinity);
760 201 : for(unsigned j=0; j<assigned_.size(); ++j) {
761 152 : rankMaximum[j]=maximumScore[j];
762 : }
763 49 : rankMaximum.back()=localFailure.kind==InvalidPairNone ? 0.0 :
764 11 : static_cast<double>(numberOfPairs-localFailure.pairIndex);
765 49 : if(!serial_ && comm.Get_size()>1) {
766 4 : comm.Max(&rankMaximum[0],static_cast<int>(rankMaximum.size()));
767 : }
768 201 : for(unsigned j=0; j<assigned_.size(); ++j) {
769 152 : maximumScore[j]=rankMaximum[j];
770 : }
771 :
772 49 : if(rankMaximum.back()>0.0) {
773 14 : const unsigned errorKey=static_cast<unsigned>(rankMaximum.back());
774 14 : const unsigned globalPairIndex=numberOfPairs-errorKey;
775 14 : const unsigned reportingRank=(!serial_ && comm.Get_size()>1) ?
776 : globalPairIndex/pairsPerRank : communicatorRank;
777 14 : InvalidPairRecord selectedFailure;
778 14 : if(communicatorRank==reportingRank) {
779 15 : for(unsigned thread=0; thread<numberOfThreads; ++thread) {
780 15 : if(threadFailures[thread].pairIndex==globalPairIndex) {
781 11 : selectedFailure=threadFailures[thread];
782 11 : break;
783 : }
784 : }
785 : }
786 :
787 14 : unsigned metadata[9]= {0,0,0,0,0,0,0,0,0};
788 14 : double details[19]= {0.0};
789 14 : if(communicatorRank==reportingRank &&
790 11 : selectedFailure.kind!=InvalidPairNone) {
791 11 : metadata[0]=static_cast<unsigned>(selectedFailure.kind);
792 11 : metadata[1]=selectedFailure.pairIndex;
793 11 : metadata[2]=selectedFailure.firstNeighborIndex;
794 11 : metadata[3]=selectedFailure.secondNeighborIndex;
795 11 : metadata[4]=selectedFailure.centerAtomSerial;
796 11 : metadata[5]=selectedFailure.assignedAtomSerial;
797 11 : metadata[6]=selectedFailure.rank;
798 11 : metadata[7]=selectedFailure.thread;
799 11 : metadata[8]=selectedFailure.pbc ? 1U : 0U;
800 44 : for(unsigned i=0; i<3; ++i) {
801 33 : details[i]=selectedFailure.centerPosition[i];
802 33 : details[3+i]=selectedFailure.assignedPosition[i];
803 33 : details[6+i]=selectedFailure.displacement[i];
804 : }
805 11 : details[9]=selectedFailure.length;
806 44 : for(unsigned i=0; i<3; ++i) {
807 132 : for(unsigned j=0; j<3; ++j) {
808 99 : details[10+3*i+j]=selectedFailure.box(i,j);
809 : }
810 : }
811 : }
812 14 : if(!serial_ && comm.Get_size()>1) {
813 4 : comm.Bcast(metadata,9,static_cast<int>(reportingRank));
814 4 : comm.Bcast(details,19,static_cast<int>(reportingRank));
815 : }
816 14 : selectedFailure.kind=static_cast<InvalidPairKind>(metadata[0]);
817 14 : selectedFailure.pairIndex=metadata[1];
818 14 : selectedFailure.firstNeighborIndex=metadata[2];
819 14 : selectedFailure.secondNeighborIndex=metadata[3];
820 14 : selectedFailure.centerAtomSerial=metadata[4];
821 14 : selectedFailure.assignedAtomSerial=metadata[5];
822 14 : selectedFailure.rank=metadata[6];
823 14 : selectedFailure.thread=metadata[7];
824 14 : selectedFailure.pbc=metadata[8]!=0;
825 56 : for(unsigned i=0; i<3; ++i) {
826 42 : selectedFailure.centerPosition[i]=details[i];
827 42 : selectedFailure.assignedPosition[i]=details[3+i];
828 42 : selectedFailure.displacement[i]=details[6+i];
829 : }
830 14 : selectedFailure.length=details[9];
831 56 : for(unsigned i=0; i<3; ++i) {
832 168 : for(unsigned j=0; j<3; ++j) {
833 126 : selectedFailure.box(i,j)=details[10+3*i+j];
834 : }
835 : }
836 14 : if(selectedFailure.kind==InvalidPairNone ||
837 : selectedFailure.pairIndex!=globalPairIndex) {
838 0 : error("internal error while selecting the first invalid reactive Voronoi pair");
839 : }
840 14 : error(invalidPairMessage(selectedFailure));
841 : }
842 :
843 35 : if(filterCandidatePairs_) {
844 2 : result.pairs.erase(
845 2 : std::remove_if(result.pairs.begin(),result.pairs.end(),
846 : [](const PairData& pair) {
847 2 : return pair.length<0.0;
848 : }),result.pairs.end());
849 : }
850 98 : for(unsigned j=0; j<assigned_.size(); ++j) {
851 63 : if(!std::isfinite(maximumScore[j])) {
852 0 : error("an ASSIGNED atom has no CENTER inside the candidate list");
853 : }
854 : }
855 :
856 : // Form stable exponential weights and normalize each assigned atom.
857 35 : std::vector<double> normalization(assigned_.size(),0.0);
858 35 : if(numberOfThreads==1) {
859 157 : for(unsigned p=0; p<result.pairs.size(); ++p) {
860 : PairData& pair=result.pairs[p];
861 122 : pair.weight=std::exp(pair.score-maximumScore[pair.assigned]);
862 122 : normalization[pair.assigned]+=pair.weight;
863 : }
864 : } else {
865 : std::vector<double> threadNormalization(
866 0 : numberOfThreads*assigned_.size(),0.0);
867 0 : #pragma omp parallel for num_threads(numberOfThreads)
868 : for(unsigned p=0; p<result.pairs.size(); ++p) {
869 : const unsigned thread=OpenMP::getThreadNum();
870 : PairData& pair=result.pairs[p];
871 : pair.weight=std::exp(pair.score-maximumScore[pair.assigned]);
872 : threadNormalization[thread*assigned_.size()+pair.assigned]+=
873 : pair.weight;
874 : }
875 0 : for(unsigned thread=0; thread<numberOfThreads; ++thread) {
876 0 : for(unsigned j=0; j<assigned_.size(); ++j) {
877 0 : normalization[j]+=
878 0 : threadNormalization[thread*assigned_.size()+j];
879 : }
880 : }
881 0 : }
882 35 : if(!serial_ && comm.Get_size()>1) {
883 0 : comm.Sum(&normalization[0],static_cast<int>(normalization.size()));
884 : }
885 98 : for(unsigned j=0; j<normalization.size(); ++j) {
886 63 : if(!std::isfinite(normalization[j]) || normalization[j]<=0.0) {
887 0 : error("soft-assignment normalization is zero or non-finite");
888 : }
889 : }
890 :
891 : // Reduce normalized weights into center occupancies.
892 35 : result.occupancy.assign(centers_.size(),0.0);
893 35 : if(numberOfThreads==1) {
894 157 : for(unsigned p=0; p<result.pairs.size(); ++p) {
895 : PairData& pair=result.pairs[p];
896 122 : pair.weight/=normalization[pair.assigned];
897 122 : result.occupancy[pair.center]+=pair.weight;
898 : }
899 : } else {
900 : std::vector<double> threadOccupancy(
901 0 : numberOfThreads*centers_.size(),0.0);
902 0 : #pragma omp parallel for num_threads(numberOfThreads)
903 : for(unsigned p=0; p<result.pairs.size(); ++p) {
904 : const unsigned thread=OpenMP::getThreadNum();
905 : PairData& pair=result.pairs[p];
906 : pair.weight/=normalization[pair.assigned];
907 : threadOccupancy[thread*centers_.size()+pair.center]+=
908 : pair.weight;
909 : }
910 0 : for(unsigned thread=0; thread<numberOfThreads; ++thread) {
911 0 : for(unsigned i=0; i<centers_.size(); ++i) {
912 0 : result.occupancy[i]+=
913 0 : threadOccupancy[thread*centers_.size()+i];
914 : }
915 : }
916 0 : }
917 35 : if(!serial_ && comm.Get_size()>1) {
918 0 : comm.Sum(&result.occupancy[0],
919 : static_cast<int>(result.occupancy.size()));
920 : }
921 35 : return result;
922 98 : }
923 :
924 35 : std::vector<double> SoftVoronoiBase::defects(
925 : const Assignment& assignment) const {
926 35 : std::vector<double> result(centers_.size(),0.0);
927 103 : for(unsigned i=0; i<centers_.size(); ++i) {
928 68 : result[i]=assignment.occupancy[i]-reference_[i];
929 : }
930 35 : return result;
931 : }
932 :
933 35 : void SoftVoronoiBase::addAssignmentDerivatives(
934 : const Assignment& assignment,
935 : const std::vector<double>& derivativeByDefect,
936 : std::vector<Vector>& derivatives, Tensor& virial) {
937 35 : if(derivativeByDefect.size()!=centers_.size()) {
938 0 : error("internal defect derivative has the wrong size");
939 : }
940 :
941 : unsigned numberOfThreads=1;
942 : #ifdef _OPENMP
943 35 : numberOfThreads=OpenMP::getGoodNumThreads(assignment.pairs);
944 : #endif
945 35 : std::vector<double> meanDerivative(assigned_.size(),0.0);
946 35 : if(numberOfThreads==1) {
947 157 : for(unsigned p=0; p<assignment.pairs.size(); ++p) {
948 : const PairData& pair=assignment.pairs[p];
949 122 : meanDerivative[pair.assigned]+=
950 122 : derivativeByDefect[pair.center]*pair.weight;
951 : }
952 : } else {
953 : std::vector<double> threadMeanDerivative(
954 0 : numberOfThreads*assigned_.size(),0.0);
955 0 : #pragma omp parallel for num_threads(numberOfThreads)
956 : for(unsigned p=0; p<assignment.pairs.size(); ++p) {
957 : const unsigned thread=OpenMP::getThreadNum();
958 : const PairData& pair=assignment.pairs[p];
959 : threadMeanDerivative[thread*assigned_.size()+pair.assigned]+=
960 : derivativeByDefect[pair.center]*pair.weight;
961 : }
962 0 : for(unsigned thread=0; thread<numberOfThreads; ++thread) {
963 0 : for(unsigned j=0; j<assigned_.size(); ++j) {
964 0 : meanDerivative[j]+=
965 0 : threadMeanDerivative[thread*assigned_.size()+j];
966 : }
967 : }
968 0 : }
969 35 : if(!serial_ && comm.Get_size()>1) {
970 0 : comm.Sum(&meanDerivative[0],static_cast<int>(meanDerivative.size()));
971 : }
972 35 : if(numberOfThreads==1) {
973 157 : for(unsigned p=0; p<assignment.pairs.size(); ++p) {
974 : const PairData& pair=assignment.pairs[p];
975 : const double radialDerivative=
976 122 : -kappa_*pair.weight*
977 122 : (derivativeByDefect[pair.center]-meanDerivative[pair.assigned]);
978 : const Vector pairDerivative=
979 122 : (radialDerivative/pair.length)*pair.distance;
980 122 : const unsigned assignedIndex=pair.assigned+centers_.size();
981 : derivatives[pair.center]-=pairDerivative;
982 122 : derivatives[assignedIndex]+=pairDerivative;
983 244 : virial-=Tensor(pairDerivative,pair.distance);
984 : }
985 : } else {
986 : std::vector<std::vector<Vector>> threadDerivatives(
987 0 : numberOfThreads,std::vector<Vector>(derivatives.size()));
988 0 : std::vector<Tensor> threadVirials(numberOfThreads);
989 0 : #pragma omp parallel for num_threads(numberOfThreads)
990 : for(unsigned p=0; p<assignment.pairs.size(); ++p) {
991 : const unsigned thread=OpenMP::getThreadNum();
992 : const PairData& pair=assignment.pairs[p];
993 : const double radialDerivative=
994 : -kappa_*pair.weight*
995 : (derivativeByDefect[pair.center]-meanDerivative[pair.assigned]);
996 : const Vector pairDerivative=
997 : (radialDerivative/pair.length)*pair.distance;
998 : const unsigned assignedIndex=pair.assigned+centers_.size();
999 : threadDerivatives[thread][pair.center]-=pairDerivative;
1000 : threadDerivatives[thread][assignedIndex]+=pairDerivative;
1001 : threadVirials[thread]-=Tensor(pairDerivative,pair.distance);
1002 : }
1003 0 : for(unsigned thread=0; thread<numberOfThreads; ++thread) {
1004 0 : for(unsigned i=0; i<derivatives.size(); ++i) {
1005 0 : derivatives[i]+=threadDerivatives[thread][i];
1006 : }
1007 0 : virial+=threadVirials[thread];
1008 : }
1009 0 : }
1010 35 : }
1011 :
1012 12 : bool SoftVoronoiBase::ownsDirectDerivatives() const {
1013 12 : return serial_ || comm.Get_rank()==0;
1014 : }
1015 :
1016 35 : void SoftVoronoiBase::finalize(const double value,
1017 : std::vector<Vector>& derivatives,
1018 : Tensor& virial) {
1019 35 : if(!serial_ && comm.Get_size()>1) {
1020 0 : if(!derivatives.empty()) {
1021 0 : comm.Sum(&derivatives[0][0],
1022 : static_cast<int>(3*derivatives.size()));
1023 : }
1024 0 : comm.Sum(virial);
1025 : }
1026 :
1027 35 : if(!std::isfinite(value)) {
1028 0 : error("reactive Voronoi value is non-finite");
1029 : }
1030 35 : if(getPntrToValue()->getNumberOfDerivatives()!=3*derivatives.size()+9) {
1031 0 : error("internal derivative storage does not match the requested atom count");
1032 : }
1033 166 : for(unsigned i=0; i<derivatives.size(); ++i) {
1034 524 : for(unsigned k=0; k<3; ++k) {
1035 393 : if(!std::isfinite(derivatives[i][k])) {
1036 0 : error("reactive Voronoi coordinate derivative is non-finite");
1037 : }
1038 : }
1039 131 : setAtomsDerivatives(i,derivatives[i]);
1040 : }
1041 140 : for(unsigned i=0; i<3; ++i) {
1042 420 : for(unsigned j=0; j<3; ++j) {
1043 315 : if(!std::isfinite(virial(i,j))) {
1044 0 : error("reactive Voronoi box derivative is non-finite");
1045 : }
1046 : }
1047 : }
1048 35 : setBoxDerivatives(virial);
1049 35 : setValue(value);
1050 35 : }
1051 :
1052 : //+PLUMEDOC COLVAR VORONOI_COORDINATION
1053 : /*
1054 : Calculate a scalar reduction of smooth Voronoi coordination defects.
1055 :
1056 : Reactive processes such as proton transfer are difficult to describe with a
1057 : fixed molecular identity because the atom that carries the proton can change.
1058 : This Action assigns every atom in ASSIGNED continuously to the atoms in
1059 : CENTERS and reduces the resulting coordination defects to one scalar. The
1060 : same assignment is shared by [VORONOI_DISTANCE](VORONOI_DISTANCE.md) and
1061 : [VORONOI_POSITION](VORONOI_POSITION.md). Together the three Actions describe
1062 : the amount, separation, and location of coordination defects without
1063 : hard-coding water, glycine, a catalyst, or an atom-list position.
1064 :
1065 : The construction follows the descriptors introduced for acid-base equilibria
1066 : \cite Grifoni2019AcidBase and condensed-phase tautomerism
1067 : \cite Grifoni2020Tautomeric. Applications to solvated glycine, interfacial
1068 : water ions, electric-field effects, and electrocatalytic nitrogen reduction
1069 : are discussed in \cite Zhang2024Glycine, \cite Zhang2025Interfaces,
1070 : \cite Zhang2025ElectricField, and \cite Zhang2026NRR.
1071 :
1072 : An illustrated bilingual guide with downloadable teaching structures,
1073 : application examples, performance data, and additional references is
1074 : available on the
1075 : [Reactive Soft-Voronoi CV project page](https://zhang-pchao.github.io/code/reactive-voronoi/).
1076 :
1077 : VORONOI_COORDINATION, [VORONOI_DISTANCE](VORONOI_DISTANCE.md), and
1078 : [VORONOI_POSITION](VORONOI_POSITION.md) are part of PLUMED's default `colvar`
1079 : module and have no external library dependency. In a PLUMED installation that
1080 : contains these Actions, use them directly in the input; no [LOAD](LOAD.md)
1081 : line or optional module is required. Bias examples using
1082 : [OPES_METAD](OPES_METAD.md) additionally require the `opes` module.
1083 :
1084 : ## Soft assignment and coordination defects
1085 :
1086 : For a center $i$ and an assigned atom $j$, let $d_{ij}$ be their
1087 : minimum-image distance. The assignment weight is
1088 :
1089 : $$
1090 : w_{ij}=\frac{\exp(-\kappa d_{ij})}
1091 : {\sum_k\exp(-\kappa d_{kj})}.
1092 : $$
1093 :
1094 : By default, $d_{ij}$ is evaluated with the minimum-image convention.
1095 : NOPBC instead uses the direct coordinate difference, so the input coordinates
1096 : must already follow a consistent unwrapped image convention.
1097 :
1098 : The denominator contains all CENTERS for the same assigned atom, so
1099 : $\sum_i w_{ij}=1$. KAPPA is positive and has inverse units of the current
1100 : PLUMED length unit. Increasing KAPPA sharpens the assignment toward the
1101 : nearest center; reducing it spreads an assigned atom over more centers. The
1102 : implementation uses a shifted softmax, which improves numerical stability
1103 : without changing the mathematical value.
1104 :
1105 : For monitoring or structural diagnosis, a larger KAPPA can provide sharper
1106 : state labels once the relevant basins are already known. For biased sampling,
1107 : a smaller KAPPA usually gives smoother assignment and force changes when an
1108 : assigned atom switches between nearby centers, although a value that is too
1109 : small can blur distinct chemical basins. Scan KAPPA on representative
1110 : reactant, transition, product, and host-switching frames before applying a
1111 : bias. Values such as 5 for smoother sampling or 100 for sharper diagnosis are
1112 : application examples, not transferable defaults.
1113 :
1114 : The smooth occupancy and coordination defect of center $i$ are
1115 :
1116 : $$
1117 : n_i=\sum_j w_{ij}, \qquad q_i=n_i-\nu_i .
1118 : $$
1119 :
1120 : REFERENCE supplies $\nu_i$. A single value is broadcast to all CENTERS;
1121 : otherwise provide exactly one value per center in the same order as CENTERS.
1122 : For example, `REFERENCE=2` with two water O centers is the vector `(2,2)`, not
1123 : a one-center calculation. Fractional entries are valid when the model
1124 : deliberately shares one reference occupancy over symmetry-related centers;
1125 : they should not be introduced by averaging chemically nonequivalent sites.
1126 : The identities
1127 :
1128 : $$
1129 : \sum_i n_i=N_{\mathrm{assigned}}, \qquad
1130 : \sum_i q_i=N_{\mathrm{assigned}}-\sum_i\nu_i
1131 : $$
1132 :
1133 : provide useful checks on a new chemical mapping. The Action does not infer
1134 : elements, molecules, water, or a special reactive site from atom order.
1135 : CENTERS and ASSIGNED must be nonempty, internally unique, and disjoint.
1136 :
1137 : VORONOI_COORDINATION returns
1138 :
1139 : $$
1140 : Q_p=\sum_{i\in S} a_i q_i^p ,
1141 : $$
1142 :
1143 : where SELECT defines $S$ and COEFFICIENTS supplies $a_i$. SELECT
1144 : defaults to all CENTERS and COEFFICIENTS defaults to one. SIGN restricts the
1145 : sum to $q_i>0$ or $q_i<0$ when POSITIVE or NEGATIVE is selected.
1146 :
1147 : - POWER=1 preserves the signed defect. With sign filtering it is
1148 : non-differentiable exactly at $q_i=0$.
1149 : - POWER=2 measures defect activity without cancellation. With sign
1150 : filtering its value and first derivative are continuous at $q_i=0$,
1151 : but its second derivative has a cusp there.
1152 : - COEFFICIENTS can distinguish chemically different selected centers or
1153 : reproduce a published scalar. Coefficients do not change the assignment.
1154 :
1155 : Analytical coordinate and box derivatives are provided. The derivative with
1156 : respect to a defect is $a_i$ for POWER=1 and $2a_iq_i$ for POWER=2
1157 : inside the selected sign branch, and zero outside it. These derivatives are
1158 : propagated through every assignment weight. A CV intended for biasing should
1159 : avoid a POWER=1 sign boundary, a neighbor-list membership change, or another
1160 : non-smooth surface discussed below.
1161 :
1162 : The geometric defects are not formal electronic charges. Their physical
1163 : meaning comes from the chosen atom sets and reference occupancies and should
1164 : always be checked for the system of interest.
1165 :
1166 : ## Translating chemistry into keywords
1167 :
1168 : Build the input from chemistry rather than from atom-list positions:
1169 :
1170 : 1. Put atoms that can receive an assigned atom in CENTERS. For proton
1171 : transfer these are commonly O and N atoms.
1172 : 2. Put only the transferable atoms in ASSIGNED. Hydrogen atoms that cannot
1173 : participate in the process need not be included.
1174 : 3. Give every center its neutral or intended occupancy in REFERENCE. Typical
1175 : examples are 2 for a water O and a model-dependent value for a reactive O
1176 : or N.
1177 : 4. Start with the exact full-pair calculation. Inspect representative
1178 : neutral, product, transition, and multi-defect configurations before
1179 : choosing SELECT, POWER, SIGN, or COEFFICIENTS.
1180 : 5. Add [VORONOI_DISTANCE](VORONOI_DISTANCE.md) only when separation is needed,
1181 : and [VORONOI_POSITION](VORONOI_POSITION.md) only when a fixed spatial frame
1182 : is physically meaningful.
1183 :
1184 : REFERENCE follows CENTERS order. COEFFICIENTS follows SELECT order. If the
1185 : centers are reordered, reorder the corresponding numeric vector as well.
1186 : Atom-valued SELECT, GROUP1, and GROUP2 lists use absolute atom numbers and do
1187 : not rely on a center being last, first, or one of a fixed number of species.
1188 :
1189 : ## Exact and neighbor-list calculations
1190 :
1191 : Without NLIST, every CENTER-ASSIGNED pair is included and the finite-system
1192 : definition above is exact. NLIST truncates the candidate centers and
1193 : renormalizes the weights over the retained candidates, so it is an
1194 : approximation rather than an algebraically exact acceleration. Every
1195 : ASSIGNED atom must retain at least one CENTER or the calculation stops with an
1196 : error.
1197 :
1198 : Omitting NLIST therefore has no hidden cutoff. When NLIST is enabled,
1199 : NL_CUTOFF is an absolute CENTER-ASSIGNED distance in the active PLUMED length
1200 : units, not a bond cutoff or a water-specific constant. If $d_{\min}$ is
1201 : the nearest-center distance and $R$ is a trial cutoff, the screening
1202 : relation
1203 :
1204 : $$
1205 : \exp[-\kappa(R-d_{\min})] \le \epsilon
1206 : $$
1207 :
1208 : can provide an initial estimate for making a single omitted score small.
1209 : Smaller KAPPA generally requires a larger cutoff. This estimate does not
1210 : replace convergence because omitted scores, normalization errors, and force
1211 : derivatives can accumulate over many centers.
1212 :
1213 : NL_CUTOFF and NL_STRIDE are required whenever NLIST is present. NL_SKIN is
1214 : optional, defaults to zero, and can be used only when NL_STRIDE is greater
1215 : than one.
1216 :
1217 : Before using NLIST in production:
1218 :
1219 : 1. evaluate representative configurations with the exact full-pair form;
1220 : 2. increase NL_CUTOFF until values and forces agree within the required
1221 : tolerance;
1222 : 3. start with NL_STRIDE=1, then increase it only after checking consecutive
1223 : MD steps. With NL_STRIDE greater than one, NL_SKIN can add a Verlet buffer.
1224 : The list is rebuilt early if any requested atom moves by more than half the
1225 : skin or if the periodic box changes, while the evaluated pairs are still
1226 : filtered at the true NL_CUTOFF on every step.
1227 :
1228 : For replica exchange, every exchange step must also be a neighbor-list update
1229 : step. Choose the update schedule accordingly; normally NL_STRIDE should
1230 : divide the exchange stride.
1231 :
1232 : For a fixed box, the half-skin displacement check ensures that a pair cannot
1233 : cross NL_CUTOFF before it was present in the buffered candidate list. A
1234 : larger skin reduces rebuilds but retains more candidate pairs, so both the
1235 : skin and stride should be benchmarked. On PLUMED 2.11 and later, sufficiently
1236 : large periodic candidate lists use the built-in link-cell broad phase;
1237 : smaller lists retain the threaded pair scan because it is faster there.
1238 : Runtime plugins built against older PLUMED versions keep the compatible pair
1239 : scan, but can still amortize its rebuild cost with NL_SKIN.
1240 :
1241 : A cutoff copied from another system is not a convergence test. Pair-list
1242 : changes can introduce small discontinuities because the retained weights are
1243 : renormalized. Prefer exact mode for derivative validation and for small or
1244 : moderate systems. Use NLIST only after a value-and-force convergence scan
1245 : demonstrates a useful speed/accuracy tradeoff for the target system.
1246 :
1247 : ## CPU parallelism and practical scaling
1248 :
1249 : The exact calculation evaluates $N_{\mathrm{centers}}
1250 : N_{\mathrm{assigned}}$ pairs per step. For water with every O in CENTERS
1251 : and every H in ASSIGNED, this is $2N_{\mathrm{water}}^2$. Exact mode
1252 : therefore remains quadratic even when it is parallel: additional CPU workers
1253 : reduce elapsed time but do not change the asymptotic cost.
1254 :
1255 : The pair, normalization, occupancy, and derivative loops use OpenMP when it is
1256 : available. Set the PLUMED and OpenMP thread counts to the CPU cores allocated
1257 : to each molecular-dynamics rank, for example:
1258 :
1259 : ```bash
1260 : export PLUMED_NUM_THREADS=4
1261 : export OMP_NUM_THREADS=4
1262 : plumed driver --plumed plumed.dat --ixyz trajectory.xyz --box 3.0,3.0,3.0
1263 : ```
1264 :
1265 : Benchmark 1, 2, 4, and 8 threads on the same frames because small candidate
1266 : lists can spend more time entering parallel regions than doing pair work, and
1267 : large lists can become limited by memory bandwidth. Avoid CPU
1268 : oversubscription. In a GPU molecular-dynamics run these Actions still execute
1269 : on CPUs, so a GPU allocation alone does not accelerate them. MPI also
1270 : partitions pairs unless SERIAL is present; SERIAL is a debugging mode that
1271 : repeats the work on every rank, not a performance option.
1272 :
1273 : For large systems, a converged NLIST is the only option here that reduces the
1274 : number of retained assignment pairs. NL_CUTOFF controls truncation accuracy,
1275 : whereas NL_STRIDE controls how often membership is rebuilt. Start with
1276 : NL_STRIDE=1, converge NL_CUTOFF against exact values and derivatives, and only
1277 : then test a larger stride against the maximum atomic displacement between
1278 : updates. Report CV time per call as well as whole-simulation throughput so
1279 : that GPU force-model time and output time are not mistaken for CV cost.
1280 :
1281 : ## Worked example 1: water autoionization
1282 :
1283 : For water, oxygen atoms can be used as CENTERS, hydrogen atoms as ASSIGNED,
1284 : and the neutral reference occupancy is two. In a configuration containing
1285 : one hydronium and one hydroxide, the corresponding defects approach +1 and
1286 : -1.
1287 :
1288 : ```plumed
1289 : UNITS LENGTH=A
1290 : WaterO: GROUP ATOMS=1-4
1291 : WaterH: GROUP ATOMS=5-12
1292 :
1293 : ionization: VORONOI_COORDINATION CENTERS=WaterO ASSIGNED=WaterH KAPPA=5 REFERENCE=2 POWER=2
1294 : positive_amount: VORONOI_COORDINATION CENTERS=WaterO ASSIGNED=WaterH KAPPA=5 REFERENCE=2 POWER=2 SIGN=POSITIVE
1295 : negative_amount: VORONOI_COORDINATION CENTERS=WaterO ASSIGNED=WaterH KAPPA=5 REFERENCE=2 POWER=2 SIGN=NEGATIVE
1296 : positive_signed: VORONOI_COORDINATION CENTERS=WaterO ASSIGNED=WaterH KAPPA=5 REFERENCE=2 POWER=1 SIGN=POSITIVE
1297 : negative_signed: VORONOI_COORDINATION CENTERS=WaterO ASSIGNED=WaterH KAPPA=5 REFERENCE=2 POWER=1 SIGN=NEGATIVE
1298 :
1299 : PRINT ARG=ionization,positive_amount,negative_amount,positive_signed,negative_signed FILE=COLVAR
1300 : ```
1301 :
1302 : For an isolated ion pair, `ionization` approaches 2, the two squared branches
1303 : approach 1, and the signed branches approach +1 and -1. In a neutral frame
1304 : all five values approach zero. Soft values between these limits are expected
1305 : during proton transfer.
1306 :
1307 : ## Worked example 2: exact-to-NLIST acceleration
1308 :
1309 : The numeric cutoff, skin, and stride below are only input examples, not
1310 : transferable recommendations. Use the three stages in order: establish an
1311 : exact reference, converge NL_CUTOFF while rebuilding every step, and only
1312 : then test whether a displacement-safe skin and a longer stride improve
1313 : performance without changing values or forces.
1314 :
1315 : ```plumed
1316 : UNITS LENGTH=A
1317 : WaterO: GROUP ATOMS=1-4
1318 : WaterH: GROUP ATOMS=5-12
1319 :
1320 : # 1. Exact full-pair reference
1321 : exact: VORONOI_COORDINATION CENTERS=WaterO ASSIGNED=WaterH KAPPA=5 REFERENCE=2 POWER=2
1322 :
1323 : # 2. Converge NL_CUTOFF with a rebuild every step
1324 : trial: VORONOI_COORDINATION CENTERS=WaterO ASSIGNED=WaterH KAPPA=5 REFERENCE=2 POWER=2 NLIST NL_CUTOFF=8.0 NL_STRIDE=1
1325 :
1326 : # 3. Only after stage 2, amortize rebuilds with a tested skin and stride
1327 : fast: VORONOI_COORDINATION CENTERS=WaterO ASSIGNED=WaterH KAPPA=5 REFERENCE=2 POWER=2 NLIST NL_CUTOFF=8.0 NL_SKIN=1.0 NL_STRIDE=10
1328 :
1329 : PRINT ARG=exact,trial,fast FILE=COLVAR
1330 : DUMPDERIVATIVES ARG=exact,trial,fast FILE=DERIVATIVES STRIDE=1
1331 : ```
1332 :
1333 : ## Worked example 3: applying a bias
1334 :
1335 : First monitor the unbiased CV and check its scale. The following compact
1336 : input then biases the total ionization activity with OPES. PACE, BARRIER,
1337 : TEMP, and all production settings are system-dependent and must be justified
1338 : for the simulation being run.
1339 :
1340 : ```plumed
1341 : UNITS LENGTH=A ENERGY=kj/mol
1342 : WaterO: GROUP ATOMS=1-4
1343 : WaterH: GROUP ATOMS=5-12
1344 : ionization: VORONOI_COORDINATION CENTERS=WaterO ASSIGNED=WaterH KAPPA=5 REFERENCE=2 POWER=2
1345 : opes: OPES_METAD ARG=ionization PACE=500 BARRIER=40 TEMP=300
1346 : PRINT ARG=ionization,opes.bias FILE=COLVAR STRIDE=10
1347 : ```
1348 :
1349 : Passing a one-frame `plumed driver` test only establishes parsing and local
1350 : evaluation. Before production biasing, verify analytical derivatives,
1351 : energy/force units, restart behavior, the unbiased CV distribution, and a
1352 : short molecular-dynamics force-path run.
1353 :
1354 : ## Validation checklist
1355 :
1356 : For every new chemical system:
1357 :
1358 : 1. Check CENTERS, ASSIGNED, REFERENCE, SELECT, and COEFFICIENTS against a
1359 : labeled structure rather than a topology-order assumption.
1360 : 2. Confirm the occupancy and defect conservation identities on neutral and
1361 : reactive frames.
1362 : 3. Compare analytical and numerical derivatives away from sign cusps,
1363 : zero-distance configurations, and periodic branch cuts.
1364 : 4. Translate the whole system by a lattice vector and verify the same value.
1365 : 5. Reorder CENTERS and ASSIGNED, reorder numeric vectors consistently, and
1366 : verify the same value and forces.
1367 : 6. If NLIST is requested, converge both values and forces against exact mode.
1368 : 7. Compare serial and intended MPI/OpenMP execution on the same frames.
1369 : 8. Run a short, fixed-seed molecular-dynamics smoke test before enhanced
1370 : sampling, then inspect finite COLVAR, bias, force, and restart output.
1371 :
1372 : ## Troubleshooting and production cautions
1373 :
1374 : - `Action VORONOI_COORDINATION is not known`: use a PLUMED installation that
1375 : contains these Actions and verify that the simulation is loading the same
1376 : PLUMED kernel as `plumed driver`.
1377 : - `REFERENCE must contain ... values`: provide one value for broadcast or one
1378 : value per CENTER, in CENTERS order.
1379 : - `every atom in SELECT must also be present in CENTERS`: SELECT is a subset,
1380 : not an independent chemical group.
1381 : - `an ASSIGNED atom has no CENTER inside the candidate list`: NL_CUTOFF is too
1382 : small for that configuration or the atom sets are wrong. The calculation
1383 : intentionally stops instead of returning a biased normalization.
1384 : - Large changes after enabling NLIST indicate an unconverged truncation, not a
1385 : harmless implementation detail.
1386 : - Very large KAPPA makes the assignment nearly discrete and can produce sharp
1387 : force changes near equidistant centers. Choose it from physical and
1388 : numerical validation, not only from endpoint values.
1389 : - POWER=1 sign filtering, ABSOLUTE positions, periodic position branch cuts,
1390 : and neighbor-list membership changes are not globally smooth bias
1391 : coordinates. Choose a smooth path or constrain the sampled region when
1392 : these features are unavoidable.
1393 : */
1394 : //+ENDPLUMEDOC
1395 :
1396 : class VoronoiCoordination : public SoftVoronoiBase {
1397 : int power_;
1398 : std::string sign_;
1399 : std::vector<unsigned> selected_;
1400 : std::vector<double> coefficients_;
1401 :
1402 : public:
1403 : explicit VoronoiCoordination(const ActionOptions&);
1404 : static void registerKeywords(Keywords&);
1405 : void calculate() override;
1406 : };
1407 :
1408 : PLUMED_REGISTER_ACTION(VoronoiCoordination,"VORONOI_COORDINATION")
1409 :
1410 31 : void VoronoiCoordination::registerKeywords(Keywords& keys) {
1411 31 : Colvar::registerKeywords(keys);
1412 31 : registerCommonKeywords(keys);
1413 31 : keys.add("atoms","SELECT","Subset of CENTERS included in the scalar reduction; the default is all centers");
1414 31 : keys.add("optional","COEFFICIENTS","One coefficient, or one value per atom in SELECT");
1415 31 : keys.add("compulsory","POWER","1","Power of the selected occupancy defects; supported values are 1 and 2");
1416 31 : keys.add("compulsory","SIGN","ALL","Use ALL, POSITIVE, or NEGATIVE defects");
1417 31 : setScalarDescription(keys,"the selected reduction of the smooth occupancy defects");
1418 31 : }
1419 :
1420 27 : VoronoiCoordination::VoronoiCoordination(const ActionOptions& ao):
1421 : Action(ao),
1422 : SoftVoronoiBase(ao),
1423 27 : power_(1),
1424 27 : sign_("ALL") {
1425 27 : std::vector<AtomNumber> selectedAtoms;
1426 27 : parseAtomList("SELECT",selectedAtoms);
1427 27 : selected_=mapSelection(selectedAtoms,"SELECT",true);
1428 :
1429 54 : parseVector("COEFFICIENTS",coefficients_);
1430 27 : if(coefficients_.empty()) {
1431 20 : coefficients_.assign(selected_.size(),1.0);
1432 : } else {
1433 14 : broadcastOrCheck(coefficients_,selected_.size(),"COEFFICIENTS");
1434 : }
1435 :
1436 27 : parse("POWER",power_);
1437 27 : if(power_!=1 && power_!=2) {
1438 0 : error("POWER must be 1 or 2");
1439 : }
1440 27 : parse("SIGN",sign_);
1441 27 : sign_=normalizedSign(sign_);
1442 :
1443 27 : finishSetup("VORONOI_COORDINATION");
1444 27 : log.printf(" reducing %u selected centers with POWER=%d SIGN=%s\n",
1445 : static_cast<unsigned>(selected_.size()),power_,sign_.c_str());
1446 27 : }
1447 :
1448 37 : void VoronoiCoordination::calculate() {
1449 37 : const Assignment assignment=calculateAssignment();
1450 23 : const std::vector<double> defect=defects(assignment);
1451 23 : std::vector<double> derivativeByDefect(centers_.size(),0.0);
1452 : double value=0.0;
1453 66 : for(unsigned s=0; s<selected_.size(); ++s) {
1454 43 : const unsigned i=selected_[s];
1455 43 : const double coefficient=coefficients_[s];
1456 : const bool included=
1457 51 : sign_=="ALL" ||
1458 43 : (sign_=="POSITIVE" && defect[i]>0.0) ||
1459 6 : (sign_=="NEGATIVE" && defect[i]<0.0);
1460 4 : if(!included) {
1461 4 : continue;
1462 : }
1463 39 : if(power_==2) {
1464 14 : value+=coefficient*defect[i]*defect[i];
1465 14 : derivativeByDefect[i]+=2.0*coefficient*defect[i];
1466 : } else {
1467 25 : value+=coefficient*defect[i];
1468 25 : derivativeByDefect[i]+=coefficient;
1469 : }
1470 : }
1471 :
1472 23 : std::vector<Vector> derivatives(getNumberOfAtoms());
1473 23 : Tensor virial;
1474 23 : addAssignmentDerivatives(assignment,derivativeByDefect,
1475 : derivatives,virial);
1476 23 : finalize(value,derivatives,virial);
1477 46 : }
1478 :
1479 : //+PLUMEDOC COLVAR VORONOI_DISTANCE
1480 : /*
1481 : Calculate a distance-weighted product of smooth coordination defects.
1482 :
1483 : This Action combines the continuously changing identities defined by
1484 : [VORONOI_COORDINATION](VORONOI_COORDINATION.md) with physical center-center
1485 : distances. It is useful when the progress of a reactive event depends not
1486 : only on whether defects exist, but also on how far apart the corresponding
1487 : sites are.
1488 :
1489 : ## Definition and group semantics
1490 :
1491 : Let $q_i$ be the smooth coordination defect of center $i$. When both
1492 : GROUP1 and GROUP2 are supplied, the Action returns
1493 :
1494 : $$
1495 : D(G_1,G_2)=-\sum_{i\in G_1}\sum_{k\in G_2}d_{ik}q_iq_k .
1496 : $$
1497 :
1498 : GROUP1 and GROUP2 are explicit, nonempty, disjoint subsets of CENTERS. When
1499 : GROUP2 is omitted, unique pairs within GROUP1 are used:
1500 :
1501 : $$
1502 : D(G_1)=-\sum_{\substack{i,k\in G_1\\i<k}}d_{ik}q_iq_k .
1503 : $$
1504 :
1505 : The center-center distance $d_{ik}$ uses the minimum image unless NOPBC
1506 : is specified. CENTERS, ASSIGNED, KAPPA, REFERENCE, NOPBC, and NLIST have the
1507 : same meaning as on the VORONOI_COORDINATION page. NLIST truncates only the
1508 : CENTER-ASSIGNED candidates used to construct the defects; it never removes
1509 : GROUP1/GROUP2 center-center reduction pairs.
1510 :
1511 : The product $q_iq_k$ emphasizes pairs of centers carrying correlated or
1512 : complementary defects without assigning a permanent ion identity. For one
1513 : positive and one negative defect, the overall minus sign makes the leading
1514 : term positive and approximately equal to their separation when the defects
1515 : approach +1 and -1. Multiple defects generate a sum over all requested
1516 : pairs; the output should not then be interpreted automatically as one unique
1517 : ion-ion distance.
1518 :
1519 : ## Derivatives, periodicity, and bias suitability
1520 :
1521 : Both contributions to the derivative are included: the direct derivative of
1522 : $d_{ik}$ and the chain-rule derivative of every $q_i$ through the
1523 : soft assignment. Analytical box derivatives are also provided using the
1524 : same minimum-image vectors. The Action stops if a requested center-center
1525 : distance is zero or non-finite because its radial derivative is then
1526 : undefined.
1527 :
1528 : The exact full-pair assignment is smooth away from coincident atoms and the
1529 : usual minimum-image branch surfaces. NLIST uses the approximate truncated
1530 : normalization described on the VORONOI_COORDINATION page and may add
1531 : membership discontinuities. Validate the exact Action with numerical
1532 : derivatives before biasing it, and converge values and forces separately if
1533 : NLIST is required.
1534 :
1535 : VORONOI_DISTANCE is not normalized by the amount of defect. It naturally
1536 : approaches zero in a neutral configuration, but its magnitude also changes as
1537 : defects form or disappear. This coupling is often desired in a reaction
1538 : coordinate; if a pure conditional distance is intended, its low-defect limit
1539 : must be defined explicitly rather than obtained by dividing by a nearly zero
1540 : weight.
1541 :
1542 : ## Choosing GROUP1 and GROUP2
1543 :
1544 : - Use GROUP1 without GROUP2 for separation among possible sites in one pool,
1545 : such as water O atoms that can host a hydronium or hydroxide defect.
1546 : - Use both groups for cross separation between two chemical pools, such as
1547 : water O atoms and one reactive N or O site.
1548 : - Combine multiple VORONOI_DISTANCE Actions with
1549 : [COMBINE](COMBINE.md) when a published CV contains several physically
1550 : distinct pair sets. Keeping the pair sets explicit avoids hidden
1551 : chemistry and atom-order rules.
1552 :
1553 : ## Worked example 1: water self-ion separation
1554 :
1555 : For one hydronium-hydroxide pair distributed over water oxygen atoms, unique
1556 : pairs within WaterO describe the solution ion-ion separation.
1557 :
1558 : ```plumed
1559 : UNITS LENGTH=A
1560 : WaterO: GROUP ATOMS=1-4
1561 : WaterH: GROUP ATOMS=5-12
1562 :
1563 : ionization: VORONOI_COORDINATION CENTERS=WaterO ASSIGNED=WaterH KAPPA=5 REFERENCE=2 POWER=2
1564 : solution_distance: VORONOI_DISTANCE CENTERS=WaterO ASSIGNED=WaterH KAPPA=5 REFERENCE=2 GROUP1=WaterO
1565 : PRINT ARG=ionization,solution_distance FILE=COLVAR
1566 : ```
1567 :
1568 : Check the CV on neutral water, a separated ion pair, and intermediate proton
1569 : transfer frames. A near-zero value in neutral water is expected and does not
1570 : mean the distance calculation failed.
1571 :
1572 : ## Worked example 2: one reactive O and one reactive N
1573 :
1574 : The following example shows that reactive sites are selected by atom number,
1575 : not by being the last atoms in CENTERS. The reference values are illustrative
1576 : and must be replaced by the chemically intended occupancies.
1577 :
1578 : ```plumed
1579 : UNITS LENGTH=A
1580 : WaterO: GROUP ATOMS=1-3
1581 : ReactiveO: GROUP ATOMS=4
1582 : ReactiveN: GROUP ATOMS=5
1583 : TransferableH: GROUP ATOMS=6-13
1584 : Centers: GROUP ATOMS=ReactiveN,WaterO,ReactiveO
1585 :
1586 : activity: VORONOI_COORDINATION CENTERS=Centers ASSIGNED=TransferableH KAPPA=5 REFERENCE=1,2,2,2,1 POWER=2
1587 : water_to_O: VORONOI_DISTANCE CENTERS=Centers ASSIGNED=TransferableH KAPPA=5 REFERENCE=1,2,2,2,1 GROUP1=WaterO GROUP2=ReactiveO
1588 : water_to_N: VORONOI_DISTANCE CENTERS=Centers ASSIGNED=TransferableH KAPPA=5 REFERENCE=1,2,2,2,1 GROUP1=WaterO GROUP2=ReactiveN
1589 : PRINT ARG=activity,water_to_O,water_to_N FILE=COLVAR
1590 : ```
1591 :
1592 : If only one site is present, remove the other site from CENTERS and remove its
1593 : matching REFERENCE entry. No C++ change is required.
1594 :
1595 : ## Worked example 3: solvated glycine
1596 :
1597 : The glycine distance used in \cite Zhang2024Glycine contains a water-glycine
1598 : cross term and an internal N-O term. Writing them separately removes the
1599 : legacy `NRX` and last-three-atoms convention. The numbered fixture below
1600 : contains zwitterionic [Z] glycine, one complete water near its ammonium group,
1601 : and one complete water near its carboxylate group. Glycine C atoms and C-H
1602 : atoms remain in the coordinate file to make the topology readable, but the
1603 : non-transferable C-H atoms are excluded from ASSIGNED.
1604 :
1605 : ```plumed
1606 : UNITS LENGTH=A
1607 : WaterO: GROUP ATOMS=1,4
1608 : WaterH: GROUP ATOMS=2,3,5,6
1609 : GlyN: GROUP ATOMS=7
1610 : GlyH: GROUP ATOMS=8,9,10
1611 : GlyO1: GROUP ATOMS=15
1612 : GlyO2: GROUP ATOMS=16
1613 : AllH: GROUP ATOMS=WaterH,GlyH
1614 : Centers: GROUP ATOMS=WaterO,GlyN,GlyO1,GlyO2
1615 :
1616 : sp: VORONOI_COORDINATION CENTERS=Centers ASSIGNED=AllH KAPPA=5 REFERENCE=2,2,2,0.5,0.5 POWER=1 COEFFICIENTS=1,1,2,2,2
1617 : sd_water: VORONOI_DISTANCE CENTERS=Centers ASSIGNED=AllH KAPPA=5 REFERENCE=2,2,2,0.5,0.5 GROUP1=WaterO GROUP2=GlyN,GlyO1,GlyO2
1618 : sd_internal: VORONOI_DISTANCE CENTERS=Centers ASSIGNED=AllH KAPPA=5 REFERENCE=2,2,2,0.5,0.5 GROUP1=GlyN GROUP2=GlyO1,GlyO2
1619 : sd: COMBINE ARG=sd_water,sd_internal COEFFICIENTS=1,1 PERIODIC=NO
1620 :
1621 : PRINT ARG=sp,sd_water,sd_internal,sd FILE=COLVAR
1622 : ```
1623 :
1624 : CENTERS expands to atoms `(1,4,7,15,16)`, so the five REFERENCE entries map to
1625 : two water O atoms, glycine N, and two symmetry-related glycine O atoms. The
1626 : vector `(2,2,2,0.5,0.5)` declares the neutral [N] chemical origin even though
1627 : the fixture is [Z]. An ideal [Z] frame therefore has defects near
1628 : `(0,0,+1,-0.5,-0.5)`, and the seven references sum to the seven atoms in AllH.
1629 : The complete full-pair fixture gives `sp=0.0073`, `sd_water=0.0354`
1630 : $\mathrm{\AA}$, `sd_internal=2.9789` $\mathrm{\AA}$, and
1631 : `sd=3.0143` $\mathrm{\AA}$.
1632 :
1633 : An alternative reference partition can accidentally give the same linear
1634 : `sp` scalar when coefficients and total references cancel. That does not
1635 : make it equivalent for VORONOI_DISTANCE, POWER=2, or sign filtering. Use one
1636 : physically declared reference vector consistently across the coupled Actions.
1637 : The coordinate files, runnable inputs, commands, and expected outputs are in
1638 : `user-doc/tutorials/others/reactive-voronoi`.
1639 : The complete application is available in the
1640 : [GlycineTautomerism repository](https://github.com/Zhang-pchao/GlycineTautomerism/tree/main/Enhanced_Sampling).
1641 :
1642 : ## Worked example 4: nitrogen-reduction proton transfer
1643 :
1644 : The published Ru single-atom nitrogen-reduction input in
1645 : \cite Zhang2026NRR uses water oxygen atoms and one reactive nitrogen as the
1646 : active Voronoi centers. Catalyst and other adsorbate groups may appear in the
1647 : full simulation input, but they are not automatically active in these Actions.
1648 :
1649 : ```plumed
1650 : UNITS LENGTH=A
1651 : WaterO: GROUP ATOMS=1-3
1652 : ReactiveN: GROUP ATOMS=4
1653 : TransferableH: GROUP ATOMS=5-11
1654 : Centers: GROUP ATOMS=WaterO,ReactiveN
1655 :
1656 : solvation: VORONOI_COORDINATION CENTERS=Centers ASSIGNED=TransferableH KAPPA=5 REFERENCE=2,2,2,1 POWER=1 COEFFICIENTS=1,1,1,2
1657 : solution_ion_distance: VORONOI_DISTANCE CENTERS=Centers ASSIGNED=TransferableH KAPPA=5 REFERENCE=2,2,2,1 GROUP1=WaterO
1658 : site_ion_distance: VORONOI_DISTANCE CENTERS=Centers ASSIGNED=TransferableH KAPPA=5 REFERENCE=2,2,2,1 GROUP1=WaterO GROUP2=ReactiveN
1659 : PRINT ARG=solvation,solution_ion_distance,site_ion_distance FILE=COLVAR
1660 : ```
1661 :
1662 : The corresponding files are available in the
1663 : [OPES-DPMD-NRR repository](https://github.com/Zhang-pchao/research/tree/main/OPES-DPMD-NRR).
1664 : If the full input uses OPES, enable the optional module as described on the
1665 : [VORONOI_COORDINATION](VORONOI_COORDINATION.md) page.
1666 :
1667 : ## Troubleshooting
1668 :
1669 : - `GROUP1 must contain at least two atoms when GROUP2 is omitted`: a
1670 : within-group distance requires at least one unique pair. For a single
1671 : reactive site, use a nonempty disjoint GROUP2.
1672 : - `GROUP1 and GROUP2 must be disjoint`: a self-pair has no valid radial
1673 : derivative and is intentionally rejected.
1674 : - A distance is unexpectedly negative: same-sign defect pairs contribute
1675 : negatively because of the leading minus sign. Inspect the individual
1676 : defects and pair-set definition before changing the sign convention.
1677 : - A value grows with the number of defects: the Action is a pair sum, not a
1678 : normalized nearest-pair distance. Reduce the groups or define a different
1679 : mathematical observable if that behavior is not intended.
1680 : - A periodic trajectory jumps: inspect minimum-image branch crossings and
1681 : whether an unwrapped or externally referenced coordinate is actually needed.
1682 : */
1683 : //+ENDPLUMEDOC
1684 :
1685 : class VoronoiDistance : public SoftVoronoiBase {
1686 : std::vector<std::pair<unsigned,unsigned> > reductionPairs_;
1687 :
1688 : public:
1689 : explicit VoronoiDistance(const ActionOptions&);
1690 : static void registerKeywords(Keywords&);
1691 : void calculate() override;
1692 : };
1693 :
1694 : PLUMED_REGISTER_ACTION(VoronoiDistance,"VORONOI_DISTANCE")
1695 :
1696 6 : void VoronoiDistance::registerKeywords(Keywords& keys) {
1697 6 : Colvar::registerKeywords(keys);
1698 6 : registerCommonKeywords(keys);
1699 6 : keys.add("atoms","GROUP1","First explicit subset of CENTERS");
1700 6 : keys.add("atoms","GROUP2","Optional disjoint second subset of CENTERS; if omitted, unique pairs within GROUP1 are used");
1701 6 : setScalarDescription(keys,"the distance-weighted product of selected smooth occupancy defects");
1702 6 : }
1703 :
1704 2 : VoronoiDistance::VoronoiDistance(const ActionOptions& ao):
1705 : Action(ao),
1706 2 : SoftVoronoiBase(ao) {
1707 2 : std::vector<AtomNumber> group1Atoms;
1708 2 : std::vector<AtomNumber> group2Atoms;
1709 2 : parseAtomList("GROUP1",group1Atoms);
1710 2 : parseAtomList("GROUP2",group2Atoms);
1711 : const std::vector<unsigned> group1=
1712 4 : mapSelection(group1Atoms,"GROUP1",false);
1713 :
1714 2 : if(group2Atoms.empty()) {
1715 1 : if(group1.size()<2) {
1716 0 : error("GROUP1 must contain at least two atoms when GROUP2 is omitted");
1717 : }
1718 3 : for(unsigned i=0; i<group1.size(); ++i) {
1719 3 : for(unsigned j=i+1; j<group1.size(); ++j) {
1720 1 : reductionPairs_.push_back(
1721 1 : std::make_pair(group1[i],group1[j]));
1722 : }
1723 : }
1724 : } else {
1725 : const std::vector<unsigned> group2=
1726 2 : mapSelection(group2Atoms,"GROUP2",false);
1727 1 : std::set<unsigned> group1Ids(group1.begin(),group1.end());
1728 2 : for(unsigned j=0; j<group2.size(); ++j) {
1729 : if(group1Ids.count(group2[j])>0) {
1730 0 : error("GROUP1 and GROUP2 must be disjoint");
1731 : }
1732 : }
1733 2 : for(unsigned i=0; i<group1.size(); ++i) {
1734 2 : for(unsigned j=0; j<group2.size(); ++j) {
1735 1 : reductionPairs_.push_back(
1736 1 : std::make_pair(group1[i],group2[j]));
1737 : }
1738 : }
1739 1 : }
1740 :
1741 2 : finishSetup("VORONOI_DISTANCE");
1742 2 : log.printf(" reducing %u explicit center pairs\n",
1743 : static_cast<unsigned>(reductionPairs_.size()));
1744 2 : }
1745 :
1746 4 : void VoronoiDistance::calculate() {
1747 4 : const Assignment assignment=calculateAssignment();
1748 4 : const std::vector<double> defect=defects(assignment);
1749 4 : std::vector<double> derivativeByDefect(centers_.size(),0.0);
1750 4 : std::vector<Vector> derivatives(getNumberOfAtoms());
1751 4 : Tensor virial;
1752 : double value=0.0;
1753 :
1754 : // The reducer is small for the intended explicit groups. Evaluate its
1755 : // scalar and q derivatives redundantly; only rank zero owns direct forces.
1756 4 : const bool directOwner=ownsDirectDerivatives();
1757 8 : for(unsigned p=0; p<reductionPairs_.size(); ++p) {
1758 4 : const unsigned i=reductionPairs_[p].first;
1759 4 : const unsigned k=reductionPairs_[p].second;
1760 4 : const Vector distance=pbc_ ?
1761 2 : pbcDistance(getPosition(i),getPosition(k)) :
1762 2 : delta(getPosition(i),getPosition(k));
1763 4 : const double length=distance.modulo();
1764 4 : if(!std::isfinite(length) ||
1765 : length<=std::numeric_limits<double>::epsilon()) {
1766 0 : error("a GROUP1-GROUP2 center distance is zero or non-finite");
1767 : }
1768 :
1769 4 : value-=length*defect[i]*defect[k];
1770 4 : derivativeByDefect[i]-=length*defect[k];
1771 4 : derivativeByDefect[k]-=length*defect[i];
1772 4 : if(directOwner) {
1773 : const Vector pairDerivative=
1774 4 : (-defect[i]*defect[k]/length)*distance;
1775 : derivatives[i]-=pairDerivative;
1776 : derivatives[k]+=pairDerivative;
1777 4 : virial-=Tensor(pairDerivative,distance);
1778 : }
1779 : }
1780 :
1781 4 : addAssignmentDerivatives(assignment,derivativeByDefect,
1782 : derivatives,virial);
1783 4 : finalize(value,derivatives,virial);
1784 12 : }
1785 :
1786 : //+PLUMEDOC COLVAR VORONOI_POSITION
1787 : /*
1788 : Calculate a defect-weighted Cartesian position relative to a fixed origin.
1789 :
1790 : Molecular identity alone does not say where a hydronium, hydroxide, or another
1791 : coordination defect is located. This Action weights the coordinates of the
1792 : possible host centers by the smooth defects defined in
1793 : [VORONOI_COORDINATION](VORONOI_COORDINATION.md). It replaces atom-list-index
1794 : moments with a physical Cartesian coordinate and provides analytical
1795 : coordinate and box derivatives.
1796 :
1797 : Legacy index-weighted observables such as $\sum_i i q_i^2$ depend on the
1798 : order in which centers are listed: reordering chemically identical input can
1799 : change the value even when the configuration is unchanged. Such an index is
1800 : useful at most as an internal diagnostic and is a poor general bias coordinate.
1801 : It is intentionally not reproduced here. Use SELECT and a physical
1802 : VORONOI_POSITION coordinate when the spatial location of a defect is needed.
1803 :
1804 : ## Definition
1805 :
1806 : For selected centers $S$, Cartesian axis $\alpha$, and fixed origin
1807 : $x_0$, define the displacement
1808 :
1809 : $$
1810 : u_i=\operatorname{minimage}(x_{i,\alpha}-x_0)
1811 : $$
1812 :
1813 : when periodic boundaries are active. NOPBC replaces it by the direct
1814 : Cartesian difference. The unnormalized Action is
1815 :
1816 : $$
1817 : P=\sum_{i\in S}g(q_i)f(u_i),
1818 : $$
1819 :
1820 : where $g(q_i)=q_i^2$ for defects retained by SIGN and zero for excluded
1821 : defects. The spatial function is $f(u)=u$ by default and
1822 : $f(u)=|u|$ with ABSOLUTE. SELECT defaults to all CENTERS.
1823 :
1824 : With NORMALIZE, the returned value is
1825 :
1826 : $$
1827 : \bar P=\frac{\sum_i g(q_i)f(u_i)}{\sum_i g(q_i)}.
1828 : $$
1829 :
1830 : The unnormalized form naturally approaches zero in a neutral configuration
1831 : and measures position multiplied by defect activity. The normalized form is
1832 : a defect-weighted mean position, but is defined only when a selected defect is
1833 : present. If its denominator is no larger than TOLERANCE, the Action stops
1834 : instead of returning a NaN or an arbitrarily amplified coordinate.
1835 :
1836 : ## Selecting the observable
1837 :
1838 : - SIGN=POSITIVE tracks over-coordinated centers; SIGN=NEGATIVE tracks
1839 : under-coordinated centers; SIGN=ALL includes both.
1840 : - ABSOLUTE measures distance from the fixed origin. Without ABSOLUTE, the
1841 : sign of the selected Cartesian side is retained.
1842 : - NORMALIZE removes the magnitude of the total squared defect. Use it only
1843 : when the relevant defect is guaranteed to exist throughout the sampled
1844 : region.
1845 : - SELECT restricts possible hosts, for example to water O atoms while leaving
1846 : a reactive molecular site in CENTERS for the shared assignment.
1847 :
1848 : The sign-gated squared weight has a continuous first derivative at
1849 : $q_i=0$ but a discontinuous second derivative. ABSOLUTE is
1850 : non-differentiable at $u_i=0$. Without ABSOLUTE, the direct Cartesian
1851 : derivative is smooth away from the periodic branch cut.
1852 :
1853 : ## Derivatives and bias suitability
1854 :
1855 : The implementation differentiates both the defect weight and the selected
1856 : center coordinate. With NORMALIZE, the quotient rule is evaluated
1857 : analytically. Box derivatives use the same periodic images as the value.
1858 : This makes the Action usable as a bias coordinate within its declared smooth
1859 : region, but the following surfaces require special care:
1860 :
1861 : 1. the minimum-image branch at half the periodic cell length;
1862 : 2. the ABSOLUTE cusp at ORIGIN;
1863 : 3. the second-derivative cusp where a sign-selected defect changes sign;
1864 : 4. an NLIST membership change;
1865 : 5. the low-weight boundary of NORMALIZE.
1866 :
1867 : A numerical derivative check must use a frame away from these surfaces. A
1868 : successful check at one frame does not remove a branch cut elsewhere in the
1869 : sampled domain.
1870 :
1871 : ## Periodic boundaries and reference frames
1872 :
1873 : ORIGIN is a fixed Cartesian coordinate, not an atom, a center of mass, or an
1874 : automatically detected interface. With periodic boundaries, the displacement
1875 : is minimum-image and therefore has a branch cut. This is appropriate only
1876 : when the cell and origin define a reproducible frame.
1877 :
1878 : For a drifting slab, moving droplet, flexible pore, or fluctuating interface,
1879 : first construct a physically justified external reference and keep all
1880 : coordinates in one consistent image convention. VORONOI_POSITION does not
1881 : unwrap trajectories or locate an instantaneous interface. NOPBC removes the
1882 : minimum-image mapping but does not by itself make wrapped molecular-dynamics
1883 : coordinates continuous.
1884 :
1885 : NLIST has the same approximate normalization and convergence requirements as
1886 : [VORONOI_COORDINATION](VORONOI_COORDINATION.md). Position values can magnify
1887 : small assignment errors through the distance from ORIGIN, so converge both
1888 : the CV and its derivatives against full-pair mode.
1889 :
1890 : ## Worked example 1: water self-ions at an interface
1891 :
1892 : For a slab normal to $z$, positive and negative defects can be monitored
1893 : relative to a fixed reference plane. The unnormalized form below remains
1894 : defined in neutral configurations.
1895 :
1896 : ```plumed
1897 : UNITS LENGTH=A
1898 : WaterO: GROUP ATOMS=1-4
1899 : WaterH: GROUP ATOMS=5-12
1900 :
1901 : h3o_z: VORONOI_POSITION CENTERS=WaterO ASSIGNED=WaterH KAPPA=5 REFERENCE=2 AXIS=Z ORIGIN=5 SIGN=POSITIVE ABSOLUTE
1902 : oh_z: VORONOI_POSITION CENTERS=WaterO ASSIGNED=WaterH KAPPA=5 REFERENCE=2 AXIS=Z ORIGIN=5 SIGN=NEGATIVE ABSOLUTE
1903 : PRINT ARG=h3o_z,oh_z FILE=COLVAR
1904 : ```
1905 :
1906 : The two values combine defect amount and distance from the plane. They are
1907 : well suited to monitoring neutral frames because both become zero as their
1908 : selected defect disappears.
1909 :
1910 : ## Worked example 2: normalized ion location
1911 :
1912 : When a positive defect is known to remain present, NORMALIZE returns its
1913 : defect-weighted mean signed location. This input intentionally omits ABSOLUTE
1914 : so the two sides of the reference plane remain distinguishable.
1915 :
1916 : ```plumed
1917 : UNITS LENGTH=A
1918 : WaterO: GROUP ATOMS=1-4
1919 : WaterH: GROUP ATOMS=5-12
1920 :
1921 : h3o_signed_z: VORONOI_POSITION CENTERS=WaterO ASSIGNED=WaterH KAPPA=5 REFERENCE=2 AXIS=Z ORIGIN=5 SIGN=POSITIVE NORMALIZE TOLERANCE=1e-10
1922 : PRINT ARG=h3o_signed_z FILE=COLVAR
1923 : ```
1924 :
1925 : Do not use this normalized form in a trajectory that can become neutral. In
1926 : that case the fail-closed TOLERANCE error is expected; use the unnormalized
1927 : form or a separately justified state-dependent workflow.
1928 :
1929 : ## Worked example 3: restrict hosts while retaining a reactive site
1930 :
1931 : All centers participate in assigning the transferable hydrogens, while SELECT
1932 : restricts the reported location to water oxygen atoms. This is useful when a
1933 : reactive molecular or catalytic site must compete for protons but should not
1934 : be interpreted as an interfacial water ion.
1935 :
1936 : ```plumed
1937 : UNITS LENGTH=A
1938 : WaterO: GROUP ATOMS=1-3
1939 : ReactiveN: GROUP ATOMS=4
1940 : TransferableH: GROUP ATOMS=5-11
1941 : Centers: GROUP ATOMS=WaterO,ReactiveN
1942 :
1943 : water_positive_z: VORONOI_POSITION CENTERS=Centers ASSIGNED=TransferableH SELECT=WaterO KAPPA=5 REFERENCE=2,2,2,1 AXIS=Z ORIGIN=5 SIGN=POSITIVE ABSOLUTE
1944 : PRINT ARG=water_positive_z FILE=COLVAR
1945 : ```
1946 :
1947 : The air-water and oil-water applications are available in the
1948 : [OilWaterInterface repository](https://github.com/Zhang-pchao/OilWaterInterface/tree/main)
1949 : and are described in \cite Zhang2025Interfaces. Historical NL_CUTOFF,
1950 : NL_STRIDE, and ORIGIN values are system-specific and should not be copied
1951 : without new full-pair and reference-frame validation.
1952 :
1953 : ## Validation checklist and troubleshooting
1954 :
1955 : - Visualize the selected centers and ORIGIN in representative periodic frames.
1956 : - Translate the system by a lattice vector and verify the same minimum-image
1957 : value.
1958 : - Place a test center on both sides of a boundary and map the expected branch
1959 : explicitly.
1960 : - Compare analytical and numerical coordinate and box derivatives away from
1961 : the branch cut, ABSOLUTE cusp, and sign boundary.
1962 : - Confirm that NORMALIZE remains safely above TOLERANCE over the intended
1963 : trajectory.
1964 : - If the sign is reversed, inspect AXIS, ORIGIN, image convention, and whether
1965 : ABSOLUTE was intended.
1966 : - If the value jumps by about a box length, the trajectory crossed the
1967 : minimum-image branch; this is a coordinate-design problem, not a force
1968 : precision problem.
1969 : - If a drifting interface is the true reference, define and validate that
1970 : moving reference before applying a bias. Changing ORIGIN to a guessed
1971 : constant only hides the drift.
1972 : */
1973 : //+ENDPLUMEDOC
1974 :
1975 : class VoronoiPosition : public SoftVoronoiBase {
1976 : unsigned axis_;
1977 : double origin_;
1978 : double tolerance_;
1979 : bool absolute_;
1980 : bool normalize_;
1981 : std::string sign_;
1982 : std::vector<unsigned> selected_;
1983 :
1984 : public:
1985 : explicit VoronoiPosition(const ActionOptions&);
1986 : static void registerKeywords(Keywords&);
1987 : void calculate() override;
1988 : };
1989 :
1990 : PLUMED_REGISTER_ACTION(VoronoiPosition,"VORONOI_POSITION")
1991 :
1992 8 : void VoronoiPosition::registerKeywords(Keywords& keys) {
1993 8 : Colvar::registerKeywords(keys);
1994 8 : registerCommonKeywords(keys);
1995 8 : keys.add("atoms","SELECT","Subset of CENTERS included in the position; the default is all centers");
1996 8 : keys.add("compulsory","AXIS","Cartesian axis X, Y, or Z");
1997 8 : keys.add("compulsory","ORIGIN","Fixed Cartesian origin coordinate");
1998 8 : keys.add("compulsory","SIGN","ALL","Use ALL, POSITIVE, or NEGATIVE defects");
1999 8 : keys.add("compulsory","TOLERANCE","1e-12","Minimum total weight accepted by NORMALIZE");
2000 8 : keys.addFlag("ABSOLUTE",false,"Use the absolute displacement from ORIGIN");
2001 8 : keys.addFlag("NORMALIZE",false,"Divide by the total selected defect weight");
2002 8 : setScalarDescription(keys,"the selected defect-weighted Cartesian position relative to ORIGIN");
2003 8 : }
2004 :
2005 4 : VoronoiPosition::VoronoiPosition(const ActionOptions& ao):
2006 : Action(ao),
2007 : SoftVoronoiBase(ao),
2008 4 : axis_(0),
2009 4 : origin_(0.0),
2010 4 : tolerance_(1e-12),
2011 4 : absolute_(false),
2012 4 : normalize_(false),
2013 4 : sign_("ALL") {
2014 4 : std::vector<AtomNumber> selectedAtoms;
2015 4 : parseAtomList("SELECT",selectedAtoms);
2016 8 : selected_=mapSelection(selectedAtoms,"SELECT",true);
2017 :
2018 : std::string axis;
2019 8 : parse("AXIS",axis);
2020 4 : std::transform(axis.begin(),axis.end(),axis.begin(),
2021 : [](const char value) {
2022 4 : return static_cast<char>(std::toupper(static_cast<unsigned char>(value)));
2023 : });
2024 4 : if(axis=="X") {
2025 4 : axis_=0;
2026 0 : } else if(axis=="Y") {
2027 0 : axis_=1;
2028 0 : } else if(axis=="Z") {
2029 0 : axis_=2;
2030 : } else {
2031 0 : error("AXIS must be X, Y, or Z");
2032 : }
2033 :
2034 4 : parse("ORIGIN",origin_);
2035 4 : if(!std::isfinite(origin_)) {
2036 0 : error("ORIGIN must be finite");
2037 : }
2038 4 : parse("SIGN",sign_);
2039 4 : sign_=normalizedSign(sign_);
2040 4 : parse("TOLERANCE",tolerance_);
2041 4 : if(!std::isfinite(tolerance_) || tolerance_<=0.0) {
2042 0 : error("TOLERANCE must be finite and positive");
2043 : }
2044 4 : parseFlag("ABSOLUTE",absolute_);
2045 4 : parseFlag("NORMALIZE",normalize_);
2046 :
2047 4 : finishSetup("VORONOI_POSITION");
2048 4 : log.printf(" reducing %u selected centers along axis %u from origin %g\n",
2049 4 : static_cast<unsigned>(selected_.size()),axis_+1,origin_);
2050 4 : }
2051 :
2052 8 : void VoronoiPosition::calculate() {
2053 8 : const Assignment assignment=calculateAssignment();
2054 8 : const std::vector<double> defect=defects(assignment);
2055 8 : std::vector<double> displacement(selected_.size(),0.0);
2056 8 : std::vector<double> function(selected_.size(),0.0);
2057 8 : std::vector<double> weight(selected_.size(),0.0);
2058 8 : std::vector<double> weightDerivative(selected_.size(),0.0);
2059 8 : std::vector<double> functionDerivative(selected_.size(),0.0);
2060 :
2061 : double numerator=0.0;
2062 : double denominator=0.0;
2063 22 : for(unsigned s=0; s<selected_.size(); ++s) {
2064 14 : const unsigned i=selected_[s];
2065 14 : Vector reference=getPosition(i);
2066 14 : reference[axis_]=origin_;
2067 14 : const Vector relative=pbc_ ?
2068 : pbcDistance(reference,getPosition(i)) :
2069 : delta(reference,getPosition(i));
2070 14 : displacement[s]=relative[axis_];
2071 14 : if(!std::isfinite(displacement[s])) {
2072 0 : error("VORONOI_POSITION displacement is non-finite");
2073 : }
2074 14 : function[s]=absolute_ ?
2075 : std::fabs(displacement[s]) : displacement[s];
2076 14 : functionDerivative[s]=absolute_ ?
2077 8 : (displacement[s]>0.0 ? 1.0 :
2078 4 : (displacement[s]<0.0 ? -1.0 : 0.0)) : 1.0;
2079 :
2080 : const bool included=
2081 26 : sign_=="ALL" ||
2082 14 : (sign_=="POSITIVE" && defect[i]>0.0) ||
2083 8 : (sign_=="NEGATIVE" && defect[i]<0.0);
2084 : if(included) {
2085 8 : weight[s]=defect[i]*defect[i];
2086 8 : weightDerivative[s]=2.0*defect[i];
2087 8 : numerator+=weight[s]*function[s];
2088 8 : denominator+=weight[s];
2089 : }
2090 : }
2091 :
2092 8 : if(normalize_ && (!std::isfinite(denominator) ||
2093 2 : denominator<=tolerance_)) {
2094 0 : error("VORONOI_POSITION total weight is below TOLERANCE");
2095 : }
2096 8 : const double value=normalize_ ? numerator/denominator : numerator;
2097 8 : std::vector<double> derivativeByDefect(centers_.size(),0.0);
2098 8 : std::vector<Vector> derivatives(getNumberOfAtoms());
2099 8 : Tensor virial;
2100 8 : const bool directOwner=ownsDirectDerivatives();
2101 22 : for(unsigned s=0; s<selected_.size(); ++s) {
2102 14 : const unsigned i=selected_[s];
2103 14 : if(weight[s]==0.0 && weightDerivative[s]==0.0) {
2104 6 : continue;
2105 : }
2106 8 : derivativeByDefect[i]+=
2107 8 : normalize_ ?
2108 2 : weightDerivative[s]*(function[s]-value)/denominator :
2109 6 : weightDerivative[s]*function[s];
2110 :
2111 8 : if(directOwner) {
2112 : const double directScalar=
2113 8 : (normalize_ ? weight[s]/denominator : weight[s])*
2114 8 : functionDerivative[s];
2115 8 : Vector directDerivative;
2116 8 : directDerivative[axis_]=directScalar;
2117 : derivatives[i]+=directDerivative;
2118 8 : Vector imagePosition=getPosition(i);
2119 8 : if(pbc_) {
2120 2 : imagePosition[axis_]=origin_+displacement[s];
2121 : }
2122 16 : virial-=Tensor(imagePosition,directDerivative);
2123 : }
2124 : }
2125 :
2126 8 : addAssignmentDerivatives(assignment,derivativeByDefect,
2127 : derivatives,virial);
2128 8 : finalize(value,derivatives,virial);
2129 16 : }
2130 :
2131 : }
2132 : }
|
