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 "tools/Matrix.h"
26 : #include "core/ActionRegister.h"
27 :
28 : #include <string>
29 : #include <cmath>
30 :
31 : //+PLUMEDOC MATRIX TOPOLOGY_MATRIX
32 : /*
33 : Adjacency matrix in which two atoms are adjacent if they are connected topologically
34 :
35 : The functionality in this action was developed as part of a project that attempted to study
36 : the nucleation of bubbles. The idea was to develop a criterion that would allow one to determine
37 : if two gas atoms $i$ and $j$ are part of the same bubble or not. This criterion would then be used
38 : to construct a adjacency matrix, which could be used in the same way that [CONTACT_MATRIX](CONTACT_MATRIX.md) is used in other
39 : methods.
40 :
41 : The criterion that was developed to determine whether atom $i$ and $j$ are connected in this way works by
42 : determining if the density within a cylinder that is centered on the vector connecting atoms $i$ and $j$ is
43 : less than a certain threshold value. To make this determination we first determine whether any given atom $k$
44 : is within the cylinder centered on the vector connecting atoms $i$ and $j$ by using the following expression
45 :
46 : $$
47 : 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)
48 : $$
49 :
50 : 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$.
51 :
52 : 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:
53 :
54 : $$
55 : \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)
56 : $$
57 :
58 : 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$.
59 :
60 : The final value that is stored in the $i, j$ element of the output matrix is:
61 :
62 : $$
63 : T_{ij} = s_3(|\mathbf{r}_{ij}|)s_4(\rho_{ij})
64 : $$
65 :
66 : where $s_3$ and $s_4$ are switching functions.
67 :
68 : 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).
69 : However, if you would like to try this method an example input that uses this action would look like this:
70 :
71 : ```plumed
72 : mat: TOPOLOGY_MATRIX ...
73 : GROUP=1-85 BACKGROUND_ATOMS=86-210
74 : BIN_SIZE=1.02 SIGMA=0.17 KERNEL=triangular
75 : CYLINDER_SWITCH={RATIONAL R_0=0.5 D_MAX=1.0}
76 : SWITCH={RATIONAL D_0=30 R_0=0.5 D_MAX=32}
77 : RADIUS={RATIONAL D_0=0.375 R_0=0.1 D_MAX=0.43}
78 : DENSITY_THRESHOLD={RATIONAL R_0=0.1 D_MAX=0.5}
79 : ...
80 : ```
81 :
82 : 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.
83 : 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
84 : from the method as we are determining the connectivity between the atoms specified via the `GROUP` keyword.
85 :
86 : The above example assumes that you want to calculate the connectivity within a single group of atoms. If you would calculate connectivity between two different groups
87 : of atoms you use the GROUPA and GROUPB keywords as shown below:
88 :
89 : ```plumed
90 : mat: TOPOLOGY_MATRIX ...
91 : GROUPA=1-20 GROUPB=21-85 BACKGROUND_ATOMS=86-210
92 : BIN_SIZE=1.02 SIGMA=0.17 KERNEL=triangular
93 : CYLINDER_SWITCH={RATIONAL R_0=0.5 D_MAX=1.0}
94 : SWITCH={RATIONAL D_0=30 R_0=0.5 D_MAX=32}
95 : RADIUS={RATIONAL D_0=0.375 R_0=0.1 D_MAX=0.43}
96 : DENSITY_THRESHOLD={RATIONAL R_0=0.1 D_MAX=0.5}
97 : COMPONENTS NOPBC
98 : ...
99 : ```
100 :
101 : Notice that we have also added the NOPBC and COMPONENTS keywords in this input. The action above thus outputs four matrices with the labels
102 : `mat.w`, `mat.x`, `mat.y` and `mat.z.` The matrix with the label `mat.w` is the adjacency matrix
103 : that would be output if you had not added the COMPONENTS flag. The $i,j$ component of the matrices `mat.x`, `mat.y` and `mat.z` contain the $x$, $y$ and $z$
104 : components of the vector connecting atoms $i$ and $k$. Importantly, however, the components of these vectors are only stored in `mat.x`, `mat.y` and `mat.z`
105 : if the elements of `mat.w` are non-zero. Using the COMPONENTS flag in this way ensures that you can use HBOND_MATRIX in tandem with many of the functionalities
106 : that are part of the [symfunc module](module_symfunc.md).
107 :
108 : The NOPBC flag, meanwhile, ensures that all distances are calculated in a way that __does not__ take the periodic boundary conditions into account. By default,
109 : distances are calculated in a way that takes periodic boundary conditions into account.
110 :
111 : ## The MASK keyword
112 :
113 : You use the MASK keyword with TOPOLOGY_MATRIX in the same way that is used in [CONTACT_MATRIX](CONTACT_MATRIX.md). This keyword thus expects a vector in input,
114 : which tells TOPOLOGY_MATRIX that it is safe to not calculate certain rows of the output matrix. An example where this keyword is used is shown below:
115 :
116 : ```plumed
117 : # Fixed virtual atom which serves as the probe volume's center (pos. in nm)
118 : center: FIXEDATOM AT=2.5,2.5,2.5
119 : # Vector in which element i is one if atom i is in sphere of interest and zero otherwise
120 : sphere: INSPHERE ATOMS=1-85 CENTER=center RADIUS={GAUSSIAN D_0=0.5 R_0=0.01 D_MAX=0.52}
121 : # Calculates cooordination numbers
122 : cmap: TOPOLOGY_MATRIX ...
123 : GROUP=1-85 BACKGROUND_ATOMS=86-210
124 : BIN_SIZE=1.02 SIGMA=0.17 KERNEL=triangular
125 : CYLINDER_SWITCH={RATIONAL R_0=0.5 D_MAX=1.0}
126 : SWITCH={RATIONAL D_0=30 R_0=0.5 D_MAX=32}
127 : RADIUS={RATIONAL D_0=0.375 R_0=0.1 D_MAX=0.43}
128 : DENSITY_THRESHOLD={RATIONAL R_0=0.1 D_MAX=0.5}
129 : MASK=sphere
130 : ...
131 : ones: ONES SIZE=85
132 : cc: MATRIX_VECTOR_PRODUCT ARG=cmap,ones
133 : # Multiply coordination numbers by sphere vector
134 : prod: CUSTOM ARG=cc,sphere FUNC=x*y PERIODIC=NO
135 : # Sum of coordination numbers for atoms that are in the sphere of interest
136 : numer: SUM ARG=prod PERIODIC=NO
137 : # Number of atoms that are in sphere of interest
138 : denom: SUM ARG=sphere PERIODIC=NO
139 : # Average coordination number for atoms in sphere of interest
140 : av: CUSTOM ARG=prod,sphere FUNC=x/y PERIODIC=NO
141 : # And print out final CV to a file
142 : PRINT ARG=av FILE=colvar STRIDE=1
143 : ```
144 :
145 : This input calculates the average number of topological connections there are for each of the atoms that are within a spherical region
146 : that is centered on the point $(2.5,2.5,2.5)$.
147 :
148 :
149 : */
150 : //+ENDPLUMEDOC
151 :
152 : namespace PLMD {
153 : namespace adjmat {
154 :
155 : class TopologyMatrix {
156 : public:
157 : double sigma;
158 : HistogramBead::KernelType kerneltype;
159 : /// The maximum number of bins that will be used
160 : /// This is calculated based on the dmax of the switching functions
161 : unsigned maxbins;
162 : /// The volume of the cells
163 : double cell_volume;
164 : double binw_mat;
165 : /// switching functions
166 : SwitchingFunction switchingFunction;
167 : SwitchingFunction cylinder_sw;
168 : SwitchingFunction low_sf;
169 : SwitchingFunction threshold_switch;
170 : static void registerKeywords( Keywords& keys );
171 : void parseInput( AdjacencyMatrixBase<TopologyMatrix>* action );
172 : static void calculateWeight( const TopologyMatrix& data,
173 : const AdjacencyMatrixInput& input,
174 : MatrixOutput output );
175 : };
176 :
177 : typedef AdjacencyMatrixBase<TopologyMatrix> tmap;
178 : PLUMED_REGISTER_ACTION(tmap,"TOPOLOGY_MATRIX")
179 :
180 10 : void TopologyMatrix::registerKeywords( Keywords& keys ) {
181 10 : keys.add("atoms","BACKGROUND_ATOMS","the list of atoms that should be considered as part of the background density");
182 10 : keys.add("compulsory","SWITCH","This keyword is used if you want to employ an alternative to the continuous swiching function defined above. "
183 : "The following provides information on the \\ref switchingfunction that are available. "
184 : "When this keyword is present you no longer need the NN, MM, D_0 and R_0 keywords.");
185 20 : keys.linkActionInDocs("SWITCH","LESS_THAN");
186 10 : keys.add("compulsory","RADIUS","swtiching function that acts upon the radial distance of the atom from the center of the cylinder");
187 20 : keys.linkActionInDocs("RADIUS","LESS_THAN");
188 10 : keys.add("compulsory","CYLINDER_SWITCH","a switching function on ( r_ij . r_ik - 1 )/r_ij");
189 20 : keys.linkActionInDocs("CYLINDER_SWITCH","LESS_THAN");
190 10 : keys.add("compulsory","BIN_SIZE","the size to use for the bins");
191 10 : keys.add("compulsory","DENSITY_THRESHOLD","a switching function that acts upon the maximum density in the cylinder");
192 20 : keys.linkActionInDocs("DENSITY_THRESHOLD","LESS_THAN");
193 10 : keys.add("compulsory","SIGMA","the width of the function to be used for kernel density estimation");
194 10 : keys.add("compulsory","KERNEL","gaussian","the type of kernel function to be used");
195 10 : }
196 :
197 8 : void TopologyMatrix::parseInput( AdjacencyMatrixBase<TopologyMatrix>* action ) {
198 : std::string errors;
199 :
200 : std::string sfinput;
201 16 : action->parse("SWITCH",sfinput);
202 8 : if( sfinput.length()==0 ) {
203 0 : action->error("could not find SWITCH keyword");
204 : }
205 8 : switchingFunction.set(sfinput,errors);
206 8 : if( errors.length()!=0 ) {
207 0 : action->error("problem reading SWITCH keyword : " + errors );
208 : }
209 :
210 : std::string lowsfinput;
211 16 : action->parse("CYLINDER_SWITCH",lowsfinput);
212 8 : if( lowsfinput.length()==0 ) {
213 0 : action->error("could not find CYLINDER_SWITCH keyword");
214 : }
215 8 : low_sf.set(lowsfinput,errors);
216 8 : if( errors.length()!=0 ) {
217 0 : action->error("problem reading CYLINDER_SWITCH keyword : " + errors );
218 : }
219 :
220 : std::string cyinput;
221 16 : action->parse("RADIUS",cyinput);
222 8 : if( cyinput.length()==0 ) {
223 0 : action->error("could not find RADIUS keyword");
224 : }
225 8 : cylinder_sw.set(cyinput,errors);
226 8 : if( errors.length()!=0 ) {
227 0 : action->error("problem reading RADIUS keyword : " + errors );
228 : }
229 :
230 : std::string thresholdinput;
231 16 : action->parse("DENSITY_THRESHOLD",thresholdinput);
232 8 : if( thresholdinput.length()==0 ) {
233 0 : action->error("could not find DENSITY_THRESHOLD keyword");
234 : }
235 8 : threshold_switch.set(thresholdinput,errors);
236 8 : if( errors.length()!=0 ) {
237 0 : action->error("problem reading DENSITY_THRESHOLD keyword : " + errors );
238 : }
239 : // Read in stuff for grid
240 16 : action->parse("SIGMA",sigma);
241 : std::string mykerneltype;
242 8 : action->parse("KERNEL",mykerneltype);
243 8 : kerneltype = HistogramBead::getKernelType(mykerneltype);
244 8 : action->parse("BIN_SIZE",binw_mat);
245 :
246 : // Set the link cell cutoff
247 8 : action->setLinkCellCutoff( true, switchingFunction.get_dmax(),
248 : std::numeric_limits<double>::max() );
249 : // Set the number of bins
250 8 : maxbins = std::floor( switchingFunction.get_dmax() / binw_mat ) + 1;
251 : // Set the cell volume
252 8 : double r=cylinder_sw.get_d0() + cylinder_sw.get_r0();
253 8 : cell_volume=binw_mat*PLMD::pi*r*r;
254 8 : }
255 :
256 76368 : void TopologyMatrix::calculateWeight( const TopologyMatrix& data,
257 : const AdjacencyMatrixInput& input,
258 : MatrixOutput output ) {
259 : // Compute switching function on distance between atoms
260 76368 : Vector distance = input.pos;
261 : double len2 = distance.modulo2();
262 76368 : if( len2>data.switchingFunction.get_dmax2() ) {
263 71910 : return;
264 : }
265 : double dfuncl;
266 33550 : double sw = data.switchingFunction.calculateSqr( len2, dfuncl );
267 :
268 : // Now run through all sea atoms
269 33550 : HistogramBead bead( data.kerneltype, 0.0, data.binw_mat, data.sigma );
270 : Vector g1derivf,g2derivf,lderivf;
271 : Tensor vir;
272 33550 : double binlength = data.maxbins * data.binw_mat;
273 33550 : std::vector<double> tvals( data.maxbins, 0 );
274 33550 : Matrix<double> tvals_derivs( data.maxbins, 6 + 3*input.natoms + 9 );
275 : tvals_derivs = 0;
276 : // tvals.resize( data.maxbins, 6 + 3*input.natoms + 9, 0 );
277 237642622 : for(unsigned i=0; i<input.natoms; ++i) {
278 : // Position of sea atom (this will be the origin)
279 : Vector d2(input.extra_positions[i][0],
280 : input.extra_positions[i][1],
281 237609072 : input.extra_positions[i][2]);
282 : // Vector connecting sea atom and first in bond taking pbc into account
283 237609072 : Vector d20 = input.pbc->distance( d2, Vector(0,0,0) );
284 : // Vector connecting sea atom and second in bond taking pbc into account
285 237609072 : Vector d21 = input.pbc->distance( d2, input.pos );
286 : // Now length of bond modulus and so on -- no pbc here as we want sea atom in middle
287 : Vector d1 = delta( d20, d21 );
288 237609072 : double d1_len = d1.modulo();
289 237609072 : d1 = d1 / d1_len;
290 : // Switching function on distance between nodes
291 237609072 : if( d1_len>data.switchingFunction.get_dmax() ) {
292 187842012 : continue ;
293 : }
294 : // Ensure that the center of the bins are on the center of the bond connecting the two atoms
295 184481984 : double start2atom = 0.5*(binlength-d1_len);
296 : Vector dstart = d20 - start2atom*d1;
297 : // Now calculate projection of axis of cylinder
298 : double proj=dotProduct(-dstart,d1);
299 : // Calculate length of vector connecting start of cylinder to first atom
300 : // And now work out projection on vector connecting start and end of cylinder
301 184481984 : double proj_between = proj - start2atom;
302 : // This tells us if we are outside the end of the cylinder
303 184481984 : double excess = proj_between - d1_len;
304 : // Return if we are outside of the cylinder as calculated based on excess
305 184481984 : if( excess>data.low_sf.get_dmax() || -proj_between>data.low_sf.get_dmax() ) {
306 134714924 : continue;
307 : }
308 : // Calculate the excess swiching functions
309 : double edf1;
310 49767060 : double eval1 = data.low_sf.calculate( excess, edf1 );
311 : double edf2;
312 49767060 : double eval2 = data.low_sf.calculate( -proj_between, edf2 );
313 : // Calculate the projection on the perpendicular distance from the center of the tube
314 49767060 : double cm = dstart.modulo2() - proj*proj;
315 :
316 : // Now calculate the density in the cylinder
317 49767060 : if( cm>0 && cm<data.cylinder_sw.get_dmax2() ) {
318 : double dfuncr;
319 316634 : double val = data.cylinder_sw.calculateSqr( cm, dfuncr );
320 : Vector dc1, dc2, dc3, dd1, dd2, dd3, de1, de2, de3;
321 316634 : if( !input.noderiv ) {
322 : Tensor d1_a1;
323 : // Derivative of director connecting atom1 - atom2 wrt the position of atom 1
324 28284 : d1_a1(0,0) = ( -(d1[1]*d1[1]+d1[2]*d1[2])/d1_len ); // dx/dx
325 28284 : d1_a1(0,1) = ( d1[0]*d1[1]/d1_len ); // dx/dy
326 28284 : d1_a1(0,2) = ( d1[0]*d1[2]/d1_len ); // dx/dz
327 28284 : d1_a1(1,0) = ( d1[1]*d1[0]/d1_len ); // dy/dx
328 28284 : d1_a1(1,1) = ( -(d1[0]*d1[0]+d1[2]*d1[2])/d1_len ); // dy/dy
329 28284 : d1_a1(1,2) = ( d1[1]*d1[2]/d1_len );
330 28284 : d1_a1(2,0) = ( d1[2]*d1[0]/d1_len );
331 28284 : d1_a1(2,1) = ( d1[2]*d1[1]/d1_len );
332 28284 : d1_a1(2,2) = ( -(d1[1]*d1[1]+d1[0]*d1[0])/d1_len );
333 :
334 : // Calculate derivatives of dot product
335 28284 : dd1 = matmul(-dstart, d1_a1) - 0.5*d1;
336 28284 : dd2 = matmul(-dstart, -d1_a1) - 0.5*d1;
337 : dd3 = d1;
338 :
339 : // Calculate derivatives of cross product
340 28284 : Vector der( -0.5*binlength*matmul( d1_a1,dstart ) );
341 28284 : dc1 = dfuncr*( 0.5*dstart + der - proj*dd1 );
342 : dc2 = dfuncr*( 0.5*dstart - der - proj*dd2 );
343 : dc3 = dfuncr*( -dstart - proj*dd3 );
344 :
345 : // Calculate derivatives of excess
346 28284 : de1 = eval2*edf1*excess*(dd1 + 0.5*d1 )
347 28284 : + eval1*edf2*proj_between*(dd1 - 0.5*d1);
348 : de2 = eval2*edf1*excess*(dd2 - 0.5*d1 )
349 : + eval1*edf2*proj_between*(dd2 + 0.5*d1);
350 28284 : de3 = ( eval2*edf1*excess + eval1*edf2*proj_between )*dd3;
351 : }
352 3151486 : for(unsigned bin=0; bin<data.maxbins; ++bin) {
353 2834852 : bead.set( bin*data.binw_mat, (bin+1)*data.binw_mat, data.sigma );
354 2834852 : if( proj<(bin*data.binw_mat-bead.getCutoff())
355 2834852 : || proj>data.binw_mat*(bin+1)+bead.getCutoff() ) {
356 2409724 : continue;
357 : }
358 : double der;
359 425128 : double contr=bead.calculateWithCutoff( proj, der ) / data.cell_volume;
360 425128 : der /= data.cell_volume;
361 425128 : tvals[bin] += contr*val*eval1*eval2;
362 :
363 425128 : if( !input.noderiv ) {
364 38132 : g1derivf=contr*eval1*eval2*dc1 + val*eval1*eval2*der*dd1 + contr*val*de1;
365 38132 : tvals_derivs[bin][0] += g1derivf[0];
366 38132 : tvals_derivs[bin][1] += g1derivf[1];
367 38132 : tvals_derivs[bin][2] += g1derivf[2];
368 38132 : g2derivf=contr*eval1*eval2*dc2 + val*eval1*eval2*der*dd2 + contr*val*de2;
369 38132 : tvals_derivs[bin][3] += g2derivf[0];
370 38132 : tvals_derivs[bin][4] += g2derivf[1];
371 38132 : tvals_derivs[bin][5] += g2derivf[2];
372 38132 : lderivf=contr*eval1*eval2*dc3 + val*eval1*eval2*der*dd3 + contr*val*de3;
373 38132 : tvals_derivs[bin][6 + 3*i+0] += lderivf[0];
374 38132 : tvals_derivs[bin][6 + 3*i+1] += lderivf[1];
375 38132 : tvals_derivs[bin][6 + 3*i+2] += lderivf[2];
376 : // Virial
377 38132 : vir = - Tensor( d20, g1derivf ) - Tensor( d21, g2derivf );
378 38132 : unsigned nbase = 6 + 3*input.natoms;
379 38132 : tvals_derivs[bin][nbase+0] += vir(0,0);
380 38132 : tvals_derivs[bin][nbase+1] += vir(0,1);
381 38132 : tvals_derivs[bin][nbase+2] += vir(0,2);
382 38132 : tvals_derivs[bin][nbase+3] += vir(1,0);
383 38132 : tvals_derivs[bin][nbase+4] += vir(1,1);
384 38132 : tvals_derivs[bin][nbase+5] += vir(1,2);
385 38132 : tvals_derivs[bin][nbase+6] += vir(2,0);
386 38132 : tvals_derivs[bin][nbase+7] += vir(2,1);
387 38132 : tvals_derivs[bin][nbase+8] += vir(2,2);
388 : }
389 : }
390 : }
391 : }
392 : // Find maximum density
393 33550 : double max = tvals[0];
394 : unsigned vout = 0;
395 534514 : for(unsigned i=1; i<data.maxbins; ++i) {
396 500964 : if( tvals[i]>max ) {
397 : max=tvals[i];
398 : vout=i;
399 : }
400 : }
401 : // Transform the density
402 : double df;
403 33550 : double tsw = data.threshold_switch.calculate( max, df );
404 33550 : if( fabs(sw*tsw)<epsilon ) {
405 : return;
406 : }
407 :
408 4458 : if( !input.noderiv ) {
409 942 : Vector ddd = tsw*dfuncl*distance;
410 942 : output.deriv[0] = sw*df*max*tvals_derivs[vout][0] - ddd[0];
411 942 : output.deriv[1] = sw*df*max*tvals_derivs[vout][1] - ddd[1];
412 942 : output.deriv[2] = sw*df*max*tvals_derivs[vout][2] - ddd[2];
413 942 : output.deriv[3] = sw*df*max*tvals_derivs[vout][3] + ddd[0];
414 942 : output.deriv[4] = sw*df*max*tvals_derivs[vout][4] + ddd[1];
415 942 : output.deriv[5] = sw*df*max*tvals_derivs[vout][5] + ddd[2];
416 109488 : for(unsigned i=0; i<input.natoms; ++i) {
417 108546 : output.deriv[6 + 3*i + 0] = sw*df*max*tvals_derivs[vout][6 + 3*i + 0];
418 108546 : output.deriv[6 + 3*i + 1] = sw*df*max*tvals_derivs[vout][6 + 3*i + 1];
419 108546 : output.deriv[6 + 3*i + 2] = sw*df*max*tvals_derivs[vout][6 + 3*i + 2];
420 : }
421 942 : unsigned nbase = 6 + 3*input.natoms;
422 942 : Tensor vird(ddd,distance);
423 942 : output.deriv[nbase + 0] = sw*df*max*tvals_derivs[vout][nbase+0] - vird(0,0);
424 942 : output.deriv[nbase + 1] = sw*df*max*tvals_derivs[vout][nbase+1] - vird(0,1);
425 942 : output.deriv[nbase + 2] = sw*df*max*tvals_derivs[vout][nbase+2] - vird(0,2);
426 942 : output.deriv[nbase + 3] = sw*df*max*tvals_derivs[vout][nbase+3] - vird(1,0);
427 942 : output.deriv[nbase + 4] = sw*df*max*tvals_derivs[vout][nbase+4] - vird(1,1);
428 942 : output.deriv[nbase + 5] = sw*df*max*tvals_derivs[vout][nbase+5] - vird(1,2);
429 942 : output.deriv[nbase + 6] = sw*df*max*tvals_derivs[vout][nbase+6] - vird(2,0);
430 942 : output.deriv[nbase + 7] = sw*df*max*tvals_derivs[vout][nbase+7] - vird(2,1);
431 942 : output.deriv[nbase + 8] = sw*df*max*tvals_derivs[vout][nbase+8] - vird(2,2);
432 : }
433 4458 : output.val[0] = sw*tsw;
434 : }
435 :
436 : }
437 : }
|
