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 35274 : 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, ldCPC, max_gau, max_nsgf, n_cont_a, &
129 : n_cont_b, nat, natom, nimages, nkind, nl_table_num_slots, nsoctot, nspins, slot_num
130 : INTEGER(KIND=int_8) :: nl_table_key
131 35274 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
132 : INTEGER, DIMENSION(3) :: cell_b
133 35274 : INTEGER, DIMENSION(:), POINTER :: atom_list, list_a, list_b
134 35274 : INTEGER, DIMENSION(:, :, :), POINTER :: cell_to_index
135 : LOGICAL :: dista, distb, found, is_entry_null, &
136 : is_task_valid, my_tddft, paw_atom, &
137 : use_virial
138 35274 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: has_intac, paw_kind
139 35274 : 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 35274 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: C_int_h, C_int_s, coc
144 35274 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: dCPC_h, dCPC_s
145 35274 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: PC_h, PC_s
146 : REAL(KIND=dp), DIMENSION(2) :: force_fac
147 : REAL(KIND=dp), DIMENSION(3, 3) :: pv_virial_thread
148 35274 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: C_coeff_hh_a, C_coeff_hh_b, &
149 35274 : C_coeff_ss_a, C_coeff_ss_b
150 : TYPE(alist_type), POINTER :: alist_ac, alist_bc
151 35274 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
152 35274 : TYPE(cp_2d_r_p_type), DIMENSION(:), POINTER :: mat_h, mat_p
153 : TYPE(dft_control_type), POINTER :: dft_control
154 35274 : TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
155 : TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
156 : TYPE(kpoint_type), POINTER :: kpoints
157 : TYPE(mp_para_env_type), POINTER :: para_env
158 : TYPE(neighbor_list_iterator_p_type), &
159 35274 : DIMENSION(:), POINTER :: nl_iterator
160 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
161 35274 : POINTER :: my_sab
162 : TYPE(neighbor_list_task_type), POINTER :: next_task, task
163 : TYPE(nl_hash_table_obj) :: nl_hash_table
164 : TYPE(oce_matrix_type), POINTER :: my_oce
165 35274 : TYPE(qs_force_type), DIMENSION(:), POINTER :: force
166 35274 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: my_kind_set
167 35274 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: int_local_h, int_local_s
168 35274 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: my_rho_atom
169 : TYPE(rho_atom_type), POINTER :: rho_at
170 : TYPE(virial_type), POINTER :: virial
171 :
172 35274 : !$ INTEGER(kind=omp_lock_kind), ALLOCATABLE, DIMENSION(:) :: locks
173 : !$ INTEGER :: lock_num
174 :
175 35274 : CALL timeset(routineN, handle)
176 :
177 35274 : NULLIFY (my_kind_set, atomic_kind_set, force, my_oce, para_env, my_rho_atom, my_sab)
178 35274 : NULLIFY (mat_h, mat_p, dft_control)
179 :
180 : CALL get_qs_env(qs_env=qs_env, &
181 : qs_kind_set=my_kind_set, &
182 : atomic_kind_set=atomic_kind_set, &
183 : force=force, &
184 : oce=my_oce, &
185 : para_env=para_env, &
186 : rho_atom_set=my_rho_atom, &
187 : virial=virial, &
188 : sab_orb=my_sab, &
189 35274 : dft_control=dft_control)
190 :
191 35274 : nspins = dft_control%nspins
192 35274 : nimages = dft_control%nimages
193 :
194 35274 : factor1 = 1.0_dp
195 35274 : factor2 = 1.0_dp
196 :
197 : !deal with externals
198 35274 : my_tddft = .FALSE.
199 35274 : IF (PRESENT(tddft)) my_tddft = tddft
200 10758 : IF (my_tddft) THEN
201 5588 : IF (nspins == 1) factor1 = 2.0_dp
202 5588 : CPASSERT(nimages == 1)
203 : END IF
204 35274 : IF (PRESENT(kscale)) THEN
205 128 : factor1 = factor1*kscale
206 128 : factor2 = factor2*kscale
207 : END IF
208 35274 : IF (PRESENT(kintegral)) factor1 = kintegral
209 35274 : IF (PRESENT(kforce)) factor2 = kforce
210 105822 : force_fac = 1.0_dp
211 35274 : IF (PRESENT(fscale)) force_fac(:) = fscale(:)
212 :
213 35274 : IF (PRESENT(rho_atom_external)) my_rho_atom => rho_atom_external
214 35274 : IF (PRESENT(kind_set_external)) my_kind_set => kind_set_external
215 35274 : IF (PRESENT(oce_external)) my_oce => oce_external
216 35274 : IF (PRESENT(sab_external)) my_sab => sab_external
217 :
218 : ! kpoint images
219 35274 : NULLIFY (cell_to_index)
220 35274 : IF (nimages > 1) THEN
221 404 : CALL get_qs_env(qs_env=qs_env, kpoints=kpoints)
222 404 : CALL get_kpoint_info(kpoint=kpoints, cell_to_index=cell_to_index)
223 : END IF
224 :
225 35274 : eps_cpc = dft_control%qs_control%gapw_control%eps_cpc
226 :
227 35274 : CALL get_atomic_kind_set(atomic_kind_set, natom=natom)
228 35274 : CALL get_qs_kind_set(my_kind_set, maxsgf=max_nsgf, maxgtops=max_gau, basis_type="GAPW_1C")
229 :
230 35274 : IF (forces) THEN
231 1904 : CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind)
232 1904 : ldCPC = max_gau
233 1904 : use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
234 : ELSE
235 : use_virial = .FALSE.
236 : END IF
237 :
238 35274 : pv_virial_thread(:, :) = 0.0_dp ! always initialize to avoid floating point exceptions in OMP REDUCTION
239 :
240 35274 : nkind = SIZE(my_kind_set, 1)
241 : ! Collect the local potential contributions from all the processors
242 : ASSOCIATE (mepos => para_env%mepos, num_pe => para_env%num_pe)
243 105812 : DO ikind = 1, nkind
244 70538 : NULLIFY (atom_list)
245 70538 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=nat)
246 70538 : CALL get_qs_kind(my_kind_set(ikind), paw_atom=paw_atom)
247 105812 : IF (paw_atom) THEN
248 : ! gather the atomic block integrals in a more compressed format
249 64514 : bo = get_limit(nat, num_pe, mepos)
250 114252 : DO iat = bo(1), bo(2)
251 49738 : iatom = atom_list(iat)
252 170164 : DO ispin = 1, nspins
253 : CALL prj_gather(my_rho_atom(iatom)%ga_Vlocal_gb_h(ispin)%r_coef, &
254 55912 : my_rho_atom(iatom)%int_scr_h(ispin)%r_coef, my_kind_set(ikind))
255 : CALL prj_gather(my_rho_atom(iatom)%ga_Vlocal_gb_s(ispin)%r_coef, &
256 105650 : my_rho_atom(iatom)%int_scr_s(ispin)%r_coef, my_kind_set(ikind))
257 : END DO
258 : END DO
259 : ! broadcast the CPC arrays to all processors (replicated data)
260 193542 : DO ip = 0, num_pe - 1
261 129028 : bo = get_limit(nat, num_pe, ip)
262 293018 : DO iat = bo(1), bo(2)
263 99476 : iatom = atom_list(iat)
264 340328 : DO ispin = 1, nspins
265 106040592 : CALL para_env%bcast(my_rho_atom(iatom)%int_scr_h(ispin)%r_coef, ip)
266 106140068 : CALL para_env%bcast(my_rho_atom(iatom)%int_scr_s(ispin)%r_coef, ip)
267 : END DO
268 : END DO
269 : END DO
270 : END IF
271 : END DO
272 : END ASSOCIATE
273 :
274 176360 : ALLOCATE (basis_set_list(nkind))
275 176370 : ALLOCATE (paw_kind(nkind), has_intac(nkind*nkind))
276 105812 : paw_kind(:) = .FALSE.
277 188444 : has_intac(:) = .FALSE.
278 105812 : DO ikind = 1, nkind
279 70538 : CALL get_qs_kind(my_kind_set(ikind), basis_set=basis_set_a, paw_atom=paw_kind(ikind))
280 105812 : IF (ASSOCIATED(basis_set_a)) THEN
281 70538 : basis_set_list(ikind)%gto_basis_set => basis_set_a
282 : ELSE
283 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
284 : END IF
285 : END DO
286 188444 : DO ikind = 1, nkind*nkind
287 188444 : has_intac(ikind) = ASSOCIATED(my_oce%intac(ikind)%alist)
288 : END DO
289 :
290 : ! build the hash table in serial...
291 : ! ... creation ...
292 35274 : CALL neighbor_list_iterator_create(nl_iterator, my_sab)
293 35274 : nl_table_num_slots = natom*natom/2 ! this is probably not optimal, but it seems a good start
294 35274 : CALL nl_hash_table_create(nl_hash_table, nmax=nl_table_num_slots)
295 : ! ... and population
296 1064333 : DO WHILE (neighbor_list_iterate(nl_iterator) == 0) ! build hash table in serial, so don't pass mepos
297 9261531 : ALLOCATE (task) ! They must be deallocated before the nl_hash_table is released
298 1029059 : CALL get_iterator_task(nl_iterator, task) ! build hash table in serial, so don't pass mepos
299 : ! tasks with the same key access the same blocks of H & P
300 1029059 : IF (task%iatom <= task%jatom) THEN
301 611871 : nl_table_key = natom*task%iatom + task%jatom
302 : ELSE
303 417188 : nl_table_key = natom*task%jatom + task%iatom
304 : END IF
305 1029059 : CALL nl_hash_table_add(nl_hash_table, nl_table_key, task)
306 : END DO
307 35274 : CALL neighbor_list_iterator_release(nl_iterator)
308 :
309 : ! Get the total number of (possibly empty) entries in the table
310 35274 : CALL nl_hash_table_status(nl_hash_table, nmax=nl_table_num_slots)
311 :
312 : !$OMP PARALLEL DEFAULT(NONE) &
313 : !$OMP SHARED(nl_table_num_slots, nl_hash_table &
314 : !$OMP , max_gau, max_nsgf, nspins, forces &
315 : !$OMP , basis_set_list, nimages, cell_to_index &
316 : !$OMP , ksmat, pmat, natom, nkind, my_kind_set, my_oce &
317 : !$OMP , my_rho_atom, factor1, factor2, use_virial &
318 : !$OMP , atom_of_kind, ldCPC, force, locks, force_fac &
319 : !$OMP , paw_kind, has_intac &
320 : !$OMP ) &
321 : !$OMP PRIVATE( slot_num, is_entry_null, TASK, is_task_valid &
322 : !$OMP , C_int_h, C_int_s, coc, a_matrix, p_matrix &
323 : !$OMP , mat_h, mat_p, dCPC_h, dCPC_s, PC_h, PC_s &
324 : !$OMP , int_local_h, int_local_s &
325 : !$OMP , ikind, jkind, iatom, jatom, cell_b &
326 : !$OMP , basis_set_a, basis_set_b, img &
327 : !$OMP , found, next_task &
328 : !$OMP , kkind, lock_num &
329 : !$OMP , iac, alist_ac, kac, n_cont_a, list_a &
330 : !$OMP , ibc, alist_bc, kbc, n_cont_b, list_b &
331 : !$OMP , katom, rho_at, nsoctot &
332 : !$OMP , C_coeff_hh_a, C_coeff_ss_a, dista, rac &
333 : !$OMP , C_coeff_hh_b, C_coeff_ss_b, distb, rbc &
334 : !$OMP , ia_kind, ja_kind, ka_kind, force_tmp &
335 : !$OMP ) &
336 : !$OMP REDUCTION(+:pv_virial_thread &
337 35274 : !$OMP )
338 :
339 : ALLOCATE (C_int_h(max_gau*max_nsgf), C_int_s(max_gau*max_nsgf), coc(max_gau*max_gau), &
340 : a_matrix(max_gau, max_gau), p_matrix(max_nsgf, max_nsgf))
341 :
342 : ALLOCATE (mat_h(nspins), mat_p(nspins))
343 : DO ispin = 1, nspins
344 : NULLIFY (mat_h(ispin)%array, mat_p(ispin)%array)
345 : END DO
346 :
347 : IF (forces) THEN
348 : ALLOCATE (dCPC_h(max_gau, max_gau), dCPC_s(max_gau, max_gau), &
349 : PC_h(max_nsgf, max_gau, nspins), PC_s(max_nsgf, max_gau, nspins))
350 : !$OMP SINGLE
351 : !$ ALLOCATE (locks(natom*nkind))
352 : !$OMP END SINGLE
353 :
354 : !$OMP DO
355 : !$ do lock_num = 1, natom*nkind
356 : !$ call omp_init_lock(locks(lock_num))
357 : !$ end do
358 : !$OMP END DO
359 : END IF
360 :
361 : ! Dynamic schedule to take account of the fact that some slots may be empty
362 : ! or contain 1 or more tasks
363 : !$OMP DO SCHEDULE(DYNAMIC,5)
364 : DO slot_num = 1, nl_table_num_slots
365 : CALL nl_hash_table_is_null(nl_hash_table, slot_num, is_entry_null)
366 :
367 : IF (.NOT. is_entry_null) THEN
368 : CALL nl_hash_table_get_from_index(nl_hash_table, slot_num, task)
369 :
370 : is_task_valid = ASSOCIATED(task)
371 : DO WHILE (is_task_valid)
372 :
373 : ikind = task%ikind
374 : jkind = task%jkind
375 : iatom = task%iatom
376 : jatom = task%jatom
377 : cell_b(:) = task%cell(:)
378 :
379 : basis_set_a => basis_set_list(ikind)%gto_basis_set
380 : IF (.NOT. ASSOCIATED(basis_set_a)) CYCLE
381 : basis_set_b => basis_set_list(jkind)%gto_basis_set
382 : IF (.NOT. ASSOCIATED(basis_set_b)) CYCLE
383 :
384 : IF (nimages > 1) THEN
385 : img = cell_to_index(cell_b(1), cell_b(2), cell_b(3))
386 : CPASSERT(img > 0)
387 : ELSE
388 : img = 1
389 : END IF
390 :
391 : DO ispin = 1, nspins
392 : NULLIFY (mat_h(ispin)%array, mat_p(ispin)%array)
393 :
394 : found = .FALSE.
395 : IF (iatom <= jatom) THEN
396 : CALL dbcsr_get_block_p(matrix=ksmat(nspins*(img - 1) + ispin)%matrix, &
397 : row=iatom, col=jatom, &
398 : BLOCK=mat_h(ispin)%array, found=found)
399 : ELSE
400 : CALL dbcsr_get_block_p(matrix=ksmat(nspins*(img - 1) + ispin)%matrix, &
401 : row=jatom, col=iatom, &
402 : BLOCK=mat_h(ispin)%array, found=found)
403 : END IF
404 : CPASSERT(found)
405 :
406 : IF (forces) THEN
407 : found = .FALSE.
408 : IF (iatom <= jatom) THEN
409 : CALL dbcsr_get_block_p(matrix=pmat(nspins*(img - 1) + ispin)%matrix, &
410 : row=iatom, col=jatom, &
411 : BLOCK=mat_p(ispin)%array, found=found)
412 : ELSE
413 : CALL dbcsr_get_block_p(matrix=pmat(nspins*(img - 1) + ispin)%matrix, &
414 : row=jatom, col=iatom, &
415 : BLOCK=mat_p(ispin)%array, found=found)
416 : END IF
417 : CPASSERT(found)
418 : END IF
419 : END DO
420 :
421 : DO kkind = 1, nkind
422 : IF (.NOT. paw_kind(kkind)) CYCLE
423 :
424 : iac = ikind + nkind*(kkind - 1)
425 : ibc = jkind + nkind*(kkind - 1)
426 : IF (.NOT. has_intac(iac)) CYCLE
427 : IF (.NOT. has_intac(ibc)) CYCLE
428 :
429 : CALL get_alist(my_oce%intac(iac), alist_ac, iatom)
430 : CALL get_alist(my_oce%intac(ibc), alist_bc, jatom)
431 : IF (.NOT. ASSOCIATED(alist_ac)) CYCLE
432 : IF (.NOT. ASSOCIATED(alist_bc)) CYCLE
433 :
434 : DO kac = 1, alist_ac%nclist
435 : DO kbc = 1, alist_bc%nclist
436 : IF (alist_ac%clist(kac)%catom /= alist_bc%clist(kbc)%catom) CYCLE
437 :
438 : IF (ALL(cell_b + alist_bc%clist(kbc)%cell - alist_ac%clist(kac)%cell == 0)) THEN
439 : n_cont_a = alist_ac%clist(kac)%nsgf_cnt
440 : n_cont_b = alist_bc%clist(kbc)%nsgf_cnt
441 : IF (n_cont_a == 0 .OR. n_cont_b == 0) CYCLE
442 :
443 : list_a => alist_ac%clist(kac)%sgf_list
444 : list_b => alist_bc%clist(kbc)%sgf_list
445 :
446 : katom = alist_ac%clist(kac)%catom
447 :
448 : IF (iatom == katom .AND. ALL(alist_ac%clist(kac)%cell == 0)) THEN
449 : C_coeff_hh_a => alist_ac%clist(kac)%achint
450 : C_coeff_ss_a => alist_ac%clist(kac)%acint
451 : dista = .FALSE.
452 : ELSE
453 : C_coeff_hh_a => alist_ac%clist(kac)%acint
454 : C_coeff_ss_a => alist_ac%clist(kac)%acint
455 : dista = .TRUE.
456 : END IF
457 :
458 : IF (jatom == katom .AND. ALL(alist_bc%clist(kbc)%cell == 0)) THEN
459 : C_coeff_hh_b => alist_bc%clist(kbc)%achint
460 : C_coeff_ss_b => alist_bc%clist(kbc)%acint
461 : distb = .FALSE.
462 : ELSE
463 : C_coeff_hh_b => alist_bc%clist(kbc)%acint
464 : C_coeff_ss_b => alist_bc%clist(kbc)%acint
465 : distb = .TRUE.
466 : END IF
467 :
468 : rho_at => my_rho_atom(katom)
469 : nsoctot = SIZE(C_coeff_ss_a, 2)
470 :
471 : CALL get_rho_atom(rho_atom=rho_at, int_scr_h=int_local_h, int_scr_s=int_local_s)
472 : CALL add_vhxca_to_ks(mat_h, C_coeff_hh_a, C_coeff_hh_b, C_coeff_ss_a, C_coeff_ss_b, &
473 : int_local_h, int_local_s, nspins, iatom, jatom, nsoctot, factor1, &
474 : list_a, n_cont_a, list_b, n_cont_b, C_int_h, C_int_s, a_matrix, dista, distb, coc)
475 :
476 : IF (forces) THEN
477 : IF (use_virial) THEN
478 : rac = alist_ac%clist(kac)%rac
479 : rbc = alist_bc%clist(kbc)%rac
480 : END IF
481 : ia_kind = atom_of_kind(iatom)
482 : ja_kind = atom_of_kind(jatom)
483 : ka_kind = atom_of_kind(katom)
484 : rho_at => my_rho_atom(katom)
485 : force_tmp(1:3, 1:3) = 0.0_dp
486 :
487 : CALL get_rho_atom(rho_atom=rho_at, int_scr_h=int_local_h, int_scr_s=int_local_s)
488 : IF (iatom <= jatom) THEN
489 : CALL add_vhxca_forces(mat_p, C_coeff_hh_a, C_coeff_hh_b, C_coeff_ss_a, C_coeff_ss_b, &
490 : int_local_h, int_local_s, force_tmp, nspins, iatom, jatom, nsoctot, &
491 : list_a, n_cont_a, list_b, n_cont_b, dCPC_h, dCPC_s, ldCPC, &
492 : PC_h, PC_s, p_matrix, force_fac)
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 : ELSE
508 : CALL add_vhxca_forces(mat_p, C_coeff_hh_b, C_coeff_hh_a, C_coeff_ss_b, C_coeff_ss_a, &
509 : int_local_h, int_local_s, force_tmp, nspins, jatom, iatom, nsoctot, &
510 : list_b, n_cont_b, list_a, n_cont_a, dCPC_h, dCPC_s, ldCPC, &
511 : PC_h, PC_s, p_matrix, force_fac)
512 : force_tmp = factor2*force_tmp
513 : !$ CALL omp_set_lock(locks((ka_kind - 1)*nkind + kkind))
514 : force(kkind)%vhxc_atom(1:3, ka_kind) = force(kkind)%vhxc_atom(1:3, ka_kind) + force_tmp(1:3, 3)
515 : !$ CALL omp_unset_lock(locks((ka_kind - 1)*nkind + kkind))
516 : !$ CALL omp_set_lock(locks((ia_kind - 1)*nkind + ikind))
517 : force(ikind)%vhxc_atom(1:3, ia_kind) = force(ikind)%vhxc_atom(1:3, ia_kind) + force_tmp(1:3, 2)
518 : !$ CALL omp_unset_lock(locks((ia_kind - 1)*nkind + ikind))
519 : !$ CALL omp_set_lock(locks((ja_kind - 1)*nkind + jkind))
520 : force(jkind)%vhxc_atom(1:3, ja_kind) = force(jkind)%vhxc_atom(1:3, ja_kind) + force_tmp(1:3, 1)
521 : !$ CALL omp_unset_lock(locks((ja_kind - 1)*nkind + jkind))
522 : IF (use_virial) THEN
523 : CALL virial_pair_force(pv_virial_thread, 1._dp, force_tmp(1:3, 2), rac)
524 : CALL virial_pair_force(pv_virial_thread, 1._dp, force_tmp(1:3, 1), rbc)
525 : END IF
526 : END IF
527 :
528 : END IF
529 : EXIT ! search loop over jatom-katom list
530 : END IF
531 : END DO ! kbc
532 : END DO ! kac
533 :
534 : END DO ! kkind
535 :
536 : next_task => task%next
537 : ! We are done with this task, we can deallocate it
538 : DEALLOCATE (task)
539 : is_task_valid = ASSOCIATED(next_task)
540 : IF (is_task_valid) task => next_task
541 :
542 : END DO
543 :
544 : ELSE
545 : ! NO KEY/VALUE
546 : END IF
547 : END DO
548 : !$OMP END DO
549 :
550 : DO ispin = 1, nspins
551 : NULLIFY (mat_h(ispin)%array, mat_p(ispin)%array)
552 : END DO
553 : DEALLOCATE (mat_h, mat_p, C_int_h, C_int_s, a_matrix, p_matrix, coc)
554 :
555 : IF (forces) THEN
556 : DEALLOCATE (dCPC_h, dCPC_s, PC_h, PC_s)
557 :
558 : ! Implicit barrier at end of OMP DO, so locks can be freed
559 : !$OMP DO
560 : !$ DO lock_num = 1, natom*nkind
561 : !$ call omp_destroy_lock(locks(lock_num))
562 : !$ END DO
563 : !$OMP END DO
564 :
565 : !$OMP SINGLE
566 : !$ DEALLOCATE (locks)
567 : !$OMP END SINGLE NOWAIT
568 : END IF
569 :
570 : !$OMP END PARALLEL
571 :
572 35274 : IF (use_virial) THEN
573 1014 : virial%pv_gapw(1:3, 1:3) = virial%pv_gapw(1:3, 1:3) + pv_virial_thread(1:3, 1:3)
574 1014 : virial%pv_virial(1:3, 1:3) = virial%pv_virial(1:3, 1:3) + pv_virial_thread(1:3, 1:3)
575 : END IF
576 :
577 35274 : CALL nl_hash_table_release(nl_hash_table)
578 :
579 35274 : DEALLOCATE (basis_set_list, paw_kind, has_intac)
580 :
581 35274 : CALL timestop(handle)
582 :
583 141096 : END SUBROUTINE update_ks_atom
584 :
585 : ! **************************************************************************************************
586 : !> \brief ...
587 : !> \param mat_h ...
588 : !> \param C_hh_a ...
589 : !> \param C_hh_b ...
590 : !> \param C_ss_a ...
591 : !> \param C_ss_b ...
592 : !> \param int_local_h ...
593 : !> \param int_local_s ...
594 : !> \param nspins ...
595 : !> \param ia ...
596 : !> \param ja ...
597 : !> \param nsp ...
598 : !> \param factor ...
599 : !> \param lista ...
600 : !> \param nconta ...
601 : !> \param listb ...
602 : !> \param ncontb ...
603 : !> \param C_int_h ...
604 : !> \param C_int_s ...
605 : !> \param a_matrix ...
606 : !> \param dista ...
607 : !> \param distb ...
608 : !> \param coc ...
609 : ! **************************************************************************************************
610 8109497 : SUBROUTINE add_vhxca_to_ks(mat_h, C_hh_a, C_hh_b, C_ss_a, C_ss_b, &
611 : int_local_h, int_local_s, &
612 8109497 : nspins, ia, ja, nsp, factor, lista, nconta, listb, ncontb, &
613 16218994 : C_int_h, C_int_s, a_matrix, dista, distb, coc)
614 : TYPE(cp_2d_r_p_type), DIMENSION(:), POINTER :: mat_h
615 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: C_hh_a, C_hh_b, C_ss_a, C_ss_b
616 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: int_local_h, int_local_s
617 : INTEGER, INTENT(IN) :: nspins, ia, ja, nsp
618 : REAL(KIND=dp), INTENT(IN) :: factor
619 : INTEGER, DIMENSION(:), INTENT(IN) :: lista
620 : INTEGER, INTENT(IN) :: nconta
621 : INTEGER, DIMENSION(:), INTENT(IN) :: listb
622 : INTEGER, INTENT(IN) :: ncontb
623 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: C_int_h, C_int_s
624 : REAL(dp), DIMENSION(:, :) :: a_matrix
625 : LOGICAL, INTENT(IN) :: dista, distb
626 : REAL(dp), DIMENSION(:) :: coc
627 :
628 : INTEGER :: i, ispin, j, k
629 8109497 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: h_block, int_hard, int_soft
630 :
631 8109497 : NULLIFY (int_hard, int_soft)
632 :
633 16556134 : DO ispin = 1, nspins
634 8446637 : h_block => mat_h(ispin)%array
635 : !
636 8446637 : int_hard => int_local_h(ispin)%r_coef
637 8446637 : int_soft => int_local_s(ispin)%r_coef
638 : !
639 16556134 : IF (ia <= ja) THEN
640 4895941 : IF (dista .AND. distb) THEN
641 4662936 : k = 0
642 112612907 : DO i = 1, nsp
643 3217413632 : DO j = 1, nsp
644 3104800725 : k = k + 1
645 3212750696 : coc(k) = int_hard(j, i) - int_soft(j, i)
646 : END DO
647 : END DO
648 : CALL DGEMM('N', 'T', nsp, ncontb, nsp, 1.0_dp, coc, nsp, C_hh_b(:, :, 1), SIZE(C_hh_b, 1), &
649 4662936 : 0.0_dp, C_int_h, nsp)
650 : CALL DGEMM('N', 'N', nconta, ncontb, nsp, factor, C_hh_a(:, :, 1), SIZE(C_hh_a, 1), &
651 4662936 : C_int_h, nsp, 0.0_dp, a_matrix, SIZE(a_matrix, 1))
652 233005 : ELSEIF (dista) THEN
653 : CALL DGEMM('N', 'T', nsp, ncontb, nsp, 1.0_dp, int_hard, SIZE(int_hard, 1), &
654 82301 : C_hh_b(:, :, 1), SIZE(C_hh_b, 1), 0.0_dp, coc, nsp)
655 : CALL DGEMM('N', 'T', nsp, ncontb, nsp, -1.0_dp, int_soft, SIZE(int_soft, 1), &
656 82301 : C_ss_b(:, :, 1), SIZE(C_ss_b, 1), 1.0_dp, coc, nsp)
657 : CALL DGEMM('N', 'N', nconta, ncontb, nsp, factor, C_hh_a(:, :, 1), SIZE(C_hh_a, 1), &
658 82301 : coc, nsp, 0.0_dp, a_matrix, SIZE(a_matrix, 1))
659 150704 : ELSEIF (distb) THEN
660 : CALL DGEMM('N', 'N', nconta, nsp, nsp, factor, C_hh_a(:, :, 1), SIZE(C_hh_a, 1), &
661 94792 : int_hard, SIZE(int_hard, 1), 0.0_dp, coc, nconta)
662 : CALL DGEMM('N', 'N', nconta, nsp, nsp, -factor, C_ss_a(:, :, 1), SIZE(C_ss_a, 1), &
663 94792 : int_soft, SIZE(int_soft, 1), 1.0_dp, coc, nconta)
664 : CALL DGEMM('N', 'T', nconta, ncontb, nsp, 1.0_dp, coc, nconta, &
665 94792 : C_hh_b(:, :, 1), SIZE(C_hh_b, 1), 0.0_dp, a_matrix, SIZE(a_matrix, 1))
666 : ELSE
667 : CALL DGEMM('N', 'T', nsp, ncontb, nsp, 1.0_dp, int_hard, SIZE(int_hard, 1), &
668 : C_hh_b(:, :, 1), SIZE(C_hh_b, 1), &
669 55912 : 0.0_dp, C_int_h, nsp)
670 : CALL DGEMM('N', 'T', nsp, ncontb, nsp, 1.0_dp, int_soft, SIZE(int_soft, 1), &
671 : C_ss_b(:, :, 1), SIZE(C_ss_b, 1), &
672 55912 : 0.0_dp, C_int_s, nsp)
673 : CALL DGEMM('N', 'N', nconta, ncontb, nsp, factor, C_hh_a(:, :, 1), SIZE(C_hh_a, 1), &
674 : C_int_h, nsp, &
675 55912 : 0.0_dp, a_matrix, SIZE(a_matrix, 1))
676 : CALL DGEMM('N', 'N', nconta, ncontb, nsp, -factor, C_ss_a(:, :, 1), SIZE(C_ss_a, 1), &
677 : C_int_s, nsp, &
678 55912 : 1.0_dp, a_matrix, SIZE(a_matrix, 1))
679 : END IF
680 : !
681 : CALL alist_post_align_blk(a_matrix, SIZE(a_matrix, 1), h_block, SIZE(h_block, 1), &
682 4895941 : lista, nconta, listb, ncontb)
683 : ELSE
684 3550696 : IF (dista .AND. distb) THEN
685 3365408 : k = 0
686 80206288 : DO i = 1, nsp
687 2159582418 : DO j = 1, nsp
688 2079376130 : k = k + 1
689 2156217010 : coc(k) = int_hard(j, i) - int_soft(j, i)
690 : END DO
691 : END DO
692 3365408 : CALL DGEMM('N', 'T', nsp, nconta, nsp, 1.0_dp, coc, nsp, C_hh_a(:, :, 1), SIZE(C_hh_a, 1), 0.0_dp, C_int_h, nsp)
693 : CALL DGEMM('N', 'N', ncontb, nconta, nsp, factor, C_hh_b(:, :, 1), SIZE(C_hh_b, 1), &
694 3365408 : C_int_h, nsp, 0.0_dp, a_matrix, SIZE(a_matrix, 1))
695 185288 : ELSEIF (dista) THEN
696 : CALL DGEMM('N', 'N', ncontb, nsp, nsp, factor, C_hh_b(:, :, 1), SIZE(C_hh_b, 1), &
697 99883 : int_hard, SIZE(int_hard, 1), 0.0_dp, coc, ncontb)
698 : CALL DGEMM('N', 'N', ncontb, nsp, nsp, -factor, C_ss_b(:, :, 1), SIZE(C_ss_b, 1), &
699 99883 : int_soft, SIZE(int_soft, 1), 1.0_dp, coc, ncontb)
700 : CALL DGEMM('N', 'T', ncontb, nconta, nsp, 1.0_dp, coc, ncontb, &
701 99883 : C_hh_a(:, :, 1), SIZE(C_hh_a, 1), 0.0_dp, a_matrix, SIZE(a_matrix, 1))
702 85405 : ELSEIF (distb) THEN
703 : CALL DGEMM('N', 'T', nsp, nconta, nsp, 1.0_dp, int_hard, SIZE(int_hard, 1), &
704 85405 : C_hh_a(:, :, 1), SIZE(C_hh_a, 1), 0.0_dp, coc, nsp)
705 : CALL DGEMM('N', 'T', nsp, nconta, nsp, -1.0_dp, int_soft, SIZE(int_soft, 1), &
706 85405 : C_ss_a(:, :, 1), SIZE(C_ss_a, 1), 1.0_dp, coc, nsp)
707 : CALL DGEMM('N', 'N', ncontb, nconta, nsp, factor, C_hh_b(:, :, 1), SIZE(C_hh_b, 1), &
708 85405 : coc, nsp, 0.0_dp, a_matrix, SIZE(a_matrix, 1))
709 : ELSE
710 : CALL DGEMM('N', 'T', nsp, nconta, nsp, 1.0_dp, int_hard, SIZE(int_hard, 1), &
711 : C_hh_a(:, :, 1), SIZE(C_hh_a, 1), &
712 0 : 0.0_dp, C_int_h, nsp)
713 : CALL DGEMM('N', 'T', nsp, nconta, nsp, 1.0_dp, int_soft, SIZE(int_soft, 1), &
714 : C_ss_a(:, :, 1), SIZE(C_ss_a, 1), &
715 0 : 0.0_dp, C_int_s, nsp)
716 : CALL DGEMM('N', 'N', ncontb, nconta, nsp, factor, C_hh_b(:, :, 1), SIZE(C_hh_b, 1), &
717 : C_int_h, nsp, &
718 0 : 0.0_dp, a_matrix, SIZE(a_matrix, 1))
719 : CALL DGEMM('N', 'N', ncontb, nconta, nsp, -factor, C_ss_b(:, :, 1), SIZE(C_ss_b, 1), &
720 : C_int_s, nsp, &
721 0 : 1.0_dp, a_matrix, SIZE(a_matrix, 1))
722 : END IF
723 : !
724 : CALL alist_post_align_blk(a_matrix, SIZE(a_matrix, 1), h_block, SIZE(h_block, 1), &
725 3550696 : listb, ncontb, lista, nconta)
726 : END IF
727 : END DO
728 :
729 8109497 : END SUBROUTINE add_vhxca_to_ks
730 :
731 : ! **************************************************************************************************
732 : !> \brief ...
733 : !> \param mat_p ...
734 : !> \param C_hh_a ...
735 : !> \param C_hh_b ...
736 : !> \param C_ss_a ...
737 : !> \param C_ss_b ...
738 : !> \param int_local_h ...
739 : !> \param int_local_s ...
740 : !> \param force ...
741 : !> \param nspins ...
742 : !> \param ia ...
743 : !> \param ja ...
744 : !> \param nsp ...
745 : !> \param lista ...
746 : !> \param nconta ...
747 : !> \param listb ...
748 : !> \param ncontb ...
749 : !> \param dCPC_h ...
750 : !> \param dCPC_s ...
751 : !> \param ldCPC ...
752 : !> \param PC_h ...
753 : !> \param PC_s ...
754 : !> \param p_matrix ...
755 : !> \param force_scaling ...
756 : ! **************************************************************************************************
757 6994494 : SUBROUTINE add_vhxca_forces(mat_p, C_hh_a, C_hh_b, C_ss_a, C_ss_b, &
758 : int_local_h, int_local_s, &
759 777166 : force, nspins, ia, ja, nsp, lista, nconta, listb, ncontb, &
760 1554332 : dCPC_h, dCPC_s, ldCPC, PC_h, PC_s, p_matrix, force_scaling)
761 : TYPE(cp_2d_r_p_type), DIMENSION(:), INTENT(IN), &
762 : POINTER :: mat_p
763 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(IN) :: C_hh_a, C_hh_b, C_ss_a, C_ss_b
764 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: int_local_h, int_local_s
765 : REAL(dp), DIMENSION(3, 3), INTENT(INOUT) :: force
766 : INTEGER, INTENT(IN) :: nspins, ia, ja, nsp
767 : INTEGER, DIMENSION(:), INTENT(IN) :: lista
768 : INTEGER, INTENT(IN) :: nconta
769 : INTEGER, DIMENSION(:), INTENT(IN) :: listb
770 : INTEGER, INTENT(IN) :: ncontb
771 : REAL(KIND=dp), DIMENSION(:, :) :: dCPC_h, dCPC_s
772 : INTEGER, INTENT(IN) :: ldCPC
773 : REAL(KIND=dp), DIMENSION(:, :, :) :: PC_h, PC_s
774 : REAL(KIND=dp), DIMENSION(:, :) :: p_matrix
775 : REAL(KIND=dp), DIMENSION(2), INTENT(IN) :: force_scaling
776 :
777 : INTEGER :: dir, ispin
778 777166 : REAL(dp), DIMENSION(:, :), POINTER :: int_hard, int_soft
779 : REAL(KIND=dp) :: ieqj, trace
780 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: p_block
781 :
782 : ! if(dista.and.distb) we could also here use ChPCh = CsPCs
783 : ! *** factor 2 because only half of the pairs with ia =/ ja are considered
784 :
785 777166 : ieqj = 2.0_dp
786 208460 : IF (ia == ja) ieqj = 1.0_dp
787 :
788 777166 : NULLIFY (int_hard, int_soft)
789 :
790 1557266 : DO ispin = 1, nspins
791 780100 : p_block => mat_p(ispin)%array
792 :
793 : CALL alist_pre_align_blk(p_block, SIZE(p_block, 1), p_matrix, SIZE(p_matrix, 1), &
794 780100 : lista, nconta, listb, ncontb)
795 :
796 780100 : int_hard => int_local_h(ispin)%r_coef
797 780100 : int_soft => int_local_s(ispin)%r_coef
798 :
799 : CALL DGEMM('N', 'N', nconta, nsp, ncontb, 1.0_dp, p_matrix, SIZE(p_matrix, 1), &
800 : C_hh_b(:, :, 1), SIZE(C_hh_b, 1), &
801 780100 : 0.0_dp, PC_h(:, :, ispin), SIZE(PC_h, 1))
802 : CALL DGEMM('N', 'N', nconta, nsp, ncontb, 1.0_dp, p_matrix(:, :), SIZE(p_matrix, 1), &
803 : C_ss_b(:, :, 1), SIZE(C_ss_b, 1), &
804 780100 : 0.0_dp, PC_s(:, :, ispin), SIZE(PC_s, 1))
805 :
806 3120400 : DO dir = 2, 4
807 : CALL DGEMM('T', 'N', nsp, nsp, nconta, 1.0_dp, PC_h(:, :, ispin), SIZE(PC_h, 1), &
808 : C_hh_a(:, :, dir), SIZE(C_hh_a, 1), &
809 2340300 : 0.0_dp, dCPC_h, SIZE(dCPC_h, 1))
810 2340300 : trace = trace_r_AxB(dCPC_h, ldCPC, int_hard, nsp, nsp, nsp)
811 2340300 : force(dir - 1, 3) = force(dir - 1, 3) + ieqj*trace*force_scaling(ispin)
812 2340300 : force(dir - 1, 1) = force(dir - 1, 1) - ieqj*trace*force_scaling(ispin)
813 :
814 : CALL DGEMM('T', 'N', nsp, nsp, nconta, 1.0_dp, PC_s(:, :, ispin), SIZE(PC_s, 1), &
815 : C_ss_a(:, :, dir), SIZE(C_ss_a, 1), &
816 2340300 : 0.0_dp, dCPC_s, SIZE(dCPC_s, 1))
817 2340300 : trace = trace_r_AxB(dCPC_s, ldCPC, int_soft, nsp, nsp, nsp)
818 2340300 : force(dir - 1, 3) = force(dir - 1, 3) - ieqj*trace*force_scaling(ispin)
819 3120400 : force(dir - 1, 1) = force(dir - 1, 1) + ieqj*trace*force_scaling(ispin)
820 : END DO
821 :
822 : ! j-k contributions
823 : CALL DGEMM('T', 'N', ncontb, nsp, nconta, 1.0_dp, p_matrix, SIZE(p_matrix, 1), &
824 : C_hh_a(:, :, 1), SIZE(C_hh_a, 1), &
825 780100 : 0.0_dp, PC_h(:, :, ispin), SIZE(PC_h, 1))
826 : CALL DGEMM('T', 'N', ncontb, nsp, nconta, 1.0_dp, p_matrix, SIZE(p_matrix, 1), &
827 : C_ss_a(:, :, 1), SIZE(C_ss_a, 1), &
828 780100 : 0.0_dp, PC_s(:, :, ispin), SIZE(PC_s, 1))
829 :
830 3897566 : DO dir = 2, 4
831 : CALL DGEMM('T', 'N', nsp, nsp, ncontb, 1.0_dp, PC_h(:, :, ispin), SIZE(PC_h, 1), &
832 : C_hh_b(:, :, dir), SIZE(C_hh_b, 1), &
833 2340300 : 0.0_dp, dCPC_h, SIZE(dCPC_h, 1))
834 2340300 : trace = trace_r_AxB(dCPC_h, ldCPC, int_hard, nsp, nsp, nsp)
835 2340300 : force(dir - 1, 3) = force(dir - 1, 3) + ieqj*trace*force_scaling(ispin)
836 2340300 : force(dir - 1, 2) = force(dir - 1, 2) - ieqj*trace*force_scaling(ispin)
837 :
838 : CALL DGEMM('T', 'N', nsp, nsp, ncontb, 1.0_dp, PC_s(:, :, ispin), SIZE(PC_s, 1), &
839 : C_ss_b(:, :, dir), SIZE(C_ss_b, 1), &
840 2340300 : 0.0_dp, dCPC_s, SIZE(dCPC_s, 1))
841 2340300 : trace = trace_r_AxB(dCPC_s, ldCPC, int_soft, nsp, nsp, nsp)
842 2340300 : force(dir - 1, 3) = force(dir - 1, 3) - ieqj*trace*force_scaling(ispin)
843 3120400 : force(dir - 1, 2) = force(dir - 1, 2) + ieqj*trace*force_scaling(ispin)
844 : END DO
845 :
846 : END DO !ispin
847 :
848 777166 : END SUBROUTINE add_vhxca_forces
849 :
850 : END MODULE qs_ks_atom
|
