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_ELECTRONIC_COLUMNBUNDLEOPERATORS_INTERNAL_H
21 #define JDFTX_ELECTRONIC_COLUMNBUNDLEOPERATORS_INTERNAL_H
22
23 #include <core/matrix3.h>
24
25 //! @cond
26
reducedL_calc(int j,int nbasis,int ncols,const complex * Y,complex * LY,const matrix3<> GGT,const vector3<int> * iGarr,const vector3<> & k,double detR)27 __hostanddev__ void reducedL_calc(int j, int nbasis, int ncols, const complex* Y, complex* LY,
28 const matrix3<> GGT, const vector3<int>* iGarr, const vector3<>& k, double detR)
29 { for (int i=0; i < ncols; i++)
30 LY[nbasis*i+j] = (-detR * GGT.metric_length_squared(iGarr[j]+k)) * Y[nbasis*i+j];
31 }
32
reducedLinv_calc(int j,int nbasis,int ncols,const complex * Y,complex * LinvY,const matrix3<> GGT,const vector3<int> * iGarr,const vector3<> & k,double detR)33 __hostanddev__ void reducedLinv_calc(int j, int nbasis, int ncols, const complex* Y, complex* LinvY,
34 const matrix3<> GGT, const vector3<int>* iGarr, const vector3<>& k, double detR)
35 { for (int i=0; i < ncols; i++)
36 { double G2 = GGT.metric_length_squared(iGarr[j]+k);
37 LinvY[nbasis*i+j] = (G2 ? -1./(detR*G2) : 0.) * Y[nbasis*i+j];
38 }
39 }
40
41 //Kernel for lattice derivative of Tr[Y^LYF] calculation
reducedLstress_calc(int j,int nbasis,int ncols,const complex * Y,const double * F,const vector3<int> * iGarr,const vector3<> & k,symmetricMatrix3<> * result)42 __hostanddev__ void reducedLstress_calc(int j, int nbasis, int ncols, const complex* Y, const double* F,
43 const vector3<int>* iGarr, const vector3<>& k, symmetricMatrix3<>* result)
44 { //Sum wave function norm at each iG over bands / columns:
45 double bandSum = 0.;
46 for (int i=0; i<ncols; i++)
47 bandSum += F[i] * Y[nbasis*i+j].norm();
48 //Collect weighted iG outer product (subsequently reduced over iG):
49 vector3<> iGk = iGarr[j]+k;
50 result[j] = outer(iGk) * bandSum;
51 }
52
53
precond_inv_kinetic_calc(int j,int nbasis,int ncols,complex * Ydata,double KErollover,const matrix3<> GGT,const vector3<int> * iGarr,const vector3<> k,double invdetR)54 __hostanddev__ void precond_inv_kinetic_calc(int j, int nbasis, int ncols, complex* Ydata,
55 double KErollover, const matrix3<> GGT, const vector3<int>* iGarr, const vector3<> k, double invdetR)
56 {
57 double x = 0.5*GGT.metric_length_squared(iGarr[j]+k)/KErollover;
58 double precondFactor = 1.+x*(1.+x*(1.+x*(1.+x*(1.+x*(1.+x*(1.+x*(1.+x)))))));
59 precondFactor = precondFactor*invdetR/(1.+x*precondFactor);
60 for(int i=0; i < ncols; i++)
61 Ydata[nbasis*i+j] *= precondFactor;
62 }
63
precond_inv_kinetic_band_calc(int j,int nbasis,int ncols,complex * Ydata,const double * KEref,const matrix3<> & GGT,const vector3<int> * iGarr,const vector3<> & k)64 __hostanddev__ void precond_inv_kinetic_band_calc(int j, int nbasis, int ncols, complex* Ydata, const double* KEref,
65 const matrix3<>& GGT, const vector3<int>* iGarr, const vector3<>& k)
66 {
67 double KE = 0.5*GGT.metric_length_squared(iGarr[j]+k);
68 for(int i=0; i<ncols; i++)
69 { double x = KE/KEref[i];
70 Ydata[nbasis*i+j] *= (27.+x*(18.+x*(12.+x*8.))) / (27.+x*(18.+x*(12.+x*(8.+x*16))));
71 }
72 }
73
74 __hostanddev__
translate_calc(int j,int nbasis,int ncols,complex * Y,const vector3<int> * iGarr,const vector3<> & k,const vector3<> & dr)75 void translate_calc(int j, int nbasis, int ncols, complex* Y, const vector3<int>* iGarr, const vector3<>& k, const vector3<>& dr)
76 { complex tFactor = cis(-2*M_PI*dot(iGarr[j]+k,dr));
77 for(int i=0; i<ncols; i++)
78 Y[nbasis*i+j] *= tFactor;
79 }
80
81 __hostanddev__
translateColumns_calc(int j,int nbasis,int ncols,complex * Y,const vector3<int> * iGarr,const vector3<> & k,const vector3<> * dr)82 void translateColumns_calc(int j, int nbasis, int ncols, complex* Y, const vector3<int>* iGarr, const vector3<>& k, const vector3<>* dr)
83 { for(int i=0; i<ncols; i++)
84 Y[nbasis*i+j] *= cis(-2*M_PI*dot(iGarr[j]+k,dr[i]));
85 }
86
reducedD_calc(int j,int nbasis,int ncols,const complex * Ydata,complex * DYdata,const vector3<int> * iGarr,double kdotGe,const vector3<> Ge)87 __hostanddev__ void reducedD_calc(int j, int nbasis, int ncols, const complex* Ydata, complex* DYdata,
88 const vector3<int>* iGarr, double kdotGe, const vector3<> Ge)
89 {
90 complex Di(0, kdotGe+dot(iGarr[j],Ge)); // i (k+G).e where e is the cartesian direction of interest
91 for(int i=0; i < ncols; i++)
92 DYdata[nbasis*i+j] = Di*Ydata[nbasis*i+j];
93 }
94
reducedDD_calc(int j,int nbasis,int ncols,const complex * Ydata,complex * DYdata,const vector3<int> * iGarr,double kdotGe1,double kdotGe2,const vector3<> Ge1,const vector3<> Ge2)95 __hostanddev__ void reducedDD_calc(int j, int nbasis, int ncols, const complex* Ydata, complex* DYdata,
96 const vector3<int>* iGarr, double kdotGe1, double kdotGe2, const vector3<> Ge1, const vector3<> Ge2)
97 {
98 complex Di(0, kdotGe1+dot(iGarr[j],Ge1)); // i (k+G).ei
99 complex Dj(0, kdotGe2+dot(iGarr[j],Ge2)); // i (k+G).ej
100 for(int i=0; i < ncols; i++)
101 DYdata[nbasis*i+j] = Di*Dj*Ydata[nbasis*i+j];
102 }
103
104 //! @endcond
105 #endif // JDFTX_ELECTRONIC_COLUMNBUNDLEOPERATORS_INTERNAL_H
106