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 { return true; }
      52             : /// Get information on how to calculate dissimilarities
      53             :   std::string getDissimilarityInstruction() const override;
      54             : /// Get the squared dissimilarity between two reference configurations
      55             :   double getDissimilarity( const unsigned& i, const unsigned& j ) override;
      56             : /// This is just to deal with ActionWithVessel
      57           0 :   void performTask( const unsigned&, const unsigned&, MultiValue& ) const override { plumed_error(); }
      58             : };
      59             : 
      60       10445 : PLUMED_REGISTER_ACTION(EuclideanDissimilarityMatrix,"EUCLIDEAN_DISSIMILARITIES")
      61             : 
      62          14 : void EuclideanDissimilarityMatrix::registerKeywords( Keywords& keys ) {
      63          42 :   AnalysisBase::registerKeywords( keys ); keys.use("ARG"); keys.reset_style("ARG","optional");
      64          28 :   keys.add("compulsory","METRIC","EUCLIDEAN","the method that you are going to use to measure the distances between points");
      65          28 :   keys.add("atoms","ATOMS","the list of atoms that you are going to use in the measure of distance that you are using");
      66          14 : }
      67             : 
      68          13 : EuclideanDissimilarityMatrix::EuclideanDissimilarityMatrix( const ActionOptions& ao ):
      69             :   Action(ao),
      70          13 :   AnalysisBase(ao)
      71             : {
      72          26 :   parse("METRIC",mtype); std::vector<AtomNumber> atoms;
      73          13 :   if( my_input_data->getNumberOfAtoms()>0 ) {
      74           6 :     parseAtomList("ATOMS",atoms);
      75           3 :     if( atoms.size()!=0 ) {
      76           0 :       mypdb.setAtomNumbers( atoms );
      77           0 :       for(unsigned i=0; i<atoms.size(); ++i) {
      78             :         bool found=false;
      79           0 :         for(unsigned j=0; j<my_input_data->getAtomIndexes().size(); ++j) {
      80           0 :           if( my_input_data->getAtomIndexes()[j]==atoms[i] ) { found=true; break; }
      81             :         }
      82           0 :         if( !found ) {
      83           0 :           std::string num; Tools::convert( atoms[i].serial(), num );
      84           0 :           error("atom number " + num + " is not stored in any action that has been input");
      85             :         }
      86             :       }
      87           0 :       mypdb.addBlockEnd( atoms.size() );
      88           3 :     } else if( getNumberOfArguments()==0 ) {
      89           1 :       mypdb.setAtomNumbers( my_input_data->getAtomIndexes() );
      90           1 :       mypdb.addBlockEnd( my_input_data->getAtomIndexes().size() );
      91           1 :       if( mtype=="EUCLIDEAN" ) mtype="OPTIMAL";
      92             :     }
      93             :   }
      94          13 :   log.printf("  measuring distances using %s metric \n",mtype.c_str() );
      95          13 :   if( my_input_data->getArgumentNames().size()>0 ) {
      96          12 :     if( getNumberOfArguments()==0 && atoms.size()==0 ) {
      97          10 :       std::vector<std::string> argnames( my_input_data->getArgumentNames() );
      98          10 :       mypdb.setArgumentNames( argnames ); requestArguments( my_input_data->getArgumentList() );
      99          10 :     } else {
     100           2 :       std::vector<Value*> myargs( getArguments() );
     101           2 :       std::vector<std::string> inargnames( my_input_data->getArgumentNames() );
     102           2 :       std::vector<std::string> argnames( myargs.size() );
     103           6 :       for(unsigned i=0; i<myargs.size(); ++i) {
     104           4 :         argnames[i]=myargs[i]->getName();
     105             :         bool found=false;
     106           6 :         for(unsigned j=0; j<inargnames.size(); ++j) {
     107           6 :           if( argnames[i]==inargnames[j] ) { found=true; break; }
     108             :         }
     109           4 :         if( !found ) error("input named " + my_input_data->getLabel() + " does not store/calculate quantity named " + argnames[i] );
     110             :       }
     111           2 :       mypdb.setArgumentNames( argnames ); requestArguments( myargs );
     112           2 :     }
     113             :   }
     114          13 : }
     115             : 
     116          20 : void EuclideanDissimilarityMatrix::performAnalysis() {
     117             :   // Resize dissimilarities matrix and set all elements to zero
     118          20 :   if( !usingLowMem() ) {
     119          20 :     dissimilarities.resize( getNumberOfDataPoints(), getNumberOfDataPoints() ); dissimilarities=0;
     120             :   }
     121          20 : }
     122             : 
     123         413 : std::string EuclideanDissimilarityMatrix::getDissimilarityInstruction() const {
     124         413 :   return "TYPE=" + mtype;
     125             : }
     126             : 
     127      611659 : double EuclideanDissimilarityMatrix::getDissimilarity( const unsigned& iframe, const unsigned& jframe ) {
     128             :   plumed_dbg_assert( iframe<getNumberOfDataPoints() && jframe<getNumberOfDataPoints() );
     129      611659 :   if( !usingLowMem() ) {
     130      611659 :     if( dissimilarities(iframe,jframe)>0. ) { return dissimilarities(iframe,jframe); }
     131             :   }
     132      256165 :   if( iframe!=jframe ) {
     133             :     double dd;
     134      255613 :     getStoredData( iframe, true ).transferDataToPDB( mypdb );
     135      255613 :     auto myref1=metricRegister().create<ReferenceConfiguration>(mtype, mypdb);
     136      255613 :     getStoredData( jframe, true ).transferDataToPDB( mypdb );
     137      255613 :     auto myref2=metricRegister().create<ReferenceConfiguration>(mtype, mypdb);
     138      255613 :     if( !usingLowMem() ) dd=dissimilarities(iframe,jframe) = dissimilarities(jframe,iframe) = distance( getPbc(), getArguments(), myref1.get(), myref2.get(), true );
     139           0 :     else dd=distance( getPbc(), getArguments(), myref1.get(), myref2.get(), true );
     140             :     return dd;
     141      255613 :   }
     142             :   return 0.0;
     143             : }
     144             : 
     145             : }
     146             : }