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 "ClusterAnalysisBase.h"
      23             : #include "AdjacencyMatrixVessel.h"
      24             : #include "core/ActionRegister.h"
      25             : #include "tools/SwitchingFunction.h"
      26             : 
      27             : //+PLUMEDOC CONCOMP CLUSTER_DISTRIBUTION
      28             : /*
      29             : Calculate functions of the distribution of properties in your connected components.
      30             : 
      31             : This collective variable was developed for looking at nucleation phenomena, where you are
      32             : interested in using studying the behavior of atoms in small aggregates or nuclei.  In these sorts of
      33             : problems you might be interested in the distribution of the sizes of the clusters in your system.
      34             : A detailed description of this CV can be found in \cite tribello-clustering.
      35             : 
      36             : \par Examples
      37             : 
      38             : The input provided below calculates the local q6 Steinhardt parameter on each atom.  The coordination number
      39             : that atoms with a high value for the local q6 Steinhardt parameter have with other atoms that have a high
      40             : value for the local q6 Steinhardt parameter is then computed.  A contact matrix is then computed that measures
      41             : whether atoms atoms \f$i\f$ and \f$j\f$ have a high value for this coordination number and if they are within
      42             : 3.6 nm of each other.  The connected components of this matrix are then found using a depth first clustering
      43             : algorithm on the corresponding graph. The number of components in this graph that contain more than 27 atoms is then computed.
      44             : As discussed in \cite tribello-clustering an input similar to this one was used to analyze the formation of a polycrystal of GeTe from amorphous
      45             : GeTe.
      46             : 
      47             : \plumedfile
      48             : q6: Q6 SPECIES=1-300 SWITCH={GAUSSIAN D_0=5.29 R_0=0.01 D_MAX=5.3} LOWMEM
      49             : lq6: LOCAL_Q6 SPECIES=q6 SWITCH={GAUSSIAN D_0=5.29 R_0=0.01 D_MAX=5.3} LOWMEM
      50             : flq6: MFILTER_MORE DATA=lq6 SWITCH={GAUSSIAN D_0=0.19 R_0=0.01 D_MAX=0.2}
      51             : cc: COORDINATIONNUMBER SPECIES=flq6 SWITCH={GAUSSIAN D_0=3.59 R_0=0.01 D_MAX=3.6}
      52             : fcc: MFILTER_MORE DATA=cc SWITCH={GAUSSIAN D_0=5.99 R_0=0.01 D_MAX=6.0}
      53             : mat: CONTACT_MATRIX ATOMS=fcc SWITCH={GAUSSIAN D_0=3.59 R_0=0.01 D_MAX=3.6}
      54             : dfs: DFSCLUSTERING MATRIX=mat
      55             : nclust: CLUSTER_DISTRIBUTION CLUSTERS=dfs TRANSFORM={GAUSSIAN D_0=5.99 R_0=0.01 D_MAX=6.0} MORE_THAN={GAUSSIAN D_0=26.99 R_0=0.01 D_MAX=27}
      56             : PRINT ARG=nclust.* FILE=colvar
      57             : \endplumedfile
      58             : 
      59             : */
      60             : //+ENDPLUMEDOC
      61             : 
      62             : namespace PLMD {
      63             : namespace adjmat {
      64             : 
      65             : class ClusterDistribution : public ClusterAnalysisBase {
      66             : private:
      67             :   unsigned nderivatives;
      68             : ///
      69             :   bool use_switch, inverse;
      70             : //
      71             :   SwitchingFunction sf;
      72             : public:
      73             : /// Create manual
      74             :   static void registerKeywords( Keywords& keys );
      75             : /// Constructor
      76             :   explicit ClusterDistribution(const ActionOptions&);
      77             : /// Do the calculation
      78             :   void calculate() override;
      79             : /// We can use ActionWithVessel to run all the calculation
      80             :   void performTask( const unsigned&, const unsigned&, MultiValue& ) const override;
      81             : };
      82             : 
      83       13787 : PLUMED_REGISTER_ACTION(ClusterDistribution,"CLUSTER_DISTRIBUTION")
      84             : 
      85           5 : void ClusterDistribution::registerKeywords( Keywords& keys ) {
      86           5 :   ClusterAnalysisBase::registerKeywords( keys );
      87          10 :   keys.add("compulsory","TRANSFORM","none","the switching function to use to convert the crystallinity parameter to a number between zero and one");
      88          10 :   keys.addFlag("INVERSE_TRANSFORM",false,"when TRANSFORM appears alone the input symmetry functions, \\f$x\\f$ are transformed used \\f$1-s(x)\\f$ "
      89             :                "where \\f$s(x)\\f$ is a switching function.  When this option is used you instead transform using \\f$s(x)\\f$ only.");
      90           5 :   keys.use("MORE_THAN");
      91           5 :   keys.use("LESS_THAN");
      92           5 :   keys.use("BETWEEN");
      93           5 :   keys.use("HISTOGRAM");
      94           5 :   keys.use("ALT_MIN");
      95           5 :   keys.use("MIN");
      96           5 :   keys.use("MAX");
      97           5 : }
      98             : 
      99           1 : ClusterDistribution::ClusterDistribution(const ActionOptions&ao):
     100             :   Action(ao),
     101             :   ClusterAnalysisBase(ao),
     102           1 :   nderivatives(0) {
     103           1 :   use_switch=false;
     104             :   std::string input, errors;
     105           2 :   parse("TRANSFORM",input);
     106           1 :   if( input!="none" ) {
     107           0 :     use_switch=true;
     108           0 :     sf.set( input, errors );
     109           0 :     if( errors.length()!=0 ) {
     110           0 :       error("problem reading SWITCH keyword : " + errors );
     111             :     }
     112             :   }
     113           1 :   parseFlag("INVERSE_TRANSFORM",inverse);
     114           1 :   if( inverse && !use_switch ) {
     115           0 :     error("INVERSE_TRANSFORM option was specified but no TRANSOFRM switching function was given");
     116             :   }
     117             : 
     118             :   // Create all tasks by copying those from underlying DFS object (which is actually MultiColvar)
     119           7 :   for(unsigned i=0; i<getNumberOfNodes(); ++i) {
     120           6 :     addTaskToList(i);
     121             :   }
     122             : 
     123             :   // And now finish the setup of everything in the base
     124             :   std::vector<AtomNumber> fake_atoms;
     125           1 :   setupMultiColvarBase( fake_atoms );
     126           1 : }
     127             : 
     128           1 : void ClusterDistribution::calculate() {
     129             :   // Activate the relevant tasks
     130           1 :   nderivatives = getNumberOfDerivatives();
     131           1 :   deactivateAllTasks();
     132           4 :   for(unsigned i=0; i<getNumberOfClusters(); ++i) {
     133           3 :     taskFlags[i]=1;
     134             :   }
     135           1 :   lockContributors();
     136             :   // Now do the calculation
     137           1 :   runAllTasks();
     138           1 : }
     139             : 
     140           3 : void ClusterDistribution::performTask( const unsigned& task_index, const unsigned& current, MultiValue& myvals ) const {
     141             :   std::vector<unsigned> myatoms;
     142           3 :   retrieveAtomsInCluster( current+1, myatoms );
     143             :   // This deals with filters
     144           3 :   if( myatoms.size()==1 && !nodeIsActive(myatoms[0]) ) {
     145             :     return ;
     146             :   }
     147             : 
     148           3 :   std::vector<double> vals( getNumberOfQuantities() );
     149           3 :   MultiValue tvals( getNumberOfQuantities(), nderivatives );
     150             : 
     151             :   // And this builds everything for this particular atom
     152             :   double vv, df, tval=0;
     153           9 :   for(unsigned j=0; j<myatoms.size(); ++j) {
     154           6 :     unsigned i=myatoms[j];
     155           6 :     getPropertiesOfNode( i, vals );
     156           6 :     if( use_switch && !inverse ) {
     157           0 :       vv = 1.0 - sf.calculate( vals[1], df );
     158           0 :       tval += vals[0]*vv;
     159           0 :       df=-df*vals[1];
     160           6 :     } else if( use_switch ) {
     161           0 :       vv = sf.calculate( vals[1], df );
     162           0 :       tval += vals[0]*vv;
     163           0 :       df=df*vals[1];
     164             :     } else {
     165           6 :       tval += vals[0]*vals[1];
     166           6 :       df=1.;
     167             :       vv=vals[1];
     168             :     }
     169           6 :     if( !doNotCalculateDerivatives() ) {
     170           6 :       getNodePropertyDerivatives( i, tvals );
     171         168 :       for(unsigned k=0; k<tvals.getNumberActive(); ++k) {
     172         162 :         unsigned kat=tvals.getActiveIndex(k);
     173         162 :         myvals.addDerivative( 1, kat, vals[0]*df*tvals.getDerivative(1,kat) + vv*tvals.getDerivative(0,kat) );
     174             :       }
     175           6 :       tvals.clearAll();
     176             :     }
     177             :   }
     178             :   myvals.setValue( 0, 1.0 );
     179             :   myvals.setValue( 1, tval );
     180           3 : }
     181             : 
     182             : }
     183             : }