Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2011-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             : #ifndef __PLUMED_matrixtools_MatrixDissimilarities_h
      23             : #define __PLUMED_matrixtools_MatrixDissimilarities_h
      24             : #include "MatrixTimesMatrix.h"
      25             : 
      26             : namespace PLMD {
      27             : namespace matrixtools {
      28             : struct Dissimilarities {
      29             :   typedef void isDissimilarities;
      30             :   bool squared{false};
      31             :   bool periodic{false};
      32             :   double min{0};
      33             :   double max{0};
      34             :   double max_minus_min{0};
      35             :   double inv_max_minus_min{0};
      36             :   static void registerKeywords( Keywords& keys );
      37             :   template <typename PTM>
      38             :   void setup( MatrixTimesMatrix<Dissimilarities,PTM>* action,
      39             :               const Value* myval );
      40             :   static void calculate( bool noderiv,
      41             :                          const Dissimilarities& actdata,
      42             :                          const InputVectors& vectors,
      43             :                          MatrixElementOutput& output );
      44             : #ifdef __PLUMED_HAS_OPENACC
      45             :   void toACCDevice() const {
      46             : #pragma acc enter data copyin(this[0:1], squared, periodic, min, max, \
      47             :                               max_minus_min, inv_max_minus_min)
      48             : 
      49             :   }
      50             :   void removeFromACCDevice() const {
      51             : #pragma acc exit data delete(inv_max_minus_min, max_minus_min, max, \
      52             :                              min, periodic, squared, this[0:1])
      53             :   }
      54             : #endif //__PLUMED_HAS_OPENACC
      55             : };
      56             : 
      57             : 
      58          51 : void Dissimilarities::registerKeywords( Keywords& keys ) {
      59         102 :   keys.setDisplayName("DISSIMILARITY");
      60          51 :   keys.addFlag("SQUARED",false,"calculate the squares of the dissimilarities (this option cannot be used with MATRIX_PRODUCT)");
      61         102 :   keys.setValueDescription("matrix","the dissimilarity matrix");
      62          51 : }
      63             : 
      64             : template <typename PTM>
      65          12 : void Dissimilarities::setup( MatrixTimesMatrix<Dissimilarities, PTM>* action,
      66             :                              const Value* myval ) {
      67          12 :   action->parseFlag("SQUARED",squared);
      68          12 :   if( squared ) {
      69           7 :     action->log.printf("  calculating the squares of the dissimilarities \n");
      70             :   }
      71          12 :   periodic = myval->isPeriodic();
      72          12 :   if( periodic ) {
      73             :     std::string strmin, strmax;
      74           3 :     myval->getDomain( strmin, strmax );
      75           3 :     Tools::convert( strmin, min );
      76           3 :     Tools::convert( strmax, max );
      77           3 :     max_minus_min = max - min;
      78           3 :     plumed_assert( max_minus_min>0 );
      79           3 :     inv_max_minus_min=1.0/max_minus_min;
      80             :   }
      81          12 : }
      82             : 
      83      545753 : void Dissimilarities::calculate( const bool noderiv,
      84             :                                  const Dissimilarities& actdata,
      85             :                                  const InputVectors& vectors,
      86             :                                  MatrixElementOutput& output ) {
      87             :   // same layout of the derivatives as in MatrixProduct::calculate
      88      545753 :   const std::size_t pp = vectors.nelem;
      89     2175346 :   for(unsigned i=0; i<vectors.nelem; ++i) {
      90     1629593 :     double tmp1 = vectors.arg1[i] - vectors.arg2[i];
      91     1629593 :     if( actdata.periodic ) {
      92     1530000 :       tmp1 = tmp1*actdata.inv_max_minus_min;
      93     1530000 :       tmp1 = Tools::pbc(tmp1);
      94     1530000 :       tmp1 = tmp1*actdata.max_minus_min;
      95             :     }
      96     1629593 :     output.values[0] += tmp1*tmp1;
      97     1629593 :     output.derivs[0][i] = 2*tmp1;
      98     1629593 :     output.derivs[0][pp+i] = -2*tmp1;
      99             :   }
     100      545753 :   if( actdata.squared ) {
     101             :     return ;
     102             :   }
     103         388 :   output.values[0] = sqrt( output.values[0] );
     104             : 
     105        1151 :   for(unsigned i=0; i<vectors.nelem; ++i) {
     106         763 :     output.derivs[0][i] = output.derivs[0][i] / (2*output.values[0]);
     107         763 :     output.derivs[0][pp+i] = output.derivs[0][pp+i] / (2*output.values[0]);
     108             :   }
     109             : }
     110             : } //namespace matrixtools
     111             : } //PLMD
     112             : #endif //__PLUMED_matrixtools_MatrixDissimilarities_h