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             : {
      38          23 :   function=dynamic_cast<AdjacencyMatrixBase*>( getAction() );
      39          23 :   plumed_assert( function );
      40          46 :   parseFlag("SYMMETRIC",symmetric); parseFlag("HBONDS",hbonds);
      41          23 :   if( symmetric && hbonds ) error("matrix should be either symmetric or hbonds");
      42          23 :   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          23 : }
      45             : 
      46         775 : bool AdjacencyMatrixVessel::isSymmetric() const {
      47         775 :   return symmetric;
      48             : }
      49             : 
      50       72251 : bool AdjacencyMatrixVessel::undirectedGraph() const {
      51       72251 :   return ( symmetric || hbonds );
      52             : }
      53             : 
      54        1596 : unsigned AdjacencyMatrixVessel::getNumberOfRows() const {
      55        1596 :   return function->ablocks[0].size();
      56             : }
      57             : 
      58       39153 : unsigned AdjacencyMatrixVessel::getNumberOfColumns() const {
      59       39153 :   return function->ablocks[1].size();
      60             : }
      61             : 
      62       47892 : bool AdjacencyMatrixVessel::matrixElementIsActive( const unsigned& ielem, const unsigned& jelem ) const {
      63       47892 :   return StoreDataVessel::storedValueIsActive( getStoreIndexFromMatrixIndices( ielem, jelem ) );
      64             : }
      65             : 
      66       74692 : unsigned AdjacencyMatrixVessel::getStoreIndexFromMatrixIndices( const unsigned& ielem, const unsigned& jelem ) const {
      67       74692 :   if( !symmetric && !hbonds ) return (function->ablocks[1].size())*ielem + jelem;
      68       72692 :   if( !symmetric ) {
      69             :     plumed_dbg_assert( ielem!=jelem );
      70       32256 :     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          17 : AdjacencyMatrixBase* AdjacencyMatrixVessel::getMatrixAction() {
      78          17 :   return function;
      79             : }
      80             : 
      81       24491 : void AdjacencyMatrixVessel::getMatrixIndices( const unsigned& code, unsigned& i, unsigned& j ) const {
      82       24491 :   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          17 :   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        1547 :     unsigned j, k; getMatrixIndices( function->getPositionInFullTaskList(i), k, j );
      95             : 
      96        1547 :     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        8217 :   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       15013 :     unsigned j, k; getMatrixIndices( function->getPositionInFullTaskList(i), k, j );
     115             : 
     116       15013 :     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       15013 :     adj_list(k,nneigh[k])=j; nneigh[k]++;
     120       15013 :     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           4 :   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        6384 :     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             : }
     154             :