Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief 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_api, ONLY: dbcsr_add,&
27 : dbcsr_distribute,&
28 : dbcsr_get_block_diag,&
29 : dbcsr_get_block_p,&
30 : dbcsr_p_type,&
31 : dbcsr_replicate_all,&
32 : dbcsr_set,&
33 : dbcsr_sum_replicated
34 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set,&
35 : dbcsr_deallocate_matrix_set
36 : USE cp_log_handling, ONLY: cp_get_default_logger,&
37 : cp_logger_type
38 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
39 : cp_print_key_unit_nr
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 39152 : 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 39152 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, se_natorb
132 : LOGICAL :: anag, atener, check, defined, found, &
133 : switch, use_virial
134 39152 : 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 39152 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: ksa_block_a, ksa_block_b, ksb_block_a, &
141 39152 : ksb_block_b, pa_block_a, pa_block_b, &
142 39152 : pb_block_a, pb_block_b
143 39152 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
144 : TYPE(atprop_type), POINTER :: atprop
145 : TYPE(cell_type), POINTER :: cell
146 39152 : 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 39152 : DIMENSION(:), POINTER :: nl_iterator
153 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
154 39152 : POINTER :: sab_se
155 39152 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
156 39152 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
157 39152 : 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 39152 : 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 39152 : CALL timeset(routineN, handle)
166 :
167 39152 : NULLIFY (dft_control, cell, force, particle_set, diagmat_ks, diagmat_p, &
168 39152 : 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 39152 : qs_kind_set=qs_kind_set, particle_set=particle_set, virial=virial)
173 :
174 : ! Parameters
175 39152 : CALL initialize_se_taper(se_taper, coulomb=.TRUE.)
176 39152 : se_control => dft_control%qs_control%se_control
177 39152 : anag = se_control%analytical_gradients
178 39152 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
179 39152 : 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 76568 : max_multipole=se_control%max_multipole, pc_coulomb_int=.FALSE.)
185 39152 : IF (se_control%do_ewald_gks) THEN
186 106 : CALL get_qs_env(qs_env=qs_env, ewald_env=ewald_env, ewald_pw=ewald_pw)
187 106 : 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 106 : dg=se_int_control%ewald_gks%dg)
190 : END IF
191 :
192 39152 : nspins = dft_control%nspins
193 39152 : CPASSERT(ASSOCIATED(matrix_p))
194 39152 : CPASSERT(SIZE(ks_matrix) > 0)
195 :
196 39152 : nkind = SIZE(atomic_kind_set)
197 39152 : IF (calculate_forces) THEN
198 3002 : CALL get_qs_env(qs_env=qs_env, force=force)
199 3002 : delta = se_control%delta
200 3002 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
201 : END IF
202 :
203 39152 : CALL dbcsr_allocate_matrix_set(diagmat_ks, nspins)
204 39152 : CALL dbcsr_allocate_matrix_set(diagmat_p, nspins)
205 :
206 80868 : DO ispin = 1, nspins
207 : ! Allocate diagonal block matrices
208 41716 : ALLOCATE (diagmat_p(ispin)%matrix, diagmat_ks(ispin)%matrix) !sm->dbcsr
209 41716 : CALL dbcsr_get_block_diag(matrix_p(ispin)%matrix, diagmat_p(ispin)%matrix)
210 41716 : CALL dbcsr_get_block_diag(ks_matrix(ispin)%matrix, diagmat_ks(ispin)%matrix)
211 41716 : CALL dbcsr_set(diagmat_ks(ispin)%matrix, 0.0_dp)
212 41716 : CALL dbcsr_replicate_all(diagmat_p(ispin)%matrix)
213 80868 : CALL dbcsr_replicate_all(diagmat_ks(ispin)%matrix)
214 : END DO
215 :
216 39152 : ecore2 = 0.0_dp
217 39152 : itype = get_se_type(dft_control%qs_control%method_id)
218 331394 : ALLOCATE (se_defined(nkind), se_kind_list(nkind), se_natorb(nkind))
219 135634 : DO ikind = 1, nkind
220 96482 : CALL get_qs_kind(qs_kind_set(ikind), se_parameter=se_kind_a)
221 96482 : se_kind_list(ikind)%se_param => se_kind_a
222 96482 : CALL get_se_param(se_kind_a, defined=defined, natorb=natorb_a)
223 96482 : se_defined(ikind) = (defined .AND. natorb_a >= 1)
224 232116 : se_natorb(ikind) = natorb_a
225 : END DO
226 :
227 39152 : CALL neighbor_list_iterator_create(nl_iterator, sab_se)
228 9898746 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
229 9859594 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, inode=inode, r=rij)
230 9859594 : IF (.NOT. se_defined(ikind)) CYCLE
231 9859594 : IF (.NOT. se_defined(jkind)) CYCLE
232 9859594 : se_kind_a => se_kind_list(ikind)%se_param
233 9859594 : se_kind_b => se_kind_list(jkind)%se_param
234 9859594 : natorb_a = se_natorb(ikind)
235 9859594 : natorb_b = se_natorb(jkind)
236 9859594 : natorb_a2 = natorb_a**2
237 9859594 : natorb_b2 = natorb_b**2
238 :
239 9859594 : IF (inode == 1) THEN
240 : CALL dbcsr_get_block_p(matrix=diagmat_p(1)%matrix, &
241 620889 : row=iatom, col=iatom, BLOCK=pa_block_a, found=found)
242 620889 : CPASSERT(ASSOCIATED(pa_block_a))
243 620889 : check = (SIZE(pa_block_a, 1) == natorb_a) .AND. (SIZE(pa_block_a, 2) == natorb_a)
244 0 : CPASSERT(check)
245 : CALL dbcsr_get_block_p(matrix=diagmat_ks(1)%matrix, &
246 620889 : row=iatom, col=iatom, BLOCK=ksa_block_a, found=found)
247 620889 : CPASSERT(ASSOCIATED(ksa_block_a))
248 9733845 : p_block_tot_a(1:natorb_a, 1:natorb_a) = 2.0_dp*pa_block_a
249 1241778 : pa_a(1:natorb_a2) = RESHAPE(pa_block_a, (/natorb_a2/))
250 620889 : IF (nspins == 2) THEN
251 : CALL dbcsr_get_block_p(matrix=diagmat_p(2)%matrix, &
252 8389 : row=iatom, col=iatom, BLOCK=pa_block_b, found=found)
253 8389 : CPASSERT(ASSOCIATED(pa_block_b))
254 8389 : check = (SIZE(pa_block_b, 1) == natorb_a) .AND. (SIZE(pa_block_b, 2) == natorb_a)
255 0 : CPASSERT(check)
256 : CALL dbcsr_get_block_p(matrix=diagmat_ks(2)%matrix, &
257 8389 : row=iatom, col=iatom, BLOCK=ksa_block_b, found=found)
258 8389 : CPASSERT(ASSOCIATED(ksa_block_b))
259 117595 : p_block_tot_a(1:natorb_a, 1:natorb_a) = pa_block_a + pa_block_b
260 16778 : pa_b(1:natorb_a2) = RESHAPE(pa_block_b, (/natorb_a2/))
261 : END IF
262 :
263 : END IF
264 :
265 39438376 : dr1 = DOT_PRODUCT(rij, rij)
266 9898746 : IF (dr1 > rij_threshold) THEN
267 : ! Determine the order of the atoms, and in case switch them..
268 9670609 : IF (iatom <= jatom) THEN
269 5083252 : switch = .FALSE.
270 : ELSE
271 4587357 : switch = .TRUE.
272 : END IF
273 : ! Retrieve blocks for KS and P
274 : CALL dbcsr_get_block_p(matrix=diagmat_p(1)%matrix, &
275 9670609 : row=jatom, col=jatom, BLOCK=pb_block_a, found=found)
276 9670609 : CPASSERT(ASSOCIATED(pb_block_a))
277 9670609 : check = (SIZE(pb_block_a, 1) == natorb_b) .AND. (SIZE(pb_block_a, 2) == natorb_b)
278 0 : CPASSERT(check)
279 : CALL dbcsr_get_block_p(matrix=diagmat_ks(1)%matrix, &
280 9670609 : row=jatom, col=jatom, BLOCK=ksb_block_a, found=found)
281 9670609 : CPASSERT(ASSOCIATED(ksb_block_a))
282 147643123 : p_block_tot_b(1:natorb_b, 1:natorb_b) = 2.0_dp*pb_block_a
283 19341218 : pb_a(1:natorb_b2) = RESHAPE(pb_block_a, (/natorb_b2/))
284 : ! Handle more than one configuration
285 9670609 : IF (nspins == 2) THEN
286 : CALL dbcsr_get_block_p(matrix=diagmat_p(2)%matrix, &
287 225262 : row=jatom, col=jatom, BLOCK=pb_block_b, found=found)
288 225262 : CPASSERT(ASSOCIATED(pb_block_b))
289 225262 : check = (SIZE(pb_block_b, 1) == natorb_b) .AND. (SIZE(pb_block_b, 2) == natorb_b)
290 0 : CPASSERT(check)
291 : CALL dbcsr_get_block_p(matrix=diagmat_ks(2)%matrix, &
292 225262 : row=jatom, col=jatom, BLOCK=ksb_block_b, found=found)
293 225262 : CPASSERT(ASSOCIATED(ksb_block_b))
294 225262 : check = (SIZE(pb_block_a, 1) == SIZE(pb_block_b, 1)) .AND. (SIZE(pb_block_a, 2) == SIZE(pb_block_b, 2))
295 0 : CPASSERT(check)
296 2764564 : p_block_tot_b(1:natorb_b, 1:natorb_b) = pb_block_a + pb_block_b
297 450524 : pb_b(1:natorb_b2) = RESHAPE(pb_block_b, (/natorb_b2/))
298 : END IF
299 :
300 19341218 : 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 9670609 : IF (nspins == 1) THEN
306 9445347 : 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 9445347 : se_taper=se_taper, store_int_env=store_int_env)
310 9445347 : ecore2 = ecore2 + ecab(1) + ecab(2)
311 225262 : ELSE IF (nspins == 2) THEN
312 225262 : 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 225262 : 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 225262 : se_taper=se_taper, store_int_env=store_int_env)
319 225262 : ecore2 = ecore2 + ecab(1) + ecab(2)
320 : END IF
321 9670609 : IF (atener) THEN
322 765 : atprop%atecoul(iatom) = atprop%atecoul(iatom) + ecab(1)
323 765 : atprop%atecoul(jatom) = atprop%atecoul(jatom) + ecab(2)
324 : END IF
325 : ! Coulomb Terms
326 9670609 : 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 9445347 : store_int_env=store_int_env)
330 225262 : 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 225262 : 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 225262 : store_int_env=store_int_env)
338 : END IF
339 :
340 9670609 : IF (calculate_forces) THEN
341 302154 : atom_a = atom_of_kind(iatom)
342 302154 : atom_b = atom_of_kind(jatom)
343 :
344 : ! Derivatives of the Two-centre One-electron terms
345 302154 : force_ab = 0.0_dp
346 302154 : 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 301993 : delta=delta)
350 161 : 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 161 : 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 161 : se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta)
355 : END IF
356 302154 : IF (use_virial) THEN
357 18728 : CALL virial_pair_force(virial%pv_virial, -1.0_dp, force_ab, rij)
358 : END IF
359 :
360 : ! Sum up force components
361 302154 : force(ikind)%all_potential(1, atom_a) = force(ikind)%all_potential(1, atom_a) - force_ab(1)
362 302154 : force(jkind)%all_potential(1, atom_b) = force(jkind)%all_potential(1, atom_b) + force_ab(1)
363 :
364 302154 : force(ikind)%all_potential(2, atom_a) = force(ikind)%all_potential(2, atom_a) - force_ab(2)
365 302154 : force(jkind)%all_potential(2, atom_b) = force(jkind)%all_potential(2, atom_b) + force_ab(2)
366 :
367 302154 : force(ikind)%all_potential(3, atom_a) = force(ikind)%all_potential(3, atom_a) - force_ab(3)
368 302154 : 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 302154 : force_ab = 0._dp
372 302154 : 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 301993 : anag=anag, se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta)
375 161 : 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 161 : 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 161 : anag=anag, se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta)
381 : END IF
382 302154 : IF (switch) THEN
383 149436 : force_ab(1) = -force_ab(1)
384 149436 : force_ab(2) = -force_ab(2)
385 149436 : force_ab(3) = -force_ab(3)
386 : END IF
387 302154 : IF (use_virial) THEN
388 18728 : CALL virial_pair_force(virial%pv_virial, -1.0_dp, force_ab, rij)
389 : END IF
390 : ! Sum up force components
391 302154 : force(ikind)%rho_elec(1, atom_a) = force(ikind)%rho_elec(1, atom_a) - force_ab(1)
392 302154 : force(jkind)%rho_elec(1, atom_b) = force(jkind)%rho_elec(1, atom_b) + force_ab(1)
393 :
394 302154 : force(ikind)%rho_elec(2, atom_a) = force(ikind)%rho_elec(2, atom_a) - force_ab(2)
395 302154 : force(jkind)%rho_elec(2, atom_b) = force(jkind)%rho_elec(2, atom_b) + force_ab(2)
396 :
397 302154 : force(ikind)%rho_elec(3, atom_a) = force(ikind)%rho_elec(3, atom_a) - force_ab(3)
398 302154 : force(jkind)%rho_elec(3, atom_b) = force(jkind)%rho_elec(3, atom_b) + force_ab(3)
399 : END IF
400 : CASE DEFAULT
401 9670609 : CPABORT("")
402 : END SELECT
403 : ELSE
404 188985 : IF (se_int_control%do_ewald_gks) THEN
405 159 : CPASSERT(iatom == jatom)
406 : ! Two-centers One-electron terms
407 159 : ecores = 0._dp
408 159 : 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 159 : se_taper=se_taper, store_int_env=store_int_env)
412 0 : 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 0 : 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 0 : se_taper=se_taper, store_int_env=store_int_env)
419 : END IF
420 159 : ecore2 = ecore2 + ecores
421 159 : IF (atener) THEN
422 0 : atprop%atecoul(iatom) = atprop%atecoul(iatom) + ecores
423 : END IF
424 : ! Coulomb Terms
425 159 : 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 159 : store_int_env=store_int_env)
429 0 : 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 0 : 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 0 : store_int_env=store_int_env)
436 : END IF
437 : END IF
438 : END IF
439 : END DO
440 39152 : CALL neighbor_list_iterator_release(nl_iterator)
441 :
442 39152 : DEALLOCATE (se_kind_list, se_defined, se_natorb)
443 :
444 80868 : DO ispin = 1, nspins
445 41716 : CALL dbcsr_sum_replicated(diagmat_ks(ispin)%matrix)
446 41716 : CALL dbcsr_distribute(diagmat_ks(ispin)%matrix)
447 : CALL dbcsr_add(ks_matrix(ispin)%matrix, diagmat_ks(ispin)%matrix, &
448 80868 : 1.0_dp, 1.0_dp)
449 : END DO
450 39152 : CALL dbcsr_deallocate_matrix_set(diagmat_p)
451 39152 : CALL dbcsr_deallocate_matrix_set(diagmat_ks)
452 :
453 : ! Two-centers one-electron terms
454 39152 : CALL para_env%sum(ecore2)
455 39152 : energy%hartree = ecore2 - energy%core
456 : ! WRITE ( *, * ) 'IN SE_F_COUL', ecore2, energy%core
457 :
458 39152 : CALL finalize_se_taper(se_taper)
459 :
460 39152 : CALL timestop(handle)
461 78304 : 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 1630 : 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 1630 : 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 1630 : 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 1630 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: ksa_block_a, pa_block_a
498 1630 : 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 1630 : 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 1630 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
511 1630 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
512 1630 : 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 1630 : CALL timeset(routineN, handle)
520 :
521 1630 : NULLIFY (dft_control, cell, force, particle_set, diagmat_ks, diagmat_p, local_particles, &
522 1630 : se_control, ewald_env, ewald_pw, se_nddo_mpole, se_nonbond_env, se_section, mpole, &
523 1630 : logger, virial, atprop)
524 :
525 1630 : 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 1630 : se_nonbond_env=se_nonbond_env, virial=virial, atprop=atprop)
531 :
532 4538 : nlocal_particles = SUM(local_particles%n_el(:))
533 1630 : natoms = SIZE(particle_set)
534 1630 : CALL ewald_env_get(ewald_env, ewald_type=ewald_type)
535 : SELECT CASE (ewald_type)
536 : CASE (do_ewald_ewald)
537 0 : forces_g_size = nlocal_particles
538 : CASE DEFAULT
539 1630 : CPABORT("Periodic SE implemented only for standard EWALD sums.")
540 : END SELECT
541 :
542 : ! Parameters
543 1630 : se_section => section_vals_get_subs_vals(qs_env%input, "DFT%QS%SE")
544 1630 : se_control => dft_control%qs_control%se_control
545 1630 : anag = se_control%analytical_gradients
546 1630 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) .AND. calculate_forces
547 1630 : atener = atprop%energy
548 :
549 1630 : nspins = dft_control%nspins
550 1630 : CPASSERT(ASSOCIATED(matrix_p))
551 1630 : CPASSERT(SIZE(ks_matrix) > 0)
552 :
553 1630 : CALL dbcsr_allocate_matrix_set(diagmat_ks, nspins)
554 1630 : CALL dbcsr_allocate_matrix_set(diagmat_p, nspins)
555 :
556 1630 : nkind = SIZE(atomic_kind_set)
557 3260 : DO ispin = 1, nspins
558 : ! Allocate diagonal block matrices
559 1630 : ALLOCATE (diagmat_p(ispin)%matrix, diagmat_ks(ispin)%matrix) !sm->dbcsr
560 1630 : CALL dbcsr_get_block_diag(matrix_p(ispin)%matrix, diagmat_p(ispin)%matrix)
561 1630 : CALL dbcsr_get_block_diag(ks_matrix(ispin)%matrix, diagmat_ks(ispin)%matrix)
562 1630 : CALL dbcsr_set(diagmat_ks(ispin)%matrix, 0.0_dp)
563 1630 : CALL dbcsr_replicate_all(diagmat_p(ispin)%matrix)
564 3260 : CALL dbcsr_replicate_all(diagmat_ks(ispin)%matrix)
565 : END DO
566 :
567 : ! Check for implemented SE methods
568 1630 : 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 0 : itype = get_se_type(dft_control%qs_control%method_id)
572 : CASE DEFAULT
573 1630 : CPABORT("")
574 : END SELECT
575 :
576 : ! Zero arrays and possibly build neighbor lists
577 1630 : energy_local = 0.0_dp
578 1630 : energy_glob = 0.0_dp
579 1630 : e_neut = 0.0_dp
580 1630 : e_self = 0.0_dp
581 1630 : task = .FALSE.
582 1630 : SELECT CASE (se_control%max_multipole)
583 : CASE (do_multipole_none)
584 : ! Do Nothing
585 : CASE (do_multipole_charge)
586 0 : task(1) = .TRUE.
587 : CASE (do_multipole_dipole)
588 0 : task = .TRUE.
589 0 : task(3) = .FALSE.
590 : CASE (do_multipole_quadrupole)
591 4832 : task = .TRUE.
592 : CASE DEFAULT
593 1630 : 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 1630 : cell, se_nonbond_env, para_env, se_section)
599 :
600 1630 : enuc = 0.0_dp
601 1630 : energy%core_overlap = 0.0_dp
602 16610 : se_nddo_mpole%charge = 0.0_dp
603 61550 : se_nddo_mpole%dipole = 0.0_dp
604 196370 : se_nddo_mpole%quadrupole = 0.0_dp
605 :
606 3260 : DO ispin = 1, nspins
607 : ! Compute the NDDO mpole expansion
608 6168 : DO ikind = 1, nkind
609 2908 : CALL get_qs_kind(qs_kind_set(ikind), se_parameter=se_kind_a)
610 2908 : CALL get_se_param(se_kind_a, defined=defined, natorb=natorb_a)
611 :
612 2908 : IF (.NOT. defined .OR. natorb_a < 1) CYCLE
613 :
614 2908 : nparticle_local = local_particles%n_el(ikind)
615 14936 : DO ilist = 1, nparticle_local
616 7490 : iatom = local_particles%list(ikind)%array(ilist)
617 : CALL dbcsr_get_block_p(matrix=diagmat_p(ispin)%matrix, &
618 7490 : row=iatom, col=iatom, BLOCK=pa_block_a, found=found)
619 7490 : CPASSERT(ASSOCIATED(pa_block_a))
620 : ! Nuclei
621 7490 : IF (task(1) .AND. ispin == 1) se_nddo_mpole%charge(iatom) = se_kind_a%zeff
622 : ! Electrons
623 7490 : size_1c_int = SIZE(se_kind_a%w_mpole)
624 70294 : DO jint = 1, size_1c_int
625 62804 : mpole => se_kind_a%w_mpole(jint)%mpole
626 62804 : indi = se_orbital_pointer(mpole%indi)
627 62804 : indj = se_orbital_pointer(mpole%indj)
628 62804 : fac = 1.0_dp
629 62804 : IF (indi /= indj) fac = 2.0_dp
630 :
631 : ! Charge
632 62804 : IF (mpole%task(1) .AND. task(1)) THEN
633 : se_nddo_mpole%charge(iatom) = se_nddo_mpole%charge(iatom) + &
634 18092 : fac*pa_block_a(indi, indj)*mpole%c
635 : END IF
636 :
637 : ! Dipole
638 62804 : IF (mpole%task(2) .AND. task(2)) THEN
639 : se_nddo_mpole%dipole(:, iatom) = se_nddo_mpole%dipole(:, iatom) + &
640 89019 : fac*pa_block_a(indi, indj)*mpole%d(:)
641 : END IF
642 :
643 : ! Quadrupole
644 70294 : IF (mpole%task(3) .AND. task(3)) THEN
645 151740 : qsph = fac*mpole%qs*pa_block_a(indi, indj)
646 25290 : CALL quadrupole_sph_to_cart(qcart, qsph)
647 : se_nddo_mpole%quadrupole(:, :, iatom) = se_nddo_mpole%quadrupole(:, :, iatom) + &
648 328770 : qcart
649 : END IF
650 : END DO
651 : ! Print some info about charge, dipole and quadrupole (debug purpose only)
652 10398 : 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 31590 : CALL para_env%sum(se_nddo_mpole%charge)
660 121470 : CALL para_env%sum(se_nddo_mpole%dipole)
661 391110 : CALL para_env%sum(se_nddo_mpole%quadrupole)
662 :
663 : ! Initialize for virial
664 1630 : IF (use_virial) THEN
665 2 : pv_glob = 0.0_dp
666 2 : pv_local = 0.0_dp
667 : END IF
668 :
669 : ! Ewald Multipoles Sum
670 1630 : iw = cp_print_key_unit_nr(logger, se_section, "PRINT%EWALD_INFO", extension=".seLog")
671 1630 : IF (calculate_forces) THEN
672 20 : CALL get_qs_env(qs_env=qs_env, force=force)
673 20 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
674 :
675 : ! Allocate and zeroing arrays
676 59 : ALLOCATE (forces_g(3, forces_g_size))
677 60 : ALLOCATE (forces_r(3, natoms))
678 668 : forces_g = 0.0_dp
679 1316 : 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 20 : do_debug=.TRUE.)
688 : ! Only SR force have to be summed up.. the one in g-space are already fully local..
689 20 : 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 1610 : iw=iw, do_debug=.TRUE.)
698 : END IF
699 1630 : 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 1630 : 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 1400 : pv_glob)
706 : END IF
707 :
708 : ! Virial for the long-range part and correction
709 1630 : IF (use_virial) THEN
710 : ! Sum up contribution of pv_glob on each thread and keep only one copy of pv_local
711 26 : virial%pv_virial = virial%pv_virial + pv_glob
712 2 : IF (para_env%is_source()) THEN
713 13 : 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 4538 : DO ikind = 1, nkind
727 2908 : CALL get_qs_kind(qs_kind_set(ikind), se_parameter=se_kind_a)
728 2908 : CALL get_se_param(se_kind_a, defined=defined, natorb=natorb_a)
729 :
730 2908 : IF (.NOT. defined .OR. natorb_a < 1) CYCLE
731 :
732 2908 : nparticle_local = local_particles%n_el(ikind)
733 14936 : DO ilist = 1, nparticle_local
734 7490 : node = node + 1
735 7490 : iatom = local_particles%list(ikind)%array(ilist)
736 14980 : DO ispin = 1, nspins
737 : CALL dbcsr_get_block_p(matrix=diagmat_ks(ispin)%matrix, &
738 7490 : row=iatom, col=iatom, BLOCK=ksa_block_a, found=found)
739 7490 : CPASSERT(ASSOCIATED(ksa_block_a))
740 :
741 : ! Modify Hamiltonian Matrix accordingly potential, field and electric field gradient
742 7490 : size_1c_int = SIZE(se_kind_a%w_mpole)
743 85274 : DO jint = 1, size_1c_int
744 62804 : tmp = 0.0_dp
745 62804 : mpole => se_kind_a%w_mpole(jint)%mpole
746 62804 : indi = se_orbital_pointer(mpole%indi)
747 62804 : indj = se_orbital_pointer(mpole%indj)
748 :
749 : ! Charge
750 62804 : IF (mpole%task(1) .AND. task(1)) THEN
751 18092 : tmp = tmp + mpole%c*se_nddo_mpole%efield0(iatom)
752 : END IF
753 :
754 : ! Dipole
755 62804 : IF (mpole%task(2) .AND. task(2)) THEN
756 50868 : tmp = tmp - DOT_PRODUCT(mpole%d, se_nddo_mpole%efield1(:, iatom))
757 : END IF
758 :
759 : ! Quadrupole
760 62804 : IF (mpole%task(3) .AND. task(3)) THEN
761 328770 : tmp = tmp - (1.0_dp/3.0_dp)*SUM(mpole%qc*RESHAPE(se_nddo_mpole%efield2(:, iatom), (/3, 3/)))
762 : END IF
763 :
764 62804 : ksa_block_a(indi, indj) = ksa_block_a(indi, indj) + tmp
765 70294 : ksa_block_a(indj, indi) = ksa_block_a(indi, indj)
766 : END DO
767 : END DO
768 :
769 : ! Nuclear term and forces
770 7490 : IF (task(1)) enuc = enuc + se_kind_a%zeff*se_nddo_mpole%efield0(iatom)
771 7490 : IF (atener) THEN
772 : atprop%atecoul(iatom) = atprop%atecoul(iatom) + &
773 0 : 0.5_dp*se_kind_a%zeff*se_nddo_mpole%efield0(iatom)
774 : END IF
775 10398 : IF (calculate_forces) THEN
776 162 : atom_a = atom_of_kind(iatom)
777 648 : force_a = forces_r(1:3, iatom) + forces_g(1:3, node)
778 : ! Derivatives of the periodic Coulomb Terms
779 648 : 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 1630 : CALL para_env%sum(enuc)
785 1630 : 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 3260 : DO ispin = 1, nspins
793 1630 : CALL dbcsr_sum_replicated(diagmat_ks(ispin)%matrix)
794 1630 : CALL dbcsr_distribute(diagmat_ks(ispin)%matrix)
795 : CALL dbcsr_add(ks_matrix(ispin)%matrix, diagmat_ks(ispin)%matrix, &
796 3260 : 1.0_dp, 1.0_dp)
797 : END DO
798 1630 : CALL dbcsr_deallocate_matrix_set(diagmat_p)
799 1630 : CALL dbcsr_deallocate_matrix_set(diagmat_ks)
800 :
801 : ! Set the Fock matrix contribution to SCP
802 1630 : IF (calculate_forces) THEN
803 20 : DEALLOCATE (forces_g)
804 20 : DEALLOCATE (forces_r)
805 : END IF
806 :
807 1630 : CALL timestop(handle)
808 3260 : 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 1400 : 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 1400 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, se_natorb
848 : LOGICAL :: anag, atener, check, defined, found, &
849 : switch, use_virial
850 1400 : 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 1400 : REAL(KIND=dp), DIMENSION(:), POINTER :: efield0
860 1400 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: efield1, efield2, ksa_block_a, ksa_block_b, &
861 1400 : ksb_block_a, ksb_block_b, pa_block_a, pa_block_b, pb_block_a, pb_block_b
862 1400 : 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 1400 : DIMENSION(:), POINTER :: nl_iterator
869 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
870 1400 : POINTER :: sab_lrc
871 1400 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
872 1400 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
873 1400 : 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 1400 : TYPE(semi_empirical_p_type), DIMENSION(:), POINTER :: se_kind_list
878 : TYPE(semi_empirical_type), POINTER :: se_kind_a, se_kind_b
879 :
880 1400 : CALL timeset(routineN, handle)
881 1400 : NULLIFY (dft_control, cell, force, particle_set, se_control, se_taper, &
882 1400 : 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 1400 : qs_kind_set=qs_kind_set, particle_set=particle_set, atprop=atprop)
887 :
888 : ! Parameters
889 1400 : CALL initialize_se_taper(se_taper, lr_corr=.TRUE.)
890 1400 : se_control => dft_control%qs_control%se_control
891 1400 : anag = se_control%analytical_gradients
892 1400 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer) .AND. calculate_forces
893 1400 : 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 1400 : pc_coulomb_int=.TRUE.)
899 :
900 1400 : nspins = dft_control%nspins
901 1400 : CPASSERT(ASSOCIATED(matrix_p))
902 1400 : CPASSERT(SIZE(ks_matrix) > 0)
903 1400 : CPASSERT(ASSOCIATED(diagmat_p))
904 1400 : CPASSERT(ASSOCIATED(diagmat_ks))
905 : MARK_USED(ks_matrix)
906 : MARK_USED(matrix_p)
907 :
908 1400 : nkind = SIZE(atomic_kind_set)
909 1400 : IF (calculate_forces) THEN
910 20 : CALL get_qs_env(qs_env=qs_env, force=force)
911 20 : delta = se_control%delta
912 20 : 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 1400 : size1 = SIZE(se_nddo_mpole%efield0)
917 4200 : ALLOCATE (efield0(size1))
918 12582 : efield0 = 0.0_dp
919 1400 : size1 = SIZE(se_nddo_mpole%efield1, 1)
920 1400 : size2 = SIZE(se_nddo_mpole%efield1, 2)
921 5600 : ALLOCATE (efield1(size1, size2))
922 46128 : efield1 = 0.0_dp
923 1400 : size1 = SIZE(se_nddo_mpole%efield2, 1)
924 1400 : size2 = SIZE(se_nddo_mpole%efield2, 2)
925 5600 : ALLOCATE (efield2(size1, size2))
926 113220 : efield2 = 0.0_dp
927 :
928 : ! Initialize if virial is requested
929 1400 : IF (use_virial) THEN
930 2 : ptens11 = 0.0_dp; ptens12 = 0.0_dp; ptens13 = 0.0_dp
931 2 : ptens21 = 0.0_dp; ptens22 = 0.0_dp; ptens23 = 0.0_dp
932 2 : 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 1400 : ecore2 = 0.0_dp
937 1400 : enuclear = 0.0_dp
938 1400 : itype = get_se_type(dft_control%qs_control%method_id)
939 :
940 10848 : ALLOCATE (se_defined(nkind), se_kind_list(nkind), se_natorb(nkind))
941 3848 : DO ikind = 1, nkind
942 2448 : CALL get_qs_kind(qs_kind_set(ikind), se_parameter=se_kind_a)
943 2448 : se_kind_list(ikind)%se_param => se_kind_a
944 2448 : CALL get_se_param(se_kind_a, defined=defined, natorb=natorb_a)
945 2448 : se_defined(ikind) = (defined .AND. natorb_a >= 1)
946 6296 : se_natorb(ikind) = natorb_a
947 : END DO
948 :
949 1400 : CALL neighbor_list_iterator_create(nl_iterator, sab_lrc)
950 1432964 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
951 1431564 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, inode=inode, r=rij)
952 1431564 : IF (.NOT. se_defined(ikind)) CYCLE
953 1431564 : IF (.NOT. se_defined(jkind)) CYCLE
954 1431564 : se_kind_a => se_kind_list(ikind)%se_param
955 1431564 : se_kind_b => se_kind_list(jkind)%se_param
956 1431564 : natorb_a = se_natorb(ikind)
957 1431564 : natorb_b = se_natorb(jkind)
958 1431564 : natorb_a2 = natorb_a**2
959 1431564 : natorb_b2 = natorb_b**2
960 :
961 1431564 : IF (inode == 1) THEN
962 : CALL dbcsr_get_block_p(matrix=diagmat_p(1)%matrix, &
963 9998 : row=iatom, col=iatom, BLOCK=pa_block_a, found=found)
964 9998 : CPASSERT(ASSOCIATED(pa_block_a))
965 9998 : check = (SIZE(pa_block_a, 1) == natorb_a) .AND. (SIZE(pa_block_a, 2) == natorb_a)
966 0 : CPASSERT(check)
967 : CALL dbcsr_get_block_p(matrix=diagmat_ks(1)%matrix, &
968 9998 : row=iatom, col=iatom, BLOCK=ksa_block_a, found=found)
969 9998 : CPASSERT(ASSOCIATED(ksa_block_a))
970 166686 : p_block_tot_a(1:natorb_a, 1:natorb_a) = 2.0_dp*pa_block_a
971 19996 : pa_a(1:natorb_a2) = RESHAPE(pa_block_a, (/natorb_a2/))
972 9998 : IF (nspins == 2) THEN
973 : CALL dbcsr_get_block_p(matrix=diagmat_p(2)%matrix, &
974 0 : row=iatom, col=iatom, BLOCK=pa_block_b, found=found)
975 0 : CPASSERT(ASSOCIATED(pa_block_b))
976 0 : check = (SIZE(pa_block_b, 1) == natorb_a) .AND. (SIZE(pa_block_b, 2) == natorb_a)
977 0 : CPASSERT(check)
978 : CALL dbcsr_get_block_p(matrix=diagmat_ks(2)%matrix, &
979 0 : row=iatom, col=iatom, BLOCK=ksa_block_b, found=found)
980 0 : CPASSERT(ASSOCIATED(ksa_block_b))
981 0 : p_block_tot_a(1:natorb_a, 1:natorb_a) = pa_block_a + pa_block_b
982 0 : pa_b(1:natorb_a2) = RESHAPE(pa_block_b, (/natorb_a2/))
983 : END IF
984 :
985 : END IF
986 :
987 5726256 : dr1 = DOT_PRODUCT(rij, rij)
988 1432964 : IF (dr1 > rij_threshold) THEN
989 : ! Determine the order of the atoms, and in case switch them..
990 1425973 : IF (iatom <= jatom) THEN
991 767729 : switch = .FALSE.
992 : ELSE
993 658244 : 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 1425973 : ptens32=ptens32, ptens33=ptens33)
1004 :
1005 : ! Retrieve blocks for KS and P
1006 : CALL dbcsr_get_block_p(matrix=diagmat_p(1)%matrix, &
1007 1425973 : row=jatom, col=jatom, BLOCK=pb_block_a, found=found)
1008 1425973 : CPASSERT(ASSOCIATED(pb_block_a))
1009 1425973 : check = (SIZE(pb_block_a, 1) == natorb_b) .AND. (SIZE(pb_block_a, 2) == natorb_b)
1010 0 : CPASSERT(check)
1011 : CALL dbcsr_get_block_p(matrix=diagmat_ks(1)%matrix, &
1012 1425973 : row=jatom, col=jatom, BLOCK=ksb_block_a, found=found)
1013 1425973 : CPASSERT(ASSOCIATED(ksb_block_a))
1014 23202075 : p_block_tot_b(1:natorb_b, 1:natorb_b) = 2.0_dp*pb_block_a
1015 2851946 : pb_a(1:natorb_b2) = RESHAPE(pb_block_a, (/natorb_b2/))
1016 : ! Handle more than one configuration
1017 1425973 : IF (nspins == 2) THEN
1018 : CALL dbcsr_get_block_p(matrix=diagmat_p(2)%matrix, &
1019 0 : row=jatom, col=jatom, BLOCK=pb_block_b, found=found)
1020 0 : CPASSERT(ASSOCIATED(pb_block_b))
1021 0 : check = (SIZE(pb_block_b, 1) == natorb_b) .AND. (SIZE(pb_block_b, 2) == natorb_b)
1022 0 : CPASSERT(check)
1023 : CALL dbcsr_get_block_p(matrix=diagmat_ks(2)%matrix, &
1024 0 : row=jatom, col=jatom, BLOCK=ksb_block_b, found=found)
1025 0 : CPASSERT(ASSOCIATED(ksb_block_b))
1026 0 : check = (SIZE(pb_block_a, 1) == SIZE(pb_block_b, 1)) .AND. (SIZE(pb_block_a, 2) == SIZE(pb_block_b, 2))
1027 0 : CPASSERT(check)
1028 0 : p_block_tot_b(1:natorb_b, 1:natorb_b) = pb_block_a + pb_block_b
1029 0 : pb_b(1:natorb_b2) = RESHAPE(pb_block_b, (/natorb_b2/))
1030 : END IF
1031 :
1032 2851946 : 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 1425973 : se_int_control=se_int_control, se_taper=se_taper)
1038 1425973 : enuclear = enuclear + enuc
1039 :
1040 : ! Two-centers One-electron terms
1041 1425973 : IF (nspins == 1) THEN
1042 1425973 : 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 1425973 : se_taper=se_taper, store_int_env=store_int_env)
1046 1425973 : ecore2 = ecore2 + ecab(1) + ecab(2)
1047 0 : ELSE IF (nspins == 2) THEN
1048 0 : 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 0 : 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 0 : se_taper=se_taper, store_int_env=store_int_env)
1055 0 : ecore2 = ecore2 + ecab(1) + ecab(2)
1056 : END IF
1057 1425973 : IF (atener) THEN
1058 0 : atprop%atecoul(iatom) = atprop%atecoul(iatom) + ecab(1) + 0.5_dp*enuc
1059 0 : atprop%atecoul(jatom) = atprop%atecoul(jatom) + ecab(2) + 0.5_dp*enuc
1060 : END IF
1061 : ! Coulomb Terms
1062 1425973 : 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 1425973 : store_int_env=store_int_env)
1066 0 : 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 0 : 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 0 : store_int_env=store_int_env)
1074 : END IF
1075 :
1076 1425973 : IF (calculate_forces) THEN
1077 43876 : atom_a = atom_of_kind(iatom)
1078 43876 : 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 43876 : anag=anag, se_int_control=se_int_control, se_taper=se_taper)
1083 :
1084 : ! Derivatives of the Two-centre One-electron terms
1085 43876 : 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 43876 : delta=delta)
1089 0 : 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 0 : 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 0 : se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta)
1094 : END IF
1095 43876 : IF (use_virial) THEN
1096 2376 : CALL virial_pair_force(virial%pv_virial, -1.0_dp, force_ab, rij)
1097 : END IF
1098 175504 : force_ab = force_ab + force_ab0
1099 :
1100 : ! Sum up force components
1101 43876 : force(ikind)%all_potential(1, atom_a) = force(ikind)%all_potential(1, atom_a) - force_ab(1)
1102 43876 : force(jkind)%all_potential(1, atom_b) = force(jkind)%all_potential(1, atom_b) + force_ab(1)
1103 :
1104 43876 : force(ikind)%all_potential(2, atom_a) = force(ikind)%all_potential(2, atom_a) - force_ab(2)
1105 43876 : force(jkind)%all_potential(2, atom_b) = force(jkind)%all_potential(2, atom_b) + force_ab(2)
1106 :
1107 43876 : force(ikind)%all_potential(3, atom_a) = force(ikind)%all_potential(3, atom_a) - force_ab(3)
1108 43876 : 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 43876 : force_ab = 0._dp
1112 43876 : 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 43876 : anag=anag, se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta)
1115 0 : 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 0 : 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 0 : anag=anag, se_int_control=se_int_control, se_taper=se_taper, force=force_ab, delta=delta)
1121 : END IF
1122 43876 : IF (switch) THEN
1123 21260 : force_ab(1) = -force_ab(1)
1124 21260 : force_ab(2) = -force_ab(2)
1125 21260 : force_ab(3) = -force_ab(3)
1126 : END IF
1127 43876 : IF (use_virial) THEN
1128 2376 : CALL virial_pair_force(virial%pv_virial, -1.0_dp, force_ab, rij)
1129 : END IF
1130 :
1131 : ! Sum up force components
1132 43876 : force(ikind)%rho_elec(1, atom_a) = force(ikind)%rho_elec(1, atom_a) - force_ab(1)
1133 43876 : force(jkind)%rho_elec(1, atom_b) = force(jkind)%rho_elec(1, atom_b) + force_ab(1)
1134 :
1135 43876 : force(ikind)%rho_elec(2, atom_a) = force(ikind)%rho_elec(2, atom_a) - force_ab(2)
1136 43876 : force(jkind)%rho_elec(2, atom_b) = force(jkind)%rho_elec(2, atom_b) + force_ab(2)
1137 :
1138 43876 : force(ikind)%rho_elec(3, atom_a) = force(ikind)%rho_elec(3, atom_a) - force_ab(3)
1139 43876 : force(jkind)%rho_elec(3, atom_b) = force(jkind)%rho_elec(3, atom_b) + force_ab(3)
1140 : END IF
1141 : CASE DEFAULT
1142 1425973 : CPABORT("")
1143 : END SELECT
1144 : END IF
1145 : END DO
1146 1400 : CALL neighbor_list_iterator_release(nl_iterator)
1147 :
1148 1400 : DEALLOCATE (se_kind_list, se_defined, se_natorb)
1149 :
1150 : ! Sum-up Virial constribution (long-range correction)
1151 1400 : IF (use_virial) THEN
1152 2 : pv_glob(1, 1) = pv_glob(1, 1) + ptens11
1153 2 : pv_glob(1, 2) = pv_glob(1, 2) + (ptens12 + ptens21)*0.5_dp
1154 2 : pv_glob(1, 3) = pv_glob(1, 3) + (ptens13 + ptens31)*0.5_dp
1155 2 : pv_glob(2, 1) = pv_glob(1, 2)
1156 2 : pv_glob(2, 2) = pv_glob(2, 2) + ptens22
1157 2 : pv_glob(2, 3) = pv_glob(2, 3) + (ptens23 + ptens32)*0.5_dp
1158 2 : pv_glob(3, 1) = pv_glob(1, 3)
1159 2 : pv_glob(3, 2) = pv_glob(2, 3)
1160 2 : pv_glob(3, 3) = pv_glob(3, 3) + ptens33
1161 : END IF
1162 :
1163 : ! collect information on potential, field, field gradient
1164 23764 : CALL para_env%sum(efield0)
1165 90856 : CALL para_env%sum(efield1)
1166 225040 : CALL para_env%sum(efield2)
1167 23764 : se_nddo_mpole%efield0 = se_nddo_mpole%efield0 - efield0
1168 90856 : se_nddo_mpole%efield1 = se_nddo_mpole%efield1 - efield1
1169 225040 : se_nddo_mpole%efield2 = se_nddo_mpole%efield2 - efield2
1170 : ! deallocate working arrays
1171 1400 : DEALLOCATE (efield0)
1172 1400 : DEALLOCATE (efield1)
1173 1400 : DEALLOCATE (efield2)
1174 :
1175 : ! Corrections for two-centers one-electron terms + nuclear terms
1176 1400 : CALL para_env%sum(enuclear)
1177 1400 : CALL para_env%sum(ecore2)
1178 1400 : energy%hartree = energy%hartree + ecore2
1179 1400 : energy%core_overlap = enuclear
1180 1400 : CALL finalize_se_taper(se_taper)
1181 1400 : CALL timestop(handle)
1182 2800 : 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 0 : 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 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, se_natorb
1207 : LOGICAL :: anag, atener, check, defined, found, &
1208 : switch
1209 0 : 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 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: ksa_block_a, ksa_block_b, ksb_block_a, &
1216 0 : ksb_block_b, pa_block_a, pa_block_b, &
1217 0 : pb_block_a, pb_block_b
1218 0 : 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 0 : 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 0 : DIMENSION(:), POINTER :: nl_iterator
1232 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
1233 0 : POINTER :: sab_orb
1234 0 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
1235 0 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
1236 0 : 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 0 : TYPE(semi_empirical_p_type), DIMENSION(:), POINTER :: se_kind_list
1240 : TYPE(semi_empirical_type), POINTER :: se_kind_a, se_kind_b
1241 :
1242 0 : CALL timeset(routineN, handle)
1243 :
1244 0 : CALL get_qs_env(qs_env=qs_env) !sm->dbcsr
1245 :
1246 0 : NULLIFY (dft_control, cell, force, particle_set, diagmat_ks, &
1247 0 : diagmat_p, local_particles, se_control, ewald_env, ewald_pw, &
1248 0 : se_nddo_mpole, se_nonbond_env, se_section, sab_orb, logger, atprop)
1249 :
1250 0 : 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 0 : se_nonbond_env=se_nonbond_env, sab_orb=sab_orb)
1256 :
1257 0 : nlocal_particles = SUM(local_particles%n_el(:))
1258 0 : natoms = SIZE(particle_set)
1259 0 : CALL ewald_env_get(ewald_env, ewald_type=ewald_type)
1260 0 : SELECT CASE (ewald_type)
1261 : CASE (do_ewald_ewald)
1262 : ! Do Nothing
1263 : CASE DEFAULT
1264 0 : CPABORT("Periodic SE implemented only for standard EWALD sums.")
1265 : END SELECT
1266 :
1267 : ! Parameters
1268 0 : se_section => section_vals_get_subs_vals(qs_env%input, "DFT%QS%SE")
1269 0 : se_control => dft_control%qs_control%se_control
1270 0 : anag = se_control%analytical_gradients
1271 0 : atener = atprop%energy
1272 :
1273 0 : nspins = dft_control%nspins
1274 0 : CPASSERT(ASSOCIATED(matrix_p))
1275 0 : CPASSERT(SIZE(ks_matrix) > 0)
1276 :
1277 0 : CALL dbcsr_allocate_matrix_set(diagmat_ks, nspins)
1278 0 : CALL dbcsr_allocate_matrix_set(diagmat_p, nspins)
1279 :
1280 0 : nkind = SIZE(atomic_kind_set)
1281 0 : DO ispin = 1, nspins
1282 : ! Allocate diagonal block matrices
1283 0 : ALLOCATE (diagmat_p(ispin)%matrix, diagmat_ks(ispin)%matrix) !sm->dbcsr
1284 0 : CALL dbcsr_get_block_diag(matrix_p(ispin)%matrix, diagmat_p(ispin)%matrix)
1285 0 : CALL dbcsr_get_block_diag(ks_matrix(ispin)%matrix, diagmat_ks(ispin)%matrix)
1286 0 : CALL dbcsr_set(diagmat_ks(ispin)%matrix, 0.0_dp)
1287 0 : CALL dbcsr_replicate_all(diagmat_p(ispin)%matrix)
1288 0 : CALL dbcsr_replicate_all(diagmat_ks(ispin)%matrix)
1289 : END DO
1290 :
1291 : ! Possibly compute forces
1292 0 : IF (calculate_forces) THEN
1293 0 : CALL get_qs_env(qs_env=qs_env, force=force)
1294 0 : 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 0 : ecore2 = 0.0_dp
1299 : ! Real space part of the 1/R^3 sum
1300 0 : ALLOCATE (se_defined(nkind), se_kind_list(nkind), se_natorb(nkind))
1301 0 : DO ikind = 1, nkind
1302 0 : CALL get_qs_kind(qs_kind_set(ikind), se_parameter=se_kind_a)
1303 0 : se_kind_list(ikind)%se_param => se_kind_a
1304 0 : CALL get_se_param(se_kind_a, defined=defined, natorb=natorb_a)
1305 0 : se_defined(ikind) = (defined .AND. natorb_a >= 1)
1306 0 : se_natorb(ikind) = natorb_a
1307 : END DO
1308 :
1309 0 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
1310 0 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
1311 0 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, inode=inode, r=rij)
1312 0 : IF (.NOT. se_defined(ikind)) CYCLE
1313 0 : IF (.NOT. se_defined(jkind)) CYCLE
1314 0 : se_kind_a => se_kind_list(ikind)%se_param
1315 0 : se_kind_b => se_kind_list(jkind)%se_param
1316 0 : natorb_a = se_natorb(ikind)
1317 0 : natorb_b = se_natorb(jkind)
1318 0 : natorb_a2 = natorb_a**2
1319 0 : natorb_b2 = natorb_b**2
1320 :
1321 0 : IF (inode == 1) THEN
1322 : CALL dbcsr_get_block_p(matrix=diagmat_p(1)%matrix, &
1323 0 : row=iatom, col=iatom, BLOCK=pa_block_a, found=found)
1324 0 : CPASSERT(ASSOCIATED(pa_block_a))
1325 0 : check = (SIZE(pa_block_a, 1) == natorb_a) .AND. (SIZE(pa_block_a, 2) == natorb_a)
1326 0 : CPASSERT(check)
1327 : CALL dbcsr_get_block_p(matrix=diagmat_ks(1)%matrix, &
1328 0 : row=iatom, col=iatom, BLOCK=ksa_block_a, found=found)
1329 0 : CPASSERT(ASSOCIATED(ksa_block_a))
1330 0 : p_block_tot_a(1:natorb_a, 1:natorb_a) = 2.0_dp*pa_block_a
1331 0 : pa_a(1:natorb_a2) = RESHAPE(pa_block_a, (/natorb_a2/))
1332 0 : IF (nspins == 2) THEN
1333 : CALL dbcsr_get_block_p(matrix=diagmat_p(2)%matrix, &
1334 0 : row=iatom, col=iatom, BLOCK=pa_block_b, found=found)
1335 0 : CPASSERT(ASSOCIATED(pa_block_b))
1336 0 : check = (SIZE(pa_block_b, 1) == natorb_a) .AND. (SIZE(pa_block_b, 2) == natorb_a)
1337 0 : CPASSERT(check)
1338 : CALL dbcsr_get_block_p(matrix=diagmat_ks(2)%matrix, &
1339 0 : row=iatom, col=iatom, BLOCK=ksa_block_b, found=found)
1340 0 : CPASSERT(ASSOCIATED(ksa_block_b))
1341 0 : p_block_tot_a(1:natorb_a, 1:natorb_a) = pa_block_a + pa_block_b
1342 0 : pa_b(1:natorb_a2) = RESHAPE(pa_block_b, (/natorb_a2/))
1343 : END IF
1344 : END IF
1345 :
1346 0 : dr1 = DOT_PRODUCT(rij, rij)
1347 0 : IF (dr1 > rij_threshold) THEN
1348 : ! Determine the order of the atoms, and in case switch them..
1349 0 : IF (iatom <= jatom) THEN
1350 0 : switch = .FALSE.
1351 : ELSE
1352 0 : switch = .TRUE.
1353 : END IF
1354 : ! Retrieve blocks for KS and P
1355 : CALL dbcsr_get_block_p(matrix=diagmat_p(1)%matrix, &
1356 0 : row=jatom, col=jatom, BLOCK=pb_block_a, found=found)
1357 0 : CPASSERT(ASSOCIATED(pb_block_a))
1358 0 : check = (SIZE(pb_block_a, 1) == natorb_b) .AND. (SIZE(pb_block_a, 2) == natorb_b)
1359 0 : CPASSERT(check)
1360 : CALL dbcsr_get_block_p(matrix=diagmat_ks(1)%matrix, &
1361 0 : row=jatom, col=jatom, BLOCK=ksb_block_a, found=found)
1362 0 : CPASSERT(ASSOCIATED(ksb_block_a))
1363 0 : p_block_tot_b(1:natorb_b, 1:natorb_b) = 2.0_dp*pb_block_a
1364 0 : pb_a(1:natorb_b2) = RESHAPE(pb_block_a, (/natorb_b2/))
1365 : ! Handle more than one configuration
1366 0 : IF (nspins == 2) THEN
1367 : CALL dbcsr_get_block_p(matrix=diagmat_p(2)%matrix, &
1368 0 : row=jatom, col=jatom, BLOCK=pb_block_b, found=found)
1369 0 : CPASSERT(ASSOCIATED(pb_block_b))
1370 0 : check = (SIZE(pb_block_b, 1) == natorb_b) .AND. (SIZE(pb_block_b, 2) == natorb_b)
1371 0 : CPASSERT(check)
1372 : CALL dbcsr_get_block_p(matrix=diagmat_ks(2)%matrix, &
1373 0 : row=jatom, col=jatom, BLOCK=ksb_block_b, found=found)
1374 0 : CPASSERT(ASSOCIATED(ksb_block_b))
1375 0 : check = (SIZE(pb_block_a, 1) == SIZE(pb_block_b, 1)) .AND. (SIZE(pb_block_a, 2) == SIZE(pb_block_b, 2))
1376 0 : CPASSERT(check)
1377 0 : p_block_tot_b(1:natorb_b, 1:natorb_b) = pb_block_a + pb_block_b
1378 0 : pb_b(1:natorb_b2) = RESHAPE(pb_block_b, (/natorb_b2/))
1379 : END IF
1380 :
1381 0 : 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 0 : r2inv = 1.0_dp/dr1
1387 0 : rinv = SQRT(r2inv)
1388 0 : r3inv = rinv**3
1389 :
1390 : ! Two-centers One-electron terms
1391 0 : IF (nspins == 1) THEN
1392 0 : 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 0 : e2a=se_kind_a%expns3_int(jkind)%expns3%e2a, rp=r3inv)
1396 0 : ecore2 = ecore2 + ecab(1) + ecab(2)
1397 0 : ELSE IF (nspins == 2) THEN
1398 0 : 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 0 : 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 0 : e2a=se_kind_a%expns3_int(jkind)%expns3%e2a, rp=r3inv)
1406 0 : ecore2 = ecore2 + ecab(1) + ecab(2)
1407 : END IF
1408 0 : IF (atener) THEN
1409 0 : atprop%atecoul(iatom) = atprop%atecoul(iatom) + ecab(1)
1410 0 : atprop%atecoul(jatom) = atprop%atecoul(jatom) + ecab(2)
1411 : END IF
1412 : ! Coulomb Terms
1413 0 : 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 0 : ksb_block_a, factor=0.5_dp, w=se_kind_a%expns3_int(jkind)%expns3%w, rp=r3inv)
1416 0 : 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 0 : 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 0 : 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 0 : IF (calculate_forces) THEN
1426 0 : dr3inv = -3.0_dp*rij*r3inv*r2inv
1427 0 : atom_a = atom_of_kind(iatom)
1428 0 : atom_b = atom_of_kind(jatom)
1429 :
1430 0 : force_ab = 0.0_dp
1431 : ! Derivatives of the One-centre One-electron terms
1432 0 : IF (nspins == 1) THEN
1433 : CALL dfock2_1el_r3(se_kind_a, se_kind_b, dr3inv, pa_a, pb_a, force_ab, &
1434 0 : se_kind_a%expns3_int(jkind)%expns3%e1b, se_kind_a%expns3_int(jkind)%expns3%e2a)
1435 0 : 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 0 : 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 0 : 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 0 : force(ikind)%all_potential(1, atom_a) = force(ikind)%all_potential(1, atom_a) - force_ab(1)
1445 0 : force(jkind)%all_potential(1, atom_b) = force(jkind)%all_potential(1, atom_b) + force_ab(1)
1446 :
1447 0 : force(ikind)%all_potential(2, atom_a) = force(ikind)%all_potential(2, atom_a) - force_ab(2)
1448 0 : force(jkind)%all_potential(2, atom_b) = force(jkind)%all_potential(2, atom_b) + force_ab(2)
1449 :
1450 0 : force(ikind)%all_potential(3, atom_a) = force(ikind)%all_potential(3, atom_a) - force_ab(3)
1451 0 : 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 0 : force_ab = 0.0_dp
1455 0 : 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 0 : w=se_kind_a%expns3_int(jkind)%expns3%w, drp=dr3inv, force=force_ab)
1458 0 : 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 0 : 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 0 : w=se_kind_a%expns3_int(jkind)%expns3%w, drp=dr3inv, force=force_ab)
1464 : END IF
1465 :
1466 : ! Sum up force components
1467 0 : force(ikind)%rho_elec(1, atom_a) = force(ikind)%rho_elec(1, atom_a) - force_ab(1)
1468 0 : force(jkind)%rho_elec(1, atom_b) = force(jkind)%rho_elec(1, atom_b) + force_ab(1)
1469 :
1470 0 : force(ikind)%rho_elec(2, atom_a) = force(ikind)%rho_elec(2, atom_a) - force_ab(2)
1471 0 : force(jkind)%rho_elec(2, atom_b) = force(jkind)%rho_elec(2, atom_b) + force_ab(2)
1472 :
1473 0 : force(ikind)%rho_elec(3, atom_a) = force(ikind)%rho_elec(3, atom_a) - force_ab(3)
1474 0 : force(jkind)%rho_elec(3, atom_b) = force(jkind)%rho_elec(3, atom_b) + force_ab(3)
1475 : END IF
1476 : CASE DEFAULT
1477 0 : CPABORT("")
1478 : END SELECT
1479 : END IF
1480 : END DO
1481 0 : CALL neighbor_list_iterator_release(nl_iterator)
1482 :
1483 0 : DEALLOCATE (se_kind_list, se_defined, se_natorb)
1484 :
1485 0 : DO ispin = 1, nspins
1486 0 : CALL dbcsr_sum_replicated(diagmat_ks(ispin)%matrix)
1487 0 : CALL dbcsr_distribute(diagmat_ks(ispin)%matrix)
1488 : CALL dbcsr_add(ks_matrix(ispin)%matrix, diagmat_ks(ispin)%matrix, &
1489 0 : 1.0_dp, 1.0_dp)
1490 : END DO
1491 0 : CALL dbcsr_deallocate_matrix_set(diagmat_p)
1492 0 : CALL dbcsr_deallocate_matrix_set(diagmat_ks)
1493 :
1494 : ! Two-centers one-electron terms
1495 0 : CALL para_env%sum(ecore2)
1496 0 : energy%hartree = energy%hartree + ecore2
1497 :
1498 0 : CALL timestop(handle)
1499 0 : END SUBROUTINE build_fock_matrix_coul_lr_r3
1500 :
1501 : END MODULE se_fock_matrix_coulomb
1502 :
|
