Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2023 of Nicholas Herringer and Siva Dasetty.
3 :
4 : The PINES module is free software: you can redistribute it and/or modify
5 : it under the terms of the GNU Lesser General Public License as published by
6 : the Free Software Foundation, either version 3 of the License, or
7 : (at your option) any later version.
8 :
9 : The PINES module is distributed in the hope that it will be useful,
10 : but WITHOUT ANY WARRANTY; without even the implied warranty of
11 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 : GNU Lesser General Public License for more details.
13 :
14 : You should have received a copy of the GNU Lesser General Public License
15 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
16 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
17 : #include <string>
18 : #include <cmath>
19 : #include <iostream>
20 : #include <fstream>
21 : #include <cstdio>
22 : #include <unordered_map>
23 : #include <set>
24 : #include <numeric>
25 : #ifdef _OPENMP
26 : #include <omp.h>
27 : #endif
28 : #include <unordered_set>
29 : #include <utility>
30 : #include <functional>
31 : #include <vector>
32 : #include <tuple>
33 : #include <cstddef>
34 :
35 : #include "core/ActionWithVirtualAtom.h"
36 : #include "core/ActionWithValue.h"
37 : #include "core/PlumedMain.h"
38 : #include "tools/Stopwatch.h"
39 : #include "tools/PDB.h"
40 : #include "tools/Vector.h"
41 : #include "tools/AtomNumber.h"
42 : #include "tools/SwitchingFunction.h"
43 : #include "tools/Communicator.h"
44 : #include "tools/Units.h"
45 :
46 : #include "colvar/Colvar.h"
47 : #include "core/ActionRegister.h"
48 :
49 :
50 :
51 : namespace PLMD {
52 : namespace pines {
53 :
54 : struct AtomNumberLess {
55 : bool operator()(const AtomNumber& a, const AtomNumber& b) const {
56 888612 : return a.index() < b.index();
57 : }
58 : };
59 :
60 : struct AtomPair {
61 : AtomNumber alpha_group;
62 : AtomNumber beta_group;
63 :
64 : // equality
65 : bool operator==(const AtomPair& o) const noexcept {
66 0 : return alpha_group == o.alpha_group && beta_group == o.beta_group;
67 : }
68 :
69 : // lexicographic ordering for std::set/std::map and sorting
70 1925144 : bool operator<(const AtomPair& o) const noexcept {
71 1925144 : return std::tie(alpha_group, beta_group) < std::tie(o.alpha_group, o.beta_group);
72 : }
73 : };
74 :
75 : // non-member convenience operators
76 : inline bool operator!=(const AtomPair& a, const AtomPair& b) noexcept {
77 : return !(a == b);
78 : }
79 : inline bool operator>(const AtomPair& a, const AtomPair& b) noexcept {
80 : return b < a;
81 : }
82 : inline bool operator<=(const AtomPair& a, const AtomPair& b) noexcept {
83 : return !(b < a);
84 : }
85 : inline bool operator>=(const AtomPair& a, const AtomPair& b) noexcept {
86 : return !(a < b);
87 : }
88 :
89 :
90 : class PINES : public PLMD::colvar::Colvar {
91 : private:
92 : PLMD::Stopwatch timer;
93 : int N_Blocks;
94 : int total_PIV_length;
95 : int inited_step;
96 : bool driver_mode;
97 : bool freeze_selection;
98 : std::vector<int> steps_since_update;
99 : std::vector<int> nstride;
100 : std::string ref_file;
101 : std::vector<SwitchingFunction> sfs;
102 : std::vector<std::string> sw;
103 : std::vector<double> r00;
104 : std::vector<std:: vector<double> > PIV;
105 : std::vector<std:: vector<Vector> > ann_deriv;
106 :
107 : std::vector<std:: vector<AtomNumber> > listall;
108 : std::vector<std:: vector<AtomNumber> > listreduced;
109 : std::set<AtomNumber, AtomNumberLess> listreducedall;
110 : std::vector<AtomNumber> listreducedall_vec;
111 : std::unordered_map<int,int> atom_ind_hashmap;
112 :
113 : std::vector<bool> stale_tolerance;
114 : PDB mypdb;
115 : std::vector<std::string> block_params;
116 : std::vector<std::array<std::vector<AtomNumber>,2> > block_groups_atom_list;
117 : std::vector<int> block_lengths;
118 : std::vector<int> Buffer_Pairs;
119 : std::vector<int> tot_num_pairs;
120 : std::vector<std::array<std::array<bool,3>,2> > input_filters;
121 : std::vector<double> delta_pd;
122 : std::vector<double> r_tolerance;
123 : std::vector<std::vector<Vector> > PL_atoms_ref_coords;
124 : std::vector<std::vector<AtomPair > > Exclude_Pairs;
125 : std::vector<std::array<std::vector<std::string>, 2>> Name_list;
126 : std::vector<std::array<std::vector<AtomNumber>, 2>> ID_list;
127 : std::vector<std::array<std::vector<int>, 2>> ResID_list;
128 : std::vector<char> all_g1g2_pairs;
129 : std::vector<std::vector<std::pair<double, AtomPair > > > vecMaxHeapVecs;
130 : std::vector<std::vector<AtomPair > > latched_pairs;
131 : std::vector<int> preupdated_block;
132 : std::vector<bool> isFirstBuild;
133 :
134 : bool atomMatchesFilters(int n, int g, AtomNumber ind, int resid, const std::string& atom_name);
135 : void buildMaxHeapVecBlock(int n, const PDB& refpdb, std::vector<std::pair<double, AtomPair>>& heap);
136 : void updateBlockPairList(int n, std::vector<std::pair<double, AtomPair>>& heap);
137 : double calculateDistance(int n, const AtomNumber& ind0, const AtomNumber& ind1, const PDB& refpdb);
138 : std::ofstream log;
139 : bool ensureBlockUpdated(int n);
140 : void latchFromCurrentHeaps();
141 : void logMsg(const std::string& msg, const std::string& section);
142 : void logMsg(const Vector& vec, const std::string& section);
143 : void resizeAllContainers(int N);
144 :
145 : public:
146 : static void registerKeywords( Keywords& keys );
147 : explicit PINES(const ActionOptions&);
148 : //~PINES();
149 : // active methods:
150 : struct MinCompareDist {
151 : bool operator()(const std::pair<double, AtomPair>& p1, const std::pair<double, AtomPair>& p2) {
152 908703 : return p1.first < p2.first; // Min heap
153 : }
154 : };
155 : struct MaxCompareDist {
156 : bool operator()(const std::pair<double, AtomPair>& p1, const std::pair<double, AtomPair>& p2) {
157 4632 : return p1.first > p2.first; // Max heap
158 : }
159 : };
160 :
161 : void calculate() override;
162 0 : void checkFieldsAllowed() override {}
163 : // -- SD prepare to requestAtoms during simulation
164 : void prepare() override;
165 : };
166 : //+PLUMEDOC PINES
167 : //Documentation to be added.
168 : //+ENDPLUMEDOC PINES
169 : PLUMED_REGISTER_ACTION(PINES, "PINES")
170 : constexpr int CACHE_LINE_SIZE = 64;
171 :
172 : template<typename T>
173 : struct alignas(CACHE_LINE_SIZE) Padded {
174 : T value;
175 :
176 : // Optional constructors
177 : Padded() : value() {}
178 : explicit Padded(const T& val) : value(val) {}
179 : };
180 :
181 5 : void PINES::registerKeywords(Keywords &keys) {
182 5 : Colvar::registerKeywords(keys);
183 5 : keys.add("numbered", "SWITCH", "The switching functions parameter. You must specify a Switching function for all PINES blocks."
184 : "Details of the various switching functions you can use are provided in the documentation for LESS_THAN.");
185 5 : keys.add("numbered", "BLOCK", "Each block of the PIV");
186 5 : keys.add("compulsory", "REF_FILE", "PDB file name that contains the information about system connectivity and labels.");
187 5 : keys.add("compulsory", "N_BLOCKS", "Number of blocks in PIV");
188 5 : keys.add("compulsory", "BUFFER", "Number of additional pairwise distances to include as a buffer for each PIV block");
189 5 : keys.add("compulsory", "PL_REFRESH", "Upper limit refresh rate for each block pair list");
190 5 : keys.add("compulsory", "SIZE", "Length of each PIV Block");
191 5 : keys.add("optional", "ATOMID", "AtomIDs");
192 5 : keys.add("optional", "RESID", "ResIDs");
193 5 : keys.add("optional", "NAME", "Atom Names");
194 5 : keys.add("optional", "EXCLUDE_PAIRS", "Excluded pairs");
195 5 : keys.addFlag("DRIVERMODE", false, "Use when post-processing a trajectory. This will force the pair list to be rebuilt at every step (i.e. no pair list)");
196 5 : keys.addFlag("FREEZE_SELECTION", false, "Freeze top-K membership during derivative checks");
197 10 : keys.addOutputComponent("ELEMENT", "default", "scalar", "Elements of the PINES block");
198 10 : keys.reset_style("SWITCH", "compulsory");
199 5 : }
200 :
201 9156 : bool PINES::atomMatchesFilters(int n, int g, AtomNumber ind, int resid, const std::string& atom_name) {
202 9156 : bool id_check = !input_filters[n][g][0];
203 9156 : bool res_check = !input_filters[n][g][1];
204 9156 : bool name_check = !input_filters[n][g][2];
205 :
206 9156 : if (input_filters[n][g][0]) {
207 0 : for (auto& atomID : ID_list[n][g]) {
208 0 : if (ind == atomID) {
209 : id_check = true;
210 : break;
211 : }
212 : }
213 : }
214 9156 : if (input_filters[n][g][1]) {
215 0 : for (auto& rID : ResID_list[n][g]) {
216 0 : if (resid == rID) {
217 : res_check = true;
218 : break;
219 : }
220 : }
221 : }
222 9156 : if (input_filters[n][g][2]) {
223 17802 : for (auto& name : Name_list[n][g]) {
224 10174 : if (atom_name == name) {
225 : name_check = true;
226 : break;
227 : }
228 : }
229 : }
230 9156 : return id_check && res_check && name_check;
231 : }
232 :
233 106 : void PINES::buildMaxHeapVecBlock(int n, const PDB& refpdb, std::vector<std::pair<double, AtomPair>>& heap) {
234 : #ifdef _OPENMP
235 106 : int nthreads = omp_get_max_threads();
236 : #else
237 : int nthreads = 1;
238 : #endif
239 106 : std::vector<std::vector<std::pair<double, AtomPair>>> thread_heaps(nthreads);
240 :
241 :
242 : // Rough estimate, tweak as needed
243 106 : int estimated_pairs_per_thread = block_groups_atom_list[n][0].size() * block_groups_atom_list[n][1].size() / nthreads;
244 :
245 530 : for (int tid = 0; tid < nthreads; ++tid) {
246 424 : thread_heaps[tid].reserve(estimated_pairs_per_thread);
247 : }
248 :
249 106 : heap.clear();
250 : #ifdef _OPENMP
251 106 : #pragma omp parallel for schedule(dynamic)
252 : #endif
253 : for (int i = 0; i < static_cast<int>(block_groups_atom_list[n][0].size()); i++) {
254 : AtomNumber ind0 = block_groups_atom_list[n][0][i];
255 : Vector Pos0 = isFirstBuild[n] ? refpdb.getPosition(ind0) : getPosition(atom_ind_hashmap[ind0.index()]);
256 : #ifdef _OPENMP
257 : int tid = omp_get_thread_num();
258 : #else
259 : int tid = 0;
260 : #endif
261 : auto& local_heap = thread_heaps[tid];
262 : std::set<AtomPair> local_unique_pairs;
263 :
264 :
265 : for (int j = 0; j < static_cast<int>(block_groups_atom_list[n][1].size()); j++) {
266 : AtomNumber ind1 = block_groups_atom_list[n][1][j];
267 : if (ind1 == ind0) {
268 : continue;
269 : }
270 :
271 : AtomPair test_pair{ind0, ind1};
272 : AtomPair reverse_pair{ind1, ind0};
273 :
274 : if (std::find(Exclude_Pairs[n].begin(), Exclude_Pairs[n].end(), test_pair) != Exclude_Pairs[n].end() ||
275 : std::find(Exclude_Pairs[n].begin(), Exclude_Pairs[n].end(), reverse_pair) != Exclude_Pairs[n].end()) {
276 : continue;
277 : }
278 : if (local_unique_pairs.count(reverse_pair)) {
279 : continue;
280 : }
281 : local_unique_pairs.insert(test_pair);
282 :
283 : Vector Pos1 = isFirstBuild[n] ? refpdb.getPosition(ind1) : getPosition(atom_ind_hashmap[ind1.index()]);
284 :
285 : double mag = pbcDistance(Pos0, Pos1).modulo();
286 : local_heap.emplace_back(mag, test_pair);
287 : }
288 : }
289 :
290 : // Flatten all thread-local heaps into one vector
291 106 : std::vector<std::pair<double, AtomPair>> all_pairs;
292 530 : for (auto& th : thread_heaps) {
293 424 : all_pairs.insert(all_pairs.end(), th.begin(), th.end());
294 : }
295 :
296 : // Build max-heap from top-K shortest pairs
297 106 : std::sort(all_pairs.begin(), all_pairs.end(), MinCompareDist());
298 106 : size_t K = std::min(static_cast<size_t>(tot_num_pairs[n]), all_pairs.size());
299 106 : heap.assign(all_pairs.begin(), all_pairs.begin() + K);
300 :
301 106 : std::sort(heap.begin(), heap.end(), MaxCompareDist());
302 : int chk1 = 0;
303 106 : int chk2 = Buffer_Pairs[n];
304 106 : delta_pd[n] = heap[chk1].first - heap[chk2].first;
305 106 : r_tolerance[n] = delta_pd[n]/4;
306 106 : listreduced[n].clear();
307 : std::set<AtomNumber, AtomNumberLess> uniqueAtoms;
308 1824 : for (const auto& pair : heap) {
309 1718 : uniqueAtoms.insert(pair.second.alpha_group); // Atom1
310 1718 : uniqueAtoms.insert(pair.second.beta_group); // Atom2
311 : }
312 1930 : for (const auto& atom : uniqueAtoms) {
313 1824 : listreduced[n].push_back(atom);
314 : }
315 106 : if(driver_mode || isFirstBuild[n]) {
316 3 : listreduced[n] = listall[n];
317 : }
318 106 : PL_atoms_ref_coords[n].resize(listreduced[n].size());
319 3421 : for (int i=0; i<static_cast<int>(listreduced[n].size()); i++) {
320 3315 : PL_atoms_ref_coords[n][i].zero();
321 3315 : PL_atoms_ref_coords[n][i] = isFirstBuild[n] ? refpdb.getPosition(listreduced[n][i]) : getPosition(atom_ind_hashmap[listreduced[n][i].index()]);
322 : }
323 : isFirstBuild[n] = false;
324 :
325 106 : }
326 :
327 303 : void PINES::updateBlockPairList(int n, std::vector<std::pair<double, AtomPair>>& heap) {
328 : // This is separate from staleness, which triggers a refresh instead of an update.
329 : // This is simply to update the pair distances in the MaxHeap/Pairlist with the new atom positions
330 : // and rearrange the ordering in the MaxHeap/Pairlisat
331 :
332 303 : plumed_massert((int)heap.size() >= tot_num_pairs[n], "heap too small");
333 303 : plumed_massert(!listreducedall_vec.empty(), "positions requested before requestAtoms");
334 :
335 : #ifdef _OPENMP
336 303 : #pragma omp parallel for
337 : #endif
338 : for (int i = 0; i < tot_num_pairs[n]; i++) {
339 : AtomNumber ind0 = heap[i].second.alpha_group;
340 : AtomNumber ind1 = heap[i].second.beta_group;
341 : auto it0 = atom_ind_hashmap.find(ind0.index());
342 : auto it1 = atom_ind_hashmap.find(ind1.index());
343 : plumed_massert(it0 != atom_ind_hashmap.end() && it1 != atom_ind_hashmap.end(),"updateBlockPairList: atom index not requested");
344 : heap[i].first = calculateDistance(n,ind0,ind1,mypdb);
345 : }
346 303 : std::nth_element(heap.begin(), heap.begin() + block_lengths[n], heap.end(), MinCompareDist());
347 303 : std::sort(heap.begin(), heap.begin() + block_lengths[n], MaxCompareDist());
348 303 : }
349 :
350 3434 : double PINES::calculateDistance(int n, const AtomNumber& ind0, const AtomNumber& ind1, const PDB& refpdb) {
351 3434 : Vector Pos0 = isFirstBuild[n] ? refpdb.getPosition(ind0) : getPosition(atom_ind_hashmap[ind0.index()]);
352 3434 : Vector Pos1 = isFirstBuild[n] ? refpdb.getPosition(ind1) : getPosition(atom_ind_hashmap[ind1.index()]);
353 3434 : return pbcDistance(Pos0, Pos1).modulo();
354 : }
355 :
356 1 : void PINES::resizeAllContainers(int N) {
357 : // Outer containers
358 1 : nstride.assign(N,1);
359 1 : steps_since_update.assign(N, 0);
360 1 : block_params.resize(N);
361 1 : block_groups_atom_list.assign(N, std::array<std::vector<PLMD::AtomNumber>, 2> { { {}, {} } });
362 1 : block_lengths.assign(N, 0);
363 1 : Buffer_Pairs.assign(N, 0);
364 1 : tot_num_pairs.assign(N, 0);
365 1 : Exclude_Pairs.resize(N);
366 1 : all_g1g2_pairs.assign(N,false);
367 1 : vecMaxHeapVecs.resize(N);
368 1 : PIV.assign(N, std::vector<double> {});
369 1 : listall.resize(N);
370 1 : listreduced.resize(N);
371 1 : stale_tolerance.assign(N,false);
372 1 : r00.assign(N, 0.0);
373 1 : sw.resize(N);
374 1 : sfs.resize(N);
375 1 : delta_pd.assign(N, 0.0);
376 1 : r_tolerance.assign(N, 0.0);
377 1 : PL_atoms_ref_coords.resize(N);
378 1 : input_filters.assign(N, std::array<std::array<bool,3>, 2> { { {false,false,false}, {false,false,false} } });
379 1 : ID_list.assign(N, std::array<std::vector<PLMD::AtomNumber>, 2> { { {}, {} } });
380 1 : ResID_list.assign(N, std::array<std::vector<int>, 2> { { {}, {} } });
381 2 : Name_list.assign(N, std::array<std::vector<std::string>, 2> { { {}, {} } });
382 : atom_ind_hashmap.clear();
383 1 : latched_pairs.resize(N);
384 1 : isFirstBuild.assign(N,true);
385 :
386 : // Per-block inner structures
387 4 : for (int n = 0; n < N; n++) {
388 3 : PL_atoms_ref_coords[n].resize(0); // will be filled per atom if needed
389 : }
390 :
391 : // Final 3D input_filters init
392 4 : for (int n = 0; n < N; n++) {
393 9 : for (int g = 0; g < 2; g++) {
394 6 : input_filters[n][g].fill(false); // [ID, ResID, Name]
395 : }
396 : }
397 1 : }
398 :
399 : // member state
400 : // static long inited_step = -1;
401 : // static std::vector<char> preupdated_block; // size N_Blocks, resets each step
402 :
403 606 : bool PINES::ensureBlockUpdated(int n) {
404 : const bool do_update = true; //(plumed.getStep() == 0 || driver_mode);
405 : if (!do_update) {
406 : return false; // not an update step
407 : }
408 606 : if ((int)vecMaxHeapVecs[n].size() != tot_num_pairs[n]) {
409 0 : buildMaxHeapVecBlock(n, mypdb, vecMaxHeapVecs[n]); // safe: PDB-based build
410 : }
411 606 : if (preupdated_block[n]) {
412 : return false; // already updated this step
413 : }
414 303 : updateBlockPairList(n, vecMaxHeapVecs[n]); // uses getPosition(); safe in calculate()
415 303 : preupdated_block[n] = 1;
416 303 : return true;
417 : }
418 :
419 0 : void PINES::latchFromCurrentHeaps() {
420 0 : if ((int)latched_pairs.size() != N_Blocks) {
421 0 : latched_pairs.resize(N_Blocks);
422 : }
423 0 : for (int n = 0; n < N_Blocks; ++n) {
424 0 : latched_pairs[n].resize(block_lengths[n]);
425 : const bool do_update = true; // (plumed.getStep() == 0 || steps_since_update[n] != 1 || driver_mode);
426 : if (do_update) {
427 : // top-K prefix (already sorted by ensureBlockUpdated if it ran)
428 0 : for (int i = 0; i < block_lengths[n]; ++i) {
429 0 : latched_pairs[n][i] = vecMaxHeapVecs[n][i].second;
430 : }
431 : } else {
432 : // buffered slice
433 : for (int i = 0; i < block_lengths[n]; ++i) {
434 : latched_pairs[n][i] = vecMaxHeapVecs[n][i + Buffer_Pairs[n]].second;
435 : }
436 : }
437 : }
438 0 : }
439 :
440 3737 : void PINES::logMsg(const std::string& msg, const std::string& section) {
441 3737 : log << "[" << plumed.getStep() << "] "
442 7474 : << "[" << section << "] " << msg << std::endl;
443 3737 : }
444 :
445 0 : void PINES::logMsg(const Vector& vec, const std::string& section) {
446 0 : log << "[" << plumed.getStep() << "] "
447 0 : << "[" << section << "] Vector = ("
448 0 : << vec[0] << ", " << vec[1] << ", " << vec[2] << ")" << std::endl;
449 0 : }
450 :
451 1 : PINES::PINES(const ActionOptions &ao) : PLUMED_COLVAR_INIT(ao),
452 1 : N_Blocks(),
453 1 : total_PIV_length(),
454 1 : inited_step(-1),
455 1 : driver_mode(),
456 1 : freeze_selection(),
457 : steps_since_update(),
458 : nstride(),
459 : ref_file(),
460 : sfs(),
461 : sw(),
462 : r00(),
463 : PIV(),
464 : ann_deriv(),
465 : listall(),
466 : listreduced(),
467 : listreducedall(),
468 : listreducedall_vec(),
469 1 : atom_ind_hashmap(),
470 : stale_tolerance(),
471 1 : mypdb(),
472 : block_params(),
473 : block_groups_atom_list(),
474 : block_lengths(),
475 : Buffer_Pairs(),
476 : tot_num_pairs(),
477 : input_filters(),
478 : delta_pd(),
479 : r_tolerance(),
480 : PL_atoms_ref_coords(),
481 : Exclude_Pairs(),
482 : Name_list(),
483 : ID_list(),
484 : ResID_list(),
485 : all_g1g2_pairs(),
486 : vecMaxHeapVecs(),
487 : latched_pairs(),
488 : preupdated_block(),
489 1 : isFirstBuild() {
490 1 : log.open("pines_debug.log", std::ios::out);
491 1 : if (!log.is_open()) {
492 0 : error("Failed to open debug log file.");
493 : }
494 :
495 : // Reference PDB file from which atom names, types, ids, and initial positions are determined
496 2 : parse("REF_FILE", ref_file);
497 1 : FILE *fp = fopen(ref_file.c_str(), "r");
498 1 : if (fp != NULL) {
499 1 : mypdb.readFromFilepointer(fp, usingNaturalUnits(), 0.1 / getUnits().getLength());
500 1 : fclose(fp);
501 : } else {
502 0 : error("Error in reference PDB file");
503 : }
504 :
505 : // Create variable to get number of blocks
506 1 : parse("N_BLOCKS", N_Blocks);
507 1 : parseFlag("DRIVERMODE", driver_mode);
508 1 : parseFlag("FREEZE_SELECTION", freeze_selection);
509 1 : resizeAllContainers(N_Blocks);
510 :
511 :
512 : // Check that the correct number of Blocks are specified
513 4 : for (int n = 0; n < N_Blocks; n++) {
514 6 : if (!parseNumbered("BLOCK", n + 1, block_params[n])) {
515 : break;
516 : }
517 : }
518 :
519 : // Parse blocks for keywords
520 4 : for (int n = 0; n < N_Blocks; n++) {
521 : std::string block_length;
522 3 : std::vector<std::string> ex_pairs_n;
523 3 : std::vector<std::string> g1_ids;
524 3 : std::vector<std::string> g2_ids;
525 3 : std::vector<std::string> g1_resids;
526 3 : std::vector<std::string> g2_resids;
527 3 : std::vector<std::string> g1_names;
528 3 : std::vector<std::string> g2_names;
529 :
530 3 : std::vector<std::string> block_data = Tools::getWords(block_params[n]);
531 :
532 3 : std::vector<std::string> G1_data;
533 3 : Tools::parseVector(block_data, "G1", G1_data, -1);
534 :
535 3 : std::vector<std::string> G2_data;
536 3 : Tools::parseVector(block_data, "G2", G2_data, -1);
537 :
538 3 : Tools::parseVector(G1_data, "ATOMID", g1_ids, -1);
539 :
540 3 : for (int i = 0; i < static_cast<int>(g1_ids.size()); i++) {
541 : AtomNumber g1i;
542 0 : g1i.setIndex(std::stoi(g1_ids[i]));
543 0 : ID_list[n][0].push_back(g1i);
544 : }
545 :
546 3 : Tools::parseVector(G2_data, "ATOMID", g2_ids, -1);
547 :
548 3 : for (int i = 0; i < static_cast<int>(g2_ids.size()); i++) {
549 : AtomNumber g2i;
550 0 : g2i.setIndex(std::stoi(g2_ids[i]));
551 0 : ID_list[n][1].push_back(g2i);
552 : }
553 :
554 3 : Tools::parseVector(G1_data, "RESID", g1_resids, -1);
555 3 : for (int i = 0; i < static_cast<int>(g1_resids.size()); i++) {
556 0 : ResID_list[n][0].push_back(std::stoi(g1_resids[i]));
557 : }
558 3 : Tools::parseVector(G2_data, "RESID", g2_resids, -1);
559 3 : for (int i = 0; i < static_cast<int>(g2_resids.size()); i++) {
560 0 : ResID_list[n][1].push_back(std::stoi(g2_resids[i]));
561 : }
562 3 : Tools::parseVector(G1_data, "NAME", g1_names, -1);
563 6 : for (int i = 0; i < static_cast<int>(g1_names.size()); i++) {
564 3 : Name_list[n][0].push_back(g1_names[i]);
565 : }
566 3 : Tools::parseVector(G2_data, "NAME", g2_names, -1);
567 7 : for (int i = 0; i < static_cast<int>(g2_names.size()); i++) {
568 4 : Name_list[n][1].push_back(g2_names[i]);
569 : }
570 :
571 6 : Tools::parse(block_data, "SIZE", block_length);
572 3 : block_lengths[n] = std::stoi(block_length);
573 :
574 6 : Tools::parseVector(block_data, "EXCLUDE_PAIRS", ex_pairs_n, -1);
575 3 : if (!ex_pairs_n.empty()) {
576 0 : for (int i = 0; i < static_cast<int>(ex_pairs_n.size() - 1); i+=2) {
577 : AtomNumber atom1;
578 0 : atom1.setIndex(std::stoi(ex_pairs_n[i]));
579 : AtomNumber atom2;
580 0 : atom2.setIndex(std::stoi(ex_pairs_n[i+1]));
581 :
582 : AtomPair excluded_pair;
583 0 : excluded_pair = {atom1, atom2};
584 0 : Exclude_Pairs[n].push_back(excluded_pair);
585 : }
586 : }
587 : std::string buffer_pairs;
588 : std::string pl_refresh;
589 3 : Tools::parse(block_data, "BUFFER", buffer_pairs);
590 6 : Tools::parse(block_data, "PL_REFRESH", pl_refresh);
591 :
592 3 : if (!buffer_pairs.empty()) {
593 3 : Buffer_Pairs[n] = std::stoi(buffer_pairs);
594 : } else {
595 0 : Buffer_Pairs[n] = 0;
596 : }
597 :
598 3 : if (!pl_refresh.empty()) {
599 3 : nstride[n] = std::stoi(pl_refresh);
600 : } else {
601 0 : nstride[n] = 25;
602 : }
603 3 : tot_num_pairs[n] = block_lengths[n] + Buffer_Pairs[n];
604 3 : }
605 :
606 1 : total_PIV_length = std::accumulate(block_lengths.begin(), block_lengths.end(), 0);
607 4 : for (int n = 0; n < N_Blocks; n++) {
608 : // pseudo-code ish
609 9 : for (int g = 0; g < 2; g++) {
610 6 : if (ID_list[n][g].size() > 0) {
611 0 : input_filters[n][g][0] = true;
612 : }
613 6 : if (ResID_list[n][g].size() > 0) {
614 0 : input_filters[n][g][1] = true;
615 : }
616 6 : if (Name_list[n][g].size() > 0) {
617 6 : input_filters[n][g][2] = true;
618 : }
619 : }
620 : }
621 :
622 4 : for (int n=0; n < N_Blocks; n++) {
623 3 : listall[n].clear();
624 : }
625 :
626 1527 : for (int i=0; i < static_cast<int>(mypdb.getAtomNumbers().size()); i++) {
627 1526 : AtomNumber ind = mypdb.getAtomNumbers()[i];
628 1526 : int resid = mypdb.getResidueNumber(ind);
629 1526 : std::string atom_name = mypdb.getAtomName(ind);
630 6104 : for (int n = 0; n < N_Blocks; n++) {
631 : bool atom_added = false;
632 13734 : for (int g = 0; g < 2; g++) {
633 9156 : if (atomMatchesFilters(n, g, ind, resid, atom_name)) {
634 1528 : block_groups_atom_list[n][g].push_back(ind);
635 : atom_added = true;
636 : }
637 : }
638 4578 : if (atom_added) {
639 1528 : listall[n].push_back(ind);
640 : }
641 : }
642 : }
643 :
644 : int g1g2_pairs;
645 4 : for (int n = 0; n < N_Blocks; n++) {
646 3 : g1g2_pairs = block_groups_atom_list[n][0].size() * block_groups_atom_list[n][1].size() - Exclude_Pairs[n].size();
647 3 : if (g1g2_pairs == tot_num_pairs[n]) {
648 1 : all_g1g2_pairs[n] = 1;
649 : } else {
650 2 : all_g1g2_pairs[n] = 0;
651 : }
652 : }
653 :
654 4 : for (int n = 0; n < N_Blocks; n++)
655 6 : if (!parseNumbered("SWITCH", n + 1, sw[n])) {
656 : break;
657 : }
658 :
659 : std::string errors;
660 4 : for (int n = 0; n < N_Blocks; n++) {
661 3 : sfs[n].set(sw[n], errors);
662 : std::string num;
663 3 : Tools::convert(n + 1, num);
664 3 : if (errors.length() != 0) {
665 0 : error("problem reading SWITCH" + num + " keyword : " + errors);
666 : }
667 3 : r00[n] = sfs[n].get_r0();
668 : }
669 1 : checkRead();
670 :
671 : int total_count = 0;
672 4 : for (int n = 0; n < N_Blocks; n++) {
673 7 : for (int i = 0; i < block_lengths[n]; i++) {
674 4 : std::string comp = "ELEMENT-" + std::to_string(total_count);
675 4 : addComponentWithDerivatives(comp);
676 4 : componentIsNotPeriodic(comp);
677 4 : total_count += 1;
678 : }
679 : }
680 1 : }
681 :
682 101 : void PINES::prepare() {
683 : bool heap_refreshed = false;
684 : bool stale_refresh_prep = false;
685 404 : for (int n = 0; n < N_Blocks; n++) {
686 303 : if (steps_since_update[n] == 0 || plumed.getStep() == 0 || driver_mode || (int)vecMaxHeapVecs[n].size() != tot_num_pairs[n]) {
687 103 : buildMaxHeapVecBlock(n, mypdb, vecMaxHeapVecs[n]);
688 : heap_refreshed = true;
689 200 : } else if((steps_since_update[n] >= nstride[n] || stale_tolerance[n]) && !all_g1g2_pairs[n]) {
690 100 : listreduced[n] = listall[n];
691 : stale_refresh_prep = true;
692 100 : steps_since_update[n] = -1;
693 : stale_tolerance[n] = false;
694 : }
695 303 : steps_since_update[n]+=1;
696 : }
697 101 : if (heap_refreshed || stale_refresh_prep) {
698 : // Collate all atoms from all lists
699 : listreducedall.clear();
700 404 : for (const auto& blockList : listreduced) { // Loop over each block's vector
701 303 : listreducedall.insert(blockList.begin(), blockList.end()); // Insert all atoms into set
702 : }
703 : atom_ind_hashmap.clear();
704 101 : listreducedall_vec = std::vector<AtomNumber>(listreducedall.begin(), listreducedall.end());
705 79677 : for (int i=0; i < static_cast<int>(listreducedall_vec.size()); i++) {
706 79576 : atom_ind_hashmap[listreducedall_vec[i].index()] = i;
707 : }
708 101 : requestAtoms(listreducedall_vec);
709 101 : ann_deriv.resize(listreducedall_vec.size());
710 :
711 79677 : for (int i=0; i < static_cast<int>(ann_deriv.size()); i++) {
712 79576 : ann_deriv[i].resize(total_PIV_length);
713 : }
714 :
715 : }
716 101 : }
717 :
718 101 : void PINES::calculate() {
719 101 : const long step = plumed.getStep();
720 101 : if (inited_step == -1) {
721 4 : for (int n =0; n < N_Blocks; n++) {
722 3 : buildMaxHeapVecBlock(n, mypdb, vecMaxHeapVecs[n]);
723 : }
724 : }
725 101 : if (step != inited_step) {
726 101 : if ((int)preupdated_block.size() != N_Blocks) {
727 1 : preupdated_block.assign(N_Blocks, 0);
728 : } else {
729 : std::fill(preupdated_block.begin(), preupdated_block.end(), 0);
730 : }
731 :
732 : // One pass: ensure any block that should update does so now (once per step)
733 404 : for (int n = 0; n < N_Blocks; n++) {
734 303 : ensureBlockUpdated(n);
735 : }
736 :
737 : // Optional: latch here (membership/order frozen for this step)
738 101 : if (freeze_selection) {
739 0 : latchFromCurrentHeaps();
740 : }
741 :
742 404 : for (int n = 0; n < N_Blocks; n++) {
743 909 : logMsg("_________________N: "+std::to_string(n)+"____________________\n","Calculate");
744 3737 : for (int p = 0; p < static_cast<int>(vecMaxHeapVecs[n].size()); p++) {
745 10302 : logMsg("P: "+std::to_string(vecMaxHeapVecs[n][p].first)+"\n","Calculate");
746 : }
747 : }
748 101 : inited_step = step;
749 : }
750 :
751 101 : Vector ref_xyz, step_xyz;
752 : AtomNumber aID;
753 : double delta_r;
754 :
755 :
756 404 : for(int n = 0; n < N_Blocks; n++) {
757 303 : if(steps_since_update[n] > 0 && steps_since_update[n] < nstride[n] && r_tolerance[n] > 0.0 && !all_g1g2_pairs[n]) {
758 0 : for(int i = 0; i < static_cast<int>(listreduced[n].size()); i++) {
759 0 : aID = listreduced[n][i];
760 0 : ref_xyz = PL_atoms_ref_coords[n][i];
761 0 : step_xyz = getPosition(atom_ind_hashmap[aID.index()]);
762 0 : delta_r = pbcDistance(ref_xyz, step_xyz).modulo();
763 0 : if(delta_r >= r_tolerance[n]) {
764 : stale_tolerance[n] = true;
765 : }
766 : }
767 : }
768 : }
769 101 : const Vector zeroVec(0., 0., 0.);
770 : #ifdef _OPENMP
771 101 : #pragma omp parallel for schedule(static)
772 : #endif
773 : for (unsigned j = 0; j < ann_deriv.size(); j++) {
774 : for (unsigned i = 0; i < ann_deriv[j].size(); i++) {
775 : ann_deriv[j][i] = zeroVec;
776 : }
777 : }
778 :
779 : #ifdef _OPENMP
780 101 : #pragma omp parallel for schedule(static)
781 : #endif
782 : for (int n = 0; n < N_Blocks; n++) {
783 : PIV[n].resize(block_lengths[n]);
784 : for (int i = 0; i < block_lengths[n]; i++) {
785 : PIV[n][i] = 0.;
786 : }
787 : }
788 :
789 :
790 :
791 :
792 :
793 : int PINES_element = 0;
794 101 : if (freeze_selection) {
795 0 : for (int n = 0; n < N_Blocks; ++n) {
796 0 : logMsg("***************N: "+std::to_string(n)+"***************\n","bigCalc");
797 0 : for (int i = 0; i < block_lengths[n]; ++i) {
798 0 : const auto& pr = latched_pairs[n][i];
799 0 : int a0 = atom_ind_hashmap[pr.alpha_group.index()];
800 0 : int a1 = atom_ind_hashmap[pr.beta_group.index()];
801 0 : Vector r0 = getPosition(a0);
802 0 : Vector r1 = getPosition(a1);
803 0 : double mod_dist = pbcDistance(r0, r1).modulo();
804 0 : Vector dr = pbcDistance(r0, r1) / mod_dist;
805 0 : logMsg("P: "+std::to_string(vecMaxHeapVecs[n][i].first)+"\n","bigCalc");
806 :
807 0 : double dfunc = 0.0;
808 0 : PIV[n][i] = sfs[n].calculate(mod_dist, dfunc);
809 0 : double ds = dfunc * mod_dist;
810 0 : ann_deriv[a0][PINES_element] = -ds * dr;
811 0 : ann_deriv[a1][PINES_element] = ds * dr;
812 0 : PINES_element += 1;
813 : }
814 : }
815 : } else {
816 404 : for (int n = 0; n < N_Blocks; n++) {
817 :
818 303 : ensureBlockUpdated(n);
819 : bool updated = true;
820 707 : for (int i=0; i < block_lengths[n]; i++) {
821 : double ds_element = 0.;
822 404 : double dfunc = 0.;
823 : int local_aid0, local_aid1;
824 : int piv_ind;
825 : if (updated) {
826 : piv_ind = i;
827 : } else {
828 : piv_ind = i + Buffer_Pairs[n];
829 : }
830 404 : PIV[n][i] = sfs[n].calculate(vecMaxHeapVecs[n][piv_ind].first, dfunc);
831 404 : local_aid0 = atom_ind_hashmap[vecMaxHeapVecs[n][piv_ind].second.alpha_group.index()];
832 404 : local_aid1 = atom_ind_hashmap[vecMaxHeapVecs[n][piv_ind].second.beta_group.index()];
833 404 : Vector Pos0 = getPosition(local_aid0);
834 404 : Vector Pos1 = getPosition(local_aid1);
835 404 : Vector dr_dcoord = pbcDistance(Pos0, Pos1) / vecMaxHeapVecs[n][piv_ind].first;
836 404 : ds_element = dfunc * vecMaxHeapVecs[n][piv_ind].first;
837 404 : ann_deriv[local_aid0][PINES_element] = -ds_element * dr_dcoord;
838 404 : ann_deriv[local_aid1][PINES_element] = ds_element * dr_dcoord;
839 404 : PINES_element += 1;
840 : }
841 : }
842 : }
843 :
844 :
845 101 : if (comm.initialized()) {
846 : // Flatten PIV
847 0 : std::vector<double> flat_PIV;
848 0 : std::vector<int> piv_sizes;
849 0 : for (const auto& vec : PIV) {
850 0 : piv_sizes.push_back(vec.size());
851 0 : flat_PIV.insert(flat_PIV.end(), vec.begin(), vec.end());
852 : }
853 :
854 : // Flatten ann_deriv
855 0 : std::vector<double> flat_deriv;
856 0 : std::vector<std::pair<int, int>> deriv_shapes;
857 0 : for (const auto& atom_deriv : ann_deriv) {
858 0 : deriv_shapes.emplace_back(atom_deriv.size(), 3);
859 0 : for (const auto& val : atom_deriv) {
860 0 : flat_deriv.insert(flat_deriv.end(), val[0]);
861 0 : flat_deriv.insert(flat_deriv.end(), val[1]);
862 0 : flat_deriv.insert(flat_deriv.end(), val[2]);
863 : }
864 : }
865 :
866 : // Perform global summation ONCE per array
867 0 : comm.Sum(flat_PIV);
868 0 : comm.Sum(flat_deriv);
869 :
870 : // Unflatten PIV
871 : {
872 : size_t idx = 0;
873 0 : for (int j = 0; j < N_Blocks; ++j) {
874 0 : for (int i = 0; i < piv_sizes[j]; ++i) {
875 0 : PIV[j][i] = flat_PIV[idx++];
876 : }
877 : }
878 : }
879 :
880 : // Unflatten ann_deriv
881 : {
882 : size_t idx = 0;
883 0 : for (size_t i = 0; i < ann_deriv.size(); ++i) {
884 0 : int ncols = deriv_shapes[i].first;
885 0 : for (int j = 0; j < ncols; ++j) {
886 0 : ann_deriv[i][j][0] = flat_deriv[idx++];
887 0 : ann_deriv[i][j][1] = flat_deriv[idx++];
888 0 : ann_deriv[i][j][2] = flat_deriv[idx++];
889 : }
890 : }
891 : }
892 0 : }
893 :
894 :
895 :
896 :
897 :
898 : //timer.start("valuePass");
899 : // Pass values and derivates to next stage
900 : unsigned total_count = 0;
901 404 : for (int j = 0; j < N_Blocks; j++) {
902 707 : for (int i = 0; i < block_lengths[j]; i++) {
903 404 : std::string comp = "ELEMENT-" + std::to_string(total_count);
904 404 : Value *valueNew = getPntrToComponent(comp);
905 404 : valueNew->set(PIV[j][i]);
906 318708 : for (unsigned k = 0; k < ann_deriv.size(); k++) {
907 318304 : setAtomsDerivatives(valueNew, k, ann_deriv[k][total_count]);
908 : }
909 404 : total_count += 1;
910 : }
911 : }
912 :
913 101 : const int Natoms = static_cast<int>(ann_deriv.size());
914 :
915 : // Precompute current box and fractional (scaled) coords once per step
916 101 : Tensor box = getPbc().getBox();
917 101 : std::vector<Vector> scaled(Natoms);
918 79677 : for (int a = 0; a < Natoms; ++a) {
919 79576 : scaled[a] = getPbc().realToScaled(getPosition(a));
920 : }
921 :
922 : // Build one virial tensor per component (column)
923 101 : std::vector<Tensor> virials(total_PIV_length);
924 505 : for (int col = 0; col < total_PIV_length; ++col) {
925 404 : Tensor dSdB;
926 : dSdB.zero(); // ∂S/∂B at fixed scaled coords
927 :
928 : // dSdB = sum_i s_i ⊗ g_i with s_i = scaled (fractional), g_i = ∂S/∂r_i
929 318708 : for (int a = 0; a < Natoms; ++a) {
930 318304 : const Vector& gi = ann_deriv[a][col]; // atomic gradient for this component
931 318304 : if (gi[0]==0.0 && gi[1]==0.0 && gi[2]==0.0) {
932 317496 : continue;
933 : }
934 : const Vector& si = scaled[a];
935 :
936 : // Outer product: (∂S/∂B)_{ik} += s_i[i] * g_i[k]
937 808 : dSdB(0,0) += si[0]*gi[0];
938 808 : dSdB(0,1) += si[0]*gi[1];
939 808 : dSdB(0,2) += si[0]*gi[2];
940 808 : dSdB(1,0) += si[1]*gi[0];
941 808 : dSdB(1,1) += si[1]*gi[1];
942 808 : dSdB(1,2) += si[1]*gi[2];
943 808 : dSdB(2,0) += si[2]*gi[0];
944 808 : dSdB(2,1) += si[2]*gi[1];
945 808 : dSdB(2,2) += si[2]*gi[2];
946 : }
947 :
948 : // PLUMED numerical uses: virial = - B^T * (∂S/∂B)
949 404 : virials[col] = - matmul(box.transpose(), dSdB);
950 : }
951 :
952 505 : for (int col = 0; col < total_PIV_length; ++col) {
953 : // If you already have the Value* via getPntrToComponent("ELEMENT-<col>"), reuse it.
954 : // Otherwise copyOutput(col) is fine here.
955 404 : Value* v = copyOutput(col);
956 :
957 1616 : for (int i = 0; i < 3; ++i) { // i = column
958 4848 : for (int k = 0; k < 3; ++k) { // k = row
959 3636 : v->addDerivative(3*Natoms + 3*k + i, virials[col](k,i));
960 : }
961 : }
962 : }
963 101 : }
964 : }
965 : }
|
