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 47 : void AdjacencyMatrixBase::registerKeywords( Keywords& keys ) {
33 47 : multicolvar::MultiColvarBase::registerKeywords( keys );
34 47 : keys.remove("LOWMEM");
35 47 : keys.use("HIGHMEM");
36 47 : }
37 :
38 23 : AdjacencyMatrixBase::AdjacencyMatrixBase(const ActionOptions& ao):
39 : Action(ao),
40 : MultiColvarBase(ao),
41 23 : connect_id(0),
42 23 : no_third_dim_accum(true),
43 23 : mat(NULL) {
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);
52 44 : parse(key,sw[0]);
53 44 : if(sw[0].length()==0) {
54 0 : error("could not find " + key + " keyword");
55 : }
56 : } else {
57 : std::string input;
58 2 : for(int i=1;; i++) {
59 3 : if( !parseNumbered(key, i, input ) ) {
60 : break;
61 : }
62 2 : sw.push_back( input );
63 : }
64 : }
65 45 : setupConnector( connect_id, 0, 0, sw );
66 45 : } else {
67 1 : if( multiple ) {
68 0 : error("keyword " + key + " does not work with multiple input strings");
69 : }
70 : unsigned nr, nc;
71 1 : if( nrow_t==0 ) {
72 1 : nr=nc=getNumberOfNodeTypes();
73 : } else {
74 : nr=nrow_t;
75 0 : nc = getNumberOfNodeTypes() - nr;
76 : }
77 3 : for(unsigned i=0; i<nr; ++i) {
78 : // Retrieve the base number
79 : unsigned ibase;
80 2 : if( nc<10 ) {
81 2 : ibase=(i+1)*10;
82 0 : } else if ( nc<100 ) {
83 0 : ibase=(i+1)*100;
84 : } else {
85 0 : error("wow this is an error I never would have expected");
86 : }
87 :
88 5 : for(unsigned j=i; j<nc; ++j) {
89 3 : std::vector<std::string> sw(1);
90 3 : parseNumbered(key,ibase+j+1,sw[0]);
91 3 : if(sw[0].length()==0) {
92 : std::string num;
93 0 : Tools::convert(ibase+j+1,num);
94 0 : error("could not find " + key + num + " keyword. Need one " + key + " keyword for each distinct base-multicolvar-pair type");
95 : }
96 3 : setupConnector( connect_id, i, j, sw );
97 3 : }
98 : }
99 : }
100 46 : connect_id++;
101 46 : }
102 :
103 5807 : unsigned AdjacencyMatrixBase::getSizeOfInputVectors() const {
104 5807 : if( mybasemulticolvars.size()==0 ) {
105 : return 2;
106 : }
107 :
108 5807 : unsigned nq = mybasemulticolvars[0]->getNumberOfQuantities();
109 5807 : for(unsigned i=1; i<mybasemulticolvars.size(); ++i) {
110 0 : if( mybasemulticolvars[i]->getNumberOfQuantities()!=nq ) {
111 0 : error("mismatch between vectors in base colvars");
112 : }
113 : }
114 : return nq;
115 : }
116 :
117 132 : unsigned AdjacencyMatrixBase::getNumberOfNodeTypes() const {
118 132 : unsigned size=mybasemulticolvars.size();
119 : if( size==0 ) {
120 : return 1;
121 : }
122 : return size;
123 : }
124 :
125 24 : void AdjacencyMatrixBase::retrieveTypeDimensions( unsigned& nrows, unsigned& ncols, unsigned& ntype ) const {
126 24 : bool allsame=(ablocks[0].size()==ablocks[1].size());
127 24 : if( allsame ) {
128 6952 : for(unsigned i=0; i<ablocks[0].size(); ++i) {
129 6929 : if( ablocks[0][i]!=ablocks[1][i] ) {
130 : allsame=false;
131 : }
132 : }
133 : }
134 :
135 24 : if( allsame ) {
136 23 : std::vector<unsigned> types(1);
137 23 : types[0]=atom_lab[ablocks[0][0]].first;
138 6929 : for(unsigned i=1; i<ablocks[0].size(); ++i) {
139 : bool found = false;
140 6913 : for(unsigned j=0; j<types.size(); ++j) {
141 6912 : if( atom_lab[ablocks[0][i]].first==types[j] ) {
142 : found=true;
143 : break;
144 : }
145 : }
146 6906 : if( !found ) {
147 1 : types.push_back( atom_lab[ablocks[0][i]].first );
148 : }
149 : }
150 23 : ntype=0;
151 23 : nrows=ncols=types.size();
152 : } else {
153 1 : std::vector<unsigned> types(1);
154 1 : types[0]=atom_lab[ablocks[0][0]].first;
155 5 : for(unsigned i=1; i<ablocks[0].size(); ++i) {
156 : bool found = false;
157 4 : for(unsigned j=0; j<types.size(); ++j) {
158 4 : if( atom_lab[ablocks[0][i]].first==types[j] ) {
159 : found=true;
160 : break;
161 : }
162 : }
163 4 : if( !found ) {
164 0 : types.push_back( atom_lab[ablocks[0][i]].first );
165 : }
166 : }
167 1 : nrows=ntype=types.size();
168 11 : for(unsigned i=0; i<ablocks[1].size(); ++i) {
169 : bool found = false;
170 10 : for(unsigned j=0; j<types.size(); ++j) {
171 10 : if( atom_lab[ablocks[1][i]].first==types[j] ) {
172 : found=true;
173 : break;
174 : }
175 : }
176 10 : if( !found ) {
177 0 : types.push_back( atom_lab[ablocks[1][i]].first );
178 : }
179 : }
180 1 : if( types.size()==nrows ) {
181 1 : ntype=0;
182 1 : ncols=1;
183 1 : plumed_assert( types.size()==1 && atom_lab[ablocks[0][0]].first==0 );
184 : } else {
185 0 : ncols = types.size() - ntype;
186 : }
187 : }
188 24 : }
189 :
190 23 : void AdjacencyMatrixBase::finishMatrixSetup( const bool& symmetric, const std::vector<AtomNumber>& all_atoms ) {
191 : std::string param;
192 23 : if( symmetric && ablocks[0].size()==ablocks[1].size() ) {
193 : param="SYMMETRIC";
194 : }
195 23 : if( !symmetric ) {
196 4 : bool usehbonds=( ablocks[0].size()==ablocks[1].size() );
197 4 : if( usehbonds ) {
198 138 : for(unsigned i=0; i<ablocks[0].size(); ++i) {
199 134 : if( ablocks[0][i]!=ablocks[1][i] ) {
200 : usehbonds = false;
201 : break;
202 : }
203 : }
204 4 : if( usehbonds ) {
205 : param="HBONDS";
206 : }
207 : }
208 : }
209 :
210 46 : vesselbase::VesselOptions da("","",0,param,this);
211 23 : Keywords keys;
212 23 : AdjacencyMatrixVessel::registerKeywords( keys );
213 23 : vesselbase::VesselOptions da2(da,keys);
214 23 : auto ves=Tools::make_unique<AdjacencyMatrixVessel>(da2);
215 23 : addVessel( std::move( ves ) );
216 23 : setupMultiColvarBase( all_atoms );
217 46 : }
218 :
219 13 : void AdjacencyMatrixBase::readMaxTwoSpeciesMatrix( const std::string& key0, const std::string& key1, const std::string& key2, const bool& symmetric ) {
220 : std::vector<AtomNumber> all_atoms;
221 13 : readTwoGroups( key0, key1, key2, all_atoms );
222 13 : finishMatrixSetup( symmetric, all_atoms );
223 13 : }
224 :
225 10 : void AdjacencyMatrixBase::readMaxThreeSpeciesMatrix( const std::string& key0, const std::string& key1, const std::string& key2, const std::string& keym, const bool& symmetric ) {
226 : std::vector<AtomNumber> all_atoms;
227 10 : readGroupKeywords( key0, key1, key2, keym, true, symmetric, all_atoms );
228 10 : finishMatrixSetup( symmetric, all_atoms );
229 10 : }
230 :
231 : // Maybe put this back GAT to check that it is returning an atom number that is one of the nodes
232 : // and not a hydrogen if we are doing HBPAMM
233 : // AtomNumber AdjacencyMatrixBase::getAbsoluteIndexOfCentralAtom(const unsigned& i) const {
234 : // plumed_dbg_assert( i<myinputdata.getFullNumberOfBaseTasks() );
235 : // return myinputdata.getAtomicIndex( i );
236 : // }
237 :
238 0 : void AdjacencyMatrixBase::recalculateMatrixElement( const unsigned& myelem, MultiValue& myvals ) {
239 : std::vector<unsigned> myatoms;
240 0 : decodeIndexToAtoms( getTaskCode( myelem ), myatoms );
241 0 : unsigned i=myatoms[0], j=myatoms[1];
242 0 : for(unsigned k=bookeeping(i,j).first; k<bookeeping(i,j).second; ++k) {
243 0 : if( !taskIsCurrentlyActive(k) ) {
244 0 : continue;
245 : }
246 0 : performTask( k, getTaskCode(k), myvals ); // This may not accumulate as we would like GAT
247 : }
248 0 : }
249 :
250 : }
251 : }
|
