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