1 /*-------------------------------------------------------------------
2 Copyright 2011 Ravishankar Sundararaman
3
4 This file is part of JDFTx.
5
6 JDFTx is free software: you can redistribute it and/or modify
7 it under the terms of the GNU General Public License as published by
8 the Free Software Foundation, either version 3 of the License, or
9 (at your option) any later version.
10
11 JDFTx is distributed in the hope that it will be useful,
12 but WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 GNU General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with JDFTx. If not, see <http://www.gnu.org/licenses/>.
18 -------------------------------------------------------------------*/
19
20 #ifndef JDFTX_CORE_OPERATORS_INTERNAL_H
21 #define JDFTX_CORE_OPERATORS_INTERNAL_H
22
23 #include <core/matrix3.h>
24 #include <core/tensor3.h>
25 #include <core/LoopMacros.h>
26 #include <core/SphericalHarmonics.h>
27 #include <core/RadialFunction.h>
28
29 //! @cond
30
31 //! Compute the full G-space indices corresponding to a half Gspace index and its mirror image:
32 #define COMPUTE_fullRefIndices \
33 vector3<int> iv(iG), ivRef; \
34 for(int k=0; k<3; k++) \
35 { if(iv[k]<0) iv[k] += S[k]; \
36 ivRef[k] = iv[k] ? S[k] - iv[k] : 0; \
37 } \
38 int iFull = iv[2] + S[2]*(iv[1] + S[1]*iv[0]); \
39 int iFullRef = ivRef[2] + S[2]*(ivRef[1] + S[1]*ivRef[0]);
40
41
42 __hostanddev__
RealG_calc(int iHalf,const vector3<int> iG,const vector3<int> S,const complex * vFull,complex * vHalf,double scaleFac)43 void RealG_calc(int iHalf, const vector3<int> iG, const vector3<int> S,
44 const complex* vFull, complex* vHalf, double scaleFac)
45 {
46 COMPUTE_fullRefIndices
47 vHalf[iHalf] = scaleFac*0.5*(vFull[iFull] + vFull[iFullRef].conj());
48 }
49
50 __hostanddev__
ImagG_calc(int iHalf,const vector3<int> iG,const vector3<int> S,const complex * vFull,complex * vHalf,double scaleFac)51 void ImagG_calc(int iHalf, const vector3<int> iG, const vector3<int> S,
52 const complex* vFull, complex* vHalf, double scaleFac)
53 {
54 COMPUTE_fullRefIndices
55 vHalf[iHalf] = complex(0,-scaleFac*0.5)*(vFull[iFull] - vFull[iFullRef].conj());
56 }
57
58 __hostanddev__
ComplexG_calc(int iHalf,const vector3<int> iG,const vector3<int> S,const complex * vHalf,complex * vFull,double scaleFac)59 void ComplexG_calc(int iHalf, const vector3<int> iG, const vector3<int> S,
60 const complex* vHalf, complex *vFull, double scaleFac)
61 {
62 COMPUTE_fullRefIndices
63 complex temp = scaleFac*vHalf[iHalf];
64 vFull[iFull] = temp; //Copy value into corresponding location in full-space
65 vFull[iFullRef] = temp.conj(); //Also set complex conjugate into the mirror location
66 }
67
68 __hostanddev__
changeGrid_calc(const vector3<int> & iG,const vector3<int> & Sin,const vector3<int> & Sout,const complex * in,complex * out)69 void changeGrid_calc(const vector3<int>& iG, const vector3<int>& Sin, const vector3<int>& Sout, const complex* in, complex* out)
70 { //Compute index:
71 #define COMPUTE_index(suffix) \
72 int i##suffix = 0; \
73 for(int k=0; k<2; k++) \
74 { if(2*iG[k]<1-S##suffix[k] || 2*iG[k]>S##suffix[k]) return; \
75 i##suffix = i##suffix * S##suffix[k] + (iG[k]<0 ? (iG[k]+S##suffix[k]) : iG[k]); \
76 } \
77 if(2*iG[2]>S##suffix[2]) return; \
78 else i##suffix = i##suffix*(1+S##suffix[2]/2) + iG[2];
79 COMPUTE_index(in)
80 COMPUTE_index(out)
81 #undef COMPUTE_index
82 out[iout] = in[iin];
83 }
84 __hostanddev__
changeGridFull_calc(const vector3<int> & iG,const vector3<int> & Sin,const vector3<int> & Sout,const complex * in,complex * out)85 void changeGridFull_calc(const vector3<int>& iG, const vector3<int>& Sin, const vector3<int>& Sout, const complex* in, complex* out)
86 { //Compute index:
87 #define COMPUTE_index(suffix) \
88 int i##suffix = 0; \
89 for(int k=0; k<3; k++) \
90 { if(2*iG[k]<1-S##suffix[k] || 2*iG[k]>S##suffix[k]) return; \
91 i##suffix = i##suffix * S##suffix[k] + (iG[k]<0 ? (iG[k]+S##suffix[k]) : iG[k]); \
92 }
93 COMPUTE_index(in)
94 COMPUTE_index(out)
95 #undef COMPUTE_index
96 out[iout] = in[iin];
97 }
98
D_calc(int i,const vector3<int> & iG,const complex * in,complex * out,const vector3<> & Ge)99 __hostanddev__ void D_calc(int i, const vector3<int>& iG, const complex* in, complex* out,
100 const vector3<>& Ge)
101 { out[i] = in[i] * complex(0, dot(iG,Ge));
102 }
DD_calc(int i,const vector3<int> & iG,const complex * in,complex * out,const vector3<> & Ge1,const vector3<> & Ge2)103 __hostanddev__ void DD_calc(int i, const vector3<int>& iG, const complex* in, complex* out,
104 const vector3<>& Ge1, const vector3<>& Ge2)
105 { out[i] = in[i] * (-dot(iG,Ge1) * dot(iG,Ge2));
106 }
107
108
109 //! Switch a function templated over l for all supported l with parenthesis enclosed argument list argList
110 #define SwitchTemplate_l(l,fTemplate,argList) \
111 switch(l) \
112 { case 0: fTemplate<0> argList; break; \
113 case 1: fTemplate<1> argList; break; \
114 case 2: fTemplate<2> argList; break; \
115 case 3: fTemplate<3> argList; break; \
116 case 4: fTemplate<4> argList; break; \
117 case 5: fTemplate<5> argList; break; \
118 case 6: fTemplate<6> argList; break; \
119 }
120
121 template<int l, int lpm> struct lGradient_staticLoop
setlGradient_staticLoop122 { static __hostanddev__ void set(int i, const vector3<>& g, const complex& phasedIn, const array<complex*,2*l+1>& out)
123 { out[lpm][i] = phasedIn * Ylm<l,lpm-l>(g);
124 lGradient_staticLoop<l,lpm-1>::set(i, g, phasedIn, out);
125 }
126 };
127 template<int l> struct lGradient_staticLoop<l,-1> { static __hostanddev__ void set(int i, const vector3<>& g, const complex& phasedIn, const array<complex*,2*l+1>& out) {} }; //end recursion
128
129 template<int l> __hostanddev__ void lGradient_calc(int i, const vector3<int>& iG, bool isNyq, const complex* in, const array<complex*,2*l+1>& out, const matrix3<>& G)
130 { const complex phase = cis(l*0.5*M_PI); // iota^l (computable at compile time)
131 lGradient_staticLoop<l,l+l>::set(i, iG*G, isNyq ? 0. : phase * in[i], out);
132 }
133
134 template<int l, int lpm> struct lDivergence_staticLoop
135 { static __hostanddev__ complex get(int i, const vector3<>& g, const array<const complex*,2*l+1>& in)
136 { return in[lpm][i] * Ylm<l,lpm-l>(g) + lDivergence_staticLoop<l,lpm-1>::get(i, g, in);
137 }
138 };
139 template<int l> struct lDivergence_staticLoop<l,-1> { static __hostanddev__ complex get(int i, const vector3<>& g, const array<const complex*,2*l+1>& in) { return complex(); } }; //end recursion
140
141 template<int l> __hostanddev__ void lDivergence_calc(int i, const vector3<int>& iG, bool isNyq, const array<const complex*,2*l+1>& in, complex* out, const matrix3<>& G)
142 { const complex phase = cis(l*0.5*M_PI); // iota^l (computable at compile time)
143 out[i] = (isNyq ? 0. : phase) * lDivergence_staticLoop<l,l+l>::get(i, iG*G, in);
144 }
145
146
147 template<int lm> __hostanddev__ symmetricMatrix3<> lGradientStress_calc(const vector3<int>& iG, const matrix3<>& G, const RadialFunctionG& w)
148 { vector3<> qvec = iG * G; //k+G in cartesian coordinates
149 double q = qvec.length();
150 vector3<> qhat = qvec * (q ? 1.0/q : 0.0); //unit vector || qvec (set to 0 for q=0 (doesn't matter))
151 //q-gradient of w * Y: (Note that iota^l phase handled outside)
152 vector3<> wYprime = w(q) * YlmPrime<lm>(qvec) + qhat * (w.deriv(q) * Ylm<lm>(qvec));
153 return symmetricMatrix3<>(qvec[0]*wYprime[0], qvec[1]*wYprime[1], qvec[2]*wYprime[2],
154 qvec[1]*wYprime[2], qvec[2]*wYprime[0], qvec[0]*wYprime[1]);
155 }
156
157 __hostanddev__ complex blochPhase_calc(const vector3<int>& iv, const vector3<>& invS, const vector3<>& k)
158 { return cis(2*M_PI*dot(k, vector3<>(iv[0]*invS[0], iv[1]*invS[1], iv[2]*invS[2])));
159 }
160
161 __hostanddev__ complex radialFunction_calc(const vector3<int>& iG, const matrix3<>& GGT, const RadialFunctionG& f, const vector3<>& r0)
162 { return f(sqrt(GGT.metric_length_squared(iG))) * cis(-2*M_PI*dot(iG,r0));
163 }
164
165
166 //----------- Functions from VectorField.h -------
167
168 __hostanddev__ void gradient_calc(int i, const vector3<int> iG, bool nyq, const matrix3<> G,
169 const complex* Xtilde, vector3<complex*>& gradTilde)
170 { complex iota(0.0, nyq ? 0.0 : 1.0); //zero nyquist frequencies
171 storeVector((iG*G) * (iota*Xtilde[i]), gradTilde, i);
172 }
173
174 __hostanddev__ void divergence_calc(int i, const vector3<int> iG, bool nyq, const matrix3<> G,
175 vector3<const complex*>& Vtilde, complex* divTilde)
176 { complex iota(0.0, nyq ? 0.0 : 1.0); //zero nyquist frequencies
177 divTilde[i] = iota * dot(iG*G, loadVector(Vtilde,i));
178 }
179
180
181 __hostanddev__ void tensorGradient_calc(int i, const vector3<int> iG, bool nyq, const matrix3<> G,
182 const complex* Xtilde, tensor3<complex*>& gradTilde)
183 {
184 complex minus_Xtilde = nyq ? complex(0,0) : -Xtilde[i]; //zero nyquist frequencies
185 vector3<> Gvec = iG*G;
186 double Gsq = Gvec.length_squared();
187 gradTilde.xy()[i] = minus_Xtilde*Gvec.x()*Gvec.y();
188 gradTilde.yz()[i] = minus_Xtilde*Gvec.y()*Gvec.z();
189 gradTilde.zx()[i] = minus_Xtilde*Gvec.z()*Gvec.x();
190 gradTilde.xxr()[i] = minus_Xtilde*(Gvec.x()*Gvec.x() - (1.0/3)*Gsq);
191 gradTilde.yyr()[i] = minus_Xtilde*(Gvec.y()*Gvec.y() - (1.0/3)*Gsq);
192 }
193
194 __hostanddev__ void tensorDivergence_calc(int i, const vector3<int> iG, bool nyq, const matrix3<> G,
195 tensor3<const complex*>& Vtilde, complex* divTilde)
196 {
197 complex temp = complex(0,0);
198 if(!nyq)
199 { vector3<> Gvec = iG*G;
200 temp += Vtilde.xy()[i]*( 2*Gvec.x()*Gvec.y() );
201 temp += Vtilde.yz()[i]*( 2*Gvec.y()*Gvec.z() );
202 temp += Vtilde.zx()[i]*( 2*Gvec.z()*Gvec.x() );
203 temp += Vtilde.xxr()[i]*( Gvec.x()*Gvec.x() - Gvec.z()*Gvec.z() );
204 temp += Vtilde.yyr()[i]*( Gvec.y()*Gvec.y() - Gvec.z()*Gvec.z() );
205 }
206 divTilde[i] = -temp;
207 }
208
209 //! @endcond
210 #endif //JDFTX_CORE_OPERATORS_INTERNAL_H
211