(git:71c3ab0)
Loading...
Searching...
No Matches
qs_cneo_types.F
Go to the documentation of this file.
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! **************************************************************************************************
16 USE kinds, ONLY: dp
17 USE periodic_table, ONLY: ptable
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
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
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
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
89
90! Public Subroutine
91
94
95CONTAINS
96
97! **************************************************************************************************
98!> \brief ...
99!> \param rhoz_cneo_set ...
100!> \param natom ...
101! **************************************************************************************************
102 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 IF (ASSOCIATED(rhoz_cneo_set)) THEN
108 CALL deallocate_rhoz_cneo_set(rhoz_cneo_set)
109 END IF
110
111 ALLOCATE (rhoz_cneo_set(natom))
112
113 END SUBROUTINE allocate_rhoz_cneo_set
114
115! **************************************************************************************************
116!> \brief ...
117!> \param rhoz_cneo ...
118! **************************************************************************************************
119 SUBROUTINE deallocate_rhoz_cneo(rhoz_cneo)
120
121 TYPE(rhoz_cneo_type), POINTER :: rhoz_cneo
122
123 IF (ASSOCIATED(rhoz_cneo)) THEN
124 IF (ASSOCIATED(rhoz_cneo%pmat)) THEN
125 DEALLOCATE (rhoz_cneo%pmat)
126 END IF
127 IF (ASSOCIATED(rhoz_cneo%core)) THEN
128 DEALLOCATE (rhoz_cneo%core)
129 END IF
130 IF (ASSOCIATED(rhoz_cneo%vmat)) THEN
131 DEALLOCATE (rhoz_cneo%vmat)
132 END IF
133 IF (ASSOCIATED(rhoz_cneo%fmat)) THEN
134 DEALLOCATE (rhoz_cneo%fmat)
135 END IF
136 IF (ASSOCIATED(rhoz_cneo%wfn)) THEN
137 DEALLOCATE (rhoz_cneo%wfn)
138 END IF
139 IF (ASSOCIATED(rhoz_cneo%cpc_h)) THEN
140 DEALLOCATE (rhoz_cneo%cpc_h)
141 END IF
142 IF (ASSOCIATED(rhoz_cneo%cpc_s)) THEN
143 DEALLOCATE (rhoz_cneo%cpc_s)
144 END IF
145 IF (ASSOCIATED(rhoz_cneo%rho_rad_h)) THEN
146 DEALLOCATE (rhoz_cneo%rho_rad_h)
147 END IF
148 IF (ASSOCIATED(rhoz_cneo%rho_rad_s)) THEN
149 DEALLOCATE (rhoz_cneo%rho_rad_s)
150 END IF
151 IF (ASSOCIATED(rhoz_cneo%vrho_rad_h)) THEN
152 DEALLOCATE (rhoz_cneo%vrho_rad_h)
153 END IF
154 IF (ASSOCIATED(rhoz_cneo%vrho_rad_s)) THEN
155 DEALLOCATE (rhoz_cneo%vrho_rad_s)
156 END IF
157 IF (ASSOCIATED(rhoz_cneo%ga_Vlocal_gb_h)) THEN
158 DEALLOCATE (rhoz_cneo%ga_Vlocal_gb_h)
159 END IF
160 IF (ASSOCIATED(rhoz_cneo%ga_Vlocal_gb_s)) THEN
161 DEALLOCATE (rhoz_cneo%ga_Vlocal_gb_s)
162 END IF
163 END IF
164
165 END SUBROUTINE deallocate_rhoz_cneo
166
167! **************************************************************************************************
168!> \brief ...
169!> \param rhoz_cneo_set ...
170! **************************************************************************************************
171 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 IF (ASSOCIATED(rhoz_cneo_set)) THEN
179 natom = SIZE(rhoz_cneo_set)
180 DO iat = 1, natom
181 rhoz_cneo => rhoz_cneo_set(iat)
182 CALL deallocate_rhoz_cneo(rhoz_cneo)
183 END DO
184 DEALLOCATE (rhoz_cneo_set)
185 END IF
186
187 END SUBROUTINE deallocate_rhoz_cneo_set
188
189! **************************************************************************************************
190!> \brief ...
191!> \param potential ...
192! **************************************************************************************************
193 SUBROUTINE allocate_cneo_potential(potential)
194
195 TYPE(cneo_potential_type), POINTER :: potential
196
197 IF (ASSOCIATED(potential)) THEN
198 CALL deallocate_cneo_potential(potential)
199 END IF
200
201 ALLOCATE (potential)
202
203 END SUBROUTINE allocate_cneo_potential
204
205! **************************************************************************************************
206!> \brief ...
207!> \param potential ...
208! **************************************************************************************************
209 SUBROUTINE deallocate_cneo_potential(potential)
210
211 TYPE(cneo_potential_type), POINTER :: potential
212
213 IF (ASSOCIATED(potential)) THEN
214 IF (ASSOCIATED(potential%elec_conf)) THEN
215 DEALLOCATE (potential%elec_conf)
216 END IF
217 IF (ASSOCIATED(potential%my_gcc_h)) THEN
218 DEALLOCATE (potential%my_gcc_h)
219 END IF
220 IF (ASSOCIATED(potential%my_gcc_s)) THEN
221 DEALLOCATE (potential%my_gcc_s)
222 END IF
223 IF (ASSOCIATED(potential%ovlp)) THEN
224 DEALLOCATE (potential%ovlp)
225 END IF
226 IF (ASSOCIATED(potential%kin)) THEN
227 DEALLOCATE (potential%kin)
228 END IF
229 IF (ASSOCIATED(potential%utrans)) THEN
230 DEALLOCATE (potential%utrans)
231 END IF
232 IF (ASSOCIATED(potential%distance)) THEN
233 DEALLOCATE (potential%distance)
234 END IF
235 IF (ASSOCIATED(potential%harmonics)) THEN
236 CALL deallocate_harmonics_atom(potential%harmonics)
237 END IF
238 IF (ASSOCIATED(potential%Qlm_gg)) THEN
239 DEALLOCATE (potential%Qlm_gg)
240 END IF
241 IF (ASSOCIATED(potential%gg)) THEN
242 DEALLOCATE (potential%gg)
243 END IF
244 IF (ASSOCIATED(potential%vgg)) THEN
245 DEALLOCATE (potential%vgg)
246 END IF
247 IF (ASSOCIATED(potential%n2oindex)) THEN
248 DEALLOCATE (potential%n2oindex)
249 END IF
250 IF (ASSOCIATED(potential%o2nindex)) THEN
251 DEALLOCATE (potential%o2nindex)
252 END IF
253 IF (ASSOCIATED(potential%rad2l)) THEN
254 DEALLOCATE (potential%rad2l)
255 END IF
256 IF (ASSOCIATED(potential%oorad2l)) THEN
257 DEALLOCATE (potential%oorad2l)
258 END IF
259 DEALLOCATE (potential)
260 END IF
261
262 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 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 IF (ASSOCIATED(potential)) THEN
308
309 IF (PRESENT(z)) z = potential%z
310 IF (PRESENT(zeff)) zeff = potential%zeff
311 IF (PRESENT(mass)) mass = potential%mass
312 IF (PRESENT(elec_conf)) elec_conf => potential%elec_conf
313 IF (PRESENT(nsgf)) nsgf = potential%nsgf
314 IF (PRESENT(nne)) nne = potential%nne
315 IF (PRESENT(npsgf)) npsgf = potential%npsgf
316 IF (PRESENT(nsotot)) nsotot = potential%nsotot
317 IF (PRESENT(my_gcc_h)) my_gcc_h => potential%my_gcc_h
318 IF (PRESENT(my_gcc_s)) my_gcc_s => potential%my_gcc_s
319 IF (PRESENT(ovlp)) ovlp => potential%ovlp
320 IF (PRESENT(kin)) kin => potential%kin
321 IF (PRESENT(ovlp)) ovlp => potential%ovlp
322 IF (PRESENT(utrans)) utrans => potential%utrans
323 IF (PRESENT(distance)) distance => potential%distance
324 IF (PRESENT(harmonics)) harmonics => potential%harmonics
325 IF (PRESENT(qlm_gg)) qlm_gg => potential%Qlm_gg
326 IF (PRESENT(gg)) gg => potential%gg
327 IF (PRESENT(vgg)) vgg => potential%vgg
328 IF (PRESENT(n2oindex)) n2oindex => potential%n2oindex
329 IF (PRESENT(o2nindex)) o2nindex => potential%o2nindex
330 IF (PRESENT(rad2l)) rad2l => potential%rad2l
331 IF (PRESENT(oorad2l)) oorad2l => potential%oorad2l
332
333 ELSE
334
335 cpabort("The pointer potential is not associated.")
336
337 END IF
338
339 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 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 IF (ASSOCIATED(potential)) THEN
384
385 IF (PRESENT(z)) THEN
386 potential%z = z
387 potential%zeff = real(z, dp)
388 IF (ASSOCIATED(potential%elec_conf)) THEN
389 cpabort("elec_conf is already associated")
390 END IF
391 ALLOCATE (potential%elec_conf(0:3))
392 potential%elec_conf(0:3) = ptable(z)%e_conv(0:3)
393 cpassert(potential%mass == 0.0_dp)
394 IF (z == 1) THEN
395 ! Hydrogen is 1.007825, not 1.00794
396 ! subtract the electron mass to get the proton mass
397 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 potential%mass = ptable(z)%amass - 0.000548579909_dp*real(z, dp)
402 END IF
403 END IF
404 IF (PRESENT(mass)) THEN
405 potential%mass = mass
406 END IF
407 IF (PRESENT(elec_conf)) THEN
408 IF (ASSOCIATED(potential%elec_conf)) THEN
409 DEALLOCATE (potential%elec_conf)
410 END IF
411 ALLOCATE (potential%elec_conf(0:SIZE(elec_conf) - 1))
412 potential%elec_conf(:) = elec_conf(:)
413 END IF
414 IF (PRESENT(nsgf)) potential%nsgf = nsgf
415 IF (PRESENT(nne)) potential%nne = nne
416 IF (PRESENT(npsgf)) potential%npsgf = npsgf
417 IF (PRESENT(nsotot)) potential%nsotot = nsotot
418 IF (PRESENT(my_gcc_h)) potential%my_gcc_h => my_gcc_h
419 IF (PRESENT(my_gcc_s)) potential%my_gcc_s => my_gcc_s
420 IF (PRESENT(ovlp)) potential%ovlp => ovlp
421 IF (PRESENT(kin)) potential%kin => kin
422 IF (PRESENT(utrans)) potential%utrans => utrans
423 IF (PRESENT(distance)) potential%distance => distance
424 IF (PRESENT(harmonics)) potential%harmonics => harmonics
425 IF (PRESENT(qlm_gg)) potential%Qlm_gg => qlm_gg
426 IF (PRESENT(gg)) potential%gg => gg
427 IF (PRESENT(vgg)) potential%vgg => vgg
428 IF (PRESENT(n2oindex)) potential%n2oindex => n2oindex
429 IF (PRESENT(o2nindex)) potential%o2nindex => o2nindex
430 IF (PRESENT(rad2l)) potential%rad2l => rad2l
431 IF (PRESENT(oorad2l)) potential%oorad2l => oorad2l
432
433 ELSE
434
435 cpabort("The pointer potential is not associated")
436
437 END IF
438
439 END SUBROUTINE set_cneo_potential
440
441! **************************************************************************************************
442!> \brief ...
443!> \param potential ...
444!> \param output_unit ...
445! **************************************************************************************************
446 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 IF (output_unit > 0 .AND. ASSOCIATED(potential)) THEN
454 WRITE (unit=output_unit, fmt="(/,T6,A,/)") &
455 "CNEO Potential information"
456 WRITE (unit=output_unit, fmt="(T8,A,T41,A,I4,A,F11.6)") &
457 "Description: ", "Z =", potential%z, &
458 ", nuclear mass =", potential%mass
459 WRITE (unit=string, fmt="(5I4)") potential%elec_conf
460 WRITE (unit=output_unit, fmt="(T8,A,T61,A20)") &
461 "Electronic configuration (s p d ...):", &
462 adjustr(trim(string))
463 END IF
464
465 END SUBROUTINE write_cneo_potential
466
467END MODULE qs_cneo_types
Defines the basic variable types.
Definition kinds.F:23
integer, parameter, public dp
Definition kinds.F:34
Periodic Table related data definitions.
type(atom), dimension(0:nelem), public ptable
Types used by CNEO-DFT (see J. Chem. Theory Comput. 2025, 21, 16, 7865–7877)
subroutine, public allocate_cneo_potential(potential)
...
subroutine, public write_cneo_potential(potential, output_unit)
...
subroutine, public deallocate_rhoz_cneo_set(rhoz_cneo_set)
...
subroutine, public deallocate_cneo_potential(potential)
...
subroutine, public 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)
...
subroutine, public 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)
...
subroutine, public allocate_rhoz_cneo_set(rhoz_cneo_set, natom)
...
subroutine, public deallocate_harmonics_atom(harmonics)
Deallocate the spherical harmonics set for the atom grid.