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