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             : #ifndef __PLUMED_adjmat_ContactMatrix_h
      23             : #define __PLUMED_adjmat_ContactMatrix_h
      24             : 
      25             : #include "AdjacencyMatrixBase.h"
      26             : #include "tools/SwitchingFunction.h"
      27             : namespace PLMD {
      28             : namespace adjmat {
      29             : 
      30             : struct ContactMatrix {
      31             : #ifdef __PLUMED_HAS_OPENACC
      32             :   SwitchingFunctionAccelerable switchingFunction;
      33             : #else
      34             :   SwitchingFunction switchingFunction;
      35             : #endif
      36             :   static void registerKeywords( Keywords& keys );
      37             :   template <typename myPTM>
      38             :   void parseInput( AdjacencyMatrixBase<ContactMatrix,myPTM>* action );
      39             :   static void calculateWeight( const ContactMatrix& data,
      40             :                                const AdjacencyMatrixInput& input,
      41             :                                MatrixOutput output );
      42             : #ifdef __PLUMED_HAS_OPENACC
      43             :   void toACCDevice() const {
      44             : #pragma acc enter data copyin(this[0:1])
      45             :     switchingFunction.toACCDevice();
      46             :   }
      47             :   void removeFromACCDevice() const {
      48             :     switchingFunction.removeFromACCDevice();
      49             : #pragma acc exit data delete(this[0:1])
      50             :   }
      51             : #endif //__PLUMED_HAS_OPENACC
      52             : };
      53             : 
      54         998 : void ContactMatrix::registerKeywords( Keywords& keys ) {
      55        1996 :   keys.setDisplayName("CONTACT_MATRIX");
      56         998 :   keys.add("compulsory","NN","6","The n parameter of the switching function ");
      57         998 :   keys.add("compulsory","MM","0","The m parameter of the switching function; 0 implies 2*NN");
      58         998 :   keys.add("compulsory","D_0","0.0","The d_0 parameter of the switching function");
      59         998 :   keys.add("compulsory","R_0","The r_0 parameter of the switching function");
      60         998 :   keys.add("optional","SWITCH","This keyword is used if you want to employ an alternative to the continuous switching function defined above. "
      61             :            "The following provides information on the switching functions that are available. "
      62             :            "When this keyword is present you no longer need the NN, MM, D_0 and R_0 keywords.");
      63        1996 :   keys.linkActionInDocs("SWITCH","LESS_THAN");
      64         998 : }
      65             : 
      66             : template <typename myPTM>
      67         278 : void ContactMatrix::parseInput( AdjacencyMatrixBase<ContactMatrix,myPTM>* action ) {
      68             :   std::string errors;
      69             :   std::string swinput;
      70         556 :   action->parse("SWITCH",swinput);
      71         278 :   if( swinput.length()>0 ) {
      72         250 :     switchingFunction.set( swinput, errors );
      73         250 :     if( errors.length()!=0 ) {
      74           0 :       action->error("problem reading switching function description " + errors);
      75             :     }
      76             :   } else {
      77          28 :     int nn=0;
      78          28 :     int mm=0;
      79          28 :     double r_0=-1.0;
      80          28 :     double d_0= 0.0;
      81          28 :     action->parse("NN",nn);
      82          28 :     action->parse("MM",mm);
      83          28 :     action->parse("R_0",r_0);
      84          28 :     action->parse("D_0",d_0);
      85          28 :     if( r_0<0.0 ) {
      86           0 :       action->error("you must set a value for R_0");
      87             :     }
      88          28 :     switchingFunction.set(nn,mm,r_0,d_0);
      89             :   }
      90             :   // And set the link cell cutoff
      91         278 :   action->log.printf("  switching function cutoff is %s \n",switchingFunction.description().c_str() );
      92         278 :   action->setLinkCellCutoff( true, switchingFunction.get_dmax(), -1.0, true );
      93         278 : }
      94             : 
      95    14314724 : void ContactMatrix::calculateWeight( const ContactMatrix& data,
      96             :                                      const AdjacencyMatrixInput& input,
      97             :                                      MatrixOutput output ) {
      98             :   const double mod2 = input.pos.modulo2();
      99    14314724 :   if( mod2<epsilon ) {
     100        1871 :     output.val[0]=0;
     101    10119487 :     return;  // Atoms can't be bonded to themselves
     102             :   }
     103             :   double dfunc;
     104    14312853 :   output.val[0] = data.switchingFunction.calculateSqr( mod2, dfunc );
     105    14312853 :   if( fabs(output.val[0])<epsilon ) {
     106      838040 :     output.val[0] = 0.0;
     107      838040 :     return;
     108             :   }
     109    13474813 :   if( input.noderiv ) {
     110             :     return;
     111             :   }
     112     4195237 :   const Vector v { (-dfunc)*input.pos[0],
     113             :                    (-dfunc)*input.pos[1],
     114     4195237 :                    (-dfunc)*input.pos[2] };
     115     4195237 :   output.deriv[0] = v[0];
     116     4195237 :   output.deriv[1] = v[1];
     117     4195237 :   output.deriv[2] = v[2];
     118     4195237 :   output.deriv[3] =-v[0];
     119     4195237 :   output.deriv[4] =-v[1];
     120     4195237 :   output.deriv[5] =-v[2];
     121             : 
     122     4195237 :   output.assignOuterProduct( 6, v, input.pos);
     123             : }
     124             : }
     125             : }
     126             : #endif