From c6b901030d6fbdf22a8e13daa534f8efe42d8439 Mon Sep 17 00:00:00 2001 From: zdy Date: Thu, 22 Apr 2021 22:30:42 +0800 Subject: [PATCH] update soc class for non-local pseudopotential initial --- ABACUS.develop/source/src_lcao/ORB_read.cpp | 14 ++-- .../source/src_pw/pseudopot_cell_vnl.cpp | 14 ++-- ABACUS.develop/source/src_pw/soc.cpp | 77 ++++++++----------- ABACUS.develop/source/src_pw/soc.h | 26 +++++-- 4 files changed, 63 insertions(+), 68 deletions(-) diff --git a/ABACUS.develop/source/src_lcao/ORB_read.cpp b/ABACUS.develop/source/src_lcao/ORB_read.cpp index 6377f76957..a0fb4b561e 100644 --- a/ABACUS.develop/source/src_lcao/ORB_read.cpp +++ b/ABACUS.develop/source/src_lcao/ORB_read.cpp @@ -376,15 +376,11 @@ void LCAO_Orbitals::Set_NonLocal(const int &it, int &n_projectors) { for(int is2=0;is2<2;is2++) { - complex coeff = complex(0.0,0.0); - for(int m=-l1-1;mdion(p1,p2) * soc.fcoef(it, is1, is2, ip1, ip2); if(p1 != p2) soc.fcoef(it, is1, is2, ip1, ip2) = complex(0.0,0.0); } diff --git a/ABACUS.develop/source/src_pw/pseudopot_cell_vnl.cpp b/ABACUS.develop/source/src_pw/pseudopot_cell_vnl.cpp index fd83eb7201..67af7a0a6f 100644 --- a/ABACUS.develop/source/src_pw/pseudopot_cell_vnl.cpp +++ b/ABACUS.develop/source/src_pw/pseudopot_cell_vnl.cpp @@ -301,15 +301,11 @@ void pseudopot_cell_vnl::init_vnl(UnitCell_pseudo &cell) { for(int is2=0;is2<2;is2++) { - complex coeff = complex(0.0,0.0); - for(int m=-l1-1;mlmaxkb ; - const int mj = soc.sph_ind(l2,j2,m,is2) + this->lmaxkb ; - coeff += soc.rotylm(m1,mi) * soc.spinor(l1,j1,m,is1)* - conj(soc.rotylm(m2,mj))*soc.spinor(l2,j2,m,is2); - } - soc.fcoef(it,is1,is2,ip,ip2) = coeff; + soc.set_fcoef(l1, l2, + is1, is2, + m1, m2, + j1, j2, + it, ip, ip2); } } } diff --git a/ABACUS.develop/source/src_pw/soc.cpp b/ABACUS.develop/source/src_pw/soc.cpp index 0b0f5cf938..2a236e921a 100644 --- a/ABACUS.develop/source/src_pw/soc.cpp +++ b/ABACUS.develop/source/src_pw/soc.cpp @@ -51,11 +51,6 @@ Soc::~Soc() delete[] p_rot; } -void Soc::init() -{ - return; -} - double Soc::spinor(const int l, const double j, const int m, const int spin) { @@ -101,55 +96,27 @@ double Soc::spinor(const int l, const double j, const int m, const int spin) void Soc::rot_ylm(const int lmax) { int m,n,l; - // this->rotylm = new complex* [2*lmax +1]; - // for(int i=0;i< 2 * lmax + 1; i++){ - // this->rotylm[i] = new complex [2*lmax +1]; - // for(int j = 0;j< 2 * lmax + 1;j++){ - // this->rotylm[i][j] = complex(0.0,0.0) ; - // } - // } - this->l_max = 2 * lmax + 1; + l_max_ = lmax; + l2plus1_ = 2 * l_max_ + 1; delete[] p_rot; - this->p_rot = new complex [this->l_max * this->l_max]; - ZEROS(p_rot, this->l_max*this->l_max); - this->rotylm(0,lmax) = complex(1.0, 0.0); + this->p_rot = new complex [l2plus1_ * l2plus1_]; + ZEROS(p_rot, l2plus1_ * l2plus1_); + this->p_rot[l_max_] = complex(1.0, 0.0); - l = lmax; + l = l_max_; for(int i=1;i<2*l+1;i += 2) { m = (i+1)/2; n = l-m; - this->rotylm(i,n) = complex(pow(-1.0 , m) / sqrt(2) , 0.0); - this->rotylm(i+1,n) = complex(0.0,-pow(-1.0 , m)/sqrt(2)); + this->p_rot[l2plus1_*i + n] = complex(pow(-1.0 , m) / sqrt(2) , 0.0); + this->p_rot[l2plus1_*(i+1) + n] = complex(0.0,-pow(-1.0 , m)/sqrt(2)); n = l+m; - this->rotylm(i,n) = complex(1.0/sqrt(2), 0.0); - this->rotylm(i+1,n) = complex(0.0, 1.0/sqrt(2)); + this->p_rot[l2plus1_*i + n] = complex(1.0/sqrt(2), 0.0); + this->p_rot[l2plus1_*(i+1) + n] = complex(0.0, 1.0/sqrt(2)); } return; } -/*void Soc::allocate_fcoef(const int nhm,const int nt) - { -//allocate -fcoef = new complex ****[nt]; -for(int i=0;i ***[2]; -for(int j=0;j<2;j++){ -fcoef[i][j] = new complex **[2]; -for(int k=0;k *[nhm]; -for(int l=0;l [nhm]; -for(int m=0;m(0.0,0.0); -} -} -} -} -} -return; -}*/ - int Soc::sph_ind(const int l, const double j, const int m, const int spin) { // This function calculates the m index of the spherical harmonic @@ -255,3 +222,27 @@ bool Soc::judge_parallel(double a[3], Vector3 b) jp = (fabs(cross)<1e-6); return jp; } + +void Soc::set_fcoef( + const int &l1, + const int &l2, + const int &is1, + const int &is2, + const int &m1, + const int &m2, + const double &j1, + const double &j2, + const int &it, + const int &ip1, + const int &ip2) +{ + complex coeff = complex(0.0,0.0); + for(int m=-l1-1;mfcoef(it,is1,is2,ip1,ip2) = coeff; +} \ No newline at end of file diff --git a/ABACUS.develop/source/src_pw/soc.h b/ABACUS.develop/source/src_pw/soc.h index 953aef7255..15964af682 100644 --- a/ABACUS.develop/source/src_pw/soc.h +++ b/ABACUS.develop/source/src_pw/soc.h @@ -41,8 +41,6 @@ class Soc Soc(); ~Soc(); - void init(); - double spinor(const int l, const double j, const int m, const int spin); int sph_ind(const int l, const double j, const int m, const int spin); @@ -50,12 +48,23 @@ class Soc void rot_ylm(const int lmax); // complex **rotylm; - complex *p_rot; - - complex & rotylm(const int &i1,const int &i2) - { return p_rot[l_max*i1 + i2]; } + const complex & rotylm(const int &i1,const int &i2) const + { return p_rot[l2plus1_*i1 + i2]; } Fcoef fcoef; + void set_fcoef( + const int &l1, + const int &l2, + const int &is1, + const int &is2, + const int &m1, + const int &m2, + const double &j1, + const double &j2, + const int &it, + const int &ip1, + const int &ip2 + ); Vector3 *m_loc; //magnetization for each atom axis double *angle1; @@ -70,6 +79,9 @@ class Soc bool judge_parallel(double a[3],Vector3 b); - int l_max; + int l_max_; + int l2plus1_; + complex *p_rot; + }; #endif