Line data    Source code

       1             : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
       2             :    Copyright (c) 2015-2020 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           3 : 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             : /// The data is stored in a grid
     103             :   vesselbase::StoreDataVessel* mydata;
     104             : public:
     105             :   static void registerKeywords( Keywords& keys );
     106             :   explicit FindContour(const ActionOptions&ao);
     107           0 :   bool checkAllActive() const { return gbuffer==0; }
     108             :   void prepareForAveraging();
     109           0 :   bool isPeriodic() { return false; }
     110          14 :   unsigned getNumberOfQuantities() const { return 1 + ingrid->getDimension(); }
     111             :   void compute( const unsigned& current, MultiValue& myvals ) const ;
     112             :   void finishAveraging();
     113             : };
     114             : 
     115        7358 : PLUMED_REGISTER_ACTION(FindContour,"FIND_CONTOUR")
     116             : 
     117           2 : void FindContour::registerKeywords( Keywords& keys ) {
     118           2 :   ContourFindingBase::registerKeywords( keys );
     119             : // We want a better way of doing this bit
     120          10 :   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");
     121           8 :   keys.add("compulsory","FILE","file on which to output coordinates");
     122          10 :   keys.add("compulsory","UNITS","PLUMED","the units in which to print out the coordinates. PLUMED means internal PLUMED units");
     123           8 :   keys.add("optional", "PRECISION","The number of digits in trajectory file");
     124           2 : }
     125             : 
     126           1 : FindContour::FindContour(const ActionOptions&ao):
     127             :   Action(ao),
     128             :   ContourFindingBase(ao),
     129           1 :   firsttime(true)
     130             : {
     131             : 
     132           2 :   parse("BUFFER",gbuffer);
     133           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);
     134             : 
     135           2 :   std::string file; parse("FILE",file);
     136           1 :   if( file.length()==0 ) error("name out output file was not specified");
     137           1 :   std::string type=Tools::extension(file);
     138           1 :   log<<"  file name "<<file<<"\n";
     139           1 :   if(type!="xyz") error("can only print xyz file type with contour finding");
     140             : 
     141             :   fmt_xyz="%f";
     142           2 :   std::string precision; parse("PRECISION",precision);
     143           1 :   if(precision.length()>0) {
     144           1 :     int p; Tools::convert(precision,p);
     145           1 :     log<<"  with precision "<<p<<"\n";
     146             :     std::string a,b;
     147           1 :     Tools::convert(p+5,a);
     148           1 :     Tools::convert(p,b);
     149           5 :     fmt_xyz="%"+a+"."+b+"f";
     150             :   }
     151           2 :   std::string unitname; parse("UNITS",unitname);
     152           1 :   if(unitname!="PLUMED") {
     153           0 :     Units myunit; myunit.setLength(unitname);
     154           0 :     lenunit=plumed.getAtoms().getUnits().getLength()/myunit.getLength();
     155             :   }
     156           1 :   else lenunit=1.0;
     157           1 :   of.link(*this); of.open(file);
     158           1 :   checkRead(); mydata=buildDataStashes( NULL );
     159           1 : }
     160             : 
     161           2 : void FindContour::prepareForAveraging() {
     162             :   // Create a task list if first time
     163           2 :   if( firsttime ) {
     164       16466 :     for(unsigned i=0; i<ingrid->getDimension()*ingrid->getNumberOfPoints(); ++i) addTaskToList( i );
     165             :   }
     166           2 :   firsttime=false; deactivateAllTasks();
     167             : 
     168             :   // We now need to identify the grid points that we need to search through
     169           2 :   std::vector<unsigned> nbin( ingrid->getNbin() );
     170           4 :   std::vector<unsigned> ind( ingrid->getDimension() );
     171           4 :   std::vector<unsigned> ones( ingrid->getDimension(), 1 );
     172             :   unsigned num_neighbours; std::vector<unsigned> neighbours;
     173       21956 :   for(unsigned i=0; i<ingrid->getNumberOfPoints(); ++i) {
     174             :     // Ensure inactive grid points are ignored
     175       10976 :     if( ingrid->inactive(i) ) continue;
     176             : 
     177             :     // Get the index of the current grid point
     178       10976 :     ingrid->getIndices( i, ind );
     179       10976 :     ingrid->getNeighbors( ind, ones, num_neighbours, neighbours );
     180             :     bool cycle=false;
     181      603680 :     for(unsigned j=0; j<num_neighbours; ++j) {
     182      592704 :       if( ingrid->inactive( neighbours[j]) ) { cycle=true; break; }
     183             :     }
     184       10976 :     if( cycle ) continue;
     185             : 
     186             :     // Get the value of a point on the grid
     187       10976 :     double val1=getFunctionValue( i ) - contour;
     188             :     bool edge=false;
     189      120736 :     for(unsigned j=0; j<ingrid->getDimension(); ++j) {
     190             :       // Make sure we don't search at the edge of the grid
     191       32928 :       if( !ingrid->isPeriodic(j) && (ind[j]+1)==nbin[j] ) continue;
     192       65856 :       else if( (ind[j]+1)==nbin[j] ) { edge=true; ind[j]=0; }
     193       30968 :       else ind[j]+=1;
     194       32928 :       double val2=getFunctionValue( ind ) - contour;
     195       35144 :       if( val1*val2<0 ) taskFlags[ ingrid->getDimension()*i + j ] = 1;
     196       69776 :       if( ingrid->isPeriodic(j) && edge ) { edge=false; ind[j]=nbin[j]-1; }
     197       30968 :       else ind[j]-=1;
     198             :     }
     199             :   }
     200           2 :   lockContributors();
     201           2 : }
     202             : 
     203         554 : void FindContour::compute( const unsigned& current, MultiValue& myvals ) const {
     204             :   // Retrieve the initial grid point coordinates
     205        1108 :   unsigned gpoint = std::floor( current / ingrid->getDimension() );
     206         554 :   std::vector<double> point( ingrid->getDimension() );
     207         554 :   ingrid->getGridPointCoordinates( gpoint, point );
     208             : 
     209             :   // Retrieve the direction we are searching for the contour
     210        1108 :   unsigned gdir = current%(ingrid->getDimension() );
     211         554 :   std::vector<double> direction( ingrid->getDimension(), 0 );
     212        1662 :   direction[gdir] = 0.999999999*ingrid->getGridSpacing()[gdir];
     213             : 
     214             :   // Now find the contour
     215             :   findContour( direction, point );
     216             :   // And transfer to the store data vessel
     217        7756 :   for(unsigned i=0; i<ingrid->getDimension(); ++i) myvals.setValue( 1+i, point[i] );
     218         554 : }
     219             : 
     220           1 : void FindContour::finishAveraging() {
     221             :   // And update the list of active grid points
     222           1 :   if( gbuffer>0 ) {
     223             :     std::vector<unsigned> neighbours; unsigned num_neighbours;
     224           0 :     std::vector<unsigned> ugrid_indices( ingrid->getDimension() );
     225           0 :     std::vector<bool> active( ingrid->getNumberOfPoints(), false );
     226           0 :     std::vector<unsigned> gbuffer_vec( ingrid->getDimension(), gbuffer );
     227           0 :     for(unsigned i=0; i<getCurrentNumberOfActiveTasks(); ++i) {
     228             :       // Get the point we are operating on
     229           0 :       unsigned ipoint = std::floor( getActiveTask(i) / ingrid->getDimension() );
     230             :       // Get the indices of this point
     231           0 :       ingrid->getIndices( ipoint, ugrid_indices );
     232             :       // Now activate buffer region
     233           0 :       ingrid->getNeighbors( ugrid_indices, gbuffer_vec, num_neighbours, neighbours );
     234           0 :       for(unsigned n=0; n<num_neighbours; ++n) active[ neighbours[n] ]=true;
     235             :     }
     236           0 :     ingrid->activateThesePoints( active );
     237             :   }
     238           2 :   std::vector<double> point( 1 + ingrid->getDimension() );
     239           1 :   of.printf("%u\n",mydata->getNumberOfStoredValues());
     240           1 :   of.printf("Points found on isocontour\n");
     241         555 :   for(unsigned i=0; i<mydata->getNumberOfStoredValues(); ++i) {
     242         554 :     mydata->retrieveSequentialValue( i, false, point ); of.printf("X");
     243        9418 :     for(unsigned j=0; j<ingrid->getDimension(); ++j) of.printf( (" " + fmt_xyz).c_str(), lenunit*point[1+j] );
     244         554 :     of.printf("\n");
     245             :   }
     246           1 : }
     247             : 
     248             : }
     249        5517 : }