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 used by CNEO-DFT
10 : !> (see J. Chem. Theory Comput. 2025, 21, 16, 7865–7877)
11 : !> \par History
12 : !> 08.2025 created [zc62]
13 : !> \author Zehua Chen
14 : ! **************************************************************************************************
15 : MODULE qs_cneo_types
16 : USE kinds, ONLY: dp
17 : USE periodic_table, ONLY: ptable
18 : USE qs_harmonics_atom, ONLY: deallocate_harmonics_atom,&
19 : harmonics_atom_type
20 : #include "./base/base_uses.f90"
21 :
22 : IMPLICIT NONE
23 :
24 : PRIVATE
25 :
26 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_cneo_types'
27 :
28 : ! Essential matrices, density and potential for each quantum nucleus
29 : TYPE rhoz_cneo_type
30 : LOGICAL :: ready = .FALSE. ! if pmat is ready. Useful in first iter
31 : REAL(dp), DIMENSION(:, :), &
32 : POINTER :: pmat => Null(), & ! nuclear density matrix
33 : core => Null(), & ! nuclear core Hamiltonian
34 : vmat => Null(), & ! nuclear Hartree from soft basis
35 : fmat => Null(), & ! Fock = core + Hartree
36 : wfn => Null() ! nuclear orbital coefficients
37 : REAL(dp), DIMENSION(3) &
38 : :: f = [0.0_dp, 0.0_dp, 0.0_dp] ! Lagrange multiplier for CNEO
39 : REAL(dp) :: e_core = 0.0_dp ! nuclear core energy
40 : REAL(dp), DIMENSION(:, :), &
41 : POINTER :: cpc_h => Null(), & ! decontracted density matrix
42 : cpc_s => Null(), & ! decontracted density matrix, soft tail
43 : rho_rad_h => Null(), & ! density on radial grid
44 : rho_rad_s => Null(), & ! density on radial grid, soft tail
45 : vrho_rad_h => Null(), & ! potential on radial grid
46 : vrho_rad_s => Null(), & ! potential on radial grid, soft tail
47 : ga_Vlocal_gb_h => Null(), & ! local Hartree integral
48 : ga_Vlocal_gb_s => Null() ! local Hartree integral, soft tail
49 : END TYPE rhoz_cneo_type
50 :
51 : ! S, T, transformation matrices and distance are shared by the same kind,
52 : ! since they are not affected by the position of basis center
53 : TYPE cneo_potential_type
54 : INTEGER :: z = 0 ! atomic number
55 : REAL(dp) :: zeff = 0.0_dp, & ! zeff = REAL(z)
56 : mass = 0.0_dp ! atomic mass in dalton
57 : INTEGER, DIMENSION(:), &
58 : POINTER :: elec_conf => Null()
59 : INTEGER :: nsgf = 0, & ! nuclear basis set is usually uncontracted
60 : nne = 0, & ! number of linear-independent basis functions
61 : npsgf = 0, & ! number of primitive SGFs
62 : nsotot = 0 ! maxso * nset
63 : REAL(dp), DIMENSION(:, :), &
64 : POINTER :: my_gcc_h => Null(), & ! 3D-normalized contraction coefficients
65 : my_gcc_s => Null(), & ! contraction coefficients for the soft tail
66 : ovlp => Null(), & ! nuclear basis overlap matrix (unused)
67 : kin => Null(), & ! nuclear kinetic energy matrix
68 : utrans => Null() ! nuclear basis transformation matrix
69 : REAL(dp), DIMENSION(:, :, :), &
70 : POINTER :: distance => Null() ! distance from nuclear basis center
71 : TYPE(harmonics_atom_type), &
72 : POINTER :: harmonics => Null() ! most of the data will be missing
73 : REAL(dp), DIMENSION(:, :, :), &
74 : POINTER :: Qlm_gg => Null(), & ! multipole expansion of nuclear gg
75 : gg => Null() ! precompute and store gg on radial grid
76 : REAL(dp), DIMENSION(:, :, :, :), &
77 : POINTER :: vgg => Null() ! precompute and store vgg on grid
78 : INTEGER, DIMENSION(:), &
79 : POINTER :: n2oindex => Null(), & ! new to old index
80 : o2nindex => Null() ! old to new index
81 : REAL(dp), DIMENSION(:, :), &
82 : POINTER :: rad2l => Null(), & ! store my own rad2l
83 : oorad2l => Null() ! store my own oorad2l
84 : END TYPE cneo_potential_type
85 :
86 : ! Public Types
87 :
88 : PUBLIC :: cneo_potential_type, rhoz_cneo_type
89 :
90 : ! Public Subroutine
91 :
92 : PUBLIC :: allocate_cneo_potential, allocate_rhoz_cneo_set, deallocate_cneo_potential, &
93 : deallocate_rhoz_cneo_set, get_cneo_potential, set_cneo_potential, write_cneo_potential
94 :
95 : CONTAINS
96 :
97 : ! **************************************************************************************************
98 : !> \brief ...
99 : !> \param rhoz_cneo_set ...
100 : !> \param natom ...
101 : ! **************************************************************************************************
102 8 : SUBROUTINE allocate_rhoz_cneo_set(rhoz_cneo_set, natom)
103 :
104 : TYPE(rhoz_cneo_type), DIMENSION(:), POINTER :: rhoz_cneo_set
105 : INTEGER, INTENT(IN) :: natom
106 :
107 8 : IF (ASSOCIATED(rhoz_cneo_set)) THEN
108 0 : CALL deallocate_rhoz_cneo_set(rhoz_cneo_set)
109 : END IF
110 :
111 68 : ALLOCATE (rhoz_cneo_set(natom))
112 :
113 8 : END SUBROUTINE allocate_rhoz_cneo_set
114 :
115 : ! **************************************************************************************************
116 : !> \brief ...
117 : !> \param rhoz_cneo ...
118 : ! **************************************************************************************************
119 20 : SUBROUTINE deallocate_rhoz_cneo(rhoz_cneo)
120 :
121 : TYPE(rhoz_cneo_type), POINTER :: rhoz_cneo
122 :
123 20 : IF (ASSOCIATED(rhoz_cneo)) THEN
124 20 : IF (ASSOCIATED(rhoz_cneo%pmat)) THEN
125 14 : DEALLOCATE (rhoz_cneo%pmat)
126 : END IF
127 20 : IF (ASSOCIATED(rhoz_cneo%core)) THEN
128 14 : DEALLOCATE (rhoz_cneo%core)
129 : END IF
130 20 : IF (ASSOCIATED(rhoz_cneo%vmat)) THEN
131 14 : DEALLOCATE (rhoz_cneo%vmat)
132 : END IF
133 20 : IF (ASSOCIATED(rhoz_cneo%fmat)) THEN
134 7 : DEALLOCATE (rhoz_cneo%fmat)
135 : END IF
136 20 : IF (ASSOCIATED(rhoz_cneo%wfn)) THEN
137 7 : DEALLOCATE (rhoz_cneo%wfn)
138 : END IF
139 20 : IF (ASSOCIATED(rhoz_cneo%cpc_h)) THEN
140 14 : DEALLOCATE (rhoz_cneo%cpc_h)
141 : END IF
142 20 : IF (ASSOCIATED(rhoz_cneo%cpc_s)) THEN
143 14 : DEALLOCATE (rhoz_cneo%cpc_s)
144 : END IF
145 20 : IF (ASSOCIATED(rhoz_cneo%rho_rad_h)) THEN
146 7 : DEALLOCATE (rhoz_cneo%rho_rad_h)
147 : END IF
148 20 : IF (ASSOCIATED(rhoz_cneo%rho_rad_s)) THEN
149 7 : DEALLOCATE (rhoz_cneo%rho_rad_s)
150 : END IF
151 20 : IF (ASSOCIATED(rhoz_cneo%vrho_rad_h)) THEN
152 7 : DEALLOCATE (rhoz_cneo%vrho_rad_h)
153 : END IF
154 20 : IF (ASSOCIATED(rhoz_cneo%vrho_rad_s)) THEN
155 7 : DEALLOCATE (rhoz_cneo%vrho_rad_s)
156 : END IF
157 20 : IF (ASSOCIATED(rhoz_cneo%ga_Vlocal_gb_h)) THEN
158 7 : DEALLOCATE (rhoz_cneo%ga_Vlocal_gb_h)
159 : END IF
160 20 : IF (ASSOCIATED(rhoz_cneo%ga_Vlocal_gb_s)) THEN
161 7 : DEALLOCATE (rhoz_cneo%ga_Vlocal_gb_s)
162 : END IF
163 : END IF
164 :
165 20 : END SUBROUTINE deallocate_rhoz_cneo
166 :
167 : ! **************************************************************************************************
168 : !> \brief ...
169 : !> \param rhoz_cneo_set ...
170 : ! **************************************************************************************************
171 8 : SUBROUTINE deallocate_rhoz_cneo_set(rhoz_cneo_set)
172 :
173 : TYPE(rhoz_cneo_type), DIMENSION(:), POINTER :: rhoz_cneo_set
174 :
175 : INTEGER :: iat, natom
176 : TYPE(rhoz_cneo_type), POINTER :: rhoz_cneo
177 :
178 8 : IF (ASSOCIATED(rhoz_cneo_set)) THEN
179 8 : natom = SIZE(rhoz_cneo_set)
180 28 : DO iat = 1, natom
181 20 : rhoz_cneo => rhoz_cneo_set(iat)
182 28 : CALL deallocate_rhoz_cneo(rhoz_cneo)
183 : END DO
184 8 : DEALLOCATE (rhoz_cneo_set)
185 : END IF
186 :
187 8 : END SUBROUTINE deallocate_rhoz_cneo_set
188 :
189 : ! **************************************************************************************************
190 : !> \brief ...
191 : !> \param potential ...
192 : ! **************************************************************************************************
193 8 : SUBROUTINE allocate_cneo_potential(potential)
194 :
195 : TYPE(cneo_potential_type), POINTER :: potential
196 :
197 8 : IF (ASSOCIATED(potential)) THEN
198 0 : CALL deallocate_cneo_potential(potential)
199 : END IF
200 :
201 8 : ALLOCATE (potential)
202 :
203 8 : END SUBROUTINE allocate_cneo_potential
204 :
205 : ! **************************************************************************************************
206 : !> \brief ...
207 : !> \param potential ...
208 : ! **************************************************************************************************
209 8 : SUBROUTINE deallocate_cneo_potential(potential)
210 :
211 : TYPE(cneo_potential_type), POINTER :: potential
212 :
213 8 : IF (ASSOCIATED(potential)) THEN
214 8 : IF (ASSOCIATED(potential%elec_conf)) THEN
215 8 : DEALLOCATE (potential%elec_conf)
216 : END IF
217 8 : IF (ASSOCIATED(potential%my_gcc_h)) THEN
218 8 : DEALLOCATE (potential%my_gcc_h)
219 : END IF
220 8 : IF (ASSOCIATED(potential%my_gcc_s)) THEN
221 8 : DEALLOCATE (potential%my_gcc_s)
222 : END IF
223 8 : IF (ASSOCIATED(potential%ovlp)) THEN
224 8 : DEALLOCATE (potential%ovlp)
225 : END IF
226 8 : IF (ASSOCIATED(potential%kin)) THEN
227 8 : DEALLOCATE (potential%kin)
228 : END IF
229 8 : IF (ASSOCIATED(potential%utrans)) THEN
230 8 : DEALLOCATE (potential%utrans)
231 : END IF
232 8 : IF (ASSOCIATED(potential%distance)) THEN
233 8 : DEALLOCATE (potential%distance)
234 : END IF
235 8 : IF (ASSOCIATED(potential%harmonics)) THEN
236 8 : CALL deallocate_harmonics_atom(potential%harmonics)
237 : END IF
238 8 : IF (ASSOCIATED(potential%Qlm_gg)) THEN
239 8 : DEALLOCATE (potential%Qlm_gg)
240 : END IF
241 8 : IF (ASSOCIATED(potential%gg)) THEN
242 8 : DEALLOCATE (potential%gg)
243 : END IF
244 8 : IF (ASSOCIATED(potential%vgg)) THEN
245 8 : DEALLOCATE (potential%vgg)
246 : END IF
247 8 : IF (ASSOCIATED(potential%n2oindex)) THEN
248 8 : DEALLOCATE (potential%n2oindex)
249 : END IF
250 8 : IF (ASSOCIATED(potential%o2nindex)) THEN
251 8 : DEALLOCATE (potential%o2nindex)
252 : END IF
253 8 : IF (ASSOCIATED(potential%rad2l)) THEN
254 8 : DEALLOCATE (potential%rad2l)
255 : END IF
256 8 : IF (ASSOCIATED(potential%oorad2l)) THEN
257 8 : DEALLOCATE (potential%oorad2l)
258 : END IF
259 8 : DEALLOCATE (potential)
260 : END IF
261 :
262 8 : END SUBROUTINE deallocate_cneo_potential
263 :
264 : ! **************************************************************************************************
265 : !> \brief ...
266 : !> \param potential ...
267 : !> \param z ...
268 : !> \param zeff ...
269 : !> \param mass ...
270 : !> \param elec_conf ...
271 : !> \param nsgf ...
272 : !> \param nne ...
273 : !> \param npsgf ...
274 : !> \param nsotot ...
275 : !> \param my_gcc_h ...
276 : !> \param my_gcc_s ...
277 : !> \param ovlp ...
278 : !> \param kin ...
279 : !> \param utrans ...
280 : !> \param distance ...
281 : !> \param harmonics ...
282 : !> \param Qlm_gg ...
283 : !> \param gg ...
284 : !> \param vgg ...
285 : !> \param n2oindex ...
286 : !> \param o2nindex ...
287 : !> \param rad2l ...
288 : !> \param oorad2l ...
289 : ! **************************************************************************************************
290 433 : SUBROUTINE get_cneo_potential(potential, z, zeff, mass, elec_conf, nsgf, nne, npsgf, &
291 : nsotot, my_gcc_h, my_gcc_s, ovlp, kin, utrans, distance, &
292 : harmonics, Qlm_gg, gg, vgg, n2oindex, o2nindex, rad2l, oorad2l)
293 :
294 : TYPE(cneo_potential_type), POINTER :: potential
295 : INTEGER, INTENT(OUT), OPTIONAL :: z
296 : REAL(dp), INTENT(OUT), OPTIONAL :: zeff, mass
297 : INTEGER, DIMENSION(:), OPTIONAL, POINTER :: elec_conf
298 : INTEGER, INTENT(OUT), OPTIONAL :: nsgf, nne, npsgf, nsotot
299 : REAL(dp), DIMENSION(:, :), OPTIONAL, POINTER :: my_gcc_h, my_gcc_s, ovlp, kin, utrans
300 : REAL(dp), DIMENSION(:, :, :), OPTIONAL, POINTER :: distance
301 : TYPE(harmonics_atom_type), OPTIONAL, POINTER :: harmonics
302 : REAL(dp), DIMENSION(:, :, :), OPTIONAL, POINTER :: Qlm_gg, gg
303 : REAL(dp), DIMENSION(:, :, :, :), OPTIONAL, POINTER :: vgg
304 : INTEGER, DIMENSION(:), OPTIONAL, POINTER :: n2oindex, o2nindex
305 : REAL(dp), DIMENSION(:, :), OPTIONAL, POINTER :: rad2l, oorad2l
306 :
307 433 : IF (ASSOCIATED(potential)) THEN
308 :
309 433 : IF (PRESENT(z)) z = potential%z
310 433 : IF (PRESENT(zeff)) zeff = potential%zeff
311 433 : IF (PRESENT(mass)) mass = potential%mass
312 433 : IF (PRESENT(elec_conf)) elec_conf => potential%elec_conf
313 433 : IF (PRESENT(nsgf)) nsgf = potential%nsgf
314 433 : IF (PRESENT(nne)) nne = potential%nne
315 433 : IF (PRESENT(npsgf)) npsgf = potential%npsgf
316 433 : IF (PRESENT(nsotot)) nsotot = potential%nsotot
317 433 : IF (PRESENT(my_gcc_h)) my_gcc_h => potential%my_gcc_h
318 433 : IF (PRESENT(my_gcc_s)) my_gcc_s => potential%my_gcc_s
319 433 : IF (PRESENT(ovlp)) ovlp => potential%ovlp
320 433 : IF (PRESENT(kin)) kin => potential%kin
321 433 : IF (PRESENT(ovlp)) ovlp => potential%ovlp
322 433 : IF (PRESENT(utrans)) utrans => potential%utrans
323 433 : IF (PRESENT(distance)) distance => potential%distance
324 433 : IF (PRESENT(harmonics)) harmonics => potential%harmonics
325 433 : IF (PRESENT(Qlm_gg)) Qlm_gg => potential%Qlm_gg
326 433 : IF (PRESENT(gg)) gg => potential%gg
327 433 : IF (PRESENT(vgg)) vgg => potential%vgg
328 433 : IF (PRESENT(n2oindex)) n2oindex => potential%n2oindex
329 433 : IF (PRESENT(o2nindex)) o2nindex => potential%o2nindex
330 433 : IF (PRESENT(rad2l)) rad2l => potential%rad2l
331 433 : IF (PRESENT(oorad2l)) oorad2l => potential%oorad2l
332 :
333 : ELSE
334 :
335 0 : CPABORT("The pointer potential is not associated.")
336 :
337 : END IF
338 :
339 433 : END SUBROUTINE get_cneo_potential
340 :
341 : ! **************************************************************************************************
342 : !> \brief ...
343 : !> \param potential ...
344 : !> \param z ...
345 : !> \param mass ...
346 : !> \param elec_conf ...
347 : !> \param nsgf ...
348 : !> \param nne ...
349 : !> \param npsgf ...
350 : !> \param nsotot ...
351 : !> \param my_gcc_h ...
352 : !> \param my_gcc_s ...
353 : !> \param ovlp ...
354 : !> \param kin ...
355 : !> \param utrans ...
356 : !> \param distance ...
357 : !> \param harmonics ...
358 : !> \param Qlm_gg ...
359 : !> \param gg ...
360 : !> \param vgg ...
361 : !> \param n2oindex ...
362 : !> \param o2nindex ...
363 : !> \param rad2l ...
364 : !> \param oorad2l ...
365 : ! **************************************************************************************************
366 50 : SUBROUTINE set_cneo_potential(potential, z, mass, elec_conf, nsgf, nne, npsgf, &
367 : nsotot, my_gcc_h, my_gcc_s, ovlp, kin, utrans, distance, &
368 : harmonics, Qlm_gg, gg, vgg, n2oindex, o2nindex, rad2l, oorad2l)
369 :
370 : TYPE(cneo_potential_type), POINTER :: potential
371 : INTEGER, INTENT(IN), OPTIONAL :: z
372 : REAL(dp), INTENT(IN), OPTIONAL :: mass
373 : INTEGER, DIMENSION(:), OPTIONAL, POINTER :: elec_conf
374 : INTEGER, INTENT(IN), OPTIONAL :: nsgf, nne, npsgf, nsotot
375 : REAL(dp), DIMENSION(:, :), OPTIONAL, POINTER :: my_gcc_h, my_gcc_s, ovlp, kin, utrans
376 : REAL(dp), DIMENSION(:, :, :), OPTIONAL, POINTER :: distance
377 : TYPE(harmonics_atom_type), OPTIONAL, POINTER :: harmonics
378 : REAL(dp), DIMENSION(:, :, :), OPTIONAL, POINTER :: Qlm_gg, gg
379 : REAL(dp), DIMENSION(:, :, :, :), OPTIONAL, POINTER :: vgg
380 : INTEGER, DIMENSION(:), OPTIONAL, POINTER :: n2oindex, o2nindex
381 : REAL(dp), DIMENSION(:, :), OPTIONAL, POINTER :: rad2l, oorad2l
382 :
383 50 : IF (ASSOCIATED(potential)) THEN
384 :
385 50 : IF (PRESENT(z)) THEN
386 8 : potential%z = z
387 8 : potential%zeff = REAL(z, dp)
388 8 : IF (ASSOCIATED(potential%elec_conf)) THEN
389 0 : CPABORT("elec_conf is already associated")
390 : END IF
391 8 : ALLOCATE (potential%elec_conf(0:3))
392 40 : potential%elec_conf(0:3) = ptable(z)%e_conv(0:3)
393 8 : CPASSERT(potential%mass == 0.0_dp)
394 8 : IF (z == 1) THEN
395 : ! Hydrogen is 1.007825, not 1.00794
396 : ! subtract the electron mass to get the proton mass
397 8 : potential%mass = 1.007825_dp - 0.000548579909_dp
398 : ELSE
399 : ! In principle, the most abundant pure isotope mass
400 : ! should be used, but no such data is available in ptable
401 0 : potential%mass = ptable(z)%amass - 0.000548579909_dp*REAL(z, dp)
402 : END IF
403 : END IF
404 50 : IF (PRESENT(mass)) THEN
405 2 : potential%mass = mass
406 : END IF
407 50 : IF (PRESENT(elec_conf)) THEN
408 0 : IF (ASSOCIATED(potential%elec_conf)) THEN
409 0 : DEALLOCATE (potential%elec_conf)
410 : END IF
411 0 : ALLOCATE (potential%elec_conf(0:SIZE(elec_conf) - 1))
412 0 : potential%elec_conf(:) = elec_conf(:)
413 : END IF
414 50 : IF (PRESENT(nsgf)) potential%nsgf = nsgf
415 50 : IF (PRESENT(nne)) potential%nne = nne
416 50 : IF (PRESENT(npsgf)) potential%npsgf = npsgf
417 50 : IF (PRESENT(nsotot)) potential%nsotot = nsotot
418 50 : IF (PRESENT(my_gcc_h)) potential%my_gcc_h => my_gcc_h
419 50 : IF (PRESENT(my_gcc_s)) potential%my_gcc_s => my_gcc_s
420 50 : IF (PRESENT(ovlp)) potential%ovlp => ovlp
421 50 : IF (PRESENT(kin)) potential%kin => kin
422 50 : IF (PRESENT(utrans)) potential%utrans => utrans
423 50 : IF (PRESENT(distance)) potential%distance => distance
424 50 : IF (PRESENT(harmonics)) potential%harmonics => harmonics
425 50 : IF (PRESENT(Qlm_gg)) potential%Qlm_gg => Qlm_gg
426 50 : IF (PRESENT(gg)) potential%gg => gg
427 50 : IF (PRESENT(vgg)) potential%vgg => vgg
428 50 : IF (PRESENT(n2oindex)) potential%n2oindex => n2oindex
429 50 : IF (PRESENT(o2nindex)) potential%o2nindex => o2nindex
430 50 : IF (PRESENT(rad2l)) potential%rad2l => rad2l
431 50 : IF (PRESENT(oorad2l)) potential%oorad2l => oorad2l
432 :
433 : ELSE
434 :
435 0 : CPABORT("The pointer potential is not associated")
436 :
437 : END IF
438 :
439 50 : END SUBROUTINE set_cneo_potential
440 :
441 : ! **************************************************************************************************
442 : !> \brief ...
443 : !> \param potential ...
444 : !> \param output_unit ...
445 : ! **************************************************************************************************
446 4 : SUBROUTINE write_cneo_potential(potential, output_unit)
447 :
448 : TYPE(cneo_potential_type), POINTER :: potential
449 : INTEGER, INTENT(IN) :: output_unit
450 :
451 : CHARACTER(LEN=20) :: string
452 :
453 4 : IF (output_unit > 0 .AND. ASSOCIATED(potential)) THEN
454 : WRITE (UNIT=output_unit, FMT="(/,T6,A,/)") &
455 4 : "CNEO Potential information"
456 : WRITE (UNIT=output_unit, FMT="(T8,A,T41,A,I4,A,F11.6)") &
457 4 : "Description: ", "Z =", potential%z, &
458 8 : ", nuclear mass =", potential%mass
459 20 : WRITE (UNIT=string, FMT="(5I4)") potential%elec_conf
460 : WRITE (UNIT=output_unit, FMT="(T8,A,T61,A20)") &
461 4 : "Electronic configuration (s p d ...):", &
462 8 : ADJUSTR(TRIM(string))
463 : END IF
464 :
465 4 : END SUBROUTINE write_cneo_potential
466 :
467 0 : END MODULE qs_cneo_types
|
