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 : #include "function/FunctionTemplateBase.h"
23 : #include "function/FunctionShortcut.h"
24 : #include "function/FunctionOfMatrix.h"
25 : #include "core/ActionRegister.h"
26 :
27 : #include <complex>
28 :
29 : namespace PLMD {
30 : namespace symfunc {
31 :
32 : //+PLUMEDOC MCOLVAR SPHERICAL_HARMONIC
33 : /*
34 : Calculate the values of all the spherical harmonic funtions for a particular value of l.
35 :
36 : This action allows you to the all the [spherical harmonic](https://en.wikipedia.org/wiki/Spherical_harmonics) functions for a particular input
37 : $L$ value. As discussed in more detail in the article provided above the spherical harmonics this action calculates have the following general form:
38 :
39 : $$
40 : Y_l^m(\theta,\phi) = wNe^{im\phi} P_l^m(\cos\theta)
41 : $$
42 :
43 : 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.
44 : $e^{i\phi}$ and $\cos(\theta) are computed from the other input arguments, $x, y$ and $z$ as follows:
45 :
46 : $$
47 : 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}
48 : $$
49 :
50 : 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.
51 : The arguments must all have the same shape as the two output components will also be matrices that are
52 : 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
53 : the real and imaginary parts of the $Y_l^m$ functions with the the $2l+1$ possible $m$ values.
54 :
55 : The following intput provides an example that demonstrates how this function is used:
56 :
57 : ```plumed
58 : d: DISTANCE_MATRIX GROUP=1-10 COMPONENTS
59 : c: SPHERICAL_HARMONIC L=1 ARG=d.x,d.y,d.z
60 : PRINT ARG=c.rm-n1 FILE=real_part_m-1
61 : PRINT ARG=c.im-n1 FILE=imaginary_part_m-1
62 : PRINT ARG=c.rm-0 FILE=real_part_m0
63 : PRINT ARG=c.im-0 FILE=imaginary_part_m0
64 : PRINT ARG=c.rm-p1 FILE=real_part_m+1
65 : PRINT ARG=c.im-p1 FILE=imaginary_part_m+1
66 : ```
67 :
68 : 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,
69 : 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$
70 : matrices are then output to six separate files.
71 :
72 : 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
73 : above.
74 :
75 : ```plumed
76 : s: CONTACT_MATRIX GROUP=1-10 SWITCH={RATIONAL R_0=1.0}
77 : sc: CONTACT_MATRIX GROUP=1-10 SWITCH={RATIONAL R_0=1.0} COMPONENTS
78 : c: SPHERICAL_HARMONIC L=1 ARG=sc.x,sc.y,sc.z,s
79 : PRINT ARG=c.rm-n1 FILE=real_part_m-1
80 : PRINT ARG=c.im-n1 FILE=imaginary_part_m-1
81 : PRINT ARG=c.rm-0 FILE=real_part_m0
82 : PRINT ARG=c.im-0 FILE=imaginary_part_m0
83 : PRINT ARG=c.rm-p1 FILE=real_part_m+1
84 : PRINT ARG=c.im-p1 FILE=imaginary_part_m+1
85 : ```
86 :
87 : 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).
88 :
89 : */
90 : //+ENDPLUMEDOC
91 :
92 : //+PLUMEDOC MCOLVAR SPHERICAL_HARMONIC_MATRIX
93 : /*
94 : Calculate the values of all the spherical harmonic funtions for a particular value of l for all the elements of a set of three input matrices
95 :
96 : \par Examples
97 :
98 :
99 : */
100 : //+ENDPLUMEDOC
101 :
102 299 : class SphericalHarmonic : public function::FunctionTemplateBase {
103 : private:
104 : int tmom;
105 : std::vector<double> coeff_poly;
106 : std::vector<double> normaliz;
107 : unsigned factorial( const unsigned& n ) const ;
108 : double deriv_poly( const unsigned& m, const double& val, double& df ) const ;
109 : void addVectorDerivatives( const unsigned& ival, const Vector& der, Matrix<double>& derivatives ) const ;
110 : public:
111 : void registerKeywords( Keywords& keys ) override;
112 : void read( ActionWithArguments* action ) override;
113 : std::vector<std::string> getComponentsPerLabel() const override;
114 : void setPeriodicityForOutputs( ActionWithValue* action ) override;
115 : void calc( const ActionWithArguments* action, const std::vector<double>& args, std::vector<double>& vals, Matrix<double>& derivatives ) const override;
116 : };
117 :
118 : typedef function::FunctionShortcut<SphericalHarmonic> SpHarmShortcut;
119 : PLUMED_REGISTER_ACTION(SpHarmShortcut,"SPHERICAL_HARMONIC")
120 : typedef function::FunctionOfMatrix<SphericalHarmonic> MatrixSpHarm;
121 : PLUMED_REGISTER_ACTION(MatrixSpHarm,"SPHERICAL_HARMONIC_MATRIX")
122 :
123 215 : void SphericalHarmonic::registerKeywords( Keywords& keys ) {
124 215 : keys.add("compulsory","L","the value of the angular momentum");
125 430 : keys.addOutputComponent("rm","default","matrix","the real parts of the spherical harmonic values with the m value given");
126 430 : keys.addOutputComponent("im","default","matrix","the real parts of the spherical harmonic values with the m value given");
127 215 : }
128 :
129 1858 : unsigned SphericalHarmonic::factorial( const unsigned& n ) const {
130 1858 : return (n == 1 || n == 0) ? 1 : factorial(n - 1) * n;
131 : }
132 :
133 42 : void SphericalHarmonic::read( ActionWithArguments* action ) {
134 42 : parse(action,"L",tmom);
135 42 : action->log.printf(" calculating %dth order spherical harmonic with %s, %s and %s as input \n", tmom, action->getPntrToArgument(0)->getName().c_str(), action->getPntrToArgument(1)->getName().c_str(), action->getPntrToArgument(2)->getName().c_str() );
136 42 : if( action->getNumberOfArguments()==4 ) {
137 42 : action->log.printf(" multiplying cylindrical harmonic by weight from %s \n", action->getPntrToArgument(3)->getName().c_str() );
138 : }
139 :
140 42 : normaliz.resize( tmom+1 );
141 232 : for(unsigned i=0; i<=tmom; ++i) {
142 190 : normaliz[i] = sqrt( (2*tmom+1)*factorial(tmom-i)/(4*pi*factorial(tmom+i)) );
143 190 : if( i%2==1 ) {
144 84 : normaliz[i]*=-1;
145 : }
146 : }
147 :
148 42 : coeff_poly.resize( tmom+1 );
149 42 : if( tmom==1 ) {
150 : // Legendre polynomial coefficients of order one
151 19 : coeff_poly[0]=0;
152 19 : coeff_poly[1]=1.0;
153 : } else if( tmom==2 ) {
154 : // Legendre polynomial coefficients of order two
155 0 : coeff_poly[0]=-0.5;
156 0 : coeff_poly[1]=0.0;
157 0 : coeff_poly[2]=1.5;
158 : } else if( tmom==3 ) {
159 : // Legendre polynomial coefficients of order three
160 1 : coeff_poly[0]=0.0;
161 1 : coeff_poly[1]=-1.5;
162 1 : coeff_poly[2]=0.0;
163 1 : coeff_poly[3]=2.5;
164 : } else if( tmom==4 ) {
165 : // Legendre polynomial coefficients of order four
166 3 : coeff_poly[0]=0.375;
167 3 : coeff_poly[1]=0.0;
168 3 : coeff_poly[2]=-3.75;
169 3 : coeff_poly[3]=0.0;
170 3 : coeff_poly[4]=4.375;
171 : } else if( tmom==5 ) {
172 : // Legendre polynomial coefficients of order five
173 0 : coeff_poly[0]=0.0;
174 0 : coeff_poly[1]=1.875;
175 0 : coeff_poly[2]=0.0;
176 0 : coeff_poly[3]=-8.75;
177 0 : coeff_poly[4]=0.0;
178 0 : coeff_poly[5]=7.875;
179 : } else if( tmom==6 ) {
180 : // Legendre polynomial coefficients of order six
181 19 : coeff_poly[0]=-0.3125;
182 19 : coeff_poly[1]=0.0;
183 19 : coeff_poly[2]=6.5625;
184 19 : coeff_poly[3]=0.0;
185 19 : coeff_poly[4]=-19.6875;
186 19 : coeff_poly[5]=0.0;
187 19 : coeff_poly[6]=14.4375;
188 : } else {
189 0 : action->error("Insert Legendre polynomial coefficients into SphericalHarmonics code");
190 : }
191 42 : }
192 :
193 84 : std::vector<std::string> SphericalHarmonic::getComponentsPerLabel() const {
194 : std::vector<std::string> comp;
195 : std::string num;
196 760 : for(int i=-tmom; i<=tmom; ++i) {
197 676 : Tools::convert(fabs(i),num);
198 676 : if( i<0 ) {
199 592 : comp.push_back( "-n" + num );
200 380 : } else if( i>0 ) {
201 592 : comp.push_back( "-p" + num );
202 : } else {
203 168 : comp.push_back( "-0");
204 : }
205 : }
206 84 : return comp;
207 0 : }
208 :
209 42 : void SphericalHarmonic::setPeriodicityForOutputs( ActionWithValue* action ) {
210 42 : std::vector<std::string> comp( getComponentsPerLabel() );
211 380 : for(unsigned i=0; i<comp.size(); ++i) {
212 676 : action->componentIsNotPeriodic("rm" + comp[i]);
213 676 : action->componentIsNotPeriodic("im" + comp[i]);
214 : }
215 42 : }
216 :
217 1742755 : void SphericalHarmonic::calc( const ActionWithArguments* action, const std::vector<double>& args, std::vector<double>& vals, Matrix<double>& derivatives ) const {
218 : double weight=1;
219 1742755 : if( args.size()==4 ) {
220 1742755 : weight = args[3];
221 : }
222 1742755 : if( weight<epsilon ) {
223 1468550 : return;
224 : }
225 :
226 274205 : double dlen2 = args[0]*args[0]+args[1]*args[1]+args[2]*args[2];
227 274205 : double dlen = sqrt( dlen2 );
228 274205 : double dlen3 = dlen2*dlen;
229 274205 : double dpoly_ass, poly_ass=deriv_poly( 0, args[2]/dlen, dpoly_ass );
230 : // Derivatives of z/r wrt x, y, z
231 274205 : Vector dz;
232 274205 : dz[0] = -( args[2] / dlen3 )*args[0];
233 274205 : dz[1] = -( args[2] / dlen3 )*args[1];
234 274205 : dz[2] = -( args[2] / dlen3 )*args[2] + (1.0 / dlen);
235 : // Accumulate for m=0
236 274205 : vals[tmom] = weight*poly_ass;
237 274205 : addVectorDerivatives( tmom, weight*dpoly_ass*dz, derivatives );
238 274205 : if( args.size()==4 ) {
239 274205 : derivatives(tmom, 3) = poly_ass;
240 : }
241 :
242 : // The complex number of which we have to take powers
243 274205 : std::complex<double> com1( args[0]/dlen,args[1]/dlen ), dp_x, dp_y, dp_z;
244 : std::complex<double> powered = std::complex<double>(1.0,0.0);
245 : std::complex<double> ii( 0.0, 1.0 );
246 274205 : Vector myrealvec, myimagvec, real_dz, imag_dz;
247 : // Do stuff for all other m values
248 1794050 : for(unsigned m=1; m<=tmom; ++m) {
249 : // Calculate Legendre Polynomial
250 1519845 : poly_ass=deriv_poly( m, args[2]/dlen, dpoly_ass );
251 : // Real and imaginary parts of z
252 : double real_z = real(com1*powered), imag_z = imag(com1*powered);
253 :
254 : // Calculate steinhardt parameter
255 1519845 : double tq6=poly_ass*real_z; // Real part of steinhardt parameter
256 1519845 : double itq6=poly_ass*imag_z; // Imaginary part of steinhardt parameter
257 :
258 : // Derivatives wrt ( x/r + iy )^m
259 1519845 : double md=static_cast<double>(m);
260 1519845 : dp_x = md*powered*( (1.0/dlen)-(args[0]*args[0])/dlen3-ii*(args[0]*args[1])/dlen3 );
261 1519845 : dp_y = md*powered*( ii*(1.0/dlen)-(args[0]*args[1])/dlen3-ii*(args[1]*args[1])/dlen3 );
262 1519845 : dp_z = md*powered*( -(args[0]*args[2])/dlen3-ii*(args[1]*args[2])/dlen3 );
263 :
264 : // Derivatives of real and imaginary parts of above
265 1519845 : real_dz[0] = real( dp_x );
266 1519845 : real_dz[1] = real( dp_y );
267 1519845 : real_dz[2] = real( dp_z );
268 1519845 : imag_dz[0] = imag( dp_x );
269 1519845 : imag_dz[1] = imag( dp_y );
270 1519845 : imag_dz[2] = imag( dp_z );
271 :
272 : // Complete derivative of steinhardt parameter
273 1519845 : myrealvec = weight*dpoly_ass*real_z*dz + weight*poly_ass*real_dz;
274 1519845 : myimagvec = weight*dpoly_ass*imag_z*dz + weight*poly_ass*imag_dz;
275 :
276 : // Real part
277 1519845 : vals[tmom+m] = weight*tq6;
278 1519845 : addVectorDerivatives( tmom+m, myrealvec, derivatives );
279 : // Imaginary part
280 1519845 : vals[3*tmom+1+m] = weight*itq6;
281 1519845 : addVectorDerivatives( 3*tmom+1+m, myimagvec, derivatives );
282 : // Store -m part of vector
283 1519845 : double pref=pow(-1.0,m);
284 : // -m part of vector is just +m part multiplied by (-1.0)**m and multiplied by complex
285 : // conjugate of Legendre polynomial
286 : // Real part
287 1519845 : vals[tmom-m] = pref*weight*tq6;
288 1519845 : addVectorDerivatives( tmom-m, pref*myrealvec, derivatives );
289 : // Imaginary part
290 1519845 : vals[3*tmom+1-m] = -pref*weight*itq6;
291 1519845 : addVectorDerivatives( 3*tmom+1-m, -pref*myimagvec, derivatives );
292 1519845 : if( args.size()==4 ) {
293 1519845 : derivatives(tmom+m,3)=tq6;
294 1519845 : derivatives(3*tmom+1+m, 3)=itq6;
295 1519845 : derivatives(tmom-m,3)=pref*tq6;
296 1519845 : derivatives(3*tmom+1-m, 3)=-pref*itq6;
297 : }
298 : // Calculate next power of complex number
299 : powered *= com1;
300 : }
301 : }
302 :
303 1794050 : double SphericalHarmonic::deriv_poly( const unsigned& m, const double& val, double& df ) const {
304 : double fact=1.0;
305 7004030 : for(unsigned j=1; j<=m; ++j) {
306 5209980 : fact=fact*j;
307 : }
308 1794050 : double res=coeff_poly[m]*fact;
309 :
310 1794050 : double pow=1.0, xi=val, dxi=1.0;
311 1794050 : df=0.0;
312 7004030 : for(int i=m+1; i<=tmom; ++i) {
313 : fact=1.0;
314 13759570 : for(unsigned j=i-m+1; j<=i; ++j) {
315 8549590 : fact=fact*j;
316 : }
317 5209980 : res=res+coeff_poly[i]*fact*xi;
318 5209980 : df = df + pow*coeff_poly[i]*fact*dxi;
319 5209980 : xi=xi*val;
320 5209980 : dxi=dxi*val;
321 5209980 : pow+=1.0;
322 : }
323 1794050 : df = df*normaliz[m];
324 1794050 : return normaliz[m]*res;
325 : }
326 :
327 6353585 : void SphericalHarmonic::addVectorDerivatives( const unsigned& ival, const Vector& der, Matrix<double>& derivatives ) const {
328 25414340 : for(unsigned j=0; j<3; ++j) {
329 19060755 : derivatives(ival,j) = der[j];
330 : }
331 6353585 : }
332 :
333 : }
334 : }
335 :
|
