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 : \par Examples
37 :
38 :
39 : */
40 : //+ENDPLUMEDOC
41 :
42 : //+PLUMEDOC MCOLVAR SPHERICAL_HARMONIC_MATRIX
43 : /*
44 : 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
45 :
46 : \par Examples
47 :
48 :
49 : */
50 : //+ENDPLUMEDOC
51 :
52 299 : class SphericalHarmonic : public function::FunctionTemplateBase {
53 : private:
54 : int tmom;
55 : std::vector<double> coeff_poly;
56 : std::vector<double> normaliz;
57 : unsigned factorial( const unsigned& n ) const ;
58 : double deriv_poly( const unsigned& m, const double& val, double& df ) const ;
59 : void addVectorDerivatives( const unsigned& ival, const Vector& der, Matrix<double>& derivatives ) const ;
60 : public:
61 : void registerKeywords( Keywords& keys ) override;
62 : void read( ActionWithArguments* action ) override;
63 : std::vector<std::string> getComponentsPerLabel() const override;
64 : void setPeriodicityForOutputs( ActionWithValue* action ) override;
65 : void calc( const ActionWithArguments* action, const std::vector<double>& args, std::vector<double>& vals, Matrix<double>& derivatives ) const override;
66 : };
67 :
68 : typedef function::FunctionShortcut<SphericalHarmonic> SpHarmShortcut;
69 : PLUMED_REGISTER_ACTION(SpHarmShortcut,"SPHERICAL_HARMONIC")
70 : typedef function::FunctionOfMatrix<SphericalHarmonic> MatrixSpHarm;
71 : PLUMED_REGISTER_ACTION(MatrixSpHarm,"SPHERICAL_HARMONIC_MATRIX")
72 :
73 215 : void SphericalHarmonic::registerKeywords( Keywords& keys ) {
74 430 : keys.add("compulsory","L","the value of the angular momentum");
75 430 : keys.addOutputComponent("rm","default","the real parts of the spherical harmonic values with the m value given");
76 430 : keys.addOutputComponent("im","default","the real parts of the spherical harmonic values with the m value given");
77 215 : }
78 :
79 1858 : unsigned SphericalHarmonic::factorial( const unsigned& n ) const {
80 1858 : return (n == 1 || n == 0) ? 1 : factorial(n - 1) * n;
81 : }
82 :
83 42 : void SphericalHarmonic::read( ActionWithArguments* action ) {
84 42 : parse(action,"L",tmom);
85 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() );
86 42 : if( action->getNumberOfArguments()==4 ) {
87 42 : action->log.printf(" multiplying cylindrical harmonic by weight from %s \n", action->getPntrToArgument(3)->getName().c_str() );
88 : }
89 :
90 42 : normaliz.resize( tmom+1 );
91 232 : for(unsigned i=0; i<=tmom; ++i) {
92 190 : normaliz[i] = sqrt( (2*tmom+1)*factorial(tmom-i)/(4*pi*factorial(tmom+i)) );
93 190 : if( i%2==1 ) {
94 84 : normaliz[i]*=-1;
95 : }
96 : }
97 :
98 42 : coeff_poly.resize( tmom+1 );
99 42 : if( tmom==1 ) {
100 : // Legendre polynomial coefficients of order one
101 19 : coeff_poly[0]=0;
102 19 : coeff_poly[1]=1.0;
103 : } else if( tmom==2 ) {
104 : // Legendre polynomial coefficients of order two
105 0 : coeff_poly[0]=-0.5;
106 0 : coeff_poly[1]=0.0;
107 0 : coeff_poly[2]=1.5;
108 : } else if( tmom==3 ) {
109 : // Legendre polynomial coefficients of order three
110 1 : coeff_poly[0]=0.0;
111 1 : coeff_poly[1]=-1.5;
112 1 : coeff_poly[2]=0.0;
113 1 : coeff_poly[3]=2.5;
114 : } else if( tmom==4 ) {
115 : // Legendre polynomial coefficients of order four
116 3 : coeff_poly[0]=0.375;
117 3 : coeff_poly[1]=0.0;
118 3 : coeff_poly[2]=-3.75;
119 3 : coeff_poly[3]=0.0;
120 3 : coeff_poly[4]=4.375;
121 : } else if( tmom==5 ) {
122 : // Legendre polynomial coefficients of order five
123 0 : coeff_poly[0]=0.0;
124 0 : coeff_poly[1]=1.875;
125 0 : coeff_poly[2]=0.0;
126 0 : coeff_poly[3]=-8.75;
127 0 : coeff_poly[4]=0.0;
128 0 : coeff_poly[5]=7.875;
129 : } else if( tmom==6 ) {
130 : // Legendre polynomial coefficients of order six
131 19 : coeff_poly[0]=-0.3125;
132 19 : coeff_poly[1]=0.0;
133 19 : coeff_poly[2]=6.5625;
134 19 : coeff_poly[3]=0.0;
135 19 : coeff_poly[4]=-19.6875;
136 19 : coeff_poly[5]=0.0;
137 19 : coeff_poly[6]=14.4375;
138 : } else {
139 0 : action->error("Insert Legendre polynomial coefficients into SphericalHarmonics code");
140 : }
141 42 : }
142 :
143 84 : std::vector<std::string> SphericalHarmonic::getComponentsPerLabel() const {
144 : std::vector<std::string> comp;
145 : std::string num;
146 760 : for(int i=-tmom; i<=tmom; ++i) {
147 676 : Tools::convert(fabs(i),num);
148 676 : if( i<0 ) {
149 592 : comp.push_back( "-n" + num );
150 380 : } else if( i>0 ) {
151 592 : comp.push_back( "-p" + num );
152 : } else {
153 168 : comp.push_back( "-0");
154 : }
155 : }
156 84 : return comp;
157 0 : }
158 :
159 42 : void SphericalHarmonic::setPeriodicityForOutputs( ActionWithValue* action ) {
160 42 : std::vector<std::string> comp( getComponentsPerLabel() );
161 380 : for(unsigned i=0; i<comp.size(); ++i) {
162 676 : action->componentIsNotPeriodic("rm" + comp[i]);
163 676 : action->componentIsNotPeriodic("im" + comp[i]);
164 : }
165 42 : }
166 :
167 1742755 : void SphericalHarmonic::calc( const ActionWithArguments* action, const std::vector<double>& args, std::vector<double>& vals, Matrix<double>& derivatives ) const {
168 : double weight=1;
169 1742755 : if( args.size()==4 ) {
170 1742755 : weight = args[3];
171 : }
172 1742755 : if( weight<epsilon ) {
173 1468550 : return;
174 : }
175 :
176 274205 : double dlen2 = args[0]*args[0]+args[1]*args[1]+args[2]*args[2];
177 274205 : double dlen = sqrt( dlen2 );
178 274205 : double dlen3 = dlen2*dlen;
179 274205 : double dpoly_ass, poly_ass=deriv_poly( 0, args[2]/dlen, dpoly_ass );
180 : // Derivatives of z/r wrt x, y, z
181 274205 : Vector dz;
182 274205 : dz[0] = -( args[2] / dlen3 )*args[0];
183 274205 : dz[1] = -( args[2] / dlen3 )*args[1];
184 274205 : dz[2] = -( args[2] / dlen3 )*args[2] + (1.0 / dlen);
185 : // Accumulate for m=0
186 274205 : vals[tmom] = weight*poly_ass;
187 274205 : addVectorDerivatives( tmom, weight*dpoly_ass*dz, derivatives );
188 274205 : if( args.size()==4 ) {
189 274205 : derivatives(tmom, 3) = poly_ass;
190 : }
191 :
192 : // The complex number of which we have to take powers
193 274205 : std::complex<double> com1( args[0]/dlen,args[1]/dlen ), dp_x, dp_y, dp_z;
194 : std::complex<double> powered = std::complex<double>(1.0,0.0);
195 : std::complex<double> ii( 0.0, 1.0 );
196 274205 : Vector myrealvec, myimagvec, real_dz, imag_dz;
197 : // Do stuff for all other m values
198 1794050 : for(unsigned m=1; m<=tmom; ++m) {
199 : // Calculate Legendre Polynomial
200 1519845 : poly_ass=deriv_poly( m, args[2]/dlen, dpoly_ass );
201 : // Real and imaginary parts of z
202 : double real_z = real(com1*powered), imag_z = imag(com1*powered);
203 :
204 : // Calculate steinhardt parameter
205 1519845 : double tq6=poly_ass*real_z; // Real part of steinhardt parameter
206 1519845 : double itq6=poly_ass*imag_z; // Imaginary part of steinhardt parameter
207 :
208 : // Derivatives wrt ( x/r + iy )^m
209 1519845 : double md=static_cast<double>(m);
210 1519845 : dp_x = md*powered*( (1.0/dlen)-(args[0]*args[0])/dlen3-ii*(args[0]*args[1])/dlen3 );
211 1519845 : dp_y = md*powered*( ii*(1.0/dlen)-(args[0]*args[1])/dlen3-ii*(args[1]*args[1])/dlen3 );
212 1519845 : dp_z = md*powered*( -(args[0]*args[2])/dlen3-ii*(args[1]*args[2])/dlen3 );
213 :
214 : // Derivatives of real and imaginary parts of above
215 1519845 : real_dz[0] = real( dp_x );
216 1519845 : real_dz[1] = real( dp_y );
217 1519845 : real_dz[2] = real( dp_z );
218 1519845 : imag_dz[0] = imag( dp_x );
219 1519845 : imag_dz[1] = imag( dp_y );
220 1519845 : imag_dz[2] = imag( dp_z );
221 :
222 : // Complete derivative of steinhardt parameter
223 1519845 : myrealvec = weight*dpoly_ass*real_z*dz + weight*poly_ass*real_dz;
224 1519845 : myimagvec = weight*dpoly_ass*imag_z*dz + weight*poly_ass*imag_dz;
225 :
226 : // Real part
227 1519845 : vals[tmom+m] = weight*tq6;
228 1519845 : addVectorDerivatives( tmom+m, myrealvec, derivatives );
229 : // Imaginary part
230 1519845 : vals[3*tmom+1+m] = weight*itq6;
231 1519845 : addVectorDerivatives( 3*tmom+1+m, myimagvec, derivatives );
232 : // Store -m part of vector
233 1519845 : double pref=pow(-1.0,m);
234 : // -m part of vector is just +m part multiplied by (-1.0)**m and multiplied by complex
235 : // conjugate of Legendre polynomial
236 : // Real part
237 1519845 : vals[tmom-m] = pref*weight*tq6;
238 1519845 : addVectorDerivatives( tmom-m, pref*myrealvec, derivatives );
239 : // Imaginary part
240 1519845 : vals[3*tmom+1-m] = -pref*weight*itq6;
241 1519845 : addVectorDerivatives( 3*tmom+1-m, -pref*myimagvec, derivatives );
242 1519845 : if( args.size()==4 ) {
243 1519845 : derivatives(tmom+m,3)=tq6;
244 1519845 : derivatives(3*tmom+1+m, 3)=itq6;
245 1519845 : derivatives(tmom-m,3)=pref*tq6;
246 1519845 : derivatives(3*tmom+1-m, 3)=-pref*itq6;
247 : }
248 : // Calculate next power of complex number
249 : powered *= com1;
250 : }
251 : }
252 :
253 1794050 : double SphericalHarmonic::deriv_poly( const unsigned& m, const double& val, double& df ) const {
254 : double fact=1.0;
255 7004030 : for(unsigned j=1; j<=m; ++j) {
256 5209980 : fact=fact*j;
257 : }
258 1794050 : double res=coeff_poly[m]*fact;
259 :
260 1794050 : double pow=1.0, xi=val, dxi=1.0;
261 1794050 : df=0.0;
262 7004030 : for(int i=m+1; i<=tmom; ++i) {
263 : fact=1.0;
264 13759570 : for(unsigned j=i-m+1; j<=i; ++j) {
265 8549590 : fact=fact*j;
266 : }
267 5209980 : res=res+coeff_poly[i]*fact*xi;
268 5209980 : df = df + pow*coeff_poly[i]*fact*dxi;
269 5209980 : xi=xi*val;
270 5209980 : dxi=dxi*val;
271 5209980 : pow+=1.0;
272 : }
273 1794050 : df = df*normaliz[m];
274 1794050 : return normaliz[m]*res;
275 : }
276 :
277 6353585 : void SphericalHarmonic::addVectorDerivatives( const unsigned& ival, const Vector& der, Matrix<double>& derivatives ) const {
278 25414340 : for(unsigned j=0; j<3; ++j) {
279 19060755 : derivatives(ival,j) = der[j];
280 : }
281 6353585 : }
282 :
283 : }
284 : }
285 :
|
