Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2013-2019 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 21 : void AdjacencyMatrixVessel::registerKeywords( Keywords& keys ) {
30 21 : StoreDataVessel::registerKeywords(keys);
31 63 : keys.addFlag("SYMMETRIC",false,"is the matrix symmetric");
32 63 : keys.addFlag("HBONDS",false,"can we think of the matrix as a undirected graph");
33 21 : }
34 :
35 21 : AdjacencyMatrixVessel::AdjacencyMatrixVessel( const vesselbase::VesselOptions& da ):
36 21 : StoreDataVessel(da)
37 : {
38 21 : function=dynamic_cast<AdjacencyMatrixBase*>( getAction() );
39 21 : plumed_assert( function );
40 63 : parseFlag("SYMMETRIC",symmetric); parseFlag("HBONDS",hbonds);
41 21 : if( symmetric && hbonds ) error("matrix should be either symmetric or hbonds");
42 39 : if( symmetric && function->ablocks[0].size()!=function->ablocks[1].size() ) error("matrix is supposed to be symmetric but nrows!=ncols");
43 23 : if( hbonds && function->ablocks[0].size()!=function->ablocks[1].size() ) error("matrix is supposed to be hbonds but nrows!=ncols");
44 21 : }
45 :
46 775 : bool AdjacencyMatrixVessel::isSymmetric() const {
47 775 : return symmetric;
48 : }
49 :
50 64056 : bool AdjacencyMatrixVessel::undirectedGraph() const {
51 64056 : return ( symmetric || hbonds );
52 : }
53 :
54 1531 : unsigned AdjacencyMatrixVessel::getNumberOfRows() const {
55 3062 : return function->ablocks[0].size();
56 : }
57 :
58 34992 : unsigned AdjacencyMatrixVessel::getNumberOfColumns() const {
59 69984 : return function->ablocks[1].size();
60 : }
61 :
62 39828 : bool AdjacencyMatrixVessel::matrixElementIsActive( const unsigned& ielem, const unsigned& jelem ) const {
63 39828 : return StoreDataVessel::storedValueIsActive( getStoreIndexFromMatrixIndices( ielem, jelem ) );
64 : }
65 :
66 58564 : unsigned AdjacencyMatrixVessel::getStoreIndexFromMatrixIndices( const unsigned& ielem, const unsigned& jelem ) const {
67 60564 : if( !symmetric && !hbonds ) return (function->ablocks[1].size())*ielem + jelem;
68 56564 : if( !symmetric ) {
69 : plumed_dbg_assert( ielem!=jelem );
70 24192 : if( jelem<ielem ) return (function->ablocks[1].size()-1)*ielem + jelem;
71 16128 : return (function->ablocks[1].size()-1)*ielem + jelem - 1;
72 : }
73 40436 : if( ielem>jelem ) return 0.5*ielem*(ielem-1)+jelem;
74 17363 : return 0.5*jelem*(jelem-1) + ielem;
75 : }
76 :
77 15 : AdjacencyMatrixBase* AdjacencyMatrixVessel::getMatrixAction() {
78 15 : return function;
79 : }
80 :
81 24491 : void AdjacencyMatrixVessel::getMatrixIndices( const unsigned& code, unsigned& i, unsigned& j ) const {
82 48982 : std::vector<unsigned> myatoms; function->decodeIndexToAtoms( function->getTaskCode(code), myatoms );
83 24491 : i=myatoms[0]; j=myatoms[1];
84 24491 : if( !undirectedGraph() ) j -= function->ablocks[0].size(); // Have to remove number of columns as returns number in ablocks[1]
85 24491 : }
86 :
87 17 : void AdjacencyMatrixVessel::retrieveMatrix( DynamicList<unsigned>& myactive_elements, Matrix<double>& mymatrix ) {
88 34 : myactive_elements.deactivateAll(); std::vector<double> vals( getNumberOfComponents() );
89 1564 : for(unsigned i=0; i<getNumberOfStoredValues(); ++i) {
90 1547 : retrieveSequentialValue( i, false, vals );
91 1547 : if( vals[0]<epsilon ) continue ;
92 :
93 1547 : myactive_elements.activate(i);
94 3094 : unsigned j, k; getMatrixIndices( function->getPositionInFullTaskList(i), k, j );
95 :
96 6188 : if( symmetric ) mymatrix(k,j)=mymatrix(j,k)=vals[0]*vals[1];
97 0 : else mymatrix(k,j)=vals[0]*vals[1];
98 : }
99 17 : myactive_elements.updateActiveMembers();
100 17 : }
101 :
102 24 : void AdjacencyMatrixVessel::retrieveAdjacencyLists( std::vector<unsigned>& nneigh, Matrix<unsigned>& adj_list ) {
103 : plumed_dbg_assert( undirectedGraph() );
104 : // Currently everything has zero neighbors
105 24627 : for(unsigned i=0; i<nneigh.size(); ++i) nneigh[i]=0;
106 :
107 : // And set up the adjacency list
108 24 : std::vector<double> myvals( getNumberOfComponents() );
109 214258 : for(unsigned i=0; i<getNumberOfStoredValues(); ++i) {
110 : // Check if atoms are connected
111 214234 : retrieveSequentialValue( i, false, myvals );
112 214234 : if( myvals[0]<epsilon || myvals[1]<epsilon ) continue ;
113 :
114 30026 : unsigned j, k; getMatrixIndices( function->getPositionInFullTaskList(i), k, j );
115 :
116 60052 : if( nneigh[j]>=adj_list.ncols() || nneigh[k]>=adj_list.ncols() ) error("adjacency lists are not large enough, increase maxconnections");
117 : // Store if atoms are connected
118 : // unsigned j, k; getMatrixIndices( i, k, j );
119 45039 : adj_list(k,nneigh[k])=j; nneigh[k]++;
120 45039 : adj_list(j,nneigh[j])=k; nneigh[j]++;
121 : }
122 24 : }
123 :
124 4 : void AdjacencyMatrixVessel::retrieveEdgeList( unsigned& nedge, std::vector<std::pair<unsigned,unsigned> >& edge_list ) {
125 4 : plumed_dbg_assert( undirectedGraph() ); nedge=0;
126 4 : std::vector<double> myvals( getNumberOfComponents() );
127 8 : if( getNumberOfStoredValues()>edge_list.size() ) error("adjacency lists are not large enough, increase maxconnections");
128 :
129 124576 : for(unsigned i=0; i<getNumberOfStoredValues(); ++i) {
130 : // Check if atoms are connected
131 124572 : retrieveSequentialValue( i, false, myvals );
132 124572 : if( myvals[0]<epsilon || myvals[1]<epsilon ) continue ;
133 :
134 19152 : getMatrixIndices( function->getPositionInFullTaskList(i), edge_list[nedge].first, edge_list[nedge].second );
135 6384 : nedge++;
136 : }
137 4 : }
138 :
139 0 : bool AdjacencyMatrixVessel::nodesAreConnected( const unsigned& iatom, const unsigned& jatom ) const {
140 0 : if( !matrixElementIsActive( iatom, jatom ) ) return false;
141 0 : unsigned ind=getStoreIndexFromMatrixIndices( iatom, jatom );
142 :
143 0 : std::vector<double> myvals( getNumberOfComponents() );
144 0 : retrieveValueWithIndex( ind, false, myvals );
145 0 : return ( myvals[0]>epsilon && myvals[1]>epsilon );
146 : }
147 :
148 8 : double AdjacencyMatrixVessel::getCutoffForConnection() const {
149 8 : return function->getLinkCellCutoff();
150 : }
151 :
152 : }
153 4839 : }
154 :
|