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 "AdjacencyMatrixBase.h"
23 : #include "tools/SwitchingFunction.h"
24 : #include "tools/HistogramBead.h"
25 : #include "core/ActionRegister.h"
26 :
27 : #include <string>
28 : #include <cmath>
29 :
30 : //+PLUMEDOC MATRIX TOPOLOGY_MATRIX
31 : /*
32 : Adjacency matrix in which two atoms are adjacent if they are connected topologically
33 :
34 : The functionality in this action was developed as part of a project that attempted to study
35 : the nucleation of bubbles. The idea was to develop a criterion that would allow one to determine
36 : if two gas atoms $i$ and $j$ are part of the same bubble or not. This criterion would then be used
37 : to construct a adjacency matrix, which could be used in the same way that [CONTACT_MATRIX](CONTACT_MATRIX.md) is used in other
38 : methods.
39 :
40 : The criterion that was developed to determine whether atom $i$ and $j$ are connected in this way works by
41 : determining if the density within a cylinder that is centered on the vector connecting atoms $i$ and $j$ is
42 : less than a certain threshold value. To make this determination we first determine whether any given atom $k$
43 : is within the cylinder centered on the vector connecting atoms $i$ and $j$ by using the following expression
44 :
45 : $$
46 : f(\mathbf{r}_{ik}, \mathbf{r}_{ij}) = s_1\left( \sqrt{ |\mathbf{r}_{ij}|^2 - \left( \frac{\mathbf{r}_{ij} \cdot \mathbf{r}_{ik}}{|\mathbf{r}_{ij} |} \right)^2} \right)s_2\left( -\frac{\mathbf{r}_{ij} \cdot \mathbf{r}_{ik}}{|\mathbf{r}_{ij} |} \right) s_2\left( \frac{\mathbf{r}_{ij} \cdot \mathbf{r}_{ik}}{|\mathbf{r}_{ij} |} - |\mathbf{r}_{ij}| \right)
47 : $$
48 :
49 : In this expression $s_1$ and $s_2$ are switching functions, while $\mathbf{r}_{lm}$ is used to indicate the vector connecting atoms $l$ and $m$.
50 :
51 : We then calculate the density for a grid of $M$ points along the vector connecting atom $i$ and atom $j$ using and find the maximum density on this grid using:
52 :
53 : $$
54 : \rho_{ij} = \max_{m \in M} \left[ \frac{M}{d_\textrm{max}} \right] \sum_k f(r_{ik}, r_{ij}) \int_{(m-1)d_{\textrm{max}}/M}^{ md_{\textrm{max}} /M } \textrm{d}x \quad K\left( \frac{x - r_{ks} \cdot r_{ij} }{ | r_{ks} | }\right)
55 : $$
56 :
57 : where $d_\textrm{max}$ is the `D_MAX` parameter of the switching function $s_3$ that appears in the next equation, $K$ is a kernal function and $s$ is used to represent a point in space that is $d_\textrm{max}$ from atom $j$ along the vector connecting atom $j$ to atom $i$.
58 :
59 : The final value that is stored in the $i, j$ element of the output matrix is:
60 :
61 : $$
62 : T_{ij} = s_3(|\mathbf{r}_{ij}|)s_4(\rho_{ij})
63 : $$
64 :
65 : where $s_3$ and $s_4$ are switching functions.
66 :
67 : We ended up abandoning this method and the project (if you want drive bubble formation you are probably better off adding a bias on the volume of the simulation cell).
68 : However, if you would like to try this method an example input that uses this action would look like this:
69 :
70 : ```plumed
71 : mat: TOPOLOGY_MATRIX ...
72 : GROUP=1-85 BACKGROUND_ATOMS=86-210
73 : BIN_SIZE=1.02 SIGMA=0.17 KERNEL=triangular
74 : CYLINDER_SWITCH={RATIONAL R_0=0.5 D_MAX=1.0}
75 : SWITCH={RATIONAL D_0=30 R_0=0.5 D_MAX=32}
76 : RADIUS={RATIONAL D_0=0.375 R_0=0.1 D_MAX=0.43}
77 : DENSITY_THRESHOLD={RATIONAL R_0=0.1 D_MAX=0.5}
78 : ...
79 : ```
80 :
81 : The switching functions $s_1$, $s_2$, $s_3$ and $s_4$ are specified using the `RADIUS`, `CYLINDER_SWITCH`, `SWITCH` and `DENSITY_THRESHOLD` keywords respectively.
82 : We loop over the atoms in the group specified using the `BACKGROUND_ATOMS` keyword when looping over $k$ in the formulas above. An $85 \times 85$ matrix is output
83 : from the method as we are determining the connectivity between the atoms specified via the `GROUP` keyword.
84 :
85 : */
86 : //+ENDPLUMEDOC
87 :
88 : namespace PLMD {
89 : namespace adjmat {
90 :
91 : class TopologyMatrix : public AdjacencyMatrixBase {
92 : private:
93 : /// The width to use for the kernel density estimation and the
94 : /// sizes of the bins to be used in kernel density estimation
95 : double sigma;
96 : std::string kerneltype;
97 : /// The maximum number of bins that will be used
98 : /// This is calculated based on the dmax of the switching functions
99 : unsigned maxbins;
100 : /// The volume of the cells
101 : double cell_volume;
102 : /// switching function
103 : SwitchingFunction switchingFunction;
104 : SwitchingFunction cylinder_sw;
105 : SwitchingFunction low_sf;
106 : double binw_mat;
107 : SwitchingFunction threshold_switch;
108 : public:
109 : static void registerKeywords( Keywords& keys );
110 : explicit TopologyMatrix(const ActionOptions&);
111 : // active methods:
112 : double calculateWeight( const Vector& pos1, const Vector& pos2, const unsigned& natoms, MultiValue& myvals ) const override;
113 : };
114 :
115 : PLUMED_REGISTER_ACTION(TopologyMatrix,"TOPOLOGY_MATRIX")
116 :
117 10 : void TopologyMatrix::registerKeywords( Keywords& keys ) {
118 10 : AdjacencyMatrixBase::registerKeywords( keys );
119 10 : keys.add("atoms","BACKGROUND_ATOMS","the list of atoms that should be considered as part of the background density");
120 10 : keys.add("compulsory","SWITCH","This keyword is used if you want to employ an alternative to the continuous swiching function defined above. "
121 : "The following provides information on the \\ref switchingfunction that are available. "
122 : "When this keyword is present you no longer need the NN, MM, D_0 and R_0 keywords.");
123 20 : keys.linkActionInDocs("SWITCH","LESS_THAN");
124 10 : keys.add("compulsory","RADIUS","swtiching function that acts upon the radial distance of the atom from the center of the cylinder");
125 20 : keys.linkActionInDocs("RADIUS","LESS_THAN");
126 10 : keys.add("compulsory","CYLINDER_SWITCH","a switching function on ( r_ij . r_ik - 1 )/r_ij");
127 20 : keys.linkActionInDocs("CYLINDER_SWITCH","LESS_THAN");
128 10 : keys.add("compulsory","BIN_SIZE","the size to use for the bins");
129 10 : keys.add("compulsory","DENSITY_THRESHOLD","a switching function that acts upon the maximum density in the cylinder");
130 20 : keys.linkActionInDocs("DENSITY_THRESHOLD","LESS_THAN");
131 10 : keys.add("compulsory","SIGMA","the width of the function to be used for kernel density estimation");
132 10 : keys.add("compulsory","KERNEL","gaussian","the type of kernel function to be used");
133 10 : }
134 :
135 8 : TopologyMatrix::TopologyMatrix(const ActionOptions&ao):
136 : Action(ao),
137 8 : AdjacencyMatrixBase(ao) {
138 : std::string sfinput,errors;
139 16 : parse("SWITCH",sfinput);
140 8 : if( sfinput.length()==0 ) {
141 0 : error("could not find SWITCH keyword");
142 : }
143 8 : switchingFunction.set(sfinput,errors);
144 8 : if( errors.length()!=0 ) {
145 0 : error("problem reading SWITCH keyword : " + errors );
146 : }
147 :
148 : std::string hsfinput;
149 16 : parse("CYLINDER_SWITCH",hsfinput);
150 8 : if( hsfinput.length()==0 ) {
151 0 : error("could not find CYLINDER_SWITCH keyword");
152 : }
153 8 : low_sf.set(hsfinput,errors);
154 8 : if( errors.length()!=0 ) {
155 0 : error("problem reading CYLINDER_SWITCH keyword : " + errors );
156 : }
157 :
158 : std::string asfinput;
159 16 : parse("RADIUS",asfinput);
160 8 : if( asfinput.length()==0 ) {
161 0 : error("could not find RADIUS keyword");
162 : }
163 8 : cylinder_sw.set(asfinput,errors);
164 8 : if( errors.length()!=0 ) {
165 0 : error("problem reading RADIUS keyword : " + errors );
166 : }
167 :
168 : std::string tsfinput;
169 16 : parse("DENSITY_THRESHOLD",tsfinput);
170 8 : if( tsfinput.length()==0 ) {
171 0 : error("could not find DENSITY_THRESHOLD keyword");
172 : }
173 8 : threshold_switch.set(tsfinput,errors);
174 8 : if( errors.length()!=0 ) {
175 0 : error("problem reading DENSITY_THRESHOLD keyword : " + errors );
176 : }
177 : // Read in stuff for grid
178 8 : parse("SIGMA",sigma);
179 8 : parse("KERNEL",kerneltype);
180 8 : parse("BIN_SIZE",binw_mat);
181 :
182 : // Set the link cell cutoff
183 8 : setLinkCellCutoff( true, switchingFunction.get_dmax(), std::numeric_limits<double>::max() );
184 : // Set the number of bins
185 8 : maxbins = std::floor( switchingFunction.get_dmax() / binw_mat ) + 1;
186 : // Set the cell volume
187 8 : double r=cylinder_sw.get_d0() + cylinder_sw.get_r0();
188 8 : cell_volume=binw_mat*pi*r*r;
189 :
190 : // And check everything has been read in correctly
191 8 : checkRead();
192 8 : }
193 :
194 69150 : double TopologyMatrix::calculateWeight( const Vector& pos1, const Vector& pos2, const unsigned& natoms, MultiValue& myvals ) const {
195 : // Compute switching function on distance between atoms
196 : Vector distance = pbcDistance( pos1, pos2 );
197 69150 : double len2 = distance.modulo2();
198 69150 : if( len2>switchingFunction.get_dmax2() ) {
199 : return 0.0;
200 : }
201 26332 : double dfuncl, sw = switchingFunction.calculateSqr( len2, dfuncl );
202 :
203 : // Now run through all sea atoms
204 26332 : HistogramBead bead;
205 : bead.isNotPeriodic();
206 26332 : bead.setKernelType( kerneltype );
207 26332 : Vector g1derivf,g2derivf,lderivf;
208 26332 : Tensor vir;
209 26332 : double binlength = maxbins * binw_mat;
210 26332 : MultiValue tvals( maxbins, myvals.getNumberOfDerivatives() );
211 236742358 : for(unsigned i=0; i<natoms; ++i) {
212 : // Position of sea atom (this will be the origin)
213 : Vector d2 = getPosition(i,myvals);
214 : // Vector connecting sea atom and first in bond taking pbc into account
215 : Vector d20 = pbcDistance( d2, pos1 );
216 : // Vector connecting sea atom and second in bond taking pbc into account
217 : Vector d21 = pbcDistance( d2, pos2 );
218 : // Now length of bond modulus and so on -- no pbc here as we want sea atom in middle
219 236716026 : Vector d1 = delta( d20, d21 );
220 236716026 : double d1_len = d1.modulo();
221 236716026 : d1 = d1 / d1_len;
222 : // Switching function on distance between nodes
223 236716026 : if( d1_len>switchingFunction.get_dmax() ) {
224 187754296 : continue ;
225 : }
226 : // Ensure that the center of the bins are on the center of the bond connecting the two atoms
227 183588938 : double start2atom = 0.5*(binlength-d1_len);
228 183588938 : Vector dstart = d20 - start2atom*d1;
229 : // Now calculate projection of axis of cylinder
230 183588938 : double proj=dotProduct(-dstart,d1);
231 : // Calculate length of vector connecting start of cylinder to first atom
232 : // And now work out projection on vector connecting start and end of cylinder
233 183588938 : double proj_between = proj - start2atom;
234 : // This tells us if we are outside the end of the cylinder
235 183588938 : double excess = proj_between - d1_len;
236 : // Return if we are outside of the cylinder as calculated based on excess
237 183588938 : if( excess>low_sf.get_dmax() || -proj_between>low_sf.get_dmax() ) {
238 134627208 : continue;
239 : }
240 : // Calculate the excess swiching functions
241 48961730 : double edf1, eval1 = low_sf.calculate( excess, edf1 );
242 48961730 : double edf2, eval2 = low_sf.calculate( -proj_between, edf2 );
243 : // Calculate the projection on the perpendicular distance from the center of the tube
244 48961730 : double cm = dstart.modulo2() - proj*proj;
245 :
246 : // Now calculate the density in the cylinder
247 48961730 : if( cm>0 && cm<cylinder_sw.get_dmax2() ) {
248 288350 : double dfuncr, val = cylinder_sw.calculateSqr( cm, dfuncr );
249 288350 : Vector dc1, dc2, dc3, dd1, dd2, dd3, de1, de2, de3;
250 288350 : if( !doNotCalculateDerivatives() ) {
251 28284 : Tensor d1_a1;
252 : // Derivative of director connecting atom1 - atom2 wrt the position of atom 1
253 28284 : d1_a1(0,0) = ( -(d1[1]*d1[1]+d1[2]*d1[2])/d1_len ); // dx/dx
254 28284 : d1_a1(0,1) = ( d1[0]*d1[1]/d1_len ); // dx/dy
255 28284 : d1_a1(0,2) = ( d1[0]*d1[2]/d1_len ); // dx/dz
256 28284 : d1_a1(1,0) = ( d1[1]*d1[0]/d1_len ); // dy/dx
257 28284 : d1_a1(1,1) = ( -(d1[0]*d1[0]+d1[2]*d1[2])/d1_len ); // dy/dy
258 28284 : d1_a1(1,2) = ( d1[1]*d1[2]/d1_len );
259 28284 : d1_a1(2,0) = ( d1[2]*d1[0]/d1_len );
260 28284 : d1_a1(2,1) = ( d1[2]*d1[1]/d1_len );
261 28284 : d1_a1(2,2) = ( -(d1[1]*d1[1]+d1[0]*d1[0])/d1_len );
262 :
263 : // Calculate derivatives of dot product
264 28284 : dd1 = matmul(-dstart, d1_a1) - 0.5*d1;
265 28284 : dd2 = matmul(-dstart, -d1_a1) - 0.5*d1;
266 28284 : dd3 = d1;
267 :
268 : // Calculate derivatives of cross product
269 28284 : Vector der( -0.5*binlength*matmul( d1_a1,dstart ) );
270 28284 : dc1 = dfuncr*( 0.5*dstart + der - proj*dd1 );
271 28284 : dc2 = dfuncr*( 0.5*dstart - der - proj*dd2 );
272 28284 : dc3 = dfuncr*( -dstart - proj*dd3 );
273 :
274 : // Calculate derivatives of excess
275 28284 : de1 = eval2*edf1*excess*(dd1 + 0.5*d1 ) + eval1*edf2*proj_between*(dd1 - 0.5*d1);
276 28284 : de2 = eval2*edf1*excess*(dd2 - 0.5*d1 ) + eval1*edf2*proj_between*(dd2 + 0.5*d1);
277 28284 : de3 = ( eval2*edf1*excess + eval1*edf2*proj_between )*dd3;
278 : }
279 2225914 : for(unsigned bin=0; bin<maxbins; ++bin) {
280 1937564 : bead.set( bin*binw_mat, (bin+1)*binw_mat, sigma );
281 1937564 : if( proj<(bin*binw_mat-bead.getCutoff()) || proj>binw_mat*(bin+1)+bead.getCutoff() ) {
282 1550568 : continue;
283 : }
284 386996 : double der, contr=bead.calculateWithCutoff( proj, der ) / cell_volume;
285 386996 : der /= cell_volume;
286 386996 : tvals.addValue( bin, contr*val*eval1*eval2 );
287 :
288 386996 : if( !doNotCalculateDerivatives() ) {
289 38132 : g1derivf=contr*eval1*eval2*dc1 + val*eval1*eval2*der*dd1 + contr*val*de1;
290 38132 : tvals.addDerivative( bin, 3*myvals.getTaskIndex()+0, g1derivf[0] );
291 38132 : tvals.addDerivative( bin, 3*myvals.getTaskIndex()+1, g1derivf[1] );
292 38132 : tvals.addDerivative( bin, 3*myvals.getTaskIndex()+2, g1derivf[2] );
293 38132 : g2derivf=contr*eval1*eval2*dc2 + val*eval1*eval2*der*dd2 + contr*val*de2;
294 38132 : tvals.addDerivative( bin, 3*myvals.getSecondTaskIndex()+0, g2derivf[0] );
295 38132 : tvals.addDerivative( bin, 3*myvals.getSecondTaskIndex()+1, g2derivf[1] );
296 38132 : tvals.addDerivative( bin, 3*myvals.getSecondTaskIndex()+2, g2derivf[2] );
297 38132 : lderivf=contr*eval1*eval2*dc3 + val*eval1*eval2*der*dd3 + contr*val*de3;
298 38132 : unsigned tindex = myvals.getIndices()[ i + myvals.getSplitIndex() ];
299 38132 : tvals.addDerivative( bin, 3*tindex+0, lderivf[0] );
300 38132 : tvals.addDerivative( bin, 3*tindex+1, lderivf[1] );
301 38132 : tvals.addDerivative( bin, 3*tindex+2, lderivf[2] );
302 : // Virial
303 38132 : vir = - Tensor( d20, g1derivf ) - Tensor( d21, g2derivf );
304 38132 : unsigned nbase = 3*getNumberOfAtoms();
305 38132 : tvals.addDerivative( bin, nbase+0, vir(0,0) );
306 38132 : tvals.addDerivative( bin, nbase+1, vir(0,1) );
307 38132 : tvals.addDerivative( bin, nbase+2, vir(0,2) );
308 38132 : tvals.addDerivative( bin, nbase+3, vir(1,0) );
309 38132 : tvals.addDerivative( bin, nbase+4, vir(1,1) );
310 38132 : tvals.addDerivative( bin, nbase+5, vir(1,2) );
311 38132 : tvals.addDerivative( bin, nbase+6, vir(2,0) );
312 38132 : tvals.addDerivative( bin, nbase+7, vir(2,1) );
313 38132 : tvals.addDerivative( bin, nbase+8, vir(2,2) );
314 : }
315 : }
316 : }
317 : }
318 : // Find maximum density
319 : double max = tvals.get(0);
320 : unsigned vout = 0;
321 305488 : for(unsigned i=1; i<maxbins; ++i) {
322 279156 : if( tvals.get(i)>max ) {
323 : max=tvals.get(i);
324 : vout=i;
325 : }
326 : }
327 : // Transform the density
328 26332 : double df, tsw = threshold_switch.calculate( max, df );
329 26332 : if( fabs(sw*tsw)<epsilon ) {
330 : return 0;
331 : }
332 :
333 3516 : if( !doNotCalculateDerivatives() ) {
334 942 : Vector ader;
335 942 : Tensor vir;
336 942 : Vector ddd = tsw*dfuncl*distance;
337 942 : ader[0] = tvals.getDerivative( vout, 3*myvals.getTaskIndex()+0 );
338 942 : ader[1] = tvals.getDerivative( vout, 3*myvals.getTaskIndex()+1 );
339 942 : ader[2] = tvals.getDerivative( vout, 3*myvals.getTaskIndex()+2 );
340 942 : addAtomDerivatives( 0, sw*df*max*ader - ddd, myvals );
341 942 : ader[0] = tvals.getDerivative( vout, 3*myvals.getSecondTaskIndex()+0 );
342 942 : ader[1] = tvals.getDerivative( vout, 3*myvals.getSecondTaskIndex()+1 );
343 942 : ader[2] = tvals.getDerivative( vout, 3*myvals.getSecondTaskIndex()+2 );
344 942 : addAtomDerivatives( 1, sw*df*max*ader + ddd, myvals );
345 109488 : for(unsigned i=0; i<natoms; ++i) {
346 108546 : unsigned tindex = myvals.getIndices()[ i + myvals.getSplitIndex() ];
347 108546 : ader[0] = tvals.getDerivative( vout, 3*tindex+0 );
348 108546 : ader[1] = tvals.getDerivative( vout, 3*tindex+1 );
349 108546 : ader[2] = tvals.getDerivative( vout, 3*tindex+2 );
350 108546 : addThirdAtomDerivatives( i, sw*df*max*ader, myvals );
351 : }
352 942 : unsigned nbase = 3*getNumberOfAtoms();
353 942 : Tensor vird(ddd,distance);
354 942 : vir(0,0) = sw*df*max*tvals.getDerivative( vout, nbase+0 ) - vird(0,0);
355 942 : vir(0,1) = sw*df*max*tvals.getDerivative( vout, nbase+1 ) - vird(0,1);
356 942 : vir(0,2) = sw*df*max*tvals.getDerivative( vout, nbase+2 ) - vird(0,2);
357 942 : vir(1,0) = sw*df*max*tvals.getDerivative( vout, nbase+3 ) - vird(1,0);
358 942 : vir(1,1) = sw*df*max*tvals.getDerivative( vout, nbase+4 ) - vird(1,1);
359 942 : vir(1,2) = sw*df*max*tvals.getDerivative( vout, nbase+5 ) - vird(1,2);
360 942 : vir(2,0) = sw*df*max*tvals.getDerivative( vout, nbase+6 ) - vird(2,0);
361 942 : vir(2,1) = sw*df*max*tvals.getDerivative( vout, nbase+7 ) - vird(2,1);
362 942 : vir(2,2) = sw*df*max*tvals.getDerivative( vout, nbase+8 ) - vird(2,2);
363 942 : addBoxDerivatives( vir, myvals );
364 : }
365 : return sw*tsw;
366 26332 : }
367 :
368 : }
369 : }
|
