Vxc for SIC
88 views
Skip to first unread message
Xiaoming Wang
Dec 27, 2019, 11:50:58 AM12/27/19
to cp2k
Hello all,
I'm trying to implement the SIC scheme proposed in PRB 99, 235139 (2019) by Giustino et al.
The new SIC scheme is a minor modification on the one proposed by Mauri [PRB 71, 205210 (2005)].
Since the Mauri SIC is implemented in CP2K, incorporation of the Giustino scheme should not need
much coding work. The Hartree self-interaction part is similar for all the schemes, so I'm only interested
in the XC part. The implementation of the XC part of the Mauri scheme is coded in qs_vxc.F.
There are two types of SIC functionals for the Mauri scheme: Mauri_SPZ and Mauri_US.
The corresponding SIC corrrected XC functionals are:
Mauri_SPZ: Exc = Exc [ alpha, beta ] - Exc [ alpha - beta, 0 ]
Mauri_US: Exc = Exc [ alpha, beta ] - Exc [ alpha, beta ] + Exc [ beta, beta ]
And the XC potentials are:
Mauri_SPZ: Vxc_up = Vxc_up [ alpha, beta ] - Vxc_up [ alpha - beta, 0 ]
Vxc_dn = Vxc_dn [ alpha, beta ] + Vxc_up [ alpha - beta, 0 ]
Mauri_US: Vxc_up = 0
Vxc_dn = Vxc_up [ beta, beta ] + Vxc_dn [ beta, beta ]
For the Giustino scheme, the XC functional is:
Exc = 0.5*Exc [ alpha, beta ] + Exc [ beta, beta ] - 0.5* Exc [ beta - m, beta ]
where m = alpha - beta. Based on my understanding, the XC potentials are:
Vxc_up = 0.5*Vxc_up [ alpha, beta ] + 0.5*Vxc_up [ beta - m, beta ]
Vxc_dn = 0.5*Vxc_dn [ alpha, beta ] + Vxc_up [ beta, beta ] + Vxc_dn [ beta, beta ]
-0.5*Vxc_dn [ beta - m, beta ] - Vxc_up [ beta - m, beta ]
Please can anyone correct me if my understanding above is wrong.
I implemented the Giustino SIC based on the above equations by slightly modifying the US scheme
as already in CP2K. However, I never get convergence for the SCF calculations (see part of the log
file below). For the same system, the SPZ and US schemes both can easily get convergence.
So I'm wondering if anyone could give me any comments on this?
Best,
Xiaoming
--
1 OT CG 0.15E+00 14.2 0.00014039 -8153.9988425820 -8.15E+03
2 OT LS 0.10E+00 5.2 -8154.0035499679
3 OT CG 0.10E+00 10.6 0.00020232 -8154.0032891957 -4.45E-03
4 OT LS 0.46E-01 6.0 -8154.0015689085
5 OT CG 0.46E-01 10.8 0.00009816 -8154.0041296811 -8.40E-04
6 OT LS 0.31E-01 5.1 -8154.0048751897
7 OT CG 0.31E-01 13.1 0.00007670 -8154.0046628243 -5.33E-04
8 OT LS 0.26E-01 5.2 -8154.0051281540
9 OT CG 0.26E-01 12.0 0.00007395 -8154.0050568034 -3.94E-04
10 OT LS 0.21E-01 6.2 -8154.0053917653
11 OT CG 0.21E-01 11.2 0.00007210 -8154.0053264707 -2.70E-04
12 OT LS 0.16E-01 5.3 -8154.0055646935
13 OT CG 0.16E-01 10.4 0.00005042 -8154.0055112362 -1.85E-04
14 OT LS 0.15E-01 5.3 -8154.0056327656
15 OT CG 0.15E-01 11.4 0.00004942 -8154.0056270301 -1.16E-04
16 OT LS 0.14E-01 5.3 -8154.0057343493
17 OT CG 0.14E-01 11.5 0.00004857 -8154.0057262676 -9.92E-05
18 OT LS 0.13E-01 5.4 -8154.0058189369
19 OT CG 0.13E-01 10.7 0.00004785 -8154.0058097707 -8.35E-05
20 OT LS 0.11E-01 5.0 -8154.0058883901
----------------------------------- OT ---------------------------------------
----------------------------------- OT ---------------------------------------
1 OT CG 0.15E+00 21.8 0.00004727 -8154.0058790748 -6.93E-05
2 OT LS 0.11E+00 4.7 -8154.0065750586
3 OT CG 0.11E+00 10.9 0.00004074 -8154.0064374937 -5.58E-04
4 OT LS 0.13E+00 5.4 -8154.0072008089
5 OT CG 0.13E+00 11.3 0.00002783 -8154.0072885701 -8.51E-04
6 OT LS 0.51E+00 5.5 -8154.0078327851
7 OT CG 0.51E+00 11.1 0.00007874 -8154.0089341020 -1.65E-03
8 OT LS 0.21E+00 5.2 -8154.0051598661
9 OT CG 0.21E+00 11.3 0.00017344 -8154.0086760074 2.58E-04
10 OT LS 0.98E-01 5.3 -8154.0049564999
11 OT CG 0.98E-01 11.2 0.00006510 -8154.0082092588 4.67E-04
12 OT LS 0.58E-01 5.7 -8154.0086320114
13 OT CG 0.58E-01 10.8 0.00003585 -8154.0085028738 -2.94E-04
14 OT LS 0.10E+00 5.3 -8154.0088394530
15 OT CG 0.10E+00 12.0 0.00007085 -8154.0090772735 -5.74E-04
16 OT LS 0.44E-01 5.4 -8154.0084763907
17 OT CG 0.44E-01 11.1 0.00003408 -8154.0088708613 2.06E-04
18 OT LS 0.72E-01 5.3 -8154.0090917387
19 OT CG 0.72E-01 11.4 0.00005874 -8154.0092230991 -3.52E-04
20 OT LS 0.30E-01 5.5 -8154.0089057509
----------------------------------- OT ---------------------------------------
----------------------------------- OT ---------------------------------------
1 OT CG 0.15E+00 22.0 0.00001421 -8154.0091066057 1.16E-04
2 OT LS 0.18E+00 5.9 -8154.0092180076
3 OT CG 0.18E+00 10.7 0.00005904 -8154.0092398262 -1.33E-04
4 OT LS 0.96E-01 5.5 -8154.0094065127
5 OT CG 0.96E-01 11.2 0.00007661 -8154.0094075743 -1.68E-04
6 OT LS 0.42E-01 5.6 -8154.0088928557
7 OT CG 0.42E-01 10.9 0.00001444 -8154.0092323658 1.75E-04
8 OT LS 0.57E-01 5.1 -8154.0092666992
9 OT CG 0.57E-01 11.2 0.00001390 -8154.0092787816 -4.64E-05
10 OT LS 0.68E-01 5.2 -8154.0093186446
11 OT CG 0.68E-01 10.5 0.00002940 -8154.0093259578 -4.72E-05
12 OT LS 0.26E-01 5.0 -8154.0092084126
13 OT CG 0.26E-01 11.0 0.00001386 -8154.0092843540 4.16E-05
14 OT LS 0.13E-01 5.4 -8154.0092839420
15 OT CG 0.13E-01 11.0 0.00001387 -8154.0092841905 1.64E-07
16 OT LS 0.16E-01 5.4 -8154.0092934733
17 OT CG 0.16E-01 11.0 0.00001375 -8154.0092962332 -1.20E-05
18 OT LS 0.20E-01 5.4 -8154.0093077912
19 OT CG 0.20E-01 10.4 0.00001362 -8154.0093106170 -1.44E-05
20 OT LS 0.24E-01 5.3 -8154.0093242773
----------------------------------- OT ---------------------------------------
----------------------------------- OT ---------------------------------------
1 OT CG 0.15E+00 22.4 0.00001343 -8154.0093267803 -1.62E-05
2 OT LS 0.15E+00 5.4 -8154.0094121812
3 OT CG 0.15E+00 11.6 0.00007460 -8154.0094132278 -8.64E-05
4 OT LS 0.74E-01 5.3 -8154.0093027278
5 OT CG 0.74E-01 10.9 0.00007707 -8154.0094484130 -3.52E-05
6 OT LS 0.33E-01 5.6 -8154.0090697156
7 OT CG 0.33E-01 14.3 0.00001494 -8154.0093119312 1.36E-04
8 OT LS 0.46E-01 6.5 -8154.0093411677
--
Xiaoming Wang
Jan 1, 2020, 12:02:45 AM1/1/20
to cp...@googlegroups.com
Happy new year!
Really appreciate if anyone could confirm my derivation of the XC potential.
Best,
Xiaoming
Thomas Kühne
Jan 1, 2020, 5:02:49 AM1/1/20
to cp...@googlegroups.com
Dear Xiaoming,
happy new year as well! Though I can’t say anything definite, but having a quick glance
into Feliciano’s PRB paper I don’t see any apparent mistake. Regarding your convergence
issue (which shouldn’t be taken as a test to prove correctness or the opposite!) I want to
comment that the excessive usage of the OUTER_LOOP and hence too short inner loop
is not the smartest thing to do and may also lead to the behavior you observe. This is
particularly true when using the CG minimizer in which case you are only computing 10
gradients before restarting the inner loop with the default STEPSIZE of 0.15, which in my
experience is basically always too large except for trivially cases such as water …
Cheers,
Thomas
Xiaoming
--
You received this message because you are subscribed to the Google Groups "cp2k" group.
To unsubscribe from this group and stop receiving emails from it, send an email to cp2k+uns...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/a2bed8c8-339a-493f-8c48-27809a74d1a6%40googlegroups.com.
==============================
Thomas D. Kühne
Dynamics of Condensed Matter
Chair of Theoretical Chemistry
University of Paderborn
Warburger Str. 100
D-33098 Paderborn
Germany
Xiaoming Wang
Jan 1, 2020, 2:16:41 PM1/1/20
to cp2k
Hi Thomas,
Thanks for your comments. I tried different setups for the SCF loop and also tried the
diagonalization method but with no success. I have a concern about the XC potential
because with US or SPZ if the XC potentials were modified slightly different from the
default ones I got the similar convergence issue. However, with the default correct XC
potentials the US and SPZ scheme could get convergence fairly easily. Below is my
implementation which is coded in qs_vxc.F.
---
! The Sio-Giustino scheme
! E(rho_alpha,rho_beta)-b/2*(E(rho_alpha,rho_beta) + E(rho_beta-m,rho_beta) -2*E(rho_beta,rho_beta))
IF (dft_control%sic_method_id .EQ. sic_mauri_giustino .AND. .NOT. sic_scaling_b_zero) THEN
rho_r(1)%pw => rho_struct_r(2)%pw
rho_r(2)%pw => rho_struct_r(2)%pw
IF (rho_g_valid) THEN
rho_g(1)%pw => rho_struct_g(2)%pw
rho_g(2)%pw => rho_struct_g(2)%pw
ENDIF
IF (my_just_energy) THEN
exc_dd = xc_exc_calc(rho_r=rho_r, tau=tau, &
rho_g=rho_g, xc_section=xc_section, &
pw_pool=xc_pw_pool)
ELSE
! virial untested
CPASSERT(.NOT. compute_virial)
CALL xc_vxc_pw_create1(vxc_rho=my_vxc_rho, vxc_tau=my_vxc_tau, rho_r=rho_r, &
rho_g=rho_g, tau=tau, exc=exc_dd, &
xc_section=xc_section, &
pw_pool=xc_pw_pool, &
compute_virial=.FALSE., &
virial_xc=virial_xc_tmp)
END IF
! and take care of the potential
IF (.NOT. my_just_energy) THEN
!
vxc_rho(2)%pw%cr3d = vxc_rho(2)%pw%cr3d/2.0_dp + 2.0_dp*my_vxc_rho(1)%pw%cr3d
CALL pw_release(my_vxc_rho(1)%pw)
CALL pw_release(my_vxc_rho(2)%pw)
DEALLOCATE (my_vxc_rho)
ENDIF
ALLOCATE (rho_m_rspace(2), rho_m_gspace(2))
CALL pw_pool_create_pw(xc_pw_pool, rho_m_gspace(1)%pw, &
use_data=COMPLEXDATA1D, &
in_space=RECIPROCALSPACE)
CALL pw_pool_create_pw(xc_pw_pool, rho_m_rspace(1)%pw, &
use_data=REALDATA3D, &
in_space=REALSPACE)
CALL pw_copy(rho_struct_r(1)%pw, rho_m_rspace(1)%pw)
CALL pw_axpy(rho_struct_r(2)%pw, rho_m_rspace(1)%pw, alpha=-1._dp)
CALL pw_copy(rho_struct_g(1)%pw, rho_m_gspace(1)%pw)
CALL pw_axpy(rho_struct_g(2)%pw, rho_m_gspace(1)%pw, alpha=-1._dp)
CALL pw_pool_create_pw(xc_pw_pool, rho_m_gspace(2)%pw, &
use_data=COMPLEXDATA1D, &
in_space=RECIPROCALSPACE)
CALL pw_pool_create_pw(xc_pw_pool, rho_m_rspace(2)%pw, &
use_data=REALDATA3D, &
in_space=REALSPACE)
CALL pw_copy(rho_struct_r(2)%pw, rho_m_rspace(2)%pw)
CALL pw_copy(rho_struct_g(2)%pw, rho_m_gspace(2)%pw)
CALL pw_axpy(rho_m_rspace(1)%pw, rho_m_rspace(2)%pw, alpha=-1._dp)
CALL pw_axpy(rho_m_gspace(1)%pw, rho_m_gspace(2)%pw, alpha=-1._dp)
rho_r(1)%pw => rho_m_rspace(2)%pw
rho_r(2)%pw => rho_struct_r(2)%pw
IF (rho_g_valid) THEN
rho_g(1)%pw => rho_m_gspace(2)%pw
rho_g(2)%pw => rho_struct_g(2)%pw
ENDIF
IF (my_just_energy) THEN
exc_dm = xc_exc_calc(rho_r=rho_r, tau=tau, &
rho_g=rho_g, xc_section=xc_section, &
pw_pool=xc_pw_pool)
ELSE
! virial untested
CPASSERT(.NOT. compute_virial)
CALL xc_vxc_pw_create1(vxc_rho=my_vxc_rho, vxc_tau=my_vxc_tau, rho_r=rho_r, &
rho_g=rho_g, tau=tau, exc=exc_dm, &
xc_section=xc_section, &
pw_pool=xc_pw_pool, &
compute_virial=.FALSE., &
virial_xc=virial_xc_tmp)
END IF
IF (.NOT. my_just_energy) THEN
vxc_rho(1)%pw%cr3d = vxc_rho(1)%pw%cr3d/2.0_dp + my_vxc_rho(1)%pw%cr3d/2.0_dp
vxc_rho(2)%pw%cr3d = vxc_rho(2)%pw%cr3d - my_vxc_rho(2)%pw%cr3d/2.0_dp - my_vxc_rho(1)%pw%cr3d
CALL pw_release(my_vxc_rho(1)%pw)
CALL pw_release(my_vxc_rho(2)%pw)
DEALLOCATE (my_vxc_rho)
ENDIF
exc = exc/2.0_dp + exc_dd - exc_dm/2.0_dp
DO ispin = 1, 2
CALL pw_pool_give_back_pw(xc_pw_pool, rho_m_rspace(ispin)%pw)
CALL pw_pool_give_back_pw(xc_pw_pool, rho_m_gspace(ispin)%pw)
ENDDO
DEALLOCATE (rho_m_rspace)
DEALLOCATE (rho_m_gspace)
ENDIF
---
Best,
Xiaoming
To unsubscribe from this group and stop receiving emails from it, send an email to cp...@googlegroups.com.
To view this discussion on the web visit https://groups.google.com/d/msgid/cp2k/a2bed8c8-339a-493f-8c48-27809a74d1a6%40googlegroups.com.
Xiaoming Wang
Jan 14, 2020, 10:25:55 AM1/14/20
to cp2k
In my implementation, I dropped the scaling so both sic_scaling_a and sic_scaling_b are 1.
Reply all
Reply to author
Forward
0 new messages