Line data Source code
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 : ! **************************************************************************************************
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)) &
125 14 : DEALLOCATE (rhoz_cneo%pmat)
126 20 : IF (ASSOCIATED(rhoz_cneo%core)) &
127 14 : DEALLOCATE (rhoz_cneo%core)
128 20 : IF (ASSOCIATED(rhoz_cneo%vmat)) &
129 14 : DEALLOCATE (rhoz_cneo%vmat)
130 20 : IF (ASSOCIATED(rhoz_cneo%fmat)) &
131 7 : DEALLOCATE (rhoz_cneo%fmat)
132 20 : IF (ASSOCIATED(rhoz_cneo%wfn)) &
133 7 : DEALLOCATE (rhoz_cneo%wfn)
134 20 : IF (ASSOCIATED(rhoz_cneo%cpc_h)) &
135 14 : DEALLOCATE (rhoz_cneo%cpc_h)
136 20 : IF (ASSOCIATED(rhoz_cneo%cpc_s)) &
137 14 : DEALLOCATE (rhoz_cneo%cpc_s)
138 20 : IF (ASSOCIATED(rhoz_cneo%rho_rad_h)) &
139 7 : DEALLOCATE (rhoz_cneo%rho_rad_h)
140 20 : IF (ASSOCIATED(rhoz_cneo%rho_rad_s)) &
141 7 : DEALLOCATE (rhoz_cneo%rho_rad_s)
142 20 : IF (ASSOCIATED(rhoz_cneo%vrho_rad_h)) &
143 7 : DEALLOCATE (rhoz_cneo%vrho_rad_h)
144 20 : IF (ASSOCIATED(rhoz_cneo%vrho_rad_s)) &
145 7 : DEALLOCATE (rhoz_cneo%vrho_rad_s)
146 20 : IF (ASSOCIATED(rhoz_cneo%ga_Vlocal_gb_h)) &
147 7 : DEALLOCATE (rhoz_cneo%ga_Vlocal_gb_h)
148 20 : IF (ASSOCIATED(rhoz_cneo%ga_Vlocal_gb_s)) &
149 7 : DEALLOCATE (rhoz_cneo%ga_Vlocal_gb_s)
150 : END IF
151 :
152 20 : END SUBROUTINE deallocate_rhoz_cneo
153 :
154 : ! **************************************************************************************************
155 : !> \brief ...
156 : !> \param rhoz_cneo_set ...
157 : ! **************************************************************************************************
158 8 : 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 8 : IF (ASSOCIATED(rhoz_cneo_set)) THEN
166 8 : natom = SIZE(rhoz_cneo_set)
167 28 : DO iat = 1, natom
168 20 : rhoz_cneo => rhoz_cneo_set(iat)
169 28 : CALL deallocate_rhoz_cneo(rhoz_cneo)
170 : END DO
171 8 : DEALLOCATE (rhoz_cneo_set)
172 : END IF
173 :
174 8 : END SUBROUTINE deallocate_rhoz_cneo_set
175 :
176 : ! **************************************************************************************************
177 : !> \brief ...
178 : !> \param potential ...
179 : ! **************************************************************************************************
180 8 : SUBROUTINE allocate_cneo_potential(potential)
181 :
182 : TYPE(cneo_potential_type), POINTER :: potential
183 :
184 8 : IF (ASSOCIATED(potential)) &
185 0 : CALL deallocate_cneo_potential(potential)
186 :
187 8 : ALLOCATE (potential)
188 :
189 8 : END SUBROUTINE allocate_cneo_potential
190 :
191 : ! **************************************************************************************************
192 : !> \brief ...
193 : !> \param potential ...
194 : ! **************************************************************************************************
195 8 : SUBROUTINE deallocate_cneo_potential(potential)
196 :
197 : TYPE(cneo_potential_type), POINTER :: potential
198 :
199 8 : IF (ASSOCIATED(potential)) THEN
200 8 : IF (ASSOCIATED(potential%elec_conf)) &
201 8 : DEALLOCATE (potential%elec_conf)
202 8 : IF (ASSOCIATED(potential%my_gcc_h)) &
203 8 : DEALLOCATE (potential%my_gcc_h)
204 8 : IF (ASSOCIATED(potential%my_gcc_s)) &
205 8 : DEALLOCATE (potential%my_gcc_s)
206 8 : IF (ASSOCIATED(potential%ovlp)) &
207 8 : DEALLOCATE (potential%ovlp)
208 8 : IF (ASSOCIATED(potential%kin)) &
209 8 : DEALLOCATE (potential%kin)
210 8 : IF (ASSOCIATED(potential%utrans)) &
211 8 : DEALLOCATE (potential%utrans)
212 8 : IF (ASSOCIATED(potential%distance)) &
213 8 : DEALLOCATE (potential%distance)
214 8 : IF (ASSOCIATED(potential%harmonics)) &
215 8 : CALL deallocate_harmonics_atom(potential%harmonics)
216 8 : IF (ASSOCIATED(potential%Qlm_gg)) &
217 8 : DEALLOCATE (potential%Qlm_gg)
218 8 : IF (ASSOCIATED(potential%gg)) &
219 8 : DEALLOCATE (potential%gg)
220 8 : IF (ASSOCIATED(potential%vgg)) &
221 8 : DEALLOCATE (potential%vgg)
222 8 : IF (ASSOCIATED(potential%n2oindex)) &
223 8 : DEALLOCATE (potential%n2oindex)
224 8 : IF (ASSOCIATED(potential%o2nindex)) &
225 8 : DEALLOCATE (potential%o2nindex)
226 8 : IF (ASSOCIATED(potential%rad2l)) &
227 8 : DEALLOCATE (potential%rad2l)
228 8 : IF (ASSOCIATED(potential%oorad2l)) &
229 8 : DEALLOCATE (potential%oorad2l)
230 8 : DEALLOCATE (potential)
231 : END IF
232 :
233 8 : 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 385 : 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 385 : IF (ASSOCIATED(potential)) THEN
279 :
280 385 : IF (PRESENT(z)) z = potential%z
281 385 : IF (PRESENT(zeff)) zeff = potential%zeff
282 385 : IF (PRESENT(mass)) mass = potential%mass
283 385 : IF (PRESENT(elec_conf)) elec_conf => potential%elec_conf
284 385 : IF (PRESENT(nsgf)) nsgf = potential%nsgf
285 385 : IF (PRESENT(nne)) nne = potential%nne
286 385 : IF (PRESENT(npsgf)) npsgf = potential%npsgf
287 385 : IF (PRESENT(nsotot)) nsotot = potential%nsotot
288 385 : IF (PRESENT(my_gcc_h)) my_gcc_h => potential%my_gcc_h
289 385 : IF (PRESENT(my_gcc_s)) my_gcc_s => potential%my_gcc_s
290 385 : IF (PRESENT(ovlp)) ovlp => potential%ovlp
291 385 : IF (PRESENT(kin)) kin => potential%kin
292 385 : IF (PRESENT(ovlp)) ovlp => potential%ovlp
293 385 : IF (PRESENT(utrans)) utrans => potential%utrans
294 385 : IF (PRESENT(distance)) distance => potential%distance
295 385 : IF (PRESENT(harmonics)) harmonics => potential%harmonics
296 385 : IF (PRESENT(Qlm_gg)) Qlm_gg => potential%Qlm_gg
297 385 : IF (PRESENT(gg)) gg => potential%gg
298 385 : IF (PRESENT(vgg)) vgg => potential%vgg
299 385 : IF (PRESENT(n2oindex)) n2oindex => potential%n2oindex
300 385 : IF (PRESENT(o2nindex)) o2nindex => potential%o2nindex
301 385 : IF (PRESENT(rad2l)) rad2l => potential%rad2l
302 385 : IF (PRESENT(oorad2l)) oorad2l => potential%oorad2l
303 :
304 : ELSE
305 :
306 0 : CPABORT("The pointer potential is not associated.")
307 :
308 : END IF
309 :
310 385 : 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 50 : 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 50 : IF (ASSOCIATED(potential)) THEN
355 :
356 50 : IF (PRESENT(z)) THEN
357 8 : potential%z = z
358 8 : potential%zeff = REAL(z, dp)
359 8 : IF (ASSOCIATED(potential%elec_conf)) &
360 0 : CPABORT("elec_conf is already associated")
361 8 : ALLOCATE (potential%elec_conf(0:3))
362 40 : potential%elec_conf(0:3) = ptable(z)%e_conv(0:3)
363 8 : CPASSERT(potential%mass == 0.0_dp)
364 8 : IF (z == 1) THEN
365 : ! Hydrogen is 1.007825, not 1.00794
366 : ! subtract the electron mass to get the proton mass
367 8 : 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 0 : potential%mass = ptable(z)%amass - 0.000548579909_dp*REAL(z, dp)
372 : END IF
373 : END IF
374 50 : IF (PRESENT(mass)) THEN
375 2 : potential%mass = mass
376 : END IF
377 50 : IF (PRESENT(elec_conf)) THEN
378 0 : IF (ASSOCIATED(potential%elec_conf)) THEN
379 0 : DEALLOCATE (potential%elec_conf)
380 : END IF
381 0 : ALLOCATE (potential%elec_conf(0:SIZE(elec_conf) - 1))
382 0 : potential%elec_conf(:) = elec_conf(:)
383 : END IF
384 50 : IF (PRESENT(nsgf)) potential%nsgf = nsgf
385 50 : IF (PRESENT(nne)) potential%nne = nne
386 50 : IF (PRESENT(npsgf)) potential%npsgf = npsgf
387 50 : IF (PRESENT(nsotot)) potential%nsotot = nsotot
388 50 : IF (PRESENT(my_gcc_h)) potential%my_gcc_h => my_gcc_h
389 50 : IF (PRESENT(my_gcc_s)) potential%my_gcc_s => my_gcc_s
390 50 : IF (PRESENT(ovlp)) potential%ovlp => ovlp
391 50 : IF (PRESENT(kin)) potential%kin => kin
392 50 : IF (PRESENT(utrans)) potential%utrans => utrans
393 50 : IF (PRESENT(distance)) potential%distance => distance
394 50 : IF (PRESENT(harmonics)) potential%harmonics => harmonics
395 50 : IF (PRESENT(Qlm_gg)) potential%Qlm_gg => Qlm_gg
396 50 : IF (PRESENT(gg)) potential%gg => gg
397 50 : IF (PRESENT(vgg)) potential%vgg => vgg
398 50 : IF (PRESENT(n2oindex)) potential%n2oindex => n2oindex
399 50 : IF (PRESENT(o2nindex)) potential%o2nindex => o2nindex
400 50 : IF (PRESENT(rad2l)) potential%rad2l => rad2l
401 50 : IF (PRESENT(oorad2l)) potential%oorad2l => oorad2l
402 :
403 : ELSE
404 :
405 0 : CPABORT("The pointer potential is not associated")
406 :
407 : END IF
408 :
409 50 : END SUBROUTINE set_cneo_potential
410 :
411 : ! **************************************************************************************************
412 : !> \brief ...
413 : !> \param potential ...
414 : !> \param output_unit ...
415 : ! **************************************************************************************************
416 4 : 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 4 : IF (output_unit > 0 .AND. ASSOCIATED(potential)) THEN
424 : WRITE (UNIT=output_unit, FMT="(/,T6,A,/)") &
425 4 : "CNEO Potential information"
426 : WRITE (UNIT=output_unit, FMT="(T8,A,T41,A,I4,A,F11.6)") &
427 4 : "Description: ", "Z =", potential%z, &
428 8 : ", nuclear mass =", potential%mass
429 20 : WRITE (UNIT=string, FMT="(5I4)") potential%elec_conf
430 : WRITE (UNIT=output_unit, FMT="(T8,A,T61,A20)") &
431 4 : "Electronic configuration (s p d ...):", &
432 8 : ADJUSTR(TRIM(string))
433 : END IF
434 :
435 4 : END SUBROUTINE write_cneo_potential
436 :
437 0 : END MODULE qs_cneo_types
|
