Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-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 "AdjacencyMatrixVessel.h"
23 : #include "AdjacencyMatrixBase.h"
24 : #include "vesselbase/ActionWithVessel.h"
25 :
26 : namespace PLMD {
27 : namespace adjmat {
28 :
29 23 : void AdjacencyMatrixVessel::registerKeywords( Keywords& keys ) {
30 23 : StoreDataVessel::registerKeywords(keys);
31 46 : keys.addFlag("SYMMETRIC",false,"is the matrix symmetric");
32 46 : keys.addFlag("HBONDS",false,"can we think of the matrix as a undirected graph");
33 23 : }
34 :
35 23 : AdjacencyMatrixVessel::AdjacencyMatrixVessel( const vesselbase::VesselOptions& da ):
36 23 : StoreDataVessel(da) {
37 23 : function=dynamic_cast<AdjacencyMatrixBase*>( getAction() );
38 23 : plumed_assert( function );
39 23 : parseFlag("SYMMETRIC",symmetric);
40 23 : parseFlag("HBONDS",hbonds);
41 23 : if( symmetric && hbonds ) {
42 0 : error("matrix should be either symmetric or hbonds");
43 : }
44 23 : if( symmetric && function->ablocks[0].size()!=function->ablocks[1].size() ) {
45 0 : error("matrix is supposed to be symmetric but nrows!=ncols");
46 : }
47 23 : if( hbonds && function->ablocks[0].size()!=function->ablocks[1].size() ) {
48 0 : error("matrix is supposed to be hbonds but nrows!=ncols");
49 : }
50 23 : }
51 :
52 775 : bool AdjacencyMatrixVessel::isSymmetric() const {
53 775 : return symmetric;
54 : }
55 :
56 72251 : bool AdjacencyMatrixVessel::undirectedGraph() const {
57 72251 : return ( symmetric || hbonds );
58 : }
59 :
60 1596 : unsigned AdjacencyMatrixVessel::getNumberOfRows() const {
61 1596 : return function->ablocks[0].size();
62 : }
63 :
64 39153 : unsigned AdjacencyMatrixVessel::getNumberOfColumns() const {
65 39153 : return function->ablocks[1].size();
66 : }
67 :
68 47892 : bool AdjacencyMatrixVessel::matrixElementIsActive( const unsigned& ielem, const unsigned& jelem ) const {
69 47892 : return StoreDataVessel::storedValueIsActive( getStoreIndexFromMatrixIndices( ielem, jelem ) );
70 : }
71 :
72 74692 : unsigned AdjacencyMatrixVessel::getStoreIndexFromMatrixIndices( const unsigned& ielem, const unsigned& jelem ) const {
73 74692 : if( !symmetric && !hbonds ) {
74 2000 : return (function->ablocks[1].size())*ielem + jelem;
75 : }
76 72692 : if( !symmetric ) {
77 : plumed_dbg_assert( ielem!=jelem );
78 32256 : if( jelem<ielem ) {
79 16128 : return (function->ablocks[1].size()-1)*ielem + jelem;
80 : }
81 16128 : return (function->ablocks[1].size()-1)*ielem + jelem - 1;
82 : }
83 40436 : if( ielem>jelem ) {
84 23073 : return 0.5*ielem*(ielem-1)+jelem;
85 : }
86 17363 : return 0.5*jelem*(jelem-1) + ielem;
87 : }
88 :
89 17 : AdjacencyMatrixBase* AdjacencyMatrixVessel::getMatrixAction() {
90 17 : return function;
91 : }
92 :
93 24491 : void AdjacencyMatrixVessel::getMatrixIndices( const unsigned& code, unsigned& i, unsigned& j ) const {
94 : std::vector<unsigned> myatoms;
95 24491 : function->decodeIndexToAtoms( function->getTaskCode(code), myatoms );
96 24491 : i=myatoms[0];
97 24491 : j=myatoms[1];
98 24491 : if( !undirectedGraph() ) {
99 0 : j -= function->ablocks[0].size(); // Have to remove number of columns as returns number in ablocks[1]
100 : }
101 24491 : }
102 :
103 17 : void AdjacencyMatrixVessel::retrieveMatrix( DynamicList<unsigned>& myactive_elements, Matrix<double>& mymatrix ) {
104 17 : myactive_elements.deactivateAll();
105 17 : std::vector<double> vals( getNumberOfComponents() );
106 1564 : for(unsigned i=0; i<getNumberOfStoredValues(); ++i) {
107 1547 : retrieveSequentialValue( i, false, vals );
108 1547 : if( vals[0]<epsilon ) {
109 0 : continue ;
110 : }
111 :
112 1547 : myactive_elements.activate(i);
113 : unsigned j, k;
114 1547 : getMatrixIndices( function->getPositionInFullTaskList(i), k, j );
115 :
116 1547 : if( symmetric ) {
117 1547 : mymatrix(k,j)=mymatrix(j,k)=vals[0]*vals[1];
118 : } else {
119 0 : mymatrix(k,j)=vals[0]*vals[1];
120 : }
121 : }
122 17 : myactive_elements.updateActiveMembers();
123 17 : }
124 :
125 24 : void AdjacencyMatrixVessel::retrieveAdjacencyLists( std::vector<unsigned>& nneigh, Matrix<unsigned>& adj_list ) {
126 : plumed_dbg_assert( undirectedGraph() );
127 : // Currently everything has zero neighbors
128 8217 : for(unsigned i=0; i<nneigh.size(); ++i) {
129 8193 : nneigh[i]=0;
130 : }
131 :
132 : // And set up the adjacency list
133 24 : std::vector<double> myvals( getNumberOfComponents() );
134 214258 : for(unsigned i=0; i<getNumberOfStoredValues(); ++i) {
135 : // Check if atoms are connected
136 214234 : retrieveSequentialValue( i, false, myvals );
137 214234 : if( myvals[0]<epsilon || myvals[1]<epsilon ) {
138 199221 : continue ;
139 : }
140 :
141 : unsigned j, k;
142 15013 : getMatrixIndices( function->getPositionInFullTaskList(i), k, j );
143 :
144 15013 : if( nneigh[j]>=adj_list.ncols() || nneigh[k]>=adj_list.ncols() ) {
145 0 : error("adjacency lists are not large enough, increase maxconnections");
146 : }
147 : // Store if atoms are connected
148 : // unsigned j, k; getMatrixIndices( i, k, j );
149 15013 : adj_list(k,nneigh[k])=j;
150 15013 : nneigh[k]++;
151 15013 : adj_list(j,nneigh[j])=k;
152 15013 : nneigh[j]++;
153 : }
154 24 : }
155 :
156 4 : void AdjacencyMatrixVessel::retrieveEdgeList( unsigned& nedge, std::vector<std::pair<unsigned,unsigned> >& edge_list ) {
157 : plumed_dbg_assert( undirectedGraph() );
158 4 : nedge=0;
159 4 : std::vector<double> myvals( getNumberOfComponents() );
160 4 : if( getNumberOfStoredValues()>edge_list.size() ) {
161 0 : error("adjacency lists are not large enough, increase maxconnections");
162 : }
163 :
164 124576 : for(unsigned i=0; i<getNumberOfStoredValues(); ++i) {
165 : // Check if atoms are connected
166 124572 : retrieveSequentialValue( i, false, myvals );
167 124572 : if( myvals[0]<epsilon || myvals[1]<epsilon ) {
168 118188 : continue ;
169 : }
170 :
171 6384 : getMatrixIndices( function->getPositionInFullTaskList(i), edge_list[nedge].first, edge_list[nedge].second );
172 6384 : nedge++;
173 : }
174 4 : }
175 :
176 0 : bool AdjacencyMatrixVessel::nodesAreConnected( const unsigned& iatom, const unsigned& jatom ) const {
177 0 : if( !matrixElementIsActive( iatom, jatom ) ) {
178 : return false;
179 : }
180 0 : unsigned ind=getStoreIndexFromMatrixIndices( iatom, jatom );
181 :
182 0 : std::vector<double> myvals( getNumberOfComponents() );
183 0 : retrieveValueWithIndex( ind, false, myvals );
184 0 : return ( myvals[0]>epsilon && myvals[1]>epsilon );
185 : }
186 :
187 8 : double AdjacencyMatrixVessel::getCutoffForConnection() const {
188 8 : return function->getLinkCellCutoff();
189 : }
190 :
191 : }
192 : }
193 :
|
