(git:374b731)
Loading...
Searching...
No Matches
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! **************************************************************************************************
19MODULE 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
106CONTAINS
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
1501END 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
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
atprop_types::atprop_type
type for the atomic properties
Definition atprop_types.F:31
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:55
cp_control_types::dft_control_type
Definition cp_control_types.F:559
cp_control_types::semi_empirical_control_type
Definition cp_control_types.F:173
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
distribution_1d_types::distribution_1d_type
structure to store local (to a processor) ordered lists of integers.
Definition distribution_1d_types.F:62
ewald_environment_types::ewald_environment_type
to build arrays of pointers
Definition ewald_environment_types.F:55
ewald_pw_types::ewald_pw_type
Definition ewald_pw_types.F:64
fist_nonbond_env_types::fist_nonbond_env_type
Definition fist_nonbond_env_types.F:81
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:126
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:763
particle_types::particle_type
Definition particle_types.F:35
qs_energy_types::qs_energy_type
Definition qs_energy_types.F:25
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:215
qs_force_types::qs_force_type
Definition qs_force_types.F:28
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:171
qs_neighbor_list_types::neighbor_list_iterator_p_type
Definition qs_neighbor_list_types.F:117
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
semi_empirical_mpole_types::nddo_mpole_type
Global Multipolar NDDO information type.
Definition semi_empirical_mpole_types.F:54
semi_empirical_mpole_types::semi_empirical_mpole_type
Semi-empirical integral multipole expansion type.
Definition semi_empirical_mpole_types.F:29
semi_empirical_store_int_types::semi_empirical_si_type
Semi-empirical store integrals type.
Definition semi_empirical_store_int_types.F:43
semi_empirical_types::se_int_control_type
Definition semi_empirical_types.F:136
semi_empirical_types::se_taper_type
Taper type use in semi-empirical calculations.
Definition semi_empirical_types.F:157
semi_empirical_types::semi_empirical_p_type
Definition semi_empirical_types.F:103
semi_empirical_types::semi_empirical_type
Semi-empirical type.
Definition semi_empirical_types.F:57
virial_types::virial_type
Definition virial_types.F:26