Skip to content

Commit a6864fd

Browse files
dyzhengjiyuang
dyzheng
and
jiyuang
authored
Feature: add constrained DFT on PW and LCAO code. (#1530)
* add __LCAO * Revert "add __LCAO" This reverts commit 9f865cb. * add __LCAO * add example for constrain-DFT * Feature: add constrained DFT on PW and LCAO code. * Fix: update CASES * Fix: delete swp file Co-authored-by: jiyuang <[email protected]>
1 parent 641e99a commit a6864fd

File tree

15 files changed

+238
-2
lines changed

15 files changed

+238
-2
lines changed

‎docs/advanced/input_files/input-main.md‎

Lines changed: 4 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1910,7 +1910,10 @@ This part of variables are used to control berry phase and wannier90 interfacae
19101910
### ocp
19111911

19121912
- **Type**: Boolean
1913-
- **Description**: option for choose whether calcualting constrained DFT or not. Only used for TDDFT.
1913+
- **Description**: option for choose whether calcualting constrained DFT or not.
1914+
- For PW and LCAO code. if set to 1, occupations of bands will be setting of "ocp_set".
1915+
- For TDDFT in LCAO code. if set to 1, occupations will be constrained since second ionic step.
1916+
- For OFDFT, this feature can't be used.
19141917
- **Default**:0
19151918

19161919
### ocp_set

‎examples/fixed_occupations/INPUT‎

Lines changed: 26 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,26 @@
1+
INPUT_PARAMETERS
2+
ntype 2
3+
ecutwfc 100
4+
nbands 220
5+
scf_nmax 400
6+
scf_thr 1e-08
7+
basis_type lcao
8+
ks_solver genelpa
9+
mixing_type pulay
10+
mixing_beta 0.4
11+
mixing_gg0 1.5
12+
pseudo_dir ../../tests/PP_ORB
13+
orbital_dir ../../tests/PP_ORB
14+
nspin 2
15+
nelec 254
16+
gamma_only 1
17+
relax_nmax 50
18+
cal_force 1
19+
force_thr_ev 0.01
20+
relax_method cg
21+
out_force 1
22+
out_stru 1
23+
ocp 1
24+
ocp_set 128*1 92*0 125*1 0 0.5 0.5 92*0
25+
suffix autotest
26+
calculation relax

‎examples/fixed_occupations/KPT‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,4 @@
1+
K_POINTS
2+
0
3+
Gamma
4+
1 1 1 0 0 0

‎examples/fixed_occupations/README.md‎

Lines changed: 20 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,20 @@
1+
For a spin-polarized calculation (where `nspin=2` in the INPUT), let's say you want to excite a spin-up (spin component 1) or spin-down (spin component 2) electron from the ground state to the excited state.
2+
3+
To do this, you will need to do the following:
4+
5+
- Find out how many bands do you have. You can obtain that by search for `NBANDS` in the running_*.log file
6+
- Find out how many k points do you have. Let's call that number `N_K`.
7+
- Find out how many electrons are occupying spin components 1 and 2. Let's say there are `O_UP` electrons that occupy the spin up, `O_DN` electrons that occupy the spin down.
8+
- Next, calculate how many states are empty in both spins: empty states in the spin up = `U_UP` where `U_UP = NBANDS - O_UP`. Same applies for `U_DN`.
9+
10+
Finally, add the `ocp_set` as follows:
11+
12+
- `ocp_set [O_UP-1]*1.0 1*0.0 1*1.0 [U_UP-1]*0.0 <repeated N_K times> [O_DN]*1.0 [U_DN]*0.0 <repeated N_K times>`
13+
14+
Note the `<repeated N_K times>` above. It means: repeat `[O_UP-1]*1.0 1*0.0 1*1.0 [U_UP-1]*0.0` for `N_K` times.
15+
16+
Let's have an example. Let's say we have 8 kpoints in both spins, 216 bands, where the spin up electrons occupy 144 bands and the spin down occupy 142 bands. Then, here are the two tags:
17+
18+
- `ocp_set 143*1.0 1*0.0 1*1.0 71*0.0 143*1.0 1*0.0 1*1.0 71*0.0 143*1.0 1*0.0 1*1.0 71*0.0 143*1.0 1*0.0 1*1.0 71*0.0 143*1.0 1*0.0 1*1.0 71*0.0 143*1.0 1*0.0 1*1.0 71*0.0 143*1.0 1*0.0 1*1.0 71*0.0 143*1.0 1*0.0 1*1.0 71*0.0 142*1.0 74*0.0 142*1.0 74*0.0 142*1.0 74*0.0 142*1.0 74*0.0 142*1.0 74*0.0 142*1.0 74*0.0 142*1.0 74*0.0 142*1.0 74*0.0`
19+
20+
I hope that helps!

‎examples/fixed_occupations/STRU‎

Lines changed: 89 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,89 @@
1+
ATOMIC_SPECIES
2+
C 12.011 C_ONCV_PBE-1.0.upf
3+
N 14.007 N_ONCV_PBE-1.0.upf
4+
5+
NUMERICAL_ORBITAL
6+
C_gga_6au_100Ry_2s2p1d.orb
7+
N_gga_6au_100Ry_2s2p1d.orb
8+
9+
LATTICE_CONSTANT
10+
1.89035917
11+
12+
LATTICE_VECTORS
13+
7.13366 0.0 0.0
14+
0.0 7.13366 0.0
15+
0.0 0.0 7.13366
16+
17+
ATOMIC_POSITIONS
18+
Direct
19+
20+
C
21+
0.5
22+
62
23+
0.0582537569646 0.0584159945857 0.058415926138 1 1 1
24+
0.0584525979788 0.0582639495068 0.55830375247 1 1 1
25+
0.0584525787913 0.558303683638 0.0582639936653 1 1 1
26+
0.058396396539 0.558224625038 0.558224715196 1 1 1
27+
0.558382169405 0.058268770834 0.0582688599633 1 1 1
28+
0.558463495114 0.0582929217445 0.558232256517 1 1 1
29+
0.558463577944 0.558232219482 0.0582928425108 1 1 1
30+
0.0569183250919 0.808357457714 0.80835749987 1 1 1
31+
0.0594265208712 0.807932319168 0.309284859406 1 1 1
32+
0.059426577572 0.309284887397 0.807932327286 1 1 1
33+
0.0582708962134 0.308419898129 0.308419946959 1 1 1
34+
0.552210490811 0.808130335617 0.808130449403 1 1 1
35+
0.56007898337 0.805702513766 0.310376795781 1 1 1
36+
0.560078966382 0.310376785976 0.805702612021 1 1 1
37+
0.557023067703 0.30867815375 0.308678282438 1 1 1
38+
0.807389926447 0.0572791102987 0.807923278007 1 1 1
39+
0.808263165087 0.0584172339073 0.308437501942 1 1 1
40+
0.806315992388 0.556601918514 0.805714215285 1 1 1
41+
0.807993747385 0.559665043596 0.308694792822 1 1 1
42+
0.308318690294 0.0597626500343 0.808381454095 1 1 1
43+
0.308748283443 0.0572708057831 0.309266803587 1 1 1
44+
0.308559910935 0.564440597291 0.808110931472 1 1 1
45+
0.310971582241 0.556642745167 0.31034723917 1 1 1
46+
0.807389883243 0.807923230389 0.057279115659 1 1 1
47+
0.80631597543 0.8057142675 0.556601902675 1 1 1
48+
0.808263172696 0.308437396592 0.0584173597563 1 1 1
49+
0.807993775679 0.308694756167 0.559665150002 1 1 1
50+
0.308318745883 0.808381485917 0.0597626931394 1 1 1
51+
0.308559854579 0.808111002217 0.564440572799 1 1 1
52+
0.308748257913 0.309266717408 0.0572707489807 1 1 1
53+
0.310971599545 0.310347183879 0.556642802022 1 1 1
54+
0.680983978304 0.931152211061 0.681896897317 1 1 1
55+
0.682933887974 0.933778836822 0.182489476871 1 1 1
56+
0.693780582343 0.422913530255 0.690760004182 1 1 1
57+
0.68319883513 0.433488872734 0.183707455502 1 1 1
58+
0.181804381008 0.934890383796 0.68364685784 1 1 1
59+
0.183671586489 0.932999768842 0.184195770119 1 1 1
60+
0.185531923621 0.435690529465 0.681903564319 1 1 1
61+
0.182902486978 0.433739164197 0.18247500625 1 1 1
62+
0.680984011378 0.681896813139 0.931152163924 1 1 1
63+
0.693780567717 0.690759882798 0.422913626245 1 1 1
64+
0.682933887465 0.18248939249 0.93377894919 1 1 1
65+
0.683198846736 0.183707388858 0.433488918408 1 1 1
66+
0.181804340788 0.683646881687 0.934890437489 1 1 1
67+
0.185531780867 0.681903603326 0.435690784688 1 1 1
68+
0.183671487857 0.1841957994 0.932999755315 1 1 1
69+
0.18290246445 0.182475050719 0.433739324283 1 1 1
70+
0.932436869302 0.932997344069 0.932997363485 1 1 1
71+
0.934262246895 0.933811117962 0.433764310377 1 1 1
72+
0.934262256629 0.433764373413 0.933811190486 1 1 1
73+
0.932979530876 0.433485845208 0.433485766269 1 1 1
74+
0.433019856016 0.934927952578 0.934928009068 1 1 1
75+
0.43480141594 0.931137243218 0.4357074411 1 1 1
76+
0.434801393697 0.43570751419 0.931137256772 1 1 1
77+
0.425932008914 0.422912441828 0.422912456519 1 1 1
78+
0.932202849631 0.682608265065 0.682608293532 1 1 1
79+
0.934700498036 0.682314731818 0.181991381118 1 1 1
80+
0.934700509495 0.181991414051 0.682314829938 1 1 1
81+
0.933292316004 0.183377148465 0.183377232291 1 1 1
82+
0.434077020534 0.682626645492 0.18446069663 1 1 1
83+
0.434076992467 0.184460817095 0.682626623521 1 1 1
84+
0.434365571398 0.181982838292 0.18198289758 1 1 1
85+
86+
N
87+
0.5
88+
1
89+
0.420690650676 0.69596553718 0.695965578765 1 1 1

‎source/module_elecstate/elecstate.cpp‎

Lines changed: 16 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -17,9 +17,25 @@ const double* ElecState::getRho(int spin) const
1717
return &(this->charge->rho[spin][0]);
1818
}
1919

20+
void ElecState::fixed_weights(const double * const ocp_kb)
21+
{
22+
for (int ik = 0; ik < this->wg.nr; ik++)
23+
{
24+
for (int ib = 0; ib < this->wg.nc; ib++)
25+
{
26+
this->wg(ik, ib) = ocp_kb[ik * this->wg.nc + ib];
27+
}
28+
}
29+
this->skip_weights = true;
30+
}
31+
2032
void ElecState::calculate_weights()
2133
{
2234
ModuleBase::TITLE("ElecState", "calculate_weights");
35+
if(this->skip_weights)
36+
{
37+
return;
38+
}
2339

2440
// for test
2541
// std::cout << " gaussian_broadening = " << use_gaussian_broadening << std::endl;

‎source/module_elecstate/elecstate.h‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -62,6 +62,8 @@ class ElecState
6262

6363
// calculate wg from ekb
6464
virtual void calculate_weights();
65+
// use occupied weights from INPUT and skip calculate_weights
66+
void fixed_weights(const double * const ocp_kb);
6567

6668
virtual void print_psi(const psi::Psi<double>& psi_in)
6769
{
@@ -100,6 +102,9 @@ class ElecState
100102
protected:
101103
// calculate ebands for all k points and all occupied bands
102104
void calEBand();
105+
106+
private:
107+
bool skip_weights = false;
103108
};
104109

105110
} // namespace elecstate

‎source/module_esolver/esolver_ks_lcao.cpp‎

Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -206,6 +206,12 @@ void ESolver_KS_LCAO::Init(Input& inp, UnitCell& ucell)
206206
// this function belongs to ions LOOP
207207
int ion_step = 0;
208208
}
209+
210+
//Fix pelec->wg by ocp_kb
211+
if(GlobalV::ocp)
212+
{
213+
this->pelec->fixed_weights(GlobalV::ocp_kb.data());
214+
}
209215
}
210216

211217
void ESolver_KS_LCAO::cal_Energy(double& etot)

‎source/module_esolver/esolver_ks_pw.cpp‎

Lines changed: 7 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -149,7 +149,13 @@ namespace ModuleESolver
149149
}
150150

151151
//temporary
152-
this->Init_GlobalC(inp,ucell);
152+
this->Init_GlobalC(inp,ucell);
153+
154+
//Fix pelec->wg by ocp_kb
155+
if(GlobalV::ocp)
156+
{
157+
this->pelec->fixed_weights(GlobalV::ocp_kb.data());
158+
}
153159
}
154160

155161
void ESolver_KS_PW::beforescf(int istep)

‎tests/integrate/304_NO_GO_ocp/INPUT‎

Lines changed: 30 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,30 @@
1+
INPUT_PARAMETERS
2+
#Parameters (1.General)
3+
suffix autotest
4+
calculation scf
5+
ntype 1
6+
nbands 8
7+
symmetry 1
8+
pseudo_dir ../../PP_ORB
9+
orbital_dir ../../PP_ORB
10+
gamma_only 1
11+
12+
#Parameters (2.Iteration)
13+
ecutwfc 20
14+
scf_thr 1e-8
15+
scf_nmax 100
16+
17+
nspin 2
18+
#Parameters (3.Basis)
19+
basis_type lcao
20+
21+
#Parameters (4.Smearing)
22+
smearing_method gauss
23+
smearing_sigma 0.002
24+
25+
#Parameters (5.Mixing)
26+
mixing_type plain
27+
mixing_beta 0.7
28+
29+
ocp 1
30+
ocp_set 3*1 0 1 3*0 3*1 0.5 0.5 3*0

0 commit comments

Comments
 (0)