Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2016-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 "ContourFindingObject.h"
24 : #include "gridtools/ActionWithGrid.h"
25 : #include "gridtools/EvaluateGridFunction.h"
26 : #include "core/ParallelTaskManager.h"
27 : #include "tools/Random.h"
28 :
29 : //+PLUMEDOC GRIDANALYSIS FIND_SPHERICAL_CONTOUR
30 : /*
31 : Find an isocontour in a three dimensional grid by searching over a Fibonacci sphere.
32 :
33 : As discussed in the documentation for the [gridtools](module_gridtools.md), PLUMED contains a number of tools that allow you to calculate
34 : a function on a grid. The function on this grid might be a [HISTOGRAM](HISTOGRAM.md) or it might be one of the phase fields that are
35 : discussed [here](module_contour.md). If this function has one or two input
36 : arguments it is relatively straightforward to plot the function. If by contrast the data has a three dimensions it can be
37 : difficult to visualize.
38 :
39 : This action provides one tool for visualizing these functions. It can be used to search for a set of points on a contour
40 : where the function takes a particular value. In other words, for the function $f(x,y,z)$ this action would find a set
41 : of points $\{x_c,y_c,z_c\}$ that have:
42 :
43 : $$
44 : f(x_c,y_c,z_c) - c = 0
45 : $$
46 :
47 : where $c$ is some constant value that is specified by the user. The points on this contour are find by searching along a
48 : set of equally spaced radii of a sphere that centered at on particular, user-specified atom or virtual atom. To ensure that
49 : these search radii are equally spaced on the surface of the sphere the search directions are generated by using a Fibonacci
50 : spiral projected on a sphere. In other words, the search directions are given by:
51 :
52 : $$
53 : \mathbf{r}_i = \left(
54 : \begin{matrix}
55 : \sqrt{1 - y^2} \cos(\phi) \\
56 : \frac{2i}{n} - 1 + \frac{1}{n} \\
57 : \sqrt{1 - y^2} \sin(\phi)
58 : \end{matrix}
59 : \right)
60 : $$
61 :
62 : where $y$ is the second component of the vector defined above, $n$ is the number of directions to look in and $\phi$ is
63 :
64 : $$
65 : \phi = \mod(i + R, n) \pi ( 3 - \sqrt{5} )
66 : $$
67 :
68 : where $R$ is a random variable between 0 and $n-1$ that is generated during the read in of the input file and that is fixed during
69 : the whole calculation.
70 :
71 : It is important to note that this action can only be used to detect contours in three dimensional functions. In addition, this action will fail to
72 : find the full set of contour points if the contour does not have the same topology as a sphere. If you are uncertain that the isocontours in your
73 : function have a spherical topology you should use [FIND_CONTOUR](FIND_CONTOUR.md) instead.
74 :
75 : ## Examples
76 :
77 : The following input demonstrates how this action can be used. The input here is used to study the shape of a droplet that has been formed during the
78 : condensation of Lennard Jones from the vapor. The input below achieves this by calculating the coordination numbers, $c_i$, of all the atoms within the gas.
79 : Obviously, those atoms within the droplet will have a large value for the coordination number while the isolated atoms in the gas will have a low value.
80 :
81 : We can detect the sizes of the droplets by constructing a matrix whose $i,j$ element tells us whether atom $i$ and atom $j$ are within 6 nm of each other and
82 : both have coordination numbers that are greater that two. The atoms within the various droplets within the system can then be found by performing a
83 : [DFSCLUSTERING](DFSCLUSTERING.md) on this matrix to detect the connected components. We can take the largest of these connected components and find the center of the droplet
84 : by exploiting the functionality within [CENTER](CENTER.md). We can then construct a phase field based on the positions of the atoms in the largest
85 : cluster and the values of the coordination numbers of these atoms as follows:
86 :
87 : $$
88 : \rho'(x,y,z) = \frac{ \sum_{i=1}^N f(s_i) K\left(\frac{x-x_i}{\lambda}, \frac{y-y_i}{\lambda}, \frac{z-z_i}{\lambda}\right) }{ \sum_{i=1}^N K\left(\frac{x-x_i}{\lambda}, \frac{y-y_i}{\lambda}, \frac{z-z_i}{\lambda}\right) }
89 : $$
90 :
91 : where $f$ is a switching function. The final line in the input then finds the a set of points on the dividing surface that separates
92 : the droplet from the surrounding gas. The value of the phase field on this isocontour is equal to 0.75.
93 :
94 : ```plumed
95 : # Calculate coordination numbers
96 : c1: COORDINATIONNUMBER SPECIES=1-512 SWITCH={EXP D_0=4.0 R_0=0.5 D_MAX=6.0}
97 : # Select coordination numbers that are more than 2.0
98 : cf: MORE_THAN ARG=c1 SWITCH={RATIONAL D_0=2.0 R_0=0.1}
99 : # Build a contact matrix
100 : c1_mat2: CONTACT_MATRIX GROUP=1-512 SWITCH={EXP D_0=4.0 R_0=0.5 D_MAX=6.0}
101 : dp_mat: OUTER_PRODUCT ARG=cf,cf
102 : # Build the final matrix
103 : mat: CUSTOM ARG=c1_mat2,dp_mat FUNC=x*y PERIODIC=NO
104 : # Find largest cluster
105 : dfs: DFSCLUSTERING ARG=mat
106 : clust1: CLUSTER_WEIGHTS CLUSTERS=dfs CLUSTER=1
107 : # Find center of largest cluster
108 : trans1: CUSTOM ARG=cf,clust1 FUNC=x*x*y PERIODIC=NO
109 : cent: CENTER ATOMS=1-512 WEIGHTS=trans1 PHASES
110 : # Calculate the phase field of the coordination
111 : dens_dist: DISTANCES ORIGIN=cent ATOMS=c1 COMPONENTS
112 : dens_numer: KDE ...
113 : VOLUMES=trans1 ARG=dens_dist.x,dens_dist.y,dens_dist.z
114 : GRID_BIN=30,30,30 BANDWIDTH=2.0,2.0,2.0
115 : ...
116 : dens_denom: KDE ...
117 : VOLUMES=clust1 ARG=dens_dist.x,dens_dist.y,dens_dist.z
118 : GRID_BIN=30,30,30 BANDWIDTH=2.0,2.0,2.0
119 : ...
120 : dens: CUSTOM ARG=dens_numer,dens_denom FUNC=x/y PERIODIC=NO
121 : # Find the isocontour around the nucleus
122 : sc: FIND_SPHERICAL_CONTOUR ARG=dens CONTOUR=0.85 INNER_RADIUS=10.0 OUTER_RADIUS=40.0 NPOINTS=100
123 : # And print the grid to a file
124 : DUMPGRID ARG=sc PRINT_XYZ FILE=mysurface.xyz STRIDE=1
125 : ```
126 :
127 : */
128 : //+ENDPLUMEDOC
129 :
130 : namespace PLMD {
131 : namespace contour {
132 :
133 1 : class FindSphericalContourObject {
134 : public:
135 : unsigned nbins;
136 : ContourFindingObject<gridtools::EvaluateGridFunction> cf;
137 3 : static void registerKeywords( Keywords& keys ) {
138 3 : ContourFindingObject<gridtools::EvaluateGridFunction>::registerKeywords( keys );
139 3 : keys.add("compulsory","NBINS","1","the number of discrete sections in which to divide the distance between the inner and outer radius when searching for a contour");
140 3 : }
141 1 : static void read( FindSphericalContourObject& func, ActionWithArguments* action, function::FunctionOptions& options ) {
142 1 : action->parse("NBINS",func.nbins);
143 1 : ContourFindingObject<gridtools::EvaluateGridFunction>::read( func.cf, action, options );
144 1 : }
145 : };
146 :
147 : class FindSphericalContour : public gridtools::ActionWithGrid {
148 : public:
149 : using input_type = FindSphericalContourObject;
150 : using PTM = ParallelTaskManager<FindSphericalContour>;
151 : private:
152 : /// The parallel task manager
153 : PTM taskmanager;
154 : unsigned npoints;
155 : double min, max;
156 : gridtools::GridCoordinatesObject gridcoords;
157 : public:
158 : static void registerKeywords( Keywords& keys );
159 : explicit FindSphericalContour(const ActionOptions&ao);
160 : unsigned getNumberOfDerivatives() override ;
161 : std::vector<std::string> getGridCoordinateNames() const override ;
162 : const gridtools::GridCoordinatesObject& getGridCoordinatesObject() const override ;
163 : void calculate() override ;
164 : void getInputData( std::vector<double>& inputdata ) const override;
165 : static void performTask( std::size_t task_index,
166 : const FindSphericalContourObject& actiondata,
167 : ParallelActionsInput& input,
168 : ParallelActionsOutput& output );
169 : };
170 :
171 : PLUMED_REGISTER_ACTION(FindSphericalContour,"FIND_SPHERICAL_CONTOUR")
172 :
173 3 : void FindSphericalContour::registerKeywords( Keywords& keys ) {
174 3 : ActionWithGrid::registerKeywords( keys );
175 6 : keys.addInputKeyword("compulsory","ARG","grid","the labels of the grid in which the contour will be found");
176 3 : keys.add("compulsory","NPOINTS","the number of points for which we are looking for the contour");
177 3 : keys.add("compulsory","INNER_RADIUS","the minimum radius on which to look for the contour");
178 3 : keys.add("compulsory","OUTER_RADIUS","the outer radius on which to look for the contour");
179 3 : FindSphericalContourObject::registerKeywords( keys );
180 6 : keys.setValueDescription("grid","a grid on a Fibonacci sphere that describes the radial distance from the origin for the points on the Willard-Chandler surface");
181 3 : keys.addDOI("10.1063/1.5134461");
182 3 : PTM::registerKeywords( keys );
183 3 : }
184 :
185 1 : FindSphericalContour::FindSphericalContour(const ActionOptions&ao):
186 : Action(ao),
187 : ActionWithGrid(ao),
188 1 : taskmanager(this) {
189 1 : if( getPntrToArgument(0)->getRank()!=3 ) {
190 0 : error("input grid must be three dimensional");
191 : }
192 :
193 1 : parse("NPOINTS",npoints);
194 1 : log.printf(" searching for %u points on dividing surface \n",npoints);
195 1 : parse("INNER_RADIUS",min);
196 1 : parse("OUTER_RADIUS",max);
197 1 : log.printf(" expecting to find dividing surface at radii between %f and %f \n",min,max);
198 : // Set this here so the same set of grid points are used on every turn
199 1 : std::vector<bool> ipbc( 3, false );
200 1 : gridcoords.setup( "fibonacci", ipbc, npoints, 0.0 );
201 : // Now create a value
202 1 : std::vector<std::size_t> shape( 3 );
203 1 : shape[0]=npoints;
204 1 : shape[1]=shape[2]=1;
205 1 : addValueWithDerivatives( shape );
206 1 : setNotPeriodic();
207 : function::FunctionOptions options;
208 1 : FindSphericalContourObject::read( taskmanager.getActionInput(), this, options );
209 1 : log.printf(" looking for contour in windows of length %f \n", (max-min)/taskmanager.getActionInput().nbins);
210 1 : log.printf(" calculating dividing surface along which function equals %f \n", taskmanager.getActionInput().cf.contour);
211 1 : taskmanager.setupParallelTaskManager( 0, 0 );
212 1 : }
213 :
214 5 : unsigned FindSphericalContour::getNumberOfDerivatives() {
215 5 : return gridcoords.getDimension();
216 : }
217 :
218 0 : std::vector<std::string> FindSphericalContour::getGridCoordinateNames() const {
219 0 : gridtools::ActionWithGrid* ag=dynamic_cast<gridtools::ActionWithGrid*>( getPntrToArgument(0)->getPntrToAction() );
220 0 : plumed_assert( ag );
221 0 : return ag->getGridCoordinateNames();
222 : }
223 :
224 3 : const gridtools::GridCoordinatesObject& FindSphericalContour::getGridCoordinatesObject() const {
225 3 : return gridcoords;
226 : }
227 :
228 2 : void FindSphericalContour::calculate() {
229 2 : taskmanager.runAllTasks();
230 2 : }
231 :
232 2 : void FindSphericalContour::getInputData( std::vector<double>& inputdata ) const {
233 2 : std::size_t ndata = gridcoords.getNumberOfPoints();
234 2 : if( inputdata.size()!=6*ndata ) {
235 1 : inputdata.resize( 6*ndata );
236 : }
237 2 : std::vector<double> direction(3);
238 202 : for(unsigned i=0; i<ndata; ++i) {
239 200 : gridcoords.getGridPointCoordinates( i, direction );
240 800 : for(unsigned j=0; j<3; ++j) {
241 600 : inputdata[6*i+j] = min*direction[j];
242 600 : inputdata[6*i+3+j] = (max-min)*direction[j] / static_cast<double>(taskmanager.getActionInput().nbins);
243 : }
244 : }
245 2 : }
246 :
247 200 : void FindSphericalContour::performTask( std::size_t task_index,
248 : const FindSphericalContourObject& actiondata,
249 : ParallelActionsInput& input,
250 : ParallelActionsOutput& output ) {
251 : bool found=false;
252 200 : View<const double> cp( input.inputdata + 6*task_index, 3 );
253 200 : std::vector<double> direction(3), contour_point(3), der(3), tmp(3);
254 200 : contour_point[0] = cp[0];
255 200 : contour_point[1] = cp[1];
256 200 : contour_point[2] = cp[2];
257 200 : View<const double> d( input.inputdata + 6*task_index+3, 3);
258 200 : direction[0] = d[0];
259 200 : direction[1] = d[1];
260 200 : direction[2] = d[2];
261 200 : for(unsigned k=0; k<actiondata.nbins; ++k) {
262 800 : for(unsigned j=0; j<3; ++j) {
263 600 : tmp[j] = contour_point[j] + direction[j];
264 : }
265 200 : double val1 = actiondata.cf.getDifferenceFromContour( contour_point, der );
266 200 : double val2 = actiondata.cf.getDifferenceFromContour( tmp, der );
267 200 : if( val1*val2<0 ) {
268 200 : ContourFindingObject<gridtools::EvaluateGridFunction>::findContour( actiondata.cf, direction, contour_point );
269 : double norm=0;
270 800 : for(unsigned j=0; j<3; ++j) {
271 600 : norm += contour_point[j]*contour_point[j];
272 : }
273 200 : output.values[0] = sqrt(norm);
274 : found=true;
275 : break;
276 : }
277 0 : for(unsigned j=0; j<3; ++j) {
278 0 : contour_point[j] = tmp[j];
279 : }
280 : }
281 : if( !found ) {
282 0 : plumed_merror("range does not bracket the dividing surface");
283 : }
284 200 : }
285 :
286 : }
287 : }
|
