comp-chem/cp2k
1
0
Fork 0
mirror of https://github.com/cp2k/cp2k.git synced 2026-07-21 06:25:15 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions

CNEO: CNEO-DFT energy and force

Browse source 
Energy and force for constrained nuclear-electronic orbital density
functional theory (J. Chem. Theory Comput. 2025, 21, 16, 7865–7877).

`POTENTIAL CNEO` triggers CNEO calculations for certain atom kinds
(hydrogen is the most common, but no limit in principle) and nuclear
basis functions are supplied by `BASIS_SET NUC basis_name`. PB series
basis sets for proton are provided in `data/NUCLEAR_BASIS_SETS`.

The `MASS` can be tuned to simulate isotopes such as deuterium, however,
one should be careful that basis functions need to be changed
accordingly. One simple modification to PB4-D with scaled exponents that
is suitable for deuterium is also added to `data/NUCLEAR_BASIS_SETS`

Note that GAPW is mandatory for atoms with quantum nuclei in CNEO-DFT
calculations, and pseudopotentials or even `GPW_TYPE T` can be used for
atom kinds with conventional classical nuclei.

For the electronic part, k-points can be enabled. However, this k-points
setting will not affect the quantum nuclear part. This is because even
quantum nuclei are delocalized as compared to conventional classical
point charges, they are still relatively localized. Therefore, a
distinguishable-particle approximation is used for quantum nuclei. This
is similar to the Gamma-point approximation for indistinguishable
quantum nuclei, with the added benefit that each quantum nucleus only
involves a small Kohn-Sham equation with localized basis functions, and
multiple quantum nuclei create multiple coupled Kohn-Sham equations,
instead of a giant one in the indistinguishable fermion case.
This commit is contained in:
zc62 2025-08-31 12:49:45 -05:00 • committed by Ole Schütt
parent 59b61dd70d
commit 15c1bfc0d2
56 changed files with 4744 additions and 234 deletions
Download patch file Download diff file Expand all files Collapse all files

393
data/NUCLEAR_BASIS_SETS Normal file
View file

@ -0,0 +1,393 @@
# PB series nuclear basis sets are from
# J. Chem. Phys. 152, 244123 (2020); doi: 10.1063/5.0009233
# Table 1. Table S2 for PB4-F2a.
H PB4-D
9
1 0 0 1 1
1.9570000000 1.0000000000
1 0 0 1 1
8.7340000000 1.0000000000
1 0 0 1 1
16.0100000000 1.0000000000
1 0 0 1 1
31.9970000000 1.0000000000
1 1 1 1 1
9.4380000000 1.0000000000
1 1 1 1 1
13.7950000000 1.0000000000
1 1 1 1 1
24.0280000000 1.0000000000
1 2 2 1 1
10.5240000000 1.0000000000
1 2 2 1 1
19.0160000000 1.0000000000
#
# scale PB4-D's exponents by sqrt(m_{deuteron}/m_{proton})
# use this for deuterium
H PB4-D_D
9
1 0 0 1 1
2.7669291965160787 1.0000000000
1 0 0 1 1
12.34867634255055 1.0000000000
1 0 0 1 1
22.635940948504047 1.0000000000
1 0 0 1 1
45.23936305617014 1.0000000000
1 1 1 1 1
13.344035644720874 1.0000000000
1 1 1 1 1
19.504235189544865 1.0000000000
1 1 1 1 1
33.97229163714273 1.0000000000
1 2 2 1 1
14.879490477330203 1.0000000000
1 2 2 1 1
26.886012059759704 1.0000000000
#
# scale PB4-D's exponents by sqrt(m_{antimuon}/m_{proton})
# use this for muonium (experimental)
H PB4-D_Mu
9
1 0 0 1 1
0.656717216271848 1.0000000000
1 0 0 1 1
2.9308983990384876 1.0000000000
1 0 0 1 1
5.372530726884153 1.0000000000
1 0 0 1 1
10.737343264716566 1.0000000000
1 1 1 1 1
3.1671420987090966 1.0000000000
1 1 1 1 1
4.6292355638580185 1.0000000000
1 1 1 1 1
8.063158545007646 1.0000000000
1 2 2 1 1
3.5315748513259724 1.0000000000
1 2 2 1 1
6.381264478602688 1.0000000000
#
H PB4-F1
10
1 0 0 1 1
5.9730000000 1.0000000000
1 0 0 1 1
10.6450000000 1.0000000000
1 0 0 1 1
17.9430000000 1.0000000000
1 0 0 1 1
28.9500000000 1.0000000000
1 1 1 1 1
7.6040000000 1.0000000000
1 1 1 1 1
14.7010000000 1.0000000000
1 1 1 1 1
23.3080000000 1.0000000000
1 2 2 1 1
9.0110000000 1.0000000000
1 2 2 1 1
19.7870000000 1.0000000000
1 3 3 1 1
10.9140000000 1.0000000000
#
H PB4-F2
11
1 0 0 1 1
5.9730000000 1.0000000000
1 0 0 1 1
10.6450000000 1.0000000000
1 0 0 1 1
17.9430000000 1.0000000000
1 0 0 1 1
28.9500000000 1.0000000000
1 1 1 1 1
7.6040000000 1.0000000000
1 1 1 1 1
14.7010000000 1.0000000000
1 1 1 1 1
23.3080000000 1.0000000000
1 2 2 1 1
9.0110000000 1.0000000000
1 2 2 1 1
19.7870000000 1.0000000000
1 3 3 1 1
10.9140000000 1.0000000000
1 3 3 1 1
20.9850000000 1.0000000000
#
H PB4-F2a
11
1 0 0 1 1
5.2290000000 1.0000000000
1 0 0 1 1
9.3370000000 1.0000000000
1 0 0 1 1
17.0180000000 1.0000000000
1 0 0 1 1
30.5080000000 1.0000000000
1 1 1 1 1
8.4840000000 1.0000000000
1 1 1 1 1
13.9930000000 1.0000000000
1 1 1 1 1
23.8090000000 1.0000000000
1 2 2 1 1
9.3170000000 1.0000000000
1 2 2 1 1
19.2670000000 1.0000000000
1 3 3 1 1
10.5360000000 1.0000000000
1 3 3 1 1
19.5740000000 1.0000000000
#
H PB5-D
12
1 0 0 1 1
1.9080000000 1.0000000000
1 0 0 1 1
9.0510000000 1.0000000000
1 0 0 1 1
15.0510000000 1.0000000000
1 0 0 1 1
29.7660000000 1.0000000000
1 0 0 1 1
40.1350000000 1.0000000000
1 1 1 1 1
4.9070000000 1.0000000000
1 1 1 1 1
10.0880000000 1.0000000000
1 1 1 1 1
15.8930000000 1.0000000000
1 1 1 1 1
25.7740000000 1.0000000000
1 2 2 1 1
10.3520000000 1.0000000000
1 2 2 1 1
18.3580000000 1.0000000000
1 2 2 1 1
36.3660000000 1.0000000000
#
H PB5-F
14
1 0 0 1 1
4.1890000000 1.0000000000
1 0 0 1 1
6.2310000000 1.0000000000
1 0 0 1 1
14.6240000000 1.0000000000
1 0 0 1 1
20.4810000000 1.0000000000
1 0 0 1 1
48.5090000000 1.0000000000
1 1 1 1 1
2.3490000000 1.0000000000
1 1 1 1 1
7.5970000000 1.0000000000
1 1 1 1 1
18.5210000000 1.0000000000
1 1 1 1 1
30.5960000000 1.0000000000
1 2 2 1 1
8.9710000000 1.0000000000
1 2 2 1 1
17.9560000000 1.0000000000
1 2 2 1 1
21.2990000000 1.0000000000
1 3 3 1 1
10.3210000000 1.0000000000
1 3 3 1 1
26.9100000000 1.0000000000
#
H PB5-G
15
1 0 0 1 1
3.2830000000 1.0000000000
1 0 0 1 1
7.6130000000 1.0000000000
1 0 0 1 1
16.0770000000 1.0000000000
1 0 0 1 1
20.4570000000 1.0000000000
1 0 0 1 1
47.0110000000 1.0000000000
1 1 1 1 1
3.1670000000 1.0000000000
1 1 1 1 1
7.7850000000 1.0000000000
1 1 1 1 1
18.9850000000 1.0000000000
1 1 1 1 1
29.4250000000 1.0000000000
1 2 2 1 1
10.2390000000 1.0000000000
1 2 2 1 1
17.2980000000 1.0000000000
1 2 2 1 1
20.6820000000 1.0000000000
1 3 3 1 1
10.9300000000 1.0000000000
1 3 3 1 1
25.9720000000 1.0000000000
1 4 4 1 1
10.5130000000 1.0000000000
#
H PB6-D
15
1 0 0 1 1
2.5130000000 1.0000000000
1 0 0 1 1
4.8400000000 1.0000000000
1 0 0 1 1
9.0880000000 1.0000000000
1 0 0 1 1
16.2310000000 1.0000000000
1 0 0 1 1
31.1100000000 1.0000000000
1 0 0 1 1
38.7540000000 1.0000000000
1 1 1 1 1
2.1740000000 1.0000000000
1 1 1 1 1
5.8970000000 1.0000000000
1 1 1 1 1
9.9960000000 1.0000000000
1 1 1 1 1
16.3470000000 1.0000000000
1 1 1 1 1
25.7160000000 1.0000000000
1 2 2 1 1
2.9300000000 1.0000000000
1 2 2 1 1
11.0150000000 1.0000000000
1 2 2 1 1
17.2560000000 1.0000000000
1 2 2 1 1
25.5410000000 1.0000000000
#
H PB6-F
18
1 0 0 1 1
2.8120000000 1.0000000000
1 0 0 1 1
5.7030000000 1.0000000000
1 0 0 1 1
10.2020000000 1.0000000000
1 0 0 1 1
15.4500000000 1.0000000000
1 0 0 1 1
29.1350000000 1.0000000000
1 0 0 1 1
39.7530000000 1.0000000000
1 1 1 1 1
1.4170000000 1.0000000000
1 1 1 1 1
6.8840000000 1.0000000000
1 1 1 1 1
10.4500000000 1.0000000000
1 1 1 1 1
16.1870000000 1.0000000000
1 1 1 1 1
26.4470000000 1.0000000000
1 2 2 1 1
4.4830000000 1.0000000000
1 2 2 1 1
9.8700000000 1.0000000000
1 2 2 1 1
18.3130000000 1.0000000000
1 2 2 1 1
25.3470000000 1.0000000000
1 3 3 1 1
6.0010000000 1.0000000000
1 3 3 1 1
10.3220000000 1.0000000000
1 3 3 1 1
24.4090000000 1.0000000000
#
H PB6-G
20
1 0 0 1 1
2.0160000000 1.0000000000
1 0 0 1 1
3.0490000000 1.0000000000
1 0 0 1 1
6.3500000000 1.0000000000
1 0 0 1 1
8.9100000000 1.0000000000
1 0 0 1 1
19.5000000000 1.0000000000
1 0 0 1 1
29.2960000000 1.0000000000
1 1 1 1 1
2.6380000000 1.0000000000
1 1 1 1 1
3.7160000000 1.0000000000
1 1 1 1 1
5.9780000000 1.0000000000
1 1 1 1 1
16.2520000000 1.0000000000
1 1 1 1 1
25.1140000000 1.0000000000
1 2 2 1 1
2.2050000000 1.0000000000
1 2 2 1 1
3.6540000000 1.0000000000
1 2 2 1 1
9.1040000000 1.0000000000
1 2 2 1 1
24.0390000000 1.0000000000
1 3 3 1 1
3.8450000000 1.0000000000
1 3 3 1 1
9.9210000000 1.0000000000
1 3 3 1 1
22.5700000000 1.0000000000
1 4 4 1 1
10.0620000000 1.0000000000
1 4 4 1 1
24.2670000000 1.0000000000
#
H PB6-H
21
1 0 0 1 1
1.3860000000 1.0000000000
1 0 0 1 1
3.2490000000 1.0000000000
1 0 0 1 1
9.5620000000 1.0000000000
1 0 0 1 1
12.4850000000 1.0000000000
1 0 0 1 1
21.4290000000 1.0000000000
1 0 0 1 1
36.9310000000 1.0000000000
1 1 1 1 1
1.3990000000 1.0000000000
1 1 1 1 1
4.2400000000 1.0000000000
1 1 1 1 1
9.9300000000 1.0000000000
1 1 1 1 1
18.3760000000 1.0000000000
1 1 1 1 1
24.1190000000 1.0000000000
1 2 2 1 1
2.6070000000 1.0000000000
1 2 2 1 1
5.1410000000 1.0000000000
1 2 2 1 1
7.7500000000 1.0000000000
1 2 2 1 1
20.7680000000 1.0000000000
1 3 3 1 1
0.5090000000 1.0000000000
1 3 3 1 1
9.1290000000 1.0000000000
1 3 3 1 1
26.4080000000 1.0000000000
1 4 4 1 1
9.4450000000 1.0000000000
1 4 4 1 1
28.4070000000 1.0000000000
1 5 5 1 1
10.1930000000 1.0000000000

4
src/CMakeLists.txt
View file

@ -525,6 +525,10 @@ list(
qs_charge_mixing.F
qs_chargemol.F
qs_charges_types.F
qs_cneo_ggrid.F
qs_cneo_methods.F
qs_cneo_types.F
qs_cneo_utils.F
qs_collocate_density.F
qs_commutators.F
qs_condnum.F

8
src/aobasis/basis_set_container_types.F
View file

@ -46,7 +46,9 @@ MODULE basis_set_container_types
aux_opt_basis = 117, &
min_basis = 118, &
tda_k_basis = 119, &
rhoin_basis = 120
rhoin_basis = 120, &
nuclear_basis = 121, &
nuclear_soft_basis = 122
! **************************************************************************************************
TYPE basis_set_container_type
PRIVATE
@ -136,6 +138,10 @@ CONTAINS
basis_type_nr = aux_opt_basis
CASE ("RHOIN")
basis_type_nr = rhoin_basis
CASE ("NUC")
basis_type_nr = nuclear_basis
CASE ("NUC_SOFT")
basis_type_nr = nuclear_soft_basis
CASE DEFAULT
basis_type_nr = unknown_basis
END SELECT

7
src/basis_set_output.F
View file

@ -59,7 +59,7 @@ CONTAINS
INTEGER :: ikind, iunit, nkind, ounit
INTEGER, SAVE :: ncalls = 0
TYPE(cp_logger_type), POINTER :: logger
TYPE(gto_basis_set_type), POINTER :: aux_fit_basis, lri_aux_basis, orb_basis, &
TYPE(gto_basis_set_type), POINTER :: aux_fit_basis, lri_aux_basis, nuclear_basis, orb_basis, &
p_lri_aux_basis, ri_aux_basis, ri_hfx_basis, ri_hxc_basis, ri_xas_basis, tda_hfx_basis
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(qs_kind_type), POINTER :: qs_kind
@ -97,6 +97,7 @@ CONTAINS
CALL get_qs_kind(qs_kind, basis_set=aux_fit_basis, basis_type="AUX_FIT")
CALL get_qs_kind(qs_kind, basis_set=ri_xas_basis, basis_type="RI_XAS")
CALL get_qs_kind(qs_kind, basis_set=tda_hfx_basis, basis_type="TDA_HFX")
CALL get_qs_kind(qs_kind, basis_set=nuclear_basis, basis_type="NUC")
IF (ounit > 0) THEN
IF (ASSOCIATED(orb_basis)) THEN
bname = "local_orbital"
@ -134,6 +135,10 @@ CONTAINS
bname = "local_tda_hfx"
CALL basis_out(tda_hfx_basis, element_symbol, bname, iunit)
END IF
IF (ASSOCIATED(nuclear_basis)) THEN
bname = "local_nuc"
CALL basis_out(nuclear_basis, element_symbol, bname, iunit)
END IF
END IF
END DO

10
src/common/bibliography.F
View file

@ -90,7 +90,7 @@ MODULE bibliography
Blase2018, Blase2020, Bruneval2015, Golze2019, Gui2018, Jacquemin2017, Liu2020, &
Sander2015, Schreiber2008, vanSetten2015, Setyawan2010, Ahart2024, Knysh2024, &
Schambeck2024, Mewes2018, Sertcan2024, Drautz2019, Lysogorskiy2021, Bochkarev2024, &
VazdaCruz2021
VazdaCruz2021, Chen2025
CONTAINS
@ -1953,6 +1953,7 @@ CONTAINS
title="Graph Atomic Cluster Expansion for semilocal interactions beyond equivariant message passing", &
source="Phys. Rev. X", volume="14", pages="021036", &
year=2024, doi="10.1103/PhysRevX.14.021036")
CALL add_reference(key=VazdaCruz2021, &
authors=s2a("V. Vaz da Cruz", "S. Eckert", "A. Fohlisch"), &
title="TD-DFT simulations of K-edge resonant inelastic X-ray scattering within "// &
@ -1960,6 +1961,13 @@ CONTAINS
source="Phys. Chem. Chem. Phys.", volume="23", pages="1835-1848", &
year=2021, doi="10.1039/d0cp04726k")
CALL add_reference(key=Chen2025, &
authors=s2a("Z. Chen", "Y. Yang"), &
title="Periodic Constrained Nuclear-Electronic Orbital Density Functional Theory for "// &
"Nuclear Quantum Effects: Method Development and Application to Hydrogen Adsorption on Pt(111)", &
source="J. Chem. Theory Comput.", volume="21", pages="7865-7877", &
year=2025, doi="10.1021/acs.jctc.5c00837")
END SUBROUTINE add_all_references
END MODULE bibliography

2
src/core_ae.F
View file

@ -58,7 +58,7 @@ MODULE core_ae
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'core_ae'
PUBLIC :: build_core_ae, build_erfc
PUBLIC :: build_core_ae, build_erfc, verfc_force
CONTAINS

11
src/cp_ddapc_util.F
View file

@ -232,7 +232,7 @@ CONTAINS
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(pw_c1d_gs_type) :: rho_tot_g
TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
TYPE(pw_c1d_gs_type), POINTER :: rho0_s_gs, rho_core
TYPE(pw_c1d_gs_type), POINTER :: rho0_s_gs, rho_core, rhoz_cneo_s_gs
TYPE(pw_env_type), POINTER :: pw_env
TYPE(pw_pool_type), POINTER :: auxbas_pool
TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
@ -250,13 +250,14 @@ CONTAINS
need_f = .FALSE.
IF (PRESENT(calc_force)) need_f = calc_force
logger => cp_get_default_logger()
NULLIFY (dft_control, rho, rho_core, rho0_s_gs, pw_env, rho_g, rho_r, &
NULLIFY (dft_control, rho, rho_core, rho0_s_gs, rhoz_cneo_s_gs, pw_env, rho_g, rho_r, &
radii, inp_radii, particle_set, qs_kind_set, qs_charges, cp_ddapc_env)
CALL get_qs_env(qs_env=qs_env, &
dft_control=dft_control, &
rho=rho, &
rho_core=rho_core, &
rho0_s_gs=rho0_s_gs, &
rhoz_cneo_s_gs=rhoz_cneo_s_gs, &
pw_env=pw_env, &
qs_charges=qs_charges, &
particle_set=particle_set, &
@ -290,7 +291,13 @@ CONTAINS
CASE (do_full_density)
! Otherwise build the total QS density (electron+nuclei) in G-space
IF (dft_control%qs_control%gapw) THEN
IF (ASSOCIATED(rhoz_cneo_s_gs)) THEN
CALL pw_axpy(rhoz_cneo_s_gs, rho0_s_gs)
END IF
CALL pw_transfer(rho0_s_gs, rho_tot_g)
IF (ASSOCIATED(rhoz_cneo_s_gs)) THEN
CALL pw_axpy(rhoz_cneo_s_gs, rho0_s_gs, -1.0_dp)
END IF
ELSE
CALL pw_transfer(rho_core, rho_tot_g)
END IF

210
src/hartree_local_methods.F
View file

@ -27,6 +27,10 @@ MODULE hartree_local_methods
USE pw_poisson_types, ONLY: pw_poisson_periodic,&
pw_poisson_type
USE qs_charges_types, ONLY: qs_charges_type
USE qs_cneo_methods, ONLY: Vh_1c_nuc_integrals,&
calculate_rhoz_cneo
USE qs_cneo_types, ONLY: cneo_potential_type,&
rhoz_cneo_type
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
USE qs_grid_atom, ONLY: grid_atom_type
@ -221,26 +225,31 @@ CONTAINS
CHARACTER(LEN=*), PARAMETER :: routineN = 'Vh_1c_gg_integrals'
INTEGER :: bo(2), handle, iat, iatom, ikind, ipgf1, is1, iset1, iso, l_ang, llmax, lmax0, &
lmax0_2nd, lmax_0, m1, max_iso, max_iso_not0, max_s_harm, maxl, maxso, mepos, n1, nat, &
nchan_0, nkind, nr, nset, nsotot, nspins, num_pe
INTEGER :: bo(2), handle, iat, iatom, ikind, ipgf1, is1, iset1, iso, l_ang, llmax, &
llmax_nuc, lmax0, lmax0_2nd, lmax_0, m1, max_iso, max_iso_not0, max_iso_not0_nuc, &
max_s_harm, max_s_harm_nuc, maxl, maxl_nuc, maxso, maxso_nuc, mepos, n1, nat, nchan_0, &
nkind, nr, nset, nset_nuc, nsotot, nsotot_nuc, nspins, num_pe
INTEGER, ALLOCATABLE, DIMENSION(:) :: cg_n_list
INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cg_list
INTEGER, DIMENSION(:), POINTER :: atom_list, lmax, lmin, npgf
LOGICAL :: core_charge, l_2nd_local_rho, &
INTEGER, DIMENSION(:), POINTER :: atom_list, lmax, lmax_nuc, lmin, &
lmin_nuc, npgf, npgf_nuc
LOGICAL :: cneo, core_charge, l_2nd_local_rho, &
my_core_2nd, my_periodic, paw_atom
REAL(dp) :: back_ch, factor
REAL(dp) :: back_ch, ecoul_1_z_cneo, factor
REAL(dp), ALLOCATABLE, DIMENSION(:) :: gexp, sqrtwr
REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: aVh1b_00, aVh1b_hh, aVh1b_ss, g0_h_w
REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: aVh1b_00, aVh1b_00_nuc, aVh1b_hh, &
aVh1b_hh_nuc, aVh1b_ss, aVh1b_ss_nuc, &
g0_h_w
REAL(dp), DIMENSION(:), POINTER :: rrad_z, vrrad_z
REAL(dp), DIMENSION(:, :), POINTER :: g0_h, g0_h_2nd, gsph, rrad_0, Vh1_h, &
Vh1_s, vrrad_0, zet
REAL(dp), DIMENSION(:, :), POINTER :: g0_h, g0_h_2nd, gsph, gsph_nuc, rrad_0, &
Vh1_h, Vh1_s, vrrad_0, zet, zet_nuc
REAL(dp), DIMENSION(:, :, :), POINTER :: my_CG, Qlm_gg, Qlm_gg_2nd
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(cell_type), POINTER :: cell
TYPE(cneo_potential_type), POINTER :: cneo_potential
TYPE(dft_control_type), POINTER :: dft_control
TYPE(grid_atom_type), POINTER :: grid_atom
TYPE(gto_basis_set_type), POINTER :: basis_1c
TYPE(gto_basis_set_type), POINTER :: basis_1c, nuc_basis
TYPE(harmonics_atom_type), POINTER :: harmonics
TYPE(pw_env_type), POINTER :: pw_env
TYPE(pw_poisson_type), POINTER :: poisson_env
@ -250,6 +259,8 @@ CONTAINS
TYPE(rho0_mpole_type), POINTER :: rho0_mpole, rho0_mpole_2nd
TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set, rho_atom_set_2nd
TYPE(rho_atom_type), POINTER :: rho_atom
TYPE(rhoz_cneo_type), DIMENSION(:), POINTER :: rhoz_cneo_set
TYPE(rhoz_cneo_type), POINTER :: rhoz_cneo
TYPE(rhoz_type), DIMENSION(:), POINTER :: rhoz_set, rhoz_set_2nd
CALL timeset(routineN, handle)
@ -260,6 +271,7 @@ CONTAINS
NULLIFY (atom_list, grid_atom, harmonics)
NULLIFY (basis_1c, lmin, lmax, npgf, zet)
NULLIFY (gsph)
NULLIFY (rhoz_cneo, rhoz_cneo_set)
CALL get_qs_env(qs_env=qs_env, &
cell=cell, dft_control=dft_control, &
@ -273,7 +285,8 @@ CONTAINS
back_ch = qs_charges%background*cell%deth
! rhoz_set is not accessed in TDDFT
CALL get_local_rho(local_rho_set, rho_atom_set, rho0_atom_set, rho0_mpole, rhoz_set) ! for integral space
CALL get_local_rho(local_rho_set, rho_atom_set, rho0_atom_set, rho0_mpole, rhoz_set, &
rhoz_cneo_set) ! for integral space
! for forces we need a second local_rho_set
l_2nd_local_rho = .FALSE.
@ -304,18 +317,36 @@ CONTAINS
! Put to 0 the local hartree energy contribution from 1 center integrals
energy_hartree_1c = 0.0_dp
! Restore total quantum nuclear density to zero
qs_charges%total_rho1_hard_nuc = 0.0_dp
rho0_mpole%tot_rhoz_cneo_s = 0.0_dp
! Here starts the loop over all the atoms
DO ikind = 1, nkind
CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=nat)
NULLIFY (cneo_potential)
CALL get_qs_kind(qs_kind_set(ikind), &
grid_atom=grid_atom, &
harmonics=harmonics, ngrid_rad=nr, &
max_iso_not0=max_iso_not0, paw_atom=paw_atom)
max_iso_not0=max_iso_not0, paw_atom=paw_atom, &
cneo_potential=cneo_potential)
CALL get_qs_kind(qs_kind_set(ikind), &
basis_set=basis_1c, basis_type="GAPW_1C")
cneo = ASSOCIATED(cneo_potential)
IF (cneo .AND. tddft) &
CPABORT("Electronic TDDFT with CNEO quantum nuclei is not implemented.")
NULLIFY (nuc_basis)
max_iso_not0_nuc = 0
IF (cneo) THEN
CPASSERT(paw_atom)
CALL get_qs_kind(qs_kind_set(ikind), &
basis_set=nuc_basis, basis_type="NUC")
max_iso_not0_nuc = cneo_potential%harmonics%max_iso_not0
END IF
IF (paw_atom) THEN
!=========== PAW ===============
CALL get_gto_basis_set(gto_basis_set=basis_1c, lmax=lmax, lmin=lmin, &
@ -330,14 +361,9 @@ CONTAINS
ALLOCATE (gexp(nr))
ALLOCATE (sqrtwr(nr), g0_h_w(nr, 0:lmax_0))
NULLIFY (Vh1_h, Vh1_s)
ALLOCATE (Vh1_h(nr, max_iso_not0))
ALLOCATE (Vh1_s(nr, max_iso_not0))
ALLOCATE (aVh1b_hh(nsotot, nsotot))
ALLOCATE (aVh1b_ss(nsotot, nsotot))
ALLOCATE (aVh1b_00(nsotot, nsotot))
ALLOCATE (cg_list(2, nsoset(maxl)**2, max_s_harm), cg_n_list(max_s_harm))
NULLIFY (Qlm_gg, g0_h)
CALL get_rho0_mpole(rho0_mpole=rho0_mpole, ikind=ikind, &
@ -352,7 +378,38 @@ CONTAINS
END IF
nchan_0 = nsoset(lmax0)
IF (nchan_0 > max_iso_not0) CPABORT("channels for rho0 > # max of spherical harmonics")
IF (nchan_0 > MAX(max_iso_not0, max_iso_not0_nuc)) &
CPABORT("channels for rho0 > # max of spherical harmonics")
NULLIFY (Vh1_h, Vh1_s)
ALLOCATE (Vh1_h(nr, max_iso_not0))
ALLOCATE (Vh1_s(nr, MAX(max_iso_not0, max_iso_not0_nuc, nchan_0)))
NULLIFY (lmax_nuc, lmin_nuc, npgf_nuc, zet_nuc, gsph_nuc)
maxso_nuc = 0
maxl_nuc = -1
nset_nuc = 0
max_s_harm_nuc = 0
llmax_nuc = -1
nsotot_nuc = 0
IF (cneo) THEN
CALL get_gto_basis_set(gto_basis_set=nuc_basis, lmax=lmax_nuc, &
lmin=lmin_nuc, maxso=maxso_nuc, npgf=npgf_nuc, &
maxl=maxl_nuc, nset=nset_nuc, zet=zet_nuc)
max_s_harm_nuc = cneo_potential%harmonics%max_s_harm
llmax_nuc = cneo_potential%harmonics%llmax
nsotot_nuc = maxso_nuc*nset_nuc
ALLOCATE (gsph_nuc(nr, nsotot_nuc))
gsph_nuc = 0.0_dp
ALLOCATE (aVh1b_hh_nuc(nsotot_nuc, nsotot_nuc))
ALLOCATE (aVh1b_ss_nuc(nsotot_nuc, nsotot_nuc))
ALLOCATE (aVh1b_00_nuc(nsotot_nuc, nsotot_nuc))
END IF
ALLOCATE (cg_list(2, nsoset(MAX(maxl, maxl_nuc))**2, &
MAX(max_s_harm, max_s_harm_nuc)), &
cg_n_list(MAX(max_s_harm, max_s_harm_nuc)))
NULLIFY (rrad_z, my_CG)
my_CG => harmonics%my_CG
@ -382,6 +439,33 @@ CONTAINS
m1 = m1 + maxso
END DO ! iset1
IF (cneo) THEN
! initialize nuclear pmat, cpc and e_core to zero
DO iat = 1, nat
iatom = atom_list(iat)
rhoz_cneo_set(iatom)%pmat = 0.0_dp
rhoz_cneo_set(iatom)%cpc_h = 0.0_dp
rhoz_cneo_set(iatom)%cpc_s = 0.0_dp
rhoz_cneo_set(iatom)%e_core = 0.0_dp
rhoz_cneo_set(iatom)%ready = .FALSE.
END DO
! calculate nuclear gsph
m1 = 0
DO iset1 = 1, nset_nuc
n1 = nsoset(lmax_nuc(iset1))
DO ipgf1 = 1, npgf_nuc(iset1)
gexp(1:nr) = EXP(-zet_nuc(ipgf1, iset1)*grid_atom%rad2(1:nr))*sqrtwr(1:nr)
DO is1 = nsoset(lmin_nuc(iset1) - 1) + 1, nsoset(lmax_nuc(iset1))
iso = is1 + (ipgf1 - 1)*n1 + m1
l_ang = indso(1, is1)
gsph_nuc(1:nr, iso) = cneo_potential%rad2l(1:nr, l_ang)*gexp(1:nr)
END DO ! is1
END DO ! ipgf1
m1 = m1 + maxso_nuc
END DO ! iset1
END IF
! Distribute the atoms of this kind
num_pe = para_env%num_pe
mepos = para_env%mepos
@ -392,10 +476,10 @@ CONTAINS
rho_atom => rho_atom_set(iatom)
NULLIFY (rrad_z, vrrad_z, rrad_0, vrrad_0)
IF (core_charge) THEN
IF (core_charge .AND. .NOT. cneo) THEN
rrad_z => rhoz_set(ikind)%r_coef ! for density
END IF
IF (my_core_2nd) THEN
IF (my_core_2nd .AND. .NOT. cneo) THEN
IF (l_2nd_local_rho) THEN
vrrad_z => rhoz_set_2nd(ikind)%vr_coef ! for potential
ELSE
@ -415,17 +499,77 @@ CONTAINS
END IF
CALL Vh_1c_atom_potential(rho_atom, vrrad_0, &
grid_atom, my_core_2nd, vrrad_z, Vh1_h, Vh1_s, & ! core charge for potential (2nd)
grid_atom, my_core_2nd .AND. .NOT. cneo, & ! core charge for potential (2nd)
vrrad_z, Vh1_h, Vh1_s, &
nchan_0, nspins, max_iso_not0, factor)
IF (l_2nd_local_rho) rho_atom => rho_atom_set(iatom) ! rho_atom for density
ecoul_1_z_cneo = 0.0_dp
IF (cneo) THEN
rhoz_cneo => rhoz_cneo_set(iatom)
! Add the soft tail of nuclear Hartree potential to total Vh1_s first.
! vrho already contains the -zeff factor.
DO iso = 1, max_iso_not0_nuc
Vh1_s(:, iso) = Vh1_s(:, iso) + rhoz_cneo%vrho_rad_s(:, iso)
END DO
! Build nuclear 1c integrals according to Vh1_h from electronic density only
! and Vh1_s from electron density and last step (or initial guess) nuclear density
! Vh1_h_nuc = -Z*Vh1_h, Vh1_s_nuc = -Z*Vh1_s
CALL Vh_1c_nuc_integrals(rhoz_cneo, cneo_potential%zeff, &
aVh1b_hh_nuc, aVh1b_ss_nuc, aVh1b_00_nuc, Vh1_h, Vh1_s, &
max_iso_not0, max_iso_not0_nuc, &
max_s_harm_nuc, llmax_nuc, cg_list, cg_n_list, &
nset_nuc, npgf_nuc, lmin_nuc, lmax_nuc, nsotot_nuc, maxso_nuc, &
nchan_0, gsph_nuc, g0_h_w, cneo_potential%harmonics%my_CG, &
cneo_potential%Qlm_gg)
! Solve the nuclear 1c problem
CALL calculate_rhoz_cneo(rhoz_cneo, cneo_potential, cg_list, cg_n_list, nset_nuc, &
npgf_nuc, lmin_nuc, lmax_nuc, maxl_nuc, maxso_nuc)
! slm_int(iso=1) = sqrt(4*Pi), slm_int(iso>1) = 0, without using Lebedev grid
! when printing, nuclear density is positive, thus needing the minus sign
qs_charges%total_rho1_hard_nuc = qs_charges%total_rho1_hard_nuc - SQRT(fourpi) &
*SUM(rhoz_cneo%rho_rad_h(:, 1)*grid_atom%wr(:))
rho0_mpole%tot_rhoz_cneo_s = rho0_mpole%tot_rhoz_cneo_s - SQRT(fourpi) &
*SUM(rhoz_cneo%rho_rad_s(:, 1)*grid_atom%wr(:))
! Calculate the contributions to Ecoul coming from Vh1_h*rhoz_cneo_h and Vh1_s*rhoz_cneo_s.
! Self-interaction of the quantum nucleus is already removed in ecoul_1_z_cneo.
! rho already contains the -zeff factor.
DO iso = 1, MIN(max_iso_not0, max_iso_not0_nuc)
ecoul_1_z_cneo = ecoul_1_z_cneo + 0.5_dp* &
(SUM(Vh1_h(:, iso)*rhoz_cneo%rho_rad_h(:, iso)*grid_atom%wr(:)) &
- SUM(Vh1_s(:, iso)*rhoz_cneo%rho_rad_s(:, iso)*grid_atom%wr(:)))
END DO
DO iso = max_iso_not0 + 1, max_iso_not0_nuc
ecoul_1_z_cneo = ecoul_1_z_cneo - 0.5_dp* &
SUM(Vh1_s(:, iso)*rhoz_cneo%rho_rad_s(:, iso)*grid_atom%wr(:))
END DO
! Add nuclear Hartree potential to total Vh1_h after solving the nuclear 1c problem
! to avoid nuclear self-interaction.
! vrho already contains the -zeff factor.
! Here the min of two max_iso_not0's is chosen, because even when the nuclear one
! is larger, it is meaningless to let Vh1_h have higher angular momentum components
! as Vh1_h now is only used for the electronic part.
DO iso = 1, MIN(max_iso_not0, max_iso_not0_nuc)
Vh1_h(:, iso) = Vh1_h(:, iso) + rhoz_cneo%vrho_rad_h(:, iso)
END DO
END IF
CALL Vh_1c_atom_energy(energy_hartree_1c, ecoul_1c, rho_atom, rrad_0, &
grid_atom, iatom, core_charge, rrad_z, Vh1_h, Vh1_s, & ! core charge for density
nchan_0, nspins, max_iso_not0)
grid_atom, iatom, core_charge .AND. .NOT. cneo, & ! core charge for density
rrad_z, Vh1_h, Vh1_s, nchan_0, nspins, max_iso_not0)
IF (l_2nd_local_rho) rho_atom => rho_atom_set_2nd(iatom) ! rho_atom for potential (2nd)
IF (cneo) THEN
CALL set_ecoul_1c(ecoul_1c, iatom, ecoul_1_z=ecoul_1c(iatom)%ecoul_1_z + ecoul_1_z_cneo)
energy_hartree_1c = energy_hartree_1c + ecoul_1_z_cneo
END IF
CALL Vh_1c_atom_integrals(rho_atom, & ! results (int_local_h and int_local_s) written on rho_atom_2nd
! int_local_h and int_local_s are used in update_ks_atom
! on int_local_h mixed core / non-core
@ -439,12 +583,27 @@ CONTAINS
DEALLOCATE (aVh1b_hh)
DEALLOCATE (aVh1b_ss)
DEALLOCATE (aVh1b_00)
IF (ALLOCATED(aVh1b_hh_nuc)) DEALLOCATE (aVh1b_hh_nuc)
IF (ALLOCATED(aVh1b_ss_nuc)) DEALLOCATE (aVh1b_ss_nuc)
IF (ALLOCATED(aVh1b_00_nuc)) DEALLOCATE (aVh1b_00_nuc)
DEALLOCATE (Vh1_h, Vh1_s)
DEALLOCATE (cg_list, cg_n_list)
DEALLOCATE (gsph)
IF (ASSOCIATED(gsph_nuc)) DEALLOCATE (gsph_nuc)
DEALLOCATE (gexp)
DEALLOCATE (sqrtwr, g0_h_w)
IF (cneo) THEN
! broadcast nuclear pmat, cpc and e_core
DO iat = 1, nat
iatom = atom_list(iat)
CALL para_env%sum(rhoz_cneo_set(iatom)%pmat)
CALL para_env%sum(rhoz_cneo_set(iatom)%cpc_h)
CALL para_env%sum(rhoz_cneo_set(iatom)%cpc_s)
CALL para_env%sum(rhoz_cneo_set(iatom)%e_core)
rhoz_cneo_set(iatom)%ready = .TRUE.
END DO
END IF
ELSE
!=========== NO PAW ===============
! This term is taken care of using the core density as in GPW
@ -453,6 +612,8 @@ CONTAINS
END DO ! ikind
CALL para_env%sum(energy_hartree_1c)
CALL para_env%sum(qs_charges%total_rho1_hard_nuc)
CALL para_env%sum(rho0_mpole%tot_rhoz_cneo_s)
CALL timestop(handle)
@ -670,7 +831,7 @@ CONTAINS
n2 = nsoset(lmax(iset2))
DO ipgf2 = 1, npgf(iset2)
! with contributions to V1_s*rho0
DO iso = 1, nchan_0
DO iso = 1, MIN(nchan_0, max_iso_not0)
l_ang = indso(1, iso)
gVg_0 = SUM(Vh1_s(:, iso)*g0_h_w(:, l_ang))
DO icg = 1, cg_n_list(iso)
@ -732,4 +893,3 @@ CONTAINS
!%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
END MODULE hartree_local_methods

9
src/input_cp2k_print_dft.F
View file

@ -402,6 +402,15 @@ CONTAINS
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL keyword_create(keyword, __LOCATION__, &
name="sab_cneo", &
description="Activates the printing of the nuclear orbital "// &
"nuclear repulsion neighbor lists (erfc potential)", &
default_l_val=.FALSE., &
lone_keyword_l_val=.TRUE.)
CALL section_add_keyword(print_key, keyword)
CALL keyword_release(keyword)
CALL section_add_subsection(section, print_key)
CALL section_release(print_key)

14
src/pw_env_methods.F
View file

@ -91,6 +91,7 @@ MODULE pw_env_methods
pw_pool_p_type,&
pw_pool_release,&
pw_pools_dealloc
USE qs_cneo_types, ONLY: cneo_potential_type
USE qs_dispersion_types, ONLY: qs_dispersion_type
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
@ -779,10 +780,10 @@ CONTAINS
TYPE(qs_environment_type), POINTER :: qs_env
CHARACTER(len=*), PARAMETER :: routineN = 'compute_max_radius'
CHARACTER(LEN=8), DIMENSION(4), PARAMETER :: &
pbas = (/"ORB ", "AUX_FIT ", "MAO ", "HARRIS "/)
CHARACTER(LEN=8), DIMENSION(9), PARAMETER :: sbas = (/"ORB ", "AUX ", "RI_AUX ", &
"MAO ", "HARRIS ", "RI_HXC ", "RI_K ", "LRI_AUX ", "RHOIN "/)
CHARACTER(LEN=8), DIMENSION(10), PARAMETER :: sbas = (/"ORB ", "AUX ", "RI_AUX ", &
"MAO ", "HARRIS ", "RI_HXC ", "RI_K ", "LRI_AUX ", "RHOIN ", "NUC "/)
CHARACTER(LEN=8), DIMENSION(5), PARAMETER :: &
pbas = (/"ORB ", "AUX_FIT ", "MAO ", "HARRIS ", "NUC "/)
REAL(KIND=dp), PARAMETER :: safety_factor = 1.015_dp
INTEGER :: handle, ibasis_set_type, igrid_level, igrid_zet0_s, ikind, ipgf, iset, ishell, &
@ -792,6 +793,7 @@ CONTAINS
REAL(KIND=dp) :: alpha, core_charge, eps_gvg, eps_rho, &
max_rpgf0_s, maxradius, zet0_h, zetp
REAL(KIND=dp), DIMENSION(:, :), POINTER :: zeta, zetb
TYPE(cneo_potential_type), POINTER :: cneo_potential
TYPE(dft_control_type), POINTER :: dft_control
TYPE(gto_basis_set_type), POINTER :: orb_basis_set
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
@ -839,6 +841,9 @@ CONTAINS
! this should, at a give point be changed
! so that also for the core a multigrid is used
DO ikind = 1, nkind
NULLIFY (cneo_potential)
CALL get_qs_kind(qs_kind_set(ikind), cneo_potential=cneo_potential)
IF (ASSOCIATED(cneo_potential)) CYCLE
CALL get_qs_kind(qs_kind_set(ikind), &
alpha_core_charge=alpha, ccore_charge=core_charge)
IF (alpha > 0.0_dp .AND. core_charge .NE. 0.0_dp) THEN
@ -1068,4 +1073,3 @@ CONTAINS
END SUBROUTINE setup_diel_rs_grid
END MODULE pw_env_methods

19
src/qmmm_gpw_forces.F
View file

@ -113,11 +113,11 @@ CONTAINS
TYPE(cp_logger_type), POINTER :: logger
TYPE(dft_control_type), POINTER :: dft_control
TYPE(mp_para_env_type), POINTER :: para_env
TYPE(pw_c1d_gs_type), POINTER :: rho0_s_gs, rho_core
TYPE(pw_c1d_gs_type), POINTER :: rho0_s_gs, rho_core, rhoz_cneo_s_gs
TYPE(pw_env_type), POINTER :: pw_env
TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
TYPE(pw_pool_type), POINTER :: auxbas_pool
TYPE(pw_r3d_rs_type) :: rho_tot_r, rho_tot_r2
TYPE(pw_r3d_rs_type) :: rho_tot_r, rho_tot_r2, rho_tot_r3
TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r
TYPE(qs_energy_type), POINTER :: energy
TYPE(qs_ks_qmmm_env_type), POINTER :: ks_qmmm_env_loc
@ -130,7 +130,7 @@ CONTAINS
periodic = qmmm_env%periodic
IF (PRESENT(calc_force)) need_f = calc_force
NULLIFY (dft_control, ks_qmmm_env_loc, rho, pw_env, energy, Forces, &
Forces_added_charges, input_section, rho0_s_gs, rho_r)
Forces_added_charges, input_section, rho0_s_gs, rhoz_cneo_s_gs, rho_r)
CALL get_qs_env(qs_env=qs_env, &
rho=rho, &
rho_core=rho_core, &
@ -139,6 +139,7 @@ CONTAINS
para_env=para_env, &
input=input_section, &
rho0_s_gs=rho0_s_gs, &
rhoz_cneo_s_gs=rhoz_cneo_s_gs, &
dft_control=dft_control)
CALL qs_rho_get(rho, rho_r=rho_r)
@ -212,9 +213,21 @@ CONTAINS
CALL auxbas_pool%create_pw(rho_tot_r2)
CALL pw_transfer(rho0_s_gs, rho_tot_r2)
CALL pw_axpy(rho_tot_r2, rho_tot_r)
IF (ASSOCIATED(rhoz_cneo_s_gs)) THEN
CALL auxbas_pool%create_pw(rho_tot_r3)
CALL pw_transfer(rhoz_cneo_s_gs, rho_tot_r3)
CALL pw_axpy(rho_tot_r3, rho_tot_r)
CALL auxbas_pool%give_back_pw(rho_tot_r3)
END IF
CALL auxbas_pool%give_back_pw(rho_tot_r2)
ELSE
CALL pw_transfer(rho0_s_gs, rho_tot_r)
IF (ASSOCIATED(rhoz_cneo_s_gs)) THEN
CALL auxbas_pool%create_pw(rho_tot_r3)
CALL pw_transfer(rhoz_cneo_s_gs, rho_tot_r3)
CALL pw_axpy(rho_tot_r3, rho_tot_r)
CALL auxbas_pool%give_back_pw(rho_tot_r3)
END IF
!
! QM/MM Nuclear Electrostatic Potential already included through rho0
!

6
src/qs_charges_types.F
View file

@ -44,6 +44,9 @@ MODULE qs_charges_types
REAL(KIND=dp) :: total_rho_soft_gspace = -1.0_dp
REAL(KIND=dp), DIMENSION(:), POINTER :: total_rho1_hard => NULL(), &
total_rho1_soft => NULL()
REAL(KIND=dp) :: total_rho1_hard_nuc = -1.0_dp
REAL(KIND=dp) :: total_rho1_soft_nuc_rspace = -1.0_dp
REAL(KIND=dp) :: total_rho1_soft_nuc_lebedev = -1.0_dp
REAL(KIND=dp) :: background = -1.0_dp
END TYPE qs_charges_type
@ -79,6 +82,9 @@ CONTAINS
qs_charges%total_rho1_hard(:) = 0.0_dp
ALLOCATE (qs_charges%total_rho1_soft(nspins))
qs_charges%total_rho1_soft(:) = 0.0_dp
qs_charges%total_rho1_hard_nuc = 0.0_dp
qs_charges%total_rho1_soft_nuc_rspace = 0.0_dp
qs_charges%total_rho1_soft_nuc_lebedev = 0.0_dp
END SUBROUTINE qs_charges_create
! **************************************************************************************************

609
src/qs_cneo_ggrid.F Normal file
View file

@ -0,0 +1,609 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
! !
! SPDX-License-Identifier: GPL-2.0-or-later !
!--------------------------------------------------------------------------------------------------!
! **************************************************************************************************
!> \brief CNEO soft nuclear densities on the global grid
!> (see J. Chem. Theory Comput. 2025, 21, 16, 7865–7877)
!> \par History
!> 08.2025 created [zc62]
!> \author Zehua Chen
! **************************************************************************************************
MODULE qs_cneo_ggrid
USE ao_util, ONLY: exp_radius_very_extended
USE atomic_kind_types, ONLY: atomic_kind_type,&
get_atomic_kind
USE basis_set_types, ONLY: get_gto_basis_set,&
gto_basis_set_type
USE cell_types, ONLY: cell_type,&
pbc
USE cp_control_types, ONLY: dft_control_type
USE gaussian_gridlevels, ONLY: gaussian_gridlevel,&
gridlevel_info_type
USE grid_api, ONLY: GRID_FUNC_AB,&
collocate_pgf_product
USE kinds, ONLY: dp
USE message_passing, ONLY: mp_para_env_type
USE orbital_pointers, ONLY: ncoset
USE particle_types, ONLY: particle_type
USE pw_env_types, ONLY: pw_env_get,&
pw_env_type
USE pw_methods, ONLY: pw_axpy,&
pw_integrate_function,&
pw_transfer,&
pw_zero
USE pw_pool_types, ONLY: pw_pool_p_type,&
pw_pool_type,&
pw_pools_create_pws,&
pw_pools_give_back_pws
USE pw_types, ONLY: pw_c1d_gs_type,&
pw_r3d_rs_type
USE qs_cneo_types, ONLY: cneo_potential_type,&
get_cneo_potential,&
rhoz_cneo_type
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
USE qs_force_types, ONLY: qs_force_type
USE qs_integrate_potential, ONLY: integrate_pgf_product
USE qs_kind_types, ONLY: get_qs_kind,&
get_qs_kind_set,&
qs_kind_type
USE qs_rho0_types, ONLY: get_rho0_mpole,&
rho0_mpole_type
USE realspace_grid_types, ONLY: map_gaussian_here,&
realspace_grid_type,&
rs_grid_zero,&
transfer_rs2pw
USE rs_pw_interface, ONLY: potential_pw2rs
USE virial_types, ONLY: virial_type
#include "./base/base_uses.f90"
IMPLICIT NONE
PRIVATE
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_cneo_ggrid'
PUBLIC :: put_rhoz_cneo_s_on_grid, rhoz_cneo_s_grid_create, integrate_vhgg_rspace
CONTAINS
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param rho0 ...
!> \param rhoz_cneo_set ...
!> \param tot_rs_int ...
! **************************************************************************************************
SUBROUTINE put_rhoz_cneo_s_on_grid(qs_env, rho0, rhoz_cneo_set, tot_rs_int)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(rho0_mpole_type), POINTER :: rho0
TYPE(rhoz_cneo_type), DIMENSION(:), POINTER :: rhoz_cneo_set
REAL(KIND=dp), INTENT(OUT) :: tot_rs_int
CHARACTER(LEN=*), PARAMETER :: routineN = 'put_rhoz_cneo_s_on_grid'
INTEGER :: atom_a, group_size, handle, iatom, igrid_level, ikind, ipgf, iset, jpgf, jset, &
m1, maxco, maxsgf_set, my_pos, na1, natom, nb1, ncoa, ncob, npgf2, nseta, offset, sgfa, &
sgfb
INTEGER, DIMENSION(:), POINTER :: atom_list, la_max, la_min, npgfa, nsgfa
INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
LOGICAL :: paw_atom
LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: map_it2
REAL(KIND=dp) :: eps_rho_rspace, radius, scale, zeff, zetp
REAL(KIND=dp), DIMENSION(3) :: ra
REAL(KIND=dp), DIMENSION(:, :), POINTER :: p_block, pab, sphi_a, work, zeta
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(cell_type), POINTER :: cell
TYPE(cneo_potential_type), POINTER :: cneo_potential
TYPE(dft_control_type), POINTER :: dft_control
TYPE(gridlevel_info_type), POINTER :: gridlevel_info
TYPE(gto_basis_set_type), POINTER :: nuc_soft_basis
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(pw_c1d_gs_type), ALLOCATABLE, DIMENSION(:) :: mgrid_gspace
TYPE(pw_c1d_gs_type), POINTER :: rhoz_cneo_s_gs
TYPE(pw_env_type), POINTER :: pw_env
TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:) :: mgrid_rspace
TYPE(pw_r3d_rs_type), POINTER :: rhoz_cneo_s_rs
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_rho
TYPE(realspace_grid_type), POINTER :: rs_grid
CALL timeset(routineN, handle)
NULLIFY (atomic_kind_set, qs_kind_set, atom_list, cell, dft_control, &
first_sgfa, gridlevel_info, la_max, la_min, nuc_soft_basis, &
npgfa, nsgfa, p_block, pab, particle_set, pw_env, pw_pools, &
rs_grid, rs_rho, sphi_a, work, zeta)
CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
atomic_kind_set=atomic_kind_set, &
cell=cell, particle_set=particle_set, &
pw_env=pw_env, &
dft_control=dft_control)
! find maximum numbers
CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
maxco=maxco, &
maxsgf_set=maxsgf_set, &
basis_type="NUC_SOFT")
IF (maxco > 0) THEN
ALLOCATE (pab(maxco, maxco), work(maxco, maxsgf_set))
eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
! *** set up the pw multi-grids *** !
CPASSERT(ASSOCIATED(pw_env))
CALL pw_env_get(pw_env=pw_env, rs_grids=rs_rho, pw_pools=pw_pools, &
gridlevel_info=gridlevel_info)
CALL pw_pools_create_pws(pw_pools, mgrid_rspace)
CALL pw_pools_create_pws(pw_pools, mgrid_gspace)
! *** set up the rs multi-grids *** !
DO igrid_level = 1, gridlevel_info%ngrid_levels
CALL rs_grid_zero(rs_rho(igrid_level))
END DO
offset = 0
my_pos = mgrid_rspace(1)%pw_grid%para%group%mepos
group_size = mgrid_rspace(1)%pw_grid%para%group%num_pe
DO ikind = 1, SIZE(atomic_kind_set)
CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom)
IF (.NOT. paw_atom) CYCLE
NULLIFY (cneo_potential)
CALL get_qs_kind(qs_kind_set(ikind), cneo_potential=cneo_potential)
IF (.NOT. ASSOCIATED(cneo_potential)) CYCLE
NULLIFY (atom_list)
CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom, atom_list=atom_list)
NULLIFY (nuc_soft_basis)
CALL get_qs_kind(qs_kind_set(ikind), basis_set=nuc_soft_basis, basis_type="NUC_SOFT")
CALL get_gto_basis_set(gto_basis_set=nuc_soft_basis, lmax=la_max, &
lmin=la_min, zet=zeta, nset=nseta, npgf=npgfa, &
sphi=sphi_a, first_sgf=first_sgfa, nsgf_set=nsgfa)
CALL get_cneo_potential(cneo_potential, zeff=zeff)
m1 = MAXVAL(npgfa(1:nseta))
ALLOCATE (map_it2(m1, m1))
DO iatom = 1, natom
atom_a = atom_list(iatom)
IF (rhoz_cneo_set(atom_a)%ready) THEN
ra(:) = pbc(particle_set(atom_a)%r, cell)
p_block => rhoz_cneo_set(atom_a)%pmat
DO iset = 1, nseta
DO jset = 1, iset
! processor mapping
map_it2 = .FALSE.
DO ipgf = 1, npgfa(iset)
IF (jset == iset) THEN
npgf2 = ipgf
ELSE
npgf2 = npgfa(jset)
END IF
DO jpgf = 1, npgf2
zetp = zeta(ipgf, iset) + zeta(jpgf, jset)
igrid_level = gaussian_gridlevel(gridlevel_info, zetp)
rs_grid => rs_rho(igrid_level)
map_it2(ipgf, jpgf) = map_gaussian_here(rs_grid, cell%h_inv, ra, offset, group_size, my_pos)
END DO
END DO
! skip empty sets (not uncommon for soft nuclear basis)
IF (npgfa(iset) > 0 .AND. npgfa(jset) > 0) THEN
offset = offset + 1
END IF
!
IF (ANY(map_it2(1:npgfa(iset), 1:npgfa(jset)))) THEN
ncoa = npgfa(iset)*ncoset(la_max(iset))
sgfa = first_sgfa(1, iset)
ncob = npgfa(jset)*ncoset(la_max(jset))
sgfb = first_sgfa(1, jset)
! decontract density block
CALL dgemm("N", "N", ncoa, nsgfa(jset), nsgfa(iset), &
1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
p_block(sgfa, sgfb), SIZE(p_block, 1), &
0.0_dp, work(1, 1), maxco)
CALL dgemm("N", "T", ncoa, ncob, nsgfa(jset), &
1.0_dp, work(1, 1), maxco, &
sphi_a(1, sgfb), SIZE(sphi_a, 1), &
0.0_dp, pab(1, 1), maxco)
DO ipgf = 1, npgfa(iset)
IF (jset == iset) THEN
npgf2 = ipgf
ELSE
npgf2 = npgfa(jset)
END IF
DO jpgf = 1, npgf2
IF (map_it2(ipgf, jpgf)) THEN
zetp = zeta(ipgf, iset) + zeta(jpgf, jset)
igrid_level = gaussian_gridlevel(gridlevel_info, zetp)
rs_grid => rs_rho(igrid_level)
na1 = (ipgf - 1)*ncoset(la_max(iset))
nb1 = (jpgf - 1)*ncoset(la_max(jset))
radius = exp_radius_very_extended(la_min=la_min(iset), &
la_max=la_max(iset), &
lb_min=la_min(jset), &
lb_max=la_max(jset), &
ra=ra, rb=ra, rp=ra, &
zetp=zetp, eps=eps_rho_rspace, &
prefactor=zeff, cutoff=1.0_dp)
IF (jset == iset .AND. jpgf == ipgf) THEN
scale = -zeff ! nuclear charge density is positive
ELSE
scale = -2.0_dp*zeff ! symmetric density matrix
END IF
CALL collocate_pgf_product( &
la_max(iset), zeta(ipgf, iset), la_min(iset), &
la_max(jset), zeta(jpgf, jset), la_min(jset), &
ra, (/0.0_dp, 0.0_dp, 0.0_dp/), &
scale, pab, na1, nb1, rs_grid, &
radius=radius, ga_gb_function=GRID_FUNC_AB)
END IF
END DO
END DO
END IF
END DO
END DO
END IF
END DO
DEALLOCATE (map_it2)
END DO
DEALLOCATE (pab, work)
NULLIFY (rhoz_cneo_s_gs, rhoz_cneo_s_rs)
CALL get_rho0_mpole(rho0_mpole=rho0, &
rhoz_cneo_s_gs=rhoz_cneo_s_gs, &
rhoz_cneo_s_rs=rhoz_cneo_s_rs)
CALL pw_zero(rhoz_cneo_s_gs)
CALL pw_zero(rhoz_cneo_s_rs)
DO igrid_level = 1, gridlevel_info%ngrid_levels
CALL pw_zero(mgrid_rspace(igrid_level))
CALL transfer_rs2pw(rs=rs_rho(igrid_level), &
pw=mgrid_rspace(igrid_level))
END DO
DO igrid_level = 1, gridlevel_info%ngrid_levels
CALL pw_zero(mgrid_gspace(igrid_level))
CALL pw_transfer(mgrid_rspace(igrid_level), &
mgrid_gspace(igrid_level))
CALL pw_axpy(mgrid_gspace(igrid_level), rhoz_cneo_s_gs)
END DO
CALL pw_transfer(rhoz_cneo_s_gs, rhoz_cneo_s_rs)
tot_rs_int = pw_integrate_function(rhoz_cneo_s_rs, isign=-1)
! *** give back the multi-grids *** !
CALL pw_pools_give_back_pws(pw_pools, mgrid_gspace)
CALL pw_pools_give_back_pws(pw_pools, mgrid_rspace)
ELSE
tot_rs_int = 0.0_dp
END IF
CALL timestop(handle)
END SUBROUTINE put_rhoz_cneo_s_on_grid
! **************************************************************************************************
!> \brief ...
!> \param pw_env ...
!> \param rho0_mpole ...
! **************************************************************************************************
SUBROUTINE rhoz_cneo_s_grid_create(pw_env, rho0_mpole)
TYPE(pw_env_type), POINTER :: pw_env
TYPE(rho0_mpole_type), POINTER :: rho0_mpole
CHARACTER(len=*), PARAMETER :: routineN = 'rhoz_cneo_s_grid_create'
INTEGER :: handle
TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
CALL timeset(routineN, handle)
CPASSERT(ASSOCIATED(pw_env))
NULLIFY (auxbas_pw_pool)
CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
CPASSERT(ASSOCIATED(auxbas_pw_pool))
! reallocate rho0 on the global grid in real and reciprocal space
CPASSERT(ASSOCIATED(rho0_mpole))
! soft nuclear charge density in real space
IF (ASSOCIATED(rho0_mpole%rhoz_cneo_s_rs)) THEN
CALL rho0_mpole%rhoz_cneo_s_rs%release()
ELSE
ALLOCATE (rho0_mpole%rhoz_cneo_s_rs)
END IF
CALL auxbas_pw_pool%create_pw(rho0_mpole%rhoz_cneo_s_rs)
! soft nuclear charge density in reciprocal space
IF (ASSOCIATED(rho0_mpole%rhoz_cneo_s_gs)) THEN
CALL rho0_mpole%rhoz_cneo_s_gs%release()
ELSE
ALLOCATE (rho0_mpole%rhoz_cneo_s_gs)
END IF
CALL auxbas_pw_pool%create_pw(rho0_mpole%rhoz_cneo_s_gs)
CALL timestop(handle)
END SUBROUTINE rhoz_cneo_s_grid_create
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param v_rspace ...
!> \param para_env ...
!> \param calculate_forces ...
!> \param rhoz_cneo_set ...
!> \param kforce ...
! **************************************************************************************************
SUBROUTINE integrate_vhgg_rspace(qs_env, v_rspace, para_env, calculate_forces, rhoz_cneo_set, &
kforce)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(pw_r3d_rs_type), INTENT(IN) :: v_rspace
TYPE(mp_para_env_type), POINTER :: para_env
LOGICAL, INTENT(IN) :: calculate_forces
TYPE(rhoz_cneo_type), DIMENSION(:), POINTER :: rhoz_cneo_set
REAL(KIND=dp), INTENT(IN), OPTIONAL :: kforce
CHARACTER(LEN=*), PARAMETER :: routineN = 'integrate_vhgg_rspace'
INTEGER :: atom_a, group_size, handle, iatom, igrid_level, ikind, ipgf, iset, jpgf, jset, &
m1, maxco, maxsgf_set, my_pos, na1, natom, nb1, ncoa, ncob, npgf2, nseta, offset, sgfa, &
sgfb
INTEGER, DIMENSION(:), POINTER :: atom_list, la_max, la_min, npgfa, nsgfa
INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
LOGICAL :: paw_atom, use_virial
LOGICAL, ALLOCATABLE, DIMENSION(:, :) :: map_it2
REAL(KIND=dp) :: eps_rho_rspace, f0, fscale, radius, &
zeff, zetp
REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, ra
REAL(KIND=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b
REAL(KIND=dp), DIMENSION(:, :), POINTER :: hab, p_block, pab, sphi_a, vmat, work, &
zeta
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(cell_type), POINTER :: cell
TYPE(cneo_potential_type), POINTER :: cneo_potential
TYPE(dft_control_type), POINTER :: dft_control
TYPE(gridlevel_info_type), POINTER :: gridlevel_info
TYPE(gto_basis_set_type), POINTER :: nuc_soft_basis
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(pw_env_type), POINTER :: pw_env
TYPE(qs_force_type), DIMENSION(:), POINTER :: force
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(realspace_grid_type), DIMENSION(:), POINTER :: rs_rho
TYPE(realspace_grid_type), POINTER :: rs_grid
TYPE(virial_type), POINTER :: virial
CALL timeset(routineN, handle)
NULLIFY (atomic_kind_set, qs_kind_set, atom_list, cell, dft_control, &
first_sgfa, force, gridlevel_info, hab, la_max, la_min, &
nuc_soft_basis, npgfa, nsgfa, p_block, pab, particle_set, &
pw_env, rs_grid, rs_rho, sphi_a, virial, vmat, work, zeta)
CALL get_qs_env(qs_env=qs_env, &
atomic_kind_set=atomic_kind_set, &
qs_kind_set=qs_kind_set, &
cell=cell, &
dft_control=dft_control, &
force=force, pw_env=pw_env, &
particle_set=particle_set, &
virial=virial)
! find maximum numbers
CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
maxco=maxco, &
maxsgf_set=maxsgf_set, &
basis_type="NUC_SOFT")
IF (maxco > 0) THEN
fscale = 1.0_dp
IF (PRESENT(kforce)) THEN
fscale = kforce
END IF
ALLOCATE (pab(maxco, maxco), work(maxco, maxsgf_set), hab(maxco, maxco))
pab = 0.0_dp
eps_rho_rspace = dft_control%qs_control%eps_rho_rspace
use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
CPASSERT(ASSOCIATED(pw_env))
CALL pw_env_get(pw_env=pw_env, rs_grids=rs_rho, gridlevel_info=gridlevel_info)
! transform the potential on the rs_multigrids
CALL potential_pw2rs(rs_rho, v_rspace, pw_env)
offset = 0
my_pos = rs_rho(1)%desc%my_pos
group_size = rs_rho(1)%desc%group_size
DO ikind = 1, SIZE(atomic_kind_set, 1)
CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom)
IF (.NOT. paw_atom) CYCLE
NULLIFY (cneo_potential)
CALL get_qs_kind(qs_kind_set(ikind), cneo_potential=cneo_potential)
IF (.NOT. ASSOCIATED(cneo_potential)) CYCLE
NULLIFY (atom_list)
CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
NULLIFY (nuc_soft_basis)
CALL get_qs_kind(qs_kind_set(ikind), basis_set=nuc_soft_basis, basis_type="NUC_SOFT")
CALL get_gto_basis_set(gto_basis_set=nuc_soft_basis, lmax=la_max, &
lmin=la_min, zet=zeta, nset=nseta, npgf=npgfa, &
sphi=sphi_a, first_sgf=first_sgfa, nsgf_set=nsgfa)
CALL get_cneo_potential(cneo_potential, zeff=zeff)
m1 = MAXVAL(npgfa(1:nseta))
ALLOCATE (map_it2(m1, m1))
DO iatom = 1, natom
atom_a = atom_list(iatom)
ra(:) = pbc(particle_set(atom_a)%r, cell)
IF (rhoz_cneo_set(atom_a)%ready) THEN
p_block => rhoz_cneo_set(atom_a)%pmat
ELSE
NULLIFY (p_block)
END IF
vmat => rhoz_cneo_set(atom_a)%vmat
vmat = 0.0_dp
DO iset = 1, nseta
DO jset = 1, iset
! processor mapping
map_it2 = .FALSE.
DO ipgf = 1, npgfa(iset)
IF (jset == iset) THEN
npgf2 = ipgf
ELSE
npgf2 = npgfa(jset)
END IF
DO jpgf = 1, npgf2
zetp = zeta(ipgf, iset) + zeta(jpgf, jset)
igrid_level = gaussian_gridlevel(gridlevel_info, zetp)
rs_grid => rs_rho(igrid_level)
map_it2(ipgf, jpgf) = map_gaussian_here(rs_grid, cell%h_inv, ra, offset, group_size, my_pos)
END DO
END DO
! skip empty sets (not uncommon for soft nuclear basis)
IF (npgfa(iset) > 0 .AND. npgfa(jset) > 0) THEN
offset = offset + 1
END IF
!
IF (ANY(map_it2(1:npgfa(iset), 1:npgfa(jset)))) THEN
hab = 0.0_dp
IF (calculate_forces) THEN
force_a = 0.0_dp
force_b = 0.0_dp
END IF
IF (use_virial) THEN
my_virial_a = 0.0_dp
my_virial_b = 0.0_dp
END IF
ncoa = npgfa(iset)*ncoset(la_max(iset))
sgfa = first_sgfa(1, iset)
ncob = npgfa(jset)*ncoset(la_max(jset))
sgfb = first_sgfa(1, jset)
IF (calculate_forces .AND. ASSOCIATED(p_block)) THEN
! decontract density block
CALL dgemm("N", "N", ncoa, nsgfa(jset), nsgfa(iset), &
1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
p_block(sgfa, sgfb), SIZE(p_block, 1), &
0.0_dp, work(1, 1), maxco)
CALL dgemm("N", "T", ncoa, ncob, nsgfa(jset), &
1.0_dp, work(1, 1), maxco, &
sphi_a(1, sgfb), SIZE(sphi_a, 1), &
0.0_dp, pab(1, 1), maxco)
END IF
DO ipgf = 1, npgfa(iset)
IF (jset == iset) THEN
npgf2 = ipgf
ELSE
npgf2 = npgfa(jset)
END IF
DO jpgf = 1, npgf2
IF (map_it2(ipgf, jpgf)) THEN
zetp = zeta(ipgf, iset) + zeta(jpgf, jset)
igrid_level = gaussian_gridlevel(gridlevel_info, zetp)
rs_grid => rs_rho(igrid_level)
na1 = (ipgf - 1)*ncoset(la_max(iset))
nb1 = (jpgf - 1)*ncoset(la_max(jset))
radius = exp_radius_very_extended(la_min=la_min(iset), &
la_max=la_max(iset), &
lb_min=la_min(jset), &
lb_max=la_max(jset), &
ra=ra, rb=ra, rp=ra, &
zetp=zetp, eps=eps_rho_rspace, &
prefactor=zeff, cutoff=1.0_dp)
CALL integrate_pgf_product( &
la_max(iset), zeta(ipgf, iset), la_min(iset), &
la_max(jset), zeta(jpgf, jset), la_min(jset), &
ra, (/0.0_dp, 0.0_dp, 0.0_dp/), rs_grid, &
hab, pab=pab, o1=na1, o2=nb1, &
radius=radius, &
calculate_forces=calculate_forces, force_a=force_a, force_b=force_b, &
use_virial=use_virial, my_virial_a=my_virial_a, my_virial_b=my_virial_b)
f0 = 0.0_dp
IF (calculate_forces .OR. use_virial) THEN
IF (jset == iset .AND. jpgf == ipgf) THEN
f0 = -fscale*zeff
ELSE
f0 = -2.0_dp*fscale*zeff
END IF
END IF
IF (calculate_forces) THEN
force(ikind)%rho_cneo_nuc(1:3, iatom) = &
force(ikind)%rho_cneo_nuc(1:3, iatom) + f0*(force_a + force_b)
END IF
IF (use_virial) THEN
virial%pv_virial = virial%pv_virial + f0*(my_virial_a + my_virial_b)
END IF
! symmetrize
IF (iset == jset .AND. ipgf /= jpgf) THEN
hab(nb1 + 1:nb1 + ncoset(la_max(jset)), &
na1 + 1:na1 + ncoset(la_max(iset))) = &
TRANSPOSE(hab(na1 + 1:na1 + ncoset(la_max(iset)), &
nb1 + 1:nb1 + ncoset(la_max(jset))))
END IF
END IF
END DO
END DO
! contract the soft basis V_Hartree integral
work(1:ncoa, 1:nsgfa(jset)) = MATMUL(hab(1:ncoa, 1:ncob), &
sphi_a(1:ncob, sgfb:sgfb + nsgfa(jset) - 1))
vmat(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfa(jset) - 1) = &
vmat(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfa(jset) - 1) - zeff* &
MATMUL(TRANSPOSE(sphi_a(1:ncoa, sgfa:sgfa + nsgfa(iset) - 1)), &
work(1:ncoa, 1:nsgfa(jset)))
! symmetrize
IF (iset /= jset) THEN
vmat(sgfb:sgfb + nsgfa(jset) - 1, sgfa:sgfa + nsgfa(iset) - 1) = &
TRANSPOSE(vmat(sgfa:sgfa + nsgfa(iset) - 1, sgfb:sgfb + nsgfa(jset) - 1))
END IF
END IF
END DO
END DO
END DO
DEALLOCATE (map_it2)
END DO
DEALLOCATE (pab, work, hab)
DO ikind = 1, SIZE(atomic_kind_set, 1)
CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom)
IF (.NOT. paw_atom) CYCLE
NULLIFY (cneo_potential)
CALL get_qs_kind(qs_kind_set(ikind), cneo_potential=cneo_potential)
IF (.NOT. ASSOCIATED(cneo_potential)) CYCLE
NULLIFY (atom_list)
CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
DO iatom = 1, natom
atom_a = atom_list(iatom)
vmat => rhoz_cneo_set(atom_a)%vmat
CALL para_env%sum(vmat)
END DO
END DO
END IF
CALL timestop(handle)
END SUBROUTINE integrate_vhgg_rspace
END MODULE qs_cneo_ggrid

1457
src/qs_cneo_methods.F Normal file
View file

File diff suppressed because it is too large Load diff

437
src/qs_cneo_types.F Normal file
View file

@ -0,0 +1,437 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
! !
! SPDX-License-Identifier: GPL-2.0-or-later !
!--------------------------------------------------------------------------------------------------!
! **************************************************************************************************
!> \brief Types used by CNEO-DFT
!> (see J. Chem. Theory Comput. 2025, 21, 16, 7865–7877)
!> \par History
!> 08.2025 created [zc62]
!> \author Zehua Chen
! **************************************************************************************************
MODULE qs_cneo_types
USE kinds, ONLY: dp
USE periodic_table, ONLY: ptable
USE qs_harmonics_atom, ONLY: deallocate_harmonics_atom,&
harmonics_atom_type
#include "./base/base_uses.f90"
IMPLICIT NONE
PRIVATE
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_cneo_types'
! Essential matrices, density and potential for each quantum nucleus
TYPE rhoz_cneo_type
LOGICAL :: ready = .FALSE. ! if pmat is ready. Useful in first iter
REAL(dp), DIMENSION(:, :), &
POINTER :: pmat => Null(), & ! nuclear density matrix
core => Null(), & ! nuclear core Hamiltonian
vmat => Null(), & ! nuclear Hartree from soft basis
fmat => Null(), & ! Fock = core + Hartree
wfn => Null() ! nuclear orbital coefficients
REAL(dp), DIMENSION(3) &
:: f = (/0.0_dp, 0.0_dp, 0.0_dp/) ! Lagrange multiplier for CNEO
REAL(dp) :: e_core = 0.0_dp ! nuclear core energy
REAL(dp), DIMENSION(:, :), &
POINTER :: cpc_h => Null(), & ! decontracted density matrix
cpc_s => Null(), & ! decontracted density matrix, soft tail
rho_rad_h => Null(), & ! density on radial grid
rho_rad_s => Null(), & ! density on radial grid, soft tail
vrho_rad_h => Null(), & ! potential on radial grid
vrho_rad_s => Null(), & ! potential on radial grid, soft tail
ga_Vlocal_gb_h => Null(), & ! local Hartree integral
ga_Vlocal_gb_s => Null() ! local Hartree integral, soft tail
END TYPE rhoz_cneo_type
! S, T, transformation matrices and distance are shared by the same kind,
! since they are not affected by the position of basis center
TYPE cneo_potential_type
INTEGER :: z = 0 ! atomic number
REAL(dp) :: zeff = 0.0_dp, & ! zeff = REAL(z)
mass = 0.0_dp ! atomic mass in dalton
INTEGER, DIMENSION(:), &
POINTER :: elec_conf => Null()
INTEGER :: nsgf = 0, & ! nuclear basis set is usually uncontracted
nne = 0, & ! number of linear-independent basis functions
npsgf = 0, & ! number of primitive SGFs
nsotot = 0 ! maxso * nset
REAL(dp), DIMENSION(:, :), &
POINTER :: my_gcc_h => Null(), & ! 3D-normalized contraction coefficients
my_gcc_s => Null(), & ! contraction coefficients for the soft tail
ovlp => Null(), & ! nuclear basis overlap matrix (unused)
kin => Null(), & ! nuclear kinetic energy matrix
utrans => Null() ! nuclear basis transformation matrix
REAL(dp), DIMENSION(:, :, :), &
POINTER :: distance => Null() ! distance from nuclear basis center
TYPE(harmonics_atom_type), &
POINTER :: harmonics => Null() ! most of the data will be missing
REAL(dp), DIMENSION(:, :, :), &
POINTER :: Qlm_gg => Null(), & ! multipole expansion of nuclear gg
gg => Null() ! precompute and store gg on radial grid
REAL(dp), DIMENSION(:, :, :, :), &
POINTER :: vgg => Null() ! precompute and store vgg on grid
INTEGER, DIMENSION(:), &
POINTER :: n2oindex => Null(), & ! new to old index
o2nindex => Null() ! old to new index
REAL(dp), DIMENSION(:, :), &
POINTER :: rad2l => Null(), & ! store my own rad2l
oorad2l => Null() ! store my own oorad2l
END TYPE cneo_potential_type
! Public Types
PUBLIC :: cneo_potential_type, rhoz_cneo_type
! Public Subroutine
PUBLIC :: allocate_cneo_potential, allocate_rhoz_cneo_set, deallocate_cneo_potential, &
deallocate_rhoz_cneo_set, get_cneo_potential, set_cneo_potential, write_cneo_potential
CONTAINS
! **************************************************************************************************
!> \brief ...
!> \param rhoz_cneo_set ...
!> \param natom ...
! **************************************************************************************************
SUBROUTINE allocate_rhoz_cneo_set(rhoz_cneo_set, natom)
TYPE(rhoz_cneo_type), DIMENSION(:), POINTER :: rhoz_cneo_set
INTEGER, INTENT(IN) :: natom
IF (ASSOCIATED(rhoz_cneo_set)) THEN
CALL deallocate_rhoz_cneo_set(rhoz_cneo_set)
END IF
ALLOCATE (rhoz_cneo_set(natom))
END SUBROUTINE allocate_rhoz_cneo_set
! **************************************************************************************************
!> \brief ...
!> \param rhoz_cneo ...
! **************************************************************************************************
SUBROUTINE deallocate_rhoz_cneo(rhoz_cneo)
TYPE(rhoz_cneo_type), POINTER :: rhoz_cneo
IF (ASSOCIATED(rhoz_cneo)) THEN
IF (ASSOCIATED(rhoz_cneo%pmat)) &
DEALLOCATE (rhoz_cneo%pmat)
IF (ASSOCIATED(rhoz_cneo%core)) &
DEALLOCATE (rhoz_cneo%core)
IF (ASSOCIATED(rhoz_cneo%vmat)) &
DEALLOCATE (rhoz_cneo%vmat)
IF (ASSOCIATED(rhoz_cneo%fmat)) &
DEALLOCATE (rhoz_cneo%fmat)
IF (ASSOCIATED(rhoz_cneo%wfn)) &
DEALLOCATE (rhoz_cneo%wfn)
IF (ASSOCIATED(rhoz_cneo%cpc_h)) &
DEALLOCATE (rhoz_cneo%cpc_h)
IF (ASSOCIATED(rhoz_cneo%cpc_s)) &
DEALLOCATE (rhoz_cneo%cpc_s)
IF (ASSOCIATED(rhoz_cneo%rho_rad_h)) &
DEALLOCATE (rhoz_cneo%rho_rad_h)
IF (ASSOCIATED(rhoz_cneo%rho_rad_s)) &
DEALLOCATE (rhoz_cneo%rho_rad_s)
IF (ASSOCIATED(rhoz_cneo%vrho_rad_h)) &
DEALLOCATE (rhoz_cneo%vrho_rad_h)
IF (ASSOCIATED(rhoz_cneo%vrho_rad_s)) &
DEALLOCATE (rhoz_cneo%vrho_rad_s)
IF (ASSOCIATED(rhoz_cneo%ga_Vlocal_gb_h)) &
DEALLOCATE (rhoz_cneo%ga_Vlocal_gb_h)
IF (ASSOCIATED(rhoz_cneo%ga_Vlocal_gb_s)) &
DEALLOCATE (rhoz_cneo%ga_Vlocal_gb_s)
END IF
END SUBROUTINE deallocate_rhoz_cneo
! **************************************************************************************************
!> \brief ...
!> \param rhoz_cneo_set ...
! **************************************************************************************************
SUBROUTINE deallocate_rhoz_cneo_set(rhoz_cneo_set)
TYPE(rhoz_cneo_type), DIMENSION(:), POINTER :: rhoz_cneo_set
INTEGER :: iat, natom
TYPE(rhoz_cneo_type), POINTER :: rhoz_cneo
IF (ASSOCIATED(rhoz_cneo_set)) THEN
natom = SIZE(rhoz_cneo_set)
DO iat = 1, natom
rhoz_cneo => rhoz_cneo_set(iat)
CALL deallocate_rhoz_cneo(rhoz_cneo)
END DO
DEALLOCATE (rhoz_cneo_set)
END IF
END SUBROUTINE deallocate_rhoz_cneo_set
! **************************************************************************************************
!> \brief ...
!> \param potential ...
! **************************************************************************************************
SUBROUTINE allocate_cneo_potential(potential)
TYPE(cneo_potential_type), POINTER :: potential
IF (ASSOCIATED(potential)) &
CALL deallocate_cneo_potential(potential)
ALLOCATE (potential)
END SUBROUTINE allocate_cneo_potential
! **************************************************************************************************
!> \brief ...
!> \param potential ...
! **************************************************************************************************
SUBROUTINE deallocate_cneo_potential(potential)
TYPE(cneo_potential_type), POINTER :: potential
IF (ASSOCIATED(potential)) THEN
IF (ASSOCIATED(potential%elec_conf)) &
DEALLOCATE (potential%elec_conf)
IF (ASSOCIATED(potential%my_gcc_h)) &
DEALLOCATE (potential%my_gcc_h)
IF (ASSOCIATED(potential%my_gcc_s)) &
DEALLOCATE (potential%my_gcc_s)
IF (ASSOCIATED(potential%ovlp)) &
DEALLOCATE (potential%ovlp)
IF (ASSOCIATED(potential%kin)) &
DEALLOCATE (potential%kin)
IF (ASSOCIATED(potential%utrans)) &
DEALLOCATE (potential%utrans)
IF (ASSOCIATED(potential%distance)) &
DEALLOCATE (potential%distance)
IF (ASSOCIATED(potential%harmonics)) &
CALL deallocate_harmonics_atom(potential%harmonics)
IF (ASSOCIATED(potential%Qlm_gg)) &
DEALLOCATE (potential%Qlm_gg)
IF (ASSOCIATED(potential%gg)) &
DEALLOCATE (potential%gg)
IF (ASSOCIATED(potential%vgg)) &
DEALLOCATE (potential%vgg)
IF (ASSOCIATED(potential%n2oindex)) &
DEALLOCATE (potential%n2oindex)
IF (ASSOCIATED(potential%o2nindex)) &
DEALLOCATE (potential%o2nindex)
IF (ASSOCIATED(potential%rad2l)) &
DEALLOCATE (potential%rad2l)
IF (ASSOCIATED(potential%oorad2l)) &
DEALLOCATE (potential%oorad2l)
DEALLOCATE (potential)
END IF
END SUBROUTINE deallocate_cneo_potential
! **************************************************************************************************
!> \brief ...
!> \param potential ...
!> \param z ...
!> \param zeff ...
!> \param mass ...
!> \param elec_conf ...
!> \param nsgf ...
!> \param nne ...
!> \param npsgf ...
!> \param nsotot ...
!> \param my_gcc_h ...
!> \param my_gcc_s ...
!> \param ovlp ...
!> \param kin ...
!> \param utrans ...
!> \param distance ...
!> \param harmonics ...
!> \param Qlm_gg ...
!> \param gg ...
!> \param vgg ...
!> \param n2oindex ...
!> \param o2nindex ...
!> \param rad2l ...
!> \param oorad2l ...
! **************************************************************************************************
SUBROUTINE get_cneo_potential(potential, z, zeff, mass, elec_conf, nsgf, nne, npsgf, &
nsotot, my_gcc_h, my_gcc_s, ovlp, kin, utrans, distance, &
harmonics, Qlm_gg, gg, vgg, n2oindex, o2nindex, rad2l, oorad2l)
TYPE(cneo_potential_type), POINTER :: potential
INTEGER, INTENT(OUT), OPTIONAL :: z
REAL(dp), INTENT(OUT), OPTIONAL :: zeff, mass
INTEGER, DIMENSION(:), OPTIONAL, POINTER :: elec_conf
INTEGER, INTENT(OUT), OPTIONAL :: nsgf, nne, npsgf, nsotot
REAL(dp), DIMENSION(:, :), OPTIONAL, POINTER :: my_gcc_h, my_gcc_s, ovlp, kin, utrans
REAL(dp), DIMENSION(:, :, :), OPTIONAL, POINTER :: distance
TYPE(harmonics_atom_type), OPTIONAL, POINTER :: harmonics
REAL(dp), DIMENSION(:, :, :), OPTIONAL, POINTER :: Qlm_gg, gg
REAL(dp), DIMENSION(:, :, :, :), OPTIONAL, POINTER :: vgg
INTEGER, DIMENSION(:), OPTIONAL, POINTER :: n2oindex, o2nindex
REAL(dp), DIMENSION(:, :), OPTIONAL, POINTER :: rad2l, oorad2l
IF (ASSOCIATED(potential)) THEN
IF (PRESENT(z)) z = potential%z
IF (PRESENT(zeff)) zeff = potential%zeff
IF (PRESENT(mass)) mass = potential%mass
IF (PRESENT(elec_conf)) elec_conf => potential%elec_conf
IF (PRESENT(nsgf)) nsgf = potential%nsgf
IF (PRESENT(nne)) nne = potential%nne
IF (PRESENT(npsgf)) npsgf = potential%npsgf
IF (PRESENT(nsotot)) nsotot = potential%nsotot
IF (PRESENT(my_gcc_h)) my_gcc_h => potential%my_gcc_h
IF (PRESENT(my_gcc_s)) my_gcc_s => potential%my_gcc_s
IF (PRESENT(ovlp)) ovlp => potential%ovlp
IF (PRESENT(kin)) kin => potential%kin
IF (PRESENT(ovlp)) ovlp => potential%ovlp
IF (PRESENT(utrans)) utrans => potential%utrans
IF (PRESENT(distance)) distance => potential%distance
IF (PRESENT(harmonics)) harmonics => potential%harmonics
IF (PRESENT(Qlm_gg)) Qlm_gg => potential%Qlm_gg
IF (PRESENT(gg)) gg => potential%gg
IF (PRESENT(vgg)) vgg => potential%vgg
IF (PRESENT(n2oindex)) n2oindex => potential%n2oindex
IF (PRESENT(o2nindex)) o2nindex => potential%o2nindex
IF (PRESENT(rad2l)) rad2l => potential%rad2l
IF (PRESENT(oorad2l)) oorad2l => potential%oorad2l
ELSE
CPABORT("The pointer potential is not associated.")
END IF
END SUBROUTINE get_cneo_potential
! **************************************************************************************************
!> \brief ...
!> \param potential ...
!> \param z ...
!> \param mass ...
!> \param elec_conf ...
!> \param nsgf ...
!> \param nne ...
!> \param npsgf ...
!> \param nsotot ...
!> \param my_gcc_h ...
!> \param my_gcc_s ...
!> \param ovlp ...
!> \param kin ...
!> \param utrans ...
!> \param distance ...
!> \param harmonics ...
!> \param Qlm_gg ...
!> \param gg ...
!> \param vgg ...
!> \param n2oindex ...
!> \param o2nindex ...
!> \param rad2l ...
!> \param oorad2l ...
! **************************************************************************************************
SUBROUTINE set_cneo_potential(potential, z, mass, elec_conf, nsgf, nne, npsgf, &
nsotot, my_gcc_h, my_gcc_s, ovlp, kin, utrans, distance, &
harmonics, Qlm_gg, gg, vgg, n2oindex, o2nindex, rad2l, oorad2l)
TYPE(cneo_potential_type), POINTER :: potential
INTEGER, INTENT(IN), OPTIONAL :: z
REAL(dp), INTENT(IN), OPTIONAL :: mass
INTEGER, DIMENSION(:), OPTIONAL, POINTER :: elec_conf
INTEGER, INTENT(IN), OPTIONAL :: nsgf, nne, npsgf, nsotot
REAL(dp), DIMENSION(:, :), OPTIONAL, POINTER :: my_gcc_h, my_gcc_s, ovlp, kin, utrans
REAL(dp), DIMENSION(:, :, :), OPTIONAL, POINTER :: distance
TYPE(harmonics_atom_type), OPTIONAL, POINTER :: harmonics
REAL(dp), DIMENSION(:, :, :), OPTIONAL, POINTER :: Qlm_gg, gg
REAL(dp), DIMENSION(:, :, :, :), OPTIONAL, POINTER :: vgg
INTEGER, DIMENSION(:), OPTIONAL, POINTER :: n2oindex, o2nindex
REAL(dp), DIMENSION(:, :), OPTIONAL, POINTER :: rad2l, oorad2l
IF (ASSOCIATED(potential)) THEN
IF (PRESENT(z)) THEN
potential%z = z
potential%zeff = REAL(z, dp)
IF (ASSOCIATED(potential%elec_conf)) &
CPABORT("elec_conf is already associated")
ALLOCATE (potential%elec_conf(0:3))
potential%elec_conf(0:3) = ptable(z)%e_conv(0:3)
CPASSERT(potential%mass == 0.0_dp)
IF (z == 1) THEN
! Hydrogen is 1.007825, not 1.00794
! subtract the electron mass to get the proton mass
potential%mass = 1.007825_dp - 0.000548579909_dp
ELSE
! In principle, the most abundant pure isotope mass
! should be used, but no such data is available in ptable
potential%mass = ptable(z)%amass - 0.000548579909_dp*REAL(z, dp)
END IF
END IF
IF (PRESENT(mass)) THEN
potential%mass = mass
END IF
IF (PRESENT(elec_conf)) THEN
IF (ASSOCIATED(potential%elec_conf)) THEN
DEALLOCATE (potential%elec_conf)
END IF
ALLOCATE (potential%elec_conf(0:SIZE(elec_conf) - 1))
potential%elec_conf(:) = elec_conf(:)
END IF
IF (PRESENT(nsgf)) potential%nsgf = nsgf
IF (PRESENT(nne)) potential%nne = nne
IF (PRESENT(npsgf)) potential%npsgf = npsgf
IF (PRESENT(nsotot)) potential%nsotot = nsotot
IF (PRESENT(my_gcc_h)) potential%my_gcc_h => my_gcc_h
IF (PRESENT(my_gcc_s)) potential%my_gcc_s => my_gcc_s
IF (PRESENT(ovlp)) potential%ovlp => ovlp
IF (PRESENT(kin)) potential%kin => kin
IF (PRESENT(utrans)) potential%utrans => utrans
IF (PRESENT(distance)) potential%distance => distance
IF (PRESENT(harmonics)) potential%harmonics => harmonics
IF (PRESENT(Qlm_gg)) potential%Qlm_gg => Qlm_gg
IF (PRESENT(gg)) potential%gg => gg
IF (PRESENT(vgg)) potential%vgg => vgg
IF (PRESENT(n2oindex)) potential%n2oindex => n2oindex
IF (PRESENT(o2nindex)) potential%o2nindex => o2nindex
IF (PRESENT(rad2l)) potential%rad2l => rad2l
IF (PRESENT(oorad2l)) potential%oorad2l => oorad2l
ELSE
CPABORT("The pointer potential is not associated")
END IF
END SUBROUTINE set_cneo_potential
! **************************************************************************************************
!> \brief ...
!> \param potential ...
!> \param output_unit ...
! **************************************************************************************************
SUBROUTINE write_cneo_potential(potential, output_unit)
TYPE(cneo_potential_type), POINTER :: potential
INTEGER, INTENT(IN) :: output_unit
CHARACTER(LEN=20) :: string
IF (output_unit > 0 .AND. ASSOCIATED(potential)) THEN
WRITE (UNIT=output_unit, FMT="(/,T6,A,/)") &
"CNEO Potential information"
WRITE (UNIT=output_unit, FMT="(T8,A,T41,A,I4,A,F11.6)") &
"Description: ", "Z =", potential%z, &
", nuclear mass =", potential%mass
WRITE (UNIT=string, FMT="(5I4)") potential%elec_conf
WRITE (UNIT=output_unit, FMT="(T8,A,T61,A20)") &
"Electronic configuration (s p d ...):", &
ADJUSTR(TRIM(string))
END IF
END SUBROUTINE write_cneo_potential
END MODULE qs_cneo_types

299
src/qs_cneo_utils.F Normal file
View file

@ -0,0 +1,299 @@
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
! !
! SPDX-License-Identifier: GPL-2.0-or-later !
!--------------------------------------------------------------------------------------------------!
! **************************************************************************************************
!> \brief Utility functions for CNEO-DFT
!> (see J. Chem. Theory Comput. 2025, 21, 16, 7865–7877)
!> \par History
!> 08.2025 created [zc62]
!> \author Zehua Chen
! **************************************************************************************************
MODULE qs_cneo_utils
USE ao_util, ONLY: trace_r_AxB
USE basis_set_types, ONLY: get_gto_basis_set,&
gto_basis_set_type
USE kinds, ONLY: dp
USE memory_utilities, ONLY: reallocate
USE orbital_pointers, ONLY: indso,&
nsoset
USE qs_harmonics_atom, ONLY: get_none0_cg_list,&
harmonics_atom_type
USE spherical_harmonics, ONLY: clebsch_gordon,&
clebsch_gordon_deallocate,&
clebsch_gordon_init
#include "./base/base_uses.f90"
IMPLICIT NONE
PRIVATE
! *** Global parameters ***
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_cneo_utils'
! *** Public subroutines ***
PUBLIC :: atom_solve_cneo, cneo_gather, cneo_scatter, create_harmonics_atom_cneo, &
create_my_CG_cneo, get_maxl_CG_cneo
CONTAINS
! **************************************************************************************************
!> \brief Mostly copied from qs_rho_atom_methods::init_rho_atom
!> \param my_CG ...
!> \param lcleb ...
!> \param maxl ...
!> \param llmax ...
! **************************************************************************************************
SUBROUTINE create_my_CG_cneo(my_CG, lcleb, maxl, llmax)
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: my_CG
INTEGER, INTENT(IN) :: lcleb, maxl, llmax
INTEGER :: il, iso, iso1, iso2, l1, l1l2, l2, lc1, &
lc2, lp, m1, m2, mm, mp
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: rga
! *** allocate calculate the CG coefficients up to the maxl ***
CALL clebsch_gordon_init(lcleb)
ALLOCATE (rga(lcleb, 2))
DO lc1 = 0, maxl
DO iso1 = nsoset(lc1 - 1) + 1, nsoset(lc1)
l1 = indso(1, iso1)
m1 = indso(2, iso1)
DO lc2 = 0, maxl
DO iso2 = nsoset(lc2 - 1) + 1, nsoset(lc2)
l2 = indso(1, iso2)
m2 = indso(2, iso2)
CALL clebsch_gordon(l1, m1, l2, m2, rga)
IF (l1 + l2 > llmax) THEN
l1l2 = llmax
ELSE
l1l2 = l1 + l2
END IF
mp = m1 + m2
mm = m1 - m2
IF (m1*m2 < 0 .OR. (m1*m2 == 0 .AND. (m1 < 0 .OR. m2 < 0))) THEN
mp = -ABS(mp)
mm = -ABS(mm)
ELSE
mp = ABS(mp)
mm = ABS(mm)
END IF
DO lp = MOD(l1 + l2, 2), l1l2, 2
il = lp/2 + 1
IF (ABS(mp) <= lp) THEN
IF (mp >= 0) THEN
iso = nsoset(lp - 1) + lp + 1 + mp
ELSE
iso = nsoset(lp - 1) + lp + 1 - ABS(mp)
END IF
my_CG(iso1, iso2, iso) = rga(il, 1)
END IF
IF (mp /= mm .AND. ABS(mm) <= lp) THEN
IF (mm >= 0) THEN
iso = nsoset(lp - 1) + lp + 1 + mm
ELSE
iso = nsoset(lp - 1) + lp + 1 - ABS(mm)
END IF
my_CG(iso1, iso2, iso) = rga(il, 2)
END IF
END DO
END DO ! iso2
END DO ! lc2
END DO ! iso1
END DO ! lc1
DEALLOCATE (rga)
CALL clebsch_gordon_deallocate()
END SUBROUTINE create_my_CG_cneo
! **************************************************************************************************
!> \brief Mostly copied from qs_harmonics_atom::create_harmonics_atom
!> \param harmonics ...
!> \param my_CG ...
!> \param llmax ...
!> \param maxs ...
!> \param max_s_harm ...
! **************************************************************************************************
SUBROUTINE create_harmonics_atom_cneo(harmonics, my_CG, llmax, maxs, max_s_harm)
TYPE(harmonics_atom_type), POINTER :: harmonics
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: my_CG
INTEGER, INTENT(IN) :: llmax, maxs, max_s_harm
INTEGER :: i, is
CPASSERT(ASSOCIATED(harmonics))
harmonics%max_s_harm = max_s_harm
harmonics%llmax = llmax
NULLIFY (harmonics%my_CG, harmonics%my_CG_dxyz, harmonics%my_CG_dxyz_asym)
CALL reallocate(harmonics%my_CG, 1, maxs, 1, maxs, 1, max_s_harm)
DO i = 1, max_s_harm
DO is = 1, maxs
harmonics%my_CG(1:maxs, is, i) = my_CG(1:maxs, is, i)
END DO
END DO
END SUBROUTINE create_harmonics_atom_cneo
! **************************************************************************************************
!> \brief Mostly copied from qs_harmonics_atom::get_maxl_CG
!> \param harmonics ...
!> \param orb_basis ...
!> \param llmax ...
!> \param max_s_harm ...
! **************************************************************************************************
SUBROUTINE get_maxl_CG_cneo(harmonics, orb_basis, llmax, max_s_harm)
TYPE(harmonics_atom_type), POINTER :: harmonics
TYPE(gto_basis_set_type), POINTER :: orb_basis
INTEGER, INTENT(IN) :: llmax, max_s_harm
INTEGER :: is1, is2, itmp, max_iso_not0, nset
INTEGER, DIMENSION(:), POINTER :: lmax, lmin
CPASSERT(ASSOCIATED(harmonics))
CALL get_gto_basis_set(gto_basis_set=orb_basis, lmax=lmax, lmin=lmin, nset=nset)
! *** Assign indices for the non null CG coefficients ***
max_iso_not0 = 0
DO is1 = 1, nset
DO is2 = 1, nset
CALL get_none0_cg_list(harmonics%my_CG, &
lmin(is1), lmax(is1), lmin(is2), lmax(is2), &
max_s_harm, llmax, max_iso_not0=itmp)
max_iso_not0 = MAX(max_iso_not0, itmp)
END DO ! is2
END DO ! is1
harmonics%max_iso_not0 = max_iso_not0
END SUBROUTINE get_maxl_CG_cneo
! **************************************************************************************************
!> \brief Mostly copied from atom_utils::atom_solve
!> \param hmat ...
!> \param f ...
!> \param umat ...
!> \param orb ...
!> \param ener ...
!> \param pmat ...
!> \param r ...
!> \param dist ...
!> \param nb ...
!> \param nv ...
! **************************************************************************************************
SUBROUTINE atom_solve_cneo(hmat, f, umat, orb, ener, pmat, r, dist, nb, nv)
REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: hmat
REAL(KIND=dp), DIMENSION(3), INTENT(IN) :: f
REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: umat
REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: orb
REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: ener
REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: pmat
REAL(KIND=dp), DIMENSION(3), INTENT(INOUT) :: r
REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: dist
INTEGER, INTENT(IN) :: nb, nv
INTEGER :: info, lwork, m, n
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: w, work
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: a, b, h_fx
CPASSERT(nb >= nv)
orb = 0._dp
n = nb
m = nv
IF (n > 0 .AND. m > 0) THEN
lwork = MAX(n*n, n + 100)
ALLOCATE (a(m, m), b(n, m), w(m), work(lwork))
IF (DOT_PRODUCT(f, f) /= 0.0_dp) THEN
ALLOCATE (h_fx(n, n))
h_fx(1:n, 1:n) = hmat(1:n, 1:n) + f(1)*dist(1:n, 1:n, 1) + &
f(2)*dist(1:n, 1:n, 2) + f(3)*dist(1:n, 1:n, 3)
CALL dgemm("N", "N", n, m, n, 1.0_dp, h_fx, n, umat, n, 0.0_dp, b, n)
DEALLOCATE (h_fx)
ELSE
CALL dgemm("N", "N", n, m, n, 1.0_dp, hmat, n, umat, n, 0.0_dp, b, n)
END IF
CALL dgemm("T", "N", m, m, n, 1.0_dp, umat, n, b, n, 0.0_dp, a, m)
CALL dsyev("V", "U", m, a, m, w, work, lwork, info)
CALL dgemm("N", "N", n, m, m, 1.0_dp, umat, n, a, m, 0.0_dp, b, n)
m = MIN(m, SIZE(orb, 2))
orb(1:n, 1:m) = b(1:n, 1:m)
ener(1:m) = w(1:m)
DEALLOCATE (a, b, w, work)
! calculate the density matrix using the orbital with the lowest orbital energy
pmat = 0.0_dp
CALL dger(n, n, 1.0_dp, orb(:, 1), 1, orb(:, 1), 1, pmat, n)
! calculate the expectation position (basis center as the origin)
r = (/trace_r_AxB(dist(1:n, 1:n, 1), n, pmat, n, n, n), &
trace_r_AxB(dist(1:n, 1:n, 2), n, pmat, n, n, n), &
trace_r_AxB(dist(1:n, 1:n, 3), n, pmat, n, n, n)/)
END IF
END SUBROUTINE atom_solve_cneo
! **************************************************************************************************
!> \brief Mostly copied from qs_oce_methods::prj_gather
!> \param ain ...
!> \param aout ...
!> \param nbas ...
!> \param n2oindex ...
! **************************************************************************************************
SUBROUTINE cneo_gather(ain, aout, nbas, n2oindex)
REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: ain
REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: aout
INTEGER, INTENT(IN) :: nbas
INTEGER, DIMENSION(:), POINTER :: n2oindex
INTEGER :: i, ip, j, jp
DO i = 1, nbas
ip = n2oindex(i)
DO j = 1, nbas
jp = n2oindex(j)
aout(j, i) = ain(jp, ip)
END DO
END DO
END SUBROUTINE cneo_gather
! **************************************************************************************************
!> \brief Mostly copied from qs_oce_methods::prj_scatter
!> \param ain ...
!> \param aout ...
!> \param nbas ...
!> \param n2oindex ...
! **************************************************************************************************
SUBROUTINE cneo_scatter(ain, aout, nbas, n2oindex)
REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: ain
REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: aout
INTEGER, INTENT(IN) :: nbas
INTEGER, DIMENSION(:), POINTER :: n2oindex
INTEGER :: i, ip, j, jp
DO i = 1, nbas
ip = n2oindex(i)
DO j = 1, nbas
jp = n2oindex(j)
aout(jp, ip) = aout(jp, ip) + ain(j, i)
END DO
END DO
END SUBROUTINE cneo_scatter
END MODULE qs_cneo_utils

6
src/qs_collocate_density.F
View file

@ -852,8 +852,10 @@ CONTAINS
alpha = calpha(ikind)
pab(1, 1) = ccore(ikind)
ELSE
CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom, &
alpha_core_charge=alpha, ccore_charge=pab(1, 1))
CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom)
IF (my_only_nopaw .AND. paw_atom) CYCLE
CALL get_qs_kind(qs_kind_set(ikind), alpha_core_charge=alpha, &
ccore_charge=pab(1, 1))
END IF
IF (my_only_nopaw .AND. paw_atom) CYCLE

28
src/qs_core_energies.F
View file

@ -29,6 +29,7 @@ MODULE qs_core_energies
USE message_passing, ONLY: mp_comm_type,&
mp_para_env_type
USE particle_types, ONLY: particle_type
USE qs_cneo_types, ONLY: cneo_potential_type
USE qs_energy_types, ONLY: qs_energy_type
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
@ -216,6 +217,7 @@ CONTAINS
REAL(KIND=dp), DIMENSION(3, 3) :: pv_loc
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(atprop_type), POINTER :: atprop
TYPE(cneo_potential_type), POINTER :: cneo_potential
TYPE(mp_comm_type) :: group
TYPE(neighbor_list_iterator_p_type), &
DIMENSION(:), POINTER :: nl_iterator
@ -275,10 +277,19 @@ CONTAINS
END IF
DO ikind = 1, nkind
CALL get_qs_kind(qs_kind_set(ikind), &
alpha_core_charge=alpha(ikind), &
core_charge_radius=radius(ikind), &
zeff=zeff(ikind))
! cneo quantum nuclei have their core energies calculated elsewhere
NULLIFY (cneo_potential)
CALL get_qs_kind(qs_kind_set(ikind), cneo_potential=cneo_potential)
IF (ASSOCIATED(cneo_potential)) THEN
alpha(ikind) = 1.0_dp
radius(ikind) = 1.0_dp
zeff(ikind) = 0.0_dp
ELSE
CALL get_qs_kind(qs_kind_set(ikind), &
alpha_core_charge=alpha(ikind), &
core_charge_radius=radius(ikind), &
zeff=zeff(ikind))
END IF
END DO
ecore_overlap = 0.0_dp
@ -365,6 +376,7 @@ CONTAINS
REAL(KIND=dp) :: alpha_core_charge, ecore_self, es, zeff
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(atprop_type), POINTER :: atprop
TYPE(cneo_potential_type), POINTER :: cneo_potential
TYPE(distribution_1d_type), POINTER :: local_particles
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(qs_energy_type), POINTER :: energy
@ -381,6 +393,10 @@ CONTAINS
ecore_self = 0.0_dp
DO ikind = 1, SIZE(atomic_kind_set)
! nuclear density self-interaction is already removed in CNEO
NULLIFY (cneo_potential)
CALL get_qs_kind(qs_kind_set(ikind), cneo_potential=cneo_potential)
IF (ASSOCIATED(cneo_potential)) CYCLE
CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom)
CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff, alpha_core_charge=alpha_core_charge)
ecore_self = ecore_self - REAL(natom, dp)*zeff**2*SQRT(alpha_core_charge)
@ -400,6 +416,10 @@ CONTAINS
CALL atprop_array_init(atprop%ateself, natom)
DO ikind = 1, SIZE(atomic_kind_set)
! nuclear density self-interaction is already removed in CNEO
NULLIFY (cneo_potential)
CALL get_qs_kind(qs_kind_set(ikind), cneo_potential=cneo_potential)
IF (ASSOCIATED(cneo_potential)) CYCLE
nparticle_local = local_particles%n_el(ikind)
CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff, alpha_core_charge=alpha_core_charge)
es = zeff**2*SQRT(alpha_core_charge)/SQRT(twopi)

4
src/qs_core_hamiltonian.F
View file

@ -106,6 +106,7 @@ MODULE qs_core_hamiltonian
dp
USE message_passing, ONLY: mp_para_env_type
USE particle_types, ONLY: particle_type
USE qs_cneo_methods, ONLY: cneo_core_matrices
USE qs_condnum, ONLY: overlap_condnum
USE qs_core_matrices, ONLY: core_matrices,&
kinetic_energy_matrix
@ -322,6 +323,9 @@ CONTAINS
! *** core and pseudopotentials
CALL core_matrices(qs_env, matrix_h, matrix_p, calculate_forces, nder)
! *** CNEO nuclear V_core
CALL cneo_core_matrices(qs_env, calculate_forces, nder)
! *** GAPW one-center-expansion (oce) matrices
NULLIFY (sap_oce)
CALL get_qs_env(qs_env=qs_env, sap_oce=sap_oce)

2
src/qs_energy_types.F
View file

@ -27,6 +27,7 @@ MODULE qs_energy_types
core_overlap = 0.0_dp, &
core_overlap0 = 0.0_dp, &
core_self = 0.0_dp, &
core_cneo = 0.0_dp, & ! quantum nuclear kinetic + v_core
repulsive = 0.0_dp, &
dispersion = 0.0_dp, &
dispersion_sc = 0.0_dp, &
@ -160,6 +161,7 @@ CONTAINS
qs_energy%core_overlap = 0.0_dp
qs_energy%core_overlap0 = 0.0_dp
qs_energy%core_self = 0.0_dp
qs_energy%core_cneo = 0.0_dp
qs_energy%repulsive = 0.0_dp
qs_energy%dispersion = 0.0_dp
qs_energy%gcp = 0.0_dp

65
src/qs_environment.F
View file

@ -164,8 +164,9 @@ MODULE qs_environment
write_ppl_radii,&
write_ppnl_radii
USE qs_kind_types, ONLY: &
check_qs_kind_set, get_qs_kind, get_qs_kind_set, init_gapw_basis_set, init_gapw_nlcc, &
init_qs_kind_set, qs_kind_type, set_qs_kind, write_gto_basis_sets, write_qs_kind_set
check_qs_kind_set, get_qs_kind, get_qs_kind_set, init_cneo_basis_set, init_gapw_basis_set, &
init_gapw_nlcc, init_qs_kind_set, qs_kind_type, set_qs_kind, write_gto_basis_sets, &
write_qs_kind_set
USE qs_ks_types, ONLY: qs_ks_env_create,&
qs_ks_env_type,&
set_ks_env
@ -598,13 +599,13 @@ CONTAINS
CHARACTER(len=2) :: element_symbol
INTEGER :: gfn_type, handle, ikind, ispin, iw, lmax_sphere, maxl, maxlgto, maxlgto_lri, &
maxlppl, maxlppnl, method_id, multiplicity, my_ival, n_ao, n_mo_add, natom, nelectron, &
ngauss, nkind, output_unit, sort_basis, tnadd_method
maxlgto_nuc, maxlppl, maxlppnl, method_id, multiplicity, my_ival, n_ao, n_mo_add, natom, &
nelectron, ngauss, nkind, output_unit, sort_basis, tnadd_method
INTEGER, DIMENSION(2) :: n_mo, nelectron_spin
INTEGER, DIMENSION(5) :: occ
LOGICAL :: all_potential_present, be_silent, do_kpoints, do_ri_hfx, do_ri_mp2, do_ri_rpa, &
do_ri_sos_mp2, do_rpa_ri_exx, do_wfc_im_time, e1terms, has_unit_metric, lribas, &
mp2_present, orb_gradient
LOGICAL :: all_potential_present, be_silent, cneo_potential_present, do_kpoints, do_ri_hfx, &
do_ri_mp2, do_ri_rpa, do_ri_sos_mp2, do_rpa_ri_exx, do_wfc_im_time, e1terms, &
has_unit_metric, lribas, mp2_present, orb_gradient
REAL(KIND=dp) :: alpha, ccore, ewald_rcut, fxx, maxocc, &
rcut, total_zeff_corr, verlet_skin, &
zeff_correction
@ -907,6 +908,13 @@ CONTAINS
END IF
END IF
! *** Check that no cneo potential is present if not GAPW
CALL get_qs_kind_set(qs_kind_set, cneo_potential_present=cneo_potential_present)
IF (cneo_potential_present .AND. &
dft_control%qs_control%method_id /= do_method_gapw) THEN
CPABORT("CNEO calculations require GAPW method")
END IF
! DFT+U
CALL get_qs_kind_set(qs_kind_set, dft_plus_u_atom_present=dft_control%dft_plus_u)
@ -1081,6 +1089,11 @@ CONTAINS
! *** the orbital transformation matrices ***
CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto, maxlppl=maxlppl, maxlppnl=maxlppnl)
! CNEO nuclear basis contributes to GAPW rho0
IF (cneo_potential_present) THEN
CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto_nuc, basis_type="NUC")
maxlgto = MAX(maxlgto, maxlgto_nuc)
END IF
lmax_sphere = dft_control%qs_control%gapw_control%lmax_sphere
IF (lmax_sphere .LT. 0) THEN
lmax_sphere = 2*maxlgto
@ -1115,6 +1128,11 @@ CONTAINS
CALL init_gapw_basis_set(qs_kind_set, qs_control, qs_env%input)
END IF
! *** Initialise CNEO nuclear soft basis
IF (cneo_potential_present) THEN
CALL init_cneo_basis_set(qs_kind_set, qs_control)
END IF
! *** Initialize the pretabulation for the calculation of the ***
! *** incomplete Gamma function F_n(t) after McMurchie-Davidson ***
CALL get_qs_kind_set(qs_kind_set, maxlgto=maxlgto)
@ -1147,6 +1165,7 @@ CONTAINS
CALL write_pgf_orb_radii("orb", atomic_kind_set, qs_kind_set, subsys_section)
CALL write_pgf_orb_radii("aux", atomic_kind_set, qs_kind_set, subsys_section)
CALL write_pgf_orb_radii("lri", atomic_kind_set, qs_kind_set, subsys_section)
CALL write_pgf_orb_radii("nuc", atomic_kind_set, qs_kind_set, subsys_section)
CALL write_core_charge_radii(atomic_kind_set, qs_kind_set, subsys_section)
CALL write_ppl_radii(atomic_kind_set, qs_kind_set, subsys_section)
CALL write_ppnl_radii(atomic_kind_set, qs_kind_set, subsys_section)
@ -1867,9 +1886,9 @@ CONTAINS
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(section_vals_type), POINTER :: force_env_section
INTEGER :: maxlgto, maxlppl, maxlppnl, natom, ncgf, &
nkind, npgf, nset, nsgf, nshell, &
output_unit
INTEGER :: maxlgto, maxlppl, maxlppnl, natom, &
natom_q, ncgf, nkind, nkind_q, npgf, &
nset, nsgf, nshell, output_unit
TYPE(cp_logger_type), POINTER :: logger
NULLIFY (logger)
@ -2000,6 +2019,32 @@ CONTAINS
" Maximum angular momentum ", maxlgto
END IF
! NUCLEAR BASIS
CALL get_qs_kind_set(qs_kind_set, &
nkind_q=nkind_q, &
natom_q=natom_q, &
maxlgto=maxlgto, &
ncgf=ncgf, &
npgf=npgf, &
nset=nset, &
nsgf=nsgf, &
nshell=nshell, &
basis_type="NUC")
IF (nset + npgf + ncgf > 0) THEN
WRITE (UNIT=output_unit, FMT="(/,T3,A,/,T3,A,(T30,A,T71,I10))") &
"Nuclear Basis: ", &
"Total number of", &
"- Quantum atomic kinds: ", nkind_q, &
"- Quantum atoms: ", natom_q, &
"- Shell sets: ", nset, &
"- Shells: ", nshell, &
"- Primitive Cartesian functions: ", npgf, &
"- Cartesian basis functions: ", ncgf, &
"- Spherical basis functions: ", nsgf
WRITE (UNIT=output_unit, FMT="(T30,A,T75,I6)") &
" Maximum angular momentum ", maxlgto
END IF
END IF
CALL cp_print_key_finished_output(output_unit, logger, force_env_section, &
"PRINT%TOTAL_NUMBERS")

41
src/qs_environment_types.F
View file

@ -90,6 +90,7 @@ MODULE qs_environment_types
release_active_space_type
USE qs_charges_types, ONLY: qs_charges_release,&
qs_charges_type
USE qs_cneo_types, ONLY: rhoz_cneo_type
USE qs_dftb_types, ONLY: qs_dftb_pairpot_release,&
qs_dftb_pairpot_type
USE qs_dispersion_types, ONLY: qs_dispersion_release,&
@ -358,6 +359,7 @@ CONTAINS
!> \param sab_almo ...
!> \param sab_kp ...
!> \param sab_kp_nosym ...
!> \param sab_cneo ...
!> \param particle_set ...
!> \param energy ...
!> \param force ...
@ -420,9 +422,12 @@ CONTAINS
!> \param rho0_atom_set ...
!> \param rho0_mpole ...
!> \param rhoz_set ...
!> \param rhoz_cneo_set ...
!> \param ecoul_1c ...
!> \param rho0_s_rs ...
!> \param rho0_s_gs ...
!> \param rhoz_cneo_s_rs ...
!> \param rhoz_cneo_s_gs ...
!> \param do_kpoints ...
!> \param has_unit_metric ...
!> \param requires_mo_derivs ...
@ -494,7 +499,7 @@ CONTAINS
dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, &
sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, &
sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, &
sab_kp, sab_kp_nosym, particle_set, energy, force, &
sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, &
matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, &
matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, &
matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, &
@ -505,8 +510,9 @@ CONTAINS
local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, &
molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, &
task_list, task_list_soft, &
rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, &
rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, &
rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, &
rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, &
do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, &
mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, &
neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, &
outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, &
@ -533,7 +539,7 @@ CONTAINS
LOGICAL, OPTIONAL :: qmmm, qmmm_periodic
TYPE(neighbor_list_set_p_type), DIMENSION(:), OPTIONAL, POINTER :: sac_ae, sac_ppl, sac_lri, &
sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, &
sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym
sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo
TYPE(particle_type), DIMENSION(:), OPTIONAL, &
POINTER :: particle_set
TYPE(qs_energy_type), OPTIONAL, POINTER :: energy
@ -590,10 +596,14 @@ CONTAINS
POINTER :: rho0_atom_set
TYPE(rho0_mpole_type), OPTIONAL, POINTER :: rho0_mpole
TYPE(rhoz_type), DIMENSION(:), OPTIONAL, POINTER :: rhoz_set
TYPE(rhoz_cneo_type), DIMENSION(:), OPTIONAL, &
POINTER :: rhoz_cneo_set
TYPE(ecoul_1center_type), DIMENSION(:), OPTIONAL, &
POINTER :: ecoul_1c
TYPE(pw_r3d_rs_type), OPTIONAL, POINTER :: rho0_s_rs
TYPE(pw_c1d_gs_type), OPTIONAL, POINTER :: rho0_s_gs
TYPE(pw_r3d_rs_type), OPTIONAL, POINTER :: rhoz_cneo_s_rs
TYPE(pw_c1d_gs_type), OPTIONAL, POINTER :: rhoz_cneo_s_gs
LOGICAL, OPTIONAL :: do_kpoints, has_unit_metric, &
requires_mo_derivs
TYPE(dbcsr_p_type), DIMENSION(:), OPTIONAL, &
@ -736,6 +746,8 @@ CONTAINS
CALL get_local_rho(qs_env%local_rho_set, rho0_mpole=rho0_mpole)
IF (PRESENT(rhoz_set)) &
CALL get_local_rho(qs_env%local_rho_set, rhoz_set=rhoz_set)
IF (PRESENT(rhoz_cneo_set)) &
CALL get_local_rho(qs_env%local_rho_set, rhoz_cneo_set=rhoz_cneo_set)
IF (PRESENT(ecoul_1c)) &
CALL get_hartree_local(qs_env%hartree_local, ecoul_1c=ecoul_1c)
IF (PRESENT(rho0_s_rs)) THEN
@ -750,6 +762,18 @@ CONTAINS
rho0_s_gs => rho0_m%rho0_s_gs
END IF
END IF
IF (PRESENT(rhoz_cneo_s_rs)) THEN
CALL get_local_rho(qs_env%local_rho_set, rho0_mpole=rho0_m)
IF (ASSOCIATED(rho0_m)) THEN
rhoz_cneo_s_rs => rho0_m%rhoz_cneo_s_rs
END IF
END IF
IF (PRESENT(rhoz_cneo_s_gs)) THEN
CALL get_local_rho(qs_env%local_rho_set, rho0_mpole=rho0_m)
IF (ASSOCIATED(rho0_m)) THEN
rhoz_cneo_s_gs => rho0_m%rhoz_cneo_s_gs
END IF
END IF
IF (PRESENT(xas_env)) xas_env => qs_env%xas_env
IF (PRESENT(input)) input => qs_env%input
@ -815,6 +839,7 @@ CONTAINS
sab_almo=sab_almo, &
sab_kp=sab_kp, &
sab_kp_nosym=sab_kp_nosym, &
sab_cneo=sab_cneo, &
task_list=task_list, &
task_list_soft=task_list_soft, &
kpoints=kpoints, &
@ -1008,6 +1033,7 @@ CONTAINS
!> \param mo_derivs ...
!> \param mo_loc_history ...
!> \param efield ...
!> \param rhoz_cneo_set ...
!> \param linres_control ...
!> \param xas_env ...
!> \param cp_ddapc_env ...
@ -1061,7 +1087,7 @@ CONTAINS
ks_qmmm_env, wf_history, scf_env, active_space, &
input, oce, rho_atom_set, rho0_atom_set, rho0_mpole, run_rtp, rtp, &
rhoz_set, rhoz_tot, ecoul_1c, has_unit_metric, requires_mo_derivs, mo_derivs, &
mo_loc_history, efield, &
mo_loc_history, efield, rhoz_cneo_set, &
linres_control, xas_env, cp_ddapc_env, cp_ddapc_ewald, &
outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, &
se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, ls_scf_env, &
@ -1105,6 +1131,8 @@ CONTAINS
POINTER :: mo_derivs
TYPE(cp_fm_type), DIMENSION(:), OPTIONAL, POINTER :: mo_loc_history
TYPE(efield_berry_type), OPTIONAL, POINTER :: efield
TYPE(rhoz_cneo_type), DIMENSION(:), OPTIONAL, &
POINTER :: rhoz_cneo_set
TYPE(linres_control_type), OPTIONAL, POINTER :: linres_control
TYPE(xas_environment_type), OPTIONAL, POINTER :: xas_env
TYPE(cp_ddapc_type), OPTIONAL, POINTER :: cp_ddapc_env
@ -1315,6 +1343,9 @@ CONTAINS
IF (PRESENT(rhoz_set)) THEN
CALL set_local_rho(qs_env%local_rho_set, rhoz_set=rhoz_set)
END IF
IF (PRESENT(rhoz_cneo_set)) THEN
CALL set_local_rho(qs_env%local_rho_set, rhoz_cneo_set=rhoz_cneo_set)
END IF
IF (PRESENT(rhoz_tot)) qs_env%local_rho_set%rhoz_tot = rhoz_tot
IF (PRESENT(ecoul_1c)) THEN
CALL set_hartree_local(qs_env%hartree_local, ecoul_1c=ecoul_1c)

2
src/qs_force.F
View file

@ -598,10 +598,12 @@ CONTAINS
iatom, ikind, " gth_nlcc", qs_force(ikind)%gth_nlcc(1:3, i), &
iatom, ikind, " gth_ppnl", qs_force(ikind)%gth_ppnl(1:3, i), &
iatom, ikind, " all_potential", qs_force(ikind)%all_potential(1:3, i), &
iatom, ikind, "cneo_potential", qs_force(ikind)%cneo_potential(1:3, i), &
iatom, ikind, " core_overlap", qs_force(ikind)%core_overlap(1:3, i), &
iatom, ikind, " rho_core", qs_force(ikind)%rho_core(1:3, i), &
iatom, ikind, " rho_elec", qs_force(ikind)%rho_elec(1:3, i), &
iatom, ikind, " rho_lri_elec", qs_force(ikind)%rho_lri_elec(1:3, i), &
iatom, ikind, " rho_cneo_nuc", qs_force(ikind)%rho_cneo_nuc(1:3, i), &
iatom, ikind, " vhxc_atom", qs_force(ikind)%vhxc_atom(1:3, i), &
iatom, ikind, " g0s_Vh_elec", qs_force(ikind)%g0s_Vh_elec(1:3, i), &
iatom, ikind, " ch_pulay", qs_force(ikind)%ch_pulay(1:3, i), &

24
src/qs_force_types.F
View file

@ -27,6 +27,7 @@ MODULE qs_force_types
TYPE qs_force_type
REAL(KIND=dp), DIMENSION(:, :), POINTER :: all_potential => NULL(), &
cneo_potential => NULL(), &
core_overlap => NULL(), &
gth_ppl => NULL(), &
gth_nlcc => NULL(), &
@ -37,6 +38,7 @@ MODULE qs_force_types
rho_core => NULL(), &
rho_elec => NULL(), &
rho_lri_elec => NULL(), &
rho_cneo_nuc => NULL(), &
vhxc_atom => NULL(), &
g0s_Vh_elec => NULL(), &
repulsive => NULL(), &
@ -91,6 +93,7 @@ CONTAINS
DO ikind = 1, nkind
n = natom_of_kind(ikind)
ALLOCATE (qs_force(ikind)%all_potential(3, n))
ALLOCATE (qs_force(ikind)%cneo_potential(3, n))
ALLOCATE (qs_force(ikind)%core_overlap(3, n))
ALLOCATE (qs_force(ikind)%gth_ppl(3, n))
ALLOCATE (qs_force(ikind)%gth_nlcc(3, n))
@ -101,6 +104,7 @@ CONTAINS
ALLOCATE (qs_force(ikind)%rho_core(3, n))
ALLOCATE (qs_force(ikind)%rho_elec(3, n))
ALLOCATE (qs_force(ikind)%rho_lri_elec(3, n))
ALLOCATE (qs_force(ikind)%rho_cneo_nuc(3, n))
ALLOCATE (qs_force(ikind)%vhxc_atom(3, n))
ALLOCATE (qs_force(ikind)%g0s_Vh_elec(3, n))
ALLOCATE (qs_force(ikind)%repulsive(3, n))
@ -143,6 +147,10 @@ CONTAINS
DEALLOCATE (qs_force(ikind)%all_potential)
END IF
IF (ASSOCIATED(qs_force(ikind)%cneo_potential)) THEN
DEALLOCATE (qs_force(ikind)%cneo_potential)
END IF
IF (ASSOCIATED(qs_force(ikind)%core_overlap)) THEN
DEALLOCATE (qs_force(ikind)%core_overlap)
END IF
@ -182,6 +190,10 @@ CONTAINS
DEALLOCATE (qs_force(ikind)%rho_lri_elec)
END IF
IF (ASSOCIATED(qs_force(ikind)%rho_cneo_nuc)) THEN
DEALLOCATE (qs_force(ikind)%rho_cneo_nuc)
END IF
IF (ASSOCIATED(qs_force(ikind)%vhxc_atom)) THEN
DEALLOCATE (qs_force(ikind)%vhxc_atom)
END IF
@ -256,6 +268,7 @@ CONTAINS
DO ikind = 1, SIZE(qs_force)
qs_force(ikind)%all_potential(:, :) = 0.0_dp
qs_force(ikind)%cneo_potential(:, :) = 0.0_dp
qs_force(ikind)%core_overlap(:, :) = 0.0_dp
qs_force(ikind)%gth_ppl(:, :) = 0.0_dp
qs_force(ikind)%gth_nlcc(:, :) = 0.0_dp
@ -266,6 +279,7 @@ CONTAINS
qs_force(ikind)%rho_core(:, :) = 0.0_dp
qs_force(ikind)%rho_elec(:, :) = 0.0_dp
qs_force(ikind)%rho_lri_elec(:, :) = 0.0_dp
qs_force(ikind)%rho_cneo_nuc(:, :) = 0.0_dp
qs_force(ikind)%vhxc_atom(:, :) = 0.0_dp
qs_force(ikind)%g0s_Vh_elec(:, :) = 0.0_dp
qs_force(ikind)%repulsive(:, :) = 0.0_dp
@ -300,6 +314,8 @@ CONTAINS
DO ikind = 1, SIZE(qs_force_out)
qs_force_out(ikind)%all_potential(:, :) = qs_force_out(ikind)%all_potential(:, :) + &
qs_force_in(ikind)%all_potential(:, :)
qs_force_out(ikind)%cneo_potential(:, :) = qs_force_out(ikind)%cneo_potential(:, :) + &
qs_force_in(ikind)%cneo_potential(:, :)
qs_force_out(ikind)%core_overlap(:, :) = qs_force_out(ikind)%core_overlap(:, :) + &
qs_force_in(ikind)%core_overlap(:, :)
qs_force_out(ikind)%gth_ppl(:, :) = qs_force_out(ikind)%gth_ppl(:, :) + &
@ -320,6 +336,8 @@ CONTAINS
qs_force_in(ikind)%rho_elec(:, :)
qs_force_out(ikind)%rho_lri_elec(:, :) = qs_force_out(ikind)%rho_lri_elec(:, :) + &
qs_force_in(ikind)%rho_lri_elec(:, :)
qs_force_out(ikind)%rho_cneo_nuc(:, :) = qs_force_out(ikind)%rho_cneo_nuc(:, :) + &
qs_force_in(ikind)%rho_cneo_nuc(:, :)
qs_force_out(ikind)%vhxc_atom(:, :) = qs_force_out(ikind)%vhxc_atom(:, :) + &
qs_force_in(ikind)%vhxc_atom(:, :)
qs_force_out(ikind)%g0s_Vh_elec(:, :) = qs_force_out(ikind)%g0s_Vh_elec(:, :) + &
@ -372,10 +390,12 @@ CONTAINS
CALL para_env%sum(qs_force(ikind)%gth_nlcc)
CALL para_env%sum(qs_force(ikind)%gth_ppnl)
CALL para_env%sum(qs_force(ikind)%all_potential)
CALL para_env%sum(qs_force(ikind)%cneo_potential)
CALL para_env%sum(qs_force(ikind)%core_overlap)
CALL para_env%sum(qs_force(ikind)%rho_core)
CALL para_env%sum(qs_force(ikind)%rho_elec)
CALL para_env%sum(qs_force(ikind)%rho_lri_elec)
CALL para_env%sum(qs_force(ikind)%rho_cneo_nuc)
CALL para_env%sum(qs_force(ikind)%vhxc_atom)
CALL para_env%sum(qs_force(ikind)%g0s_Vh_elec)
CALL para_env%sum(qs_force(ikind)%fock_4c)
@ -391,12 +411,14 @@ CONTAINS
qs_force(ikind)%gth_nlcc(:, :) + &
qs_force(ikind)%gth_ppnl(:, :) + &
qs_force(ikind)%all_potential(:, :) + &
qs_force(ikind)%cneo_potential(:, :) + &
qs_force(ikind)%kinetic(:, :) + &
qs_force(ikind)%overlap(:, :) + &
qs_force(ikind)%overlap_admm(:, :) + &
qs_force(ikind)%rho_core(:, :) + &
qs_force(ikind)%rho_elec(:, :) + &
qs_force(ikind)%rho_lri_elec(:, :) + &
qs_force(ikind)%rho_cneo_nuc(:, :) + &
qs_force(ikind)%vhxc_atom(:, :) + &
qs_force(ikind)%g0s_Vh_elec(:, :) + &
qs_force(ikind)%fock_4c(:, :) + &
@ -558,12 +580,14 @@ CONTAINS
qs_force(ikind)%gth_nlcc(:, ia) + &
qs_force(ikind)%gth_ppnl(:, ia) + &
qs_force(ikind)%all_potential(:, ia) + &
qs_force(ikind)%cneo_potential(:, ia) + &
qs_force(ikind)%kinetic(:, ia) + &
qs_force(ikind)%overlap(:, ia) + &
qs_force(ikind)%overlap_admm(:, ia) + &
qs_force(ikind)%rho_core(:, ia) + &
qs_force(ikind)%rho_elec(:, ia) + &
qs_force(ikind)%rho_lri_elec(:, ia) + &
qs_force(ikind)%rho_cneo_nuc(:, ia) + &
qs_force(ikind)%vhxc_atom(:, ia) + &
qs_force(ikind)%g0s_Vh_elec(:, ia) + &
qs_force(ikind)%fock_4c(:, ia) + &

34
src/qs_gapw_densities.F
View file

@ -18,6 +18,8 @@ MODULE qs_gapw_densities
pw_env_type
USE pw_pool_types, ONLY: pw_pool_p_type
USE qs_charges_types, ONLY: qs_charges_type
USE qs_cneo_ggrid, ONLY: put_rhoz_cneo_s_on_grid
USE qs_cneo_types, ONLY: rhoz_cneo_type
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
USE qs_kind_types, ONLY: get_qs_kind,&
@ -81,6 +83,7 @@ CONTAINS
TYPE(rho0_atom_type), DIMENSION(:), POINTER :: rho0_atom_set
TYPE(rho0_mpole_type), POINTER :: rho0_mpole
TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
TYPE(rhoz_cneo_type), DIMENSION(:), POINTER :: rhoz_cneo_set
CALL timeset(routineN, handle)
@ -108,13 +111,15 @@ CONTAINS
atomic_kind_set=atomic_kind_set, &
rho0_mpole=rho0_mpole, &
rho_atom_set=rho_atom_set, &
rho0_atom_set=rho0_atom_set)
rho0_atom_set=rho0_atom_set, &
rhoz_cneo_set=rhoz_cneo_set)
gapw_control => dft_control%qs_control%gapw_control
! If TDDFPT%MGRID is defined, overwrite QS grid info accordingly
IF (PRESENT(local_rho_set)) THEN
rho_atom_set => local_rho_set%rho_atom_set
rhoz_cneo_set => local_rho_set%rhoz_cneo_set
IF (my_do_rho0) THEN
rho0_mpole => local_rho_set%rho0_mpole
rho0_atom_set => local_rho_set%rho0_atom_set
@ -147,8 +152,9 @@ CONTAINS
!Calculate rho0_h and rho0_s on the radial grids centered on the atomic position
IF (my_do_rho0) &
CALL calculate_rho0_atom(gapw_control, rho_atom_set, rho0_atom_set, rho0_mpole, &
atom_list, natom, ikind, my_kind_set(ikind), rho0_h_tot)
CALL calculate_rho0_atom(gapw_control, rho_atom_set, rhoz_cneo_set, rho0_atom_set, &
rho0_mpole, atom_list, natom, ikind, my_kind_set(ikind), &
rho0_h_tot)
END DO
@ -193,6 +199,28 @@ CONTAINS
END IF
qs_charges%total_rho0_soft_rspace = tot_rs_int
qs_charges%total_rho0_hard_lebedev = rho0_h_tot
IF (rho0_mpole%do_cneo) THEN
! put soft tails of quantum nuclear charge densities on the global grid
CALL put_rhoz_cneo_s_on_grid(qs_env, rho0_mpole, rhoz_cneo_set, tot_rs_int)
IF (ABS(rho0_mpole%tot_rhoz_cneo_s) .GE. 1.0E-5_dp) THEN
IF (ABS(1.0_dp - ABS(tot_rs_int/rho0_mpole%tot_rhoz_cneo_s)) .GT. 1.0E-3_dp) THEN
IF (output_unit > 0) THEN
WRITE (output_unit, '(/,72("*"))')
WRITE (output_unit, '(T2,A,T66,1E20.8)') &
"WARNING: rhoz_cneo_s calculated on the local grid is :", &
rho0_mpole%tot_rhoz_cneo_s, &
" rhoz_cneo_s calculated on the global grid is :", tot_rs_int
WRITE (output_unit, '(T2,A)') &
" bad integration"
WRITE (output_unit, '(72("*"),/)')
END IF
END IF
END IF
qs_charges%total_rho1_soft_nuc_rspace = tot_rs_int
qs_charges%total_rho1_soft_nuc_lebedev = rho0_mpole%tot_rhoz_cneo_s
ELSE
qs_charges%total_rho1_soft_nuc_rspace = 0.0_dp
END IF
ELSE
qs_charges%total_rho0_hard_lebedev = 0.0_dp
END IF

18
src/qs_initial_guess.F
View file

@ -64,6 +64,7 @@ MODULE qs_initial_guess
USE particle_methods, ONLY: get_particle_set
USE particle_types, ONLY: particle_type
USE qs_atomic_block, ONLY: calculate_atomic_block_dm
USE qs_cneo_types, ONLY: cneo_potential_type
USE qs_density_matrices, ONLY: calculate_density_matrix
USE qs_dftb_utils, ONLY: get_dftb_atom_param
USE qs_eht_guess, ONLY: calculate_eht_guess
@ -145,8 +146,9 @@ CONTAINS
INTEGER, DIMENSION(2) :: nelectron_spin
INTEGER, DIMENSION(:), POINTER :: atom_list, elec_conf, nelec_kind, &
sort_kind
LOGICAL :: did_guess, do_hfx_ri_mo, do_kpoints, do_std_diag, exist, has_unit_metric, &
natom_mismatch, need_mos, need_wm, ofgpw, owns_ortho, print_history_log, print_log
LOGICAL :: cneo_potential_present, did_guess, do_hfx_ri_mo, do_kpoints, do_std_diag, exist, &
has_unit_metric, natom_mismatch, need_mos, need_wm, ofgpw, owns_ortho, print_history_log, &
print_log
REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: buff, buff2
REAL(dp), DIMENSION(:, :), POINTER :: pdata
REAL(KIND=dp) :: checksum, eps, length, maxocc, occ, &
@ -556,6 +558,11 @@ CONTAINS
CPABORT("calculate_first_density_matrix: core_guess not implemented for k-points")
END IF
CALL get_qs_kind_set(qs_kind_set, cneo_potential_present=cneo_potential_present)
IF (cneo_potential_present) THEN
CPABORT("calculate_first_density_matrix: core_guess not implemented for CNEO")
END IF
owns_ortho = .FALSE.
IF (.NOT. ASSOCIATED(work1)) THEN
need_wm = .TRUE.
@ -1214,6 +1221,7 @@ CONTAINS
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: econf, pdiag, sdiag
REAL(KIND=dp), DIMENSION(0:3) :: edftb
TYPE(all_potential_type), POINTER :: all_potential
TYPE(cneo_potential_type), POINTER :: cneo_potential
TYPE(dbcsr_type), POINTER :: matrix_p
TYPE(gth_potential_type), POINTER :: gth_potential
TYPE(gto_basis_set_type), POINTER :: orb_basis_set
@ -1267,8 +1275,10 @@ CONTAINS
CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, &
all_potential=all_potential, &
gth_potential=gth_potential, &
sgp_potential=sgp_potential)
has_pot = ASSOCIATED(all_potential) .OR. ASSOCIATED(gth_potential) .OR. ASSOCIATED(sgp_potential)
sgp_potential=sgp_potential, &
cneo_potential=cneo_potential)
has_pot = ASSOCIATED(all_potential) .OR. ASSOCIATED(gth_potential) .OR. &
ASSOCIATED(sgp_potential) .OR. ASSOCIATED(cneo_potential)
IF (dft_control%qs_control%dftb) THEN
CALL get_dftb_atom_param(qs_kind_set(ikind)%dftb_parameter, &

7
src/qs_integrate_potential_single.F
View file

@ -555,12 +555,13 @@ CONTAINS
DO ikind = 1, SIZE(atomic_kind_set)
CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom_of_kind, atom_list=atom_list)
CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom, &
alpha_core_charge=alpha_core_charge, &
ccore_charge=ccore_charge)
CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom)
IF (dft_control%qs_control%gapw .AND. paw_atom) CYCLE
CALL get_qs_kind(qs_kind_set(ikind), alpha_core_charge=alpha_core_charge, &
ccore_charge=ccore_charge)
pab(1, 1) = -ccore_charge
IF (alpha_core_charge == 0.0_dp .OR. pab(1, 1) == 0.0_dp) CYCLE

59
src/qs_interactions.F
View file

@ -41,6 +41,7 @@ MODULE qs_interactions
USE paw_proj_set_types, ONLY: get_paw_proj_set,&
paw_proj_set_type,&
set_paw_proj_set
USE qs_cneo_types, ONLY: cneo_potential_type
USE qs_kind_types, ONLY: get_qs_kind,&
qs_kind_type
USE string_utilities, ONLY: uppercase
@ -98,20 +99,22 @@ CONTAINS
zet
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: c_nonlocal
TYPE(all_potential_type), POINTER :: all_potential
TYPE(cneo_potential_type), POINTER :: cneo_potential
TYPE(gth_potential_type), POINTER :: gth_potential
TYPE(gto_basis_set_type), POINTER :: aux_basis_set, aux_fit_basis_set, aux_gw_basis, &
aux_opt_basis_set, gapw_1c_basis, harris_basis, lri_basis, mao_basis, min_basis_set, &
orb_basis_set, p_lri_basis, rhoin_basis, ri_aux_basis_set, ri_basis, ri_xas_basis, &
soft_basis, tda_k_basis
nuc_basis_set, orb_basis_set, p_lri_basis, rhoin_basis, ri_aux_basis_set, ri_basis, &
ri_xas_basis, soft_basis, tda_k_basis
TYPE(paw_proj_set_type), POINTER :: paw_proj_set
TYPE(sgp_potential_type), POINTER :: sgp_potential
CALL timeset(routineN, handle)
NULLIFY (all_potential, gth_potential, sgp_potential)
NULLIFY (all_potential, gth_potential, sgp_potential, cneo_potential)
NULLIFY (aux_basis_set, aux_fit_basis_set, aux_gw_basis, tda_k_basis, &
harris_basis, lri_basis, mao_basis, orb_basis_set, p_lri_basis, ri_aux_basis_set, &
ri_basis, ri_xas_basis, soft_basis, gapw_1c_basis, aux_opt_basis_set, min_basis_set)
ri_basis, ri_xas_basis, soft_basis, gapw_1c_basis, aux_opt_basis_set, min_basis_set, &
nuc_basis_set)
NULLIFY (nprj_ppnl, nprj)
NULLIFY (alpha_ppnl, cexp_ppl, cprj_ppnl, zet)
@ -137,12 +140,14 @@ CONTAINS
CALL get_qs_kind(qs_kind_set(ikind), basis_set=gapw_1c_basis, basis_type="GAPW_1C")
CALL get_qs_kind(qs_kind_set(ikind), basis_set=tda_k_basis, basis_type="TDA_HFX")
CALL get_qs_kind(qs_kind_set(ikind), basis_set=rhoin_basis, basis_type="RHOIN")
CALL get_qs_kind(qs_kind_set(ikind), basis_set=nuc_basis_set, basis_type="NUC")
CALL get_qs_kind(qs_kind_set(ikind), &
paw_proj_set=paw_proj_set, &
paw_atom=paw_atom, &
all_potential=all_potential, &
gth_potential=gth_potential, &
sgp_potential=sgp_potential)
sgp_potential=sgp_potential, &
cneo_potential=cneo_potential)
! Calculate the orbital basis function radii ***
! For ALL electron this has to come before the calculation of the
@ -310,6 +315,14 @@ CONTAINS
core_charge_radius=core_charge_radius, &
ppl_radius=ppl_radius, &
ppnl_radius=ppnl_radius)
ELSE IF (ASSOCIATED(cneo_potential)) THEN
IF (ASSOCIATED(nuc_basis_set)) THEN
CALL init_interaction_radii_orb_basis(nuc_basis_set, qs_control%eps_pgf_orb/ &
SQRT(cneo_potential%zeff))
END IF
END IF
! Calculate the aux fit orbital basis function radii
@ -671,11 +684,11 @@ CONTAINS
REAL(KIND=dp), DIMENSION(:, :), POINTER :: pgf_radius
TYPE(cp_logger_type), POINTER :: logger
TYPE(gto_basis_set_type), POINTER :: aux_basis_set, lri_basis_set, &
orb_basis_set
nuc_basis_set, orb_basis_set
NULLIFY (logger)
logger => cp_get_default_logger()
NULLIFY (aux_basis_set, orb_basis_set, lri_basis_set)
NULLIFY (aux_basis_set, orb_basis_set, lri_basis_set, nuc_basis_set)
bas = " "
bas(1:3) = basis(1:3)
CALL uppercase(bas)
@ -685,6 +698,8 @@ CONTAINS
bas = "AUXILLIARY"
ELSE IF (bas(1:3) == "LRI") THEN
bas = "LOCAL RI"
ELSE IF (bas(1:3) == "NUC") THEN
bas = "NUCLEAR "
ELSE
CPABORT("Undefined basis set type")
END IF
@ -721,10 +736,18 @@ CONTAINS
CALL get_gto_basis_set(gto_basis_set=lri_basis_set, &
kind_radius=kind_radius)
END IF
ELSE IF (bas(1:3) == "NUC") THEN
CALL get_qs_kind(qs_kind_set(ikind), basis_set=nuc_basis_set, basis_type="NUC")
IF (ASSOCIATED(nuc_basis_set)) THEN
CALL get_gto_basis_set(gto_basis_set=nuc_basis_set, &
kind_radius=kind_radius, &
short_kind_radius=short_kind_radius)
END IF
ELSE
CPABORT("Undefined basis set type")
END IF
IF (ASSOCIATED(aux_basis_set) .OR. ASSOCIATED(orb_basis_set)) THEN
IF (ASSOCIATED(aux_basis_set) .OR. ASSOCIATED(orb_basis_set) .OR. &
ASSOCIATED(nuc_basis_set)) THEN
WRITE (UNIT=output_unit, FMT="(T45,I5,T53,A5,T57,F12.6,T69,F12.6)") &
ikind, name, kind_radius*conv, short_kind_radius*conv
ELSE
@ -766,10 +789,18 @@ CONTAINS
nset=nset, &
set_radius=set_radius)
END IF
ELSE IF (bas(1:3) == "NUC") THEN
CALL get_qs_kind(qs_kind_set(ikind), basis_set=nuc_basis_set, basis_type="NUC")
IF (ASSOCIATED(nuc_basis_set)) THEN
CALL get_gto_basis_set(gto_basis_set=nuc_basis_set, &
nset=nset, &
set_radius=set_radius)
END IF
ELSE
CPABORT("Undefined basis set type")
END IF
IF (ASSOCIATED(aux_basis_set) .OR. ASSOCIATED(orb_basis_set)) THEN
IF (ASSOCIATED(aux_basis_set) .OR. ASSOCIATED(orb_basis_set) .OR. &
ASSOCIATED(nuc_basis_set)) THEN
WRITE (UNIT=output_unit, FMT="(T50,I5,T57,A5,(T63,I5,T69,F12.6))") &
ikind, name, (iset, set_radius(iset)*conv, iset=1, nset)
ELSE
@ -813,11 +844,19 @@ CONTAINS
npgf=npgf, &
pgf_radius=pgf_radius)
END IF
ELSE IF (bas(1:3) == "NUC") THEN
CALL get_qs_kind(qs_kind_set(ikind), basis_set=nuc_basis_set, basis_type="NUC")
IF (ASSOCIATED(nuc_basis_set)) THEN
CALL get_gto_basis_set(gto_basis_set=nuc_basis_set, &
nset=nset, &
npgf=npgf, &
pgf_radius=pgf_radius)
END IF
ELSE
CPABORT("Undefined basis type")
END IF
IF (ASSOCIATED(aux_basis_set) .OR. ASSOCIATED(orb_basis_set) .OR. &
ASSOCIATED(lri_basis_set)) THEN
ASSOCIATED(lri_basis_set) .OR. ASSOCIATED(nuc_basis_set)) THEN
DO iset = 1, nset
WRITE (UNIT=output_unit, FMT="(T50,I5,T57,A5,T63,I5,(T69,F12.6))") &
ikind, name, iset, &

218
src/qs_kind_types.F
View file

@ -84,6 +84,12 @@ MODULE qs_kind_types
USE physcon, ONLY: angstrom,&
bohr,&
evolt
USE qs_cneo_types, ONLY: allocate_cneo_potential,&
cneo_potential_type,&
deallocate_cneo_potential,&
get_cneo_potential,&
set_cneo_potential,&
write_cneo_potential
USE qs_dftb_types, ONLY: qs_dftb_atom_type
USE qs_dftb_utils, ONLY: deallocate_dftb_atom_param,&
get_dftb_atom_param,&
@ -181,6 +187,7 @@ MODULE qs_kind_types
TYPE(xtb_atom_type), POINTER :: xtb_parameter => Null()
!
TYPE(atom_upfpot_type), POINTER :: upf_potential => Null()
TYPE(cneo_potential_type), POINTER :: cneo_potential => Null()
!
TYPE(basis_set_container_type), &
DIMENSION(20) :: basis_sets = basis_set_container_type()
@ -246,7 +253,8 @@ MODULE qs_kind_types
set_qs_kind, &
write_qs_kind_set, &
write_gto_basis_sets, &
init_atom_electronic_state, set_pseudo_state
init_atom_electronic_state, set_pseudo_state, &
init_cneo_basis_set
! Public data types
PUBLIC :: qs_kind_type, pao_potential_type, pao_descriptor_type
@ -287,6 +295,9 @@ CONTAINS
CALL atom_release_upf(qs_kind_set(ikind)%upf_potential)
DEALLOCATE (qs_kind_set(ikind)%upf_potential)
END IF
IF (ASSOCIATED(qs_kind_set(ikind)%cneo_potential)) THEN
CALL deallocate_cneo_potential(qs_kind_set(ikind)%cneo_potential)
END IF
IF (ASSOCIATED(qs_kind_set(ikind)%se_parameter)) THEN
CALL semi_empirical_release(qs_kind_set(ikind)%se_parameter)
END IF
@ -370,6 +381,7 @@ CONTAINS
!> \param gth_potential ...
!> \param sgp_potential ...
!> \param upf_potential ...
!> \param cneo_potential ...
!> \param se_parameter ...
!> \param dftb_parameter ...
!> \param xtb_parameter ...
@ -437,7 +449,7 @@ CONTAINS
SUBROUTINE get_qs_kind(qs_kind, &
basis_set, basis_type, ncgf, nsgf, &
all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, &
se_parameter, dftb_parameter, xtb_parameter, &
cneo_potential, se_parameter, dftb_parameter, xtb_parameter, &
dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, &
alpha_core_charge, ccore_charge, core_charge, core_charge_radius, &
paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, &
@ -461,6 +473,7 @@ CONTAINS
TYPE(gth_potential_type), OPTIONAL, POINTER :: gth_potential
TYPE(sgp_potential_type), OPTIONAL, POINTER :: sgp_potential
TYPE(atom_upfpot_type), OPTIONAL, POINTER :: upf_potential
TYPE(cneo_potential_type), OPTIONAL, POINTER :: cneo_potential
TYPE(semi_empirical_type), OPTIONAL, POINTER :: se_parameter
TYPE(qs_dftb_atom_type), OPTIONAL, POINTER :: dftb_parameter
TYPE(xtb_atom_type), OPTIONAL, POINTER :: xtb_parameter
@ -558,6 +571,7 @@ CONTAINS
IF (PRESENT(gth_potential)) gth_potential => qs_kind%gth_potential
IF (PRESENT(sgp_potential)) sgp_potential => qs_kind%sgp_potential
IF (PRESENT(upf_potential)) upf_potential => qs_kind%upf_potential
IF (PRESENT(cneo_potential)) cneo_potential => qs_kind%cneo_potential
IF (PRESENT(se_parameter)) se_parameter => qs_kind%se_parameter
IF (PRESENT(dftb_parameter)) dftb_parameter => qs_kind%dftb_parameter
IF (PRESENT(xtb_parameter)) xtb_parameter => qs_kind%xtb_parameter
@ -575,6 +589,8 @@ CONTAINS
ELSE IF (ASSOCIATED(qs_kind%sgp_potential)) THEN
CALL get_potential(potential=qs_kind%sgp_potential, &
alpha_core_charge=alpha_core_charge)
ELSE IF (ASSOCIATED(qs_kind%cneo_potential)) THEN
CPABORT("CNEO ALPHA CORE CHARGE NOT AVAILABLE")
ELSE
alpha_core_charge = 1.0_dp
END IF
@ -590,7 +606,9 @@ CONTAINS
CALL get_potential(potential=qs_kind%sgp_potential, &
ccore_charge=ccore_charge)
ELSE IF (ASSOCIATED(qs_kind%upf_potential)) THEN
CPABORT("UPF CCORE CHARGE RADIUS NOT AVAILABLE")
CPABORT("UPF CCORE CHARGE NOT AVAILABLE")
ELSE IF (ASSOCIATED(qs_kind%cneo_potential)) THEN
CPABORT("CNEO CCORE CHARGE NOT AVAILABLE")
ELSE
ccore_charge = 0.0_dp
END IF
@ -607,6 +625,8 @@ CONTAINS
core_charge_radius=core_charge_radius)
ELSE IF (ASSOCIATED(qs_kind%upf_potential)) THEN
CPABORT("UPF CORE CHARGE RADIUS NOT AVAILABLE")
ELSE IF (ASSOCIATED(qs_kind%cneo_potential)) THEN
CPABORT("CNEO CORE CHARGE RADIUS NOT AVAILABLE")
ELSE
core_charge_radius = 0.0_dp
END IF
@ -623,6 +643,9 @@ CONTAINS
zeff=core_charge)
ELSE IF (ASSOCIATED(qs_kind%upf_potential)) THEN
CPABORT("UPF CORE CHARGE NOT AVAILABLE")
ELSE IF (ASSOCIATED(qs_kind%cneo_potential)) THEN
CALL get_cneo_potential(potential=qs_kind%cneo_potential, &
zeff=core_charge)
ELSE
core_charge = 0.0_dp
END IF
@ -643,6 +666,8 @@ CONTAINS
CALL get_potential(potential=qs_kind%sgp_potential, zeff=zeff)
ELSE IF (ASSOCIATED(qs_kind%upf_potential)) THEN
zeff = qs_kind%upf_potential%zion
ELSE IF (ASSOCIATED(qs_kind%cneo_potential)) THEN
CALL get_cneo_potential(potential=qs_kind%cneo_potential, zeff=zeff)
ELSE
zeff = 0.0_dp
END IF
@ -867,6 +892,9 @@ CONTAINS
!> \param basis_type ...
!> \param total_zeff_corr ... [SGh]
!> \param npgf_seg total number of primitive GTOs in "segmented contraction format"
!> \param cneo_potential_present ...
!> \param nkind_q ...
!> \param natom_q ...
! **************************************************************************************************
SUBROUTINE get_qs_kind_set(qs_kind_set, &
all_potential_present, tnadd_potential_present, gth_potential_present, &
@ -875,7 +903,8 @@ CONTAINS
ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, &
nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, &
basis_rcut, &
basis_type, total_zeff_corr, npgf_seg)
basis_type, total_zeff_corr, npgf_seg, &
cneo_potential_present, nkind_q, natom_q)
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
LOGICAL, INTENT(OUT), OPTIONAL :: all_potential_present, tnadd_potential_present, &
@ -890,6 +919,8 @@ CONTAINS
CHARACTER(len=*), OPTIONAL :: basis_type
REAL(KIND=dp), INTENT(OUT), OPTIONAL :: total_zeff_corr
INTEGER, INTENT(OUT), OPTIONAL :: npgf_seg
LOGICAL, INTENT(OUT), OPTIONAL :: cneo_potential_present
INTEGER, INTENT(OUT), OPTIONAL :: nkind_q, natom_q
CHARACTER(len=default_string_length) :: my_basis_type
INTEGER :: ikind, imax, lmax_rho0_kind, &
@ -899,6 +930,7 @@ CONTAINS
LOGICAL :: dft_plus_u_atom, ecp_semi_local, paw_atom
REAL(KIND=dp) :: brcut, zeff, zeff_correction
TYPE(all_potential_type), POINTER :: all_potential
TYPE(cneo_potential_type), POINTER :: cneo_potential
TYPE(gth_potential_type), POINTER :: gth_potential
TYPE(gto_basis_set_type), POINTER :: tmp_basis_set
TYPE(local_potential_type), POINTER :: tnadd_potential
@ -939,6 +971,9 @@ CONTAINS
IF (PRESENT(tnadd_potential_present)) tnadd_potential_present = .FALSE.
IF (PRESENT(gth_potential_present)) gth_potential_present = .FALSE.
IF (PRESENT(sgp_potential_present)) sgp_potential_present = .FALSE.
IF (PRESENT(cneo_potential_present)) cneo_potential_present = .FALSE.
IF (PRESENT(nkind_q)) nkind_q = 0
IF (PRESENT(natom_q)) natom_q = 0
IF (PRESENT(paw_atom_present)) paw_atom_present = .FALSE.
IF (PRESENT(max_ngrid_rad)) max_ngrid_rad = 0
IF (PRESENT(max_sph_harm)) max_sph_harm = 0
@ -955,6 +990,7 @@ CONTAINS
tnadd_potential=tnadd_potential, &
gth_potential=gth_potential, &
sgp_potential=sgp_potential, &
cneo_potential=cneo_potential, &
paw_proj_set=paw_proj_set, &
dftb_parameter=dftb_parameter, &
ngrid_rad=ngrid_rad, &
@ -1010,8 +1046,12 @@ CONTAINS
maxppnl = MAX(maxppnl, imax)
END IF
CALL get_basis_from_container(qs_kind%basis_sets, basis_set=tmp_basis_set, &
basis_type=my_basis_type)
IF (my_basis_type(1:3) == "NUC" .AND. .NOT. ASSOCIATED(cneo_potential)) THEN
NULLIFY (tmp_basis_set)
ELSE
CALL get_basis_from_container(qs_kind%basis_sets, basis_set=tmp_basis_set, &
basis_type=my_basis_type)
END IF
IF (PRESENT(maxcgf)) THEN
IF (ASSOCIATED(tmp_basis_set)) THEN
@ -1139,6 +1179,10 @@ CONTAINS
ELSE IF (ASSOCIATED(qs_kind%sgp_potential)) THEN
CALL get_potential(potential=qs_kind%sgp_potential, &
zeff=zeff, zeff_correction=zeff_correction)
ELSE IF (ASSOCIATED(qs_kind%cneo_potential)) THEN
CALL get_cneo_potential(potential=qs_kind%cneo_potential, &
zeff=zeff)
zeff_correction = 0.0_dp
ELSE
zeff = 0.0_dp
zeff_correction = 0.0_dp
@ -1197,6 +1241,24 @@ CONTAINS
END IF
END IF
IF (PRESENT(cneo_potential_present)) THEN
IF (ASSOCIATED(cneo_potential)) THEN
cneo_potential_present = .TRUE.
END IF
END IF
IF (PRESENT(nkind_q)) THEN
IF (ASSOCIATED(cneo_potential)) THEN
nkind_q = nkind_q + 1
END IF
END IF
IF (PRESENT(natom_q)) THEN
IF (ASSOCIATED(cneo_potential)) THEN
natom_q = natom_q + qs_kind_set(ikind)%natom
END IF
END IF
IF (PRESENT(paw_atom_present)) THEN
IF (paw_atom) THEN
paw_atom_present = .TRUE.
@ -1354,9 +1416,13 @@ CONTAINS
NULLIFY (soft_basis)
CALL allocate_gto_basis_set(soft_basis)
! Quantum nuclear wave functions are very localized. Even if soft
! electronic basis is used, the atomic kind needs PAW treatment
! because the local Hartree potential is needed.
CALL create_soft_basis(orb_basis, soft_basis, &
qs_control%gapw_control%eps_fit, rc, paw_atom, &
qs_control%gapw_control%force_paw, gpw)
(qs_control%gapw_control%force_paw .OR. &
ASSOCIATED(qs_kind%cneo_potential)), gpw)
CALL add_basis_set_to_container(qs_kind%basis_sets, soft_basis, "ORB_SOFT")
CALL set_qs_kind(qs_kind=qs_kind, paw_atom=paw_atom)
@ -1494,6 +1560,51 @@ CONTAINS
END SUBROUTINE init_gapw_nlcc
! **************************************************************************************************
!> \brief ...
!> \param qs_kind_set ...
!> \param qs_control ...
! **************************************************************************************************
SUBROUTINE init_cneo_basis_set(qs_kind_set, qs_control)
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(qs_control_type), POINTER :: qs_control
INTEGER :: ikind, nkind
LOGICAL :: paw_atom
REAL(dp) :: rc
TYPE(gto_basis_set_type), POINTER :: orb_basis, soft_basis
TYPE(qs_kind_type), POINTER :: qs_kind
IF (ASSOCIATED(qs_kind_set)) THEN
nkind = SIZE(qs_kind_set)
DO ikind = 1, nkind
qs_kind => qs_kind_set(ikind)
IF (ASSOCIATED(qs_kind%cneo_potential)) THEN
CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis, basis_type="NUC", &
hard_radius=rc)
NULLIFY (soft_basis)
CALL allocate_gto_basis_set(soft_basis)
CALL create_soft_basis(orb_basis, soft_basis, &
qs_control%gapw_control%eps_fit/ &
SQRT(qs_kind%cneo_potential%zeff), &
rc, paw_atom, .TRUE., .FALSE.)
CALL add_basis_set_to_container(qs_kind%basis_sets, soft_basis, "NUC_SOFT")
END IF
END DO
ELSE
CPABORT("The pointer qs_kind_set is not associated")
END IF
END SUBROUTINE init_cneo_basis_set
! **************************************************************************************************
!> \brief Read an atomic kind data set from the input file.
!> \param qs_kind ...
@ -1538,9 +1649,10 @@ CONTAINS
nb_rep, nexp, ngauss, nlcc, nloc, nnl, norbitals, npaodesc, npaopot, nppnl, nspin, nu, z
INTEGER, DIMENSION(:), POINTER :: add_el, elec_conf, orbitals
LOGICAL :: check, ecp_semi_local, explicit, explicit_basis, explicit_J, explicit_kgpot, &
explicit_potential, explicit_U, explicit_u_m_j, nobasis, section_enabled, &
explicit_potential, explicit_U, explicit_u_m_j, nobasis, nobasis_nuc, section_enabled, &
subsection_enabled, update_input
REAL(KIND=dp) :: alpha, ccore, r, rc, zeff_correction
REAL(KIND=dp) :: alpha, ccore, mass, r, rc, &
zeff_correction
REAL(KIND=dp), DIMENSION(6) :: error
REAL(KIND=dp), DIMENSION(:), POINTER :: a_nl, aloc, anlcc, cloc, cnlcc, nelec
REAL(KIND=dp), DIMENSION(:, :), POINTER :: h_nl
@ -2212,7 +2324,7 @@ CONTAINS
! Allocate and initialise the potential data set structure
IF (potential_name /= '') THEN
SELECT CASE (TRIM(potential_type))
CASE ("ALL", "ECP")
CASE ("ALL", "ECP", "CNEO")
CALL cp_abort(__LOCATION__, &
"PW DFT calculations only with potential type UPF or GTH possible."// &
" <"//TRIM(potential_type)//"> was specified "// &
@ -2518,6 +2630,43 @@ CONTAINS
CALL set_potential(qs_kind%sgp_potential, elec_conf=elec_conf)
CALL atom_release_upf(upfpot)
CALL atom_sgp_release(sgppot)
CASE ("CNEO")
IF (zeff_correction /= 0.0_dp) &
CPABORT("CORE_CORRECTION is not compatible with CNEO")
CALL allocate_cneo_potential(qs_kind%cneo_potential)
CALL set_cneo_potential(qs_kind%cneo_potential, z=z)
mass = 0.0_dp
! Input mass is the mass of the neutral atom, not the nucleus.
! The mass of electrons will get subtracted later.
! In principle, the most abundant pure isotope mass should be used.
IF (k_rep > 0) THEN
CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
keyword_name="MASS", n_rep_val=i)
IF (i > 0) CALL section_vals_val_get(kind_section, i_rep_section=k_rep, &
keyword_name="MASS", r_val=mass)
END IF
! Remove electron mass from atomic mass to get nuclear mass
IF (mass - REAL(z, dp)*0.000548579909_dp > 0.0_dp) THEN
mass = mass - REAL(z, dp)*0.000548579909_dp
CALL set_cneo_potential(qs_kind%cneo_potential, mass=mass)
END IF
! In case the mass is not set by user, get the default mass from z
CALL get_cneo_potential(qs_kind%cneo_potential, mass=mass)
! Warn if the mass is from ptable
IF (ABS(mass + REAL(z, dp)*0.000548579909_dp - ptable(z)%amass) < 1.e-4_dp) THEN
CALL cp_warn(__LOCATION__, &
"Atomic mass of the atomic kind <"//TRIM(qs_kind%name)// &
"> is very close to its average mass. Is it a pure isotope? "// &
"Pure isotopes are preferable for CNEO. "// &
"(e.g., mass of 1H is 1.007825, not 1.00794)")
END IF
CALL get_qs_kind(qs_kind, elec_conf=elec_conf)
IF (.NOT. ASSOCIATED(elec_conf)) THEN
CALL get_cneo_potential(potential=qs_kind%cneo_potential, elec_conf=elec_conf)
CALL set_qs_kind(qs_kind, elec_conf=elec_conf)
ELSE
CALL set_cneo_potential(potential=qs_kind%cneo_potential, elec_conf=elec_conf)
END IF
CASE DEFAULT
CALL cp_abort(__LOCATION__, &
"An invalid potential type <"// &
@ -2566,6 +2715,28 @@ CONTAINS
END SELECT
END IF
END IF
! check that we have a nuclear orbital basis set if CNEO is requested
nobasis_nuc = ASSOCIATED(qs_kind%cneo_potential)
DO i = 1, SIZE(qs_kind%basis_sets)
NULLIFY (tmp_basis_set)
CALL get_basis_from_container(qs_kind%basis_sets, basis_set=tmp_basis_set, &
inumbas=i, basis_type=basis_type)
IF (basis_type == "NUC") THEN
nobasis_nuc = .FALSE.
IF (.NOT. ASSOCIATED(qs_kind%cneo_potential)) THEN
CALL cp_warn(__LOCATION__, &
"POTENTIAL is not set to CNEO, NUC type basis set for KIND <"// &
TRIM(qs_kind%name)//"> will be ignored!")
END IF
END IF
END DO
IF (nobasis_nuc) THEN
CALL cp_abort(__LOCATION__, &
"No NUC type basis set was defined for the "// &
"atomic kind <"//TRIM(qs_kind%name)//">, which is required by "// &
"POTENTIAL CNEO.")
END IF
END SELECT
CALL timestop(handle)
@ -2876,6 +3047,8 @@ CONTAINS
CALL set_potential(potential=qs_kind%gth_potential, zeff=zeff)
ELSE IF (ASSOCIATED(qs_kind%sgp_potential)) THEN
CALL set_potential(potential=qs_kind%sgp_potential, zeff=zeff)
ELSE IF (ASSOCIATED(qs_kind%cneo_potential)) THEN
CPABORT("CNEO potential ZEFF should not be manually set")
END IF
END IF
IF (PRESENT(ghost)) qs_kind%ghost = ghost
@ -2954,6 +3127,12 @@ CONTAINS
bstring = "RI XAS Basis Set"
CASE ("RI_HFX")
bstring = "RI HFX Basis Set"
CASE ("NUC")
bstring = "Nuclear Basis Set"
do_print = .FALSE.
CASE ("NUC_SOFT")
bstring = "Nuclear Soft Basis Set"
do_print = .FALSE.
END SELECT
IF (do_print) THEN
@ -3046,6 +3225,19 @@ CONTAINS
END IF
END IF
END IF
IF (ASSOCIATED(qs_kind%cneo_potential)) THEN
WRITE (UNIT=output_unit, FMT="(/,T6,A)") &
"The nuclei of this atomic kind are quantum mechanical (CNEO)"
CALL write_cneo_potential(qs_kind%cneo_potential, output_unit)
NULLIFY (tmp_basis)
CALL get_basis_from_container(qs_kind%basis_sets, basis_set=tmp_basis, &
basis_type="NUC")
CALL write_orb_basis_set(tmp_basis, output_unit, "Nuclear Basis Set")
NULLIFY (tmp_basis)
CALL get_basis_from_container(qs_kind%basis_sets, basis_set=tmp_basis, &
basis_type="NUC_SOFT")
CALL write_orb_basis_set(tmp_basis, output_unit, "Nuclear Soft Basis Set")
END IF
ELSE
CPABORT("")
END IF
@ -3162,6 +3354,12 @@ CONTAINS
bstring = "LRI Basis Set for TDDFPT"
CASE ("RI_HFX")
bstring = "RI HFX Basis Set"
CASE ("NUC")
bstring = "Nuclear Basis Set"
IF (.NOT. ASSOCIATED(qs_kind%cneo_potential)) NULLIFY (tmp_basis)
CASE ("NUC_SOFT")
bstring = "Nuclear Soft Basis Set"
IF (.NOT. ASSOCIATED(qs_kind%cneo_potential)) NULLIFY (tmp_basis)
END SELECT
IF (ASSOCIATED(tmp_basis)) CALL write_gto_basis_set(tmp_basis, output_unit, bstring)

6
src/qs_kpp1_env_methods.F
View file

@ -209,8 +209,12 @@ CONTAINS
DO ispin = 2, nspins
CALL pw_axpy(rho1_g(ispin), rho1_tot_gspace)
END DO
IF (gapw) &
IF (gapw) THEN
CALL pw_axpy(p_env%local_rho_set%rho0_mpole%rho0_s_gs, rho1_tot_gspace)
IF (ASSOCIATED(p_env%local_rho_set%rho0_mpole%rhoz_cneo_s_gs)) THEN
CALL pw_axpy(p_env%local_rho_set%rho0_mpole%rhoz_cneo_s_gs, rho1_tot_gspace)
END IF
END IF
scf_section => section_vals_get_subs_vals(input, "DFT%SCF")
IF (cp_print_key_should_output(logger%iter_info, scf_section, "PRINT%TOTAL_DENSITIES") &

45
src/qs_ks_methods.F
View file

@ -215,9 +215,8 @@ CONTAINS
TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, v_rspace_embed, v_rspace_new, &
v_rspace_new_aux_fit, v_tau_rspace, &
v_tau_rspace_aux_fit
TYPE(pw_r3d_rs_type), POINTER :: rho0_s_rs, rho_nlcc, v_hartree_rspace, &
v_sccs_rspace, v_sic_rspace, &
v_spin_ddapc_rest_r, vee, vppl_rspace
TYPE(pw_r3d_rs_type), POINTER :: rho0_s_rs, rho_nlcc, rhoz_cneo_s_rs, v_hartree_rspace, &
v_sccs_rspace, v_sic_rspace, v_spin_ddapc_rest_r, vee, vppl_rspace
TYPE(qs_energy_type), POINTER :: energy
TYPE(qs_ks_env_type), POINTER :: ks_env
TYPE(qs_rho_type), POINTER :: rho, rho_struct, rho_xc
@ -431,12 +430,6 @@ CONTAINS
CALL calc_v_sic_rspace(v_sic_rspace, energy, qs_env, dft_control, rho, poisson_env, &
just_energy, calculate_forces, auxbas_pw_pool)
IF (gapw) THEN
CALL get_qs_env(qs_env, ecoul_1c=ecoul_1c, local_rho_set=local_rho_set)
CALL Vh_1c_gg_integrals(qs_env, energy%hartree_1c, ecoul_1c, local_rho_set, para_env, tddft=.FALSE., &
core_2nd=.FALSE.)
END IF
! Check if CDFT constraint is needed
CALL qs_ks_cdft_constraint(qs_env, auxbas_pw_pool, calculate_forces, cdft_control)
@ -450,9 +443,16 @@ CONTAINS
IF (.NOT. just_energy) THEN
IF (gapw) THEN
CALL get_qs_env(qs_env=qs_env, &
rho0_s_rs=rho0_s_rs)
rho0_s_rs=rho0_s_rs, &
rhoz_cneo_s_rs=rhoz_cneo_s_rs)
CPASSERT(ASSOCIATED(rho0_s_rs))
IF (ASSOCIATED(rhoz_cneo_s_rs)) THEN
CALL pw_axpy(rhoz_cneo_s_rs, rho0_s_rs)
END IF
ee_ener = ee_ener + pw_integral_ab(rho0_s_rs, vee)
IF (ASSOCIATED(rhoz_cneo_s_rs)) THEN
CALL pw_axpy(rhoz_cneo_s_rs, rho0_s_rs, -1.0_dp)
END IF
END IF
END IF
! the sign accounts for the charge of the electrons
@ -798,6 +798,18 @@ CONTAINS
v_qmmm=vee, scale=-1.0_dp)
END IF
CALL integrate_vhg0_rspace(qs_env, v_hartree_rspace, para_env, calculate_forces)
! Place Vh_1c_gg_integrals after integrate_vhg0_rspace for CNEO calculations
! because vhg0 integral is needed to build the complete nuclear equation
CALL get_qs_env(qs_env, ecoul_1c=ecoul_1c, local_rho_set=local_rho_set)
CALL Vh_1c_gg_integrals(qs_env, energy%hartree_1c, ecoul_1c, local_rho_set, para_env, tddft=.FALSE., &
core_2nd=.FALSE.)
! CNEO quantum nuclear core energy (kinetic + Z*erfc(r)/r potential from classical nuclei)
energy%core_cneo = 0.0_dp
IF (ASSOCIATED(local_rho_set%rhoz_cneo_set)) THEN
DO iatom = 1, SIZE(local_rho_set%rhoz_cneo_set)
energy%core_cneo = energy%core_cneo + local_rho_set%rhoz_cneo_set(iatom)%e_core
END DO
END IF
END IF
IF (gapw .OR. gapw_xc) THEN
@ -907,8 +919,8 @@ CONTAINS
END IF
! Sum all energy terms to obtain the total energy
energy%total = energy%core_overlap + energy%core_self + energy%core + energy%hartree + &
energy%hartree_1c + energy%exc + energy%exc1 + energy%ex + &
energy%total = energy%core_overlap + energy%core_self + energy%core_cneo + energy%core + &
energy%hartree + energy%hartree_1c + energy%exc + energy%exc1 + energy%ex + &
energy%dispersion + energy%gcp + energy%qmmm_el + energy%mulliken + &
SUM(energy%ddapc_restraint) + energy%s2_restraint + &
energy%dft_plus_u + energy%kTS + &
@ -947,7 +959,7 @@ CONTAINS
LOGICAL :: my_skip
TYPE(dft_control_type), POINTER :: dft_control
TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
TYPE(pw_c1d_gs_type), POINTER :: rho0_s_gs, rho_core
TYPE(pw_c1d_gs_type), POINTER :: rho0_s_gs, rho_core, rhoz_cneo_s_gs
TYPE(qs_charges_type), POINTER :: qs_charges
my_skip = .FALSE.
@ -960,10 +972,13 @@ CONTAINS
NULLIFY (rho_core)
CALL get_qs_env(qs_env=qs_env, rho_core=rho_core)
IF (dft_control%qs_control%gapw) THEN
NULLIFY (rho0_s_gs)
CALL get_qs_env(qs_env=qs_env, rho0_s_gs=rho0_s_gs)
NULLIFY (rho0_s_gs, rhoz_cneo_s_gs)
CALL get_qs_env(qs_env=qs_env, rho0_s_gs=rho0_s_gs, rhoz_cneo_s_gs=rhoz_cneo_s_gs)
CPASSERT(ASSOCIATED(rho0_s_gs))
CALL pw_copy(rho0_s_gs, rho_tot_gspace)
IF (ASSOCIATED(rhoz_cneo_s_gs)) THEN
CALL pw_axpy(rhoz_cneo_s_gs, rho_tot_gspace)
END IF
IF (dft_control%qs_control%gapw_control%nopaw_as_gpw) THEN
CALL pw_axpy(rho_core, rho_tot_gspace)
END IF

11
src/qs_ks_qmmm_methods.F
View file

@ -152,7 +152,7 @@ CONTAINS
INTEGER :: handle, ispin, output_unit
TYPE(cp_logger_type), POINTER :: logger
TYPE(dft_control_type), POINTER :: dft_control
TYPE(pw_r3d_rs_type), POINTER :: rho0_s_rs
TYPE(pw_r3d_rs_type), POINTER :: rho0_s_rs, rhoz_cneo_s_rs
TYPE(section_vals_type), POINTER :: input
CALL timeset(routineN, handle)
@ -174,9 +174,16 @@ CONTAINS
END DO
IF (dft_control%qs_control%gapw) THEN
CALL get_qs_env(qs_env=qs_env, &
rho0_s_rs=rho0_s_rs)
rho0_s_rs=rho0_s_rs, &
rhoz_cneo_s_rs=rhoz_cneo_s_rs)
CPASSERT(ASSOCIATED(rho0_s_rs))
IF (ASSOCIATED(rhoz_cneo_s_rs)) THEN
CALL pw_axpy(rhoz_cneo_s_rs, rho0_s_rs)
END IF
qmmm_energy = qmmm_energy + pw_integral_ab(rho0_s_rs, v_qmmm)
IF (ASSOCIATED(rhoz_cneo_s_rs)) THEN
CALL pw_axpy(rhoz_cneo_s_rs, rho0_s_rs, -1.0_dp)
END IF
END IF
CALL cp_print_key_finished_output(output_unit, logger, input, &

20
src/qs_ks_types.F
View file

@ -118,6 +118,7 @@ MODULE qs_ks_types
!> \param sab_almo: neighbor lists to create ALMO delocalization template
!> \param sab_kp: neighbor lists to create kp image cell lists
!> \param sab_kp_nosym: neighbor lists to create kp image cell lists, non-symmetric
!> \param sab_cneo: neighbor lists for the calculation of the quantum nuclear core Hamiltonian matrix
!>
!> \param kpoints information on the kpoints used
!> \param subsys the particles, molecules,... of this environment
@ -189,7 +190,8 @@ MODULE qs_ks_types
sab_lrc => Null(), &
sab_almo => Null(), &
sab_kp => Null(), &
sab_kp_nosym => Null()
sab_kp_nosym => Null(), &
sab_cneo => Null()
TYPE(task_list_type), POINTER :: task_list => Null()
TYPE(task_list_type), POINTER :: task_list_soft => Null()
@ -277,6 +279,7 @@ CONTAINS
!> \param sab_almo ...
!> \param sab_kp ...
!> \param sab_kp_nosym ...
!> \param sab_cneo ...
!> \param task_list ...
!> \param task_list_soft ...
!> \param kpoints ...
@ -322,7 +325,7 @@ CONTAINS
neighbor_list_id, &
sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, &
sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, &
sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, &
sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, &
task_list, task_list_soft, &
kpoints, do_kpoints, &
atomic_kind_set, qs_kind_set, cell, cell_ref, use_ref_cell, &
@ -351,7 +354,8 @@ CONTAINS
INTEGER, OPTIONAL :: neighbor_list_id
TYPE(neighbor_list_set_p_type), DIMENSION(:), OPTIONAL, POINTER :: sab_orb, sab_all, sac_ae, &
sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, &
sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym
sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, &
sab_cneo
TYPE(task_list_type), OPTIONAL, POINTER :: task_list, task_list_soft
TYPE(kpoint_type), OPTIONAL, POINTER :: kpoints
LOGICAL, OPTIONAL :: do_kpoints
@ -446,6 +450,7 @@ CONTAINS
IF (PRESENT(sab_almo)) sab_almo => ks_env%sab_almo
IF (PRESENT(sab_kp)) sab_kp => ks_env%sab_kp
IF (PRESENT(sab_kp_nosym)) sab_kp_nosym => ks_env%sab_kp_nosym
IF (PRESENT(sab_cneo)) sab_cneo => ks_env%sab_cneo
IF (PRESENT(dft_control)) dft_control => ks_env%dft_control
IF (PRESENT(dbcsr_dist)) dbcsr_dist => ks_env%dbcsr_dist
IF (PRESENT(distribution_2d)) distribution_2d => ks_env%distribution_2d
@ -542,6 +547,7 @@ CONTAINS
!> \param sab_almo ...
!> \param sab_kp ...
!> \param sab_kp_nosym ...
!> \param sab_cneo ...
!> \param task_list ...
!> \param task_list_soft ...
!> \param subsys ...
@ -565,7 +571,7 @@ CONTAINS
kpoints, &
sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, &
sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, &
sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, &
sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, &
task_list, task_list_soft, &
subsys, dft_control, dbcsr_dist, distribution_2d, pw_env, &
para_env, blacs_env)
@ -590,7 +596,8 @@ CONTAINS
TYPE(kpoint_type), OPTIONAL, POINTER :: kpoints
TYPE(neighbor_list_set_p_type), DIMENSION(:), OPTIONAL, POINTER :: sab_orb, sab_all, sac_ae, &
sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, &
sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym
sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, &
sab_cneo
TYPE(task_list_type), OPTIONAL, POINTER :: task_list, task_list_soft
TYPE(qs_subsys_type), OPTIONAL, POINTER :: subsys
TYPE(dft_control_type), OPTIONAL, POINTER :: dft_control
@ -666,6 +673,7 @@ CONTAINS
IF (PRESENT(sab_almo)) ks_env%sab_almo => sab_almo
IF (PRESENT(sab_kp)) ks_env%sab_kp => sab_kp
IF (PRESENT(sab_kp_nosym)) ks_env%sab_kp_nosym => sab_kp_nosym
IF (PRESENT(sab_cneo)) ks_env%sab_cneo => sab_cneo
IF (PRESENT(task_list)) ks_env%task_list => task_list
IF (PRESENT(task_list_soft)) ks_env%task_list_soft => task_list_soft
@ -808,6 +816,7 @@ CONTAINS
CALL release_neighbor_list_sets(ks_env%sab_almo)
CALL release_neighbor_list_sets(ks_env%sab_kp)
CALL release_neighbor_list_sets(ks_env%sab_kp_nosym)
CALL release_neighbor_list_sets(ks_env%sab_cneo)
IF (ASSOCIATED(ks_env%dft_control)) THEN
CALL dft_control_release(ks_env%dft_control)
DEALLOCATE (ks_env%dft_control)
@ -902,6 +911,7 @@ CONTAINS
CALL release_neighbor_list_sets(ks_env%sab_almo)
CALL release_neighbor_list_sets(ks_env%sab_kp)
CALL release_neighbor_list_sets(ks_env%sab_kp_nosym)
CALL release_neighbor_list_sets(ks_env%sab_cneo)
CALL kpoint_release(ks_env%kpoints)
CALL pw_env_release(ks_env%pw_env, ks_env%para_env)
END SUBROUTINE qs_ks_part_release

32
src/qs_ks_utils.F
View file

@ -931,7 +931,9 @@ CONTAINS
REAL(n_electrons, dp), &
"Core density on regular grids:", &
qs_charges%total_rho_core_rspace, &
qs_charges%total_rho_core_rspace - REAL(n_electrons + dft_control%charge, dp)
qs_charges%total_rho_core_rspace + &
qs_charges%total_rho1_hard_nuc - &
REAL(n_electrons + dft_control%charge, dp)
END IF
IF (dft_control%qs_control%gapw) THEN
tot1_h = qs_charges%total_rho1_hard(1)
@ -949,12 +951,28 @@ CONTAINS
tot_rho_r + tot1_h - tot1_s, &
"Total charge density (r-space): ", &
tot_rho_r + tot1_h - tot1_s &
+ qs_charges%total_rho_core_rspace, &
"Total Rho_soft + Rho0_soft (g-space):", &
qs_charges%total_rho_gspace
+ qs_charges%total_rho_core_rspace &
+ qs_charges%total_rho1_hard_nuc
IF (qs_charges%total_rho1_hard_nuc /= 0.0_dp) THEN
WRITE (UNIT=output_unit, FMT="(T3,A,T41,F20.10)") &
"Total CNEO nuc. char. den. (Lebedev): ", &
qs_charges%total_rho1_hard_nuc, &
"Total CNEO soft char. den. (Lebedev): ", &
qs_charges%total_rho1_soft_nuc_lebedev, &
"Total CNEO soft char. den. (r-space): ", &
qs_charges%total_rho1_soft_nuc_rspace, &
"Total soft Rho_e+n+0 (g-space):", &
qs_charges%total_rho_gspace
ELSE
WRITE (UNIT=output_unit, FMT="(T3,A,T41,F20.10)") &
"Total Rho_soft + Rho0_soft (g-space):", &
qs_charges%total_rho_gspace
END IF
END IF
qs_charges%background = tot_rho_r + tot1_h - tot1_s + &
qs_charges%total_rho_core_rspace
qs_charges%total_rho_core_rspace + &
qs_charges%total_rho1_hard_nuc
! only add total_rho1_hard_nuc for gapw as cneo requires gapw
ELSE IF (dft_control%qs_control%gapw_xc) THEN
tot1_h = qs_charges%total_rho1_hard(1)
tot1_s = qs_charges%total_rho1_soft(1)
@ -1153,6 +1171,10 @@ CONTAINS
WRITE (UNIT=output_unit, FMT="(T3,A,T41,2F20.10)") &
"S2 restraint (order_p,energy) : ", s2_order_p, energy%s2_restraint
END IF
IF (energy%core_cneo /= 0.0_dp) THEN
WRITE (UNIT=output_unit, FMT="(T3,A,T61,F20.10)") &
"CNEO| quantum nuclear core energy: ", energy%core_cneo
END IF
END IF ! output_unit
CALL cp_print_key_finished_output(output_unit, logger, input, &

7
src/qs_linres_epr_nablavks.F
View file

@ -52,6 +52,7 @@ MODULE qs_linres_epr_nablavks
USE pw_pool_types, ONLY: pw_pool_type
USE pw_types, ONLY: pw_c1d_gs_type,&
pw_r3d_rs_type
USE qs_cneo_types, ONLY: cneo_potential_type
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
USE qs_gapw_densities, ONLY: prepare_gapw_den
@ -125,6 +126,7 @@ CONTAINS
TYPE(all_potential_type), POINTER :: all_potential
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(cell_type), POINTER :: cell
TYPE(cneo_potential_type), POINTER :: cneo_potential
TYPE(cp_logger_type), POINTER :: logger
TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao
TYPE(dft_control_type), POINTER :: dft_control
@ -472,6 +474,7 @@ CONTAINS
gth_potential=gth_potential, &
sgp_potential=sgp_potential, &
all_potential=all_potential, &
cneo_potential=cneo_potential, &
paw_atom=paw_atom)
IF (ASSOCIATED(gth_potential)) THEN
@ -775,6 +778,10 @@ CONTAINS
CPABORT("EPR with SGP potentials is not implemented")
ELSE IF (ASSOCIATED(cneo_potential)) THEN
CPABORT("EPR with CNEO potentials is not implemented")
ELSE IF (ASSOCIATED(all_potential)) THEN
CALL get_potential(potential=all_potential, &

6
src/qs_linres_kernel.F
View file

@ -288,8 +288,12 @@ CONTAINS
DO ispin = 2, nspins
CALL pw_axpy(rho1_g(ispin), rho1_tot_gspace)
END DO
IF (gapw) &
IF (gapw) THEN
CALL pw_axpy(p_env%local_rho_set%rho0_mpole%rho0_s_gs, rho1_tot_gspace)
IF (ASSOCIATED(p_env%local_rho_set%rho0_mpole%rhoz_cneo_s_gs)) THEN
CALL pw_axpy(p_env%local_rho_set%rho0_mpole%rhoz_cneo_s_gs, rho1_tot_gspace)
END IF
END IF
IF (.NOT. (nspins == 1 .AND. lr_triplet)) THEN
CALL pw_poisson_solve(poisson_env, rho1_tot_gspace, &

41
src/qs_local_rho_types.F
View file

@ -11,6 +11,8 @@ MODULE qs_local_rho_types
USE mathconstants, ONLY: fourpi,&
pi
USE memory_utilities, ONLY: reallocate
USE qs_cneo_types, ONLY: deallocate_rhoz_cneo_set,&
rhoz_cneo_type
USE qs_grid_atom, ONLY: grid_atom_type
USE qs_harmonics_atom, ONLY: harmonics_atom_type
USE qs_rho0_types, ONLY: deallocate_rho0_atom,&
@ -45,7 +47,9 @@ MODULE qs_local_rho_types
TYPE(rho0_mpole_type), POINTER :: rho0_mpole => NULL()
TYPE(rho0_atom_type), DIMENSION(:), POINTER :: rho0_atom_set => NULL()
TYPE(rhoz_type), DIMENSION(:), POINTER :: rhoz_set => NULL()
REAL(dp) :: rhoz_tot = -1.0_dp
TYPE(rhoz_cneo_type), DIMENSION(:), POINTER :: rhoz_cneo_set => NULL()
REAL(dp) :: rhoz_tot = -1.0_dp, &
rhoz_cneo_tot = -1.0_dp
END TYPE local_rho_type
! Public Types
@ -152,9 +156,12 @@ CONTAINS
nkind = SIZE(rhoz_set)
DO ikind = 1, nkind
DEALLOCATE (rhoz_set(ikind)%r_coef)
DEALLOCATE (rhoz_set(ikind)%dr_coef)
DEALLOCATE (rhoz_set(ikind)%vr_coef)
IF (ASSOCIATED(rhoz_set(ikind)%r_coef)) &
DEALLOCATE (rhoz_set(ikind)%r_coef)
IF (ASSOCIATED(rhoz_set(ikind)%dr_coef)) &
DEALLOCATE (rhoz_set(ikind)%dr_coef)
IF (ASSOCIATED(rhoz_set(ikind)%vr_coef)) &
DEALLOCATE (rhoz_set(ikind)%vr_coef)
END DO
DEALLOCATE (rhoz_set)
@ -168,8 +175,10 @@ CONTAINS
!> \param rho0_atom_set ...
!> \param rho0_mpole ...
!> \param rhoz_set ...
!> \param rhoz_cneo_set ...
! **************************************************************************************************
SUBROUTINE get_local_rho(local_rho_set, rho_atom_set, rho0_atom_set, rho0_mpole, rhoz_set)
SUBROUTINE get_local_rho(local_rho_set, rho_atom_set, rho0_atom_set, rho0_mpole, rhoz_set, &
rhoz_cneo_set)
TYPE(local_rho_type), POINTER :: local_rho_set
TYPE(rho_atom_type), DIMENSION(:), OPTIONAL, &
@ -178,11 +187,14 @@ CONTAINS
POINTER :: rho0_atom_set
TYPE(rho0_mpole_type), OPTIONAL, POINTER :: rho0_mpole
TYPE(rhoz_type), DIMENSION(:), OPTIONAL, POINTER :: rhoz_set
TYPE(rhoz_cneo_type), DIMENSION(:), OPTIONAL, &
POINTER :: rhoz_cneo_set
IF (PRESENT(rho_atom_set)) rho_atom_set => local_rho_set%rho_atom_set
IF (PRESENT(rho0_atom_set)) rho0_atom_set => local_rho_set%rho0_atom_set
IF (PRESENT(rho0_mpole)) rho0_mpole => local_rho_set%rho0_mpole
IF (PRESENT(rhoz_set)) rhoz_set => local_rho_set%rhoz_set
IF (PRESENT(rhoz_cneo_set)) rhoz_cneo_set => local_rho_set%rhoz_cneo_set
END SUBROUTINE get_local_rho
@ -200,8 +212,10 @@ CONTAINS
NULLIFY (local_rho_set%rho0_atom_set)
NULLIFY (local_rho_set%rho0_mpole)
NULLIFY (local_rho_set%rhoz_set)
NULLIFY (local_rho_set%rhoz_cneo_set)
local_rho_set%rhoz_tot = 0.0_dp
local_rho_set%rhoz_cneo_tot = 0.0_dp
END SUBROUTINE local_rho_set_create
@ -230,6 +244,10 @@ CONTAINS
CALL deallocate_rhoz(local_rho_set%rhoz_set)
END IF
IF (ASSOCIATED(local_rho_set%rhoz_cneo_set)) THEN
CALL deallocate_rhoz_cneo_set(local_rho_set%rhoz_cneo_set)
END IF
DEALLOCATE (local_rho_set)
END IF
@ -242,9 +260,10 @@ CONTAINS
!> \param rho0_atom_set ...
!> \param rho0_mpole ...
!> \param rhoz_set ...
!> \param rhoz_cneo_set ...
! **************************************************************************************************
SUBROUTINE set_local_rho(local_rho_set, rho_atom_set, rho0_atom_set, rho0_mpole, &
rhoz_set)
rhoz_set, rhoz_cneo_set)
TYPE(local_rho_type), POINTER :: local_rho_set
TYPE(rho_atom_type), DIMENSION(:), OPTIONAL, &
@ -253,6 +272,8 @@ CONTAINS
POINTER :: rho0_atom_set
TYPE(rho0_mpole_type), OPTIONAL, POINTER :: rho0_mpole
TYPE(rhoz_type), DIMENSION(:), OPTIONAL, POINTER :: rhoz_set
TYPE(rhoz_cneo_type), DIMENSION(:), OPTIONAL, &
POINTER :: rhoz_cneo_set
IF (PRESENT(rho_atom_set)) THEN
IF (ASSOCIATED(local_rho_set%rho_atom_set)) THEN
@ -282,7 +303,13 @@ CONTAINS
local_rho_set%rhoz_set => rhoz_set
END IF
IF (PRESENT(rhoz_cneo_set)) THEN
IF (ASSOCIATED(local_rho_set%rhoz_cneo_set)) THEN
CALL deallocate_rhoz_cneo_set(local_rho_set%rhoz_cneo_set)
END IF
local_rho_set%rhoz_cneo_set => rhoz_cneo_set
END IF
END SUBROUTINE set_local_rho
END MODULE qs_local_rho_types

73
src/qs_neighbor_lists.F
View file

@ -67,6 +67,7 @@ MODULE qs_neighbor_lists
paw_proj_set_type
USE periodic_table, ONLY: ptable
USE physcon, ONLY: bohr
USE qs_cneo_types, ONLY: cneo_potential_type
USE qs_dftb_types, ONLY: qs_dftb_atom_type
USE qs_dftb_utils, ONLY: get_dftb_atom_param
USE qs_dispersion_types, ONLY: qs_dispersion_type
@ -295,20 +296,22 @@ CONTAINS
CHARACTER(LEN=default_string_length) :: print_key_path
INTEGER :: handle, hfx_pot, ikind, ingp, iw, jkind, &
maxatom, ngp, nkind, zat
LOGICAL :: all_potential_present, almo, dftb, do_hfx, dokp, gth_potential_present, &
lri_optbas, lrigpw, mic, molecule_only, nddo, paw_atom, paw_atom_present, rigpw, &
sgp_potential_present, xtb
LOGICAL :: all_potential_present, almo, cneo_potential_present, dftb, do_hfx, dokp, &
gth_potential_present, lri_optbas, lrigpw, mic, molecule_only, nddo, paw_atom, &
paw_atom_present, rigpw, sgp_potential_present, xtb
LOGICAL, ALLOCATABLE, DIMENSION(:) :: all_present, aux_fit_present, aux_present, &
core_present, default_present, nonbond1_atom, nonbond2_atom, oce_present, orb_present, &
ppl_present, ppnl_present, ri_present, xb1_atom, xb2_atom
cneo_present, core_present, default_present, nonbond1_atom, nonbond2_atom, oce_present, &
orb_present, ppl_present, ppnl_present, ri_present, xb1_atom, xb2_atom
REAL(dp) :: almo_rcov, almo_rvdw, eps_schwarz, &
omega, pdist, rcut, roperator, subcells
REAL(dp), ALLOCATABLE, DIMENSION(:) :: all_pot_rad, aux_fit_radius, c_radius, calpha, &
core_radius, oce_radius, orb_radius, ppl_radius, ppnl_radius, ri_radius, zeff
core_radius, nuc_orb_radius, oce_radius, orb_radius, ppl_radius, ppnl_radius, ri_radius, &
zeff
REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: pair_radius, pair_radius_lb
TYPE(all_potential_type), POINTER :: all_potential
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(cell_type), POINTER :: cell
TYPE(cneo_potential_type), POINTER :: cneo_potential
TYPE(cp_logger_type), POINTER :: logger
TYPE(dft_control_type), POINTER :: dft_control
TYPE(distribution_1d_type), POINTER :: distribution_1d
@ -316,13 +319,14 @@ CONTAINS
TYPE(ewald_environment_type), POINTER :: ewald_env
TYPE(gth_potential_type), POINTER :: gth_potential
TYPE(gto_basis_set_type), POINTER :: aux_basis_set, aux_fit_basis_set, &
orb_basis_set, ri_basis_set
nuc_basis_set, orb_basis_set, &
ri_basis_set
TYPE(kpoint_type), POINTER :: kpoints
TYPE(local_atoms_type), ALLOCATABLE, DIMENSION(:) :: atom2d
TYPE(molecule_type), DIMENSION(:), POINTER :: molecule_set
TYPE(neighbor_list_set_p_type), DIMENSION(:), POINTER :: saa_list, sab_all, sab_almo, &
sab_cn, sab_core, sab_gcp, sab_kp, sab_kp_nosym, sab_lrc, sab_orb, sab_scp, sab_se, &
sab_tbe, sab_vdw, sab_xb, sab_xtb_nonbond, sab_xtb_pp, sab_xtbe, sac_ae, sac_lri, &
sab_cn, sab_cneo, sab_core, sab_gcp, sab_kp, sab_kp_nosym, sab_lrc, sab_orb, sab_scp, &
sab_se, sab_tbe, sab_vdw, sab_xb, sab_xtb_nonbond, sab_xtb_pp, sab_xtbe, sac_ae, sac_lri, &
sac_ppl, sap_oce, sap_ppnl, soa_list, soo_list
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(paw_proj_set_type), POINTER :: paw_proj
@ -365,6 +369,7 @@ CONTAINS
NULLIFY (sab_almo)
NULLIFY (sab_kp)
NULLIFY (sab_kp_nosym)
NULLIFY (sab_cneo)
CALL get_qs_env(qs_env, &
ks_env=ks_env, &
@ -405,7 +410,8 @@ CONTAINS
sab_all=sab_all, &
sab_almo=sab_almo, &
sab_kp=sab_kp, &
sab_kp_nosym=sab_kp_nosym)
sab_kp_nosym=sab_kp_nosym, &
sab_cneo=sab_cneo)
dokp = (kpoints%nkp > 0)
nddo = dft_control%qs_control%semi_empirical
@ -437,7 +443,8 @@ CONTAINS
CALL get_qs_kind_set(qs_kind_set, paw_atom_present=paw_atom_present, &
gth_potential_present=gth_potential_present, &
sgp_potential_present=sgp_potential_present, &
all_potential_present=all_potential_present)
all_potential_present=all_potential_present, &
cneo_potential_present=cneo_potential_present)
CALL section_vals_val_get(qs_env%input, "DFT%SUBCELLS", r_val=subcells)
@ -469,6 +476,10 @@ CONTAINS
ALLOCATE (all_present(nkind), all_pot_rad(nkind))
all_pot_rad = 0.0_dp
END IF
IF (cneo_potential_present) THEN
ALLOCATE (cneo_present(nkind), nuc_orb_radius(nkind))
nuc_orb_radius = 0.0_dp
END IF
! Initialize the local data structures
ALLOCATE (atom2d(nkind))
@ -482,13 +493,15 @@ CONTAINS
CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis_set, basis_type="ORB")
CALL get_qs_kind(qs_kind_set(ikind), basis_set=aux_basis_set, basis_type="AUX")
CALL get_qs_kind(qs_kind_set(ikind), basis_set=aux_fit_basis_set, basis_type="AUX_FIT")
CALL get_qs_kind(qs_kind_set(ikind), basis_set=nuc_basis_set, basis_type="NUC")
CALL get_qs_kind(qs_kind_set(ikind), &
paw_proj_set=paw_proj, &
paw_atom=paw_atom, &
all_potential=all_potential, &
gth_potential=gth_potential, &
sgp_potential=sgp_potential)
sgp_potential=sgp_potential, &
cneo_potential=cneo_potential)
IF (dftb) THEN
! Set the interaction radius for the neighbor lists (DFTB case)
@ -519,15 +532,22 @@ CONTAINS
aux_fit_present(ikind) = .FALSE.
END IF
! core overlap
CALL get_qs_kind(qs_kind_set(ikind), &
alpha_core_charge=calpha(ikind), &
core_charge_radius=core_radius(ikind), &
zeff=zeff(ikind))
IF (zeff(ikind) /= 0._dp .AND. calpha(ikind) /= 0._dp) THEN
core_present(ikind) = .TRUE.
core_present(ikind) = .FALSE.
IF (ASSOCIATED(cneo_potential) .AND. ASSOCIATED(nuc_basis_set)) THEN
cneo_present(ikind) = .TRUE.
CALL get_gto_basis_set(gto_basis_set=nuc_basis_set, kind_radius=nuc_orb_radius(ikind))
ELSE
core_present(ikind) = .FALSE.
IF (cneo_potential_present) cneo_present(ikind) = .FALSE.
! core overlap
CALL get_qs_kind(qs_kind_set(ikind), &
alpha_core_charge=calpha(ikind), &
core_charge_radius=core_radius(ikind), &
zeff=zeff(ikind))
IF (zeff(ikind) /= 0._dp .AND. calpha(ikind) /= 0._dp) THEN
core_present(ikind) = .TRUE.
ELSE
core_present(ikind) = .FALSE.
END IF
END IF
! Pseudopotentials
@ -717,6 +737,16 @@ CONTAINS
END IF
END IF
! Build quantum nuclear orbital-classical nuclear ERFC potential list for CNEO
IF (cneo_potential_present) THEN
CALL pair_radius_setup(cneo_present, core_present, nuc_orb_radius, core_radius, pair_radius)
CALL build_neighbor_lists(sab_cneo, particle_set, atom2d, cell, pair_radius, &
subcells=subcells, symmetric=.FALSE., operator_type="PP", nlname="sab_cneo")
CALL set_ks_env(ks_env=ks_env, sab_cneo=sab_cneo)
CALL write_neighbor_lists(sab_cneo, particle_set, cell, para_env, neighbor_list_section, &
"/SAB_CNEO", "sab_cneo", "NUCLEAR ORBITAL ERFC POTENTIAL")
END IF
IF (nddo) THEN
! Semi-empirical neighbor lists
default_present = .TRUE.
@ -1000,6 +1030,9 @@ CONTAINS
IF (all_potential_present .OR. sgp_potential_present) THEN
DEALLOCATE (all_present, all_pot_rad)
END IF
IF (cneo_potential_present) THEN
DEALLOCATE (cneo_present, nuc_orb_radius)
END IF
CALL timestop(handle)

116
src/qs_rho0_ggrid.F
View file

@ -38,6 +38,10 @@ MODULE qs_rho0_ggrid
pw_pool_type
USE pw_types, ONLY: pw_c1d_gs_type,&
pw_r3d_rs_type
USE qs_cneo_ggrid, ONLY: integrate_vhgg_rspace,&
rhoz_cneo_s_grid_create
USE qs_cneo_types, ONLY: cneo_potential_type,&
rhoz_cneo_type
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
USE qs_force_types, ONLY: qs_force_type
@ -314,6 +318,10 @@ CONTAINS
CALL timestop(handle)
IF (rho0_mpole%do_cneo) THEN
CALL rhoz_cneo_s_grid_create(pw_env, rho0_mpole)
END IF
END SUBROUTINE rho0_s_grid_create
! **************************************************************************************************
@ -347,24 +355,27 @@ CONTAINS
CHARACTER(LEN=*), PARAMETER :: routineN = 'integrate_vhg0_rspace'
INTEGER :: auxbas_grid, bo(2), handle, i, iat, iatom, ic, icg, ico, ig1, ig2, igrid, ii, &
ikind, ipgf1, ipgf2, is, iset1, iset2, iso, iso1, iso2, ispin, j, l0_ikind, llmax, lmax0, &
lshell, lx, ly, lz, m1, m2, max_iso_not0_local, max_s_harm, maxl, maxso, mepos, n1, n2, &
nat, nch_ik, nch_max, ncurr, nset, nsotot, nspins, num_pe
ikind, ipgf1, ipgf2, is, iset1, iset2, iso, iso1, iso2, ispin, j, l0_ikind, llmax, &
llmax_nuc, lmax0, lshell, lx, ly, lz, m1, m2, max_iso_not0_local, max_s_harm, &
max_s_harm_nuc, maxl, maxl_nuc, maxso, maxso_nuc, mepos, n1, n2, nat, nch_ik, nch_max, &
ncurr, nset, nset_nuc, nsotot, nsotot_nuc, nspins, num_pe
INTEGER, ALLOCATABLE, DIMENSION(:) :: cg_n_list
INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cg_list
INTEGER, DIMENSION(:), POINTER :: atom_list, lmax, lmin, npgf
INTEGER, DIMENSION(:), POINTER :: atom_list, lmax, lmax_nuc, lmin, &
lmin_nuc, npgf, npgf_nuc
LOGICAL :: grid_distributed, paw_atom, use_virial
REAL(KIND=dp) :: eps_rho_rspace, force_tmp(3), fscale, &
ra(3), rpgf0, zet0
REAL(KIND=dp), DIMENSION(3, 3) :: my_virial_a, my_virial_b
REAL(KIND=dp), DIMENSION(:), POINTER :: hab_sph, norm_l, Qlm
REAL(KIND=dp), DIMENSION(:, :), POINTER :: hab, hdab_sph, intloc, pab
REAL(KIND=dp), DIMENSION(:, :), POINTER :: hab, hdab_sph, intloc, intloc_nuc, pab
REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: a_hdab_sph, hdab, Qlm_gg
REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: a_hdab
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(cell_type), POINTER :: cell
TYPE(cneo_potential_type), POINTER :: cneo_potential
TYPE(dft_control_type), POINTER :: dft_control
TYPE(gto_basis_set_type), POINTER :: basis_1c_set
TYPE(gto_basis_set_type), POINTER :: basis_1c_set, nuc_basis_set
TYPE(harmonics_atom_type), POINTER :: harmonics
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(pw_c1d_gs_type) :: coeff_gaux, coeff_gspace
@ -382,6 +393,8 @@ CONTAINS
TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: int_local_h, int_local_s
TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
TYPE(rho_atom_type), POINTER :: rho_atom
TYPE(rhoz_cneo_type), DIMENSION(:), POINTER :: rhoz_cneo_set
TYPE(rhoz_cneo_type), POINTER :: rhoz_cneo
TYPE(virial_type), POINTER :: virial
CALL timeset(routineN, handle)
@ -396,7 +409,7 @@ CONTAINS
! END IF
NULLIFY (atomic_kind_set, qs_kind_set, dft_control, particle_set)
NULLIFY (cell, force, pw_env, rho0_mpole, rho_atom_set)
NULLIFY (cell, force, pw_env, rho0_mpole, rho_atom_set, rhoz_cneo_set)
CALL get_qs_env(qs_env=qs_env, &
atomic_kind_set=atomic_kind_set, &
@ -406,6 +419,7 @@ CONTAINS
force=force, pw_env=pw_env, &
rho0_mpole=rho0_mpole, &
rho_atom_set=rho_atom_set, &
rhoz_cneo_set=rhoz_cneo_set, &
particle_set=particle_set, &
virial=virial)
@ -427,7 +441,8 @@ CONTAINS
!END IF
IF (PRESENT(local_rho_set)) &
CALL get_local_rho(local_rho_set, rho0_mpole=rho0_mpole, rho_atom_set=rho_atom_set)
CALL get_local_rho(local_rho_set, rho0_mpole=rho0_mpole, rho_atom_set=rho_atom_set, &
rhoz_cneo_set=rhoz_cneo_set)
! Q from rho0_mpole of local_rho_set
! for TDDFT forces we need mixed potential / integral space
! potential stored on local_rho_set_2nd
@ -514,12 +529,13 @@ CONTAINS
END IF
DO ikind = 1, SIZE(atomic_kind_set, 1)
NULLIFY (basis_1c_set, atom_list, harmonics)
NULLIFY (basis_1c_set, atom_list, harmonics, cneo_potential)
CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=nat)
CALL get_qs_kind(qs_kind_set(ikind), &
basis_set=basis_1c_set, basis_type="GAPW_1C", &
paw_atom=paw_atom, &
harmonics=harmonics)
harmonics=harmonics, &
cneo_potential=cneo_potential)
IF (.NOT. paw_atom) CYCLE
@ -543,7 +559,29 @@ CONTAINS
max_s_harm = harmonics%max_s_harm
llmax = harmonics%llmax
ALLOCATE (cg_list(2, nsoset(maxl)**2, max_s_harm), cg_n_list(max_s_harm))
NULLIFY (intloc_nuc)
maxl_nuc = -1
max_s_harm_nuc = 0
llmax_nuc = -1
IF (ASSOCIATED(cneo_potential)) THEN
NULLIFY (nuc_basis_set)
CALL get_qs_kind(qs_kind_set(ikind), &
basis_set=nuc_basis_set, &
basis_type="NUC")
CALL get_gto_basis_set(gto_basis_set=nuc_basis_set, &
lmax=lmax_nuc, lmin=lmin_nuc, &
maxso=maxso_nuc, maxl=maxl_nuc, &
nset=nset_nuc, npgf=npgf_nuc)
nsotot_nuc = maxso_nuc*nset_nuc
ALLOCATE (intloc_nuc(nsotot_nuc, nsotot_nuc))
max_s_harm_nuc = cneo_potential%harmonics%max_s_harm
llmax_nuc = cneo_potential%harmonics%llmax
END IF
ALLOCATE (cg_list(2, nsoset(MAX(maxl, maxl_nuc))**2, &
MAX(max_s_harm, max_s_harm_nuc)), &
cg_n_list(MAX(max_s_harm, max_s_harm_nuc)))
num_pe = para_env%num_pe
mepos = para_env%mepos
@ -645,9 +683,49 @@ CONTAINS
m1 = m1 + maxso
END DO ! iset1
IF (ASSOCIATED(cneo_potential)) THEN
intloc_nuc = 0.0_dp
m1 = 0
DO iset1 = 1, nset_nuc
n1 = nsoset(lmax_nuc(iset1))
m2 = 0
DO iset2 = 1, nset_nuc
n2 = nsoset(lmax_nuc(iset2))
CALL get_none0_cg_list(cneo_potential%harmonics%my_CG, lmin_nuc(iset1), &
lmax_nuc(iset1), lmin_nuc(iset2), lmax_nuc(iset2), &
max_s_harm_nuc, llmax_nuc, cg_list, cg_n_list, &
max_iso_not0_local)
DO ipgf1 = 1, npgf_nuc(iset1)
DO ipgf2 = 1, npgf_nuc(iset2)
DO iso = 1, MIN(nsoset(l0_ikind), max_iso_not0_local)
DO icg = 1, cg_n_list(iso)
iso1 = cg_list(1, icg, iso)
iso2 = cg_list(2, icg, iso)
ig1 = iso1 + n1*(ipgf1 - 1) + m1
ig2 = iso2 + n2*(ipgf2 - 1) + m2
intloc_nuc(ig1, ig2) = intloc_nuc(ig1, ig2) - cneo_potential%zeff* &
cneo_potential%Qlm_gg(ig1, ig2, iso)*hab_sph(iso)
END DO ! icg
END DO ! iso
END DO ! ipgf2
END DO ! ipgf1
m2 = m2 + maxso_nuc
END DO ! iset2
m1 = m1 + maxso_nuc
END DO ! iset1
END IF
IF (grid_distributed) THEN
! Sum result over all processors
CALL para_env%sum(intloc)
IF (ASSOCIATED(cneo_potential)) THEN
CALL para_env%sum(intloc_nuc)
END IF
END IF
IF (j == mepos) THEN
@ -657,6 +735,11 @@ CONTAINS
int_local_h(ispin)%r_coef = int_local_h(ispin)%r_coef + intloc
int_local_s(ispin)%r_coef = int_local_s(ispin)%r_coef + intloc
END DO
IF (ASSOCIATED(cneo_potential)) THEN
rhoz_cneo => rhoz_cneo_set(iatom)
rhoz_cneo%ga_Vlocal_gb_h = rhoz_cneo%ga_Vlocal_gb_h + intloc_nuc
rhoz_cneo%ga_Vlocal_gb_s = rhoz_cneo%ga_Vlocal_gb_s + intloc_nuc
END IF
END IF
IF (PRESENT(atener)) THEN
@ -691,6 +774,7 @@ CONTAINS
END DO
DEALLOCATE (intloc)
IF (ASSOCIATED(intloc_nuc)) DEALLOCATE (intloc_nuc)
DEALLOCATE (cg_list, cg_n_list)
END DO ! ikind
@ -701,6 +785,16 @@ CONTAINS
CALL timestop(handle)
IF (rho0_mpole%do_cneo) THEN
IF (PRESENT(kforce)) THEN
CALL integrate_vhgg_rspace(qs_env, v_rspace, para_env, calculate_forces, &
rhoz_cneo_set, kforce)
ELSE
CALL integrate_vhgg_rspace(qs_env, v_rspace, para_env, calculate_forces, &
rhoz_cneo_set)
END IF
END IF
END SUBROUTINE integrate_vhg0_rspace
END MODULE qs_rho0_ggrid

243
src/qs_rho0_methods.F
View file

@ -35,12 +35,18 @@ MODULE qs_rho0_methods
nso,&
nsoset
USE orbital_transformation_matrices, ONLY: orbtramat
USE qs_cneo_methods, ONLY: allocate_rhoz_cneo_internals,&
init_cneo_potential_internals
USE qs_cneo_types, ONLY: cneo_potential_type,&
rhoz_cneo_type
USE qs_cneo_utils, ONLY: cneo_scatter
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
USE qs_grid_atom, ONLY: grid_atom_type
USE qs_harmonics_atom, ONLY: get_none0_cg_list,&
harmonics_atom_type
USE qs_kind_types, ONLY: get_qs_kind,&
get_qs_kind_set,&
qs_kind_type,&
set_qs_kind
USE qs_local_rho_types, ONLY: allocate_rhoz,&
@ -74,15 +80,15 @@ CONTAINS
! **************************************************************************************************
!> \brief ...
!> \param mp_gau ...
!> \param Qlm_gg ...
!> \param basis_1c ...
!> \param harmonics ...
!> \param nchannels ...
!> \param nsotot ...
! **************************************************************************************************
SUBROUTINE calculate_mpole_gau(mp_gau, basis_1c, harmonics, nchannels, nsotot)
SUBROUTINE calculate_mpole_gau(Qlm_gg, basis_1c, harmonics, nchannels, nsotot)
TYPE(mpole_gau_overlap) :: mp_gau
REAL(dp), DIMENSION(:, :, :), POINTER :: Qlm_gg
TYPE(gto_basis_set_type), POINTER :: basis_1c
TYPE(harmonics_atom_type), POINTER :: harmonics
INTEGER, INTENT(IN) :: nchannels, nsotot
@ -102,7 +108,8 @@ CONTAINS
NULLIFY (lmax, lmin, npgf, my_CG, zet)
CALL reallocate(mp_gau%Qlm_gg, 1, nsotot, 1, nsotot, 1, nchannels)
CALL reallocate(Qlm_gg, 1, nsotot, 1, nsotot, 1, &
MIN(nchannels, harmonics%max_iso_not0))
CALL get_gto_basis_set(gto_basis_set=basis_1c, &
lmax=lmax, lmin=lmin, maxso=maxso, &
@ -142,8 +149,8 @@ CONTAINS
ig1 = iso1 + n1*(ipgf1 - 1) + m1
ig2 = iso2 + n2*(ipgf2 - 1) + m2
mp_gau%Qlm_gg(ig1, ig2, iso) = fourpi/(2._dp*l + 1._dp)* &
my_CG(iso1, iso2, iso)*gaussint_sph(zet1 + zet2, l + l1 + l2)
Qlm_gg(ig1, ig2, iso) = fourpi/(2._dp*l + 1._dp)* &
my_CG(iso1, iso2, iso)*gaussint_sph(zet1 + zet2, l + l1 + l2)
END DO ! icg
END DO ! iso
@ -163,6 +170,7 @@ CONTAINS
!> \brief ...
!> \param gapw_control ...
!> \param rho_atom_set ...
!> \param rhoz_cneo_set ...
!> \param rho0_atom_set ...
!> \param rho0_mp ...
!> \param a_list ...
@ -171,11 +179,12 @@ CONTAINS
!> \param qs_kind ...
!> \param rho0_h_tot ...
! **************************************************************************************************
SUBROUTINE calculate_rho0_atom(gapw_control, rho_atom_set, rho0_atom_set, &
SUBROUTINE calculate_rho0_atom(gapw_control, rho_atom_set, rhoz_cneo_set, rho0_atom_set, &
rho0_mp, a_list, natom, ikind, qs_kind, rho0_h_tot)
TYPE(gapw_control_type), POINTER :: gapw_control
TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
TYPE(rhoz_cneo_type), DIMENSION(:), POINTER :: rhoz_cneo_set
TYPE(rho0_atom_type), DIMENSION(:), POINTER :: rho0_atom_set
TYPE(rho0_mpole_type), POINTER :: rho0_mp
INTEGER, DIMENSION(:), INTENT(IN) :: a_list
@ -187,12 +196,14 @@ CONTAINS
INTEGER :: handle, iat, iatom, ic, ico, ir, is, &
iso, ispin, l, lmax0, lshell, lx, ly, &
lz, nr, nsotot, nspins
lz, nr, nsotot, nsotot_nuc, nspins
LOGICAL :: paw_atom
REAL(KIND=dp) :: sum1
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: cpc_h_nuc, cpc_s_nuc
REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: cpc_ah, cpc_as
REAL(KIND=dp), DIMENSION(:), POINTER :: norm_g0l_h
REAL(KIND=dp), DIMENSION(:, :), POINTER :: g0_h, vg0_h
TYPE(cneo_potential_type), POINTER :: cneo_potential
TYPE(grid_atom_type), POINTER :: g_atom
TYPE(harmonics_atom_type), POINTER :: harmonics
TYPE(mpole_gau_overlap), POINTER :: mpole_gau
@ -206,6 +217,7 @@ CONTAINS
NULLIFY (mpole_rho)
NULLIFY (g0_h, vg0_h, g_atom)
NULLIFY (norm_g0l_h, harmonics)
NULLIFY (cneo_potential)
CALL get_rho0_mpole(rho0_mpole=rho0_mp, ikind=ikind, &
l0_ikind=lmax0, mp_gau_ikind=mpole_gau, &
@ -213,7 +225,8 @@ CONTAINS
vg0_h=vg0_h, &
norm_g0l_h=norm_g0l_h)
CALL get_qs_kind(qs_kind, harmonics=harmonics, paw_atom=paw_atom, grid_atom=g_atom)
CALL get_qs_kind(qs_kind, harmonics=harmonics, paw_atom=paw_atom, grid_atom=g_atom, &
cneo_potential=cneo_potential)
nr = g_atom%nr
@ -222,7 +235,14 @@ CONTAINS
iatom = a_list(iat)
rho0_atom_set(iatom)%rho0_rad_h%r_coef = 0.0_dp
rho0_mp%mp_rho(iatom)%Qlm_tot = 0.0_dp
rho0_mp%mp_rho(iatom)%Qlm_tot(1) = rho0_mp%mp_rho(iatom)%Qlm_z
! When no nuclear density matrix is available, use spherical zeff
IF (.NOT. ASSOCIATED(cneo_potential)) THEN
rho0_mp%mp_rho(iatom)%Qlm_tot(1) = rho0_mp%mp_rho(iatom)%Qlm_z
ELSE
IF (.NOT. rhoz_cneo_set(iatom)%ready) THEN
rho0_mp%mp_rho(iatom)%Qlm_tot(1) = rho0_mp%mp_rho(iatom)%Qlm_z
END IF
END IF
rho0_mp%mp_rho(iatom)%Q0 = 0.0_dp
rho0_mp%mp_rho(iatom)%Qlm_car = 0.0_dp
END DO
@ -246,52 +266,86 @@ CONTAINS
CALL prj_scatter(cpc_h(ispin)%r_coef, cpc_ah(:, :, ispin), qs_kind)
CALL prj_scatter(cpc_s(ispin)%r_coef, cpc_as(:, :, ispin), qs_kind)
END DO
END IF
! Total charge (hard-soft) at atom
IF (paw_atom) THEN
DO ispin = 1, nspins
mpole_rho%Q0(ispin) = (trace_r_AxB(mpole_gau%Qlm_gg(:, :, 1), nsotot, &
cpc_ah(:, :, ispin), nsotot, nsotot, nsotot) &
- trace_r_AxB(mpole_gau%Qlm_gg(:, :, 1), nsotot, &
cpc_as(:, :, ispin), nsotot, nsotot, nsotot))/SQRT(fourpi)
END DO
END IF
! Multipoles of local charge distribution
DO iso = 1, nsoset(lmax0)
l = indso(1, iso)
IF (paw_atom) THEN
mpole_rho%Qlm_h(iso) = 0.0_dp
mpole_rho%Qlm_s(iso) = 0.0_dp
DO ispin = 1, nspins
mpole_rho%Qlm_h(iso) = mpole_rho%Qlm_h(iso) + &
trace_r_AxB(mpole_gau%Qlm_gg(:, :, iso), nsotot, &
cpc_ah(:, :, ispin), nsotot, nsotot, nsotot)
mpole_rho%Qlm_s(iso) = mpole_rho%Qlm_s(iso) + &
trace_r_AxB(mpole_gau%Qlm_gg(:, :, iso), nsotot, &
cpc_as(:, :, ispin), nsotot, nsotot, nsotot)
END DO ! ispin
mpole_rho%Qlm_tot(iso) = mpole_rho%Qlm_tot(iso) + &
mpole_rho%Qlm_h(iso) - mpole_rho%Qlm_s(iso)
nsotot_nuc = 0
IF (ASSOCIATED(cneo_potential)) THEN
IF (rhoz_cneo_set(iatom)%ready) THEN
nsotot_nuc = SIZE(cneo_potential%Qlm_gg, 1)
ALLOCATE (cpc_h_nuc(nsotot_nuc, nsotot_nuc), cpc_s_nuc(nsotot_nuc, nsotot_nuc))
cpc_h_nuc = 0._dp
cpc_s_nuc = 0._dp
CALL cneo_scatter(rhoz_cneo_set(iatom)%cpc_h, cpc_h_nuc, cneo_potential%npsgf, &
cneo_potential%n2oindex)
CALL cneo_scatter(rhoz_cneo_set(iatom)%cpc_s, cpc_s_nuc, cneo_potential%npsgf, &
cneo_potential%n2oindex)
END IF
END IF
! Total charge (hard-soft) at atom
IF (paw_atom) THEN
DO ispin = 1, nspins
mpole_rho%Q0(ispin) = (trace_r_AxB(mpole_gau%Qlm_gg(:, :, 1), nsotot, &
cpc_ah(:, :, ispin), nsotot, nsotot, nsotot) &
- trace_r_AxB(mpole_gau%Qlm_gg(:, :, 1), nsotot, &
cpc_as(:, :, ispin), nsotot, nsotot, nsotot))/SQRT(fourpi)
END DO
END IF
! Multipoles of local charge distribution
DO iso = 1, MIN(nsoset(lmax0), harmonics%max_iso_not0)
IF (paw_atom) THEN
mpole_rho%Qlm_h(iso) = 0.0_dp
mpole_rho%Qlm_s(iso) = 0.0_dp
DO ispin = 1, nspins
mpole_rho%Qlm_h(iso) = mpole_rho%Qlm_h(iso) + &
trace_r_AxB(mpole_gau%Qlm_gg(:, :, iso), nsotot, &
cpc_ah(:, :, ispin), nsotot, nsotot, nsotot)
mpole_rho%Qlm_s(iso) = mpole_rho%Qlm_s(iso) + &
trace_r_AxB(mpole_gau%Qlm_gg(:, :, iso), nsotot, &
cpc_as(:, :, ispin), nsotot, nsotot, nsotot)
END DO ! ispin
mpole_rho%Qlm_tot(iso) = mpole_rho%Qlm_tot(iso) + &
mpole_rho%Qlm_h(iso) - mpole_rho%Qlm_s(iso)
END IF
END DO ! iso
! Multipoles of CNEO quantum nuclear charge distribuition
IF (ASSOCIATED(cneo_potential)) THEN
IF (rhoz_cneo_set(iatom)%ready) THEN
DO iso = 1, MIN(nsoset(lmax0), cneo_potential%harmonics%max_iso_not0)
mpole_rho%Qlm_tot(iso) = mpole_rho%Qlm_tot(iso) - cneo_potential%zeff* &
trace_r_AxB(cneo_potential%Qlm_gg(:, :, iso), nsotot_nuc, &
cpc_h_nuc - cpc_s_nuc, nsotot_nuc, &
nsotot_nuc, nsotot_nuc)
END DO ! iso
END IF
END IF
DEALLOCATE (cpc_ah, cpc_as)
IF (ALLOCATED(cpc_h_nuc)) DEALLOCATE (cpc_h_nuc)
IF (ALLOCATED(cpc_s_nuc)) DEALLOCATE (cpc_s_nuc)
END IF
DO iso = 1, nsoset(lmax0)
l = indso(1, iso)
rho0_atom_set(iatom)%rho0_rad_h%r_coef(1:nr, iso) = &
g0_h(1:nr, l)*mpole_rho%Qlm_tot(iso)
rho0_atom_set(iatom)%vrho0_rad_h%r_coef(1:nr, iso) = &
vg0_h(1:nr, l)*mpole_rho%Qlm_tot(iso)
sum1 = 0.0_dp
DO ir = 1, nr
sum1 = sum1 + g_atom%wr(ir)* &
rho0_atom_set(iatom)%rho0_rad_h%r_coef(ir, iso)
END DO
rho0_h_tot = rho0_h_tot + sum1*harmonics%slm_int(iso)
! When CNEO is enabled, it is possible for rho0 to have a higher angualr momentum
! than that of the electronic density. In that case, the nuclear density must have
! a higher angular momentum, but cneo_potential%harmonics%slm_int is not initialized.
! For higher angular momenta, simply use the fact that slm_int(iso>1)=0
IF (iso <= harmonics%max_iso_not0) THEN
sum1 = 0.0_dp
DO ir = 1, nr
sum1 = sum1 + g_atom%wr(ir)* &
rho0_atom_set(iatom)%rho0_rad_h%r_coef(ir, iso)
END DO
rho0_h_tot = rho0_h_tot + sum1*harmonics%slm_int(iso)
END IF
END DO ! iso
IF (paw_atom) THEN
DEALLOCATE (cpc_ah, cpc_as)
END IF
END DO ! iat
END IF
@ -354,21 +408,23 @@ CONTAINS
CHARACTER(len=*), PARAMETER :: routineN = 'init_rho0'
CHARACTER(LEN=default_string_length) :: unit_str
INTEGER :: handle, iat, iatom, ikind, l, l_rho1_max, laddg, lmaxg, maxl, maxnset, maxso, &
nat, natom, nchan_c, nchan_s, nkind, nr, nset, nsotot, output_unit
INTEGER :: handle, iat, iatom, ikind, l, l_rho1_max, l_rho1_max_nuc, laddg, lmaxg, maxl, &
maxl_nuc, maxnset, maxso, maxso_nuc, nat, natom, nchan_c, nchan_s, nkind, nr, nset, &
nset_nuc, nsotot, nsotot_nuc, output_unit
INTEGER, DIMENSION(:), POINTER :: atom_list
LOGICAL :: paw_atom
REAL(KIND=dp) :: alpha_core, eps_Vrho0, max_rpgf0_s, &
radius, rc_min, rc_orb, &
total_rho_core_rspace, zeff
LOGICAL :: cneo_potential_present, paw_atom
REAL(KIND=dp) :: alpha_core, eps_Vrho0, max_rpgf0_s, radius, rc_min, rc_orb, &
total_rho_core_rspace, total_rho_nuc_cneo_rspace, zeff
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(cneo_potential_type), POINTER :: cneo_potential
TYPE(cp_logger_type), POINTER :: logger
TYPE(grid_atom_type), POINTER :: grid_atom
TYPE(gto_basis_set_type), POINTER :: basis_1c
TYPE(gto_basis_set_type), POINTER :: basis_1c, nuc_basis, nuc_soft_basis
TYPE(harmonics_atom_type), POINTER :: harmonics
TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
TYPE(rho0_atom_type), DIMENSION(:), POINTER :: rho0_atom_set
TYPE(rho0_mpole_type), POINTER :: rho0_mpole
TYPE(rhoz_cneo_type), DIMENSION(:), POINTER :: rhoz_cneo_set
TYPE(rhoz_type), DIMENSION(:), POINTER :: rhoz_set
TYPE(section_vals_type), POINTER :: dft_section
@ -384,6 +440,10 @@ CONTAINS
NULLIFY (rho0_mpole)
NULLIFY (rho0_atom_set)
NULLIFY (rhoz_set)
NULLIFY (cneo_potential)
NULLIFY (nuc_basis)
NULLIFY (nuc_soft_basis)
NULLIFY (rhoz_cneo_set)
CALL get_qs_env(qs_env=qs_env, qs_kind_set=qs_kind_set, &
atomic_kind_set=atomic_kind_set)
@ -394,6 +454,7 @@ CONTAINS
! Initialize rhoz total to zero
! in gapw rhoz is calculated on local the lebedev grids
total_rho_core_rspace = 0.0_dp
total_rho_nuc_cneo_rspace = 0.0_dp
CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
@ -425,10 +486,24 @@ CONTAINS
paw_atom=paw_atom, &
hard0_radius=rc_orb, &
zeff=zeff, &
alpha_core_charge=alpha_core)
cneo_potential=cneo_potential)
CALL get_qs_kind(qs_kind_set(ikind), &
basis_set=basis_1c, basis_type="GAPW_1C")
IF (ASSOCIATED(cneo_potential)) THEN
IF (PRESENT(zcore) .AND. zcore == 0.0_dp) &
CPABORT("Electronic TDDFT with CNEO quantum nuclei is not implemented.")
CPASSERT(paw_atom)
NULLIFY (nuc_basis, nuc_soft_basis)
CALL get_qs_kind(qs_kind_set(ikind), &
basis_set=nuc_basis, basis_type="NUC")
CALL get_qs_kind(qs_kind_set(ikind), &
basis_set=nuc_soft_basis, basis_type="NUC_SOFT")
alpha_core = 1.0_dp
ELSE
CALL get_qs_kind(qs_kind_set(ikind), alpha_core_charge=alpha_core)
END IF
! Set charge distribution of ionic cores to zero when computing the response-density
IF (PRESENT(zcore)) zeff = zcore
@ -439,8 +514,25 @@ CONTAINS
maxnset = MAX(maxnset, nset)
l_rho1_max = indso(1, harmonics%max_iso_not0)
maxl_nuc = -1
maxso_nuc = 0
nset_nuc = 0
l_rho1_max_nuc = -1
IF (ASSOCIATED(cneo_potential)) THEN
CALL get_gto_basis_set(gto_basis_set=nuc_basis, &
maxl=maxl_nuc, &
maxso=maxso_nuc, nset=nset_nuc)
! Initialize CNEO potential internals
CALL init_cneo_potential_internals(cneo_potential, nuc_basis, nuc_soft_basis, &
gapw_control, grid_atom)
l_rho1_max_nuc = indso(1, cneo_potential%harmonics%max_iso_not0)
END IF
IF (paw_atom) THEN
rho0_mpole%lmax0_kind(ikind) = MIN(2*maxl, l_rho1_max, maxl + laddg, lmaxg)
rho0_mpole%lmax0_kind(ikind) = MIN(2*MAX(maxl, maxl_nuc), &
MAX(l_rho1_max, l_rho1_max_nuc), &
MAX(maxl, maxl_nuc) + laddg, lmaxg)
ELSE
rho0_mpole%lmax0_kind(ikind) = 0
END IF
@ -453,6 +545,7 @@ CONTAINS
nchan_s = nsoset(rho0_mpole%lmax0_kind(ikind))
nchan_c = ncoset(rho0_mpole%lmax0_kind(ikind))
nsotot = maxso*nset
nsotot_nuc = maxso_nuc*nset_nuc
DO iat = 1, nat
iatom = atom_list(iat)
@ -465,18 +558,36 @@ CONTAINS
IF (paw_atom) THEN
! Calculate multipoles given by the product of 2 primitives Qlm_gg
CALL calculate_mpole_gau(rho0_mpole%mp_gau(ikind), &
CALL calculate_mpole_gau(rho0_mpole%mp_gau(ikind)%Qlm_gg, &
basis_1c, harmonics, nchan_s, nsotot)
END IF
! Calculate the core density rhoz
! exp(-alpha_c**2 r**2)Z(alpha_c**2/pi)**(3/2)
! on the logarithmic radial grid
! WARNING: alpha_core_charge = alpha_c**2
CALL calculate_rhoz(rhoz_set(ikind), grid_atom, alpha_core, zeff, &
nat, total_rho_core_rspace, harmonics)
IF (ASSOCIATED(cneo_potential)) THEN
rho0_mpole%do_cneo = .TRUE.
! Calculate multipoles given by the product of two nuclear primitives Qlm_gg
CALL calculate_mpole_gau(cneo_potential%Qlm_gg, nuc_basis, &
cneo_potential%harmonics, nchan_s, nsotot_nuc)
! initial CNEO quantum nuclear charge density is a simple Zeff sum,
! but it will be calculated from numerical integration during SCF
total_rho_nuc_cneo_rspace = total_rho_nuc_cneo_rspace - zeff*nat
ELSE
! Calculate the core density rhoz
! exp(-alpha_c**2 r**2)Z(alpha_c**2/pi)**(3/2)
! on the logarithmic radial grid
! WARNING: alpha_core_charge = alpha_c**2
CALL calculate_rhoz(rhoz_set(ikind), grid_atom, alpha_core, zeff, &
nat, total_rho_core_rspace, harmonics)
END IF
END DO ! ikind
total_rho_core_rspace = -total_rho_core_rspace
total_rho_nuc_cneo_rspace = -total_rho_nuc_cneo_rspace
! Allocate internals for quantum nuclear densities, if requested
CALL get_qs_kind_set(qs_kind_set, cneo_potential_present=cneo_potential_present)
IF (cneo_potential_present) THEN
CALL allocate_rhoz_cneo_internals(rhoz_cneo_set, atomic_kind_set, &
qs_kind_set, qs_env)
END IF
IF (gapw_control%alpha0_hard_from_input) THEN
! The exponent for the compensation charge rho0_hard is read from input
@ -509,8 +620,10 @@ CONTAINS
END DO
rho0_mpole%max_rpgf0_s = max_rpgf0_s
CALL set_local_rho(local_rho_set, rho0_atom_set=rho0_atom_set, rho0_mpole=rho0_mpole, rhoz_set=rhoz_set)
CALL set_local_rho(local_rho_set, rho0_atom_set=rho0_atom_set, rho0_mpole=rho0_mpole, &
rhoz_set=rhoz_set, rhoz_cneo_set=rhoz_cneo_set)
local_rho_set%rhoz_tot = total_rho_core_rspace
local_rho_set%rhoz_cneo_tot = total_rho_nuc_cneo_rspace
dft_section => section_vals_get_subs_vals(qs_env%input, "DFT")
output_unit = cp_print_key_unit_nr(logger, dft_section, "PRINT%GAPW%RHO0_INFORMATION", &

32
src/qs_rho0_types.F
View file

@ -61,6 +61,17 @@ MODULE qs_rho0_types
INTEGER :: lmax_0 = -1, igrid_zet0_s = -1
TYPE(pw_r3d_rs_type), POINTER :: rho0_s_rs => NULL()
TYPE(pw_c1d_gs_type), POINTER :: rho0_s_gs => NULL()
! CNEO nuclear charge density related stuff.
! These are put here because it is preferred that function rho0_s_grid_create
! can initialize PW for both rho0 and nuclear rho1s. This is akin to the behavior
! that init_rho0 also initializes data structures for CNEO nuclear charge densities,
! such that code modifications are mimimal for places where init_rho0 and
! rho0_s_grid_create are called.
LOGICAL :: do_cneo = .FALSE. ! if CNEO potential is used
REAL(dp) :: tot_rhoz_cneo_s = 0.0_dp ! soft nuclear charge on the local grid
TYPE(pw_r3d_rs_type), POINTER :: rhoz_cneo_s_rs => Null() ! soft nuclear charge density, r-space
TYPE(pw_c1d_gs_type), POINTER :: rhoz_cneo_s_gs => Null() ! soft nuclear charge density, g-space
END TYPE rho0_mpole_type
! **************************************************************************************************
@ -383,8 +394,18 @@ CONTAINS
IF (ASSOCIATED(rho0%rho0_s_gs)) THEN
CALL rho0%rho0_s_gs%release()
DEALLOCATE (rho0%rho0_s_gs)
END IF
IF (ASSOCIATED(rho0%rhoz_cneo_s_rs)) THEN
CALL rho0%rhoz_cneo_s_rs%release()
DEALLOCATE (rho0%rhoz_cneo_s_rs)
END IF
IF (ASSOCIATED(rho0%rhoz_cneo_s_gs)) THEN
CALL rho0%rhoz_cneo_s_gs%release()
DEALLOCATE (rho0%rhoz_cneo_s_gs)
END IF
DEALLOCATE (rho0)
ELSE
CALL cp_abort(__LOCATION__, &
@ -416,12 +437,15 @@ CONTAINS
!> \param max_rpgf0_s ...
!> \param rho0_s_rs ...
!> \param rho0_s_gs ...
!> \param rhoz_cneo_s_rs ...
!> \param rhoz_cneo_s_gs ...
! **************************************************************************************************
SUBROUTINE get_rho0_mpole(rho0_mpole, g0_h, vg0_h, iat, ikind, lmax_0, l0_ikind, &
mp_gau_ikind, mp_rho, norm_g0l_h, &
Qlm_gg, Qlm_car, Qlm_tot, &
zet0_h, igrid_zet0_s, rpgf0_h, rpgf0_s, &
max_rpgf0_s, rho0_s_rs, rho0_s_gs)
max_rpgf0_s, rho0_s_rs, rho0_s_gs, &
rhoz_cneo_s_rs, rhoz_cneo_s_gs)
TYPE(rho0_mpole_type), POINTER :: rho0_mpole
REAL(dp), DIMENSION(:, :), OPTIONAL, POINTER :: g0_h, vg0_h
@ -438,6 +462,8 @@ CONTAINS
REAL(dp), INTENT(OUT), OPTIONAL :: rpgf0_h, rpgf0_s, max_rpgf0_s
TYPE(pw_r3d_rs_type), OPTIONAL, POINTER :: rho0_s_rs
TYPE(pw_c1d_gs_type), OPTIONAL, POINTER :: rho0_s_gs
TYPE(pw_r3d_rs_type), OPTIONAL, POINTER :: rhoz_cneo_s_rs
TYPE(pw_c1d_gs_type), OPTIONAL, POINTER :: rhoz_cneo_s_gs
IF (ASSOCIATED(rho0_mpole)) THEN
@ -449,6 +475,8 @@ CONTAINS
IF (PRESENT(max_rpgf0_s)) max_rpgf0_s = rho0_mpole%max_rpgf0_s
IF (PRESENT(rho0_s_rs)) rho0_s_rs => rho0_mpole%rho0_s_rs
IF (PRESENT(rho0_s_gs)) rho0_s_gs => rho0_mpole%rho0_s_gs
IF (PRESENT(rhoz_cneo_s_rs)) rhoz_cneo_s_rs => rho0_mpole%rhoz_cneo_s_rs
IF (PRESENT(rhoz_cneo_s_gs)) rhoz_cneo_s_gs => rho0_mpole%rhoz_cneo_s_gs
IF (PRESENT(ikind)) THEN
IF (PRESENT(l0_ikind)) l0_ikind = rho0_mpole%lmax0_kind(ikind)

29
src/qs_scf_output.F
View file

@ -620,7 +620,9 @@ CONTAINS
accurate_sum(tot_rho_r) + nelectron_total, &
"Core density on regular grids:", &
qs_charges%total_rho_core_rspace, &
qs_charges%total_rho_core_rspace - REAL(nelectron_total + dft_control%charge, dp)
qs_charges%total_rho_core_rspace + &
qs_charges%total_rho1_hard_nuc - &
REAL(nelectron_total + dft_control%charge, dp)
IF (dft_control%correct_surf_dip) THEN
WRITE (UNIT=output_unit, FMT="((T3,A,/,T3,A,T41,F20.10))") &
@ -644,9 +646,24 @@ CONTAINS
accurate_sum(tot_rho_r) + tot1_h - tot1_s, &
"Total charge density (r-space): ", &
accurate_sum(tot_rho_r) + tot1_h - tot1_s &
+ qs_charges%total_rho_core_rspace, &
"Total Rho_soft + Rho0_soft (g-space):", &
qs_charges%total_rho_gspace
+ qs_charges%total_rho_core_rspace &
+ qs_charges%total_rho1_hard_nuc
IF (qs_charges%total_rho1_hard_nuc /= 0.0_dp) THEN
WRITE (UNIT=output_unit, FMT="(T3,A,T41,F20.10)") &
"Total CNEO nuc. char. den. (Lebedev): ", &
qs_charges%total_rho1_hard_nuc, &
"Total CNEO soft char. den. (Lebedev): ", &
qs_charges%total_rho1_soft_nuc_lebedev, &
"Total CNEO soft char. den. (r-space): ", &
qs_charges%total_rho1_soft_nuc_rspace, &
"Total soft Rho_e+n+0 (g-space):", &
qs_charges%total_rho_gspace
ELSE
WRITE (UNIT=output_unit, FMT="(T3,A,T41,F20.10)") &
"Total Rho_soft + Rho0_soft (g-space):", &
qs_charges%total_rho_gspace
END IF
! only add total_rho1_hard_nuc for gapw as cneo requires gapw
ELSE
WRITE (UNIT=output_unit, FMT="(T3,A,T41,F20.10)") &
"Total charge density on r-space grids: ", &
@ -779,6 +796,10 @@ CONTAINS
WRITE (UNIT=output_unit, FMT="(/,(T3,A,T56,F25.14))") &
"GAPW_XC| Exc from hard and soft atomic rho1: ", exc1_energy
END IF
IF (energy%core_cneo /= 0.0_dp) THEN
WRITE (UNIT=output_unit, FMT="(T3,A,T56,F25.14)") &
"CNEO| quantum nuclear core energy: ", energy%core_cneo
END IF
END IF
IF (dft_control%hairy_probes .EQV. .TRUE.) THEN
WRITE (UNIT=output_unit, FMT="((T3,A,T56,F25.14))") &

18
src/qs_scf_post_gpw.F
View file

@ -1908,7 +1908,7 @@ CONTAINS
TYPE(particle_list_type), POINTER :: particles
TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
TYPE(pw_c1d_gs_type) :: aux_g, rho_elec_gspace
TYPE(pw_c1d_gs_type), POINTER :: rho0_s_gs, rho_core
TYPE(pw_c1d_gs_type), POINTER :: rho0_s_gs, rho_core, rhoz_cneo_s_gs
TYPE(pw_env_type), POINTER :: pw_env
TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
@ -1955,7 +1955,7 @@ CONTAINS
! Print the total density (electronic + core charge)
IF (BTEST(cp_print_key_should_output(logger%iter_info, input, &
"DFT%PRINT%TOT_DENSITY_CUBE"), cp_p_file)) THEN
NULLIFY (rho_core, rho0_s_gs)
NULLIFY (rho_core, rho0_s_gs, rhoz_cneo_s_gs)
append_cube = section_get_lval(input, "DFT%PRINT%TOT_DENSITY_CUBE%APPEND")
my_pos_cube = "REWIND"
IF (append_cube) THEN
@ -1963,17 +1963,29 @@ CONTAINS
END IF
CALL get_qs_env(qs_env=qs_env, pw_env=pw_env, rho_core=rho_core, &
rho0_s_gs=rho0_s_gs)
rho0_s_gs=rho0_s_gs, rhoz_cneo_s_gs=rhoz_cneo_s_gs)
CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
pw_pools=pw_pools)
CALL auxbas_pw_pool%create_pw(wf_r)
IF (dft_control%qs_control%gapw) THEN
IF (dft_control%qs_control%gapw_control%nopaw_as_gpw) THEN
CALL pw_axpy(rho_core, rho0_s_gs)
IF (ASSOCIATED(rhoz_cneo_s_gs)) THEN
CALL pw_axpy(rhoz_cneo_s_gs, rho0_s_gs)
END IF
CALL pw_transfer(rho0_s_gs, wf_r)
CALL pw_axpy(rho_core, rho0_s_gs, -1.0_dp)
IF (ASSOCIATED(rhoz_cneo_s_gs)) THEN
CALL pw_axpy(rhoz_cneo_s_gs, rho0_s_gs, -1.0_dp)
END IF
ELSE
IF (ASSOCIATED(rhoz_cneo_s_gs)) THEN
CALL pw_axpy(rhoz_cneo_s_gs, rho0_s_gs)
END IF
CALL pw_transfer(rho0_s_gs, wf_r)
IF (ASSOCIATED(rhoz_cneo_s_gs)) THEN
CALL pw_axpy(rhoz_cneo_s_gs, rho0_s_gs, -1.0_dp)
END IF
END IF
ELSE
CALL pw_transfer(rho_core, wf_r)

9
src/qs_subsys_methods.F
View file

@ -37,6 +37,8 @@ MODULE qs_subsys_methods
USE molecule_kind_types, ONLY: get_molecule_kind,&
molecule_kind_type,&
set_molecule_kind
USE qs_cneo_types, ONLY: cneo_potential_type,&
get_cneo_potential
USE qs_kind_types, ONLY: create_qs_kind_set,&
get_qs_kind,&
init_atom_electronic_state,&
@ -170,6 +172,7 @@ CONTAINS
REAL(KIND=dp), DIMENSION(0:lmat, 10, 2) :: edelta
TYPE(all_potential_type), POINTER :: all_potential
TYPE(atomic_kind_type), POINTER :: atomic_kind
TYPE(cneo_potential_type), POINTER :: cneo_potential
TYPE(gth_potential_type), POINTER :: gth_potential
TYPE(gto_basis_set_type), POINTER :: orb_basis_set
TYPE(molecule_kind_type), POINTER :: molecule_kind
@ -201,7 +204,8 @@ CONTAINS
basis_set=orb_basis_set, &
all_potential=all_potential, &
gth_potential=gth_potential, &
sgp_potential=sgp_potential)
sgp_potential=sgp_potential, &
cneo_potential=cneo_potential)
! Obtain the electronic state of the atom
! The same state is used to calculate the ATOMIC GUESS
@ -237,6 +241,9 @@ CONTAINS
ELSE IF (ASSOCIATED(sgp_potential)) THEN
CALL get_potential(potential=sgp_potential, zeff=zeff, &
zeff_correction=zeff_correction)
ELSE IF (ASSOCIATED(cneo_potential)) THEN
CALL get_cneo_potential(potential=cneo_potential, zeff=zeff)
zeff_correction = 0.0_dp
ELSE
zeff = 0.0_dp
zeff_correction = 0.0_dp

3
src/qs_tddfpt2_fhxc_forces.F
View file

@ -432,6 +432,9 @@ CONTAINS
! calculate associated hartree potential
IF (gapw) THEN
CALL pw_axpy(local_rho_set%rho0_mpole%rho0_s_gs, rhox_tot_gspace)
IF (ASSOCIATED(local_rho_set%rho0_mpole%rhoz_cneo_s_gs)) THEN
CALL pw_axpy(local_rho_set%rho0_mpole%rhoz_cneo_s_gs, rhox_tot_gspace)
END IF
END IF
CALL pw_poisson_solve(poisson_env, rhox_tot_gspace, xehartree, &
xv_hartree_gspace)

3
src/qs_tddfpt2_forces.F
View file

@ -784,6 +784,9 @@ CONTAINS
CALL get_rho0_mpole(local_rho_set%rho0_mpole, Qlm_tot=Qlm_tot)
rhotot = rhotot + local_rho_set%rho0_mpole%total_rho0_h
CALL pw_axpy(local_rho_set%rho0_mpole%rho0_s_gs, rho_tot_gspace)
IF (ASSOCIATED(local_rho_set%rho0_mpole%rhoz_cneo_s_gs)) THEN
CALL pw_axpy(local_rho_set%rho0_mpole%rhoz_cneo_s_gs, rho_tot_gspace)
END IF
END IF
IF (ABS(rhotot) > 1.e-05_dp) THEN

3
src/qs_tddfpt2_operators.F
View file

@ -218,6 +218,9 @@ CONTAINS
IF (gapw) THEN
CPASSERT(ASSOCIATED(local_rho_set))
CALL pw_axpy(local_rho_set%rho0_mpole%rho0_s_gs, rho_ia_g)
IF (ASSOCIATED(local_rho_set%rho0_mpole%rhoz_cneo_s_gs)) THEN
CALL pw_axpy(local_rho_set%rho0_mpole%rhoz_cneo_s_gs, rho_ia_g)
END IF
END IF
CALL pw_poisson_solve(poisson_env, rho_ia_g, pair_energy, work_v_gspace)

3
src/qs_update_s_mstruct.F
View file

@ -85,6 +85,9 @@ CONTAINS
CALL get_qs_env(qs_env, rho_core=rho_core)
IF (dft_control%qs_control%gapw) THEN
qs_env%qs_charges%total_rho_core_rspace = qs_env%local_rho_set%rhoz_tot
! Initial CNEO quantum nuclear charge density is a simple Zeff sum.
! Later it will be calculated from numerical integration during SCF.
qs_env%qs_charges%total_rho1_hard_nuc = qs_env%local_rho_set%rhoz_cneo_tot
IF (dft_control%qs_control%gapw_control%nopaw_as_gpw) THEN
CPASSERT(ASSOCIATED(rho_core))
CALL calculate_rho_core(rho_core, &

9
src/response_solver.F
View file

@ -1136,6 +1136,9 @@ CONTAINS
CALL get_qs_env(qs_env, natom=natom)
! add rho0 contributions to GS density (only for Coulomb) only for gapw
CALL pw_axpy(local_rho_set_gs%rho0_mpole%rho0_s_gs, rho_tot_gspace_gs)
IF (ASSOCIATED(local_rho_set_gs%rho0_mpole%rhoz_cneo_s_gs)) THEN
CALL pw_axpy(local_rho_set_gs%rho0_mpole%rhoz_cneo_s_gs, rho_tot_gspace_gs)
END IF
IF (dft_control%qs_control%gapw_control%nopaw_as_gpw) THEN
CALL get_qs_env(qs_env=qs_env, rho_core=rho_core)
CALL pw_axpy(rho_core, rho_tot_gspace_gs)
@ -1495,6 +1498,9 @@ CONTAINS
! contribution for both T and D^Z
IF (gapw) THEN
CALL pw_axpy(local_rho_set_t%rho0_mpole%rho0_s_gs, rhoz_tot_gspace)
IF (ASSOCIATED(local_rho_set_t%rho0_mpole%rhoz_cneo_s_gs)) THEN
CALL pw_axpy(local_rho_set_t%rho0_mpole%rhoz_cneo_s_gs, rhoz_tot_gspace)
END IF
END IF
CALL pw_poisson_solve(poisson_env, rhoz_tot_gspace, ehartree, zv_hartree_gspace)
@ -1590,6 +1596,9 @@ CONTAINS
! add rho0 contributions to response density (only for Coulomb) only for gapw
IF (gapw) THEN
CALL pw_axpy(local_rho_set_t%rho0_mpole%rho0_s_gs, rho_tot_gspace_t)
IF (ASSOCIATED(local_rho_set_t%rho0_mpole%rhoz_cneo_s_gs)) THEN
CALL pw_axpy(local_rho_set_t%rho0_mpole%rhoz_cneo_s_gs, rho_tot_gspace_t)
END IF
! compute response Coulomb potential
CALL auxbas_pw_pool%create_pw(v_hartree_gspace_t)
CALL auxbas_pw_pool%create_pw(v_hartree_rspace_t)

17
src/surface_dipole.F
View file

@ -72,7 +72,7 @@ CONTAINS
REAL(dp), ALLOCATABLE, DIMENSION(:) :: rhoavsurf
TYPE(cell_type), POINTER :: cell
TYPE(dft_control_type), POINTER :: dft_control
TYPE(pw_c1d_gs_type), POINTER :: rho0_s_gs, rho_core
TYPE(pw_c1d_gs_type), POINTER :: rho0_s_gs, rho_core, rhoz_cneo_s_gs
TYPE(pw_env_type), POINTER :: pw_env
TYPE(pw_pool_p_type), DIMENSION(:), POINTER :: pw_pools
TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
@ -84,13 +84,14 @@ CONTAINS
CALL timeset(routineN, handle)
NULLIFY (cell, dft_control, rho, pw_env, auxbas_pw_pool, &
pw_pools, subsys, v_hartree_rspace, rho_r)
pw_pools, subsys, v_hartree_rspace, rho_r, rhoz_cneo_s_gs)
CALL get_qs_env(qs_env, &
dft_control=dft_control, &
rho=rho, &
rho_core=rho_core, &
rho0_s_gs=rho0_s_gs, &
rhoz_cneo_s_gs=rhoz_cneo_s_gs, &
cell=cell, &
pw_env=pw_env, &
subsys=subsys, &
@ -104,10 +105,22 @@ CONTAINS
IF (dft_control%qs_control%gapw) THEN
IF (dft_control%qs_control%gapw_control%nopaw_as_gpw) THEN
CALL pw_axpy(rho_core, rho0_s_gs)
IF (ASSOCIATED(rhoz_cneo_s_gs)) THEN
CALL pw_axpy(rhoz_cneo_s_gs, rho0_s_gs)
END IF
CALL pw_transfer(rho0_s_gs, wf_r)
CALL pw_axpy(rho_core, rho0_s_gs, -1.0_dp)
IF (ASSOCIATED(rhoz_cneo_s_gs)) THEN
CALL pw_axpy(rhoz_cneo_s_gs, rho0_s_gs, -1.0_dp)
END IF
ELSE
IF (ASSOCIATED(rhoz_cneo_s_gs)) THEN
CALL pw_axpy(rhoz_cneo_s_gs, rho0_s_gs)
END IF
CALL pw_transfer(rho0_s_gs, wf_r)
IF (ASSOCIATED(rhoz_cneo_s_gs)) THEN
CALL pw_axpy(rhoz_cneo_s_gs, rho0_s_gs, -1.0_dp)
END IF
END IF
ELSE
CALL pw_transfer(rho_core, wf_r)

61
tests/QS/regtest-cneo/H-cneo.inp Normal file
View file

@ -0,0 +1,61 @@
&GLOBAL
PRINT_LEVEL MEDIUM
PROJECT H-cneo
RUN_TYPE ENERGY
&END GLOBAL
&FORCE_EVAL
METHOD Quickstep
&DFT
BASIS_SET_FILE_NAME EMSL_BASIS_SETS
BASIS_SET_FILE_NAME NUCLEAR_BASIS_SETS
LSD
&MGRID
CUTOFF 150
NGRIDS 1
&END MGRID
&QS
ALPHA0_HARD 10
EPSFIT 1.E-4
EPSISO 1.0E-12
EPSRHO0 1.E-8
EPSSVD 0.0
EPS_GVG 1.0E-6
EPS_PGF_ORB 1.0E-6
LMAXN0 2
LMAXN1 6
METHOD GAPW
QUADRATURE GC_LOG
&END QS
&SCF
EPS_SCF 1.0E-4
IGNORE_CONVERGENCE_FAILURE
MAX_SCF 2
SCF_GUESS ATOMIC
&END SCF
&XC
&XC_FUNCTIONAL
&BECKE88
&END BECKE88
&LYP
&END LYP
&END XC_FUNCTIONAL
&END XC
&END DFT
&SUBSYS
&CELL
ABC 6.0 6.0 6.0
&END CELL
&COORD
H 0.0 0.0 0.0
&END COORD
&KIND H
BASIS_SET Ahlrichs-def2-SVP
BASIS_SET NUC PB4-D
HARD_EXP_RADIUS 1.5117809071370015011
LEBEDEV_GRID 50
POTENTIAL CNEO
RADIAL_GRID 50
&END KIND
&END SUBSYS
&END FORCE_EVAL

62
tests/QS/regtest-cneo/H2O-cneo-all.inp Normal file
View file

@ -0,0 +1,62 @@
&GLOBAL
PRINT_LEVEL MEDIUM
PROJECT H2O-cneo-all
&END GLOBAL
&FORCE_EVAL
METHOD Quickstep
&DFT
BASIS_SET_FILE_NAME BASIS_SET
BASIS_SET_FILE_NAME NUCLEAR_BASIS_SETS
POTENTIAL_FILE_NAME POTENTIAL
&MGRID
CUTOFF 150
NGRIDS 1
&END MGRID
&QS
ALPHA0_H 10
EPSFIT 1.E-4
EPSISO 1.0E-12
EPSRHO0 1.E-8
EPSSVD 0.0
EPS_GVG 1.0E-6
EPS_PGF_ORB 1.0E-6
LMAXN0 2
LMAXN1 6
METHOD GAPW
QUADRATURE GC_LOG
&END QS
&SCF
EPS_SCF 1.0E-4
SCF_GUESS ATOMIC
&END SCF
&XC
&XC_FUNCTIONAL Pade
&END XC_FUNCTIONAL
&END XC
&END DFT
&SUBSYS
&CELL
ABC 6.0 6.0 6.0
&END CELL
&COORD
O 0.000000 0.000000 -0.065587
H 0.000000 -0.757136 0.520545
H 0.000000 0.757136 0.520545
&END COORD
&KIND H
BASIS_SET TZVP-ALL-PADE
BASIS_SET NUC PB4-D
HARD_EXP_RADIUS 1.5117809071370015011
LEBEDEV_GRID 50
POTENTIAL CNEO
RADIAL_GRID 50
&END KIND
&KIND O
BASIS_SET TZVP-ALL-PADE
LEBEDEV_GRID 50
POTENTIAL ALL
RADIAL_GRID 50
&END KIND
&END SUBSYS
&END FORCE_EVAL

62
tests/QS/regtest-cneo/H2O-cneo-gth.inp Normal file
View file

@ -0,0 +1,62 @@
&GLOBAL
PRINT_LEVEL MEDIUM
PROJECT H2O-cneo-gth
&END GLOBAL
&FORCE_EVAL
METHOD Quickstep
&DFT
BASIS_SET_FILE_NAME BASIS_SET
BASIS_SET_FILE_NAME NUCLEAR_BASIS_SETS
POTENTIAL_FILE_NAME POTENTIAL
&MGRID
CUTOFF 150
NGRIDS 1
&END MGRID
&QS
ALPHA0_H 10
EPSFIT 1.E-4
EPSISO 1.0E-12
EPSRHO0 1.E-8
EPSSVD 0.0
EPS_GVG 1.0E-6
EPS_PGF_ORB 1.0E-6
LMAXN0 2
LMAXN1 6
METHOD GAPW
QUADRATURE GC_LOG
&END QS
&SCF
EPS_SCF 1.0E-4
SCF_GUESS ATOMIC
&END SCF
&XC
&XC_FUNCTIONAL Pade
&END XC_FUNCTIONAL
&END XC
&END DFT
&SUBSYS
&CELL
ABC 6.0 6.0 6.0
&END CELL
&COORD
O 0.000000 0.000000 -0.065587
H 0.000000 -0.757136 0.520545
H 0.000000 0.757136 0.520545
&END COORD
&KIND H
BASIS_SET TZVP-ALL-PADE
BASIS_SET NUC PB4-D
HARD_EXP_RADIUS 1.5117809071370015011
LEBEDEV_GRID 50
POTENTIAL CNEO
RADIAL_GRID 50
&END KIND
&KIND O
BASIS_SET DZVP-GTH-PADE
LEBEDEV_GRID 50
POTENTIAL GTH-PADE-q6
RADIAL_GRID 50
&END KIND
&END SUBSYS
&END FORCE_EVAL

55
tests/QS/regtest-cneo/Mu2S-cneo.inp Normal file
View file

@ -0,0 +1,55 @@
&GLOBAL
PRINT_LEVEL MEDIUM
PROJECT Mu2S-cneo
&END GLOBAL
&FORCE_EVAL
METHOD Quickstep
&DFT
BASIS_SET_FILE_NAME GTH_BASIS_SETS
BASIS_SET_FILE_NAME EMSL_BASIS_SETS
BASIS_SET_FILE_NAME NUCLEAR_BASIS_SETS
POTENTIAL_FILE_NAME POTENTIAL
&MGRID
CUTOFF 240
&END MGRID
&QS
EPS_DEFAULT 1.0E-8
EPS_SVD 0.0
EXTRAPOLATION PS
EXTRAPOLATION_ORDER 3
METHOD GAPW
&END QS
&SCF
EPS_SCF 1.0E-5
IGNORE_CONVERGENCE_FAILURE
MAX_SCF 2
SCF_GUESS ATOMIC
&END SCF
&XC
&XC_FUNCTIONAL BLYP
&END XC_FUNCTIONAL
&END XC
&END DFT
&SUBSYS
&CELL
ABC 4.0 4.0 4.0
&END CELL
&COORD
S 0.0 0.0 0.104467
H 0.0 0.975586 -0.835734
H 0.0 -0.975586 -0.835734
&END COORD
&KIND H
BASIS_SET Ahlrichs-def2-SVP
BASIS_SET NUC PB4-D_Mu
HARD_EXP_RADIUS 1.5117809071370015011
MASS 0.1134289259
POTENTIAL CNEO
&END KIND
&KIND S
BASIS_SET TZV2P-GTH
POTENTIAL GTH-BLYP-q6
&END KIND
&END SUBSYS
&END FORCE_EVAL

8
tests/QS/regtest-cneo/TEST_FILES.toml Normal file
View file

@ -0,0 +1,8 @@
# Hydrogen atom
"H-cneo.inp" = [{matcher="E_total", tol=1e-12, ref=-0.46201823434196}]
# test water with O ALL and GTH
"H2O-cneo-all.inp" = [{matcher="E_total", tol=1e-12, ref=-75.79965995506093}]
"H2O-cneo-gth.inp" = [{matcher="E_total", tol=1e-12, ref=-17.07942644605168}]
# test with presence of soft Muonium and S
"Mu2S-cneo.inp" = [{matcher="E_total", tol=1e-12, ref=-11.07849418175051}]
#EOF

1
tests/TEST_DIRS
View file

@ -394,3 +394,4 @@ QS/regtest-eht-guess libdftd4
QS/regtest-trexio trexio
QS/regtest-trexio-2 trexio libgrpp
QS/regtest-rtbse-gxac libint greenx
QS/regtest-cneo