Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Construction of the Exchange part of the Fock Matrix
10 : !> \author Teodoro Laino [tlaino] (05.2009) - Split and module reorganization
11 : !> \par History
12 : !> Teodoro Laino (04.2008) [tlaino] - University of Zurich : d-orbitals
13 : !> Teodoro Laino (09.2008) [tlaino] - University of Zurich : Speed-up
14 : !> Teodoro Laino (09.2008) [tlaino] - University of Zurich : Periodic SE
15 : ! **************************************************************************************************
16 : MODULE se_fock_matrix_exchange
17 : USE atomic_kind_types, ONLY: atomic_kind_type,&
18 : get_atomic_kind_set
19 : USE cell_types, ONLY: cell_type
20 : USE cp_control_types, ONLY: dft_control_type,&
21 : semi_empirical_control_type
22 : USE cp_dbcsr_api, ONLY: dbcsr_get_block_p,&
23 : dbcsr_p_type
24 : USE input_constants, ONLY: do_se_IS_kdso,&
25 : do_se_IS_kdso_d
26 : USE kinds, ONLY: dp
27 : USE message_passing, ONLY: mp_para_env_type
28 : USE multipole_types, ONLY: do_multipole_none
29 : USE particle_types, ONLY: particle_type
30 : USE qs_environment_types, ONLY: get_qs_env,&
31 : qs_environment_type
32 : USE qs_force_types, ONLY: qs_force_type
33 : USE qs_kind_types, ONLY: get_qs_kind,&
34 : qs_kind_type
35 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
36 : neighbor_list_iterate,&
37 : neighbor_list_iterator_create,&
38 : neighbor_list_iterator_p_type,&
39 : neighbor_list_iterator_release,&
40 : neighbor_list_set_p_type
41 : USE se_fock_matrix_integrals, ONLY: dfock2E,&
42 : fock1_2el,&
43 : fock2E
44 : USE semi_empirical_int_arrays, ONLY: rij_threshold
45 : USE semi_empirical_store_int_types, ONLY: semi_empirical_si_type
46 : USE semi_empirical_types, ONLY: get_se_param,&
47 : se_int_control_type,&
48 : se_taper_type,&
49 : semi_empirical_p_type,&
50 : semi_empirical_type,&
51 : setup_se_int_control_type
52 : USE semi_empirical_utils, ONLY: finalize_se_taper,&
53 : initialize_se_taper
54 : USE virial_methods, ONLY: virial_pair_force
55 : USE virial_types, ONLY: virial_type
56 : #include "./base/base_uses.f90"
57 :
58 : IMPLICIT NONE
59 : PRIVATE
60 :
61 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'se_fock_matrix_exchange'
62 : LOGICAL, PARAMETER, PRIVATE :: debug_this_module = .FALSE.
63 :
64 : PUBLIC :: build_fock_matrix_exchange
65 :
66 : CONTAINS
67 :
68 : ! **************************************************************************************************
69 : !> \brief Construction of the Exchange part of the Fock matrix
70 : !> \param qs_env ...
71 : !> \param ks_matrix ...
72 : !> \param matrix_p ...
73 : !> \param calculate_forces ...
74 : !> \param store_int_env ...
75 : !> \author JGH
76 : ! **************************************************************************************************
77 41094 : SUBROUTINE build_fock_matrix_exchange(qs_env, ks_matrix, matrix_p, calculate_forces, &
78 : store_int_env)
79 :
80 : TYPE(qs_environment_type), POINTER :: qs_env
81 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_matrix, matrix_p
82 : LOGICAL, INTENT(in) :: calculate_forces
83 : TYPE(semi_empirical_si_type), POINTER :: store_int_env
84 :
85 : CHARACTER(len=*), PARAMETER :: routineN = 'build_fock_matrix_exchange'
86 :
87 : INTEGER :: atom_a, atom_b, handle, iatom, icol, &
88 : ikind, integral_screening, irow, &
89 : jatom, jkind, natorb_a, nkind, nspins
90 41094 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
91 : INTEGER, DIMENSION(2) :: size_p_block_a
92 : LOGICAL :: anag, check, defined, found, switch, &
93 : use_virial
94 41094 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: se_defined
95 : REAL(KIND=dp) :: delta, dr
96 : REAL(KIND=dp), DIMENSION(3) :: force_ab, rij
97 : REAL(KIND=dp), DIMENSION(45, 45) :: p_block_tot
98 41094 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: ks_block_a, ks_block_b, p_block_a, &
99 41094 : p_block_b
100 41094 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
101 : TYPE(cell_type), POINTER :: cell
102 : TYPE(dft_control_type), POINTER :: dft_control
103 : TYPE(mp_para_env_type), POINTER :: para_env
104 : TYPE(neighbor_list_iterator_p_type), &
105 41094 : DIMENSION(:), POINTER :: nl_iterator
106 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
107 41094 : POINTER :: sab_orb
108 41094 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
109 41094 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
110 41094 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
111 : TYPE(se_int_control_type) :: se_int_control
112 : TYPE(se_taper_type), POINTER :: se_taper
113 : TYPE(semi_empirical_control_type), POINTER :: se_control
114 41094 : TYPE(semi_empirical_p_type), DIMENSION(:), POINTER :: se_kind_list
115 : TYPE(semi_empirical_type), POINTER :: se_kind_a, se_kind_b
116 : TYPE(virial_type), POINTER :: virial
117 :
118 41094 : CALL timeset(routineN, handle)
119 :
120 41094 : NULLIFY (dft_control, cell, force, particle_set, se_control, se_taper)
121 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, cell=cell, se_taper=se_taper, &
122 41094 : para_env=para_env, virial=virial)
123 :
124 41094 : CALL initialize_se_taper(se_taper, exchange=.TRUE.)
125 41094 : se_control => dft_control%qs_control%se_control
126 41094 : anag = se_control%analytical_gradients
127 41094 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
128 41094 : nspins = dft_control%nspins
129 :
130 41094 : CPASSERT(ASSOCIATED(matrix_p))
131 41094 : CPASSERT(SIZE(ks_matrix) > 0)
132 :
133 : ! Identify proper integral screening (according user requests)
134 41094 : integral_screening = se_control%integral_screening
135 41094 : IF ((integral_screening == do_se_IS_kdso_d) .AND. (.NOT. se_control%force_kdsod_EX)) THEN
136 1462 : integral_screening = do_se_IS_kdso
137 : END IF
138 : CALL setup_se_int_control_type(se_int_control, shortrange=.FALSE., &
139 : do_ewald_r3=.FALSE., do_ewald_gks=.FALSE., integral_screening=integral_screening, &
140 41094 : max_multipole=do_multipole_none, pc_coulomb_int=.FALSE.)
141 :
142 : CALL get_qs_env(qs_env=qs_env, sab_orb=sab_orb, &
143 41094 : atomic_kind_set=atomic_kind_set, qs_kind_set=qs_kind_set)
144 :
145 41094 : nkind = SIZE(atomic_kind_set)
146 41094 : IF (calculate_forces) THEN
147 3002 : CALL get_qs_env(qs_env=qs_env, particle_set=particle_set, force=force)
148 3002 : delta = se_control%delta
149 3002 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
150 : END IF
151 :
152 347136 : ALLOCATE (se_defined(nkind), se_kind_list(nkind))
153 141666 : DO ikind = 1, nkind
154 100572 : CALL get_qs_kind(qs_kind_set(ikind), se_parameter=se_kind_a)
155 100572 : se_kind_list(ikind)%se_param => se_kind_a
156 100572 : CALL get_se_param(se_kind_a, defined=defined, natorb=natorb_a)
157 141666 : se_defined(ikind) = (defined .AND. natorb_a >= 1)
158 : END DO
159 :
160 41094 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb)
161 3821694 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0)
162 3780600 : CALL get_iterator_info(nl_iterator, ikind=ikind, jkind=jkind, iatom=iatom, jatom=jatom, r=rij)
163 3780600 : IF (.NOT. se_defined(ikind)) CYCLE
164 3780600 : IF (.NOT. se_defined(jkind)) CYCLE
165 3780600 : se_kind_a => se_kind_list(ikind)%se_param
166 3780600 : se_kind_b => se_kind_list(jkind)%se_param
167 :
168 3780600 : IF (iatom <= jatom) THEN
169 1941399 : irow = iatom
170 1941399 : icol = jatom
171 1941399 : switch = .FALSE.
172 : ELSE
173 1839201 : irow = jatom
174 1839201 : icol = iatom
175 1839201 : switch = .TRUE.
176 : END IF
177 : ! Retrieve blocks for KS and P
178 : CALL dbcsr_get_block_p(matrix=ks_matrix(1)%matrix, &
179 3780600 : row=irow, col=icol, BLOCK=ks_block_a, found=found)
180 3780600 : CPASSERT(ASSOCIATED(ks_block_a))
181 : CALL dbcsr_get_block_p(matrix=matrix_p(1)%matrix, &
182 3780600 : row=irow, col=icol, BLOCK=p_block_a, found=found)
183 3780600 : CPASSERT(ASSOCIATED(p_block_a))
184 3780600 : size_p_block_a(1) = SIZE(p_block_a, 1)
185 3780600 : size_p_block_a(2) = SIZE(p_block_a, 2)
186 43382323 : p_block_tot(1:size_p_block_a(1), 1:size_p_block_a(2)) = 2.0_dp*p_block_a
187 :
188 : ! Handle more configurations
189 3780600 : IF (nspins == 2) THEN
190 : CALL dbcsr_get_block_p(matrix=ks_matrix(2)%matrix, &
191 31652 : row=irow, col=icol, BLOCK=ks_block_b, found=found)
192 31652 : CPASSERT(ASSOCIATED(ks_block_b))
193 : CALL dbcsr_get_block_p(matrix=matrix_p(2)%matrix, &
194 31652 : row=irow, col=icol, BLOCK=p_block_b, found=found)
195 31652 : CPASSERT(ASSOCIATED(p_block_b))
196 31652 : check = (size_p_block_a(1) == SIZE(p_block_b, 1)) .AND. (size_p_block_a(2) == SIZE(p_block_b, 2))
197 0 : CPASSERT(check)
198 350981 : p_block_tot(1:SIZE(p_block_a, 1), 1:SIZE(p_block_a, 2)) = p_block_a + p_block_b
199 : END IF
200 :
201 15122400 : dr = DOT_PRODUCT(rij, rij)
202 3821694 : IF (iatom == jatom .AND. dr < rij_threshold) THEN
203 : ! Once center - Two electron Terms
204 193095 : IF (nspins == 1) THEN
205 188231 : CALL fock1_2el(se_kind_a, p_block_tot, p_block_a, ks_block_a, factor=0.5_dp)
206 4864 : ELSE IF (nspins == 2) THEN
207 4864 : CALL fock1_2el(se_kind_a, p_block_tot, p_block_a, ks_block_a, factor=1.0_dp)
208 4864 : CALL fock1_2el(se_kind_a, p_block_tot, p_block_b, ks_block_b, factor=1.0_dp)
209 : END IF
210 : ELSE
211 : ! Exchange Terms
212 3587505 : IF (nspins == 1) THEN
213 : CALL fock2E(se_kind_a, se_kind_b, rij, switch, size_p_block_a, p_block_a, ks_block_a, &
214 : factor=0.5_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper, &
215 3560717 : store_int_env=store_int_env)
216 26788 : ELSE IF (nspins == 2) THEN
217 : CALL fock2E(se_kind_a, se_kind_b, rij, switch, size_p_block_a, p_block_a, ks_block_a, &
218 : factor=1.0_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper, &
219 26788 : store_int_env=store_int_env)
220 :
221 : CALL fock2E(se_kind_a, se_kind_b, rij, switch, size_p_block_a, p_block_b, ks_block_b, &
222 : factor=1.0_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper, &
223 26788 : store_int_env=store_int_env)
224 : END IF
225 3587505 : IF (calculate_forces) THEN
226 172272 : atom_a = atom_of_kind(iatom)
227 172272 : atom_b = atom_of_kind(jatom)
228 172272 : force_ab = 0.0_dp
229 172272 : IF (nspins == 1) THEN
230 : CALL dfock2E(se_kind_a, se_kind_b, rij, switch, size_p_block_a, p_block_a, &
231 : factor=0.5_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper, force=force_ab, &
232 172111 : delta=delta)
233 161 : ELSE IF (nspins == 2) THEN
234 : CALL dfock2E(se_kind_a, se_kind_b, rij, switch, size_p_block_a, p_block_a, &
235 : factor=1.0_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper, force=force_ab, &
236 161 : delta=delta)
237 :
238 : CALL dfock2E(se_kind_a, se_kind_b, rij, switch, size_p_block_a, p_block_b, &
239 : factor=1.0_dp, anag=anag, se_int_control=se_int_control, se_taper=se_taper, force=force_ab, &
240 161 : delta=delta)
241 : END IF
242 172272 : IF (switch) THEN
243 88778 : force_ab(1) = -force_ab(1)
244 88778 : force_ab(2) = -force_ab(2)
245 88778 : force_ab(3) = -force_ab(3)
246 : END IF
247 172272 : IF (use_virial) THEN
248 0 : CALL virial_pair_force(virial%pv_virial, -1.0_dp, force_ab, rij)
249 : END IF
250 :
251 172272 : force(ikind)%rho_elec(1, atom_a) = force(ikind)%rho_elec(1, atom_a) - force_ab(1)
252 172272 : force(jkind)%rho_elec(1, atom_b) = force(jkind)%rho_elec(1, atom_b) + force_ab(1)
253 :
254 172272 : force(ikind)%rho_elec(2, atom_a) = force(ikind)%rho_elec(2, atom_a) - force_ab(2)
255 172272 : force(jkind)%rho_elec(2, atom_b) = force(jkind)%rho_elec(2, atom_b) + force_ab(2)
256 :
257 172272 : force(ikind)%rho_elec(3, atom_a) = force(ikind)%rho_elec(3, atom_a) - force_ab(3)
258 172272 : force(jkind)%rho_elec(3, atom_b) = force(jkind)%rho_elec(3, atom_b) + force_ab(3)
259 : END IF
260 : END IF
261 : END DO
262 41094 : CALL neighbor_list_iterator_release(nl_iterator)
263 :
264 41094 : DEALLOCATE (se_kind_list, se_defined)
265 :
266 41094 : CALL finalize_se_taper(se_taper)
267 :
268 41094 : CALL timestop(handle)
269 :
270 82188 : END SUBROUTINE build_fock_matrix_exchange
271 :
272 : END MODULE se_fock_matrix_exchange
273 :
|
