Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2015-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 "AnalysisBase.h"
      23             : #include "tools/PDB.h"
      24             : #include "core/ActionRegister.h"
      25             : #include "reference/MetricRegister.h"
      26             : #include "reference/ReferenceConfiguration.h"
      27             : 
      28             : //+PLUMEDOC ANALYSIS EUCLIDEAN_DISSIMILARITIES
      29             : /*
      30             : Calculate the matrix of dissimilarities between a trajectory of atomic configurations.
      31             : 
      32             : \par Examples
      33             : 
      34             : */
      35             : //+ENDPLUMEDOC
      36             : 
      37             : namespace PLMD {
      38             : namespace analysis {
      39             : 
      40             : class EuclideanDissimilarityMatrix : public AnalysisBase {
      41             : private:
      42             :   PDB mypdb;
      43             :   std::string mtype;
      44             :   Matrix<double> dissimilarities;
      45             : public:
      46             :   static void registerKeywords( Keywords& keys );
      47             :   explicit EuclideanDissimilarityMatrix( const ActionOptions& ao );
      48             : /// Do the analysis
      49             :   void performAnalysis() override;
      50             : /// This ensures that classes that use this data know that dissimilarities were set
      51         436 :   bool dissimilaritiesWereSet() const override {
      52         436 :     return true;
      53             :   }
      54             : /// Get information on how to calculate dissimilarities
      55             :   std::string getDissimilarityInstruction() const override;
      56             : /// Get the squared dissimilarity between two reference configurations
      57             :   double getDissimilarity( const unsigned& i, const unsigned& j ) override;
      58             : /// This is just to deal with ActionWithVessel
      59           0 :   void performTask( const unsigned&, const unsigned&, MultiValue& ) const override {
      60           0 :     plumed_error();
      61             :   }
      62             : };
      63             : 
      64       13811 : PLUMED_REGISTER_ACTION(EuclideanDissimilarityMatrix,"EUCLIDEAN_DISSIMILARITIES")
      65             : 
      66          17 : void EuclideanDissimilarityMatrix::registerKeywords( Keywords& keys ) {
      67          17 :   AnalysisBase::registerKeywords( keys );
      68          17 :   keys.use("ARG");
      69          34 :   keys.reset_style("ARG","optional");
      70          34 :   keys.add("compulsory","METRIC","EUCLIDEAN","the method that you are going to use to measure the distances between points");
      71          34 :   keys.add("atoms","ATOMS","the list of atoms that you are going to use in the measure of distance that you are using");
      72          17 : }
      73             : 
      74          13 : EuclideanDissimilarityMatrix::EuclideanDissimilarityMatrix( const ActionOptions& ao ):
      75             :   Action(ao),
      76          13 :   AnalysisBase(ao) {
      77          26 :   parse("METRIC",mtype);
      78             :   std::vector<AtomNumber> atoms;
      79          13 :   if( my_input_data->getNumberOfAtoms()>0 ) {
      80           6 :     parseAtomList("ATOMS",atoms);
      81           3 :     if( atoms.size()!=0 ) {
      82           0 :       mypdb.setAtomNumbers( atoms );
      83           0 :       for(unsigned i=0; i<atoms.size(); ++i) {
      84             :         bool found=false;
      85           0 :         for(unsigned j=0; j<my_input_data->getAtomIndexes().size(); ++j) {
      86           0 :           if( my_input_data->getAtomIndexes()[j]==atoms[i] ) {
      87             :             found=true;
      88             :             break;
      89             :           }
      90             :         }
      91           0 :         if( !found ) {
      92             :           std::string num;
      93           0 :           Tools::convert( atoms[i].serial(), num );
      94           0 :           error("atom number " + num + " is not stored in any action that has been input");
      95             :         }
      96             :       }
      97           0 :       mypdb.addBlockEnd( atoms.size() );
      98           3 :     } else if( getNumberOfArguments()==0 ) {
      99           1 :       mypdb.setAtomNumbers( my_input_data->getAtomIndexes() );
     100           1 :       mypdb.addBlockEnd( my_input_data->getAtomIndexes().size() );
     101           1 :       if( mtype=="EUCLIDEAN" ) {
     102             :         mtype="OPTIMAL";
     103             :       }
     104             :     }
     105             :   }
     106          13 :   log.printf("  measuring distances using %s metric \n",mtype.c_str() );
     107          13 :   if( my_input_data->getArgumentNames().size()>0 ) {
     108          12 :     if( getNumberOfArguments()==0 && atoms.size()==0 ) {
     109          10 :       std::vector<std::string> argnames( my_input_data->getArgumentNames() );
     110          10 :       mypdb.setArgumentNames( argnames );
     111          10 :       requestArguments( my_input_data->getArgumentList() );
     112          10 :     } else {
     113           2 :       std::vector<Value*> myargs( getArguments() );
     114           2 :       std::vector<std::string> inargnames( my_input_data->getArgumentNames() );
     115           2 :       std::vector<std::string> argnames( myargs.size() );
     116           6 :       for(unsigned i=0; i<myargs.size(); ++i) {
     117           4 :         argnames[i]=myargs[i]->getName();
     118             :         bool found=false;
     119           6 :         for(unsigned j=0; j<inargnames.size(); ++j) {
     120           6 :           if( argnames[i]==inargnames[j] ) {
     121             :             found=true;
     122             :             break;
     123             :           }
     124             :         }
     125           4 :         if( !found ) {
     126           0 :           error("input named " + my_input_data->getLabel() + " does not store/calculate quantity named " + argnames[i] );
     127             :         }
     128             :       }
     129           2 :       mypdb.setArgumentNames( argnames );
     130           2 :       requestArguments( myargs );
     131           2 :     }
     132             :   }
     133          13 : }
     134             : 
     135          20 : void EuclideanDissimilarityMatrix::performAnalysis() {
     136             :   // Resize dissimilarities matrix and set all elements to zero
     137          20 :   if( !usingLowMem() ) {
     138          20 :     dissimilarities.resize( getNumberOfDataPoints(), getNumberOfDataPoints() );
     139             :     dissimilarities=0;
     140             :   }
     141          20 : }
     142             : 
     143         413 : std::string EuclideanDissimilarityMatrix::getDissimilarityInstruction() const {
     144         413 :   return "TYPE=" + mtype;
     145             : }
     146             : 
     147      611659 : double EuclideanDissimilarityMatrix::getDissimilarity( const unsigned& iframe, const unsigned& jframe ) {
     148             :   plumed_dbg_assert( iframe<getNumberOfDataPoints() && jframe<getNumberOfDataPoints() );
     149      611659 :   if( !usingLowMem() ) {
     150      611659 :     if( dissimilarities(iframe,jframe)>0. ) {
     151             :       return dissimilarities(iframe,jframe);
     152             :     }
     153             :   }
     154      256165 :   if( iframe!=jframe ) {
     155             :     double dd;
     156      255613 :     getStoredData( iframe, true ).transferDataToPDB( mypdb );
     157      255613 :     auto myref1=metricRegister().create<ReferenceConfiguration>(mtype, mypdb);
     158      255613 :     getStoredData( jframe, true ).transferDataToPDB( mypdb );
     159      255613 :     auto myref2=metricRegister().create<ReferenceConfiguration>(mtype, mypdb);
     160      255613 :     if( !usingLowMem() ) {
     161      255613 :       dd=dissimilarities(iframe,jframe) = dissimilarities(jframe,iframe) = distance( getPbc(), getArguments(), myref1.get(), myref2.get(), true );
     162             :     } else {
     163           0 :       dd=distance( getPbc(), getArguments(), myref1.get(), myref2.get(), true );
     164             :     }
     165             :     return dd;
     166      255613 :   }
     167             :   return 0.0;
     168             : }
     169             : 
     170             : }
     171             : }