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             : #include "ActionVolume.h"
      23             : 
      24             : namespace PLMD {
      25             : namespace multicolvar {
      26             : 
      27          34 : void ActionVolume::registerKeywords( Keywords& keys ) {
      28          34 :   VolumeGradientBase::registerKeywords( keys );
      29         102 :   if( keys.reserved("VMEAN") ) keys.use("VMEAN");
      30         102 :   keys.use("MEAN"); keys.use("LESS_THAN"); keys.use("MORE_THAN");
      31         102 :   keys.use("BETWEEN"); keys.use("HISTOGRAM"); keys.use("SUM");
      32          68 :   keys.add("compulsory","SIGMA","the width of the function to be used for kernel density estimation");
      33          68 :   keys.add("compulsory","KERNEL","gaussian","the type of kernel function to be used");
      34          68 :   keys.addFlag("OUTSIDE",false,"calculate quantities for colvars that are on atoms outside the region of interest");
      35          34 : }
      36             : 
      37          28 : ActionVolume::ActionVolume(const ActionOptions&ao):
      38             :   Action(ao),
      39          28 :   VolumeGradientBase(ao)
      40             : {
      41             :   // Find number of quantities
      42          28 :   if( getPntrToMultiColvar()->isDensity() ) nquantities=2;                           // Value + weight
      43          17 :   else if( getPntrToMultiColvar()->getNumberOfQuantities()==2 ) nquantities=2;       // Value + weight
      44           2 :   else nquantities = 1 + getPntrToMultiColvar()->getNumberOfQuantities()-2 + 1;      // Norm  + vector + weight
      45             : 
      46             :   // Output some nice information
      47          28 :   std::string functype=getPntrToMultiColvar()->getName();
      48         279 :   std::transform( functype.begin(), functype.end(), functype.begin(), [](unsigned char c) { return std::tolower(c); } );
      49          28 :   log.printf("  calculating %s inside region of insterest\n",functype.c_str() );
      50             : 
      51          28 :   parseFlag("OUTSIDE",not_in); sigma=0.0;
      52          81 :   if( keywords.exists("SIGMA") ) parse("SIGMA",sigma);
      53          84 :   if( keywords.exists("KERNEL") ) parse("KERNEL",kerneltype);
      54             : 
      55          28 :   if( getPntrToMultiColvar()->isDensity() ) {
      56             :     std::string input;
      57          22 :     addVessel( "SUM", input, -1 );  // -1 here means that this value will be named getLabel()
      58             :   }
      59          28 :   readVesselKeywords();
      60          28 : }
      61             : 
      62       69358 : void ActionVolume::calculateAllVolumes( const unsigned& curr, MultiValue& outvals ) const {
      63       69358 :   Vector catom_pos=getPntrToMultiColvar()->getCentralAtomPos( curr );
      64             : 
      65       69358 :   double weight; Vector wdf; Tensor vir; std::vector<Vector> refders( getNumberOfAtoms() );
      66       69358 :   weight=calculateNumberInside( catom_pos, wdf, vir, refders );
      67       69358 :   if( not_in ) {
      68        4000 :     weight = 1.0 - weight; wdf *= -1.; vir *=-1;
      69        8000 :     for(unsigned i=0; i<refders.size(); ++i) refders[i]*=-1;
      70             :   }
      71       69358 :   setNumberInVolume( 0, curr, weight, wdf, vir, refders, outvals );
      72       69358 : }
      73             : 
      74      148129 : bool ActionVolume::inVolumeOfInterest( const unsigned& curr ) const {
      75      148129 :   Vector catom_pos=getPntrToMultiColvar()->getCentralAtomPos( curr );
      76      148129 :   Vector wdf; Tensor vir; std::vector<Vector> refders( getNumberOfAtoms() );
      77      148129 :   double weight=calculateNumberInside( catom_pos, wdf, vir, refders );
      78      148129 :   if( not_in ) weight = 1.0 - weight;
      79      148129 :   if( weight<getTolerance() ) return false;
      80             :   return true;
      81             : }
      82             : 
      83             : }
      84             : }