@@ -54,6 +54,11 @@ struct AtomNumberLess {
5454 }
5555};
5656
57+ struct AtomPair {
58+ AtomNumber alpha_group;
59+ AtomNumber beta_group;
60+ };
61+
5762class PINES : public PLMD ::colvar::Colvar {
5863private:
5964 PLMD::Stopwatch timer;
@@ -81,27 +86,27 @@ class PINES : public PLMD::colvar::Colvar {
8186 std::vector<bool > stale_tolerance;
8287 PDB mypdb;
8388 std::vector<std::string> block_params;
84- std::vector<std::vector <std::vector<AtomNumber> > > block_groups_atom_list;
89+ std::vector<std::array <std::vector<AtomNumber>, 2 > > block_groups_atom_list;
8590 std::vector<int > block_lengths;
8691 std::vector<int > Buffer_Pairs;
8792 std::vector<int > tot_num_pairs;
88- std::vector<std::vector <std::vector <bool > > > input_filters;
93+ std::vector<std::array <std::array <bool , 3 >, 2 > > input_filters;
8994 std::vector<double > delta_pd;
9095 std::vector<double > r_tolerance;
9196 std::vector<std::vector<Vector> > PL_atoms_ref_coords;
92- std::vector<std::vector<std::pair<AtomNumber,AtomNumber> > > Exclude_Pairs;
93- std::vector<std::vector <std::vector<std::string> > > Name_list;
94- std::vector<std::vector <std::vector<AtomNumber> > > ID_list;
95- std::vector<std::vector <std::vector<int > > > ResID_list;
97+ std::vector<std::vector<AtomPair > > Exclude_Pairs;
98+ std::vector<std::array <std::vector<std::string> >, 2 > Name_list;
99+ std::vector<std::array <std::vector<AtomNumber> >, 2 > ID_list;
100+ std::vector<std::array <std::vector<int > >, 2 > ResID_list;
96101 std::vector<char > all_g1g2_pairs;
97- std::vector<std::vector<std::pair<double , std::pair<AtomNumber,AtomNumber> > > > vecMaxHeapVecs;
98- std::vector<std::vector<std::pair<AtomNumber,AtomNumber> > > latched_pairs;
99- std::vector<char > preupdated_block;
102+ std::vector<std::vector<std::pair<double , AtomPair > > > vecMaxHeapVecs;
103+ std::vector<std::vector<AtomPair > > latched_pairs;
104+ std::vector<int > preupdated_block;
100105 std::vector<bool > isFirstBuild;
101106
102107 bool atomMatchesFilters (int n, int g, AtomNumber ind, int resid, const std::string& atom_name);
103- void buildMaxHeapVecBlock (int n, const PDB& mypdb, std::vector<std::pair<double , std::pair<AtomNumber, AtomNumber> >>& heap);
104- void updateBlockPairList (int n, std::vector<std::pair<double , std::pair<AtomNumber, AtomNumber> >>& heap);
108+ void buildMaxHeapVecBlock (int n, const PDB& mypdb, std::vector<std::pair<double , AtomPair >>& heap);
109+ void updateBlockPairList (int n, std::vector<std::pair<double , AtomPair >>& heap);
105110 double calculateDistance (int n, const AtomNumber& ind0, const AtomNumber& ind1, const PDB& mypdb);
106111 std::ofstream log;
107112 bool ensureBlockUpdated (int n);
@@ -116,12 +121,12 @@ class PINES : public PLMD::colvar::Colvar {
116121 // ~PINES();
117122 // active methods:
118123 struct MinCompareDist {
119- bool operator ()(const std::pair<double , std::pair<AtomNumber, AtomNumber>> & p1, const std::pair<double , std::pair<AtomNumber, AtomNumber> >& p2) {
124+ bool operator ()(const std::pair<double , AtomPair> & p1, const std::pair<double , AtomPair >& p2) {
120125 return p1.first < p2.first ; // Min heap
121126 }
122127 };
123128 struct MaxCompareDist {
124- bool operator ()(const std::pair<double , std::pair<AtomNumber, AtomNumber>> & p1, const std::pair<double , std::pair<AtomNumber, AtomNumber> >& p2) {
129+ bool operator ()(const std::pair<double , AtomPair> & p1, const std::pair<double , AtomPair >& p2) {
125130 return p1.first > p2.first ; // Max heap
126131 }
127132 };
@@ -199,11 +204,11 @@ bool PINES::atomMatchesFilters(int n, int g, AtomNumber ind, int resid, const st
199204 return id_check && res_check && name_check;
200205}
201206
202- void PINES::buildMaxHeapVecBlock (int n, const PDB& mypdb, std::vector<std::pair<double , std::pair<AtomNumber, AtomNumber> >>& heap) {
207+ void PINES::buildMaxHeapVecBlock (int n, const PDB& mypdb, std::vector<std::pair<double , AtomPair >>& heap) {
203208 // logMsg("Building Pair List for block: " + std::to_string(n), "buildMaxHeapVecBlock");
204209
205210 int nthreads = omp_get_max_threads ();
206- std::vector<std::vector<std::pair<double , std::pair<AtomNumber, AtomNumber> >>> thread_heaps (nthreads);
211+ std::vector<std::vector<std::pair<double , AtomPair >>> thread_heaps (nthreads);
207212
208213 // logMsg("isFirstBuild[n]? " + std::to_string(isFirstBuild[n]), "buildMaxHeapVecBlock");
209214
@@ -221,7 +226,7 @@ void PINES::buildMaxHeapVecBlock(int n, const PDB& mypdb, std::vector<std::pair<
221226 Vector Pos0 = isFirstBuild[n] ? mypdb.getPosition (ind0) : getPosition (atom_ind_hashmap[ind0.index ()]);
222227 int tid = omp_get_thread_num ();
223228 auto & local_heap = thread_heaps[tid];
224- std::unordered_set<std::pair<AtomNumber, AtomNumber> , pair_hash> local_unique_pairs;
229+ std::unordered_set<AtomPair , pair_hash> local_unique_pairs;
225230 // logMsg("ind0: " + std::to_string(ind0.index()), "buildMaxHeapVecBlock");
226231 // logMsg("Num pairs in heap: " + std::to_string(block_groups_atom_list[n][0].size()), "buildMaxHeapVecBlock");
227232 // logMsg(Pos0, "buildMaxHeapVecBlock");
@@ -255,7 +260,7 @@ void PINES::buildMaxHeapVecBlock(int n, const PDB& mypdb, std::vector<std::pair<
255260 }
256261
257262 // Flatten all thread-local heaps into one vector
258- std::vector<std::pair<double , std::pair<AtomNumber, AtomNumber> >> all_pairs;
263+ std::vector<std::pair<double , AtomPair >> all_pairs;
259264 for (auto & th : thread_heaps) {
260265 all_pairs.insert (all_pairs.end (), th.begin (), th.end ());
261266 }
@@ -276,8 +281,8 @@ void PINES::buildMaxHeapVecBlock(int n, const PDB& mypdb, std::vector<std::pair<
276281 listreduced[n].clear ();
277282 std::set<AtomNumber, AtomNumberLess> uniqueAtoms;
278283 for (const auto & pair : heap) {
279- uniqueAtoms.insert (pair.second .first ); // Atom1
280- uniqueAtoms.insert (pair.second .second ); // Atom2
284+ uniqueAtoms.insert (pair.second .alpha_group ); // Atom1
285+ uniqueAtoms.insert (pair.second .beta_group ); // Atom2
281286 }
282287 for (const auto & atom : uniqueAtoms) {
283288 listreduced[n].push_back (atom);
@@ -299,7 +304,7 @@ void PINES::buildMaxHeapVecBlock(int n, const PDB& mypdb, std::vector<std::pair<
299304
300305}
301306
302- void PINES::updateBlockPairList (int n, std::vector<std::pair<double , std::pair<AtomNumber, AtomNumber> >>& heap) {
307+ void PINES::updateBlockPairList (int n, std::vector<std::pair<double , AtomPair >>& heap) {
303308 // This is separate from staleness, which triggers a refresh instead of an update.
304309 // This is simply to update the pair distances in the MaxHeap/Pairlist with the new atom positions
305310 // and rearrange the ordering in the MaxHeap/Pairlisat
@@ -309,8 +314,8 @@ void PINES::updateBlockPairList(int n, std::vector<std::pair<double, std::pair<A
309314
310315 #pragma omp parallel for
311316 for (int i = 0 ; i < tot_num_pairs[n]; i++) {
312- AtomNumber ind0 = heap[i].second .first ;
313- AtomNumber ind1 = heap[i].second .second ;
317+ AtomNumber ind0 = heap[i].second .alpha_group ;
318+ AtomNumber ind1 = heap[i].second .beta_group ;
314319 auto it0 = atom_ind_hashmap.find (ind0.index ());
315320 auto it1 = atom_ind_hashmap.find (ind1.index ());
316321 plumed_massert (it0 != atom_ind_hashmap.end () && it1 != atom_ind_hashmap.end ()," updateBlockPairList: atom index not requested" );
@@ -335,38 +340,33 @@ void PINES::resizeAllContainers(int N) {
335340 nstride.assign (N,1 );
336341 steps_since_update.assign (N, 0 );
337342 block_params.resize (N);
338- block_groups_atom_list.resize (N );
339- block_lengths.resize (N );
340- Buffer_Pairs.resize (N );
341- tot_num_pairs.resize (N );
343+ block_groups_atom_list.assign (N, { {}, {} } );
344+ block_lengths.assign (N, 0 );
345+ Buffer_Pairs.assign (N, 0 );
346+ tot_num_pairs.assign (N, 0 );
342347 Exclude_Pairs.resize (N);
343- all_g1g2_pairs.resize (N );
348+ all_g1g2_pairs.assign (N, false );
344349 vecMaxHeapVecs.resize (N);
345- PIV.resize (N );
350+ PIV.assign (N, 0.0 );
346351 listall.resize (N);
347352 listreduced.resize (N);
348- stale_tolerance.resize (N );
349- r00.resize (N );
353+ stale_tolerance.assign (N, false );
354+ r00.assign (N, 0.0 );
350355 sw.resize (N);
351356 sfs.resize (N);
352357 delta_pd.assign (N, 0.0 );
353358 r_tolerance.assign (N, 0.0 );
354359 PL_atoms_ref_coords.resize (N);
355- input_filters.resize (N );
356- ID_list.resize (N );
357- ResID_list.resize (N );
358- Name_list.resize (N );
360+ input_filters.assign (N, { {}, {} } );
361+ ID_list.assign (N, { {}, {} } );
362+ ResID_list.assign (N, { {}, {} } );
363+ Name_list.assign (N, { {}, {} } );
359364 atom_ind_hashmap.clear ();
360365 latched_pairs.resize (N);
361366 isFirstBuild.assign (N,true );
362367
363368 // Per-block inner structures
364369 for (int n = 0 ; n < N; n++) {
365- block_groups_atom_list[n].resize (2 );
366- ID_list[n].resize (2 );
367- ResID_list[n].resize (2 );
368- Name_list[n].resize (2 );
369- input_filters[n].resize (2 );
370370 PL_atoms_ref_coords[n].resize (0 ); // will be filled per atom if needed
371371 }
372372
@@ -431,53 +431,45 @@ void PINES::logMsg(const Vector& vec, const std::string& section) {
431431}
432432
433433PINES::PINES (const ActionOptions &ao) : PLUMED_COLVAR_INIT(ao),
434- N_Blocks (1 ),
435- total_PIV_length(1 ),
436- driver_mode(false ),
437- freeze_selection(false ),
438- last_step_latched(- 1 ),
434+ N_Blocks (),
435+ total_PIV_length(),
436+ driver_mode(),
437+ freeze_selection(),
438+ last_step_latched(),
439439 inited_step(-1 ),
440- steps_since_update(std::vector< int >( 1 ) ),
441- nstride(std::vector< int >( 1 , 10 ) ),
442- ref_file(std::string() ),
440+ steps_since_update(),
441+ nstride(),
442+ ref_file(),
443443 sfs(),
444- sw(std::vector<string>( 1 ) ),
445- r00(std::vector< double >( 1 ) ),
446- PIV(std::vector<std::vector< double >>( 1 ) ),
447- ann_deriv(std::vector<std::vector<Vector> >( 1 ) ),
448- listall(std::vector<std::vector<AtomNumber> >() ),
449- listreduced(std::vector<std::vector<AtomNumber> >() ),
450- listreducedall(std::set<AtomNumber, AtomNumberLess>() ),
451- listreducedall_vec(std::vector<AtomNumber>() ),
444+ sw(),
445+ r00(),
446+ PIV(),
447+ ann_deriv(),
448+ listall(),
449+ listreduced(),
450+ listreducedall(),
451+ listreducedall_vec(),
452452 atom_ind_hashmap(),
453- stale_tolerance(std::vector< bool >( 1 , false ) ),
453+ stale_tolerance(),
454454 mypdb(),
455- block_params(std::vector<string>(1 )),
456- block_groups_atom_list(std::vector<std::vector<std::vector<AtomNumber> > >()),
457- block_lengths(std::vector<int >(1 ,1 )),
458- Buffer_Pairs(std::vector<int >(1 ,0 )),
459- tot_num_pairs(std::vector<int >(1 ,1 )),
460- input_filters(
461- std::vector<std::vector<std::vector<bool > > >(
462- 1 , // outermost dimension
463- std::vector<std::vector<bool >>(
464- 1 , // middle dimension
465- std::vector<bool >(1 , false ) // innermost dimension with value `false`
466- )
467- )
468- ),
469- delta_pd(std::vector<double >()),
470- r_tolerance(std::vector<double >()),
455+ block_params(),
456+ block_groups_atom_list(),
457+ block_lengths(),
458+ Buffer_Pairs(),
459+ tot_num_pairs(),
460+ input_filters(),
461+ delta_pd(),
462+ r_tolerance(),
471463 PL_atoms_ref_coords(),
472- Exclude_Pairs(std::vector<std::vector<std::pair<AtomNumber, AtomNumber> > >() ),
473- Name_list(std::vector<std::vector<std::vector<string> > >() ),
474- ID_list(std::vector<std::vector<std::vector<AtomNumber> > >() ),
475- ResID_list(std::vector<std::vector<std::vector< int > > >() ),
476- all_g1g2_pairs(std::vector< char >() ),
464+ Exclude_Pairs(),
465+ Name_list(),
466+ ID_list(),
467+ ResID_list(),
468+ all_g1g2_pairs(),
477469 vecMaxHeapVecs(),
478- latched_pairs(std::vector<std::vector<std::pair<AtomNumber,AtomNumber> > >() ),
479- preupdated_block(std::vector< char >() ),
480- isFirstBuild(std::vector< bool >() ) {
470+ latched_pairs(),
471+ preupdated_block(),
472+ isFirstBuild() {
481473 log.open (" pines_debug.log" , std::ios::out);
482474 if (!log.is_open ()) {
483475 error (" Failed to open debug log file." );
@@ -572,7 +564,7 @@ PINES::PINES(const ActionOptions &ao) : PLUMED_COLVAR_INIT(ao),
572564 AtomNumber atom2;
573565 atom2.setIndex (std::stoi (ex_pairs_n[i+1 ]));
574566
575- std::pair<AtomNumber, AtomNumber> excluded_pair;
567+ AtomPair excluded_pair;
576568 excluded_pair = {atom1, atom2};
577569 Exclude_Pairs[n].push_back (excluded_pair);
578570 }
@@ -654,10 +646,6 @@ PINES::PINES(const ActionOptions &ao) : PLUMED_COLVAR_INIT(ao),
654646 }
655647 }
656648
657- r00.resize (N_Blocks);
658- sw.resize (N_Blocks);
659- sfs.resize (N_Blocks);
660-
661649 for (unsigned n = 0 ; n < N_Blocks; n++)
662650 if (!parseNumbered (" SWITCH" , n + 1 , sw[n])) {
663651 break ;
@@ -811,8 +799,8 @@ void PINES::calculate() {
811799 logMsg (" ***************N: " +std::to_string (n)+" ***************\n " ," bigCalc" );
812800 for (int i = 0 ; i < block_lengths[n]; ++i) {
813801 const auto & pr = latched_pairs[n][i];
814- int a0 = atom_ind_hashmap[pr.first .index ()];
815- int a1 = atom_ind_hashmap[pr.second .index ()];
802+ int a0 = atom_ind_hashmap[pr.alpha_group .index ()];
803+ int a1 = atom_ind_hashmap[pr.beta_group .index ()];
816804 Vector r0 = getPosition (a0);
817805 Vector r1 = getPosition (a1);
818806 double mod_dist = pbcDistance (r0, r1).modulo ();
@@ -852,8 +840,8 @@ void PINES::calculate() {
852840 // timer.stop("sfs");
853841 // logMsg("PIV Values: PIV[" + std::to_string(n) + "][" + std::to_string(i) + "]= " + std::to_string(PIV[n][i]), "Calculate");
854842 // timer.start("distCalc");
855- local_aid0 = atom_ind_hashmap[vecMaxHeapVecs[n][piv_ind].second .first .index ()];
856- local_aid1 = atom_ind_hashmap[vecMaxHeapVecs[n][piv_ind].second .second .index ()];
843+ local_aid0 = atom_ind_hashmap[vecMaxHeapVecs[n][piv_ind].second .alpha_group .index ()];
844+ local_aid1 = atom_ind_hashmap[vecMaxHeapVecs[n][piv_ind].second .beta_group .index ()];
857845 Vector Pos0 = getPosition (local_aid0);
858846 Vector Pos1 = getPosition (local_aid1);
859847 Vector dr_dcoord = pbcDistance (Pos0, Pos1) / vecMaxHeapVecs[n][piv_ind].first ;
0 commit comments