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