(git:ccc2433)
se_fock_matrix_coulomb.F
Go to the documentation of this file.
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 Module that collects all Coulomb parts of the fock matrix construction
10 !>
11 !> \author Teodoro Laino (05.2009) [tlaino] - split and module reorganization
12 !> \par History
13 !> Teodoro Laino (04.2008) [tlaino] - University of Zurich : d-orbitals
14 !> Teodoro Laino (09.2008) [tlaino] - University of Zurich : Speed-up
15 !> Teodoro Laino (09.2008) [tlaino] - University of Zurich : Periodic SE
16 !> Teodoro Laino (05.2009) [tlaino] - Stress Tensor
17 !>
18 ! **************************************************************************************************
19 MODULE se_fock_matrix_coulomb
20  USE atomic_kind_types, ONLY: atomic_kind_type,&
21  get_atomic_kind_set
22  USE atprop_types, ONLY: atprop_type
23  USE cell_types, ONLY: cell_type
24  USE cp_control_types, ONLY: dft_control_type,&
25  semi_empirical_control_type
26  USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set,&
27  dbcsr_deallocate_matrix_set
28  USE cp_log_handling, ONLY: cp_get_default_logger,&
29  cp_logger_type
30  USE cp_output_handling, ONLY: cp_print_key_finished_output,&
31  cp_print_key_unit_nr
32  USE dbcsr_api, ONLY: dbcsr_add,&
33  dbcsr_distribute,&
34  dbcsr_get_block_diag,&
35  dbcsr_get_block_p,&
36  dbcsr_p_type,&
37  dbcsr_replicate_all,&
38  dbcsr_set,&
39  dbcsr_sum_replicated
40  USE distribution_1d_types, ONLY: distribution_1d_type
41  USE ewald_environment_types, ONLY: ewald_env_get,&
42  ewald_environment_type
43  USE ewald_pw_types, ONLY: ewald_pw_get,&
44  ewald_pw_type
45  USE ewalds_multipole, ONLY: ewald_multipole_evaluate
46  USE fist_neighbor_list_control, ONLY: list_control
47  USE fist_nonbond_env_types, ONLY: fist_nonbond_env_type
48  USE input_constants, ONLY: &
49  do_method_am1, do_method_mndo, do_method_mndod, do_method_pdg, do_method_pm3, &
50  do_method_pm6, do_method_pm6fm, do_method_pnnl, do_method_rm1, do_se_is_slater
51  USE input_section_types, ONLY: section_vals_get_subs_vals,&
52  section_vals_type
53  USE kinds, ONLY: dp
54  USE message_passing, ONLY: mp_para_env_type
55  USE multipole_types, ONLY: do_multipole_charge,&
56  do_multipole_dipole,&
57  do_multipole_none,&
58  do_multipole_quadrupole
59  USE particle_types, ONLY: particle_type
60  USE pw_poisson_types, ONLY: do_ewald_ewald
61  USE qs_energy_types, ONLY: qs_energy_type
62  USE qs_environment_types, ONLY: get_qs_env,&
63  qs_environment_type
64  USE qs_force_types, ONLY: qs_force_type
65  USE qs_kind_types, ONLY: get_qs_kind,&
66  qs_kind_type
67  USE qs_neighbor_list_types, ONLY: get_iterator_info,&
68  neighbor_list_iterate,&
69  neighbor_list_iterator_create,&
70  neighbor_list_iterator_p_type,&
71  neighbor_list_iterator_release,&
72  neighbor_list_set_p_type
73  USE se_fock_matrix_integrals, ONLY: &
74  dfock2c, dfock2c_r3, dfock2_1el, dfock2_1el_r3, fock2c, fock2c_ew, fock2c_r3, fock2_1el, &
75  fock2_1el_ew, fock2_1el_r3, se_coulomb_ij_interaction
76  USE semi_empirical_int_arrays, ONLY: rij_threshold,&
77  se_orbital_pointer
78  USE semi_empirical_integrals, ONLY: corecore_el,&
79  dcorecore_el
80  USE semi_empirical_mpole_methods, ONLY: quadrupole_sph_to_cart
81  USE semi_empirical_mpole_types, ONLY: nddo_mpole_type,&
82  semi_empirical_mpole_type
83  USE semi_empirical_store_int_types, ONLY: semi_empirical_si_type
84  USE semi_empirical_types, ONLY: get_se_param,&
85  se_int_control_type,&
86  se_taper_type,&
87  semi_empirical_p_type,&
88  semi_empirical_type,&
89  setup_se_int_control_type
90  USE semi_empirical_utils, ONLY: finalize_se_taper,&
91  get_se_type,&
92  initialize_se_taper
93  USE virial_methods, ONLY: virial_pair_force
94  USE virial_types, ONLY: virial_type
95 #include "./base/base_uses.f90"
96 
97  IMPLICIT NONE
98  PRIVATE
99 
100  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'se_fock_matrix_coulomb'
101  LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .false.
102 
103  PUBLIC :: build_fock_matrix_coulomb, build_fock_matrix_coulomb_lr, &
104  build_fock_matrix_coul_lr_r3
105 
106 CONTAINS
107 
108 ! **************************************************************************************************
109 !> \brief Construction of the Coulomb part of the Fock matrix
110 !> \param qs_env ...
111 !> \param ks_matrix ...
112 !> \param matrix_p ...
113 !> \param energy ...
114 !> \param calculate_forces ...
115 !> \param store_int_env ...
116 !> \author JGH
117 ! **************************************************************************************************
118  SUBROUTINE build_fock_matrix_coulomb(qs_env, ks_matrix, matrix_p, energy, calculate_forces, &
119  store_int_env)
120 
121  TYPE(qs_environment_type), POINTER :: qs_env
122  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_matrix, matrix_p
123  TYPE(qs_energy_type), POINTER :: energy
124  LOGICAL, INTENT(in) :: calculate_forces
125  TYPE(semi_empirical_si_type), POINTER :: store_int_env
126 
127  CHARACTER(len=*), PARAMETER :: routinen = 'build_fock_matrix_coulomb'
128 
129  INTEGER :: atom_a, atom_b, handle, iatom, ikind, inode, ispin, itype, jatom, jkind, &
130  natorb_a, natorb_a2, natorb_b, natorb_b2, nkind, nspins
131  INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, se_natorb
132  LOGICAL :: anag, atener, check, defined, found, &
133  switch, use_virial
134  LOGICAL, ALLOCATABLE, DIMENSION(:) :: se_defined
135  REAL(kind=dp) :: delta, dr1, ecore2, ecores
136  REAL(kind=dp), DIMENSION(2) :: ecab
137  REAL(kind=dp), DIMENSION(2025) :: pa_a, pa_b, pb_a, pb_b
138  REAL(kind=dp), DIMENSION(3) :: force_ab, rij
139  REAL(kind=dp), DIMENSION(45, 45) :: p_block_tot_a, p_block_tot_b
140  REAL(kind=dp), DIMENSION(:, :), POINTER :: ksa_block_a, ksa_block_b, ksb_block_a, &
141  ksb_block_b, pa_block_a, pa_block_b, &
142  pb_block_a, pb_block_b
143  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
144  TYPE(atprop_type), POINTER :: atprop
145  TYPE(cell_type), POINTER :: cell
146  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: diagmat_ks, diagmat_p
147  TYPE(dft_control_type), POINTER :: dft_control
148  TYPE(ewald_environment_type), POINTER :: ewald_env
149  TYPE(ewald_pw_type), POINTER :: ewald_pw
150  TYPE(mp_para_env_type), POINTER :: para_env
151  TYPE(neighbor_list_iterator_p_type), &
152  DIMENSION(:), POINTER :: nl_iterator
153  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
154  POINTER :: sab_se
155  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
156  TYPE(qs_force_type), DIMENSION(:), POINTER :: force
157  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
158  TYPE(se_int_control_type) :: se_int_control
159  TYPE(se_taper_type), POINTER :: se_taper
160  TYPE(semi_empirical_control_type), POINTER :: se_control
161  TYPE(semi_empirical_p_type), DIMENSION(:), POINTER :: se_kind_list
162  TYPE(semi_empirical_type), POINTER :: se_kind_a, se_kind_b
163  TYPE(virial_type), POINTER :: virial
164 
165  CALL timeset(routinen, handle)
166 
167  NULLIFY (dft_control, cell, force, particle_set, diagmat_ks, diagmat_p, &
168  se_control, se_taper, virial, atprop)
169 
170  CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, cell=cell, se_taper=se_taper, &
171  para_env=para_env, sab_se=sab_se, atomic_kind_set=atomic_kind_set, atprop=atprop, &
172  qs_kind_set=qs_kind_set, particle_set=particle_set, virial=virial)
173 
174  ! Parameters
175  CALL initialize_se_taper(se_taper, coulomb=.true.)
176  se_control => dft_control%qs_control%se_control
177  anag = se_control%analytical_gradients
178  use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
179  atener = atprop%energy
180 
181  CALL setup_se_int_control_type(se_int_control, do_ewald_r3=se_control%do_ewald_r3, &
182  do_ewald_gks=se_control%do_ewald_gks, integral_screening=se_control%integral_screening, &
183  shortrange=(se_control%do_ewald .OR. se_control%do_ewald_gks), &
184  max_multipole=se_control%max_multipole, pc_coulomb_int=.false.)
185  IF (se_control%do_ewald_gks) THEN
186  CALL get_qs_env(qs_env=qs_env, ewald_env=ewald_env, ewald_pw=ewald_pw)
187  CALL ewald_env_get(ewald_env, alpha=se_int_control%ewald_gks%alpha)
188  CALL ewald_pw_get(ewald_pw, pw_big_pool=se_int_control%ewald_gks%pw_pool, &
189  dg=se_int_control%ewald_gks%dg)
190  END IF
191 
192  nspins = dft_control%nspins
193  cpassert(ASSOCIATED(matrix_p))
194  cpassert(SIZE(ks_matrix) > 0)
195 
196  nkind = SIZE(atomic_kind_set)
197  IF (calculate_forces) THEN
198  CALL get_qs_env(qs_env=qs_env, force=force)
199  delta = se_control%delta
200  CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
201  END IF
202 
203  CALL dbcsr_allocate_matrix_set(diagmat_ks, nspins)
204  CALL dbcsr_allocate_matrix_set(diagmat_p, nspins)
205 
206  DO ispin = 1, nspins
207  ! Allocate diagonal block matrices
208  ALLOCATE (diagmat_p(ispin)%matrix, diagmat_ks(ispin)%matrix) !sm->dbcsr
209  CALL dbcsr_get_block_diag(matrix_p(ispin)%matrix, diagmat_p(ispin)%matrix)
210  CALL dbcsr_get_block_diag(ks_matrix(ispin)%matrix, diagmat_ks(ispin)%matrix)
211  CALL dbcsr_set(diagmat_ks(ispin)%matrix, 0.0_dp)
212  CALL dbcsr_replicate_all(diagmat_p(ispin)%matrix)
213  CALL dbcsr_replicate_all(diagmat_ks(ispin)%matrix)
214  END DO
215 
216  ecore2 = 0.0_dp
217  itype = get_se_type(dft_control%qs_control%method_id)
218  ALLOCATE (se_defined(nkind), se_kind_list(nkind), se_natorb(nkind))
219  DO ikind = 1, nkind
220  CALL get_qs_kind(qs_kind_set(ikind), se_parameter=se_kind_a)
221  se_kind_list(ikind)%se_param => se_kind_a
222  CALL get_se_param(se_kind_a, defined=defined, natorb=natorb_a)
223  se_defined(ikind) = (defined .AND. natorb_a >= 1)
224  se_natorb(ikind) = natorb_a
225  END DO
226 
227  CALL neighbor_list_iterator_create(nl_iterator, sab_se)
228  DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
229  CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, inode=inode, r=rij)
230  IF (.NOT. se_defined(ikind)) cycle
231  IF (.NOT. se_defined(jkind)) cycle
232  se_kind_a => se_kind_list(ikind)%se_param
233  se_kind_b => se_kind_list(jkind)%se_param
234  natorb_a = se_natorb(ikind)
235  natorb_b = se_natorb(jkind)
236  natorb_a2 = natorb_a**2
237  natorb_b2 = natorb_b**2
238 
239  IF (inode == 1) THEN
240  CALL dbcsr_get_block_p(matrix=diagmat_p(1)%matrix, &
241  row=iatom, col=iatom, block=pa_block_a, found=found)
242  cpassert(ASSOCIATED(pa_block_a))
243  check = (SIZE(pa_block_a, 1) == natorb_a) .AND. (SIZE(pa_block_a, 2) == natorb_a)
244  cpassert(check)
245  CALL dbcsr_get_block_p(matrix=diagmat_ks(1)%matrix, &
246  row=iatom, col=iatom, block=ksa_block_a, found=found)
247  cpassert(ASSOCIATED(ksa_block_a))
248  p_block_tot_a(1:natorb_a, 1:natorb_a) = 2.0_dp*pa_block_a
249  pa_a(1:natorb_a2) = reshape(pa_block_a, (/natorb_a2/))
250  IF (nspins == 2) THEN
251  CALL dbcsr_get_block_p(matrix=diagmat_p(2)%matrix, &
252  row=iatom, col=iatom, block=pa_block_b, found=found)
253  cpassert(ASSOCIATED(pa_block_b))
254  check = (SIZE(pa_block_b, 1) == natorb_a) .AND. (SIZE(pa_block_b, 2) == natorb_a)
255  cpassert(check)
256  CALL dbcsr_get_block_p(matrix=diagmat_ks(2)%matrix, &
257  row=iatom, col=iatom, block=ksa_block_b, found=found)
258  cpassert(ASSOCIATED(ksa_block_b))
259  p_block_tot_a(1:natorb_a, 1:natorb_a) = pa_block_a + pa_block_b
260  pa_b(1:natorb_a2) = reshape(pa_block_b, (/natorb_a2/))
261  END IF
262 
263  END IF
264 
265  dr1 = dot_product(rij, rij)
266  IF (dr1 > rij_threshold) THEN
267  ! Determine the order of the atoms, and in case switch them..
268  IF (iatom <= jatom) THEN
269  switch = .false.
270  ELSE
271  switch = .true.
272  END IF
273  ! Retrieve blocks for KS and P
274  CALL dbcsr_get_block_p(matrix=diagmat_p(1)%matrix, &
275  row=jatom, col=jatom, block=pb_block_a, found=found)
276  cpassert(ASSOCIATED(pb_block_a))
277  check = (SIZE(pb_block_a, 1) == natorb_b) .AND. (SIZE(pb_block_a, 2) == natorb_b)
278  cpassert(check)
279  CALL dbcsr_get_block_p(matrix=diagmat_ks(1)%matrix, &
280  row=jatom, col=jatom, block=ksb_block_a, found=found)
281  cpassert(ASSOCIATED(ksb_block_a))
282  p_block_tot_b(1:natorb_b, 1:natorb_b) = 2.0_dp*pb_block_a
283  pb_a(1:natorb_b2) = reshape(pb_block_a, (/natorb_b2/))
284  ! Handle more than one configuration
285  IF (nspins == 2) THEN
286  CALL dbcsr_get_block_p(matrix=diagmat_p(2)%matrix, &
287  row=jatom, col=jatom, block=pb_block_b, found=found)
288  cpassert(ASSOCIATED(pb_block_b))
289  check = (SIZE(pb_block_b, 1) == natorb_b) .AND. (SIZE(pb_block_b, 2) == natorb_b)
290  cpassert(check)
291  CALL dbcsr_get_block_p(matrix=diagmat_ks(2)%matrix, &
292  row=jatom, col=jatom, block=ksb_block_b, found=found)
293  cpassert(ASSOCIATED(ksb_block_b))
294  check = (SIZE(pb_block_a, 1) == SIZE(pb_block_b, 1)) .AND. (SIZE(pb_block_a, 2) == SIZE(pb_block_b, 2))
295  cpassert(check)
296  p_block_tot_b(1:natorb_b, 1:natorb_b) = pb_block_a + pb_block_b
297  pb_b(1:natorb_b2) = reshape(pb_block_b, (/natorb_b2/))
298  END IF
299 
300  SELECT CASE (dft_control%qs_control%method_id)
301  CASE (do_method_mndo, do_method_am1, do_method_pm3, do_method_pm6, do_method_pm6fm, do_method_pdg, &
302  do_method_rm1, do_method_mndod, do_method_pnnl)
303 
304  ! Two-centers One-electron terms
305  IF (nspins == 1) THEN
306  ecab = 0._dp
307  CALL fock2_1el(se_kind_a, se_kind_b, rij, ksa_block_a, ksb_block_a, &
308  pa_a, pb_a, ecore=ecab, itype=itype, anag=anag, se_int_control=se_int_control, &
309  se_taper=se_taper, store_int_env=store_int_env)
310  ecore2 = ecore2 + ecab(1) + ecab(2)
311  ELSE IF (nspins == 2) THEN
312  ecab = 0._dp
313  CALL fock2_1el(se_kind_a, se_kind_b, rij, ksa_block_a, ksb_block_a, &
314  pa_block_a, pb_block_a, ecore=ecab, itype=itype, anag=anag, &
315  se_int_control=se_int_control, se_taper=se_taper, store_int_env=store_int_env)
316  CALL fock2_1el(se_kind_a, se_kind_b, rij, ksa_block_b, ksb_block_b, &
317  pa_b, pb_b, ecore=ecab, itype=itype, anag=anag, se_int_control=se_int_control, &
318  se_taper=se_taper, store_int_env=store_int_env)
319  ecore2 = ecore2 + ecab(1) + ecab(2)
320  END IF
321  IF (atener) THEN
322  atprop%atecoul(iatom) = atprop%atecoul(iatom) + ecab(1)
323  atprop%atecoul(jatom) = atprop%atecoul(jatom) + ecab(2)
324  END IF
325  ! Coulomb Terms
326  IF (nspins == 1) THEN
327  CALL fock2c(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, ksa_block_a, p_block_tot_b, &
328  ksb_block_a, factor=0.5_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper, &
329  store_int_env=store_int_env)
330  ELSE IF (nspins == 2) THEN
331  CALL fock2c(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, ksa_block_a, p_block_tot_b, &
332  ksb_block_a, factor=1.0_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper, &
333  store_int_env=store_int_env)
334 
335  CALL fock2c(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, ksa_block_b, p_block_tot_b, &
336  ksb_block_b, factor=1.0_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper, &
337  store_int_env=store_int_env)
338  END IF
339 
340  IF (calculate_forces) THEN
341  atom_a = atom_of_kind(iatom)
342  atom_b = atom_of_kind(jatom)
343 
344  ! Derivatives of the Two-centre One-electron terms
345  force_ab = 0.0_dp
346  IF (nspins == 1) THEN
347  CALL dfock2_1el(se_kind_a, se_kind_b, rij, pa_a, pb_a, itype=itype, anag=anag, &
348  se_int_control=se_int_control, se_taper=se_taper, force=force_ab, &
349  delta=delta)
350  ELSE IF (nspins == 2) THEN
351  CALL dfock2_1el(se_kind_a, se_kind_b, rij, pa_block_a, pb_block_a, itype=itype, anag=anag, &
352  se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta)
353  CALL dfock2_1el(se_kind_a, se_kind_b, rij, pa_b, pb_b, itype=itype, anag=anag, &
354  se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta)
355  END IF
356  IF (use_virial) THEN
357  CALL virial_pair_force(virial%pv_virial, -1.0_dp, force_ab, rij)
358  END IF
359 
360  ! Sum up force components
361  force(ikind)%all_potential(1, atom_a) = force(ikind)%all_potential(1, atom_a) - force_ab(1)
362  force(jkind)%all_potential(1, atom_b) = force(jkind)%all_potential(1, atom_b) + force_ab(1)
363 
364  force(ikind)%all_potential(2, atom_a) = force(ikind)%all_potential(2, atom_a) - force_ab(2)
365  force(jkind)%all_potential(2, atom_b) = force(jkind)%all_potential(2, atom_b) + force_ab(2)
366 
367  force(ikind)%all_potential(3, atom_a) = force(ikind)%all_potential(3, atom_a) - force_ab(3)
368  force(jkind)%all_potential(3, atom_b) = force(jkind)%all_potential(3, atom_b) + force_ab(3)
369 
370  ! Derivatives of the Coulomb Terms
371  force_ab = 0._dp
372  IF (nspins == 1) THEN
373  CALL dfock2c(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, p_block_tot_b, factor=0.25_dp, &
374  anag=anag, se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta)
375  ELSE IF (nspins == 2) THEN
376  CALL dfock2c(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, p_block_tot_b, factor=0.50_dp, &
377  anag=anag, se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta)
378 
379  CALL dfock2c(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, p_block_tot_b, factor=0.50_dp, &
380  anag=anag, se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta)
381  END IF
382  IF (switch) THEN
383  force_ab(1) = -force_ab(1)
384  force_ab(2) = -force_ab(2)
385  force_ab(3) = -force_ab(3)
386  END IF
387  IF (use_virial) THEN
388  CALL virial_pair_force(virial%pv_virial, -1.0_dp, force_ab, rij)
389  END IF
390  ! Sum up force components
391  force(ikind)%rho_elec(1, atom_a) = force(ikind)%rho_elec(1, atom_a) - force_ab(1)
392  force(jkind)%rho_elec(1, atom_b) = force(jkind)%rho_elec(1, atom_b) + force_ab(1)
393 
394  force(ikind)%rho_elec(2, atom_a) = force(ikind)%rho_elec(2, atom_a) - force_ab(2)
395  force(jkind)%rho_elec(2, atom_b) = force(jkind)%rho_elec(2, atom_b) + force_ab(2)
396 
397  force(ikind)%rho_elec(3, atom_a) = force(ikind)%rho_elec(3, atom_a) - force_ab(3)
398  force(jkind)%rho_elec(3, atom_b) = force(jkind)%rho_elec(3, atom_b) + force_ab(3)
399  END IF
400  CASE DEFAULT
401  cpabort("")
402  END SELECT
403  ELSE
404  IF (se_int_control%do_ewald_gks) THEN
405  cpassert(iatom == jatom)
406  ! Two-centers One-electron terms
407  ecores = 0._dp
408  IF (nspins == 1) THEN
409  CALL fock2_1el_ew(se_kind_a, rij, ksa_block_a, pa_a, &
410  ecore=ecores, itype=itype, anag=anag, se_int_control=se_int_control, &
411  se_taper=se_taper, store_int_env=store_int_env)
412  ELSE IF (nspins == 2) THEN
413  CALL fock2_1el_ew(se_kind_a, rij, ksa_block_a, pa_block_a, &
414  ecore=ecores, itype=itype, anag=anag, se_int_control=se_int_control, &
415  se_taper=se_taper, store_int_env=store_int_env)
416  CALL fock2_1el_ew(se_kind_a, rij, ksa_block_b, pa_b, &
417  ecore=ecores, itype=itype, anag=anag, se_int_control=se_int_control, &
418  se_taper=se_taper, store_int_env=store_int_env)
419  END IF
420  ecore2 = ecore2 + ecores
421  IF (atener) THEN
422  atprop%atecoul(iatom) = atprop%atecoul(iatom) + ecores
423  END IF
424  ! Coulomb Terms
425  IF (nspins == 1) THEN
426  CALL fock2c_ew(se_kind_a, rij, p_block_tot_a, ksa_block_a, &
427  factor=0.5_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper, &
428  store_int_env=store_int_env)
429  ELSE IF (nspins == 2) THEN
430  CALL fock2c_ew(se_kind_a, rij, p_block_tot_a, ksa_block_a, &
431  factor=1.0_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper, &
432  store_int_env=store_int_env)
433  CALL fock2c_ew(se_kind_a, rij, p_block_tot_a, ksa_block_b, &
434  factor=1.0_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper, &
435  store_int_env=store_int_env)
436  END IF
437  END IF
438  END IF
439  END DO
440  CALL neighbor_list_iterator_release(nl_iterator)
441 
442  DEALLOCATE (se_kind_list, se_defined, se_natorb)
443 
444  DO ispin = 1, nspins
445  CALL dbcsr_sum_replicated(diagmat_ks(ispin)%matrix)
446  CALL dbcsr_distribute(diagmat_ks(ispin)%matrix)
447  CALL dbcsr_add(ks_matrix(ispin)%matrix, diagmat_ks(ispin)%matrix, &
448  1.0_dp, 1.0_dp)
449  END DO
450  CALL dbcsr_deallocate_matrix_set(diagmat_p)
451  CALL dbcsr_deallocate_matrix_set(diagmat_ks)
452 
453  ! Two-centers one-electron terms
454  CALL para_env%sum(ecore2)
455  energy%hartree = ecore2 - energy%core
456 ! WRITE ( *, * ) 'IN SE_F_COUL', ecore2, energy%core
457 
458  CALL finalize_se_taper(se_taper)
459 
460  CALL timestop(handle)
461  END SUBROUTINE build_fock_matrix_coulomb
462 
463 ! **************************************************************************************************
464 !> \brief Long-Range part for SE Coulomb interactions
465 !> \param qs_env ...
466 !> \param ks_matrix ...
467 !> \param matrix_p ...
468 !> \param energy ...
469 !> \param calculate_forces ...
470 !> \param store_int_env ...
471 !> \date 08.2008 [created]
472 !> \author Teodoro Laino [tlaino] - University of Zurich
473 ! **************************************************************************************************
474  SUBROUTINE build_fock_matrix_coulomb_lr(qs_env, ks_matrix, matrix_p, energy, &
475  calculate_forces, store_int_env)
476 
477  TYPE(qs_environment_type), POINTER :: qs_env
478  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_matrix, matrix_p
479  TYPE(qs_energy_type), POINTER :: energy
480  LOGICAL, INTENT(in) :: calculate_forces
481  TYPE(semi_empirical_si_type), POINTER :: store_int_env
482 
483  CHARACTER(len=*), PARAMETER :: routinen = 'build_fock_matrix_coulomb_lr'
484 
485  INTEGER :: atom_a, ewald_type, forces_g_size, handle, iatom, ikind, ilist, indi, indj, &
486  ispin, itype, iw, jint, natoms, natorb_a, nkind, nlocal_particles, node, nparticle_local, &
487  nspins, size_1c_int
488  INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
489  LOGICAL :: anag, atener, defined, found, use_virial
490  LOGICAL, DIMENSION(3) :: task
491  REAL(kind=dp) :: e_neut, e_self, energy_glob, &
492  energy_local, enuc, fac, tmp
493  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: forces_g, forces_r
494  REAL(kind=dp), DIMENSION(3) :: force_a
495  REAL(kind=dp), DIMENSION(3, 3) :: pv_glob, pv_local, qcart
496  REAL(kind=dp), DIMENSION(5) :: qsph
497  REAL(kind=dp), DIMENSION(:, :), POINTER :: ksa_block_a, pa_block_a
498  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
499  TYPE(atprop_type), POINTER :: atprop
500  TYPE(cell_type), POINTER :: cell
501  TYPE(cp_logger_type), POINTER :: logger
502  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: diagmat_ks, diagmat_p
503  TYPE(dft_control_type), POINTER :: dft_control
504  TYPE(distribution_1d_type), POINTER :: local_particles
505  TYPE(ewald_environment_type), POINTER :: ewald_env
506  TYPE(ewald_pw_type), POINTER :: ewald_pw
507  TYPE(fist_nonbond_env_type), POINTER :: se_nonbond_env
508  TYPE(mp_para_env_type), POINTER :: para_env
509  TYPE(nddo_mpole_type), POINTER :: se_nddo_mpole
510  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
511  TYPE(qs_force_type), DIMENSION(:), POINTER :: force
512  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
513  TYPE(section_vals_type), POINTER :: se_section
514  TYPE(semi_empirical_control_type), POINTER :: se_control
515  TYPE(semi_empirical_mpole_type), POINTER :: mpole
516  TYPE(semi_empirical_type), POINTER :: se_kind_a
517  TYPE(virial_type), POINTER :: virial
518 
519  CALL timeset(routinen, handle)
520 
521  NULLIFY (dft_control, cell, force, particle_set, diagmat_ks, diagmat_p, local_particles, &
522  se_control, ewald_env, ewald_pw, se_nddo_mpole, se_nonbond_env, se_section, mpole, &
523  logger, virial, atprop)
524 
525  logger => cp_get_default_logger()
526  CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, cell=cell, para_env=para_env, &
527  atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, particle_set=particle_set, &
528  ewald_env=ewald_env, &
529  local_particles=local_particles, ewald_pw=ewald_pw, se_nddo_mpole=se_nddo_mpole, &
530  se_nonbond_env=se_nonbond_env, virial=virial, atprop=atprop)
531 
532  nlocal_particles = sum(local_particles%n_el(:))
533  natoms = SIZE(particle_set)
534  CALL ewald_env_get(ewald_env, ewald_type=ewald_type)
535  SELECT CASE (ewald_type)
536  CASE (do_ewald_ewald)
537  forces_g_size = nlocal_particles
538  CASE DEFAULT
539  cpabort("Periodic SE implemented only for standard EWALD sums.")
540  END SELECT
541 
542  ! Parameters
543  se_section => section_vals_get_subs_vals(qs_env%input, "DFT%QS%SE")
544  se_control => dft_control%qs_control%se_control
545  anag = se_control%analytical_gradients
546  use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) .AND. calculate_forces
547  atener = atprop%energy
548 
549  nspins = dft_control%nspins
550  cpassert(ASSOCIATED(matrix_p))
551  cpassert(SIZE(ks_matrix) > 0)
552 
553  CALL dbcsr_allocate_matrix_set(diagmat_ks, nspins)
554  CALL dbcsr_allocate_matrix_set(diagmat_p, nspins)
555 
556  nkind = SIZE(atomic_kind_set)
557  DO ispin = 1, nspins
558  ! Allocate diagonal block matrices
559  ALLOCATE (diagmat_p(ispin)%matrix, diagmat_ks(ispin)%matrix) !sm->dbcsr
560  CALL dbcsr_get_block_diag(matrix_p(ispin)%matrix, diagmat_p(ispin)%matrix)
561  CALL dbcsr_get_block_diag(ks_matrix(ispin)%matrix, diagmat_ks(ispin)%matrix)
562  CALL dbcsr_set(diagmat_ks(ispin)%matrix, 0.0_dp)
563  CALL dbcsr_replicate_all(diagmat_p(ispin)%matrix)
564  CALL dbcsr_replicate_all(diagmat_ks(ispin)%matrix)
565  END DO
566 
567  ! Check for implemented SE methods
568  SELECT CASE (dft_control%qs_control%method_id)
569  CASE (do_method_mndo, do_method_am1, do_method_pm3, do_method_pm6, do_method_pm6fm, do_method_pdg, &
570  do_method_rm1, do_method_mndod, do_method_pnnl)
571  itype = get_se_type(dft_control%qs_control%method_id)
572  CASE DEFAULT
573  cpabort("")
574  END SELECT
575 
576  ! Zero arrays and possibly build neighbor lists
577  energy_local = 0.0_dp
578  energy_glob = 0.0_dp
579  e_neut = 0.0_dp
580  e_self = 0.0_dp
581  task = .false.
582  SELECT CASE (se_control%max_multipole)
583  CASE (do_multipole_none)
584  ! Do Nothing
585  CASE (do_multipole_charge)
586  task(1) = .true.
587  CASE (do_multipole_dipole)
588  task = .true.
589  task(3) = .false.
590  CASE (do_multipole_quadrupole)
591  task = .true.
592  CASE DEFAULT
593  cpabort("")
594  END SELECT
595 
596  ! Build-up neighbor lists for real-space part of Ewald multipoles
597  CALL list_control(atomic_kind_set, particle_set, local_particles, &
598  cell, se_nonbond_env, para_env, se_section)
599 
600  enuc = 0.0_dp
601  energy%core_overlap = 0.0_dp
602  se_nddo_mpole%charge = 0.0_dp
603  se_nddo_mpole%dipole = 0.0_dp
604  se_nddo_mpole%quadrupole = 0.0_dp
605 
606  DO ispin = 1, nspins
607  ! Compute the NDDO mpole expansion
608  DO ikind = 1, nkind
609  CALL get_qs_kind(qs_kind_set(ikind), se_parameter=se_kind_a)
610  CALL get_se_param(se_kind_a, defined=defined, natorb=natorb_a)
611 
612  IF (.NOT. defined .OR. natorb_a < 1) cycle
613 
614  nparticle_local = local_particles%n_el(ikind)
615  DO ilist = 1, nparticle_local
616  iatom = local_particles%list(ikind)%array(ilist)
617  CALL dbcsr_get_block_p(matrix=diagmat_p(ispin)%matrix, &
618  row=iatom, col=iatom, block=pa_block_a, found=found)
619  cpassert(ASSOCIATED(pa_block_a))
620  ! Nuclei
621  IF (task(1) .AND. ispin == 1) se_nddo_mpole%charge(iatom) = se_kind_a%zeff
622  ! Electrons
623  size_1c_int = SIZE(se_kind_a%w_mpole)
624  DO jint = 1, size_1c_int
625  mpole => se_kind_a%w_mpole(jint)%mpole
626  indi = se_orbital_pointer(mpole%indi)
627  indj = se_orbital_pointer(mpole%indj)
628  fac = 1.0_dp
629  IF (indi /= indj) fac = 2.0_dp
630 
631  ! Charge
632  IF (mpole%task(1) .AND. task(1)) THEN
633  se_nddo_mpole%charge(iatom) = se_nddo_mpole%charge(iatom) + &
634  fac*pa_block_a(indi, indj)*mpole%c
635  END IF
636 
637  ! Dipole
638  IF (mpole%task(2) .AND. task(2)) THEN
639  se_nddo_mpole%dipole(:, iatom) = se_nddo_mpole%dipole(:, iatom) + &
640  fac*pa_block_a(indi, indj)*mpole%d(:)
641  END IF
642 
643  ! Quadrupole
644  IF (mpole%task(3) .AND. task(3)) THEN
645  qsph = fac*mpole%qs*pa_block_a(indi, indj)
646  CALL quadrupole_sph_to_cart(qcart, qsph)
647  se_nddo_mpole%quadrupole(:, :, iatom) = se_nddo_mpole%quadrupole(:, :, iatom) + &
648  qcart
649  END IF
650  END DO
651  ! Print some info about charge, dipole and quadrupole (debug purpose only)
652  IF (debug_this_module) THEN
653  WRITE (*, '(I5,F12.6,5X,3F12.6,5X,9F12.6)') iatom, se_nddo_mpole%charge(iatom), &
654  se_nddo_mpole%dipole(:, iatom), se_nddo_mpole%quadrupole(:, :, iatom)
655  END IF
656  END DO
657  END DO
658  END DO
659  CALL para_env%sum(se_nddo_mpole%charge)
660  CALL para_env%sum(se_nddo_mpole%dipole)
661  CALL para_env%sum(se_nddo_mpole%quadrupole)
662 
663  ! Initialize for virial
664  IF (use_virial) THEN
665  pv_glob = 0.0_dp
666  pv_local = 0.0_dp
667  END IF
668 
669  ! Ewald Multipoles Sum
670  iw = cp_print_key_unit_nr(logger, se_section, "PRINT%EWALD_INFO", extension=".seLog")
671  IF (calculate_forces) THEN
672  CALL get_qs_env(qs_env=qs_env, force=force)
673  CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
674 
675  ! Allocate and zeroing arrays
676  ALLOCATE (forces_g(3, forces_g_size))
677  ALLOCATE (forces_r(3, natoms))
678  forces_g = 0.0_dp
679  forces_r = 0.0_dp
680  CALL ewald_multipole_evaluate( &
681  ewald_env, ewald_pw, se_nonbond_env, cell, &
682  particle_set, local_particles, energy_local, energy_glob, e_neut, e_self, task, &
683  do_correction_bonded=.false., do_forces=.true., do_stress=use_virial, do_efield=.true., &
684  charges=se_nddo_mpole%charge, dipoles=se_nddo_mpole%dipole, quadrupoles=se_nddo_mpole%quadrupole, &
685  forces_local=forces_g, forces_glob=forces_r, pv_glob=pv_glob, pv_local=pv_local, &
686  efield0=se_nddo_mpole%efield0, efield1=se_nddo_mpole%efield1, efield2=se_nddo_mpole%efield2, iw=iw, &
687  do_debug=.true.)
688  ! Only SR force have to be summed up.. the one in g-space are already fully local..
689  CALL para_env%sum(forces_r)
690  ELSE
691  CALL ewald_multipole_evaluate( &
692  ewald_env, ewald_pw, se_nonbond_env, cell, &
693  particle_set, local_particles, energy_local, energy_glob, e_neut, e_self, task, &
694  do_correction_bonded=.false., do_forces=.false., do_stress=.false., do_efield=.true., &
695  charges=se_nddo_mpole%charge, dipoles=se_nddo_mpole%dipole, quadrupoles=se_nddo_mpole%quadrupole, &
696  efield0=se_nddo_mpole%efield0, efield1=se_nddo_mpole%efield1, efield2=se_nddo_mpole%efield2, &
697  iw=iw, do_debug=.true.)
698  END IF
699  CALL cp_print_key_finished_output(iw, logger, se_section, "PRINT%EWALD_INFO")
700 
701  ! Apply correction only when the Integral Scheme is different from Slater
702  IF ((se_control%integral_screening /= do_se_is_slater) .AND. (.NOT. debug_this_module)) THEN
703  CALL build_fock_matrix_coul_lrc(qs_env, ks_matrix, matrix_p, energy, calculate_forces, &
704  store_int_env, se_nddo_mpole, task, diagmat_p, diagmat_ks, virial, &
705  pv_glob)
706  END IF
707 
708  ! Virial for the long-range part and correction
709  IF (use_virial) THEN
710  ! Sum up contribution of pv_glob on each thread and keep only one copy of pv_local
711  virial%pv_virial = virial%pv_virial + pv_glob
712  IF (para_env%is_source()) THEN
713  virial%pv_virial = virial%pv_virial + pv_local
714  END IF
715  END IF
716 
717  ! Debug Statements
718  IF (debug_this_module) THEN
719  CALL para_env%sum(energy_glob)
720  WRITE (*, *) "TOTAL ENERGY AFTER EWALD:", energy_local + energy_glob + e_neut + e_self, &
721  energy_local, energy_glob, e_neut, e_self
722  END IF
723 
724  ! Modify the KS matrix and possibly compute derivatives
725  node = 0
726  DO ikind = 1, nkind
727  CALL get_qs_kind(qs_kind_set(ikind), se_parameter=se_kind_a)
728  CALL get_se_param(se_kind_a, defined=defined, natorb=natorb_a)
729 
730  IF (.NOT. defined .OR. natorb_a < 1) cycle
731 
732  nparticle_local = local_particles%n_el(ikind)
733  DO ilist = 1, nparticle_local
734  node = node + 1
735  iatom = local_particles%list(ikind)%array(ilist)
736  DO ispin = 1, nspins
737  CALL dbcsr_get_block_p(matrix=diagmat_ks(ispin)%matrix, &
738  row=iatom, col=iatom, block=ksa_block_a, found=found)
739  cpassert(ASSOCIATED(ksa_block_a))
740 
741  ! Modify Hamiltonian Matrix accordingly potential, field and electric field gradient
742  size_1c_int = SIZE(se_kind_a%w_mpole)
743  DO jint = 1, size_1c_int
744  tmp = 0.0_dp
745  mpole => se_kind_a%w_mpole(jint)%mpole
746  indi = se_orbital_pointer(mpole%indi)
747  indj = se_orbital_pointer(mpole%indj)
748 
749  ! Charge
750  IF (mpole%task(1) .AND. task(1)) THEN
751  tmp = tmp + mpole%c*se_nddo_mpole%efield0(iatom)
752  END IF
753 
754  ! Dipole
755  IF (mpole%task(2) .AND. task(2)) THEN
756  tmp = tmp - dot_product(mpole%d, se_nddo_mpole%efield1(:, iatom))
757  END IF
758 
759  ! Quadrupole
760  IF (mpole%task(3) .AND. task(3)) THEN
761  tmp = tmp - (1.0_dp/3.0_dp)*sum(mpole%qc*reshape(se_nddo_mpole%efield2(:, iatom), (/3, 3/)))
762  END IF
763 
764  ksa_block_a(indi, indj) = ksa_block_a(indi, indj) + tmp
765  ksa_block_a(indj, indi) = ksa_block_a(indi, indj)
766  END DO
767  END DO
768 
769  ! Nuclear term and forces
770  IF (task(1)) enuc = enuc + se_kind_a%zeff*se_nddo_mpole%efield0(iatom)
771  IF (atener) THEN
772  atprop%atecoul(iatom) = atprop%atecoul(iatom) + &
773  0.5_dp*se_kind_a%zeff*se_nddo_mpole%efield0(iatom)
774  END IF
775  IF (calculate_forces) THEN
776  atom_a = atom_of_kind(iatom)
777  force_a = forces_r(1:3, iatom) + forces_g(1:3, node)
778  ! Derivatives of the periodic Coulomb Terms
779  force(ikind)%all_potential(:, atom_a) = force(ikind)%all_potential(:, atom_a) - force_a(:)
780  END IF
781  END DO
782  END DO
783  ! Sum nuclear energy contribution
784  CALL para_env%sum(enuc)
785  energy%core_overlap = energy%core_overlap + energy%core_overlap0 + 0.5_dp*enuc
786 
787  ! Debug Statements
788  IF (debug_this_module) THEN
789  WRITE (*, *) "ENUC: ", enuc*0.5_dp
790  END IF
791 
792  DO ispin = 1, nspins
793  CALL dbcsr_sum_replicated(diagmat_ks(ispin)%matrix)
794  CALL dbcsr_distribute(diagmat_ks(ispin)%matrix)
795  CALL dbcsr_add(ks_matrix(ispin)%matrix, diagmat_ks(ispin)%matrix, &
796  1.0_dp, 1.0_dp)
797  END DO
798  CALL dbcsr_deallocate_matrix_set(diagmat_p)
799  CALL dbcsr_deallocate_matrix_set(diagmat_ks)
800 
801  ! Set the Fock matrix contribution to SCP
802  IF (calculate_forces) THEN
803  DEALLOCATE (forces_g)
804  DEALLOCATE (forces_r)
805  END IF
806 
807  CALL timestop(handle)
808  END SUBROUTINE build_fock_matrix_coulomb_lr
809 
810 ! **************************************************************************************************
811 !> \brief When doing long-range SE calculation, this module computes the correction
812 !> between the mismatch of point-like multipoles and multipoles represented
813 !> with charges
814 !> \param qs_env ...
815 !> \param ks_matrix ...
816 !> \param matrix_p ...
817 !> \param energy ...
818 !> \param calculate_forces ...
819 !> \param store_int_env ...
820 !> \param se_nddo_mpole ...
821 !> \param task ...
822 !> \param diagmat_p ...
823 !> \param diagmat_ks ...
824 !> \param virial ...
825 !> \param pv_glob ...
826 !> \author Teodoro Laino [tlaino] - 05.2009
827 ! **************************************************************************************************
828  SUBROUTINE build_fock_matrix_coul_lrc(qs_env, ks_matrix, matrix_p, energy, &
829  calculate_forces, store_int_env, se_nddo_mpole, task, diagmat_p, diagmat_ks, &
830  virial, pv_glob)
831 
832  TYPE(qs_environment_type), POINTER :: qs_env
833  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_matrix, matrix_p
834  TYPE(qs_energy_type), POINTER :: energy
835  LOGICAL, INTENT(in) :: calculate_forces
836  TYPE(semi_empirical_si_type), POINTER :: store_int_env
837  TYPE(nddo_mpole_type), POINTER :: se_nddo_mpole
838  LOGICAL, DIMENSION(3), INTENT(IN) :: task
839  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: diagmat_p, diagmat_ks
840  TYPE(virial_type), POINTER :: virial
841  REAL(kind=dp), DIMENSION(3, 3), INTENT(INOUT) :: pv_glob
842 
843  CHARACTER(len=*), PARAMETER :: routinen = 'build_fock_matrix_coul_lrc'
844 
845  INTEGER :: atom_a, atom_b, handle, iatom, ikind, inode, itype, jatom, jkind, natorb_a, &
846  natorb_a2, natorb_b, natorb_b2, nkind, nspins, size1, size2
847  INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, se_natorb
848  LOGICAL :: anag, atener, check, defined, found, &
849  switch, use_virial
850  LOGICAL, ALLOCATABLE, DIMENSION(:) :: se_defined
851  REAL(kind=dp) :: delta, dr1, ecore2, enuc, enuclear, &
852  ptens11, ptens12, ptens13, ptens21, &
853  ptens22, ptens23, ptens31, ptens32, &
854  ptens33
855  REAL(kind=dp), DIMENSION(2) :: ecab
856  REAL(kind=dp), DIMENSION(2025) :: pa_a, pa_b, pb_a, pb_b
857  REAL(kind=dp), DIMENSION(3) :: force_ab, force_ab0, rij
858  REAL(kind=dp), DIMENSION(45, 45) :: p_block_tot_a, p_block_tot_b
859  REAL(kind=dp), DIMENSION(:), POINTER :: efield0
860  REAL(kind=dp), DIMENSION(:, :), POINTER :: efield1, efield2, ksa_block_a, ksa_block_b, &
861  ksb_block_a, ksb_block_b, pa_block_a, pa_block_b, pb_block_a, pb_block_b
862  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
863  TYPE(atprop_type), POINTER :: atprop
864  TYPE(cell_type), POINTER :: cell
865  TYPE(dft_control_type), POINTER :: dft_control
866  TYPE(mp_para_env_type), POINTER :: para_env
867  TYPE(neighbor_list_iterator_p_type), &
868  DIMENSION(:), POINTER :: nl_iterator
869  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
870  POINTER :: sab_lrc
871  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
872  TYPE(qs_force_type), DIMENSION(:), POINTER :: force
873  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
874  TYPE(se_int_control_type) :: se_int_control
875  TYPE(se_taper_type), POINTER :: se_taper
876  TYPE(semi_empirical_control_type), POINTER :: se_control
877  TYPE(semi_empirical_p_type), DIMENSION(:), POINTER :: se_kind_list
878  TYPE(semi_empirical_type), POINTER :: se_kind_a, se_kind_b
879 
880  CALL timeset(routinen, handle)
881  NULLIFY (dft_control, cell, force, particle_set, se_control, se_taper, &
882  efield0, efield1, efield2, atprop)
883 
884  CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, cell=cell, se_taper=se_taper, &
885  para_env=para_env, sab_lrc=sab_lrc, atomic_kind_set=atomic_kind_set, &
886  qs_kind_set=qs_kind_set, particle_set=particle_set, atprop=atprop)
887 
888  ! Parameters
889  CALL initialize_se_taper(se_taper, lr_corr=.true.)
890  se_control => dft_control%qs_control%se_control
891  anag = se_control%analytical_gradients
892  use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) .AND. calculate_forces
893  atener = atprop%energy
894 
895  CALL setup_se_int_control_type(se_int_control, do_ewald_r3=se_control%do_ewald_r3, &
896  do_ewald_gks=.false., integral_screening=se_control%integral_screening, &
897  shortrange=.false., max_multipole=se_control%max_multipole, &
898  pc_coulomb_int=.true.)
899 
900  nspins = dft_control%nspins
901  cpassert(ASSOCIATED(matrix_p))
902  cpassert(SIZE(ks_matrix) > 0)
903  cpassert(ASSOCIATED(diagmat_p))
904  cpassert(ASSOCIATED(diagmat_ks))
905  mark_used(ks_matrix)
906  mark_used(matrix_p)
907 
908  nkind = SIZE(atomic_kind_set)
909  IF (calculate_forces) THEN
910  CALL get_qs_env(qs_env=qs_env, force=force)
911  delta = se_control%delta
912  CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
913  END IF
914 
915  ! Allocate arrays for storing partial information on potential, field, field gradient
916  size1 = SIZE(se_nddo_mpole%efield0)
917  ALLOCATE (efield0(size1))
918  efield0 = 0.0_dp
919  size1 = SIZE(se_nddo_mpole%efield1, 1)
920  size2 = SIZE(se_nddo_mpole%efield1, 2)
921  ALLOCATE (efield1(size1, size2))
922  efield1 = 0.0_dp
923  size1 = SIZE(se_nddo_mpole%efield2, 1)
924  size2 = SIZE(se_nddo_mpole%efield2, 2)
925  ALLOCATE (efield2(size1, size2))
926  efield2 = 0.0_dp
927 
928  ! Initialize if virial is requested
929  IF (use_virial) THEN
930  ptens11 = 0.0_dp; ptens12 = 0.0_dp; ptens13 = 0.0_dp
931  ptens21 = 0.0_dp; ptens22 = 0.0_dp; ptens23 = 0.0_dp
932  ptens31 = 0.0_dp; ptens32 = 0.0_dp; ptens33 = 0.0_dp
933  END IF
934 
935  ! Start of the loop for the correction of the pair interactions
936  ecore2 = 0.0_dp
937  enuclear = 0.0_dp
938  itype = get_se_type(dft_control%qs_control%method_id)
939 
940  ALLOCATE (se_defined(nkind), se_kind_list(nkind), se_natorb(nkind))
941  DO ikind = 1, nkind
942  CALL get_qs_kind(qs_kind_set(ikind), se_parameter=se_kind_a)
943  se_kind_list(ikind)%se_param => se_kind_a
944  CALL get_se_param(se_kind_a, defined=defined, natorb=natorb_a)
945  se_defined(ikind) = (defined .AND. natorb_a >= 1)
946  se_natorb(ikind) = natorb_a
947  END DO
948 
949  CALL neighbor_list_iterator_create(nl_iterator, sab_lrc)
950  DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
951  CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, inode=inode, r=rij)
952  IF (.NOT. se_defined(ikind)) cycle
953  IF (.NOT. se_defined(jkind)) cycle
954  se_kind_a => se_kind_list(ikind)%se_param
955  se_kind_b => se_kind_list(jkind)%se_param
956  natorb_a = se_natorb(ikind)
957  natorb_b = se_natorb(jkind)
958  natorb_a2 = natorb_a**2
959  natorb_b2 = natorb_b**2
960 
961  IF (inode == 1) THEN
962  CALL dbcsr_get_block_p(matrix=diagmat_p(1)%matrix, &
963  row=iatom, col=iatom, block=pa_block_a, found=found)
964  cpassert(ASSOCIATED(pa_block_a))
965  check = (SIZE(pa_block_a, 1) == natorb_a) .AND. (SIZE(pa_block_a, 2) == natorb_a)
966  cpassert(check)
967  CALL dbcsr_get_block_p(matrix=diagmat_ks(1)%matrix, &
968  row=iatom, col=iatom, block=ksa_block_a, found=found)
969  cpassert(ASSOCIATED(ksa_block_a))
970  p_block_tot_a(1:natorb_a, 1:natorb_a) = 2.0_dp*pa_block_a
971  pa_a(1:natorb_a2) = reshape(pa_block_a, (/natorb_a2/))
972  IF (nspins == 2) THEN
973  CALL dbcsr_get_block_p(matrix=diagmat_p(2)%matrix, &
974  row=iatom, col=iatom, block=pa_block_b, found=found)
975  cpassert(ASSOCIATED(pa_block_b))
976  check = (SIZE(pa_block_b, 1) == natorb_a) .AND. (SIZE(pa_block_b, 2) == natorb_a)
977  cpassert(check)
978  CALL dbcsr_get_block_p(matrix=diagmat_ks(2)%matrix, &
979  row=iatom, col=iatom, block=ksa_block_b, found=found)
980  cpassert(ASSOCIATED(ksa_block_b))
981  p_block_tot_a(1:natorb_a, 1:natorb_a) = pa_block_a + pa_block_b
982  pa_b(1:natorb_a2) = reshape(pa_block_b, (/natorb_a2/))
983  END IF
984 
985  END IF
986 
987  dr1 = dot_product(rij, rij)
988  IF (dr1 > rij_threshold) THEN
989  ! Determine the order of the atoms, and in case switch them..
990  IF (iatom <= jatom) THEN
991  switch = .false.
992  ELSE
993  switch = .true.
994  END IF
995 
996  ! Point-like interaction corrections
997  CALL se_coulomb_ij_interaction(iatom, jatom, task, do_forces=calculate_forces, &
998  do_efield=.true., do_stress=use_virial, charges=se_nddo_mpole%charge, &
999  dipoles=se_nddo_mpole%dipole, quadrupoles=se_nddo_mpole%quadrupole, &
1000  force_ab=force_ab0, efield0=efield0, efield1=efield1, efield2=efield2, &
1001  rab2=dr1, rab=rij, ptens11=ptens11, ptens12=ptens12, ptens13=ptens13, &
1002  ptens21=ptens21, ptens22=ptens22, ptens23=ptens23, ptens31=ptens31, &
1003  ptens32=ptens32, ptens33=ptens33)
1004 
1005  ! Retrieve blocks for KS and P
1006  CALL dbcsr_get_block_p(matrix=diagmat_p(1)%matrix, &
1007  row=jatom, col=jatom, block=pb_block_a, found=found)
1008  cpassert(ASSOCIATED(pb_block_a))
1009  check = (SIZE(pb_block_a, 1) == natorb_b) .AND. (SIZE(pb_block_a, 2) == natorb_b)
1010  cpassert(check)
1011  CALL dbcsr_get_block_p(matrix=diagmat_ks(1)%matrix, &
1012  row=jatom, col=jatom, block=ksb_block_a, found=found)
1013  cpassert(ASSOCIATED(ksb_block_a))
1014  p_block_tot_b(1:natorb_b, 1:natorb_b) = 2.0_dp*pb_block_a
1015  pb_a(1:natorb_b2) = reshape(pb_block_a, (/natorb_b2/))
1016  ! Handle more than one configuration
1017  IF (nspins == 2) THEN
1018  CALL dbcsr_get_block_p(matrix=diagmat_p(2)%matrix, &
1019  row=jatom, col=jatom, block=pb_block_b, found=found)
1020  cpassert(ASSOCIATED(pb_block_b))
1021  check = (SIZE(pb_block_b, 1) == natorb_b) .AND. (SIZE(pb_block_b, 2) == natorb_b)
1022  cpassert(check)
1023  CALL dbcsr_get_block_p(matrix=diagmat_ks(2)%matrix, &
1024  row=jatom, col=jatom, block=ksb_block_b, found=found)
1025  cpassert(ASSOCIATED(ksb_block_b))
1026  check = (SIZE(pb_block_a, 1) == SIZE(pb_block_b, 1)) .AND. (SIZE(pb_block_a, 2) == SIZE(pb_block_b, 2))
1027  cpassert(check)
1028  p_block_tot_b(1:natorb_b, 1:natorb_b) = pb_block_a + pb_block_b
1029  pb_b(1:natorb_b2) = reshape(pb_block_b, (/natorb_b2/))
1030  END IF
1031 
1032  SELECT CASE (dft_control%qs_control%method_id)
1033  CASE (do_method_mndo, do_method_am1, do_method_pm3, do_method_pm6, do_method_pm6fm, do_method_pdg, &
1034  do_method_rm1, do_method_mndod)
1035  ! Evaluate nuclear contribution..
1036  CALL corecore_el(se_kind_a, se_kind_b, rij, enuc=enuc, itype=itype, anag=anag, &
1037  se_int_control=se_int_control, se_taper=se_taper)
1038  enuclear = enuclear + enuc
1039 
1040  ! Two-centers One-electron terms
1041  IF (nspins == 1) THEN
1042  ecab = 0._dp
1043  CALL fock2_1el(se_kind_a, se_kind_b, rij, ksa_block_a, ksb_block_a, &
1044  pa_a, pb_a, ecore=ecab, itype=itype, anag=anag, se_int_control=se_int_control, &
1045  se_taper=se_taper, store_int_env=store_int_env)
1046  ecore2 = ecore2 + ecab(1) + ecab(2)
1047  ELSE IF (nspins == 2) THEN
1048  ecab = 0._dp
1049  CALL fock2_1el(se_kind_a, se_kind_b, rij, ksa_block_a, ksb_block_a, &
1050  pa_block_a, pb_block_a, ecore=ecab, itype=itype, anag=anag, &
1051  se_int_control=se_int_control, se_taper=se_taper, store_int_env=store_int_env)
1052  CALL fock2_1el(se_kind_a, se_kind_b, rij, ksa_block_b, ksb_block_b, &
1053  pa_b, pb_b, ecore=ecab, itype=itype, anag=anag, se_int_control=se_int_control, &
1054  se_taper=se_taper, store_int_env=store_int_env)
1055  ecore2 = ecore2 + ecab(1) + ecab(2)
1056  END IF
1057  IF (atener) THEN
1058  atprop%atecoul(iatom) = atprop%atecoul(iatom) + ecab(1) + 0.5_dp*enuc
1059  atprop%atecoul(jatom) = atprop%atecoul(jatom) + ecab(2) + 0.5_dp*enuc
1060  END IF
1061  ! Coulomb Terms
1062  IF (nspins == 1) THEN
1063  CALL fock2c(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, ksa_block_a, p_block_tot_b, &
1064  ksb_block_a, factor=0.5_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper, &
1065  store_int_env=store_int_env)
1066  ELSE IF (nspins == 2) THEN
1067  CALL fock2c(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, ksa_block_a, p_block_tot_b, &
1068  ksb_block_a, factor=1.0_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper, &
1069  store_int_env=store_int_env)
1070 
1071  CALL fock2c(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, ksa_block_b, p_block_tot_b, &
1072  ksb_block_b, factor=1.0_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper, &
1073  store_int_env=store_int_env)
1074  END IF
1075 
1076  IF (calculate_forces) THEN
1077  atom_a = atom_of_kind(iatom)
1078  atom_b = atom_of_kind(jatom)
1079 
1080  ! Evaluate nuclear contribution..
1081  CALL dcorecore_el(se_kind_a, se_kind_b, rij, denuc=force_ab, itype=itype, delta=delta, &
1082  anag=anag, se_int_control=se_int_control, se_taper=se_taper)
1083 
1084  ! Derivatives of the Two-centre One-electron terms
1085  IF (nspins == 1) THEN
1086  CALL dfock2_1el(se_kind_a, se_kind_b, rij, pa_a, pb_a, itype=itype, anag=anag, &
1087  se_int_control=se_int_control, se_taper=se_taper, force=force_ab, &
1088  delta=delta)
1089  ELSE IF (nspins == 2) THEN
1090  CALL dfock2_1el(se_kind_a, se_kind_b, rij, pa_block_a, pb_block_a, itype=itype, anag=anag, &
1091  se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta)
1092  CALL dfock2_1el(se_kind_a, se_kind_b, rij, pa_b, pb_b, itype=itype, anag=anag, &
1093  se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta)
1094  END IF
1095  IF (use_virial) THEN
1096  CALL virial_pair_force(virial%pv_virial, -1.0_dp, force_ab, rij)
1097  END IF
1098  force_ab = force_ab + force_ab0
1099 
1100  ! Sum up force components
1101  force(ikind)%all_potential(1, atom_a) = force(ikind)%all_potential(1, atom_a) - force_ab(1)
1102  force(jkind)%all_potential(1, atom_b) = force(jkind)%all_potential(1, atom_b) + force_ab(1)
1103 
1104  force(ikind)%all_potential(2, atom_a) = force(ikind)%all_potential(2, atom_a) - force_ab(2)
1105  force(jkind)%all_potential(2, atom_b) = force(jkind)%all_potential(2, atom_b) + force_ab(2)
1106 
1107  force(ikind)%all_potential(3, atom_a) = force(ikind)%all_potential(3, atom_a) - force_ab(3)
1108  force(jkind)%all_potential(3, atom_b) = force(jkind)%all_potential(3, atom_b) + force_ab(3)
1109 
1110  ! Derivatives of the Coulomb Terms
1111  force_ab = 0._dp
1112  IF (nspins == 1) THEN
1113  CALL dfock2c(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, p_block_tot_b, factor=0.25_dp, &
1114  anag=anag, se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta)
1115  ELSE IF (nspins == 2) THEN
1116  CALL dfock2c(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, p_block_tot_b, factor=0.50_dp, &
1117  anag=anag, se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta)
1118 
1119  CALL dfock2c(se_kind_a, se_kind_b, rij, switch, p_block_tot_a, p_block_tot_b, factor=0.50_dp, &
1120  anag=anag, se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta)
1121  END IF
1122  IF (switch) THEN
1123  force_ab(1) = -force_ab(1)
1124  force_ab(2) = -force_ab(2)
1125  force_ab(3) = -force_ab(3)
1126  END IF
1127  IF (use_virial) THEN
1128  CALL virial_pair_force(virial%pv_virial, -1.0_dp, force_ab, rij)
1129  END IF
1130 
1131  ! Sum up force components
1132  force(ikind)%rho_elec(1, atom_a) = force(ikind)%rho_elec(1, atom_a) - force_ab(1)
1133  force(jkind)%rho_elec(1, atom_b) = force(jkind)%rho_elec(1, atom_b) + force_ab(1)
1134 
1135  force(ikind)%rho_elec(2, atom_a) = force(ikind)%rho_elec(2, atom_a) - force_ab(2)
1136  force(jkind)%rho_elec(2, atom_b) = force(jkind)%rho_elec(2, atom_b) + force_ab(2)
1137 
1138  force(ikind)%rho_elec(3, atom_a) = force(ikind)%rho_elec(3, atom_a) - force_ab(3)
1139  force(jkind)%rho_elec(3, atom_b) = force(jkind)%rho_elec(3, atom_b) + force_ab(3)
1140  END IF
1141  CASE DEFAULT
1142  cpabort("")
1143  END SELECT
1144  END IF
1145  END DO
1146  CALL neighbor_list_iterator_release(nl_iterator)
1147 
1148  DEALLOCATE (se_kind_list, se_defined, se_natorb)
1149 
1150  ! Sum-up Virial constribution (long-range correction)
1151  IF (use_virial) THEN
1152  pv_glob(1, 1) = pv_glob(1, 1) + ptens11
1153  pv_glob(1, 2) = pv_glob(1, 2) + (ptens12 + ptens21)*0.5_dp
1154  pv_glob(1, 3) = pv_glob(1, 3) + (ptens13 + ptens31)*0.5_dp
1155  pv_glob(2, 1) = pv_glob(1, 2)
1156  pv_glob(2, 2) = pv_glob(2, 2) + ptens22
1157  pv_glob(2, 3) = pv_glob(2, 3) + (ptens23 + ptens32)*0.5_dp
1158  pv_glob(3, 1) = pv_glob(1, 3)
1159  pv_glob(3, 2) = pv_glob(2, 3)
1160  pv_glob(3, 3) = pv_glob(3, 3) + ptens33
1161  END IF
1162 
1163  ! collect information on potential, field, field gradient
1164  CALL para_env%sum(efield0)
1165  CALL para_env%sum(efield1)
1166  CALL para_env%sum(efield2)
1167  se_nddo_mpole%efield0 = se_nddo_mpole%efield0 - efield0
1168  se_nddo_mpole%efield1 = se_nddo_mpole%efield1 - efield1
1169  se_nddo_mpole%efield2 = se_nddo_mpole%efield2 - efield2
1170  ! deallocate working arrays
1171  DEALLOCATE (efield0)
1172  DEALLOCATE (efield1)
1173  DEALLOCATE (efield2)
1174 
1175  ! Corrections for two-centers one-electron terms + nuclear terms
1176  CALL para_env%sum(enuclear)
1177  CALL para_env%sum(ecore2)
1178  energy%hartree = energy%hartree + ecore2
1179  energy%core_overlap = enuclear
1180  CALL finalize_se_taper(se_taper)
1181  CALL timestop(handle)
1182  END SUBROUTINE build_fock_matrix_coul_lrc
1183 
1184 ! **************************************************************************************************
1185 !> \brief Construction of the residual part (1/R^3) of the Coulomb long-range
1186 !> term of the Fock matrix
1187 !> The 1/R^3 correction works in real-space strictly on the zero-cell,
1188 !> in order to avoid more parameters to be provided in the input..
1189 !> \param qs_env ...
1190 !> \param ks_matrix ...
1191 !> \param matrix_p ...
1192 !> \param energy ...
1193 !> \param calculate_forces ...
1194 !> \author Teodoro Laino [tlaino] - 12.2008
1195 ! **************************************************************************************************
1196  SUBROUTINE build_fock_matrix_coul_lr_r3(qs_env, ks_matrix, matrix_p, energy, calculate_forces)
1197  TYPE(qs_environment_type), POINTER :: qs_env
1198  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_matrix, matrix_p
1199  TYPE(qs_energy_type), POINTER :: energy
1200  LOGICAL, INTENT(in) :: calculate_forces
1201 
1202  CHARACTER(len=*), PARAMETER :: routinen = 'build_fock_matrix_coul_lr_r3'
1203 
1204  INTEGER :: atom_a, atom_b, ewald_type, handle, iatom, ikind, inode, ispin, itype, jatom, &
1205  jkind, natoms, natorb_a, natorb_a2, natorb_b, natorb_b2, nkind, nlocal_particles, nspins
1206  INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, se_natorb
1207  LOGICAL :: anag, atener, check, defined, found, &
1208  switch
1209  LOGICAL, ALLOCATABLE, DIMENSION(:) :: se_defined
1210  REAL(kind=dp) :: dr1, ecore2, r2inv, r3inv, rinv
1211  REAL(kind=dp), DIMENSION(2) :: ecab
1212  REAL(kind=dp), DIMENSION(2025) :: pa_a, pa_b, pb_a, pb_b
1213  REAL(kind=dp), DIMENSION(3) :: dr3inv, force_ab, rij
1214  REAL(kind=dp), DIMENSION(45, 45) :: p_block_tot_a, p_block_tot_b
1215  REAL(kind=dp), DIMENSION(:, :), POINTER :: ksa_block_a, ksa_block_b, ksb_block_a, &
1216  ksb_block_b, pa_block_a, pa_block_b, &
1217  pb_block_a, pb_block_b
1218  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
1219  TYPE(atprop_type), POINTER :: atprop
1220  TYPE(cell_type), POINTER :: cell
1221  TYPE(cp_logger_type), POINTER :: logger
1222  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: diagmat_ks, diagmat_p
1223  TYPE(dft_control_type), POINTER :: dft_control
1224  TYPE(distribution_1d_type), POINTER :: local_particles
1225  TYPE(ewald_environment_type), POINTER :: ewald_env
1226  TYPE(ewald_pw_type), POINTER :: ewald_pw
1227  TYPE(fist_nonbond_env_type), POINTER :: se_nonbond_env
1228  TYPE(mp_para_env_type), POINTER :: para_env
1229  TYPE(nddo_mpole_type), POINTER :: se_nddo_mpole
1230  TYPE(neighbor_list_iterator_p_type), &
1231  DIMENSION(:), POINTER :: nl_iterator
1232  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1233  POINTER :: sab_orb
1234  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1235  TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1236  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
1237  TYPE(section_vals_type), POINTER :: se_section
1238  TYPE(semi_empirical_control_type), POINTER :: se_control
1239  TYPE(semi_empirical_p_type), DIMENSION(:), POINTER :: se_kind_list
1240  TYPE(semi_empirical_type), POINTER :: se_kind_a, se_kind_b
1241 
1242  CALL timeset(routinen, handle)
1243 
1244  CALL get_qs_env(qs_env=qs_env) !sm->dbcsr
1245 
1246  NULLIFY (dft_control, cell, force, particle_set, diagmat_ks, &
1247  diagmat_p, local_particles, se_control, ewald_env, ewald_pw, &
1248  se_nddo_mpole, se_nonbond_env, se_section, sab_orb, logger, atprop)
1249 
1250  logger => cp_get_default_logger()
1251  CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, cell=cell, para_env=para_env, &
1252  atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set, particle_set=particle_set, &
1253  ewald_env=ewald_env, atprop=atprop, &
1254  local_particles=local_particles, ewald_pw=ewald_pw, se_nddo_mpole=se_nddo_mpole, &
1255  se_nonbond_env=se_nonbond_env, sab_orb=sab_orb)
1256 
1257  nlocal_particles = sum(local_particles%n_el(:))
1258  natoms = SIZE(particle_set)
1259  CALL ewald_env_get(ewald_env, ewald_type=ewald_type)
1260  SELECT CASE (ewald_type)
1261  CASE (do_ewald_ewald)
1262  ! Do Nothing
1263  CASE DEFAULT
1264  cpabort("Periodic SE implemented only for standard EWALD sums.")
1265  END SELECT
1266 
1267  ! Parameters
1268  se_section => section_vals_get_subs_vals(qs_env%input, "DFT%QS%SE")
1269  se_control => dft_control%qs_control%se_control
1270  anag = se_control%analytical_gradients
1271  atener = atprop%energy
1272 
1273  nspins = dft_control%nspins
1274  cpassert(ASSOCIATED(matrix_p))
1275  cpassert(SIZE(ks_matrix) > 0)
1276 
1277  CALL dbcsr_allocate_matrix_set(diagmat_ks, nspins)
1278  CALL dbcsr_allocate_matrix_set(diagmat_p, nspins)
1279 
1280  nkind = SIZE(atomic_kind_set)
1281  DO ispin = 1, nspins
1282  ! Allocate diagonal block matrices
1283  ALLOCATE (diagmat_p(ispin)%matrix, diagmat_ks(ispin)%matrix) !sm->dbcsr
1284  CALL dbcsr_get_block_diag(matrix_p(ispin)%matrix, diagmat_p(ispin)%matrix)
1285  CALL dbcsr_get_block_diag(ks_matrix(ispin)%matrix, diagmat_ks(ispin)%matrix)
1286  CALL dbcsr_set(diagmat_ks(ispin)%matrix, 0.0_dp)
1287  CALL dbcsr_replicate_all(diagmat_p(ispin)%matrix)
1288  CALL dbcsr_replicate_all(diagmat_ks(ispin)%matrix)
1289  END DO
1290 
1291  ! Possibly compute forces
1292  IF (calculate_forces) THEN
1293  CALL get_qs_env(qs_env=qs_env, force=force)
1294  CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
1295  END IF
1296  itype = get_se_type(dft_control%qs_control%method_id)
1297 
1298  ecore2 = 0.0_dp
1299  ! Real space part of the 1/R^3 sum
1300  ALLOCATE (se_defined(nkind), se_kind_list(nkind), se_natorb(nkind))
1301  DO ikind = 1, nkind
1302  CALL get_qs_kind(qs_kind_set(ikind), se_parameter=se_kind_a)
1303  se_kind_list(ikind)%se_param => se_kind_a
1304  CALL get_se_param(se_kind_a, defined=defined, natorb=natorb_a)
1305  se_defined(ikind) = (defined .AND. natorb_a >= 1)
1306  se_natorb(ikind) = natorb_a
1307  END DO
1308 
1309  CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
1310  DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1311  CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, inode=inode, r=rij)
1312  IF (.NOT. se_defined(ikind)) cycle
1313  IF (.NOT. se_defined(jkind)) cycle
1314  se_kind_a => se_kind_list(ikind)%se_param
1315  se_kind_b => se_kind_list(jkind)%se_param
1316  natorb_a = se_natorb(ikind)
1317  natorb_b = se_natorb(jkind)
1318  natorb_a2 = natorb_a**2
1319  natorb_b2 = natorb_b**2
1320 
1321  IF (inode == 1) THEN
1322  CALL dbcsr_get_block_p(matrix=diagmat_p(1)%matrix, &
1323  row=iatom, col=iatom, block=pa_block_a, found=found)
1324  cpassert(ASSOCIATED(pa_block_a))
1325  check = (SIZE(pa_block_a, 1) == natorb_a) .AND. (SIZE(pa_block_a, 2) == natorb_a)
1326  cpassert(check)
1327  CALL dbcsr_get_block_p(matrix=diagmat_ks(1)%matrix, &
1328  row=iatom, col=iatom, block=ksa_block_a, found=found)
1329  cpassert(ASSOCIATED(ksa_block_a))
1330  p_block_tot_a(1:natorb_a, 1:natorb_a) = 2.0_dp*pa_block_a
1331  pa_a(1:natorb_a2) = reshape(pa_block_a, (/natorb_a2/))
1332  IF (nspins == 2) THEN
1333  CALL dbcsr_get_block_p(matrix=diagmat_p(2)%matrix, &
1334  row=iatom, col=iatom, block=pa_block_b, found=found)
1335  cpassert(ASSOCIATED(pa_block_b))
1336  check = (SIZE(pa_block_b, 1) == natorb_a) .AND. (SIZE(pa_block_b, 2) == natorb_a)
1337  cpassert(check)
1338  CALL dbcsr_get_block_p(matrix=diagmat_ks(2)%matrix, &
1339  row=iatom, col=iatom, block=ksa_block_b, found=found)
1340  cpassert(ASSOCIATED(ksa_block_b))
1341  p_block_tot_a(1:natorb_a, 1:natorb_a) = pa_block_a + pa_block_b
1342  pa_b(1:natorb_a2) = reshape(pa_block_b, (/natorb_a2/))
1343  END IF
1344  END IF
1345 
1346  dr1 = dot_product(rij, rij)
1347  IF (dr1 > rij_threshold) THEN
1348  ! Determine the order of the atoms, and in case switch them..
1349  IF (iatom <= jatom) THEN
1350  switch = .false.
1351  ELSE
1352  switch = .true.
1353  END IF
1354  ! Retrieve blocks for KS and P
1355  CALL dbcsr_get_block_p(matrix=diagmat_p(1)%matrix, &
1356  row=jatom, col=jatom, block=pb_block_a, found=found)
1357  cpassert(ASSOCIATED(pb_block_a))
1358  check = (SIZE(pb_block_a, 1) == natorb_b) .AND. (SIZE(pb_block_a, 2) == natorb_b)
1359  cpassert(check)
1360  CALL dbcsr_get_block_p(matrix=diagmat_ks(1)%matrix, &
1361  row=jatom, col=jatom, block=ksb_block_a, found=found)
1362  cpassert(ASSOCIATED(ksb_block_a))
1363  p_block_tot_b(1:natorb_b, 1:natorb_b) = 2.0_dp*pb_block_a
1364  pb_a(1:natorb_b2) = reshape(pb_block_a, (/natorb_b2/))
1365  ! Handle more than one configuration
1366  IF (nspins == 2) THEN
1367  CALL dbcsr_get_block_p(matrix=diagmat_p(2)%matrix, &
1368  row=jatom, col=jatom, block=pb_block_b, found=found)
1369  cpassert(ASSOCIATED(pb_block_b))
1370  check = (SIZE(pb_block_b, 1) == natorb_b) .AND. (SIZE(pb_block_b, 2) == natorb_b)
1371  cpassert(check)
1372  CALL dbcsr_get_block_p(matrix=diagmat_ks(2)%matrix, &
1373  row=jatom, col=jatom, block=ksb_block_b, found=found)
1374  cpassert(ASSOCIATED(ksb_block_b))
1375  check = (SIZE(pb_block_a, 1) == SIZE(pb_block_b, 1)) .AND. (SIZE(pb_block_a, 2) == SIZE(pb_block_b, 2))
1376  cpassert(check)
1377  p_block_tot_b(1:natorb_b, 1:natorb_b) = pb_block_a + pb_block_b
1378  pb_b(1:natorb_b2) = reshape(pb_block_b, (/natorb_b2/))
1379  END IF
1380 
1381  SELECT CASE (dft_control%qs_control%method_id)
1382  CASE (do_method_mndo, do_method_am1, do_method_pm3, do_method_pm6, do_method_pm6fm, do_method_pdg, &
1383  do_method_rm1, do_method_mndod, do_method_pnnl)
1384 
1385  ! Pre-compute some quantities..
1386  r2inv = 1.0_dp/dr1
1387  rinv = sqrt(r2inv)
1388  r3inv = rinv**3
1389 
1390  ! Two-centers One-electron terms
1391  IF (nspins == 1) THEN
1392  ecab = 0._dp
1393  CALL fock2_1el_r3(se_kind_a, se_kind_b, ksa_block_a, ksb_block_a, pa_a, pb_a, &
1394  ecore=ecab, e1b=se_kind_a%expns3_int(jkind)%expns3%e1b, &
1395  e2a=se_kind_a%expns3_int(jkind)%expns3%e2a, rp=r3inv)
1396  ecore2 = ecore2 + ecab(1) + ecab(2)
1397  ELSE IF (nspins == 2) THEN
1398  ecab = 0._dp
1399  CALL fock2_1el_r3(se_kind_a, se_kind_b, ksa_block_a, ksb_block_a, pa_block_a, &
1400  pb_block_a, ecore=ecab, e1b=se_kind_a%expns3_int(jkind)%expns3%e1b, &
1401  e2a=se_kind_a%expns3_int(jkind)%expns3%e2a, rp=r3inv)
1402 
1403  CALL fock2_1el_r3(se_kind_a, se_kind_b, ksa_block_b, ksb_block_b, pa_b, pb_b, &
1404  ecore=ecab, e1b=se_kind_a%expns3_int(jkind)%expns3%e1b, &
1405  e2a=se_kind_a%expns3_int(jkind)%expns3%e2a, rp=r3inv)
1406  ecore2 = ecore2 + ecab(1) + ecab(2)
1407  END IF
1408  IF (atener) THEN
1409  atprop%atecoul(iatom) = atprop%atecoul(iatom) + ecab(1)
1410  atprop%atecoul(jatom) = atprop%atecoul(jatom) + ecab(2)
1411  END IF
1412  ! Coulomb Terms
1413  IF (nspins == 1) THEN
1414  CALL fock2c_r3(se_kind_a, se_kind_b, switch, p_block_tot_a, ksa_block_a, p_block_tot_b, &
1415  ksb_block_a, factor=0.5_dp, w=se_kind_a%expns3_int(jkind)%expns3%w, rp=r3inv)
1416  ELSE IF (nspins == 2) THEN
1417  CALL fock2c_r3(se_kind_a, se_kind_b, switch, p_block_tot_a, ksa_block_a, p_block_tot_b, &
1418  ksb_block_a, factor=1.0_dp, w=se_kind_a%expns3_int(jkind)%expns3%w, rp=r3inv)
1419 
1420  CALL fock2c_r3(se_kind_a, se_kind_b, switch, p_block_tot_a, ksa_block_b, p_block_tot_b, &
1421  ksb_block_b, factor=1.0_dp, w=se_kind_a%expns3_int(jkind)%expns3%w, rp=r3inv)
1422  END IF
1423 
1424  ! Compute forces if requested
1425  IF (calculate_forces) THEN
1426  dr3inv = -3.0_dp*rij*r3inv*r2inv
1427  atom_a = atom_of_kind(iatom)
1428  atom_b = atom_of_kind(jatom)
1429 
1430  force_ab = 0.0_dp
1431  ! Derivatives of the One-centre One-electron terms
1432  IF (nspins == 1) THEN
1433  CALL dfock2_1el_r3(se_kind_a, se_kind_b, dr3inv, pa_a, pb_a, force_ab, &
1434  se_kind_a%expns3_int(jkind)%expns3%e1b, se_kind_a%expns3_int(jkind)%expns3%e2a)
1435  ELSE IF (nspins == 2) THEN
1436  CALL dfock2_1el_r3(se_kind_a, se_kind_b, dr3inv, pa_block_a, pb_block_a, force_ab, &
1437  se_kind_a%expns3_int(jkind)%expns3%e1b, se_kind_a%expns3_int(jkind)%expns3%e2a)
1438 
1439  CALL dfock2_1el_r3(se_kind_a, se_kind_b, dr3inv, pa_b, pb_b, force_ab, &
1440  se_kind_a%expns3_int(jkind)%expns3%e1b, se_kind_a%expns3_int(jkind)%expns3%e2a)
1441  END IF
1442 
1443  ! Sum up force components
1444  force(ikind)%all_potential(1, atom_a) = force(ikind)%all_potential(1, atom_a) - force_ab(1)
1445  force(jkind)%all_potential(1, atom_b) = force(jkind)%all_potential(1, atom_b) + force_ab(1)
1446 
1447  force(ikind)%all_potential(2, atom_a) = force(ikind)%all_potential(2, atom_a) - force_ab(2)
1448  force(jkind)%all_potential(2, atom_b) = force(jkind)%all_potential(2, atom_b) + force_ab(2)
1449 
1450  force(ikind)%all_potential(3, atom_a) = force(ikind)%all_potential(3, atom_a) - force_ab(3)
1451  force(jkind)%all_potential(3, atom_b) = force(jkind)%all_potential(3, atom_b) + force_ab(3)
1452 
1453  ! Derivatives of the Coulomb Terms
1454  force_ab = 0.0_dp
1455  IF (nspins == 1) THEN
1456  CALL dfock2c_r3(se_kind_a, se_kind_b, switch, p_block_tot_a, p_block_tot_b, factor=0.25_dp, &
1457  w=se_kind_a%expns3_int(jkind)%expns3%w, drp=dr3inv, force=force_ab)
1458  ELSE IF (nspins == 2) THEN
1459  CALL dfock2c_r3(se_kind_a, se_kind_b, switch, p_block_tot_a, p_block_tot_b, factor=0.50_dp, &
1460  w=se_kind_a%expns3_int(jkind)%expns3%w, drp=dr3inv, force=force_ab)
1461 
1462  CALL dfock2c_r3(se_kind_a, se_kind_b, switch, p_block_tot_a, p_block_tot_b, factor=0.50_dp, &
1463  w=se_kind_a%expns3_int(jkind)%expns3%w, drp=dr3inv, force=force_ab)
1464  END IF
1465 
1466  ! Sum up force components
1467  force(ikind)%rho_elec(1, atom_a) = force(ikind)%rho_elec(1, atom_a) - force_ab(1)
1468  force(jkind)%rho_elec(1, atom_b) = force(jkind)%rho_elec(1, atom_b) + force_ab(1)
1469 
1470  force(ikind)%rho_elec(2, atom_a) = force(ikind)%rho_elec(2, atom_a) - force_ab(2)
1471  force(jkind)%rho_elec(2, atom_b) = force(jkind)%rho_elec(2, atom_b) + force_ab(2)
1472 
1473  force(ikind)%rho_elec(3, atom_a) = force(ikind)%rho_elec(3, atom_a) - force_ab(3)
1474  force(jkind)%rho_elec(3, atom_b) = force(jkind)%rho_elec(3, atom_b) + force_ab(3)
1475  END IF
1476  CASE DEFAULT
1477  cpabort("")
1478  END SELECT
1479  END IF
1480  END DO
1481  CALL neighbor_list_iterator_release(nl_iterator)
1482 
1483  DEALLOCATE (se_kind_list, se_defined, se_natorb)
1484 
1485  DO ispin = 1, nspins
1486  CALL dbcsr_sum_replicated(diagmat_ks(ispin)%matrix)
1487  CALL dbcsr_distribute(diagmat_ks(ispin)%matrix)
1488  CALL dbcsr_add(ks_matrix(ispin)%matrix, diagmat_ks(ispin)%matrix, &
1489  1.0_dp, 1.0_dp)
1490  END DO
1491  CALL dbcsr_deallocate_matrix_set(diagmat_p)
1492  CALL dbcsr_deallocate_matrix_set(diagmat_ks)
1493 
1494  ! Two-centers one-electron terms
1495  CALL para_env%sum(ecore2)
1496  energy%hartree = energy%hartree + ecore2
1497 
1498  CALL timestop(handle)
1499  END SUBROUTINE build_fock_matrix_coul_lr_r3
1500 
1501 END MODULE se_fock_matrix_coulomb
1502 
fac
static GRID_HOST_DEVICE double fac(const int i)
Factorial function, e.g. fac(5) = 5! = 120.
Definition: grid_common.h:48
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition: cp_dbcsr_operations.F:81
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition: cp_dbcsr_operations.F:89
atomic_kind_types
Define the atomic kind types and their sub types.
Definition: atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Definition: atomic_kind_types.F:223
atprop_types
Holds information on atomic properties.
Definition: atprop_types.F:14
cell_types
Handles all functions related to the CELL.
Definition: cell_types.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition: cp_control_types.F:12
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition: cp_dbcsr_operations.F:18
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition: cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition: cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition: cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition: cp_output_handling.F:803
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition: cp_output_handling.F:1013
distribution_1d_types
stores a lists of integer that are local to a processor. The idea is that these integers represent ob...
Definition: distribution_1d_types.F:24
ewald_environment_types
Definition: ewald_environment_types.F:14
ewald_environment_types::ewald_env_get
subroutine, public ewald_env_get(ewald_env, ewald_type, alpha, eps_pol, epsilon, gmax, ns_max, o_spline, group, para_env, poisson_section, precs, rcut, do_multipoles, max_multipole, do_ipol, max_ipol_iter, interaction_cutoffs, cell_hmat)
Purpose: Get the EWALD environment.
Definition: ewald_environment_types.F:123
ewald_pw_types
pw_types
Definition: ewald_pw_types.F:12
ewald_pw_types::ewald_pw_get
subroutine, public ewald_pw_get(ewald_pw, pw_big_pool, pw_small_pool, rs_desc, poisson_env, dg)
get the ewald_pw environment to the correct program.
Definition: ewald_pw_types.F:325
ewalds_multipole
Treats the electrostatic for multipoles (up to quadrupoles)
Definition: ewalds_multipole.F:14
ewalds_multipole::ewald_multipole_evaluate
recursive subroutine, public ewald_multipole_evaluate(ewald_env, ewald_pw, nonbond_env, cell, particle_set, local_particles, energy_local, energy_glob, e_neut, e_self, task, do_correction_bonded, do_forces, do_stress, do_efield, radii, charges, dipoles, quadrupoles, forces_local, forces_glob, pv_local, pv_glob, efield0, efield1, efield2, iw, do_debug, atomic_kind_set, mm_section)
Computes the potential and the force for a lattice sum of multipoles (up to quadrupole)
Definition: ewalds_multipole.F:122
fist_neighbor_list_control
Definition: fist_neighbor_list_control.F:15
fist_neighbor_list_control::list_control
subroutine, public list_control(atomic_kind_set, particle_set, local_particles, cell, fist_nonbond_env, para_env, mm_section, shell_particle_set, core_particle_set, force_update, exclusions)
...
Definition: fist_neighbor_list_control.F:86
fist_nonbond_env_types
Definition: fist_nonbond_env_types.F:13
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition: input_constants.F:17
input_constants::do_method_pdg
integer, parameter, public do_method_pdg
Definition: input_constants.F:250
input_constants::do_method_pnnl
integer, parameter, public do_method_pnnl
Definition: input_constants.F:250
input_constants::do_method_rm1
integer, parameter, public do_method_rm1
Definition: input_constants.F:250
input_constants::do_method_pm3
integer, parameter, public do_method_pm3
Definition: input_constants.F:250
input_constants::do_method_mndo
integer, parameter, public do_method_mndo
Definition: input_constants.F:250
input_constants::do_method_mndod
integer, parameter, public do_method_mndod
Definition: input_constants.F:250
input_constants::do_method_am1
integer, parameter, public do_method_am1
Definition: input_constants.F:250
input_constants::do_se_is_slater
integer, parameter, public do_se_is_slater
Definition: input_constants.F:277
input_constants::do_method_pm6fm
integer, parameter, public do_method_pm6fm
Definition: input_constants.F:250
input_constants::do_method_pm6
integer, parameter, public do_method_pm6
Definition: input_constants.F:250
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition: input_section_types.F:15
input_section_types::section_vals_get_subs_vals
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
Definition: input_section_types.F:724
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
message_passing
Interface to the message passing library MPI.
Definition: message_passing.F:23
multipole_types
Multipole structure: for multipole (fixed and induced) in FF based MD.
Definition: multipole_types.F:12
multipole_types::do_multipole_quadrupole
integer, parameter, public do_multipole_quadrupole
Definition: multipole_types.F:31
multipole_types::do_multipole_dipole
integer, parameter, public do_multipole_dipole
Definition: multipole_types.F:31
multipole_types::do_multipole_charge
integer, parameter, public do_multipole_charge
Definition: multipole_types.F:31
multipole_types::do_multipole_none
integer, parameter, public do_multipole_none
Definition: multipole_types.F:31
particle_types
Define the data structure for the particle information.
Definition: particle_types.F:19
pw_poisson_types
functions related to the poisson solver on regular grids
Definition: pw_poisson_types.F:22
pw_poisson_types::do_ewald_ewald
integer, parameter, public do_ewald_ewald
Definition: pw_poisson_types.F:75
qs_energy_types
Definition: qs_energy_types.F:14
qs_environment_types
Definition: qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_RI_aux_kp, matrix_s, matrix_s_RI_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, ecoul_1c, rho0_s_rs, rho0_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, WannierCentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, rhs)
Get the QUICKSTEP environment.
Definition: qs_environment_types.F:502
qs_force_types
Definition: qs_force_types.F:13
qs_kind_types
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, U_of_dft_plus_u, J_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, J0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition: qs_kind_types.F:451
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition: qs_neighbor_list_types.F:17
qs_neighbor_list_types::neighbor_list_iterator_create
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
Definition: qs_neighbor_list_types.F:165
qs_neighbor_list_types::neighbor_list_iterator_release
subroutine, public neighbor_list_iterator_release(iterator_set)
...
Definition: qs_neighbor_list_types.F:251
qs_neighbor_list_types::neighbor_list_iterate
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
Definition: qs_neighbor_list_types.F:351
qs_neighbor_list_types::get_iterator_info
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
Definition: qs_neighbor_list_types.F:549
se_fock_matrix_coulomb
Module that collects all Coulomb parts of the fock matrix construction.
Definition: se_fock_matrix_coulomb.F:19
se_fock_matrix_coulomb::build_fock_matrix_coulomb
subroutine, public build_fock_matrix_coulomb(qs_env, ks_matrix, matrix_p, energy, calculate_forces, store_int_env)
Construction of the Coulomb part of the Fock matrix.
Definition: se_fock_matrix_coulomb.F:120
se_fock_matrix_coulomb::build_fock_matrix_coul_lr_r3
subroutine, public build_fock_matrix_coul_lr_r3(qs_env, ks_matrix, matrix_p, energy, calculate_forces)
Construction of the residual part (1/R^3) of the Coulomb long-range term of the Fock matrix The 1/R^3...
Definition: se_fock_matrix_coulomb.F:1197
se_fock_matrix_coulomb::build_fock_matrix_coulomb_lr
subroutine, public build_fock_matrix_coulomb_lr(qs_env, ks_matrix, matrix_p, energy, calculate_forces, store_int_env)
Long-Range part for SE Coulomb interactions.
Definition: se_fock_matrix_coulomb.F:476
se_fock_matrix_integrals
Provides the low level routines to build both the exchange and the Coulomb Fock matrices....
Definition: se_fock_matrix_integrals.F:18
se_fock_matrix_integrals::fock2c_r3
subroutine, public fock2c_r3(sepi, sepj, switch, pi_tot, fi_mat, pj_tot, fj_mat, factor, w, rp)
Construction of 2-center Fock Matrix - Coulomb Terms for the residual (1/R^3) integral part.
Definition: se_fock_matrix_integrals.F:889
se_fock_matrix_integrals::fock2_1el_ew
subroutine, public fock2_1el_ew(sep, rij, ks_block, p_block, ecore, itype, anag, se_int_control, se_taper, store_int_env)
Construction of 2-center 1-electron Fock Matrix (Ewald self term)
Definition: se_fock_matrix_integrals.F:289
se_fock_matrix_integrals::dfock2c
subroutine, public dfock2c(sepi, sepj, rij, switch, pi_tot, pj_tot, factor, anag, se_int_control, se_taper, force, delta)
Derivatives of 2-center Fock Matrix - Coulomb Terms.
Definition: se_fock_matrix_integrals.F:507
se_fock_matrix_integrals::fock2_1el_r3
subroutine, public fock2_1el_r3(sepi, sepj, ksi_block, ksj_block, pi_block, pj_block, e1b, e2a, ecore, rp)
Construction of 2-center 1-electron Fock Matrix for the residual (1/R^3) integral part.
Definition: se_fock_matrix_integrals.F:780
se_fock_matrix_integrals::dfock2c_r3
subroutine, public dfock2c_r3(sepi, sepj, switch, pi_tot, pj_tot, factor, w, drp, force)
Derivatives of 2-center Fock Matrix - Coulomb Terms for the residual (1/R^3) integral part.
Definition: se_fock_matrix_integrals.F:958
se_fock_matrix_integrals::fock2c
subroutine, public fock2c(sepi, sepj, rij, switch, pi_tot, fi_mat, pj_tot, fj_mat, factor, anag, se_int_control, se_taper, store_int_env)
Construction of 2-center Fock Matrix - Coulomb Terms.
Definition: se_fock_matrix_integrals.F:417
se_fock_matrix_integrals::se_coulomb_ij_interaction
subroutine, public se_coulomb_ij_interaction(atom_a, atom_b, my_task, do_forces, do_efield, do_stress, charges, dipoles, quadrupoles, force_ab, efield0, efield1, efield2, rab2, rab, integral_value, ptens11, ptens12, ptens13, ptens21, ptens22, ptens23, ptens31, ptens32, ptens33)
Coulomb interaction multipolar correction.
Definition: se_fock_matrix_integrals.F:1046
se_fock_matrix_integrals::dfock2_1el_r3
subroutine, public dfock2_1el_r3(sepi, sepj, drp, pi_block, pj_block, force, e1b, e2a)
Derivatives of 2-center 1-electron Fock Matrix residual (1/R^3) integral part.
Definition: se_fock_matrix_integrals.F:831
se_fock_matrix_integrals::fock2_1el
subroutine, public fock2_1el(sepi, sepj, rij, ksi_block, ksj_block, pi_block, pj_block, ecore, itype, anag, se_int_control, se_taper, store_int_env)
Construction of 2-center 1-electron Fock Matrix.
Definition: se_fock_matrix_integrals.F:65
se_fock_matrix_integrals::dfock2_1el
subroutine, public dfock2_1el(sepi, sepj, rij, pi_block, pj_block, itype, anag, se_int_control, se_taper, force, delta)
Derivatives of 2-center 1-electron Fock Matrix.
Definition: se_fock_matrix_integrals.F:143
se_fock_matrix_integrals::fock2c_ew
subroutine, public fock2c_ew(sep, rij, p_tot, f_mat, factor, anag, se_int_control, se_taper, store_int_env)
Construction of 2-center Fock Matrix - Coulomb Self Terms (Ewald)
Definition: se_fock_matrix_integrals.F:347
semi_empirical_int_arrays
Arrays of parameters used in the semi-empirical calculations \References Everywhere in this module TC...
Definition: semi_empirical_int_arrays.F:17
semi_empirical_int_arrays::rij_threshold
real(kind=dp), parameter, public rij_threshold
Definition: semi_empirical_int_arrays.F:27
semi_empirical_int_arrays::se_orbital_pointer
integer, dimension(9), public se_orbital_pointer
Definition: semi_empirical_int_arrays.F:30
semi_empirical_integrals
Set of wrappers for semi-empirical analytical/numerical Integrals routines.
Definition: semi_empirical_integrals.F:16
semi_empirical_integrals::dcorecore_el
subroutine, public dcorecore_el(sepi, sepj, rij, denuc, itype, delta, anag, se_int_control, se_taper)
wrapper for numerical/analytical routines core-core electrostatic (only) integrals derivatives
Definition: semi_empirical_integrals.F:589
semi_empirical_integrals::corecore_el
subroutine, public corecore_el(sepi, sepj, rij, enuc, itype, anag, se_int_control, se_taper)
wrapper for numerical/analytical routines core-core electrostatic (only) integrals
Definition: semi_empirical_integrals.F:437
semi_empirical_mpole_methods
Setup and Methods for semi-empirical multipole types.
Definition: semi_empirical_mpole_methods.F:12
semi_empirical_mpole_methods::quadrupole_sph_to_cart
subroutine, public quadrupole_sph_to_cart(qcart, qsph)
Transforms the quadrupole components from sphericals to cartesians.
Definition: semi_empirical_mpole_methods.F:237
semi_empirical_mpole_types
Definition of the semi empirical multipole integral expansions types.
Definition: semi_empirical_mpole_types.F:12
semi_empirical_store_int_types
Type to store integrals for semi-empirical calculations.
Definition: semi_empirical_store_int_types.F:13
semi_empirical_types
Definition of the semi empirical parameter types.
Definition: semi_empirical_types.F:12
semi_empirical_types::setup_se_int_control_type
subroutine, public setup_se_int_control_type(se_int_control, shortrange, do_ewald_r3, do_ewald_gks, integral_screening, max_multipole, pc_coulomb_int)
Setup the Semiempirical integral control type.
Definition: semi_empirical_types.F:555
semi_empirical_types::get_se_param
subroutine, public get_se_param(sep, name, typ, defined, z, zeff, natorb, eheat, beta, sto_exponents, uss, upp, udd, uff, alp, eisol, gss, gsp, gpp, gp2, acoul, nr, de, ass, asp, app, hsp, gsd, gpd, gdd, ppddg, dpddg, ngauss)
Get info from the semi-empirical type.
Definition: semi_empirical_types.F:363
semi_empirical_utils
Working with the semi empirical parameter types.
Definition: semi_empirical_utils.F:12
semi_empirical_utils::finalize_se_taper
subroutine, public finalize_se_taper(se_taper)
Finalizes the semi-empirical taper for a chunk calculation.
Definition: semi_empirical_utils.F:199
semi_empirical_utils::get_se_type
integer function, public get_se_type(se_method)
Gives back the unique semi_empirical METHOD type.
Definition: semi_empirical_utils.F:336
semi_empirical_utils::initialize_se_taper
subroutine, public initialize_se_taper(se_taper, coulomb, exchange, lr_corr)
Initializes the semi-empirical taper for a chunk calculation.
Definition: semi_empirical_utils.F:163
virial_methods
Definition: virial_methods.F:12
virial_methods::virial_pair_force
pure subroutine, public virial_pair_force(pv_virial, f0, force, rab)
Computes the contribution to the stress tensor from two-body pair-wise forces.
Definition: virial_methods.F:139
virial_types
Definition: virial_types.F:13