Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2015-2023 The plumed team
3 : (see the PEOPLE file at the root of the distribution for a list of names)
4 :
5 : See http://www.plumed.org for more information.
6 :
7 : This file is part of plumed, version 2.
8 :
9 : plumed is free software: you can redistribute it and/or modify
10 : it under the terms of the GNU Lesser General Public License as published by
11 : the Free Software Foundation, either version 3 of the License, or
12 : (at your option) any later version.
13 :
14 : plumed is distributed in the hope that it will be useful,
15 : but WITHOUT ANY WARRANTY; without even the implied warranty of
16 : MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17 : GNU Lesser General Public License for more details.
18 :
19 : You should have received a copy of the GNU Lesser General Public License
20 : along with plumed. If not, see <http://www.gnu.org/licenses/>.
21 : +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
22 : #include "AdjacencyMatrixBase.h"
23 : #include "multicolvar/BridgedMultiColvarFunction.h"
24 : #include "multicolvar/AtomValuePack.h"
25 : #include "multicolvar/CatomPack.h"
26 : #include "core/PlumedMain.h"
27 : #include "core/ActionSet.h"
28 :
29 : namespace PLMD {
30 : namespace adjmat {
31 :
32 29 : void AdjacencyMatrixBase::registerKeywords( Keywords& keys ) {
33 29 : multicolvar::MultiColvarBase::registerKeywords( keys );
34 58 : keys.remove("LOWMEM"); keys.use("HIGHMEM");
35 29 : }
36 :
37 23 : AdjacencyMatrixBase::AdjacencyMatrixBase(const ActionOptions& ao):
38 : Action(ao),
39 : MultiColvarBase(ao),
40 23 : connect_id(0),
41 23 : no_third_dim_accum(true),
42 23 : mat(NULL)
43 : {
44 46 : log<<" Bibliography "<<plumed.cite("Tribello, Giberti, Sosso, Salvalaglio and Parrinello, J. Chem. Theory Comput. 13, 1317 (2017)")<<"\n";
45 23 : }
46 :
47 46 : void AdjacencyMatrixBase::parseConnectionDescriptions( const std::string& key, const bool& multiple, const unsigned& nrow_t ) {
48 46 : if( getNumberOfNodeTypes()==1 || (getNumberOfNodeTypes()==2 && nrow_t==1) ) {
49 : std::vector<std::string> sw;
50 45 : if( !multiple ) {
51 44 : sw.resize(1); parse(key,sw[0]);
52 44 : if(sw[0].length()==0) error("could not find " + key + " keyword");
53 : } else {
54 : std::string input;
55 2 : for(int i=1;; i++) {
56 3 : if( !parseNumbered(key, i, input ) ) break;
57 2 : sw.push_back( input );
58 : }
59 : }
60 45 : setupConnector( connect_id, 0, 0, sw );
61 45 : } else {
62 1 : if( multiple ) error("keyword " + key + " does not work with multiple input strings");
63 : unsigned nr, nc;
64 1 : if( nrow_t==0 ) {
65 1 : nr=nc=getNumberOfNodeTypes();
66 : } else {
67 0 : nr=nrow_t; nc = getNumberOfNodeTypes() - nr;
68 : }
69 3 : for(unsigned i=0; i<nr; ++i) {
70 : // Retrieve the base number
71 : unsigned ibase;
72 2 : if( nc<10 ) {
73 2 : ibase=(i+1)*10;
74 0 : } else if ( nc<100 ) {
75 0 : ibase=(i+1)*100;
76 : } else {
77 0 : error("wow this is an error I never would have expected");
78 : }
79 :
80 5 : for(unsigned j=i; j<nc; ++j) {
81 3 : std::vector<std::string> sw(1); parseNumbered(key,ibase+j+1,sw[0]);
82 3 : if(sw[0].length()==0) {
83 0 : std::string num; Tools::convert(ibase+j+1,num);
84 0 : error("could not find " + key + num + " keyword. Need one " + key + " keyword for each distinct base-multicolvar-pair type");
85 : }
86 3 : setupConnector( connect_id, i, j, sw );
87 3 : }
88 : }
89 : }
90 46 : connect_id++;
91 46 : }
92 :
93 5807 : unsigned AdjacencyMatrixBase::getSizeOfInputVectors() const {
94 5807 : if( mybasemulticolvars.size()==0 ) return 2;
95 :
96 5807 : unsigned nq = mybasemulticolvars[0]->getNumberOfQuantities();
97 5807 : for(unsigned i=1; i<mybasemulticolvars.size(); ++i) {
98 0 : if( mybasemulticolvars[i]->getNumberOfQuantities()!=nq ) error("mismatch between vectors in base colvars");
99 : }
100 : return nq;
101 : }
102 :
103 132 : unsigned AdjacencyMatrixBase::getNumberOfNodeTypes() const {
104 132 : unsigned size=mybasemulticolvars.size();
105 : if( size==0 ) return 1;
106 : return size;
107 : }
108 :
109 24 : void AdjacencyMatrixBase::retrieveTypeDimensions( unsigned& nrows, unsigned& ncols, unsigned& ntype ) const {
110 24 : bool allsame=(ablocks[0].size()==ablocks[1].size());
111 24 : if( allsame ) {
112 6952 : for(unsigned i=0; i<ablocks[0].size(); ++i) {
113 6929 : if( ablocks[0][i]!=ablocks[1][i] ) allsame=false;
114 : }
115 : }
116 :
117 24 : if( allsame ) {
118 23 : std::vector<unsigned> types(1); types[0]=atom_lab[ablocks[0][0]].first;
119 6929 : for(unsigned i=1; i<ablocks[0].size(); ++i) {
120 : bool found = false;
121 6913 : for(unsigned j=0; j<types.size(); ++j) {
122 6912 : if( atom_lab[ablocks[0][i]].first==types[j] ) { found=true; break; }
123 : }
124 6906 : if( !found ) types.push_back( atom_lab[ablocks[0][i]].first );
125 : }
126 23 : ntype=0; nrows=ncols=types.size();
127 : } else {
128 1 : std::vector<unsigned> types(1); types[0]=atom_lab[ablocks[0][0]].first;
129 5 : for(unsigned i=1; i<ablocks[0].size(); ++i) {
130 : bool found = false;
131 4 : for(unsigned j=0; j<types.size(); ++j) {
132 4 : if( atom_lab[ablocks[0][i]].first==types[j] ) { found=true; break; }
133 : }
134 4 : if( !found ) types.push_back( atom_lab[ablocks[0][i]].first );
135 : }
136 1 : nrows=ntype=types.size();
137 11 : for(unsigned i=0; i<ablocks[1].size(); ++i) {
138 : bool found = false;
139 10 : for(unsigned j=0; j<types.size(); ++j) {
140 10 : if( atom_lab[ablocks[1][i]].first==types[j] ) { found=true; break; }
141 : }
142 10 : if( !found ) types.push_back( atom_lab[ablocks[1][i]].first );
143 : }
144 1 : if( types.size()==nrows ) { ntype=0; ncols=1; plumed_assert( types.size()==1 && atom_lab[ablocks[0][0]].first==0 ); }
145 0 : else ncols = types.size() - ntype;
146 : }
147 24 : }
148 :
149 23 : void AdjacencyMatrixBase::finishMatrixSetup( const bool& symmetric, const std::vector<AtomNumber>& all_atoms ) {
150 : std::string param;
151 23 : if( symmetric && ablocks[0].size()==ablocks[1].size() ) param="SYMMETRIC";
152 23 : if( !symmetric ) {
153 4 : bool usehbonds=( ablocks[0].size()==ablocks[1].size() );
154 4 : if( usehbonds ) {
155 138 : for(unsigned i=0; i<ablocks[0].size(); ++i) {
156 134 : if( ablocks[0][i]!=ablocks[1][i] ) { usehbonds = false; break; }
157 : }
158 4 : if( usehbonds ) param="HBONDS";
159 : }
160 : }
161 :
162 46 : vesselbase::VesselOptions da("","",0,param,this);
163 23 : Keywords keys; AdjacencyMatrixVessel::registerKeywords( keys );
164 23 : vesselbase::VesselOptions da2(da,keys);
165 23 : auto ves=Tools::make_unique<AdjacencyMatrixVessel>(da2);
166 23 : addVessel( std::move( ves ) );
167 23 : setupMultiColvarBase( all_atoms );
168 46 : }
169 :
170 13 : void AdjacencyMatrixBase::readMaxTwoSpeciesMatrix( const std::string& key0, const std::string& key1, const std::string& key2, const bool& symmetric ) {
171 13 : std::vector<AtomNumber> all_atoms; readTwoGroups( key0, key1, key2, all_atoms );
172 13 : finishMatrixSetup( symmetric, all_atoms );
173 13 : }
174 :
175 10 : void AdjacencyMatrixBase::readMaxThreeSpeciesMatrix( const std::string& key0, const std::string& key1, const std::string& key2, const std::string& keym, const bool& symmetric ) {
176 10 : std::vector<AtomNumber> all_atoms; readGroupKeywords( key0, key1, key2, keym, true, symmetric, all_atoms );
177 10 : finishMatrixSetup( symmetric, all_atoms );
178 10 : }
179 :
180 : // Maybe put this back GAT to check that it is returning an atom number that is one of the nodes
181 : // and not a hydrogen if we are doing HBPAMM
182 : // AtomNumber AdjacencyMatrixBase::getAbsoluteIndexOfCentralAtom(const unsigned& i) const {
183 : // plumed_dbg_assert( i<myinputdata.getFullNumberOfBaseTasks() );
184 : // return myinputdata.getAtomicIndex( i );
185 : // }
186 :
187 0 : void AdjacencyMatrixBase::recalculateMatrixElement( const unsigned& myelem, MultiValue& myvals ) {
188 0 : std::vector<unsigned> myatoms; decodeIndexToAtoms( getTaskCode( myelem ), myatoms );
189 0 : unsigned i=myatoms[0], j=myatoms[1];
190 0 : for(unsigned k=bookeeping(i,j).first; k<bookeeping(i,j).second; ++k) {
191 0 : if( !taskIsCurrentlyActive(k) ) continue;
192 0 : performTask( k, getTaskCode(k), myvals ); // This may not accumulate as we would like GAT
193 : }
194 0 : }
195 :
196 : }
197 : }
|
