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 : #ifndef __PLUMED_symfunc_SphericalHarmonic_h
23 : #define __PLUMED_symfunc_SphericalHarmonic_h
24 : #include "function/FunctionSetup.h"
25 : #include "function/FunctionOfMatrix.h"
26 :
27 : #include <complex>
28 :
29 : namespace PLMD {
30 : namespace symfunc {
31 :
32 : class SphericalHarmonic {
33 : int tmom=0;
34 : std::vector<double> coeff_poly_v{};
35 : std::vector<double> normaliz_v{};
36 : double *coeff_poly=nullptr;
37 : double *normaliz=nullptr;
38 : public:
39 : static void registerKeywords( Keywords& keys );
40 : static unsigned factorial( unsigned n );
41 : static void read( SphericalHarmonic& func,
42 : ActionWithArguments* action,
43 : function::FunctionOptions& options );
44 : static void calc( const SphericalHarmonic& func,
45 : bool noderiv,
46 : View<const double> args,
47 : function::FunctionOutput funcout );
48 : double deriv_poly( unsigned m,
49 : double val,
50 : double& df ) const;
51 : static inline void addVectorDerivatives( unsigned ival,
52 : const Vector& der,
53 : View2D<double> derivatives );
54 : void update() {
55 44 : coeff_poly = coeff_poly_v.data();
56 44 : normaliz = normaliz_v.data();
57 : }
58 90 : SphericalHarmonic() = default;
59 1 : ~SphericalHarmonic() = default;
60 : SphericalHarmonic(const SphericalHarmonic&x):
61 : tmom(x.tmom),
62 : coeff_poly_v(x.coeff_poly_v),
63 : normaliz_v(x.normaliz_v) {
64 : update();
65 : }
66 : SphericalHarmonic(SphericalHarmonic&&x):
67 : tmom(x.tmom),
68 : coeff_poly_v(std::move(x.coeff_poly_v)),
69 : normaliz_v(std::move(x.normaliz_v)) {
70 : update();
71 : }
72 : SphericalHarmonic &operator=(const SphericalHarmonic&x) {
73 : if (this!=&x) {
74 : tmom=x.tmom;
75 : coeff_poly_v=x.coeff_poly_v;
76 : normaliz_v=x.normaliz_v;
77 : update();
78 : }
79 : return *this;
80 : }
81 : SphericalHarmonic &operator=(SphericalHarmonic&&x) {
82 : if (this!=&x) {
83 : tmom=x.tmom;
84 : coeff_poly_v=std::move(x.coeff_poly_v);
85 : normaliz_v=std::move(x.normaliz_v);
86 : update();
87 : }
88 : return *this;
89 : }
90 : #ifdef __PLUMED_HAS_OPENACC
91 : void toACCDevice() const {
92 : const auto num=tmom+1;
93 : #pragma acc enter data copyin(this[0:1], \
94 : tmom, coeff_poly[0:num], normaliz[0:num])
95 : }
96 : void removeFromACCDevice() const {
97 : const auto num=tmom+1;
98 : #pragma acc exit data delete(normaliz[0:num], \
99 : coeff_poly[0:num], tmom, this[0:1])
100 : }
101 : #endif // __PLUMED_HAS_OPENACC
102 : };
103 :
104 274 : void SphericalHarmonic::registerKeywords( Keywords& keys ) {
105 274 : keys.add("compulsory","L","the value of the angular momentum");
106 548 : keys.addOutputComponent("rm","default",
107 : "matrix","the real parts of the spherical harmonic values with the m value given");
108 548 : keys.addOutputComponent("im","default",
109 : "matrix","the real parts of the spherical harmonic values with the m value given");
110 274 : keys.add("hidden","MASKED_INPUT_ALLOWED",
111 : "turns on that you are allowed to use masked inputs");
112 274 : }
113 :
114 2028 : unsigned SphericalHarmonic::factorial( unsigned n ) {
115 2028 : return (n == 1 || n == 0) ? 1 : factorial(n - 1) * n;
116 : }
117 :
118 45 : void SphericalHarmonic::read( SphericalHarmonic& func,
119 : ActionWithArguments* action,
120 : function::FunctionOptions& options ) {
121 45 : action->parse("L",func.tmom);
122 45 : if( action->getNumberOfArguments()<3 ) {
123 1 : action->error("the x, y and z components of the vectors (and optionally a weight) should be given with ARG");
124 : }
125 44 : action->log.printf(" calculating %dth order spherical harmonic with %s, %s and %s as input \n",
126 : func.tmom,
127 : action->getPntrToArgument(0)->getName().c_str(),
128 : action->getPntrToArgument(1)->getName().c_str(),
129 : action->getPntrToArgument(2)->getName().c_str() );
130 44 : if( action->getNumberOfArguments()==4 ) {
131 44 : action->log.printf(" multiplying cylindrical harmonic by weight from %s \n",
132 : action->getPntrToArgument(3)->getName().c_str() );
133 : }
134 :
135 44 : func.normaliz_v.resize( func.tmom+1 );
136 44 : func.coeff_poly_v.resize( func.tmom+1 );
137 : func.update();
138 :
139 248 : for(decltype(func.tmom) i=0; i<=func.tmom; ++i) {
140 204 : func.normaliz[i] = ( i%2==0 ? 1.0 : -1.0 )
141 612 : * sqrt( (2*func.tmom+1)
142 204 : * factorial(func.tmom-i)
143 204 : /(4*PLMD::pi*factorial(func.tmom+i)) );
144 : }
145 :
146 44 : switch (func.tmom) {
147 19 : case 1:
148 : // Legendre polynomial coefficients of order one
149 19 : func.coeff_poly[0]=0;
150 19 : func.coeff_poly[1]=1.0;
151 19 : break;
152 0 : case 2:
153 : // Legendre polynomial coefficients of order two
154 0 : func.coeff_poly[0]=-0.5;
155 0 : func.coeff_poly[1]=0.0;
156 0 : func.coeff_poly[2]=1.5;
157 0 : break;
158 1 : case 3:
159 : // Legendre polynomial coefficients of order three
160 1 : func.coeff_poly[0]=0.0;
161 1 : func.coeff_poly[1]=-1.5;
162 1 : func.coeff_poly[2]=0.0;
163 1 : func.coeff_poly[3]=2.5;
164 1 : break;
165 3 : case 4:
166 : // Legendre polynomial coefficients of order four
167 3 : func.coeff_poly[0]=0.375;
168 3 : func.coeff_poly[1]=0.0;
169 3 : func.coeff_poly[2]=-3.75;
170 3 : func.coeff_poly[3]=0.0;
171 3 : func.coeff_poly[4]=4.375;
172 3 : break;
173 0 : case 5:
174 : // Legendre polynomial coefficients of order five
175 0 : func.coeff_poly[0]=0.0;
176 0 : func.coeff_poly[1]=1.875;
177 0 : func.coeff_poly[2]=0.0;
178 0 : func.coeff_poly[3]=-8.75;
179 0 : func.coeff_poly[4]=0.0;
180 0 : func.coeff_poly[5]=7.875;
181 0 : break;
182 21 : case 6:
183 : // Legendre polynomial coefficients of order six
184 21 : func.coeff_poly[0]=-0.3125;
185 21 : func.coeff_poly[1]=0.0;
186 21 : func.coeff_poly[2]=6.5625;
187 21 : func.coeff_poly[3]=0.0;
188 21 : func.coeff_poly[4]=-19.6875;
189 21 : func.coeff_poly[5]=0.0;
190 21 : func.coeff_poly[6]=14.4375;
191 21 : break;
192 0 : default:
193 0 : action->error("Insert Legendre polynomial coefficients into SphericalHarmonics code");
194 : }
195 :
196 : std::string num;
197 408 : for(int i=-func.tmom; i<=func.tmom; ++i) {
198 364 : Tools::convert(abs(i),num);
199 364 : if( i<0 ) {
200 320 : options.multipleValuesForEachRegisteredComponent.push_back( "-n" + num );
201 204 : } else if( i>0 ) {
202 320 : options.multipleValuesForEachRegisteredComponent.push_back( "-p" + num );
203 : } else {
204 88 : options.multipleValuesForEachRegisteredComponent.push_back( "-0");
205 : }
206 : }
207 88 : options.derivativeZeroIfValueIsZero = (action->getNumberOfArguments()==4
208 44 : && action->getPntrToArgument(3)
209 : ->isDerivativeZeroWhenValueIsZero() );
210 44 : }
211 :
212 2339822 : void SphericalHarmonic::calc( const SphericalHarmonic& func,
213 : bool noderiv,
214 : const View<const double> args,
215 : function::FunctionOutput funcout ) {
216 : double weight=1;
217 2339822 : if( args.size()==4 ) {
218 2339822 : weight = args[3];
219 : }
220 2339822 : if( weight<epsilon ) {
221 10010 : if( !noderiv ) {
222 5000 : const unsigned imbase = 2*func.tmom+1;
223 69760 : for(int m=-func.tmom; m<=func.tmom; ++m) {
224 64760 : auto pos = func.tmom+m;
225 64760 : funcout.derivs[pos][0] = 0.0;
226 64760 : funcout.derivs[pos][1] = 0.0;
227 64760 : funcout.derivs[pos][2] = 0.0;
228 64760 : funcout.derivs[pos][3] = 0.0;
229 64760 : funcout.derivs[pos + imbase][0] = 0.0;
230 64760 : funcout.derivs[pos + imbase][1] = 0.0;
231 64760 : funcout.derivs[pos + imbase][2] = 0.0;
232 64760 : if( args.size()==4 ) {
233 64760 : funcout.derivs[pos + imbase][3] = 0.0;
234 : }
235 : }
236 : }
237 10010 : return;
238 : }
239 :
240 2329812 : const double dlen2 = args[0]*args[0] + args[1]*args[1] + args[2]*args[2];
241 2329812 : const double dlen = sqrt( dlen2 );
242 2329812 : const double dleninv = 1.0/dlen;
243 2329812 : const double dlen3 = dlen2*dlen;
244 2329812 : const double dlen3inv = 1.0/dlen3;
245 : double dpoly_ass;
246 2329812 : double poly_ass=func.deriv_poly( 0, args[2]*dleninv, dpoly_ass );
247 : // Accumulate for m=0
248 2329812 : funcout.values[func.tmom] = weight*poly_ass;
249 : // Derivatives of z/r wrt x, y, z
250 : Vector dz;
251 2329812 : if( !noderiv ) {
252 453494 : dz[0] = -( args[2] * dlen3inv )*args[0];
253 453494 : dz[1] = -( args[2] * dlen3inv )*args[1];
254 453494 : dz[2] = -( args[2] * dlen3inv )*args[2] + (1.0 * dleninv);
255 453494 : funcout.derivs[func.tmom] = weight*dpoly_ass*dz ;
256 453494 : if( args.size()==4 ) {
257 453494 : funcout.derivs[func.tmom][3] = poly_ass;
258 : }
259 : }
260 :
261 : // The complex number of which we have to take powers
262 2329812 : std::complex<double> com1( args[0]*dleninv, args[1]*dleninv );
263 : std::complex<double> powered(1.0, 0.0);
264 : constexpr std::complex<double> ii(0.0, 1.0);
265 : std::complex<double> dp_x;
266 : std::complex<double> dp_y;
267 : std::complex<double> dp_z;
268 : Vector myrealvec;
269 : Vector myimagvec;
270 : Vector real_dz;
271 : Vector imag_dz;
272 : // Do stuff for all other m values
273 2329812 : const unsigned imbase = 3*func.tmom+1;
274 : // The imaginary part for m=0 is always zero, but it is a component like the
275 : // others and its derivatives are read when the forces are applied: set them
276 : // here rather than relying on the buffer being zero (on the GPU it is not)
277 2329812 : funcout.values[imbase] = 0.0;
278 2329812 : if( !noderiv ) {
279 453494 : funcout.derivs[imbase][0] = 0.0;
280 453494 : funcout.derivs[imbase][1] = 0.0;
281 453494 : funcout.derivs[imbase][2] = 0.0;
282 453494 : if( args.size()==4 ) {
283 453494 : funcout.derivs[imbase][3] = 0.0;
284 : }
285 : }
286 : double pref = -1.0;
287 13869208 : for(decltype(func.tmom) m=1; m<=func.tmom; ++m) {
288 : // Calculate Legendre Polynomial
289 11539396 : poly_ass=func.deriv_poly( m, args[2]*dleninv, dpoly_ass );
290 : // Real and imaginary parts of z
291 : double real_z = real(com1*powered);
292 : double imag_z = imag(com1*powered);
293 :
294 : // Calculate steinhardt parameter
295 11539396 : double tq6 =poly_ass*real_z; // Real part of steinhardt parameter
296 11539396 : double itq6=poly_ass*imag_z; // Imaginary part of steinhardt parameter
297 :
298 : // Real part
299 11539396 : funcout.values[func.tmom + m] = weight* tq6;
300 : // Imaginary part
301 11539396 : funcout.values[imbase + m] = weight*itq6;
302 : // Store -m part of vector
303 : // -m part of vector is just +m part multiplied by (-1.0)**m and multiplied by complex
304 : // conjugate of Legendre polynomial
305 : // Real part
306 11539396 : funcout.values[func.tmom - m] = pref*funcout.values[func.tmom+m];
307 : // Imaginary part
308 11539396 : funcout.values[imbase - m] = -pref*funcout.values[imbase + m];
309 11539396 : if( !noderiv ) {
310 : // Derivatives wrt ( x/r + iy )^m
311 1501226 : double md=static_cast<double>(m);
312 1501226 : dp_x = md*powered*( (1.0*dleninv)-(args[0]*args[0])*dlen3inv-ii*(args[0]*args[1])*dlen3inv );
313 1501226 : dp_y = md*powered*( ii*(1.0*dleninv)-(args[0]*args[1])*dlen3inv-ii*(args[1]*args[1])*dlen3inv );
314 1501226 : dp_z = md*powered*( -(args[0]*args[2])*dlen3inv-ii*(args[1]*args[2])*dlen3inv );
315 : // Derivatives of real and imaginary parts of above
316 : real_dz[0] = real( dp_x );
317 : real_dz[1] = real( dp_y );
318 : real_dz[2] = real( dp_z );
319 : imag_dz[0] = imag( dp_x );
320 : imag_dz[1] = imag( dp_y );
321 : imag_dz[2] = imag( dp_z );
322 : // Complete derivative of steinhardt parameter
323 1501226 : myrealvec = weight*(dpoly_ass*real_z*dz + poly_ass*real_dz);
324 1501226 : myimagvec = weight*(dpoly_ass*imag_z*dz + poly_ass*imag_dz);
325 1501226 : funcout.derivs[func.tmom+m] = myrealvec;
326 : funcout.derivs[imbase +m] = myimagvec;
327 :
328 1501226 : funcout.derivs[func.tmom-m] = pref*myrealvec;
329 : funcout.derivs[imbase -m] = -pref*myimagvec;
330 1501226 : if( args.size()==4 ) {
331 1501226 : funcout.derivs[func.tmom + m][3]= tq6;
332 1501226 : funcout.derivs[imbase + m][3]= itq6;
333 1501226 : funcout.derivs[func.tmom - m][3]= pref* tq6;
334 1501226 : funcout.derivs[imbase - m][3]=-pref*itq6;
335 : }
336 : }
337 : // Calculate next power of complex number
338 : powered *= com1;
339 : //hopefully the compiler will optimize with bitflipping the sign here
340 : pref = -pref;
341 : }
342 : }
343 :
344 13869208 : double SphericalHarmonic::deriv_poly( const unsigned m,
345 : const double val,
346 : double& df ) const {
347 : double fact=1.0;
348 39134296 : for(unsigned j=2; j<=m; ++j) {
349 25265088 : fact *= j;
350 : }
351 13869208 : double res=coeff_poly[m]*fact;
352 :
353 : double pow=1.0;
354 : double xi=val;
355 : double dxi=1.0;
356 13869208 : df=0.0;
357 50673692 : for(int i=m+1; i<=tmom; ++i) {
358 : fact = 1.0;
359 93375964 : for(int j=i-m+1; j<=i; ++j) {
360 56571480 : fact *= j;
361 : }
362 36804484 : res += coeff_poly[i]*fact*xi;
363 36804484 : df += pow*coeff_poly[i]*fact*dxi;
364 36804484 : xi *= val;
365 36804484 : dxi *= val;
366 36804484 : pow += 1.0;
367 : }
368 13869208 : df *= normaliz[m];
369 13869208 : return normaliz[m]*res;
370 : }
371 :
372 : }
373 : }
374 :
375 : #endif //__PLUMED_symfunc_SphericalHarmonic_h
|
