Fix: OperatorDFTU for multi-k and NSPIN=2 (#1570)

dyzheng · dyzheng · web-flow · commit c614a9a93705 · 2022-11-29T10:20:44.000+08:00
* Fix: OperatorDFTU for multi-k case

* Test: add 360* case for DFTU in Gamma_only

* Fix: solve conflicts

Co-authored-by: dyzheng &lt;zhengdy@bjaisi.com&gt;
diff --git a/source/module_hamilt/hamilt_lcao.cpp b/source/module_hamilt/hamilt_lcao.cpp
@@ -144,28 +144,29 @@ HamiltLCAO<T>::HamiltLCAO(
         this->opsd->add(exx);
     }
 
-    if (GlobalV::dft_plus_u)
+#ifdef __DEEPKS
+    if (GlobalV::deepks_scf)
     {
-        Operator<double>* dftu = new OperatorDFTU<OperatorLCAO<double>>(
+        Operator<double>* deepks = new DeePKS<OperatorLCAO<double>>(
+            loc_in,
             LM_in,
             nullptr,// no explicit call yet
             &(LM_in->Hloc)
         );
-        this->opsd->add(dftu);
+        this->opsd->add(deepks);
     }
+#endif
 
-#ifdef __DEEPKS
-    if (GlobalV::deepks_scf)
+    //end node should be OperatorDFTU
+    if (GlobalV::dft_plus_u)
     {
-        Operator<double>* deepks = new DeePKS<OperatorLCAO<double>>(
-            loc_in,
+        Operator<double>* dftu = new OperatorDFTU<OperatorLCAO<double>>(
             LM_in,
             nullptr,// no explicit call yet
             &(LM_in->Hloc)
         );
-        this->opsd->add(deepks);
+        this->opsd->add(dftu);
     }
-#endif
     return;
 }
 
@@ -288,16 +289,6 @@ HamiltLCAO<T>::HamiltLCAO(
         this->ops->add(nonlocal);
     }
 
-    if (GlobalV::dft_plus_u)
-    {
-        Operator<std::complex<double>>* dftu = new OperatorDFTU<OperatorLCAO<std::complex<double>>>(
-            LM_in,
-            nullptr,// no explicit call yet
-            &(LM_in->Hloc2)
-        );
-        this->ops->add(dftu);
-    }
-
 #ifdef __DEEPKS
     if (GlobalV::deepks_scf)
     {
@@ -310,6 +301,16 @@ HamiltLCAO<T>::HamiltLCAO(
         this->ops->add(deepks);
     }
 #endif
+    //end node should be OperatorDFTU
+    if (GlobalV::dft_plus_u)
+    {
+        Operator<std::complex<double>>* dftu = new OperatorDFTU<OperatorLCAO<std::complex<double>>>(
+            LM_in,
+            nullptr,// no explicit call yet
+            &(LM_in->Hloc2)
+        );
+        this->ops->add(dftu);
+    }
 
 }
 
diff --git a/source/module_hamilt/ks_lcao/operator_lcao.cpp b/source/module_hamilt/ks_lcao/operator_lcao.cpp
@@ -155,7 +155,11 @@ void OperatorLCAO<T>::init(const int ik_in)
         {   
             //only HK should be updated when cal_type=lcao_dftu
             //in cal_type=lcao_dftu, HK only need to update from one node
+            //dftu is a special operator, it should be the last node of chain
+            //Overlap matrix for ik is used by it, do folding first and then return
+            this->folding_fixed(ik_in);
             this->contributeHk(ik_in);
+            return;
 
             break;
         }
diff --git a/tests/integrate/260_NO_15_PK_PU_AF/INPUT b/tests/integrate/260_NO_15_PK_PU_AF/INPUT
@@ -1,12 +1,12 @@
 INPUT_PARAMETERS
 suffix    autotest
-ntype    3
+ntype    2
 nbands    40
 
 calculation    scf
-ecutwfc    20
+ecutwfc    50
 scf_thr    1.0e-4
-scf_nmax    2
+scf_nmax    3
 out_chg    0
 
 #init_chg    file
@@ -18,6 +18,7 @@ smearing_method    gaussian
 smearing_sigma    0.07
 
 #cal_force     1
+#cal_stress    1
 #force_thr_ev    0.01
 #relax_method    cg
 #relax_bfgs_init    0.5
@@ -28,7 +29,7 @@ mixing_beta    0.4
 
 ks_solver    scalapack_gvx
 basis_type    lcao
-gamma_only    1
+gamma_only    0
 symmetry    0
 nspin    2
 
diff --git a/tests/integrate/260_NO_15_PK_PU_AF/KPT b/tests/integrate/260_NO_15_PK_PU_AF/KPT
@@ -1,4 +1,4 @@
 K_POINTS
 0
 Gamma
-1 1 1 0 0 0
+2 2 2 0 0 0
diff --git a/tests/integrate/260_NO_15_PK_PU_AF/STRU b/tests/integrate/260_NO_15_PK_PU_AF/STRU
@@ -1,11 +1,9 @@
 ATOMIC_SPECIES
-Fe1 1.000 Fe_ONCV_PBE-1.0.upf
-Fe2 1.000 Fe_ONCV_PBE-1.0.upf
+Fe 1.000 Fe_ONCV_PBE-1.0.upf
 O  1.000 O_ONCV_PBE-1.0.upf
 
 NUMERICAL_ORBITAL
 Fe_gga_9au_100Ry_4s2p2d1f.orb
-Fe_gga_9au_100Ry_4s2p2d1f.orb
 O_gga_7au_100Ry_2s2p1d.orb
 
 LATTICE_CONSTANT
@@ -18,15 +16,11 @@ LATTICE_VECTORS
 ATOMIC_POSITIONS
 Direct
 
-Fe1
-1.0
-1
-0.00            0.00            0.00            1        1       1
-
-Fe2
--1.0
-1
-0.50            0.50            0.50            1        1       1
+Fe
+0.0
+2
+0.00            0.00            0.00         mag  1.0
+0.50            0.50            0.50         mag  -1.0
 
 O
 0.0
diff --git a/tests/integrate/260_NO_15_PK_PU_AF/jd b/tests/integrate/260_NO_15_PK_PU_AF/jd
@@ -0,0 +1 @@
+DFTU + NSPIN2, Fe2O2, multi-k case
diff --git a/tests/integrate/260_NO_15_PK_PU_AF/result.ref b/tests/integrate/260_NO_15_PK_PU_AF/result.ref
@@ -1,3 +1,3 @@
-etotref -7250.8489104634409159
-etotperatomref -1812.7122276159
-totaltimeref 7.86288
+etotref -7301.6331988753327096
+etotperatomref -1825.4082997188
+totaltimeref 18.05934
diff --git a/tests/integrate/360_NO_15_GO_PU_AF/INPUT b/tests/integrate/360_NO_15_GO_PU_AF/INPUT
@@ -0,0 +1,42 @@
+INPUT_PARAMETERS
+suffix    autotest
+ntype    2
+nbands    40
+
+calculation    scf
+ecutwfc    50
+scf_thr    1.0e-4
+scf_nmax    3
+out_chg    0
+
+#init_chg    file
+#out_dos    1
+#dos_sigma    0.05
+#out_band    1
+
+smearing_method    gaussian
+smearing_sigma    0.07
+
+cal_force     0
+cal_stress    0
+#force_thr_ev    0.01
+#relax_method    cg
+#relax_bfgs_init    0.5
+
+mixing_type    pulay
+mixing_beta    0.4
+#mixing_gg0    1.5
+
+ks_solver    scalapack_gvx
+basis_type    lcao
+gamma_only    1
+symmetry    0
+nspin    2
+
+#Parameter DFT+U
+dft_plus_u    1
+orbital_corr    2 2 -1
+hubbard_u    5.0 5.0 0.0
+hund_j    0.0 0.0 0.0
+pseudo_dir	../../PP_ORB
+orbital_dir	../../PP_ORB
diff --git a/tests/integrate/360_NO_15_GO_PU_AF/KPT b/tests/integrate/360_NO_15_GO_PU_AF/KPT
@@ -0,0 +1,4 @@
+K_POINTS
+0
+Gamma
+1 1 1 0 0 0
diff --git a/tests/integrate/360_NO_15_GO_PU_AF/STRU b/tests/integrate/360_NO_15_GO_PU_AF/STRU
@@ -0,0 +1,29 @@
+ATOMIC_SPECIES
+Fe 1.000 Fe_ONCV_PBE-1.0.upf
+O  1.000 O_ONCV_PBE-1.0.upf
+
+NUMERICAL_ORBITAL
+Fe_gga_9au_100Ry_4s2p2d1f.orb
+O_gga_7au_100Ry_2s2p1d.orb
+
+LATTICE_CONSTANT
+8.190
+
+LATTICE_VECTORS
+ 1.00    0.50     0.50
+ 0.50    1.00     0.50
+ 0.50    0.50     1.00
+ATOMIC_POSITIONS
+Direct
+
+Fe
+0.0
+2
+0.00            0.00            0.00         mag  1.0
+0.50            0.50            0.50         mag  -1.0
+
+O
+0.0
+2
+0.25            0.25            0.25            1        1       1
+0.75            0.75            0.75            1        1       1
diff --git a/tests/integrate/360_NO_15_GO_PU_AF/jd b/tests/integrate/360_NO_15_GO_PU_AF/jd
@@ -0,0 +1 @@
+DFTU + NSPIN2, Fe2O2, gamma-only case
diff --git a/tests/integrate/360_NO_15_GO_PU_AF/result.ref b/tests/integrate/360_NO_15_GO_PU_AF/result.ref
@@ -0,0 +1,3 @@
+etotref -7297.5434108291065058
+etotperatomref -1824.3858527073
+totaltimeref 12.30296
diff --git a/tests/integrate/CASES b/tests/integrate/CASES
@@ -175,6 +175,7 @@
 320_NO_GO_MD_MSST
 320_NO_GO_MD_NHC
 345_NO_GO_BS
+360_NO_15_GO_PU_AF
 #381_NO_GO_HSE
 
 #401_NP_KP_sp

Original file line number	Diff line number	Diff line change
`@@ -155,7 +155,11 @@ void OperatorLCAO<T>::init(const int ik_in)`
`155`	`155`	`{`
`156`	`156`	`//only HK should be updated when cal_type=lcao_dftu`
`157`	`157`	`//in cal_type=lcao_dftu, HK only need to update from one node`
	`158`	`+ //dftu is a special operator, it should be the last node of chain`
	`159`	`+ //Overlap matrix for ik is used by it, do folding first and then return`
	`160`	`+ this->folding_fixed(ik_in);`
`158`	`161`	`this->contributeHk(ik_in);`
	`162`	`+ return;`
`159`	`163`
`160`	`164`	`break;`
`161`	`165`	`}`
Original file line number	Diff line number	Diff line change
`@@ -1,4 +1,4 @@`
`1`	`1`	`K_POINTS`
`2`	`2`	`0`
`3`	`3`	`Gamma`
`4`		`-1 1 1 0 0 0`
	`4`	`+2 2 2 0 0 0`
-Original file line number
+Diff line change
@@ @@ -0,0 +1,4 @@ @@
 +K_POINTS
 +0
 +Gamma
 +1 1 1 0 0 0