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","matrix","the real parts of the spherical harmonic values with the m value given");
76 430 : keys.addOutputComponent("im","default","matrix","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 ) action->log.printf(" multiplying cylindrical harmonic by weight from %s \n", action->getPntrToArgument(3)->getName().c_str() );
87 :
88 42 : normaliz.resize( tmom+1 );
89 232 : for(unsigned i=0; i<=tmom; ++i) {
90 190 : normaliz[i] = sqrt( (2*tmom+1)*factorial(tmom-i)/(4*pi*factorial(tmom+i)) );
91 190 : if( i%2==1 ) normaliz[i]*=-1;
92 : }
93 :
94 42 : coeff_poly.resize( tmom+1 );
95 42 : if( tmom==1 ) {
96 : // Legendre polynomial coefficients of order one
97 19 : coeff_poly[0]=0; coeff_poly[1]=1.0;
98 : } else if( tmom==2 ) {
99 : // Legendre polynomial coefficients of order two
100 0 : coeff_poly[0]=-0.5; coeff_poly[1]=0.0;
101 0 : coeff_poly[2]=1.5;
102 : } else if( tmom==3 ) {
103 : // Legendre polynomial coefficients of order three
104 1 : coeff_poly[0]=0.0; coeff_poly[1]=-1.5;
105 1 : coeff_poly[2]=0.0; coeff_poly[3]=2.5;
106 : } else if( tmom==4 ) {
107 : // Legendre polynomial coefficients of order four
108 3 : coeff_poly[0]=0.375; coeff_poly[1]=0.0;
109 3 : coeff_poly[2]=-3.75; coeff_poly[3]=0.0;
110 3 : coeff_poly[4]=4.375;
111 : } else if( tmom==5 ) {
112 : // Legendre polynomial coefficients of order five
113 0 : coeff_poly[0]=0.0; coeff_poly[1]=1.875;
114 0 : coeff_poly[2]=0.0; coeff_poly[3]=-8.75;
115 0 : coeff_poly[4]=0.0; coeff_poly[5]=7.875;
116 : } else if( tmom==6 ) {
117 : // Legendre polynomial coefficients of order six
118 19 : coeff_poly[0]=-0.3125; coeff_poly[1]=0.0;
119 19 : coeff_poly[2]=6.5625; coeff_poly[3]=0.0;
120 19 : coeff_poly[4]=-19.6875; coeff_poly[5]=0.0;
121 19 : coeff_poly[6]=14.4375;
122 : } else {
123 0 : action->error("Insert Legendre polynomial coefficients into SphericalHarmonics code");
124 : }
125 42 : }
126 :
127 84 : std::vector<std::string> SphericalHarmonic::getComponentsPerLabel() const {
128 : std::vector<std::string> comp; std::string num;
129 760 : for(int i=-tmom; i<=tmom; ++i) {
130 676 : Tools::convert(fabs(i),num);
131 972 : if( i<0 ) comp.push_back( "-n" + num );
132 676 : else if( i>0 ) comp.push_back( "-p" + num );
133 168 : else comp.push_back( "-0");
134 : }
135 84 : return comp;
136 0 : }
137 :
138 42 : void SphericalHarmonic::setPeriodicityForOutputs( ActionWithValue* action ) {
139 42 : std::vector<std::string> comp( getComponentsPerLabel() );
140 718 : for(unsigned i=0; i<comp.size(); ++i) { action->componentIsNotPeriodic("rm" + comp[i]); action->componentIsNotPeriodic("im" + comp[i]); }
141 42 : }
142 :
143 1742755 : void SphericalHarmonic::calc( const ActionWithArguments* action, const std::vector<double>& args, std::vector<double>& vals, Matrix<double>& derivatives ) const {
144 1742755 : double weight=1; if( args.size()==4 ) weight = args[3];
145 1742755 : if( weight<epsilon ) return;
146 :
147 274205 : double dlen2 = args[0]*args[0]+args[1]*args[1]+args[2]*args[2];
148 274205 : double dlen = sqrt( dlen2 ); double dlen3 = dlen2*dlen;
149 274205 : double dpoly_ass, poly_ass=deriv_poly( 0, args[2]/dlen, dpoly_ass );
150 : // Derivatives of z/r wrt x, y, z
151 274205 : Vector dz;
152 274205 : dz[0] = -( args[2] / dlen3 )*args[0];
153 274205 : dz[1] = -( args[2] / dlen3 )*args[1];
154 274205 : dz[2] = -( args[2] / dlen3 )*args[2] + (1.0 / dlen);
155 : // Accumulate for m=0
156 274205 : vals[tmom] = weight*poly_ass;
157 274205 : addVectorDerivatives( tmom, weight*dpoly_ass*dz, derivatives );
158 274205 : if( args.size()==4 ) derivatives(tmom, 3) = poly_ass;
159 :
160 : // The complex number of which we have to take powers
161 274205 : std::complex<double> com1( args[0]/dlen,args[1]/dlen ), dp_x, dp_y, dp_z;
162 : std::complex<double> powered = std::complex<double>(1.0,0.0); std::complex<double> ii( 0.0, 1.0 );
163 274205 : Vector myrealvec, myimagvec, real_dz, imag_dz;
164 : // Do stuff for all other m values
165 1794050 : for(unsigned m=1; m<=tmom; ++m) {
166 : // Calculate Legendre Polynomial
167 1519845 : poly_ass=deriv_poly( m, args[2]/dlen, dpoly_ass );
168 : // Real and imaginary parts of z
169 : double real_z = real(com1*powered), imag_z = imag(com1*powered);
170 :
171 : // Calculate steinhardt parameter
172 1519845 : double tq6=poly_ass*real_z; // Real part of steinhardt parameter
173 1519845 : double itq6=poly_ass*imag_z; // Imaginary part of steinhardt parameter
174 :
175 : // Derivatives wrt ( x/r + iy )^m
176 1519845 : double md=static_cast<double>(m);
177 1519845 : dp_x = md*powered*( (1.0/dlen)-(args[0]*args[0])/dlen3-ii*(args[0]*args[1])/dlen3 );
178 1519845 : dp_y = md*powered*( ii*(1.0/dlen)-(args[0]*args[1])/dlen3-ii*(args[1]*args[1])/dlen3 );
179 1519845 : dp_z = md*powered*( -(args[0]*args[2])/dlen3-ii*(args[1]*args[2])/dlen3 );
180 :
181 : // Derivatives of real and imaginary parts of above
182 1519845 : real_dz[0] = real( dp_x ); real_dz[1] = real( dp_y ); real_dz[2] = real( dp_z );
183 1519845 : imag_dz[0] = imag( dp_x ); imag_dz[1] = imag( dp_y ); imag_dz[2] = imag( dp_z );
184 :
185 : // Complete derivative of steinhardt parameter
186 1519845 : myrealvec = weight*dpoly_ass*real_z*dz + weight*poly_ass*real_dz;
187 1519845 : myimagvec = weight*dpoly_ass*imag_z*dz + weight*poly_ass*imag_dz;
188 :
189 : // Real part
190 1519845 : vals[tmom+m] = weight*tq6;
191 1519845 : addVectorDerivatives( tmom+m, myrealvec, derivatives );
192 : // Imaginary part
193 1519845 : vals[3*tmom+1+m] = weight*itq6;
194 1519845 : addVectorDerivatives( 3*tmom+1+m, myimagvec, derivatives );
195 : // Store -m part of vector
196 1519845 : double pref=pow(-1.0,m);
197 : // -m part of vector is just +m part multiplied by (-1.0)**m and multiplied by complex
198 : // conjugate of Legendre polynomial
199 : // Real part
200 1519845 : vals[tmom-m] = pref*weight*tq6;
201 1519845 : addVectorDerivatives( tmom-m, pref*myrealvec, derivatives );
202 : // Imaginary part
203 1519845 : vals[3*tmom+1-m] = -pref*weight*itq6;
204 1519845 : addVectorDerivatives( 3*tmom+1-m, -pref*myimagvec, derivatives );
205 1519845 : if( args.size()==4 ) {
206 1519845 : derivatives(tmom+m,3)=tq6; derivatives(3*tmom+1+m, 3)=itq6;
207 1519845 : derivatives(tmom-m,3)=pref*tq6; derivatives(3*tmom+1-m, 3)=-pref*itq6;
208 : }
209 : // Calculate next power of complex number
210 : powered *= com1;
211 : }
212 : }
213 :
214 1794050 : double SphericalHarmonic::deriv_poly( const unsigned& m, const double& val, double& df ) const {
215 : double fact=1.0;
216 7004030 : for(unsigned j=1; j<=m; ++j) fact=fact*j;
217 1794050 : double res=coeff_poly[m]*fact;
218 :
219 1794050 : double pow=1.0, xi=val, dxi=1.0; df=0.0;
220 7004030 : for(int i=m+1; i<=tmom; ++i) {
221 : fact=1.0;
222 13759570 : for(unsigned j=i-m+1; j<=i; ++j) fact=fact*j;
223 5209980 : res=res+coeff_poly[i]*fact*xi;
224 5209980 : df = df + pow*coeff_poly[i]*fact*dxi;
225 5209980 : xi=xi*val; dxi=dxi*val; pow+=1.0;
226 : }
227 1794050 : df = df*normaliz[m];
228 1794050 : return normaliz[m]*res;
229 : }
230 :
231 6353585 : void SphericalHarmonic::addVectorDerivatives( const unsigned& ival, const Vector& der, Matrix<double>& derivatives ) const {
232 25414340 : for(unsigned j=0; j<3; ++j) derivatives(ival,j) = der[j];
233 6353585 : }
234 :
235 : }
236 : }
237 :
|
