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