Line data Source code
1 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2 : Copyright (c) 2012-2017 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 : #ifdef __PLUMED_HAS_OPENACC
23 : #define __PLUMED_USE_OPENACC 1
24 : #endif //__PLUMED_HAS_OPENACC
25 : #include "function/FunctionSetup.h"
26 : #include "function/FunctionShortcut.h"
27 : #include "function/FunctionOfMatrix.h"
28 : #include "core/ActionRegister.h"
29 :
30 : #include <complex>
31 :
32 : namespace PLMD {
33 : namespace symfunc {
34 :
35 : //+PLUMEDOC MCOLVAR SPHERICAL_HARMONIC
36 : /*
37 : Calculate the values of all the spherical harmonic funtions for a particular value of l.
38 :
39 : This action allows you to the all the [spherical harmonic](https://en.wikipedia.org/wiki/Spherical_harmonics) functions for a particular input
40 : $L$ value. As discussed in more detail in the article provided above the spherical harmonics this action calculates have the following general form:
41 :
42 : $$
43 : Y_l^m(\theta,\phi) = wNe^{im\phi} P_l^m(\cos\theta)
44 : $$
45 :
46 : where $N$ is a normalisation constant, $P_l^m$ is tn associated Legendre Polynomial and $w$ is an optional weight that can be passed in an input argument or simply set equal to one.
47 : $e^{i\phi}$ and $\cos(\theta) are computed from the other input arguments, $x, y$ and $z$ as follows:
48 :
49 : $$
50 : e^{i\phi} = \frac{x}{r} + i\frac{y}{r} \qquad \textrm{and} \qquad \cos(\theta) = \frac{z}{r} \qquad \textrm{where} \qquad r = \sqrt{x^2 + y^2 + z^2}
51 : $$
52 :
53 : At present, the arguments for this action must be matrices. However, it would be easy to add functionality that would allow you to compute this function for scalar or vector input.
54 : The arguments must all have the same shape as the two output components will also be matrices that are
55 : calculated by applying the function above to each of the elements of the input matrix in turn. The number of components output will be equal to $2(2L+1)$ and will contain
56 : the real and imaginary parts of the $Y_l^m$ functions with the the $2l+1$ possible $m$ values.
57 :
58 : The following intput provides an example that demonstrates how this function is used:
59 :
60 : ```plumed
61 : d: DISTANCE_MATRIX GROUP=1-10 COMPONENTS
62 : c: SPHERICAL_HARMONIC L=1 ARG=d.x,d.y,d.z
63 : PRINT ARG=c.rm-n1 FILE=real_part_m-1
64 : PRINT ARG=c.im-n1 FILE=imaginary_part_m-1
65 : PRINT ARG=c.rm-0 FILE=real_part_m0
66 : PRINT ARG=c.im-0 FILE=imaginary_part_m0
67 : PRINT ARG=c.rm-p1 FILE=real_part_m+1
68 : PRINT ARG=c.im-p1 FILE=imaginary_part_m+1
69 : ```
70 :
71 : The DISTANCE_MATRIX command in the above input computes 3 $10\times10$ matrices. These 3 $10\times10$ matrices are used in the input to the sphierical harmonic command,
72 : which in turn outputs 6 $10\times10$ matrices that contain the real and imaginary parts when the three spherical harmonic functions with $l=1$ are applied element-wise to the above input. These six $10\times10$
73 : matrices are then output to six separate files.
74 :
75 : In the above example the weights for every distance is set equal to one. The following example shows how an argument can be used to set the $w$ values to use when computing the function
76 : above.
77 :
78 : ```plumed
79 : s: CONTACT_MATRIX GROUP=1-10 SWITCH={RATIONAL R_0=1.0}
80 : sc: CONTACT_MATRIX GROUP=1-10 SWITCH={RATIONAL R_0=1.0} COMPONENTS
81 : c: SPHERICAL_HARMONIC L=1 ARG=sc.x,sc.y,sc.z,s
82 : PRINT ARG=c.rm-n1 FILE=real_part_m-1
83 : PRINT ARG=c.im-n1 FILE=imaginary_part_m-1
84 : PRINT ARG=c.rm-0 FILE=real_part_m0
85 : PRINT ARG=c.im-0 FILE=imaginary_part_m0
86 : PRINT ARG=c.rm-p1 FILE=real_part_m+1
87 : PRINT ARG=c.im-p1 FILE=imaginary_part_m+1
88 : ```
89 :
90 : This function is used in the calculation of the Steinhardt order parameters, which are described in detail [here](https://www.plumed-tutorials.org/lessons/23/001/data/Steinhardt.html).
91 :
92 : */
93 : //+ENDPLUMEDOC
94 :
95 : class SphericalHarmonic {
96 : int tmom=0;
97 : std::vector<double> coeff_poly_v{};
98 : std::vector<double> normaliz_v{};
99 : double *coeff_poly=nullptr;
100 : double *normaliz=nullptr;
101 : public:
102 : static void registerKeywords( Keywords& keys );
103 : static unsigned factorial( unsigned n );
104 : static void read( SphericalHarmonic& func,
105 : ActionWithArguments* action,
106 : function::FunctionOptions& options );
107 : static void calc( const SphericalHarmonic& func,
108 : bool noderiv,
109 : View<const double> args,
110 : function::FunctionOutput funcout );
111 : double deriv_poly( unsigned m,
112 : double val,
113 : double& df ) const;
114 : static inline void addVectorDerivatives( unsigned ival,
115 : const Vector& der,
116 : View2D<double> derivatives );
117 : void update() {
118 42 : coeff_poly = coeff_poly_v.data();
119 42 : normaliz = normaliz_v.data();
120 : }
121 254 : SphericalHarmonic() = default;
122 170 : ~SphericalHarmonic() = default;
123 : SphericalHarmonic(const SphericalHarmonic&x):
124 : tmom(x.tmom),
125 : coeff_poly_v(x.coeff_poly_v),
126 : normaliz_v(x.normaliz_v) {
127 : update();
128 : }
129 : SphericalHarmonic(SphericalHarmonic&&x):
130 : tmom(x.tmom),
131 : coeff_poly_v(std::move(x.coeff_poly_v)),
132 : normaliz_v(std::move(x.normaliz_v)) {
133 : update();
134 : }
135 : SphericalHarmonic &operator=(const SphericalHarmonic&x) {
136 : if (this!=&x) {
137 : tmom=x.tmom;
138 : coeff_poly_v=x.coeff_poly_v;
139 : normaliz_v=x.normaliz_v;
140 : update();
141 : }
142 : return *this;
143 : }
144 : SphericalHarmonic &operator=(SphericalHarmonic&&x) {
145 : if (this!=&x) {
146 : tmom=x.tmom;
147 : coeff_poly_v=std::move(x.coeff_poly_v);
148 : normaliz_v=std::move(x.normaliz_v);
149 : update();
150 : }
151 : return *this;
152 : }
153 : #ifdef __PLUMED_USE_OPENACC
154 : void toACCDevice() const {
155 : const auto num=tmom+1;
156 : #pragma acc enter data copyin(this[0:1], \
157 : tmom, coeff_poly[0:num], normaliz[0:num])
158 : }
159 : void removeFromACCDevice() const {
160 : const auto num=tmom+1;
161 : #pragma acc exit data delete(normaliz[0:num], \
162 : coeff_poly[0:num], tmom, this[0:1])
163 : }
164 : #endif // __PLUMED_USE_OPENACC
165 : };
166 :
167 : typedef function::FunctionShortcut<SphericalHarmonic> SpHarmShortcut;
168 : PLUMED_REGISTER_ACTION(SpHarmShortcut,"SPHERICAL_HARMONIC")
169 : typedef function::FunctionOfMatrix<SphericalHarmonic> MatrixSpHarm;
170 : PLUMED_REGISTER_ACTION(MatrixSpHarm,"SPHERICAL_HARMONIC_MATRIX")
171 :
172 257 : void SphericalHarmonic::registerKeywords( Keywords& keys ) {
173 257 : keys.add("compulsory","L","the value of the angular momentum");
174 514 : keys.addOutputComponent("rm","default",
175 : "matrix","the real parts of the spherical harmonic values with the m value given");
176 514 : keys.addOutputComponent("im","default",
177 : "matrix","the real parts of the spherical harmonic values with the m value given");
178 257 : keys.add("hidden","MASKED_INPUT_ALLOWED",
179 : "turns on that you are allowed to use masked inputs");
180 257 : }
181 :
182 1858 : unsigned SphericalHarmonic::factorial( unsigned n ) {
183 1858 : return (n == 1 || n == 0) ? 1 : factorial(n - 1) * n;
184 : }
185 :
186 42 : void SphericalHarmonic::read( SphericalHarmonic& func,
187 : ActionWithArguments* action,
188 : function::FunctionOptions& options ) {
189 42 : action->parse("L",func.tmom);
190 42 : action->log.printf(" calculating %dth order spherical harmonic with %s, %s and %s as input \n",
191 : func.tmom,
192 : action->getPntrToArgument(0)->getName().c_str(),
193 : action->getPntrToArgument(1)->getName().c_str(),
194 : action->getPntrToArgument(2)->getName().c_str() );
195 42 : if( action->getNumberOfArguments()==4 ) {
196 42 : action->log.printf(" multiplying cylindrical harmonic by weight from %s \n",
197 : action->getPntrToArgument(3)->getName().c_str() );
198 : }
199 :
200 42 : func.normaliz_v.resize( func.tmom+1 );
201 42 : func.coeff_poly_v.resize( func.tmom+1 );
202 : func.update();
203 :
204 232 : for(int i=0; i<=func.tmom; ++i) {
205 190 : func.normaliz[i] = ( i%2==0 ? 1.0 : -1.0 )
206 570 : * sqrt( (2*func.tmom+1)
207 190 : * factorial(func.tmom-i)
208 190 : /(4*PLMD::pi*factorial(func.tmom+i)) );
209 : }
210 :
211 42 : switch (func.tmom) {
212 19 : case 1:
213 : // Legendre polynomial coefficients of order one
214 19 : func.coeff_poly[0]=0;
215 19 : func.coeff_poly[1]=1.0;
216 19 : break;
217 0 : case 2:
218 : // Legendre polynomial coefficients of order two
219 0 : func.coeff_poly[0]=-0.5;
220 0 : func.coeff_poly[1]=0.0;
221 0 : func.coeff_poly[2]=1.5;
222 0 : break;
223 1 : case 3:
224 : // Legendre polynomial coefficients of order three
225 1 : func.coeff_poly[0]=0.0;
226 1 : func.coeff_poly[1]=-1.5;
227 1 : func.coeff_poly[2]=0.0;
228 1 : func.coeff_poly[3]=2.5;
229 1 : break;
230 3 : case 4:
231 : // Legendre polynomial coefficients of order four
232 3 : func.coeff_poly[0]=0.375;
233 3 : func.coeff_poly[1]=0.0;
234 3 : func.coeff_poly[2]=-3.75;
235 3 : func.coeff_poly[3]=0.0;
236 3 : func.coeff_poly[4]=4.375;
237 3 : break;
238 0 : case 5:
239 : // Legendre polynomial coefficients of order five
240 0 : func.coeff_poly[0]=0.0;
241 0 : func.coeff_poly[1]=1.875;
242 0 : func.coeff_poly[2]=0.0;
243 0 : func.coeff_poly[3]=-8.75;
244 0 : func.coeff_poly[4]=0.0;
245 0 : func.coeff_poly[5]=7.875;
246 0 : break;
247 19 : case 6:
248 : // Legendre polynomial coefficients of order six
249 19 : func.coeff_poly[0]=-0.3125;
250 19 : func.coeff_poly[1]=0.0;
251 19 : func.coeff_poly[2]=6.5625;
252 19 : func.coeff_poly[3]=0.0;
253 19 : func.coeff_poly[4]=-19.6875;
254 19 : func.coeff_poly[5]=0.0;
255 19 : func.coeff_poly[6]=14.4375;
256 19 : break;
257 0 : default:
258 0 : action->error("Insert Legendre polynomial coefficients into SphericalHarmonics code");
259 : }
260 :
261 : std::string num;
262 380 : for(int i=-func.tmom; i<=func.tmom; ++i) {
263 338 : Tools::convert(std::abs(i),num);
264 338 : if( i<0 ) {
265 296 : options.multipleValuesForEachRegisteredComponent.push_back( "-n" + num );
266 190 : } else if( i>0 ) {
267 296 : options.multipleValuesForEachRegisteredComponent.push_back( "-p" + num );
268 : } else {
269 84 : options.multipleValuesForEachRegisteredComponent.push_back( "-0");
270 : }
271 : }
272 84 : options.derivativeZeroIfValueIsZero = (action->getNumberOfArguments()==4
273 42 : && action->getPntrToArgument(3)
274 : ->isDerivativeZeroWhenValueIsZero() );
275 42 : }
276 :
277 19212686 : void SphericalHarmonic::calc( const SphericalHarmonic& func,
278 : bool noderiv,
279 : const View<const double> args,
280 : function::FunctionOutput funcout ) {
281 : double weight=1;
282 19212686 : if( args.size()==4 ) {
283 19212686 : weight = args[3];
284 : }
285 19212686 : if( weight<epsilon ) {
286 16892014 : if( !noderiv ) {
287 3258340 : const unsigned imbase = 2*func.tmom+1;
288 27299632 : for(int m=-func.tmom; m<=func.tmom; ++m) {
289 24041292 : auto pos = func.tmom+m;
290 24041292 : funcout.derivs[pos][0] = 0.0;
291 24041292 : funcout.derivs[pos][1] = 0.0;
292 24041292 : funcout.derivs[pos][2] = 0.0;
293 24041292 : funcout.derivs[pos][3] = 0.0;
294 24041292 : funcout.derivs[pos + imbase][0] = 0.0;
295 24041292 : funcout.derivs[pos + imbase][1] = 0.0;
296 24041292 : funcout.derivs[pos + imbase][2] = 0.0;
297 24041292 : if( args.size()==4 ) {
298 24041292 : funcout.derivs[pos + imbase][3] = 0.0;
299 : }
300 : }
301 : }
302 16892014 : return;
303 : }
304 :
305 2320672 : const double dlen2 = args[0]*args[0] + args[1]*args[1] + args[2]*args[2];
306 2320672 : const double dlen = sqrt( dlen2 );
307 2320672 : const double dleninv = 1.0/dlen;
308 2320672 : const double dlen3 = dlen2*dlen;
309 2320672 : const double dlen3inv = 1.0/dlen3;
310 : double dpoly_ass;
311 2320672 : double poly_ass=func.deriv_poly( 0, args[2]*dleninv, dpoly_ass );
312 : // Accumulate for m=0
313 2320672 : funcout.values[func.tmom] = weight*poly_ass;
314 : // Derivatives of z/r wrt x, y, z
315 : Vector dz;
316 2320672 : if( !noderiv ) {
317 449462 : dz[0] = -( args[2] * dlen3inv )*args[0];
318 449462 : dz[1] = -( args[2] * dlen3inv )*args[1];
319 449462 : dz[2] = -( args[2] * dlen3inv )*args[2] + (1.0 * dleninv);
320 449462 : funcout.derivs[func.tmom] = weight*dpoly_ass*dz ;
321 449462 : if( args.size()==4 ) {
322 449462 : funcout.derivs[func.tmom][3] = poly_ass;
323 : }
324 : }
325 :
326 : // The complex number of which we have to take powers
327 2320672 : std::complex<double> com1( args[0]*dleninv, args[1]*dleninv );
328 : std::complex<double> powered(1.0, 0.0);
329 : constexpr std::complex<double> ii(0.0, 1.0);
330 : std::complex<double> dp_x;
331 : std::complex<double> dp_y;
332 : std::complex<double> dp_z;
333 : Vector myrealvec;
334 : Vector myimagvec;
335 : Vector real_dz;
336 : Vector imag_dz;
337 : // Do stuff for all other m values
338 2320672 : const unsigned imbase = 3*func.tmom+1;
339 : double pref = -1.0;
340 13805228 : for(int m=1; m<=func.tmom; ++m) {
341 : // Calculate Legendre Polynomial
342 11484556 : poly_ass=func.deriv_poly( m, args[2]*dleninv, dpoly_ass );
343 : // Real and imaginary parts of z
344 : double real_z = real(com1*powered);
345 : double imag_z = imag(com1*powered);
346 :
347 : // Calculate steinhardt parameter
348 11484556 : double tq6 =poly_ass*real_z; // Real part of steinhardt parameter
349 11484556 : double itq6=poly_ass*imag_z; // Imaginary part of steinhardt parameter
350 :
351 : // Real part
352 11484556 : funcout.values[func.tmom + m] = weight* tq6;
353 : // Imaginary part
354 11484556 : funcout.values[imbase + m] = weight*itq6;
355 : // Store -m part of vector
356 : // -m part of vector is just +m part multiplied by (-1.0)**m and multiplied by complex
357 : // conjugate of Legendre polynomial
358 : // Real part
359 11484556 : funcout.values[func.tmom - m] = pref*funcout.values[func.tmom+m];
360 : // Imaginary part
361 11484556 : funcout.values[imbase - m] = -pref*funcout.values[imbase + m];
362 11484556 : if( !noderiv ) {
363 : // Derivatives wrt ( x/r + iy )^m
364 1477034 : double md=static_cast<double>(m);
365 1477034 : dp_x = md*powered*( (1.0*dleninv)-(args[0]*args[0])*dlen3inv-ii*(args[0]*args[1])*dlen3inv );
366 1477034 : dp_y = md*powered*( ii*(1.0*dleninv)-(args[0]*args[1])*dlen3inv-ii*(args[1]*args[1])*dlen3inv );
367 1477034 : dp_z = md*powered*( -(args[0]*args[2])*dlen3inv-ii*(args[1]*args[2])*dlen3inv );
368 : // Derivatives of real and imaginary parts of above
369 : real_dz[0] = real( dp_x );
370 : real_dz[1] = real( dp_y );
371 : real_dz[2] = real( dp_z );
372 : imag_dz[0] = imag( dp_x );
373 : imag_dz[1] = imag( dp_y );
374 : imag_dz[2] = imag( dp_z );
375 : // Complete derivative of steinhardt parameter
376 1477034 : myrealvec = weight*(dpoly_ass*real_z*dz + poly_ass*real_dz);
377 1477034 : myimagvec = weight*(dpoly_ass*imag_z*dz + poly_ass*imag_dz);
378 1477034 : funcout.derivs[func.tmom+m] = myrealvec;
379 : funcout.derivs[imbase +m] = myimagvec;
380 :
381 1477034 : funcout.derivs[func.tmom-m] = pref*myrealvec;
382 : funcout.derivs[imbase -m] = -pref*myimagvec;
383 1477034 : if( args.size()==4 ) {
384 1477034 : funcout.derivs[func.tmom + m][3]= tq6;
385 1477034 : funcout.derivs[imbase + m][3]= itq6;
386 1477034 : funcout.derivs[func.tmom - m][3]= pref* tq6;
387 1477034 : funcout.derivs[imbase - m][3]=-pref*itq6;
388 : }
389 : }
390 : // Calculate next power of complex number
391 : powered *= com1;
392 : //hopefully the compiler will optimize with bitflipping the sign here
393 : pref = -pref;
394 : }
395 : }
396 :
397 13805228 : double SphericalHarmonic::deriv_poly( const unsigned m,
398 : const double val,
399 : double& df ) const {
400 : double fact=1.0;
401 38933216 : for(unsigned j=2; j<=m; ++j) {
402 25127988 : fact *= j;
403 : }
404 13805228 : double res=coeff_poly[m]*fact;
405 :
406 : double pow=1.0;
407 : double xi=val;
408 : double dxi=1.0;
409 13805228 : df=0.0;
410 50417772 : for(int i=m+1; i<=tmom; ++i) {
411 : fact = 1.0;
412 92864124 : for(int j=i-m+1; j<=i; ++j) {
413 56251580 : fact *= j;
414 : }
415 36612544 : res += coeff_poly[i]*fact*xi;
416 36612544 : df += pow*coeff_poly[i]*fact*dxi;
417 36612544 : xi *= val;
418 36612544 : dxi *= val;
419 36612544 : pow += 1.0;
420 : }
421 13805228 : df *= normaliz[m];
422 13805228 : return normaliz[m]*res;
423 : }
424 :
425 : }
426 : }
427 :
|
