Skip to content

Commit 5d1fc9b

Browse files
jinzx10jinzx10
committed
multicollinear LDA & GGA (implemented by Xiaoyu Zhang)
1 parent cead7ce commit 5d1fc9b

File tree

6 files changed

+178
-5
lines changed

6 files changed

+178
-5
lines changed

‎source/module_hamilt_general/module_xc/CMakeLists.txt‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -17,6 +17,12 @@ add_library(
1717
xc_functional_libxc_wrapper_xc.cpp
1818
xc_functional_libxc_wrapper_gcxc.cpp
1919
xc_functional_libxc_wrapper_tauxc.cpp
20+
NCLibxc/NCLibxc.cpp
21+
NCLibxc/NCLibxc.h
22+
NCLibxc/interface_to_libxc.cpp
23+
NCLibxc/interface_to_libxc.h
24+
NCLibxc/LebedevGrid.cpp
25+
xc_functional_NCLibxc_gga.cpp
2026
)
2127

2228
if(ENABLE_COVERAGE)

‎source/module_hamilt_general/module_xc/xc_functional.h‎

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -20,6 +20,10 @@
2020
#include "module_elecstate/module_charge/charge.h"
2121
#include "module_cell/unitcell.h"
2222

23+
// for multi-collinear approach
24+
using Matrix2x2 = std::array<std::array<std::complex<double>, 2>, 2>;
25+
std::vector<std::array<double, 4>> MakeAngularGrid(int grid_level);
26+
2327
class XC_Functional
2428
{
2529
public:
@@ -308,6 +312,16 @@ class XC_Functional
308312
static void hcth(const double rho, const double grho, double &sx, double &v1x, double &v2x);
309313
static void pwcorr(const double r, const double c[], double &g, double &dg);
310314

315+
//-------------------
316+
// xc_functional_NCLibxc_gga.cpp
317+
//-------------------
318+
// This file is for implementing multi-collinear appraoch for GGA functionals.
319+
static void postlibxc_gga(int xc_id, const std::vector<double>& rho_up, const std::vector<double>& rho_down,
320+
std::vector<double>& e, std::vector<double>& v1, std::vector<double>& v2,
321+
std::vector<double>& f1, std::vector<double>& f2, std::vector<double>& f3,const Charge* const chr,const double tpiba);
322+
static std::pair<std::vector<double>, std::vector<Matrix2x2>> gga_mc(int xc_id, const std::vector<double>& n,
323+
const std::vector<double>& mx, const std::vector<double>& my, const std::vector<double>& mz,const Charge* const chr,const double tpiba);
324+
311325
};
312326

313327
#endif //XC_FUNCTION_H

‎source/module_hamilt_general/module_xc/xc_functional_libxc_vxc.cpp‎

Lines changed: 85 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -9,6 +9,8 @@
99
#include "module_base/timer.h"
1010
#include "module_base/tool_title.h"
1111

12+
#include "NCLibxc/NCLibxc.h"
13+
1214
#include <xc.h>
1315

1416
#include <vector>
@@ -122,6 +124,88 @@ std::tuple<double,double,ModuleBase::matrix> XC_Functional_Libxc::v_xc_libxc( /
122124
if(4==PARAM.inp.nspin)
123125
{
124126
v = XC_Functional_Libxc::convert_v_nspin4(nrxx, chr, amag, v);
127+
if(PARAM.inp.multicolin)// added by Xiaoyu Zhang, Peking University, 2024.10.09. multicollinear method
128+
{
129+
ModuleBase::matrix v_nspin4(PARAM.inp.nspin, nrxx);
130+
etxc=0;
131+
vtxc=0;
132+
std::complex<double> twoi(0.0, 2.0);
133+
std::complex<double> two(2.0, 0.0);
134+
double e2=2.0;
135+
136+
NCLibxc::print_NCLibxc();
137+
if(!is_gga)
138+
{ //LDA
139+
#ifdef _OPENMP
140+
#pragma omp parallel for schedule(static, 1024) reduction(+:etxc) reduction(+:vtxc)
141+
#endif
142+
for(int ir = 0;ir<nrxx; ++ir){
143+
double exc = 0.0;
144+
for(int ipol=0;ipol<4;++ipol){
145+
v_nspin4(ipol, ir) = 0;
146+
}
147+
std::vector<double> n = {chr->rho[0][ir] + chr->rho_core[ir]};
148+
std::vector<double> mx = {chr->rho[1][ir]};
149+
std::vector<double> my = {chr->rho[2][ir]};
150+
std::vector<double> mz = {chr->rho[3][ir]};
151+
double amag = sqrt( pow(chr->rho[1][ir],2) + pow(chr->rho[2][ir],2) + pow(chr->rho[3][ir],2) );
152+
if (n[0] - amag <= 0.0) { //ensure the rhoup and rhodn to libxc are positive
153+
continue;
154+
}
155+
for(const int &id : func_id){
156+
auto [E_MC, V_MC] = NCLibxc::lda_mc(id, n, mx, my, mz);
157+
exc = e2*E_MC[0];
158+
v_nspin4(0, ir) += std::real(e2*(V_MC[0][0][0]+V_MC[0][1][1])/two);
159+
v_nspin4(1, ir) += std::real(e2*(V_MC[0][0][1]+V_MC[0][1][0])/two);
160+
v_nspin4(2, ir) += std::real(e2*(V_MC[0][1][0]-V_MC[0][0][1])/twoi);
161+
v_nspin4(3, ir) += std::real(e2*(V_MC[0][0][0]-V_MC[0][1][1])/two);
162+
etxc += exc * n[0];
163+
vtxc += v_nspin4(0, ir) * chr->rho[0][ir] + v_nspin4(1, ir) * mx[0] + v_nspin4(2, ir) * my[0] + v_nspin4(3, ir) * mz[0];// vtxc is used the calculation of the total energy(Ts more specifically), because abacus doesn't directly programme the kinetic operator and instead uses the sum of occupied orbital energy
164+
// the reason why i use chr->rho here is not clear. It is based on the implementation in v_xc.
165+
}
166+
}
167+
}
168+
else
169+
{ // gga
170+
// 初始化所有格点的数组
171+
std::vector<double> n(nrxx);
172+
std::vector<double> mx(nrxx);
173+
std::vector<double> my(nrxx);
174+
std::vector<double> mz(nrxx);
175+
std::vector<double> amag(nrxx);
176+
177+
for (int ir = 0; ir < nrxx; ++ir) {
178+
n[ir] = chr->rho[0][ir] + chr->rho_core[ir];
179+
mx[ir] = chr->rho[1][ir];
180+
my[ir] = chr->rho[2][ir];
181+
mz[ir] = chr->rho[3][ir];
182+
amag[ir] = sqrt(pow(chr->rho[1][ir], 2) + pow(chr->rho[2][ir], 2) + pow(chr->rho[3][ir], 2));
183+
}
184+
185+
// 确保 rhoup 和 rhodn 为正
186+
//for (int ir = 0; ir < nrxx; ++ir) {
187+
// if (n[ir] - amag[ir] <= 0.0) {
188+
// continue;
189+
// }
190+
// }
191+
192+
for (const int& id : func_id) {
193+
auto [E_MC, V_MC] = XC_Functional::gga_mc(id, n, mx, my, mz, chr, tpiba);
194+
#ifdef _OPENMP
195+
#pragma omp parallel for schedule(static, 1024) reduction(+:etxc) reduction(+:vtxc)
196+
#endif
197+
for (int ir = 0; ir < nrxx; ++ir) {
198+
double exc = e2 * E_MC[ir];
199+
v_nspin4(0, ir) += std::real(e2 * (V_MC[ir][0][0] + V_MC[ir][1][1]) / two);
200+
v_nspin4(1, ir) += std::real(e2 * (V_MC[ir][0][1] + V_MC[ir][1][0]) / two);
201+
v_nspin4(2, ir) += std::real(e2 * (V_MC[ir][1][0] - V_MC[ir][0][1]) / twoi);
202+
v_nspin4(3, ir) += std::real(e2 * (V_MC[ir][0][0] - V_MC[ir][1][1]) / two);
203+
etxc += exc * n[ir];
204+
vtxc += v_nspin4(0, ir) * chr->rho[0][ir] + v_nspin4(1, ir) * mx[ir] + v_nspin4(2, ir) * my[ir] + v_nspin4(3, ir) * mz[ir];
205+
}
206+
}
207+
}
208+
}
125209
}
126210

127211
//-------------------------------------------------
@@ -392,4 +476,4 @@ std::tuple<double,double,ModuleBase::matrix,ModuleBase::matrix> XC_Functional_Li
392476
return std::make_tuple( etxc, vtxc, std::move(v), std::move(vofk) );
393477
}
394478

395-
#endif
479+
#endif

‎source/module_hamilt_general/module_xc/xc_functional_vxc.cpp‎

Lines changed: 65 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -11,6 +11,16 @@
1111

1212
#ifdef USE_LIBXC
1313
#include "xc_functional_libxc.h"
14+
#include "NCLibxc/NCLibxc.h"
15+
#include <tuple>
16+
#include <iostream>
17+
#include <vector>
18+
#include <array>
19+
#include <cmath>
20+
#include <complex>
21+
#include <iomanip>
22+
#include <xc.h>
23+
#include <stdexcept>
1424
#endif
1525

1626
// [etxc, vtxc, v] = XC_Functional::v_xc(...)
@@ -39,6 +49,8 @@ std::tuple<double,double,ModuleBase::matrix> XC_Functional::v_xc(
3949
// the square of the e charge
4050
// in Rydeberg unit, so * 2.0.
4151
double e2 = 2.0;
52+
std::complex<double> twoi(0.0, 2.0);
53+
std::complex<double> two(2.0, 0.0);
4254

4355
double vanishing_charge = 1.0e-10;
4456

@@ -104,8 +116,8 @@ std::tuple<double,double,ModuleBase::matrix> XC_Functional::v_xc(
104116
}
105117
}
106118
}
107-
else if(PARAM.inp.nspin == 4)//noncollinear case added by zhengdy
108-
{
119+
else if(PARAM.inp.nspin == 4 && PARAM.inp.multicolin == 0)
120+
{ //noncollinear case (Kubler's locally collinear) added by zhengdy
109121
#ifdef _OPENMP
110122
#pragma omp parallel for reduction(+:etxc) reduction(+:vtxc)
111123
#endif
@@ -164,13 +176,63 @@ std::tuple<double,double,ModuleBase::matrix> XC_Functional::v_xc(
164176
}//end if
165177
// energy terms, local-density contributions
166178

179+
if (PARAM.inp.nspin == 4 && PARAM.inp.multicolin == 1
180+
&& (func_type == 0 || func_type == 1) )
181+
{ // noncollinear case added by Xiaoyu Zhang, Peking University, 2024.10.02. multicollinear method for lda Since NCLibxc needs libxc, this part codes will not be used.
182+
NCLibxc::print_NCLibxc();
183+
#ifdef _OPENMP
184+
#pragma omp parallel for reduction(+:etxc) reduction(+:vtxc)
185+
#endif
186+
for(int ir = 0;ir<nrxx; ir++)
187+
{
188+
if(!use_libxc){
189+
std::cerr << "Error: Multi-collinear approach does not support running without Libxc." << std::endl;
190+
std::exit(EXIT_FAILURE);
191+
}
192+
double exc = 0.0;
193+
for(int ipol=0;ipol<4;ipol++){
194+
v(ipol, ir) = 0;
195+
}
196+
NCLibxc nc_libxc;
197+
std::vector<double> n = {chr->rho[0][ir] + chr->rho_core[ir]};
198+
std::vector<double> mx = {chr->rho[1][ir]};
199+
std::vector<double> my = {chr->rho[2][ir]};
200+
std::vector<double> mz = {chr->rho[3][ir]};
201+
double amag = sqrt( pow(chr->rho[1][ir],2) + pow(chr->rho[2][ir],2) + pow(chr->rho[3][ir],2) );
202+
if (n[0] - amag <= 0.0) { //ensure the rhoup and rhodn to libxc are positive
203+
continue;
204+
}
205+
for(const int &id : func_id){
206+
auto [E_MC, V_MC] = NCLibxc::lda_mc(id, n, mx, my, mz);
207+
exc = e2*E_MC[0];
208+
v(0, ir) += std::real(e2*(V_MC[0][0][0]+V_MC[0][1][1])/two);
209+
v(1, ir) += std::real(e2*(V_MC[0][0][1]+V_MC[0][1][0])/two);
210+
v(2, ir) += std::real(e2*(V_MC[0][1][0]-V_MC[0][0][1])/twoi);
211+
v(3, ir) += std::real(e2*(V_MC[0][0][0]-V_MC[0][1][1])/two);
212+
etxc += exc * n[0];
213+
vtxc += v(0, ir) * chr->rho[0][ir] + v(1, ir) * chr->rho[1][ir] + v(2, ir) * chr->rho[2][ir] + v(3, ir) * chr->rho[3][ir];
214+
}
215+
}
216+
}
217+
218+
167219
// add gradient corrections (if any)
168220
// mohan modify 2009-12-15
169221

170222
// the dummy variable dum contains gradient correction to stress
171223
// which is not used here
172224
std::vector<double> dum;
173-
gradcorr(etxc, vtxc, v, chr, chr->rhopw, ucell, dum);
225+
if(PARAM.inp.nspin == 4 && PARAM.inp.multicolin == 1) {
226+
if(func_type == 2 || func_type == 3) {
227+
NCLibxc::print_NCLibxc();
228+
}
229+
if(func_type == 4 || func_type == 5) {
230+
std::cerr << "Error: Multi-collinear approach hasn't support hybrid functioanl yet" << std::endl;
231+
std::exit(EXIT_FAILURE);
232+
}
233+
} else {
234+
gradcorr(etxc, vtxc, v, chr, chr->rhopw, ucell, dum);
235+
}
174236

175237
// parallel code : collect vtxc,etxc
176238
// mohan add 2008-06-01

‎source/module_io/read_input_item_elec_stru.cpp‎

Lines changed: 7 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -280,7 +280,7 @@ void ReadInput::item_elec_stru()
280280
item.annotation = "1: single spin; 2: up and down spin; 4: noncollinear spin";
281281
read_sync_int(input.nspin);
282282
item.reset_value = [](const Input_Item& item, Parameter& para) {
283-
if (para.input.noncolin || para.input.lspinorb)
283+
if (para.input.noncolin || para.input.lspinorb || para.input.multicolin)
284284
{
285285
para.input.nspin = 4;
286286
}
@@ -590,6 +590,12 @@ void ReadInput::item_elec_stru()
590590
read_sync_bool(input.noncolin);
591591
this->add_item(item);
592592
}
593+
{
594+
Input_Item item("multicolin");
595+
item.annotation = "multi-collinear approach (in development)";
596+
read_sync_bool(input.multicolin);
597+
this->add_item(item);
598+
}
593599
{
594600
Input_Item item("soc_lambda");
595601
item.annotation = "The fraction of averaged SOC pseudopotential is "

‎source/module_parameter/input_parameter.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -121,6 +121,7 @@ struct Input_para
121121

122122
bool lspinorb = false; ///< consider the spin-orbit interaction
123123
bool noncolin = false; ///< using non-collinear-spin
124+
bool multicolin = false; ///< multi-collinear approach (in development)
124125
double soc_lambda = 1.0; ///< The fraction of averaged SOC pseudopotential
125126
///< is given by (1-soc_lambda)
126127

0 commit comments

Comments
 (0)