Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Types needed for a for a Harris model calculation
10 : !> \par History
11 : !> 2024.07 created
12 : !> \author JGH
13 : ! **************************************************************************************************
14 : MODULE qs_harris_types
15 : USE atomic_kind_types, ONLY: atomic_kind_type,&
16 : get_atomic_kind_set
17 : USE basis_set_types, ONLY: get_gto_basis_set,&
18 : gto_basis_set_type
19 : USE distribution_1d_types, ONLY: distribution_1d_type
20 : USE kinds, ONLY: default_string_length,&
21 : dp
22 : USE pw_types, ONLY: pw_r3d_rs_type
23 : USE qs_kind_types, ONLY: get_qs_kind,&
24 : qs_kind_type
25 : #include "./base/base_uses.f90"
26 :
27 : IMPLICIT NONE
28 :
29 : PRIVATE
30 :
31 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_harris_types'
32 :
33 : ! *****************************************************************************
34 : TYPE rho_vec_type
35 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: rvecs
36 : END TYPE rho_vec_type
37 :
38 : TYPE harris_rhoin_type
39 : CHARACTER(LEN=default_string_length) :: basis_type = "NDef"
40 : TYPE(rho_vec_type), ALLOCATABLE, DIMENSION(:, :) :: rhovec
41 : TYPE(rho_vec_type), ALLOCATABLE, DIMENSION(:, :) :: intvec
42 : INTEGER :: nspin = 0
43 : INTEGER :: nbas = 0
44 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: basptr
45 : LOGICAL :: frozen = .FALSE.
46 : END TYPE harris_rhoin_type
47 :
48 : TYPE harris_energy_type
49 : REAL(KIND=dp) :: eharris = 0.0_dp
50 : REAL(KIND=dp) :: eband = 0.0_dp
51 : REAL(KIND=dp) :: exc_correction = 0.0_dp
52 : REAL(KIND=dp) :: eh_correction = 0.0_dp
53 : REAL(KIND=dp) :: ewald_correction = 0.0_dp
54 : REAL(KIND=dp) :: dispersion = 0.0_dp
55 : REAL(KIND=dp) :: trial_dm = 0.0_dp
56 : REAL(KIND=dp) :: direct_harris = 0.0_dp
57 : REAL(KIND=dp) :: direct_difference = 0.0_dp
58 : END TYPE harris_energy_type
59 :
60 : ! *****************************************************************************
61 : !> \brief Contains information on the Harris method
62 : !> \par History
63 : !> 07.2024 created
64 : !> \author JGH
65 : ! *****************************************************************************
66 : TYPE harris_type
67 : INTEGER :: energy_functional = 0
68 : INTEGER :: density_source = 0
69 : INTEGER :: orbital_basis = 0
70 : CHARACTER(LEN=default_string_length) :: density_filename = ""
71 : INTEGER :: fit_max_iter = 50
72 : INTEGER :: fit_max_backtrack = 20
73 : INTEGER :: fit_method = 0
74 : REAL(KIND=dp) :: fit_eps = 1.0E-3_dp
75 : REAL(KIND=dp) :: fit_step_size = 1.0_dp
76 : REAL(KIND=dp) :: fit_temperature = 0.0_dp
77 : REAL(KIND=dp) :: fit_relative_entropy_weight = 1.0E-3_dp
78 : LOGICAL :: density_fit_ready = .FALSE.
79 : LOGICAL :: density_target_ready = .FALSE.
80 : LOGICAL :: direct_density_matrix_energy = .FALSE.
81 : !
82 : TYPE(harris_energy_type) :: energy
83 : !
84 : TYPE(harris_rhoin_type) :: rhoin
85 : !
86 : TYPE(pw_r3d_rs_type) :: vh_rspace = pw_r3d_rs_type()
87 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: vxc_rspace => Null()
88 : TYPE(pw_r3d_rs_type) :: density_fit_rspace = pw_r3d_rs_type()
89 : TYPE(pw_r3d_rs_type) :: density_target_rspace = pw_r3d_rs_type()
90 :
91 : !
92 : LOGICAL :: debug_forces = .FALSE.
93 : LOGICAL :: debug_stress = .FALSE.
94 : END TYPE harris_type
95 : ! **************************************************************************************************
96 :
97 : PUBLIC :: harris_type, harris_energy_type, harris_env_release, &
98 : harris_print_direct_energy, harris_print_energy, harris_rhoin_type, harris_rhoin_init
99 :
100 : ! **************************************************************************************************
101 :
102 : CONTAINS
103 :
104 : ! **************************************************************************************************
105 :
106 : ! **************************************************************************************************
107 : !> \brief ...
108 : !> \param iounit ...
109 : !> \param energy ...
110 : ! **************************************************************************************************
111 56 : SUBROUTINE harris_print_energy(iounit, energy)
112 : INTEGER, INTENT(IN) :: iounit
113 : TYPE(harris_energy_type) :: energy
114 :
115 56 : IF (iounit > 0) THEN
116 28 : WRITE (UNIT=iounit, FMT="(/,(T2,A))") "HARRIS MODEL ENERGY INFORMATION"
117 : WRITE (UNIT=iounit, FMT="((T3,A,T56,F25.14))") &
118 28 : "Harris model energy: ", energy%eharris, &
119 28 : "Band energy: ", energy%eband, &
120 28 : "Hartree correction energy: ", energy%eh_correction, &
121 28 : "XC correction energy: ", energy%exc_correction, &
122 28 : "Ewald sum correction energy: ", energy%ewald_correction, &
123 56 : "Dispersion energy (pair potential): ", energy%dispersion
124 : END IF
125 :
126 56 : END SUBROUTINE harris_print_energy
127 :
128 : ! **************************************************************************************************
129 : !> \brief Prints the two direct fitted-density-matrix energy evaluations.
130 : !> \param iounit Output unit
131 : !> \param energy Harris energy data
132 : ! **************************************************************************************************
133 8 : SUBROUTINE harris_print_direct_energy(iounit, energy)
134 : INTEGER, INTENT(IN) :: iounit
135 : TYPE(harris_energy_type), INTENT(IN) :: energy
136 :
137 8 : IF (iounit > 0) THEN
138 4 : WRITE (UNIT=iounit, FMT="(/,(T2,A))") "HARRIS DIRECT DENSITY MATRIX ENERGY INFORMATION"
139 : WRITE (UNIT=iounit, FMT="((T3,A,T56,F25.14))") &
140 4 : "Consistent trial-DM energy: ", energy%trial_dm, &
141 4 : "Harris-like fitted-DM/cube energy: ", energy%direct_harris, &
142 8 : "Trial-DM minus Harris-like energy: ", energy%direct_difference
143 : END IF
144 :
145 8 : END SUBROUTINE harris_print_direct_energy
146 :
147 : ! **************************************************************************************************
148 : !> \brief ...
149 : !> \param rhoin ...
150 : !> \param basis_type ...
151 : !> \param qs_kind_set ...
152 : !> \param atomic_kind_set ...
153 : !> \param local_particles ...
154 : !> \param nspin ...
155 : ! **************************************************************************************************
156 8 : SUBROUTINE harris_rhoin_init(rhoin, basis_type, qs_kind_set, atomic_kind_set, &
157 : local_particles, nspin)
158 : TYPE(harris_rhoin_type) :: rhoin
159 : CHARACTER(LEN=*) :: basis_type
160 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
161 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
162 : TYPE(distribution_1d_type), POINTER :: local_particles
163 : INTEGER, INTENT(IN) :: nspin
164 :
165 : INTEGER :: iatom, ikind, iptr, ispin, natom, nkind, &
166 : nparticle_local, nsgf
167 8 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of, nbasf
168 : TYPE(gto_basis_set_type), POINTER :: basis_set
169 : TYPE(qs_kind_type), POINTER :: qs_kind
170 :
171 8 : CALL harris_rhoin_release(rhoin)
172 :
173 8 : rhoin%basis_type = basis_type
174 8 : rhoin%nspin = nspin
175 :
176 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, &
177 8 : atom_of_kind=atom_of_kind, kind_of=kind_of)
178 8 : natom = SIZE(atom_of_kind)
179 8 : nkind = SIZE(qs_kind_set)
180 :
181 24 : ALLOCATE (nbasf(nkind))
182 30 : DO ikind = 1, nkind
183 22 : qs_kind => qs_kind_set(ikind)
184 22 : CALL get_qs_kind(qs_kind, basis_set=basis_set, basis_type=basis_type)
185 22 : CALL get_gto_basis_set(basis_set, nsgf=nsgf)
186 30 : nbasf(ikind) = nsgf
187 : END DO
188 :
189 24 : ALLOCATE (rhoin%basptr(natom, 2))
190 8 : iptr = 1
191 44 : DO iatom = 1, natom
192 36 : ikind = kind_of(iatom)
193 36 : rhoin%basptr(iatom, 1) = iptr
194 36 : iptr = iptr + nbasf(ikind)
195 44 : rhoin%basptr(iatom, 2) = iptr - 1
196 : END DO
197 8 : rhoin%nbas = iptr - 1
198 :
199 70 : ALLOCATE (rhoin%rhovec(nkind, nspin))
200 18 : DO ispin = 1, nspin
201 46 : DO ikind = 1, nkind
202 28 : nsgf = nbasf(ikind)
203 28 : nparticle_local = local_particles%n_el(ikind)
204 114 : ALLOCATE (rhoin%rhovec(ikind, ispin)%rvecs(nsgf, nparticle_local))
205 : END DO
206 : END DO
207 :
208 62 : ALLOCATE (rhoin%intvec(nkind, nspin))
209 18 : DO ispin = 1, nspin
210 46 : DO ikind = 1, nkind
211 28 : nsgf = nbasf(ikind)
212 28 : nparticle_local = local_particles%n_el(ikind)
213 114 : ALLOCATE (rhoin%intvec(ikind, ispin)%rvecs(nsgf, nparticle_local))
214 : END DO
215 : END DO
216 :
217 8 : DEALLOCATE (nbasf)
218 :
219 16 : END SUBROUTINE harris_rhoin_init
220 :
221 : ! **************************************************************************************************
222 : !> \brief ...
223 : !> \param harris_env ...
224 : ! **************************************************************************************************
225 9094 : SUBROUTINE harris_env_release(harris_env)
226 : TYPE(harris_type), POINTER :: harris_env
227 :
228 : INTEGER :: iab
229 :
230 9094 : IF (ASSOCIATED(harris_env)) THEN
231 : !
232 9094 : CALL harris_rhoin_release(harris_env%rhoin)
233 : !
234 9094 : IF (ASSOCIATED(harris_env%vh_rspace%pw_grid)) THEN
235 28 : CALL harris_env%vh_rspace%release()
236 : END IF
237 9094 : IF (ASSOCIATED(harris_env%density_fit_rspace%pw_grid)) THEN
238 16 : CALL harris_env%density_fit_rspace%release()
239 : END IF
240 9094 : IF (ASSOCIATED(harris_env%density_target_rspace%pw_grid)) THEN
241 12 : CALL harris_env%density_target_rspace%release()
242 : END IF
243 9094 : IF (ASSOCIATED(harris_env%vxc_rspace)) THEN
244 58 : DO iab = 1, SIZE(harris_env%vxc_rspace)
245 58 : CALL harris_env%vxc_rspace(iab)%release()
246 : END DO
247 28 : DEALLOCATE (harris_env%vxc_rspace)
248 : END IF
249 : !
250 9094 : DEALLOCATE (harris_env)
251 : END IF
252 :
253 9094 : NULLIFY (harris_env)
254 :
255 9094 : END SUBROUTINE harris_env_release
256 :
257 : ! **************************************************************************************************
258 : !> \brief ...
259 : !> \param rhoin ...
260 : ! **************************************************************************************************
261 9102 : SUBROUTINE harris_rhoin_release(rhoin)
262 : TYPE(harris_rhoin_type) :: rhoin
263 :
264 : INTEGER :: i, j
265 :
266 9102 : IF (ALLOCATED(rhoin%rhovec)) THEN
267 18 : DO i = 1, SIZE(rhoin%rhovec, 2)
268 46 : DO j = 1, SIZE(rhoin%rhovec, 1)
269 38 : IF (ALLOCATED(rhoin%rhovec(j, i)%rvecs)) THEN
270 28 : DEALLOCATE (rhoin%rhovec(j, i)%rvecs)
271 : END IF
272 : END DO
273 : END DO
274 36 : DEALLOCATE (rhoin%rhovec)
275 : END IF
276 9102 : IF (ALLOCATED(rhoin%intvec)) THEN
277 18 : DO i = 1, SIZE(rhoin%intvec, 2)
278 46 : DO j = 1, SIZE(rhoin%intvec, 1)
279 38 : IF (ALLOCATED(rhoin%intvec(j, i)%rvecs)) THEN
280 28 : DEALLOCATE (rhoin%intvec(j, i)%rvecs)
281 : END IF
282 : END DO
283 : END DO
284 36 : DEALLOCATE (rhoin%intvec)
285 : END IF
286 9102 : IF (ALLOCATED(rhoin%basptr)) THEN
287 8 : DEALLOCATE (rhoin%basptr)
288 : END IF
289 9102 : rhoin%basis_type = "NDef"
290 9102 : rhoin%nspin = 0
291 9102 : rhoin%nbas = 0
292 9102 : rhoin%frozen = .FALSE.
293 :
294 9102 : END SUBROUTINE harris_rhoin_release
295 :
296 0 : END MODULE qs_harris_types
|
