Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief routines that build the Kohn-Sham matrix contributions
10 : !> coming from local atomic densities
11 : ! **************************************************************************************************
12 : MODULE qs_ks_atom
13 :
14 : USE ao_util, ONLY: trace_r_AxB
15 : USE atomic_kind_types, ONLY: atomic_kind_type,&
16 : get_atomic_kind,&
17 : get_atomic_kind_set
18 : USE basis_set_types, ONLY: gto_basis_set_p_type,&
19 : gto_basis_set_type
20 : USE cp_array_utils, ONLY: cp_2d_r_p_type
21 : USE cp_control_types, ONLY: dft_control_type
22 : USE cp_dbcsr_api, ONLY: dbcsr_get_block_p,&
23 : dbcsr_p_type
24 : USE kinds, ONLY: dp,&
25 : int_8
26 : USE kpoint_types, ONLY: get_kpoint_info,&
27 : kpoint_type
28 : USE message_passing, ONLY: mp_para_env_type
29 : USE qs_environment_types, ONLY: get_qs_env,&
30 : qs_environment_type
31 : USE qs_force_types, ONLY: qs_force_type
32 : USE qs_kind_types, ONLY: get_qs_kind,&
33 : get_qs_kind_set,&
34 : qs_kind_type
35 : USE qs_neighbor_list_types, ONLY: get_iterator_task,&
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 : neighbor_list_task_type
42 : USE qs_nl_hash_table_types, ONLY: nl_hash_table_add,&
43 : nl_hash_table_create,&
44 : nl_hash_table_get_from_index,&
45 : nl_hash_table_is_null,&
46 : nl_hash_table_obj,&
47 : nl_hash_table_release,&
48 : nl_hash_table_status
49 : USE qs_oce_methods, ONLY: prj_gather
50 : USE qs_oce_types, ONLY: oce_matrix_type
51 : USE qs_rho_atom_types, ONLY: get_rho_atom,&
52 : rho_atom_coeff,&
53 : rho_atom_type
54 : USE sap_kind_types, ONLY: alist_post_align_blk,&
55 : alist_pre_align_blk,&
56 : alist_type,&
57 : get_alist
58 : USE util, ONLY: get_limit
59 : USE virial_methods, ONLY: virial_pair_force
60 : USE virial_types, ONLY: virial_type
61 :
62 : !$ USE OMP_LIB, ONLY: omp_get_max_threads, &
63 : !$ omp_get_thread_num, &
64 : !$ omp_lock_kind, &
65 : !$ omp_init_lock, omp_set_lock, &
66 : !$ omp_unset_lock, omp_destroy_lock
67 :
68 : #include "./base/base_uses.f90"
69 :
70 : IMPLICIT NONE
71 :
72 : PRIVATE
73 :
74 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_ks_atom'
75 :
76 : PUBLIC :: update_ks_atom
77 :
78 : CONTAINS
79 :
80 : ! **************************************************************************************************
81 : !> \brief The correction to the KS matrix due to the GAPW local terms to the hartree and
82 : !> XC contributions is here added. The correspondig forces contribution are also calculated
83 : !> if required. Each sparse-matrix block A-B is corrected by all the atomic contributions
84 : !> centered on atoms C for which the triplet A-C-B exists (they are close enough)
85 : !> To this end special lists are used
86 : !> \param qs_env qs environment, for the lists, the contraction coefficients and the
87 : !> pre calculated integrals of the potential with the atomic orbitals
88 : !> \param ksmat KS matrix, sparse matrix
89 : !> \param pmat density matrix, sparse matrix, needed only for the forces
90 : !> \param forces switch for the calculation of the forces on atoms
91 : !> \param tddft switch for TDDFT linear response
92 : !> \param rho_atom_external ...
93 : !> \param kind_set_external ...
94 : !> \param oce_external ...
95 : !> \param sab_external ...
96 : !> \param kscale ...
97 : !> \param kintegral ...
98 : !> \param kforce ...
99 : !> \param fscale ...
100 : !> \par History
101 : !> created [MI]
102 : !> the loop over the spins is done internally [03-05,MI]
103 : !> Rewrite using new OCE matrices [08.02.09, jhu]
104 : !> Add OpenMP [Apr 2016, EPCC]
105 : !> Allow for external kind_set, rho_atom_set, oce, sab 12.2019 (A. Bussy)
106 : ! **************************************************************************************************
107 41308 : SUBROUTINE update_ks_atom(qs_env, ksmat, pmat, forces, tddft, rho_atom_external, &
108 : kind_set_external, oce_external, sab_external, kscale, &
109 : kintegral, kforce, fscale)
110 :
111 : TYPE(qs_environment_type), POINTER :: qs_env
112 : TYPE(dbcsr_p_type), DIMENSION(*), INTENT(INOUT) :: ksmat, pmat
113 : LOGICAL, INTENT(IN) :: forces
114 : LOGICAL, INTENT(IN), OPTIONAL :: tddft
115 : TYPE(rho_atom_type), DIMENSION(:), OPTIONAL, &
116 : POINTER :: rho_atom_external
117 : TYPE(qs_kind_type), DIMENSION(:), OPTIONAL, &
118 : POINTER :: kind_set_external
119 : TYPE(oce_matrix_type), OPTIONAL, POINTER :: oce_external
120 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
121 : OPTIONAL, POINTER :: sab_external
122 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: kscale, kintegral, kforce
123 : REAL(KIND=dp), DIMENSION(2), INTENT(IN), OPTIONAL :: fscale
124 :
125 : CHARACTER(len=*), PARAMETER :: routineN = 'update_ks_atom'
126 :
127 : INTEGER :: bo(2), handle, ia_kind, iac, iat, iatom, ibc, ikind, img, ip, ispin, ja_kind, &
128 : jatom, jkind, ka_kind, kac, katom, kbc, kkind, max_gau, max_nsgf, n_cont_a, n_cont_b, &
129 : nat, natom, nimages, nkind, nl_table_num_slots, nsoctot, nspins, slot_num
130 : INTEGER(KIND=int_8) :: nl_table_key
131 41308 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
132 : INTEGER, DIMENSION(3) :: cell_b
133 41308 : INTEGER, DIMENSION(:), POINTER :: atom_list, list_a, list_b
134 41308 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
135 : LOGICAL :: dista, distb, do_forces, found, &
136 : is_entry_null, is_task_valid, &
137 : my_tddft, paw_atom, use_virial
138 41308 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: has_intac, paw_kind
139 41308 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: a_matrix, p_matrix
140 : REAL(dp), DIMENSION(3) :: rac, rbc
141 : REAL(dp), DIMENSION(3, 3) :: force_tmp
142 : REAL(kind=dp) :: eps_cpc, factor1, factor2
143 41308 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: C_int_h, C_int_s, coc
144 : REAL(KIND=dp), DIMENSION(2) :: force_fac
145 : REAL(KIND=dp), DIMENSION(3, 3) :: pv_virial_thread
146 41308 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: C_coeff_hh_a, C_coeff_hh_b, &
147 41308 : C_coeff_ss_a, C_coeff_ss_b
148 : TYPE(alist_type), POINTER :: alist_ac, alist_bc
149 41308 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
150 41308 : TYPE(cp_2d_r_p_type), DIMENSION(:), POINTER :: mat_h, mat_p
151 : TYPE(dft_control_type), POINTER :: dft_control
152 41308 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
153 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
154 : TYPE(kpoint_type), POINTER :: kpoints
155 : TYPE(mp_para_env_type), POINTER :: para_env
156 : TYPE(neighbor_list_iterator_p_type), &
157 41308 : DIMENSION(:), POINTER :: nl_iterator
158 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
159 41308 : POINTER :: my_sab
160 : TYPE(neighbor_list_task_type), POINTER :: next_task, task
161 : TYPE(nl_hash_table_obj) :: nl_hash_table
162 : TYPE(oce_matrix_type), POINTER :: my_oce
163 41308 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
164 41308 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: my_kind_set
165 41308 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: int_local_h, int_local_s
166 41308 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: my_rho_atom
167 : TYPE(rho_atom_type), POINTER :: rho_at
168 : TYPE(virial_type), POINTER :: virial
169 :
170 41308 : !$ INTEGER(kind=omp_lock_kind), ALLOCATABLE, DIMENSION(:) :: locks
171 : !$ INTEGER :: lock_num
172 :
173 41308 : CALL timeset(routineN, handle)
174 :
175 41308 : NULLIFY (my_kind_set, atomic_kind_set, force, my_oce, para_env, my_rho_atom, my_sab)
176 41308 : NULLIFY (mat_h, mat_p, dft_control)
177 :
178 : CALL get_qs_env(qs_env=qs_env, &
179 : qs_kind_set=my_kind_set, &
180 : atomic_kind_set=atomic_kind_set, &
181 : force=force, &
182 : oce=my_oce, &
183 : para_env=para_env, &
184 : rho_atom_set=my_rho_atom, &
185 : virial=virial, &
186 : sab_orb=my_sab, &
187 41308 : dft_control=dft_control)
188 :
189 41308 : nspins = dft_control%nspins
190 41308 : nimages = dft_control%nimages
191 :
192 41308 : factor1 = 1.0_dp
193 41308 : factor2 = 1.0_dp
194 :
195 : !deal with externals
196 41308 : my_tddft = .FALSE.
197 41308 : IF (PRESENT(tddft)) my_tddft = tddft
198 11766 : IF (my_tddft) THEN
199 6028 : IF (nspins == 1) factor1 = 2.0_dp
200 6028 : CPASSERT(nimages == 1)
201 : END IF
202 41308 : IF (PRESENT(kscale)) THEN
203 256 : factor1 = factor1*kscale
204 256 : factor2 = factor2*kscale
205 : END IF
206 41308 : IF (PRESENT(kintegral)) factor1 = kintegral
207 41308 : IF (PRESENT(kforce)) factor2 = kforce
208 123924 : force_fac = 1.0_dp
209 41308 : IF (PRESENT(fscale)) force_fac(:) = fscale(:)
210 :
211 41308 : IF (PRESENT(rho_atom_external)) my_rho_atom => rho_atom_external
212 41308 : IF (PRESENT(kind_set_external)) my_kind_set => kind_set_external
213 41308 : IF (PRESENT(oce_external)) my_oce => oce_external
214 41308 : IF (PRESENT(sab_external)) my_sab => sab_external
215 :
216 : ! kpoint images
217 41308 : NULLIFY (cell_to_index)
218 41308 : IF (nimages > 1) THEN
219 628 : CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
220 628 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
221 : END IF
222 :
223 41308 : eps_cpc = dft_control%qs_control%gapw_control%eps_cpc
224 :
225 41308 : CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
226 41308 : CALL get_qs_kind_set(my_kind_set, maxsgf=max_nsgf, maxgtops=max_gau, basis_type="GAPW_1C")
227 :
228 41308 : IF (forces) THEN
229 1992 : CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind)
230 1992 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
231 : ELSE
232 : use_virial = .FALSE.
233 : END IF
234 :
235 41308 : pv_virial_thread(:, :) = 0.0_dp ! always initialize to avoid floating point exceptions in OMP REDUCTION
236 :
237 41308 : nkind = SIZE(my_kind_set, 1)
238 : ! Collect the local potential contributions from all the processors
239 : ASSOCIATE (mepos => para_env%mepos, num_pe => para_env%num_pe)
240 121456 : DO ikind = 1, nkind
241 80148 : NULLIFY (atom_list)
242 80148 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=nat)
243 80148 : CALL get_qs_kind(my_kind_set(ikind), paw_atom=paw_atom)
244 121456 : IF (paw_atom) THEN
245 : ! gather the atomic block integrals in a more compressed format
246 71892 : bo = get_limit(nat, num_pe, mepos)
247 127017 : DO iat = bo(1), bo(2)
248 55125 : iatom = atom_list(iat)
249 189542 : DO ispin = 1, nspins
250 : CALL prj_gather(my_rho_atom(iatom)%ga_Vlocal_gb_h(ispin)%r_coef, &
251 62525 : my_rho_atom(iatom)%int_scr_h(ispin)%r_coef, my_kind_set(ikind))
252 : CALL prj_gather(my_rho_atom(iatom)%ga_Vlocal_gb_s(ispin)%r_coef, &
253 117650 : my_rho_atom(iatom)%int_scr_s(ispin)%r_coef, my_kind_set(ikind))
254 : END DO
255 : END DO
256 : ! broadcast the CPC arrays to all processors (replicated data)
257 215676 : DO ip = 0, num_pe - 1
258 143784 : bo = get_limit(nat, num_pe, ip)
259 325926 : DO iat = bo(1), bo(2)
260 110250 : iatom = atom_list(iat)
261 379084 : DO ispin = 1, nspins
262 120728050 : CALL para_env%bcast(my_rho_atom(iatom)%int_scr_h(ispin)%r_coef, ip)
263 120838300 : CALL para_env%bcast(my_rho_atom(iatom)%int_scr_s(ispin)%r_coef, ip)
264 : END DO
265 : END DO
266 : END DO
267 : END IF
268 : END DO
269 : END ASSOCIATE
270 :
271 204072 : ALLOCATE (basis_set_list(nkind))
272 206540 : ALLOCATE (paw_kind(nkind), has_intac(nkind*nkind))
273 41308 : paw_kind(:) = .FALSE.
274 41308 : has_intac(:) = .FALSE.
275 121456 : DO ikind = 1, nkind
276 80148 : CALL get_qs_kind(my_kind_set(ikind), basis_set=basis_set_a, paw_atom=paw_kind(ikind))
277 121456 : IF (ASSOCIATED(basis_set_a)) THEN
278 80148 : basis_set_list(ikind)%gto_basis_set => basis_set_a
279 : ELSE
280 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
281 : END IF
282 : END DO
283 211912 : DO ikind = 1, nkind*nkind
284 211912 : has_intac(ikind) = ASSOCIATED(my_oce%intac(ikind)%alist)
285 : END DO
286 :
287 : ! build the hash table in serial...
288 : ! ... creation ...
289 41308 : CALL neighbor_list_iterator_create(nl_iterator, my_sab)
290 41308 : nl_table_num_slots = natom*natom/2 ! this is probably not optimal, but it seems a good start
291 41308 : CALL nl_hash_table_create(nl_hash_table, nmax=nl_table_num_slots)
292 : ! ... and population
293 1147803 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0) ! build hash table in serial, so don't pass mepos
294 9958455 : ALLOCATE (task) ! They must be deallocated before the nl_hash_table is released
295 1106495 : CALL get_iterator_task(nl_iterator, task) ! build hash table in serial, so don't pass mepos
296 : ! tasks with the same key access the same blocks of H & P
297 1106495 : IF (task%iatom <= task%jatom) THEN
298 665606 : nl_table_key = natom*task%iatom + task%jatom
299 : ELSE
300 440889 : nl_table_key = natom*task%jatom + task%iatom
301 : END IF
302 1106495 : CALL nl_hash_table_add(nl_hash_table, nl_table_key, task)
303 : END DO
304 41308 : CALL neighbor_list_iterator_release(nl_iterator)
305 :
306 : ! Get the total number of (possibly empty) entries in the table
307 41308 : CALL nl_hash_table_status(nl_hash_table, nmax=nl_table_num_slots)
308 :
309 : !$OMP PARALLEL DEFAULT(NONE) &
310 : !$OMP SHARED(nl_table_num_slots, nl_hash_table &
311 : !$OMP , max_gau, max_nsgf, nspins, forces &
312 : !$OMP , basis_set_list, nimages, cell_to_index &
313 : !$OMP , ksmat, pmat, natom, nkind, my_kind_set, my_oce &
314 : !$OMP , my_rho_atom, factor1, factor2, use_virial &
315 : !$OMP , atom_of_kind, force, locks, force_fac &
316 : !$OMP , paw_kind, has_intac &
317 : !$OMP ) &
318 : !$OMP PRIVATE( slot_num, is_entry_null, TASK, is_task_valid &
319 : !$OMP , C_int_h, C_int_s, coc, a_matrix, p_matrix &
320 : !$OMP , mat_h, mat_p &
321 : !$OMP , int_local_h, int_local_s &
322 : !$OMP , ikind, jkind, iatom, jatom, cell_b &
323 : !$OMP , basis_set_a, basis_set_b, img &
324 : !$OMP , found, next_task &
325 : !$OMP , kkind, lock_num &
326 : !$OMP , iac, alist_ac, kac, n_cont_a, list_a &
327 : !$OMP , ibc, alist_bc, kbc, n_cont_b, list_b &
328 : !$OMP , katom, rho_at, nsoctot &
329 : !$OMP , C_coeff_hh_a, C_coeff_ss_a, dista, rac &
330 : !$OMP , C_coeff_hh_b, C_coeff_ss_b, distb, rbc &
331 : !$OMP , ia_kind, ja_kind, ka_kind, force_tmp, do_forces &
332 : !$OMP ) &
333 : !$OMP REDUCTION(+:pv_virial_thread &
334 41308 : !$OMP )
335 :
336 : ALLOCATE (C_int_h(max_gau*max_nsgf), C_int_s(max_gau*max_nsgf), coc(max_gau*max_gau), &
337 : a_matrix(max_gau, max_gau), p_matrix(max_nsgf, max_nsgf))
338 :
339 : ALLOCATE (mat_h(nspins), mat_p(nspins))
340 : DO ispin = 1, nspins
341 : NULLIFY (mat_h(ispin)%array, mat_p(ispin)%array)
342 : END DO
343 :
344 : IF (forces) THEN
345 : !$OMP SINGLE
346 : !$ ALLOCATE (locks(natom*nkind))
347 : !$OMP END SINGLE
348 :
349 : !$OMP DO
350 : !$ do lock_num = 1, natom*nkind
351 : !$ call omp_init_lock(locks(lock_num))
352 : !$ end do
353 : !$OMP END DO
354 : END IF
355 :
356 : ! Dynamic schedule to take account of the fact that some slots may be empty
357 : ! or contain 1 or more tasks
358 : !$OMP DO SCHEDULE(DYNAMIC,5)
359 : DO slot_num = 1, nl_table_num_slots
360 : CALL nl_hash_table_is_null(nl_hash_table, slot_num, is_entry_null)
361 :
362 : IF (.NOT. is_entry_null) THEN
363 : CALL nl_hash_table_get_from_index(nl_hash_table, slot_num, task)
364 :
365 : is_task_valid = ASSOCIATED(task)
366 : DO WHILE (is_task_valid)
367 :
368 : ikind = task%ikind
369 : jkind = task%jkind
370 : iatom = task%iatom
371 : jatom = task%jatom
372 : cell_b(:) = task%cell(:)
373 :
374 : basis_set_a => basis_set_list(ikind)%gto_basis_set
375 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
376 : basis_set_b => basis_set_list(jkind)%gto_basis_set
377 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
378 :
379 : IF (nimages > 1) THEN
380 : img = cell_to_index(cell_b(1), cell_b(2), cell_b(3))
381 : CPASSERT(img > 0)
382 : ELSE
383 : img = 1
384 : END IF
385 :
386 : DO ispin = 1, nspins
387 : NULLIFY (mat_h(ispin)%array, mat_p(ispin)%array)
388 :
389 : found = .FALSE.
390 : IF (iatom <= jatom) THEN
391 : CALL dbcsr_get_block_p(matrix=ksmat(nspins*(img - 1) + ispin)%matrix, &
392 : row=iatom, col=jatom, &
393 : BLOCK=mat_h(ispin)%array, found=found)
394 : ELSE
395 : CALL dbcsr_get_block_p(matrix=ksmat(nspins*(img - 1) + ispin)%matrix, &
396 : row=jatom, col=iatom, &
397 : BLOCK=mat_h(ispin)%array, found=found)
398 : END IF
399 : CPASSERT(found)
400 :
401 : IF (forces) THEN
402 : found = .FALSE.
403 : IF (iatom <= jatom) THEN
404 : CALL dbcsr_get_block_p(matrix=pmat(nspins*(img - 1) + ispin)%matrix, &
405 : row=iatom, col=jatom, &
406 : BLOCK=mat_p(ispin)%array, found=found)
407 : ELSE
408 : CALL dbcsr_get_block_p(matrix=pmat(nspins*(img - 1) + ispin)%matrix, &
409 : row=jatom, col=iatom, &
410 : BLOCK=mat_p(ispin)%array, found=found)
411 : END IF
412 : CPASSERT(found)
413 : END IF
414 : END DO
415 :
416 : DO kkind = 1, nkind
417 : IF (.NOT. paw_kind(kkind)) CYCLE
418 :
419 : iac = ikind + nkind*(kkind - 1)
420 : ibc = jkind + nkind*(kkind - 1)
421 : IF (.NOT. has_intac(iac)) CYCLE
422 : IF (.NOT. has_intac(ibc)) CYCLE
423 :
424 : CALL get_alist(my_oce%intac(iac), alist_ac, iatom)
425 : CALL get_alist(my_oce%intac(ibc), alist_bc, jatom)
426 : IF (.NOT. ASSOCIATED(alist_ac)) CYCLE
427 : IF (.NOT. ASSOCIATED(alist_bc)) CYCLE
428 :
429 : DO kac = 1, alist_ac%nclist
430 : DO kbc = 1, alist_bc%nclist
431 : IF (alist_ac%clist(kac)%catom /= alist_bc%clist(kbc)%catom) CYCLE
432 :
433 : IF (ALL(cell_b + alist_bc%clist(kbc)%cell - alist_ac%clist(kac)%cell == 0)) THEN
434 : n_cont_a = alist_ac%clist(kac)%nsgf_cnt
435 : n_cont_b = alist_bc%clist(kbc)%nsgf_cnt
436 : IF (n_cont_a == 0 .OR. n_cont_b == 0) CYCLE
437 :
438 : list_a => alist_ac%clist(kac)%sgf_list
439 : list_b => alist_bc%clist(kbc)%sgf_list
440 :
441 : katom = alist_ac%clist(kac)%catom
442 :
443 : IF (iatom == katom .AND. ALL(alist_ac%clist(kac)%cell == 0)) THEN
444 : C_coeff_hh_a => alist_ac%clist(kac)%achint
445 : C_coeff_ss_a => alist_ac%clist(kac)%acint
446 : dista = .FALSE.
447 : ELSE
448 : C_coeff_hh_a => alist_ac%clist(kac)%acint
449 : C_coeff_ss_a => alist_ac%clist(kac)%acint
450 : dista = .TRUE.
451 : END IF
452 :
453 : IF (jatom == katom .AND. ALL(alist_bc%clist(kbc)%cell == 0)) THEN
454 : C_coeff_hh_b => alist_bc%clist(kbc)%achint
455 : C_coeff_ss_b => alist_bc%clist(kbc)%acint
456 : distb = .FALSE.
457 : ELSE
458 : C_coeff_hh_b => alist_bc%clist(kbc)%acint
459 : C_coeff_ss_b => alist_bc%clist(kbc)%acint
460 : distb = .TRUE.
461 : END IF
462 :
463 : rho_at => my_rho_atom(katom)
464 : nsoctot = SIZE(C_coeff_ss_a, 2)
465 :
466 : CALL get_rho_atom(rho_atom=rho_at, int_scr_h=int_local_h, int_scr_s=int_local_s)
467 : do_forces = forces .AND. (dista .OR. distb)
468 : ! Use the stored block ordering for both KS and force contractions.
469 : IF (iatom <= jatom) THEN
470 : CALL add_vhxca(mat_h, mat_p, C_coeff_hh_a, C_coeff_hh_b, C_coeff_ss_a, C_coeff_ss_b, &
471 : int_local_h, int_local_s, nspins, iatom, jatom, nsoctot, factor1, &
472 : list_a, n_cont_a, list_b, n_cont_b, dista, distb, &
473 : C_int_h, C_int_s, max_nsgf, coc, a_matrix, p_matrix, &
474 : do_forces, force_tmp, force_fac)
475 : ELSE
476 : CALL add_vhxca(mat_h, mat_p, C_coeff_hh_b, C_coeff_hh_a, C_coeff_ss_b, C_coeff_ss_a, &
477 : int_local_h, int_local_s, nspins, jatom, iatom, nsoctot, factor1, &
478 : list_b, n_cont_b, list_a, n_cont_a, distb, dista, &
479 : C_int_h, C_int_s, max_nsgf, coc, a_matrix, p_matrix, &
480 : do_forces, force_tmp, force_fac)
481 : ! Restore the original task ordering of the two endpoint forces.
482 : force_tmp(:, 1:2) = force_tmp(:, 2:1:-1)
483 : END IF
484 :
485 : IF (do_forces) THEN
486 : IF (use_virial) THEN
487 : rac = alist_ac%clist(kac)%rac
488 : rbc = alist_bc%clist(kbc)%rac
489 : END IF
490 : ia_kind = atom_of_kind(iatom)
491 : ja_kind = atom_of_kind(jatom)
492 : ka_kind = atom_of_kind(katom)
493 : force_tmp = factor2*force_tmp
494 : !$ CALL omp_set_lock(locks((ka_kind - 1)*nkind + kkind))
495 : force(kkind)%vhxc_atom(1:3, ka_kind) = force(kkind)%vhxc_atom(1:3, ka_kind) + force_tmp(1:3, 3)
496 : !$ CALL omp_unset_lock(locks((ka_kind - 1)*nkind + kkind))
497 : !$ CALL omp_set_lock(locks((ia_kind - 1)*nkind + ikind))
498 : force(ikind)%vhxc_atom(1:3, ia_kind) = force(ikind)%vhxc_atom(1:3, ia_kind) + force_tmp(1:3, 1)
499 : !$ CALL omp_unset_lock(locks((ia_kind - 1)*nkind + ikind))
500 : !$ CALL omp_set_lock(locks((ja_kind - 1)*nkind + jkind))
501 : force(jkind)%vhxc_atom(1:3, ja_kind) = force(jkind)%vhxc_atom(1:3, ja_kind) + force_tmp(1:3, 2)
502 : !$ CALL omp_unset_lock(locks((ja_kind - 1)*nkind + jkind))
503 : IF (use_virial) THEN
504 : CALL virial_pair_force(pv_virial_thread, 1._dp, force_tmp(1:3, 1), rac)
505 : CALL virial_pair_force(pv_virial_thread, 1._dp, force_tmp(1:3, 2), rbc)
506 : END IF
507 : END IF
508 : EXIT ! search loop over jatom-katom list
509 : END IF
510 : END DO ! kbc
511 : END DO ! kac
512 :
513 : END DO ! kkind
514 :
515 : next_task => task%next
516 : ! We are done with this task, we can deallocate it
517 : DEALLOCATE (task)
518 : is_task_valid = ASSOCIATED(next_task)
519 : IF (is_task_valid) task => next_task
520 :
521 : END DO
522 :
523 : ELSE
524 : ! NO KEY/VALUE
525 : END IF
526 : END DO
527 : !$OMP END DO
528 :
529 : DO ispin = 1, nspins
530 : NULLIFY (mat_h(ispin)%array, mat_p(ispin)%array)
531 : END DO
532 : DEALLOCATE (mat_h, mat_p, C_int_h, C_int_s, a_matrix, p_matrix, coc)
533 :
534 : IF (forces) THEN
535 : ! Implicit barrier at end of OMP DO, so locks can be freed
536 : !$OMP DO
537 : !$ DO lock_num = 1, natom*nkind
538 : !$ call omp_destroy_lock(locks(lock_num))
539 : !$ END DO
540 : !$OMP END DO
541 :
542 : !$OMP SINGLE
543 : !$ DEALLOCATE (locks)
544 : !$OMP END SINGLE NOWAIT
545 : END IF
546 :
547 : !$OMP END PARALLEL
548 :
549 41308 : IF (use_virial) THEN
550 1742 : virial%pv_gapw(1:3, 1:3) = virial%pv_gapw(1:3, 1:3) + pv_virial_thread(1:3, 1:3)
551 1742 : virial%pv_virial(1:3, 1:3) = virial%pv_virial(1:3, 1:3) + pv_virial_thread(1:3, 1:3)
552 : END IF
553 :
554 41308 : CALL nl_hash_table_release(nl_hash_table)
555 :
556 41308 : DEALLOCATE (basis_set_list, paw_kind, has_intac)
557 :
558 41308 : CALL timestop(handle)
559 :
560 165232 : END SUBROUTINE update_ks_atom
561 :
562 : ! **************************************************************************************************
563 : !> \brief Add one-center hard-minus-soft KS and coefficient-derivative force contributions.
564 : !> Real symmetric one-center potentials allow the same unscaled C*V intermediates
565 : !> to be used for KS assembly and forces within each spin.
566 : !> \param mat_h KS blocks in stored (A <= B) order
567 : !> \param mat_p density blocks in the same order, used only for forces
568 : !> \param C_hh_a hard coefficients and derivatives on A
569 : !> \param C_hh_b hard coefficients and derivatives on B
570 : !> \param C_ss_a soft coefficients and derivatives on A
571 : !> \param C_ss_b soft coefficients and derivatives on B
572 : !> \param int_local_h hard one-center potentials by spin
573 : !> \param int_local_s soft one-center potentials by spin
574 : !> \param nspins number of spins
575 : !> \param ia atom A
576 : !> \param ja atom B
577 : !> \param nsp size of the one-center basis
578 : !> \param factor KS scaling, kept separate from the force scaling
579 : !> \param lista contracted basis indices on A
580 : !> \param nconta number of contracted functions on A
581 : !> \param listb contracted basis indices on B
582 : !> \param ncontb number of contracted functions on B
583 : !> \param dista true if A is offsite relative to the expansion center
584 : !> \param distb true if B is offsite relative to the expansion center
585 : !> \param work_h workspace for unscaled potential-weighted coefficients
586 : !> \param work_s workspace for soft or density-weighted coefficients
587 : !> \param ldwork leading dimension of the coefficient workspaces
588 : !> \param int_diff workspace for the hard-minus-soft potential
589 : !> \param a_matrix workspace for the contracted KS block
590 : !> \param p_matrix workspace for the aligned density block
591 : !> \param do_forces whether nonzero coefficient-derivative contributions are needed
592 : !> \param force forces on A, B, and the expansion center
593 : !> \param force_scaling force scaling by spin
594 : ! **************************************************************************************************
595 16565398 : SUBROUTINE add_vhxca(mat_h, mat_p, C_hh_a, C_hh_b, C_ss_a, C_ss_b, &
596 8282699 : int_local_h, int_local_s, nspins, ia, ja, nsp, factor, &
597 8282699 : lista, nconta, listb, ncontb, dista, distb, &
598 24848097 : work_h, work_s, ldwork, int_diff, a_matrix, p_matrix, &
599 : do_forces, force, force_scaling)
600 : TYPE(cp_2d_r_p_type), DIMENSION(:), INTENT(IN) :: mat_h, mat_p
601 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: C_hh_a, C_hh_b, C_ss_a, C_ss_b
602 : TYPE(rho_atom_coeff), DIMENSION(:), INTENT(IN) :: int_local_h, int_local_s
603 : INTEGER, INTENT(IN) :: nspins, ia, ja, nsp
604 : REAL(KIND=dp), INTENT(IN) :: factor
605 : INTEGER, DIMENSION(:), INTENT(IN) :: lista
606 : INTEGER, INTENT(IN) :: nconta
607 : INTEGER, DIMENSION(:), INTENT(IN) :: listb
608 : INTEGER, INTENT(IN) :: ncontb
609 : LOGICAL, INTENT(IN) :: dista, distb
610 : INTEGER, INTENT(IN) :: ldwork
611 : REAL(KIND=dp), INTENT(OUT) :: work_s(ldwork, nsp), &
612 : work_h(ldwork, nsp), int_diff(nsp, nsp)
613 : REAL(dp), DIMENSION(:, :), INTENT(OUT) :: a_matrix, p_matrix
614 : LOGICAL, INTENT(IN) :: do_forces
615 : REAL(dp), DIMENSION(3, 3), INTENT(OUT) :: force
616 : REAL(KIND=dp), DIMENSION(2), INTENT(IN) :: force_scaling
617 :
618 : INTEGER :: dir, ispin
619 8282699 : REAL(dp), DIMENSION(:, :), POINTER :: h_block, int_hard, int_soft, p_block
620 : REAL(KIND=dp) :: ieqj, trace
621 :
622 8282699 : force(:, :) = 0.0_dp
623 : ! Factor 2 because only half of the pairs with ia /= ja are considered.
624 8282699 : ieqj = 2.0_dp
625 2448706 : IF (ia == ja) ieqj = 1.0_dp
626 :
627 16913371 : DO ispin = 1, nspins
628 8630672 : h_block => mat_h(ispin)%array
629 8630672 : int_hard => int_local_h(ispin)%r_coef
630 8630672 : int_soft => int_local_s(ispin)%r_coef
631 :
632 4984403180 : IF (dista .AND. distb) int_diff(:, :) = int_hard(1:nsp, 1:nsp) - int_soft(1:nsp, 1:nsp)
633 :
634 8630672 : IF (distb .AND. .NOT. dista) THEN
635 : ! B is offsite: H_ab = (C_a^h V_h - C_a^s V_s) C_b^T.
636 : CALL DGEMM('N', 'N', nconta, nsp, nsp, 1.0_dp, C_hh_a(:, :, 1), SIZE(C_hh_a, 1), &
637 213182 : int_hard, SIZE(int_hard, 1), 0.0_dp, work_h, ldwork)
638 : CALL DGEMM('N', 'N', nconta, nsp, nsp, -1.0_dp, C_ss_a(:, :, 1), SIZE(C_ss_a, 1), &
639 213182 : int_soft, SIZE(int_soft, 1), 1.0_dp, work_h, ldwork)
640 : CALL DGEMM('N', 'T', nconta, ncontb, nsp, factor, work_h, ldwork, &
641 213182 : C_hh_b(:, :, 1), SIZE(C_hh_b, 1), 0.0_dp, a_matrix, SIZE(a_matrix, 1))
642 : ELSE
643 8417490 : IF (distb) THEN
644 : ! Both ends are offsite: C_h = C_s, including their derivatives.
645 : CALL DGEMM('N', 'N', ncontb, nsp, nsp, 1.0_dp, C_hh_b(:, :, 1), SIZE(C_hh_b, 1), &
646 8170750 : int_diff, nsp, 0.0_dp, work_h, ldwork)
647 : ELSE
648 : CALL DGEMM('N', 'N', ncontb, nsp, nsp, 1.0_dp, C_hh_b(:, :, 1), SIZE(C_hh_b, 1), &
649 246740 : int_hard, SIZE(int_hard, 1), 0.0_dp, work_h, ldwork)
650 246740 : IF (dista) THEN
651 : CALL DGEMM('N', 'N', ncontb, nsp, nsp, -1.0_dp, C_ss_b(:, :, 1), SIZE(C_ss_b, 1), &
652 184215 : int_soft, SIZE(int_soft, 1), 1.0_dp, work_h, ldwork)
653 : ELSE
654 : CALL DGEMM('N', 'N', ncontb, nsp, nsp, 1.0_dp, C_ss_b(:, :, 1), SIZE(C_ss_b, 1), &
655 62525 : int_soft, SIZE(int_soft, 1), 0.0_dp, work_s, ldwork)
656 : END IF
657 : END IF
658 : CALL DGEMM('N', 'T', nconta, ncontb, nsp, factor, C_hh_a(:, :, 1), SIZE(C_hh_a, 1), &
659 8417490 : work_h, ldwork, 0.0_dp, a_matrix, SIZE(a_matrix, 1))
660 8417490 : IF (.NOT. dista) THEN
661 : CALL DGEMM('N', 'T', nconta, ncontb, nsp, -factor, C_ss_a(:, :, 1), SIZE(C_ss_a, 1), &
662 62525 : work_s, ldwork, 1.0_dp, a_matrix, SIZE(a_matrix, 1))
663 : END IF
664 : END IF
665 :
666 : CALL alist_post_align_blk(a_matrix, SIZE(a_matrix, 1), h_block, SIZE(h_block, 1), &
667 8630672 : lista, nconta, listb, ncontb)
668 :
669 8630672 : IF (.NOT. do_forces) CYCLE
670 779039 : p_block => mat_p(ispin)%array
671 : CALL alist_pre_align_blk(p_block, SIZE(p_block, 1), p_matrix, SIZE(p_matrix, 1), &
672 779039 : lista, nconta, listb, ncontb)
673 :
674 : ! Onsite coefficient derivatives are zero (see build_oce_matrices).
675 : ! The offsite derivative is common to hard and soft coefficients.
676 779039 : IF (dista) THEN
677 : ! Reuse K_b = C_b^h V_h - C_b^s V_s from KS assembly.
678 : CALL DGEMM('N', 'N', nconta, nsp, ncontb, 1.0_dp, p_matrix, SIZE(p_matrix, 1), &
679 766951 : work_h, ldwork, 0.0_dp, work_s, ldwork)
680 : ! <(P C_b)^T dC_a, V>_F = <dC_a, P (C_b V)>_F.
681 3067804 : DO dir = 2, 4
682 2300853 : trace = trace_r_AxB(work_s, ldwork, C_hh_a(:, :, dir), SIZE(C_hh_a, 1), nconta, nsp)
683 2300853 : trace = ieqj*trace*force_scaling(ispin)
684 2300853 : force(dir - 1, 3) = force(dir - 1, 3) + trace
685 3067804 : force(dir - 1, 1) = force(dir - 1, 1) - trace
686 : END DO
687 : END IF
688 :
689 9061738 : IF (distb) THEN
690 766657 : IF (dista) THEN
691 : ! A force used K_b; now form K_a = C_a (V_h - V_s).
692 : CALL DGEMM('N', 'N', nconta, nsp, nsp, 1.0_dp, C_hh_a(:, :, 1), SIZE(C_hh_a, 1), &
693 754569 : int_diff, nsp, 0.0_dp, work_h, ldwork)
694 : END IF
695 : ! For an onsite A, K_a is already available from KS assembly.
696 : CALL DGEMM('T', 'N', ncontb, nsp, nconta, 1.0_dp, p_matrix, SIZE(p_matrix, 1), &
697 766657 : work_h, ldwork, 0.0_dp, work_s, ldwork)
698 3066628 : DO dir = 2, 4
699 2299971 : trace = trace_r_AxB(work_s, ldwork, C_hh_b(:, :, dir), SIZE(C_hh_b, 1), ncontb, nsp)
700 2299971 : trace = ieqj*trace*force_scaling(ispin)
701 2299971 : force(dir - 1, 3) = force(dir - 1, 3) + trace
702 3066628 : force(dir - 1, 2) = force(dir - 1, 2) - trace
703 : END DO
704 : END IF
705 : END DO
706 8282699 : END SUBROUTINE add_vhxca
707 :
708 : END MODULE qs_ks_atom
|
