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 "core/ActionRegister.h"
      23             : #include "vesselbase/StoreDataVessel.h"
      24             : #include "ContourFindingBase.h"
      25             : #include "core/PlumedMain.h"
      26             : #include "core/Atoms.h"
      27             : 
      28             : //+PLUMEDOC GRIDANALYSIS FIND_CONTOUR
      29             : /*
      30             : Find an isocontour in a smooth function.
      31             : 
      32             : As discussed in the part of the manual on \ref Analysis PLUMED contains a number of tools that allow you to calculate
      33             : a function on a grid.  The function on this grid might be a \ref HISTOGRAM as a function of a few collective variables
      34             : or it might be a phase field that has been calculated using \ref MULTICOLVARDENS.  If this function has one or two input
      35             : arguments it is relatively straightforward to plot the function.  If by contrast the data has a three or more dimensions
      36             : it can be difficult to visualize.
      37             : 
      38             : This action provides one tool for visualizing these functions.  It can be used to search for a set of points on a contour
      39             : where the function takes a particular values.  In other words, for the function \f$f(x,y)\f$ this action would find a set
      40             : of points \f$\{x_c,y_c\}\f$ that have:
      41             : 
      42             : \f[
      43             : f(x_c,y_c) - c = 0
      44             : \f]
      45             : 
      46             : where \f$c\f$ is some constant value that is specified by the user.  The points on this contour are detected using a variant
      47             : on the marching squares or marching cubes algorithm, which you can find information on here:
      48             : 
      49             : https://en.wikipedia.org/wiki/Marching_squares
      50             : https://en.wikipedia.org/wiki/Marching_cubes
      51             : 
      52             : As such, and unlike \ref FIND_CONTOUR_SURFACE or \ref FIND_SPHERICAL_CONTOUR, the function input to this action can have any dimension.
      53             : Furthermore, the topology of the contour will be determined by the algorithm and does not need to be specified by the user.
      54             : 
      55             : \par Examples
      56             : 
      57             : The input below allows you to calculate something akin to a Willard-Chandler dividing surface \cite wcsurface.
      58             : The simulation cell in this case contains a solid phase and a liquid phase.  The Willard-Chandler surface is the
      59             : surface that separates the parts of the box containing the solid from the parts containing the liquid.  To compute the position
      60             : of this surface  the \ref FCCUBIC symmetry function is calculated for each of the atoms in the system from on the geometry of the
      61             : atoms in the first coordination sphere of each of the atoms.  These quantities are then transformed using a switching function.
      62             : This procedure generates a single number for each atom in the system and this quantity has a value of one for atoms that are in
      63             : parts of the box that resemble the solid structure and zero for atoms that are in parts of the box that resemble the liquid.
      64             : The position of a virtual atom is then computed using \ref CENTER_OF_MULTICOLVAR and a phase field model is constructed using
      65             : \ref MULTICOLVARDENS.  These procedure ensures that we have a continuous function that gives a measure of the average degree of
      66             : solidness at each point in the simulation cell.  The Willard-Chandler dividing surface is calculated by finding a a set of points
      67             : at which the value of this phase field is equal to 0.5.  This set of points is output to file called mycontour.dat.  A new contour
      68             : is found on every single step for each frame that is read in.
      69             : 
      70             : \plumedfile
      71             : UNITS NATURAL
      72             : FCCUBIC ...
      73             :   SPECIES=1-96000 SWITCH={CUBIC D_0=1.2 D_MAX=1.5}
      74             :   ALPHA=27 PHI=0.0 THETA=-1.5708 PSI=-2.35619 LABEL=fcc
      75             : ... FCCUBIC
      76             : 
      77             : tfcc: MTRANSFORM_MORE DATA=fcc LOWMEM SWITCH={SMAP R_0=0.5 A=8 B=8}
      78             : center: CENTER_OF_MULTICOLVAR DATA=tfcc
      79             : 
      80             : dens: MULTICOLVARDENS ...
      81             :   DATA=tfcc ORIGIN=center DIR=xyz
      82             :   NBINS=80,80,80 BANDWIDTH=1.0,1.0,1.0 STRIDE=1 CLEAR=1
      83             : ...
      84             : 
      85             : FIND_CONTOUR GRID=dens CONTOUR=0.5 FILE=mycontour.xyz
      86             : \endplumedfile
      87             : 
      88             : */
      89             : //+ENDPLUMEDOC
      90             : 
      91             : namespace PLMD {
      92             : namespace gridtools {
      93             : 
      94             : class FindContour : public ContourFindingBase {
      95             : private:
      96             :   bool firsttime;
      97             :   unsigned gbuffer;
      98             : /// Stuff for output
      99             :   OFile of;
     100             :   double lenunit;
     101             :   std::string fmt_xyz;
     102             : public:
     103             :   static void registerKeywords( Keywords& keys );
     104             :   explicit FindContour(const ActionOptions&ao);
     105           0 :   bool checkAllActive() const override { return gbuffer==0; }
     106             :   void prepareForAveraging() override;
     107           0 :   bool isPeriodic() override { return false; }
     108           7 :   unsigned getNumberOfQuantities() const override { return 1 + ingrid->getDimension(); }
     109             :   void compute( const unsigned& current, MultiValue& myvals ) const override;
     110             :   void finishAveraging() override;
     111             : };
     112             : 
     113       10421 : PLUMED_REGISTER_ACTION(FindContour,"FIND_CONTOUR")
     114             : 
     115           2 : void FindContour::registerKeywords( Keywords& keys ) {
     116           2 :   ContourFindingBase::registerKeywords( keys );
     117             : // We want a better way of doing this bit
     118           4 :   keys.add("compulsory","BUFFER","0","number of buffer grid points around location where grid was found on last step.  If this is zero the full grid is calculated on each step");
     119           4 :   keys.add("compulsory","FILE","file on which to output coordinates");
     120           4 :   keys.add("compulsory","UNITS","PLUMED","the units in which to print out the coordinates. PLUMED means internal PLUMED units");
     121           4 :   keys.add("optional", "PRECISION","The number of digits in trajectory file");
     122           2 : }
     123             : 
     124           1 : FindContour::FindContour(const ActionOptions&ao):
     125             :   Action(ao),
     126             :   ContourFindingBase(ao),
     127           1 :   firsttime(true)
     128             : {
     129             : 
     130           1 :   parse("BUFFER",gbuffer);
     131           1 :   if( gbuffer>0 ) log.printf("  after first step a subset of only %u grid points around where the countour was found will be checked\n",gbuffer);
     132             : 
     133           2 :   std::string file; parse("FILE",file);
     134           1 :   if( file.length()==0 ) error("name out output file was not specified");
     135           1 :   std::string type=Tools::extension(file);
     136           1 :   log<<"  file name "<<file<<"\n";
     137           1 :   if(type!="xyz") error("can only print xyz file type with contour finding");
     138             : 
     139             :   fmt_xyz="%f";
     140           2 :   std::string precision; parse("PRECISION",precision);
     141           1 :   if(precision.length()>0) {
     142           1 :     int p; Tools::convert(precision,p);
     143           1 :     log<<"  with precision "<<p<<"\n";
     144             :     std::string a,b;
     145           1 :     Tools::convert(p+5,a);
     146           1 :     Tools::convert(p,b);
     147           2 :     fmt_xyz="%"+a+"."+b+"f";
     148             :   }
     149           2 :   std::string unitname; parse("UNITS",unitname);
     150           1 :   if(unitname!="PLUMED") {
     151           0 :     Units myunit; myunit.setLength(unitname);
     152           0 :     lenunit=plumed.getAtoms().getUnits().getLength()/myunit.getLength();
     153           0 :   }
     154           1 :   else lenunit=1.0;
     155           1 :   of.link(*this); of.open(file);
     156           1 :   checkRead(); mydata=buildDataStashes( NULL );
     157           1 : }
     158             : 
     159           2 : void FindContour::prepareForAveraging() {
     160             :   // Create a task list if first time
     161           2 :   if( firsttime ) {
     162       16465 :     for(unsigned i=0; i<ingrid->getDimension()*ingrid->getNumberOfPoints(); ++i) addTaskToList( i );
     163             :   }
     164           2 :   firsttime=false; deactivateAllTasks();
     165             : 
     166             :   // We now need to identify the grid points that we need to search through
     167           2 :   std::vector<unsigned> nbin( ingrid->getNbin() );
     168           2 :   std::vector<unsigned> ind( ingrid->getDimension() );
     169           2 :   std::vector<unsigned> ones( ingrid->getDimension(), 1 );
     170             :   unsigned num_neighbours; std::vector<unsigned> neighbours;
     171       10978 :   for(unsigned i=0; i<ingrid->getNumberOfPoints(); ++i) {
     172             :     // Ensure inactive grid points are ignored
     173       10976 :     if( ingrid->inactive(i) ) continue;
     174             : 
     175             :     // Get the index of the current grid point
     176       10976 :     ingrid->getIndices( i, ind );
     177       10976 :     ingrid->getNeighbors( ind, ones, num_neighbours, neighbours );
     178             :     bool cycle=false;
     179      307328 :     for(unsigned j=0; j<num_neighbours; ++j) {
     180      296352 :       if( ingrid->inactive( neighbours[j]) ) { cycle=true; break; }
     181             :     }
     182       10976 :     if( cycle ) continue;
     183             : 
     184             :     // Get the value of a point on the grid
     185       10976 :     double val1=getFunctionValue( i ) - contour;
     186             :     bool edge=false;
     187       43904 :     for(unsigned j=0; j<ingrid->getDimension(); ++j) {
     188             :       // Make sure we don't search at the edge of the grid
     189       32928 :       if( !ingrid->isPeriodic(j) && (ind[j]+1)==nbin[j] ) continue;
     190       32928 :       else if( (ind[j]+1)==nbin[j] ) { edge=true; ind[j]=0; }
     191       30968 :       else ind[j]+=1;
     192       32928 :       double val2=getFunctionValue( ind ) - contour;
     193       32928 :       if( val1*val2<0 ) taskFlags[ ingrid->getDimension()*i + j ] = 1;
     194       32928 :       if( ingrid->isPeriodic(j) && edge ) { edge=false; ind[j]=nbin[j]-1; }
     195       30968 :       else ind[j]-=1;
     196             :     }
     197             :   }
     198           2 :   lockContributors();
     199           2 : }
     200             : 
     201         554 : void FindContour::compute( const unsigned& current, MultiValue& myvals ) const {
     202             :   // Retrieve the initial grid point coordinates
     203         554 :   unsigned gpoint = std::floor( current / ingrid->getDimension() );
     204         554 :   std::vector<double> point( ingrid->getDimension() );
     205         554 :   ingrid->getGridPointCoordinates( gpoint, point );
     206             : 
     207             :   // Retrieve the direction we are searching for the contour
     208         554 :   unsigned gdir = current%(ingrid->getDimension() );
     209         554 :   std::vector<double> direction( ingrid->getDimension(), 0 );
     210         554 :   direction[gdir] = 0.999999999*ingrid->getGridSpacing()[gdir];
     211             : 
     212             :   // Now find the contour
     213             :   findContour( direction, point );
     214             :   // And transfer to the store data vessel
     215        2216 :   for(unsigned i=0; i<ingrid->getDimension(); ++i) myvals.setValue( 1+i, point[i] );
     216         554 : }
     217             : 
     218           1 : void FindContour::finishAveraging() {
     219             :   // And update the list of active grid points
     220           1 :   if( gbuffer>0 ) {
     221             :     std::vector<unsigned> neighbours; unsigned num_neighbours;
     222           0 :     std::vector<unsigned> ugrid_indices( ingrid->getDimension() );
     223           0 :     std::vector<bool> active( ingrid->getNumberOfPoints(), false );
     224           0 :     std::vector<unsigned> gbuffer_vec( ingrid->getDimension(), gbuffer );
     225           0 :     for(unsigned i=0; i<getCurrentNumberOfActiveTasks(); ++i) {
     226             :       // Get the point we are operating on
     227           0 :       unsigned ipoint = std::floor( getActiveTask(i) / ingrid->getDimension() );
     228             :       // Get the indices of this point
     229           0 :       ingrid->getIndices( ipoint, ugrid_indices );
     230             :       // Now activate buffer region
     231           0 :       ingrid->getNeighbors( ugrid_indices, gbuffer_vec, num_neighbours, neighbours );
     232           0 :       for(unsigned n=0; n<num_neighbours; ++n) active[ neighbours[n] ]=true;
     233             :     }
     234           0 :     ingrid->activateThesePoints( active );
     235             :   }
     236           1 :   std::vector<double> point( 1 + ingrid->getDimension() );
     237           1 :   of.printf("%u\n",mydata->getNumberOfStoredValues());
     238           1 :   of.printf("Points found on isocontour\n");
     239         555 :   for(unsigned i=0; i<mydata->getNumberOfStoredValues(); ++i) {
     240         554 :     mydata->retrieveSequentialValue( i, false, point ); of.printf("X");
     241        2216 :     for(unsigned j=0; j<ingrid->getDimension(); ++j) of.printf( (" " + fmt_xyz).c_str(), lenunit*point[1+j] );
     242         554 :     of.printf("\n");
     243             :   }
     244           1 : }
     245             : 
     246             : }
     247             : }