-
Notifications
You must be signed in to change notification settings - Fork 1
/
qs_gapw_densities.F
182 lines (152 loc) · 7.97 KB
/
qs_gapw_densities.F
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
!--------------------------------------------------------------------------------------------------!
! CP2K: A general program to perform molecular dynamics simulations !
! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
! !
! SPDX-License-Identifier: GPL-2.0-or-later !
!--------------------------------------------------------------------------------------------------!
! **************************************************************************************************
MODULE qs_gapw_densities
USE atomic_kind_types, ONLY: atomic_kind_type,&
get_atomic_kind
USE cp_control_types, ONLY: dft_control_type,&
gapw_control_type
USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
USE kinds, ONLY: dp
USE message_passing, ONLY: mp_para_env_type
USE qs_charges_types, ONLY: qs_charges_type
USE qs_environment_types, ONLY: get_qs_env,&
qs_environment_type
USE qs_kind_types, ONLY: get_qs_kind,&
qs_kind_type
USE qs_local_rho_types, ONLY: local_rho_type
USE qs_rho0_ggrid, ONLY: put_rho0_on_grid
USE qs_rho0_methods, ONLY: calculate_rho0_atom
USE qs_rho0_types, ONLY: rho0_atom_type,&
rho0_mpole_type
USE qs_rho_atom_methods, ONLY: calculate_rho_atom
USE qs_rho_atom_types, ONLY: rho_atom_type
#include "./base/base_uses.f90"
IMPLICIT NONE
PRIVATE
CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_gapw_densities'
PUBLIC :: prepare_gapw_den
CONTAINS
! **************************************************************************************************
!> \brief ...
!> \param qs_env ...
!> \param local_rho_set ...
!> \param do_rho0 ...
!> \param kind_set_external can be provided to use different projectors/grids/basis than the default
! **************************************************************************************************
SUBROUTINE prepare_gapw_den(qs_env, local_rho_set, do_rho0, kind_set_external)
TYPE(qs_environment_type), POINTER :: qs_env
TYPE(local_rho_type), OPTIONAL, POINTER :: local_rho_set
LOGICAL, INTENT(IN), OPTIONAL :: do_rho0
TYPE(qs_kind_type), DIMENSION(:), OPTIONAL, &
POINTER :: kind_set_external
CHARACTER(len=*), PARAMETER :: routineN = 'prepare_gapw_den'
INTEGER :: handle, ikind, ispin, natom, nspins, &
output_unit
INTEGER, DIMENSION(:), POINTER :: atom_list
LOGICAL :: extern, my_do_rho0, paw_atom
REAL(dp) :: rho0_h_tot, tot_rs_int
REAL(dp), DIMENSION(:), POINTER :: rho1_h_tot, rho1_s_tot
TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
TYPE(dft_control_type), POINTER :: dft_control
TYPE(gapw_control_type), POINTER :: gapw_control
TYPE(mp_para_env_type), POINTER :: para_env
TYPE(qs_charges_type), POINTER :: qs_charges
TYPE(qs_kind_type), DIMENSION(:), POINTER :: my_kind_set
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
CALL timeset(routineN, handle)
NULLIFY (atomic_kind_set)
NULLIFY (my_kind_set)
NULLIFY (dft_control)
NULLIFY (gapw_control)
NULLIFY (para_env)
NULLIFY (atom_list)
NULLIFY (rho0_mpole)
NULLIFY (qs_charges)
NULLIFY (rho1_h_tot, rho1_s_tot)
NULLIFY (rho_atom_set)
NULLIFY (rho0_atom_set)
my_do_rho0 = .TRUE.
IF (PRESENT(do_rho0)) my_do_rho0 = do_rho0
output_unit = cp_logger_get_default_io_unit()
CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, &
para_env=para_env, &
qs_charges=qs_charges, &
qs_kind_set=my_kind_set, &
atomic_kind_set=atomic_kind_set, &
rho0_mpole=rho0_mpole, &
rho_atom_set=rho_atom_set, &
rho0_atom_set=rho0_atom_set)
gapw_control => dft_control%qs_control%gapw_control
IF (PRESENT(local_rho_set)) THEN
rho_atom_set => local_rho_set%rho_atom_set
IF (my_do_rho0) THEN
rho0_mpole => local_rho_set%rho0_mpole
rho0_atom_set => local_rho_set%rho0_atom_set
END IF
END IF
extern = .FALSE.
IF (PRESENT(kind_set_external)) THEN
CPASSERT(ASSOCIATED(kind_set_external))
my_kind_set => kind_set_external
extern = .TRUE.
END IF
nspins = dft_control%nspins
rho0_h_tot = 0.0_dp
ALLOCATE (rho1_h_tot(1:nspins), rho1_s_tot(1:nspins))
rho1_h_tot = 0.0_dp
rho1_s_tot = 0.0_dp
DO ikind = 1, SIZE(atomic_kind_set)
CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
CALL get_qs_kind(my_kind_set(ikind), paw_atom=paw_atom)
!Calculate rho1_h and rho1_s on the radial grids centered on the atomic position
IF (paw_atom) THEN
CALL calculate_rho_atom(para_env, rho_atom_set, my_kind_set(ikind), &
atom_list, natom, nspins, rho1_h_tot, rho1_s_tot)
END IF
!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)
END DO
!Do not mess with charges if using a non-default kind_set
IF (.NOT. extern) THEN
CALL para_env%sum(rho1_h_tot)
CALL para_env%sum(rho1_s_tot)
DO ispin = 1, nspins
qs_charges%total_rho1_hard(ispin) = -rho1_h_tot(ispin)
qs_charges%total_rho1_soft(ispin) = -rho1_s_tot(ispin)
END DO
IF (my_do_rho0) THEN
rho0_mpole%total_rho0_h = -rho0_h_tot
!Put the rho0_soft on the global grid
CALL put_rho0_on_grid(qs_env, rho0_mpole, tot_rs_int)
IF (ABS(rho0_h_tot) .GE. 1.0E-5_dp) THEN
IF (ABS(1.0_dp - ABS(tot_rs_int/rho0_h_tot)) .GT. 1.0E-3_dp) THEN
IF (output_unit > 0) THEN
WRITE (output_unit, '(/,72("*"))')
WRITE (output_unit, '(T2,A,T66,1E20.8)') &
"WARNING: rho0 calculated on the local grid is :", -rho0_h_tot, &
" rho0 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_rho0_soft_rspace = tot_rs_int
qs_charges%total_rho0_hard_lebedev = rho0_h_tot
ELSE
qs_charges%total_rho0_hard_lebedev = 0.0_dp
END IF
END IF
DEALLOCATE (rho1_h_tot, rho1_s_tot)
CALL timestop(handle)
END SUBROUTINE prepare_gapw_den
END MODULE qs_gapw_densities