Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2015-2018 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 19 : void AdjacencyMatrixBase::registerKeywords( Keywords& keys ) {
33 19 : multicolvar::MultiColvarBase::registerKeywords( keys );
34 19 : keys.remove("LOWMEM"); keys.use("HIGHMEM");
35 19 : }
36 :
37 15 : AdjacencyMatrixBase::AdjacencyMatrixBase(const ActionOptions& ao):
38 : Action(ao),
39 : MultiColvarBase(ao),
40 : connect_id(0),
41 : no_third_dim_accum(true),
42 15 : mat(NULL)
43 : {
44 15 : log<<" Bibliography "<<plumed.cite("Tribello, Giberti, Sosso, Salvalaglio and Parrinello, J. Chem. Theory Comput. 13, 1317 (2017)")<<"\n";
45 15 : }
46 :
47 20 : void AdjacencyMatrixBase::parseConnectionDescriptions( const std::string& key, const bool& multiple, const unsigned& nrow_t ) {
48 20 : if( getNumberOfNodeTypes()==1 || (getNumberOfNodeTypes()==2 && nrow_t==1) ) {
49 19 : std::vector<std::string> sw;
50 19 : if( !multiple ) {
51 18 : sw.resize(1); parse(key,sw[0]);
52 18 : if(sw[0].length()==0) error("could not find " + key + " keyword");
53 : } else {
54 1 : std::string input;
55 3 : for(int i=1;; i++) {
56 3 : if( !parseNumbered(key, i, input ) ) break;
57 2 : sw.push_back( input );
58 2 : }
59 : }
60 19 : setupConnector( connect_id, 0, 0, sw );
61 : } 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 20 : connect_id++;
91 20 : }
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 5807 : return nq;
101 : }
102 :
103 34 : unsigned AdjacencyMatrixBase::getNumberOfNodeTypes() const {
104 34 : unsigned size=mybasemulticolvars.size();
105 34 : if( size==0 ) return 1;
106 18 : return size;
107 : }
108 :
109 16 : void AdjacencyMatrixBase::retrieveTypeDimensions( unsigned& nrows, unsigned& ncols, unsigned& ntype ) const {
110 16 : bool allsame=(ablocks[0].size()==ablocks[1].size());
111 16 : if( allsame ) {
112 5540 : for(unsigned i=0; i<ablocks[0].size(); ++i) {
113 5525 : if( ablocks[0][i]!=ablocks[1][i] ) allsame=false;
114 : }
115 : }
116 :
117 16 : if( allsame ) {
118 15 : std::vector<unsigned> types(1); types[0]=atom_lab[ablocks[0][0]].first;
119 5525 : for(unsigned i=1; i<ablocks[0].size(); ++i) {
120 5510 : bool found = false;
121 5517 : for(unsigned j=0; j<types.size(); ++j) {
122 5516 : if( atom_lab[ablocks[0][i]].first==types[j] ) { found=true; break; }
123 : }
124 5510 : if( !found ) types.push_back( atom_lab[ablocks[0][i]].first );
125 : }
126 15 : 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 4 : 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 10 : 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 16 : }
148 :
149 15 : void AdjacencyMatrixBase::finishMatrixSetup( const bool& symmetric, const std::vector<AtomNumber>& all_atoms ) {
150 15 : std::string param;
151 15 : if( symmetric && ablocks[0].size()==ablocks[1].size() ) param="SYMMETRIC";
152 15 : if( !symmetric ) {
153 2 : bool usehbonds=( ablocks[0].size()==ablocks[1].size() );
154 2 : if( usehbonds ) {
155 69 : for(unsigned i=0; i<ablocks[0].size(); ++i) {
156 67 : if( ablocks[0][i]!=ablocks[1][i] ) { usehbonds = false; break; }
157 : }
158 2 : if( usehbonds ) param="HBONDS";
159 : }
160 : }
161 :
162 30 : vesselbase::VesselOptions da("","",0,param,this);
163 30 : Keywords keys; AdjacencyMatrixVessel::registerKeywords( keys );
164 30 : vesselbase::VesselOptions da2(da,keys); mat = new AdjacencyMatrixVessel(da2); addVessel( mat );
165 30 : setupMultiColvarBase( all_atoms );
166 15 : }
167 :
168 13 : void AdjacencyMatrixBase::readMaxTwoSpeciesMatrix( const std::string& key0, const std::string& key1, const std::string& key2, const bool& symmetric ) {
169 13 : std::vector<AtomNumber> all_atoms; readTwoGroups( key0, key1, key2, all_atoms );
170 13 : finishMatrixSetup( symmetric, all_atoms );
171 13 : }
172 :
173 2 : void AdjacencyMatrixBase::readMaxThreeSpeciesMatrix( const std::string& key0, const std::string& key1, const std::string& key2, const std::string& keym, const bool& symmetric ) {
174 2 : std::vector<AtomNumber> all_atoms; readGroupKeywords( key0, key1, key2, keym, true, symmetric, all_atoms );
175 2 : finishMatrixSetup( symmetric, all_atoms );
176 2 : }
177 :
178 : // Maybe put this back GAT to check that it is returning an atom number that is one of the nodes
179 : // and not a hydrogen if we are doing HBPAMM
180 : // AtomNumber AdjacencyMatrixBase::getAbsoluteIndexOfCentralAtom(const unsigned& i) const {
181 : // plumed_dbg_assert( i<myinputdata.getFullNumberOfBaseTasks() );
182 : // return myinputdata.getAtomicIndex( i );
183 : // }
184 :
185 21658 : void AdjacencyMatrixBase::recalculateMatrixElement( const unsigned& myelem, MultiValue& myvals ) {
186 21658 : std::vector<unsigned> myatoms; decodeIndexToAtoms( getTaskCode( myelem ), myatoms );
187 21658 : unsigned i=myatoms[0], j=myatoms[1];
188 43316 : for(unsigned k=bookeeping(i,j).first; k<bookeeping(i,j).second; ++k) {
189 21658 : if( !taskIsCurrentlyActive(k) ) continue;
190 21658 : performTask( k, getTaskCode(k), myvals ); // This may not accumulate as we would like GAT
191 21658 : }
192 21658 : }
193 :
194 : }
195 2523 : }
|
