[CP2K-user] [CP2K:12698] Re: Vxc for SIC
Xiaoming Wang
wxia... at gmail.com
Wed Jan 1 19:16:40 UTC 2020
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
On Wednesday, January 1, 2020 at 5:02:49 AM UTC-5, tkuehne wrote:
>
> 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
>
> Am 01.01.2020 um 06:02 schrieb Xiaoming Wang <wx... at gmail.com
> <javascript:>>:
>
> Happy new year!
>
> Really appreciate if anyone could confirm my derivation of the XC
> potential.
>
> Best,
> Xiaoming
>
> On Friday, December 27, 2019 at 11:50:58 AM UTC-5, Xiaoming Wang wrote:
>>
>> 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
>>
>> --
>>
>>
>>
>>
> --
> 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 c... at googlegroups.com <javascript:>.
> To view this discussion on the web visit
> https://groups.google.com/d/msgid/cp2k/a2bed8c8-339a-493f-8c48-27809a74d1a6%40googlegroups.com
> <https://groups.google.com/d/msgid/cp2k/a2bed8c8-339a-493f-8c48-27809a74d1a6%40googlegroups.com?utm_medium=email&utm_source=footer>
> .
>
>
>
>
> ==============================
> Thomas D. Kühne
> Dynamics of Condensed Matter
> Chair of Theoretical Chemistry
> University of Paderborn
> Warburger Str. 100
> D-33098 Paderborn
> Germany
> td... at mail.upb.de <javascript:>
> +49/(0)5251/60-5726
>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://lists.cp2k.org/archives/cp2k-user/attachments/20200101/35269126/attachment.htm>
More information about the CP2K-user
mailing list