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       10421 : PLUMED_REGISTER_ACTION(ClusterDistribution,"CLUSTER_DISTRIBUTION")
      84             : 
      85           2 : void ClusterDistribution::registerKeywords( Keywords& keys ) {
      86           2 :   ClusterAnalysisBase::registerKeywords( keys );
      87           4 :   keys.add("compulsory","TRANSFORM","none","the switching function to use to convert the crystallinity parameter to a number between zero and one");
      88           4 :   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           6 :   keys.use("MORE_THAN"); keys.use("LESS_THAN"); keys.use("BETWEEN");
      91           8 :   keys.use("HISTOGRAM"); keys.use("ALT_MIN"); keys.use("MIN"); keys.use("MAX");
      92           2 : }
      93             : 
      94           1 : ClusterDistribution::ClusterDistribution(const ActionOptions&ao):
      95             :   Action(ao),
      96             :   ClusterAnalysisBase(ao),
      97           1 :   nderivatives(0)
      98             : {
      99           1 :   use_switch=false;
     100           2 :   std::string input, errors; parse("TRANSFORM",input);
     101           1 :   if( input!="none" ) {
     102           0 :     use_switch=true; sf.set( input, errors );
     103           0 :     if( errors.length()!=0 ) error("problem reading SWITCH keyword : " + errors );
     104             :   }
     105           1 :   parseFlag("INVERSE_TRANSFORM",inverse);
     106           1 :   if( inverse && !use_switch ) error("INVERSE_TRANSFORM option was specified but no TRANSOFRM switching function was given");
     107             : 
     108             :   // Create all tasks by copying those from underlying DFS object (which is actually MultiColvar)
     109           7 :   for(unsigned i=0; i<getNumberOfNodes(); ++i) addTaskToList(i);
     110             : 
     111             :   // And now finish the setup of everything in the base
     112           1 :   std::vector<AtomNumber> fake_atoms; setupMultiColvarBase( fake_atoms );
     113           1 : }
     114             : 
     115           1 : void ClusterDistribution::calculate() {
     116             :   // Activate the relevant tasks
     117           1 :   nderivatives = getNumberOfDerivatives();
     118           1 :   deactivateAllTasks();
     119           4 :   for(unsigned i=0; i<getNumberOfClusters(); ++i) taskFlags[i]=1;
     120           1 :   lockContributors();
     121             :   // Now do the calculation
     122           1 :   runAllTasks();
     123           1 : }
     124             : 
     125           3 : void ClusterDistribution::performTask( const unsigned& task_index, const unsigned& current, MultiValue& myvals ) const {
     126           3 :   std::vector<unsigned> myatoms; retrieveAtomsInCluster( current+1, myatoms );
     127             :   // This deals with filters
     128           3 :   if( myatoms.size()==1 && !nodeIsActive(myatoms[0]) ) return ;
     129             : 
     130           3 :   std::vector<double> vals( getNumberOfQuantities() );
     131           3 :   MultiValue tvals( getNumberOfQuantities(), nderivatives );
     132             : 
     133             :   // And this builds everything for this particular atom
     134             :   double vv, df, tval=0;
     135           9 :   for(unsigned j=0; j<myatoms.size(); ++j) {
     136           6 :     unsigned i=myatoms[j];
     137           6 :     getPropertiesOfNode( i, vals );
     138           6 :     if( use_switch && !inverse ) {
     139           0 :       vv = 1.0 - sf.calculate( vals[1], df );
     140           0 :       tval += vals[0]*vv; df=-df*vals[1];
     141           6 :     } else if( use_switch ) {
     142           0 :       vv = sf.calculate( vals[1], df );
     143           0 :       tval += vals[0]*vv; df=df*vals[1];
     144             :     } else {
     145           6 :       tval += vals[0]*vals[1]; df=1.; vv=vals[1];
     146             :     }
     147           6 :     if( !doNotCalculateDerivatives() ) {
     148           6 :       getNodePropertyDerivatives( i, tvals );
     149         168 :       for(unsigned k=0; k<tvals.getNumberActive(); ++k) {
     150         162 :         unsigned kat=tvals.getActiveIndex(k);
     151         162 :         myvals.addDerivative( 1, kat, vals[0]*df*tvals.getDerivative(1,kat) + vv*tvals.getDerivative(0,kat) );
     152             :       }
     153           6 :       tvals.clearAll();
     154             :     }
     155             :   }
     156             :   myvals.setValue( 0, 1.0 ); myvals.setValue( 1, tval );
     157           3 : }
     158             : 
     159             : }
     160             : }