(git:15c1bfc)
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-2025 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! **************************************************************************************************
15MODULE 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
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)) &
125 DEALLOCATE (rhoz_cneo%pmat)
126 IF (ASSOCIATED(rhoz_cneo%core)) &
127 DEALLOCATE (rhoz_cneo%core)
128 IF (ASSOCIATED(rhoz_cneo%vmat)) &
129 DEALLOCATE (rhoz_cneo%vmat)
130 IF (ASSOCIATED(rhoz_cneo%fmat)) &
131 DEALLOCATE (rhoz_cneo%fmat)
132 IF (ASSOCIATED(rhoz_cneo%wfn)) &
133 DEALLOCATE (rhoz_cneo%wfn)
134 IF (ASSOCIATED(rhoz_cneo%cpc_h)) &
135 DEALLOCATE (rhoz_cneo%cpc_h)
136 IF (ASSOCIATED(rhoz_cneo%cpc_s)) &
137 DEALLOCATE (rhoz_cneo%cpc_s)
138 IF (ASSOCIATED(rhoz_cneo%rho_rad_h)) &
139 DEALLOCATE (rhoz_cneo%rho_rad_h)
140 IF (ASSOCIATED(rhoz_cneo%rho_rad_s)) &
141 DEALLOCATE (rhoz_cneo%rho_rad_s)
142 IF (ASSOCIATED(rhoz_cneo%vrho_rad_h)) &
143 DEALLOCATE (rhoz_cneo%vrho_rad_h)
144 IF (ASSOCIATED(rhoz_cneo%vrho_rad_s)) &
145 DEALLOCATE (rhoz_cneo%vrho_rad_s)
146 IF (ASSOCIATED(rhoz_cneo%ga_Vlocal_gb_h)) &
147 DEALLOCATE (rhoz_cneo%ga_Vlocal_gb_h)
148 IF (ASSOCIATED(rhoz_cneo%ga_Vlocal_gb_s)) &
149 DEALLOCATE (rhoz_cneo%ga_Vlocal_gb_s)
150 END IF
151
152 END SUBROUTINE deallocate_rhoz_cneo
153
154! **************************************************************************************************
155!> \brief ...
156!> \param rhoz_cneo_set ...
157! **************************************************************************************************
158 SUBROUTINE deallocate_rhoz_cneo_set(rhoz_cneo_set)
159
160 TYPE(rhoz_cneo_type), DIMENSION(:), POINTER :: rhoz_cneo_set
161
162 INTEGER :: iat, natom
163 TYPE(rhoz_cneo_type), POINTER :: rhoz_cneo
164
165 IF (ASSOCIATED(rhoz_cneo_set)) THEN
166 natom = SIZE(rhoz_cneo_set)
167 DO iat = 1, natom
168 rhoz_cneo => rhoz_cneo_set(iat)
169 CALL deallocate_rhoz_cneo(rhoz_cneo)
170 END DO
171 DEALLOCATE (rhoz_cneo_set)
172 END IF
173
174 END SUBROUTINE deallocate_rhoz_cneo_set
175
176! **************************************************************************************************
177!> \brief ...
178!> \param potential ...
179! **************************************************************************************************
180 SUBROUTINE allocate_cneo_potential(potential)
181
182 TYPE(cneo_potential_type), POINTER :: potential
183
184 IF (ASSOCIATED(potential)) &
185 CALL deallocate_cneo_potential(potential)
186
187 ALLOCATE (potential)
188
189 END SUBROUTINE allocate_cneo_potential
190
191! **************************************************************************************************
192!> \brief ...
193!> \param potential ...
194! **************************************************************************************************
195 SUBROUTINE deallocate_cneo_potential(potential)
196
197 TYPE(cneo_potential_type), POINTER :: potential
198
199 IF (ASSOCIATED(potential)) THEN
200 IF (ASSOCIATED(potential%elec_conf)) &
201 DEALLOCATE (potential%elec_conf)
202 IF (ASSOCIATED(potential%my_gcc_h)) &
203 DEALLOCATE (potential%my_gcc_h)
204 IF (ASSOCIATED(potential%my_gcc_s)) &
205 DEALLOCATE (potential%my_gcc_s)
206 IF (ASSOCIATED(potential%ovlp)) &
207 DEALLOCATE (potential%ovlp)
208 IF (ASSOCIATED(potential%kin)) &
209 DEALLOCATE (potential%kin)
210 IF (ASSOCIATED(potential%utrans)) &
211 DEALLOCATE (potential%utrans)
212 IF (ASSOCIATED(potential%distance)) &
213 DEALLOCATE (potential%distance)
214 IF (ASSOCIATED(potential%harmonics)) &
215 CALL deallocate_harmonics_atom(potential%harmonics)
216 IF (ASSOCIATED(potential%Qlm_gg)) &
217 DEALLOCATE (potential%Qlm_gg)
218 IF (ASSOCIATED(potential%gg)) &
219 DEALLOCATE (potential%gg)
220 IF (ASSOCIATED(potential%vgg)) &
221 DEALLOCATE (potential%vgg)
222 IF (ASSOCIATED(potential%n2oindex)) &
223 DEALLOCATE (potential%n2oindex)
224 IF (ASSOCIATED(potential%o2nindex)) &
225 DEALLOCATE (potential%o2nindex)
226 IF (ASSOCIATED(potential%rad2l)) &
227 DEALLOCATE (potential%rad2l)
228 IF (ASSOCIATED(potential%oorad2l)) &
229 DEALLOCATE (potential%oorad2l)
230 DEALLOCATE (potential)
231 END IF
232
233 END SUBROUTINE deallocate_cneo_potential
234
235! **************************************************************************************************
236!> \brief ...
237!> \param potential ...
238!> \param z ...
239!> \param zeff ...
240!> \param mass ...
241!> \param elec_conf ...
242!> \param nsgf ...
243!> \param nne ...
244!> \param npsgf ...
245!> \param nsotot ...
246!> \param my_gcc_h ...
247!> \param my_gcc_s ...
248!> \param ovlp ...
249!> \param kin ...
250!> \param utrans ...
251!> \param distance ...
252!> \param harmonics ...
253!> \param Qlm_gg ...
254!> \param gg ...
255!> \param vgg ...
256!> \param n2oindex ...
257!> \param o2nindex ...
258!> \param rad2l ...
259!> \param oorad2l ...
260! **************************************************************************************************
261 SUBROUTINE get_cneo_potential(potential, z, zeff, mass, elec_conf, nsgf, nne, npsgf, &
262 nsotot, my_gcc_h, my_gcc_s, ovlp, kin, utrans, distance, &
263 harmonics, Qlm_gg, gg, vgg, n2oindex, o2nindex, rad2l, oorad2l)
264
265 TYPE(cneo_potential_type), POINTER :: potential
266 INTEGER, INTENT(OUT), OPTIONAL :: z
267 REAL(dp), INTENT(OUT), OPTIONAL :: zeff, mass
268 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: elec_conf
269 INTEGER, INTENT(OUT), OPTIONAL :: nsgf, nne, npsgf, nsotot
270 REAL(dp), DIMENSION(:, :), OPTIONAL, POINTER :: my_gcc_h, my_gcc_s, ovlp, kin, utrans
271 REAL(dp), DIMENSION(:, :, :), OPTIONAL, POINTER :: distance
272 TYPE(harmonics_atom_type), OPTIONAL, POINTER :: harmonics
273 REAL(dp), DIMENSION(:, :, :), OPTIONAL, POINTER :: qlm_gg, gg
274 REAL(dp), DIMENSION(:, :, :, :), OPTIONAL, POINTER :: vgg
275 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: n2oindex, o2nindex
276 REAL(dp), DIMENSION(:, :), OPTIONAL, POINTER :: rad2l, oorad2l
277
278 IF (ASSOCIATED(potential)) THEN
279
280 IF (PRESENT(z)) z = potential%z
281 IF (PRESENT(zeff)) zeff = potential%zeff
282 IF (PRESENT(mass)) mass = potential%mass
283 IF (PRESENT(elec_conf)) elec_conf => potential%elec_conf
284 IF (PRESENT(nsgf)) nsgf = potential%nsgf
285 IF (PRESENT(nne)) nne = potential%nne
286 IF (PRESENT(npsgf)) npsgf = potential%npsgf
287 IF (PRESENT(nsotot)) nsotot = potential%nsotot
288 IF (PRESENT(my_gcc_h)) my_gcc_h => potential%my_gcc_h
289 IF (PRESENT(my_gcc_s)) my_gcc_s => potential%my_gcc_s
290 IF (PRESENT(ovlp)) ovlp => potential%ovlp
291 IF (PRESENT(kin)) kin => potential%kin
292 IF (PRESENT(ovlp)) ovlp => potential%ovlp
293 IF (PRESENT(utrans)) utrans => potential%utrans
294 IF (PRESENT(distance)) distance => potential%distance
295 IF (PRESENT(harmonics)) harmonics => potential%harmonics
296 IF (PRESENT(qlm_gg)) qlm_gg => potential%Qlm_gg
297 IF (PRESENT(gg)) gg => potential%gg
298 IF (PRESENT(vgg)) vgg => potential%vgg
299 IF (PRESENT(n2oindex)) n2oindex => potential%n2oindex
300 IF (PRESENT(o2nindex)) o2nindex => potential%o2nindex
301 IF (PRESENT(rad2l)) rad2l => potential%rad2l
302 IF (PRESENT(oorad2l)) oorad2l => potential%oorad2l
303
304 ELSE
305
306 cpabort("The pointer potential is not associated.")
307
308 END IF
309
310 END SUBROUTINE get_cneo_potential
311
312! **************************************************************************************************
313!> \brief ...
314!> \param potential ...
315!> \param z ...
316!> \param mass ...
317!> \param elec_conf ...
318!> \param nsgf ...
319!> \param nne ...
320!> \param npsgf ...
321!> \param nsotot ...
322!> \param my_gcc_h ...
323!> \param my_gcc_s ...
324!> \param ovlp ...
325!> \param kin ...
326!> \param utrans ...
327!> \param distance ...
328!> \param harmonics ...
329!> \param Qlm_gg ...
330!> \param gg ...
331!> \param vgg ...
332!> \param n2oindex ...
333!> \param o2nindex ...
334!> \param rad2l ...
335!> \param oorad2l ...
336! **************************************************************************************************
337 SUBROUTINE set_cneo_potential(potential, z, mass, elec_conf, nsgf, nne, npsgf, &
338 nsotot, my_gcc_h, my_gcc_s, ovlp, kin, utrans, distance, &
339 harmonics, Qlm_gg, gg, vgg, n2oindex, o2nindex, rad2l, oorad2l)
340
341 TYPE(cneo_potential_type), POINTER :: potential
342 INTEGER, INTENT(IN), OPTIONAL :: z
343 REAL(dp), INTENT(IN), OPTIONAL :: mass
344 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: elec_conf
345 INTEGER, INTENT(IN), OPTIONAL :: nsgf, nne, npsgf, nsotot
346 REAL(dp), DIMENSION(:, :), OPTIONAL, POINTER :: my_gcc_h, my_gcc_s, ovlp, kin, utrans
347 REAL(dp), DIMENSION(:, :, :), OPTIONAL, POINTER :: distance
348 TYPE(harmonics_atom_type), OPTIONAL, POINTER :: harmonics
349 REAL(dp), DIMENSION(:, :, :), OPTIONAL, POINTER :: qlm_gg, gg
350 REAL(dp), DIMENSION(:, :, :, :), OPTIONAL, POINTER :: vgg
351 INTEGER, DIMENSION(:), OPTIONAL, POINTER :: n2oindex, o2nindex
352 REAL(dp), DIMENSION(:, :), OPTIONAL, POINTER :: rad2l, oorad2l
353
354 IF (ASSOCIATED(potential)) THEN
355
356 IF (PRESENT(z)) THEN
357 potential%z = z
358 potential%zeff = real(z, dp)
359 IF (ASSOCIATED(potential%elec_conf)) &
360 cpabort("elec_conf is already associated")
361 ALLOCATE (potential%elec_conf(0:3))
362 potential%elec_conf(0:3) = ptable(z)%e_conv(0:3)
363 cpassert(potential%mass == 0.0_dp)
364 IF (z == 1) THEN
365 ! Hydrogen is 1.007825, not 1.00794
366 ! subtract the electron mass to get the proton mass
367 potential%mass = 1.007825_dp - 0.000548579909_dp
368 ELSE
369 ! In principle, the most abundant pure isotope mass
370 ! should be used, but no such data is available in ptable
371 potential%mass = ptable(z)%amass - 0.000548579909_dp*real(z, dp)
372 END IF
373 END IF
374 IF (PRESENT(mass)) THEN
375 potential%mass = mass
376 END IF
377 IF (PRESENT(elec_conf)) THEN
378 IF (ASSOCIATED(potential%elec_conf)) THEN
379 DEALLOCATE (potential%elec_conf)
380 END IF
381 ALLOCATE (potential%elec_conf(0:SIZE(elec_conf) - 1))
382 potential%elec_conf(:) = elec_conf(:)
383 END IF
384 IF (PRESENT(nsgf)) potential%nsgf = nsgf
385 IF (PRESENT(nne)) potential%nne = nne
386 IF (PRESENT(npsgf)) potential%npsgf = npsgf
387 IF (PRESENT(nsotot)) potential%nsotot = nsotot
388 IF (PRESENT(my_gcc_h)) potential%my_gcc_h => my_gcc_h
389 IF (PRESENT(my_gcc_s)) potential%my_gcc_s => my_gcc_s
390 IF (PRESENT(ovlp)) potential%ovlp => ovlp
391 IF (PRESENT(kin)) potential%kin => kin
392 IF (PRESENT(utrans)) potential%utrans => utrans
393 IF (PRESENT(distance)) potential%distance => distance
394 IF (PRESENT(harmonics)) potential%harmonics => harmonics
395 IF (PRESENT(qlm_gg)) potential%Qlm_gg => qlm_gg
396 IF (PRESENT(gg)) potential%gg => gg
397 IF (PRESENT(vgg)) potential%vgg => vgg
398 IF (PRESENT(n2oindex)) potential%n2oindex => n2oindex
399 IF (PRESENT(o2nindex)) potential%o2nindex => o2nindex
400 IF (PRESENT(rad2l)) potential%rad2l => rad2l
401 IF (PRESENT(oorad2l)) potential%oorad2l => oorad2l
402
403 ELSE
404
405 cpabort("The pointer potential is not associated")
406
407 END IF
408
409 END SUBROUTINE set_cneo_potential
410
411! **************************************************************************************************
412!> \brief ...
413!> \param potential ...
414!> \param output_unit ...
415! **************************************************************************************************
416 SUBROUTINE write_cneo_potential(potential, output_unit)
417
418 TYPE(cneo_potential_type), POINTER :: potential
419 INTEGER, INTENT(IN) :: output_unit
420
421 CHARACTER(LEN=20) :: string
422
423 IF (output_unit > 0 .AND. ASSOCIATED(potential)) THEN
424 WRITE (unit=output_unit, fmt="(/,T6,A,/)") &
425 "CNEO Potential information"
426 WRITE (unit=output_unit, fmt="(T8,A,T41,A,I4,A,F11.6)") &
427 "Description: ", "Z =", potential%z, &
428 ", nuclear mass =", potential%mass
429 WRITE (unit=string, fmt="(5I4)") potential%elec_conf
430 WRITE (unit=output_unit, fmt="(T8,A,T61,A20)") &
431 "Electronic configuration (s p d ...):", &
432 adjustr(trim(string))
433 END IF
434
435 END SUBROUTINE write_cneo_potential
436
437END MODULE qs_cneo_types
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
periodic_table
Periodic Table related data definitions.
Definition periodic_table.F:28
periodic_table::ptable
type(atom), dimension(0:nelem), public ptable
Definition periodic_table.F:65
qs_cneo_types
Types used by CNEO-DFT (see J. Chem. Theory Comput. 2025, 21, 16, 7865–7877)
Definition qs_cneo_types.F:15
qs_cneo_types::allocate_cneo_potential
subroutine, public allocate_cneo_potential(potential)
...
Definition qs_cneo_types.F:181
qs_cneo_types::write_cneo_potential
subroutine, public write_cneo_potential(potential, output_unit)
...
Definition qs_cneo_types.F:417
qs_cneo_types::deallocate_rhoz_cneo_set
subroutine, public deallocate_rhoz_cneo_set(rhoz_cneo_set)
...
Definition qs_cneo_types.F:159
qs_cneo_types::deallocate_cneo_potential
subroutine, public deallocate_cneo_potential(potential)
...
Definition qs_cneo_types.F:196
qs_cneo_types::set_cneo_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)
...
Definition qs_cneo_types.F:340
qs_cneo_types::get_cneo_potential
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)
...
Definition qs_cneo_types.F:264
qs_cneo_types::allocate_rhoz_cneo_set
subroutine, public allocate_rhoz_cneo_set(rhoz_cneo_set, natom)
...
Definition qs_cneo_types.F:103
qs_harmonics_atom
Definition qs_harmonics_atom.F:8
qs_harmonics_atom::deallocate_harmonics_atom
subroutine, public deallocate_harmonics_atom(harmonics)
Deallocate the spherical harmonics set for the atom grid.
Definition qs_harmonics_atom.F:97
qs_cneo_types::cneo_potential_type
Definition qs_cneo_types.F:53
qs_cneo_types::rhoz_cneo_type
Definition qs_cneo_types.F:29
qs_harmonics_atom::harmonics_atom_type
Definition qs_harmonics_atom.F:32