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