Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Calculation of the energies concerning the core charge distribution
10 : !> \par History
11 : !> - Full refactoring of calculate_ecore and calculate_ecore_overlap (jhu)
12 : !> \author Matthias Krack (27.04.2001)
13 : ! **************************************************************************************************
14 : MODULE qs_core_energies
15 : USE atomic_kind_types, ONLY: atomic_kind_type,&
16 : get_atomic_kind,&
17 : get_atomic_kind_set
18 : USE atprop_types, ONLY: atprop_array_init,&
19 : atprop_type
20 : USE cell_types, ONLY: cell_type,&
21 : pbc
22 : USE cp_dbcsr_api, ONLY: dbcsr_dot,&
23 : dbcsr_p_type,&
24 : dbcsr_type
25 : USE distribution_1d_types, ONLY: distribution_1d_type
26 : USE kinds, ONLY: dp
27 : USE mathconstants, ONLY: oorootpi,&
28 : twopi
29 : USE message_passing, ONLY: mp_comm_type,&
30 : mp_para_env_type
31 : USE particle_types, ONLY: particle_type
32 : USE qs_energy_types, ONLY: qs_energy_type
33 : USE qs_environment_types, ONLY: get_qs_env,&
34 : qs_environment_type
35 : USE qs_force_types, ONLY: qs_force_type
36 : USE qs_kind_types, ONLY: get_qs_kind,&
37 : qs_kind_type
38 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
39 : neighbor_list_iterate,&
40 : neighbor_list_iterator_create,&
41 : neighbor_list_iterator_p_type,&
42 : neighbor_list_iterator_release,&
43 : neighbor_list_set_p_type
44 : USE virial_methods, ONLY: virial_pair_force
45 : USE virial_types, ONLY: virial_type
46 : #include "./base/base_uses.f90"
47 :
48 : IMPLICIT NONE
49 :
50 : PRIVATE
51 :
52 : ! *** Global parameters ***
53 :
54 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_core_energies'
55 :
56 : PUBLIC :: calculate_ptrace, &
57 : calculate_ecore_overlap, &
58 : calculate_ecore_self, calculate_ecore_alpha
59 :
60 : INTERFACE calculate_ptrace
61 : MODULE PROCEDURE calculate_ptrace_1, calculate_ptrace_gamma, calculate_ptrace_kp
62 : END INTERFACE
63 :
64 : ! **************************************************************************************************
65 :
66 : CONTAINS
67 :
68 : ! **************************************************************************************************
69 : !> \brief Calculate the trace of a operator matrix with the density matrix.
70 : !> Sum over all spin components (in P, no spin in H)
71 : !> \param hmat ...
72 : !> \param pmat ...
73 : !> \param ecore ...
74 : !> \param nspin ...
75 : !> \date 29.07.2014
76 : !> \par History
77 : !> - none
78 : !> \author JGH
79 : !> \version 1.0
80 : ! **************************************************************************************************
81 132 : SUBROUTINE calculate_ptrace_gamma(hmat, pmat, ecore, nspin)
82 :
83 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: hmat, pmat
84 : REAL(KIND=dp), INTENT(OUT) :: ecore
85 : INTEGER, INTENT(IN) :: nspin
86 :
87 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_ptrace_gamma'
88 :
89 : INTEGER :: handle, ispin
90 : REAL(KIND=dp) :: etr
91 :
92 132 : CALL timeset(routineN, handle)
93 :
94 132 : ecore = 0.0_dp
95 264 : DO ispin = 1, nspin
96 132 : etr = 0.0_dp
97 132 : CALL dbcsr_dot(hmat(1)%matrix, pmat(ispin)%matrix, etr)
98 264 : ecore = ecore + etr
99 : END DO
100 :
101 132 : CALL timestop(handle)
102 :
103 132 : END SUBROUTINE calculate_ptrace_gamma
104 :
105 : ! **************************************************************************************************
106 : !> \brief Calculate the trace of a operator matrix with the density matrix.
107 : !> Sum over all spin components (in P, no spin in H) and the real space
108 : !> coordinates
109 : !> \param hmat H matrix
110 : !> \param pmat P matrices
111 : !> \param ecore Tr(HP) output
112 : !> \param nspin Number of P matrices
113 : !> \date 29.07.2014
114 : !> \author JGH
115 : !> \version 1.0
116 : ! **************************************************************************************************
117 278266 : SUBROUTINE calculate_ptrace_kp(hmat, pmat, ecore, nspin)
118 :
119 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: hmat, pmat
120 : REAL(KIND=dp), INTENT(OUT) :: ecore
121 : INTEGER, INTENT(IN) :: nspin
122 :
123 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_ptrace_kp'
124 :
125 : INTEGER :: handle, ic, ispin, nc
126 : REAL(KIND=dp) :: etr
127 :
128 278266 : CALL timeset(routineN, handle)
129 :
130 278266 : nc = SIZE(pmat, 2)
131 :
132 278266 : ecore = 0.0_dp
133 602448 : DO ispin = 1, nspin
134 1658468 : DO ic = 1, nc
135 1056020 : etr = 0.0_dp
136 1056020 : CALL dbcsr_dot(hmat(1, ic)%matrix, pmat(ispin, ic)%matrix, etr)
137 1380202 : ecore = ecore + etr
138 : END DO
139 : END DO
140 :
141 278266 : CALL timestop(handle)
142 :
143 278266 : END SUBROUTINE calculate_ptrace_kp
144 :
145 : ! **************************************************************************************************
146 : !> \brief Calculate the core Hamiltonian energy which includes the kinetic
147 : !> and the potential energy of the electrons. It is assumed, that
148 : !> the core Hamiltonian matrix h and the density matrix p have the
149 : !> same sparse matrix structure (same atomic blocks and block
150 : !> ordering)
151 : !> \param h ...
152 : !> \param p ...
153 : !> \param ecore ...
154 : !> \date 03.05.2001
155 : !> \par History
156 : !> - simplified taking advantage of new non-redundant matrix
157 : !> structure (27.06.2003,MK)
158 : !> - simplified using DBCSR trace function (21.07.2010, jhu)
159 : !> \author MK
160 : !> \version 1.0
161 : ! **************************************************************************************************
162 28 : SUBROUTINE calculate_ptrace_1(h, p, ecore)
163 :
164 : TYPE(dbcsr_type), POINTER :: h, p
165 : REAL(KIND=dp), INTENT(OUT) :: ecore
166 :
167 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_ptrace_1'
168 :
169 : INTEGER :: handle
170 :
171 28 : CALL timeset(routineN, handle)
172 :
173 28 : ecore = 0.0_dp
174 28 : CALL dbcsr_dot(h, p, ecore)
175 :
176 28 : CALL timestop(handle)
177 :
178 28 : END SUBROUTINE calculate_ptrace_1
179 :
180 : ! **************************************************************************************************
181 : !> \brief Calculate the overlap energy of the core charge distribution.
182 : !> \param qs_env ...
183 : !> \param para_env ...
184 : !> \param calculate_forces ...
185 : !> \param molecular ...
186 : !> \param E_overlap_core ...
187 : !> \param atecc ...
188 : !> \date 30.04.2001
189 : !> \par History
190 : !> - Force calculation added (03.06.2002,MK)
191 : !> - Parallelized using a list of local atoms for rows and
192 : !> columns (19.07.2003,MK)
193 : !> - Use precomputed neighborlists (sab_core) and nl iterator (28.07.2010,jhu)
194 : !> \author MK
195 : !> \version 1.0
196 : ! **************************************************************************************************
197 18169 : SUBROUTINE calculate_ecore_overlap(qs_env, para_env, calculate_forces, molecular, &
198 18169 : E_overlap_core, atecc)
199 : TYPE(qs_environment_type), POINTER :: qs_env
200 : TYPE(mp_para_env_type), POINTER :: para_env
201 : LOGICAL, INTENT(IN) :: calculate_forces
202 : LOGICAL, INTENT(IN), OPTIONAL :: molecular
203 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: E_overlap_core
204 : REAL(KIND=dp), DIMENSION(:), OPTIONAL :: atecc
205 :
206 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_ecore_overlap'
207 :
208 : INTEGER :: atom_a, atom_b, handle, iatom, ikind, &
209 : jatom, jkind, natom, nkind
210 18169 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
211 : LOGICAL :: atenergy, only_molecule, use_virial
212 : REAL(KIND=dp) :: aab, dab, eab, ecore_overlap, f, fab, &
213 : rab2, rootaab, zab
214 18169 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: alpha, radius, zeff
215 : REAL(KIND=dp), DIMENSION(3) :: deab, rab
216 : REAL(KIND=dp), DIMENSION(3, 3) :: pv_loc
217 18169 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
218 : TYPE(atprop_type), POINTER :: atprop
219 : TYPE(mp_comm_type) :: group
220 : TYPE(neighbor_list_iterator_p_type), &
221 18169 : DIMENSION(:), POINTER :: nl_iterator
222 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
223 18169 : POINTER :: sab_core
224 18169 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
225 : TYPE(qs_energy_type), POINTER :: energy
226 18169 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
227 18169 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
228 : TYPE(virial_type), POINTER :: virial
229 :
230 18169 : CALL timeset(routineN, handle)
231 :
232 18169 : NULLIFY (atomic_kind_set)
233 18169 : NULLIFY (qs_kind_set)
234 18169 : NULLIFY (energy)
235 18169 : NULLIFY (atprop)
236 18169 : NULLIFY (force)
237 18169 : NULLIFY (particle_set)
238 :
239 18169 : group = para_env
240 :
241 18169 : only_molecule = .FALSE.
242 18169 : IF (PRESENT(molecular)) only_molecule = molecular
243 :
244 : CALL get_qs_env(qs_env=qs_env, &
245 : atomic_kind_set=atomic_kind_set, &
246 : qs_kind_set=qs_kind_set, &
247 : particle_set=particle_set, &
248 : energy=energy, &
249 : force=force, &
250 : sab_core=sab_core, &
251 : atprop=atprop, &
252 18169 : virial=virial)
253 :
254 : ! Allocate work storage
255 18169 : nkind = SIZE(atomic_kind_set)
256 18169 : natom = SIZE(particle_set)
257 :
258 18169 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
259 :
260 90845 : ALLOCATE (alpha(nkind), radius(nkind), zeff(nkind))
261 50448 : alpha(:) = 0.0_dp
262 50448 : radius(:) = 0.0_dp
263 50448 : zeff(:) = 0.0_dp
264 :
265 18169 : IF (calculate_forces) THEN
266 6157 : CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind)
267 : END IF
268 :
269 18169 : atenergy = .FALSE.
270 18169 : IF (ASSOCIATED(atprop)) THEN
271 18169 : IF (atprop%energy) THEN
272 58 : atenergy = .TRUE.
273 58 : CALL atprop_array_init(atprop%atecc, natom)
274 : END IF
275 : END IF
276 :
277 50448 : DO ikind = 1, nkind
278 : CALL get_qs_kind(qs_kind_set(ikind), &
279 : alpha_core_charge=alpha(ikind), &
280 : core_charge_radius=radius(ikind), &
281 50448 : zeff=zeff(ikind))
282 : END DO
283 :
284 18169 : ecore_overlap = 0.0_dp
285 18169 : pv_loc = 0.0_dp
286 :
287 18169 : CALL neighbor_list_iterator_create(nl_iterator, sab_core)
288 88174 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
289 70005 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, r=rab)
290 70005 : zab = zeff(ikind)*zeff(jkind)
291 70005 : aab = alpha(ikind)*alpha(jkind)/(alpha(ikind) + alpha(jkind))
292 70005 : rootaab = SQRT(aab)
293 70005 : fab = 2.0_dp*oorootpi*zab*rootaab
294 70005 : rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
295 88174 : IF (rab2 > 1.e-8_dp) THEN
296 35743 : IF (ikind == jkind .AND. iatom == jatom) THEN
297 : f = 0.5_dp
298 : ELSE
299 35527 : f = 1.0_dp
300 : END IF
301 35743 : dab = SQRT(rab2)
302 35743 : eab = zab*erfc(rootaab*dab)/dab
303 35743 : ecore_overlap = ecore_overlap + f*eab
304 35743 : IF (atenergy) THEN
305 98 : atprop%atecc(iatom) = atprop%atecc(iatom) + 0.5_dp*f*eab
306 98 : atprop%atecc(jatom) = atprop%atecc(jatom) + 0.5_dp*f*eab
307 : END IF
308 35743 : IF (PRESENT(atecc)) THEN
309 55 : atecc(iatom) = atecc(iatom) + 0.5_dp*f*eab
310 55 : atecc(jatom) = atecc(jatom) + 0.5_dp*f*eab
311 : END IF
312 35743 : IF (calculate_forces) THEN
313 50308 : deab(:) = rab(:)*f*(eab + fab*EXP(-aab*rab2))/rab2
314 12577 : atom_a = atom_of_kind(iatom)
315 12577 : atom_b = atom_of_kind(jatom)
316 50308 : force(ikind)%core_overlap(:, atom_a) = force(ikind)%core_overlap(:, atom_a) + deab(:)
317 50308 : force(jkind)%core_overlap(:, atom_b) = force(jkind)%core_overlap(:, atom_b) - deab(:)
318 12577 : IF (use_virial) THEN
319 969 : CALL virial_pair_force(pv_loc, 1._dp, deab, rab)
320 : END IF
321 : END IF
322 : END IF
323 : END DO
324 18169 : CALL neighbor_list_iterator_release(nl_iterator)
325 :
326 18169 : DEALLOCATE (alpha, radius, zeff)
327 18169 : IF (calculate_forces) THEN
328 6157 : DEALLOCATE (atom_of_kind)
329 : END IF
330 18169 : IF (calculate_forces .AND. use_virial) THEN
331 8502 : virial%pv_ecore_overlap = virial%pv_ecore_overlap + pv_loc
332 8502 : virial%pv_virial = virial%pv_virial + pv_loc
333 : END IF
334 :
335 18169 : CALL group%sum(ecore_overlap)
336 :
337 18169 : energy%core_overlap = ecore_overlap
338 :
339 18169 : IF (PRESENT(E_overlap_core)) THEN
340 2056 : E_overlap_core = energy%core_overlap
341 : END IF
342 :
343 18169 : CALL timestop(handle)
344 :
345 36338 : END SUBROUTINE calculate_ecore_overlap
346 :
347 : ! **************************************************************************************************
348 : !> \brief Calculate the self energy of the core charge distribution.
349 : !> \param qs_env ...
350 : !> \param E_self_core ...
351 : !> \param atecc ...
352 : !> \date 27.04.2001
353 : !> \author MK
354 : !> \version 1.0
355 : ! **************************************************************************************************
356 17737 : SUBROUTINE calculate_ecore_self(qs_env, E_self_core, atecc)
357 : TYPE(qs_environment_type), POINTER :: qs_env
358 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: E_self_core
359 : REAL(KIND=dp), DIMENSION(:), OPTIONAL :: atecc
360 :
361 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_ecore_self'
362 :
363 : INTEGER :: handle, iatom, ikind, iparticle_local, &
364 : natom, nparticle_local
365 : REAL(KIND=dp) :: alpha_core_charge, ecore_self, es, zeff
366 17737 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
367 : TYPE(atprop_type), POINTER :: atprop
368 : TYPE(distribution_1d_type), POINTER :: local_particles
369 17737 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
370 : TYPE(qs_energy_type), POINTER :: energy
371 17737 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
372 :
373 : ! -------------------------------------------------------------------------
374 :
375 17737 : NULLIFY (atprop)
376 17737 : CALL timeset(routineN, handle)
377 :
378 : CALL get_qs_env(qs_env=qs_env, atomic_kind_set=atomic_kind_set, &
379 17737 : qs_kind_set=qs_kind_set, energy=energy, atprop=atprop)
380 :
381 17737 : ecore_self = 0.0_dp
382 :
383 49534 : DO ikind = 1, SIZE(atomic_kind_set)
384 31797 : CALL get_atomic_kind(atomic_kind_set(ikind), natom=natom)
385 31797 : CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff, alpha_core_charge=alpha_core_charge)
386 49534 : ecore_self = ecore_self - REAL(natom, dp)*zeff**2*SQRT(alpha_core_charge)
387 : END DO
388 :
389 17737 : energy%core_self = ecore_self/SQRT(twopi)
390 17737 : IF (PRESENT(E_self_core)) THEN
391 1624 : E_self_core = energy%core_self
392 : END IF
393 :
394 17737 : IF (ASSOCIATED(atprop)) THEN
395 17737 : IF (atprop%energy) THEN
396 : ! atomic energy
397 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, &
398 58 : local_particles=local_particles)
399 58 : natom = SIZE(particle_set)
400 58 : CALL atprop_array_init(atprop%ateself, natom)
401 :
402 172 : DO ikind = 1, SIZE(atomic_kind_set)
403 114 : nparticle_local = local_particles%n_el(ikind)
404 114 : CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff, alpha_core_charge=alpha_core_charge)
405 114 : es = zeff**2*SQRT(alpha_core_charge)/SQRT(twopi)
406 270 : DO iparticle_local = 1, nparticle_local
407 98 : iatom = local_particles%list(ikind)%array(iparticle_local)
408 212 : atprop%ateself(iatom) = atprop%ateself(iatom) - es
409 : END DO
410 : END DO
411 : END IF
412 : END IF
413 17737 : IF (PRESENT(atecc)) THEN
414 : ! atomic energy
415 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, &
416 30 : local_particles=local_particles)
417 30 : natom = SIZE(particle_set)
418 88 : DO ikind = 1, SIZE(atomic_kind_set)
419 58 : nparticle_local = local_particles%n_el(ikind)
420 58 : CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff, alpha_core_charge=alpha_core_charge)
421 58 : es = zeff**2*SQRT(alpha_core_charge)/SQRT(twopi)
422 141 : DO iparticle_local = 1, nparticle_local
423 53 : iatom = local_particles%list(ikind)%array(iparticle_local)
424 111 : atecc(iatom) = atecc(iatom) - es
425 : END DO
426 : END DO
427 : END IF
428 :
429 17737 : CALL timestop(handle)
430 :
431 17737 : END SUBROUTINE calculate_ecore_self
432 :
433 : ! **************************************************************************************************
434 : !> \brief Calculate the overlap and self energy of the core charge distribution for a given alpha
435 : !> Use a minimum image convention and double loop over all atoms
436 : !> \param qs_env ...
437 : !> \param alpha ...
438 : !> \param atecc ...
439 : !> \author JGH
440 : !> \version 1.0
441 : ! **************************************************************************************************
442 2 : SUBROUTINE calculate_ecore_alpha(qs_env, alpha, atecc)
443 : TYPE(qs_environment_type), POINTER :: qs_env
444 : REAL(KIND=dp), INTENT(IN) :: alpha
445 : REAL(KIND=dp), DIMENSION(:) :: atecc
446 :
447 : CHARACTER(len=*), PARAMETER :: routineN = 'calculate_ecore_alpha'
448 :
449 : INTEGER :: handle, iatom, ikind, jatom, jkind, &
450 : natom, nkind
451 2 : INTEGER, ALLOCATABLE, DIMENSION(:) :: kind_of
452 : REAL(KIND=dp) :: dab, eab, fab, rootaab, zab
453 2 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: zeff
454 : REAL(KIND=dp), DIMENSION(3) :: rab
455 2 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
456 : TYPE(cell_type), POINTER :: cell
457 2 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
458 2 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
459 :
460 2 : CALL timeset(routineN, handle)
461 :
462 : CALL get_qs_env(qs_env=qs_env, &
463 : cell=cell, &
464 : atomic_kind_set=atomic_kind_set, &
465 : qs_kind_set=qs_kind_set, &
466 2 : particle_set=particle_set)
467 2 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, kind_of=kind_of)
468 : !
469 2 : nkind = SIZE(atomic_kind_set)
470 2 : natom = SIZE(particle_set)
471 6 : ALLOCATE (zeff(nkind))
472 6 : zeff(:) = 0.0_dp
473 6 : DO ikind = 1, nkind
474 6 : CALL get_qs_kind(qs_kind_set(ikind), zeff=zeff(ikind))
475 : END DO
476 :
477 2 : rootaab = SQRT(0.5_dp*alpha)
478 8 : DO iatom = 1, natom
479 6 : ikind = kind_of(iatom)
480 6 : atecc(iatom) = atecc(iatom) - zeff(ikind)**2*SQRT(alpha/twopi)
481 14 : DO jatom = iatom + 1, natom
482 6 : jkind = kind_of(jatom)
483 6 : zab = zeff(ikind)*zeff(jkind)
484 6 : fab = 2.0_dp*oorootpi*zab*rootaab
485 24 : rab = particle_set(iatom)%r - particle_set(jatom)%r
486 24 : rab = pbc(rab, cell)
487 24 : dab = SQRT(SUM(rab(:)**2))
488 6 : eab = zab*erfc(rootaab*dab)/dab
489 6 : atecc(iatom) = atecc(iatom) + 0.5_dp*eab
490 12 : atecc(jatom) = atecc(jatom) + 0.5_dp*eab
491 : END DO
492 : END DO
493 :
494 2 : DEALLOCATE (zeff)
495 :
496 2 : CALL timestop(handle)
497 :
498 4 : END SUBROUTINE calculate_ecore_alpha
499 :
500 : END MODULE qs_core_energies
|
