(git:ccc2433)
kg_tnadd_mat.F
Go to the documentation of this file.
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 !> \brief Calculation of the local potential contribution of the nonadditive kinetic energy
9 !> <a|V(local)|b> = <a|Sum e^a*rc**2|b>
10 !> \par History
11 !> - adapted from core_ppl
12 ! **************************************************************************************************
13 MODULE kg_tnadd_mat
14  USE ai_overlap_ppl, ONLY: ppl_integral
15  USE atomic_kind_types, ONLY: atomic_kind_type,&
16  get_atomic_kind_set
17  USE basis_set_types, ONLY: gto_basis_set_p_type,&
18  gto_basis_set_type
19  USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
20  USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set,&
21  dbcsr_deallocate_matrix_set
22  USE dbcsr_api, ONLY: dbcsr_add,&
23  dbcsr_create,&
24  dbcsr_distribution_type,&
25  dbcsr_finalize,&
26  dbcsr_get_block_p,&
27  dbcsr_p_type,&
28  dbcsr_set,&
29  dbcsr_type_symmetric
30  USE external_potential_types, ONLY: get_potential,&
31  local_potential_type
32  USE kg_environment_types, ONLY: kg_environment_type
33  USE kinds, ONLY: dp
34  USE orbital_pointers, ONLY: init_orbital_pointers,&
35  ncoset
36  USE particle_methods, ONLY: get_particle_set
37  USE particle_types, ONLY: particle_type
38  USE qs_force_types, ONLY: qs_force_type
39  USE qs_kind_types, ONLY: get_qs_kind,&
40  get_qs_kind_set,&
41  qs_kind_type
42  USE qs_neighbor_list_types, ONLY: get_iterator_info,&
43  neighbor_list_iterate,&
44  neighbor_list_iterator_create,&
45  neighbor_list_iterator_p_type,&
46  neighbor_list_iterator_release,&
47  neighbor_list_set_p_type,&
48  nl_set_sub_iterator,&
49  nl_sub_iterate
50  USE virial_methods, ONLY: virial_pair_force
51  USE virial_types, ONLY: virial_type
52 
53 !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num, omp_get_num_threads
54 
55 #include "./base/base_uses.f90"
56 
57  IMPLICIT NONE
58 
59  PRIVATE
60 
61  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'kg_tnadd_mat'
62 
63  PUBLIC :: build_tnadd_mat
64 
65 CONTAINS
66 
67 !==========================================================================================================
68 
69 ! **************************************************************************************************
70 !> \brief ...
71 !> \param kg_env ...
72 !> \param matrix_p ...
73 !> \param force ...
74 !> \param virial ...
75 !> \param calculate_forces ...
76 !> \param use_virial ...
77 !> \param qs_kind_set ...
78 !> \param atomic_kind_set ...
79 !> \param particle_set ...
80 !> \param sab_orb ...
81 !> \param dbcsr_dist ...
82 ! **************************************************************************************************
83  SUBROUTINE build_tnadd_mat(kg_env, matrix_p, force, virial, calculate_forces, use_virial, &
84  qs_kind_set, atomic_kind_set, particle_set, sab_orb, dbcsr_dist)
85 
86  TYPE(kg_environment_type), POINTER :: kg_env
87  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: matrix_p
88  TYPE(qs_force_type), DIMENSION(:), POINTER :: force
89  TYPE(virial_type), POINTER :: virial
90  LOGICAL, INTENT(IN) :: calculate_forces
91  LOGICAL :: use_virial
92  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
93  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
94  TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
95  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
96  POINTER :: sab_orb
97  TYPE(dbcsr_distribution_type), POINTER :: dbcsr_dist
98 
99  CHARACTER(LEN=*), PARAMETER :: routinen = 'build_tnadd_mat'
100  INTEGER, PARAMETER :: nexp_max = 10
101 
102  INTEGER :: atom_a, atom_b, atom_c, handle, i, iatom, icol, ikind, img, imol, inode, irow, &
103  iset, jatom, jkind, jmol, jset, katom, kkind, kmol, last_iatom, last_jatom, ldai, ldsab, &
104  maxco, maxder, maxl, maxlgto, maxnset, maxpol, maxsgf, mepos, natom, ncoa, ncob, nder, &
105  ngau, nkind, npol, nseta, nsetb, nthread, sgfa, sgfb
106  INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
107  INTEGER, DIMENSION(:), POINTER :: atom_to_molecule, col_blk_sizes, la_max, &
108  la_min, lb_max, lb_min, npgfa, npgfb, &
109  nsgfa, nsgfb, row_blk_sizes
110  INTEGER, DIMENSION(:, :), POINTER :: first_sgfa, first_sgfb
111  INTEGER, DIMENSION(nexp_max) :: nct
112  LOGICAL :: found, new_atom_pair
113  REAL(kind=dp) :: dab, dac, dbc, f0, radius
114  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: work
115  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: ppl_fwork, ppl_work
116  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: hab, pab
117  REAL(kind=dp), DIMENSION(3) :: force_a, force_b, rab, rac, rbc
118  REAL(kind=dp), DIMENSION(:), POINTER :: alpha, set_radius_a, set_radius_b
119  REAL(kind=dp), DIMENSION(:, :), POINTER :: ccval, cval, h_block, p_block, rpgfa, &
120  rpgfb, sphi_a, sphi_b, zeta, zetb
121  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_kg
122  TYPE(gto_basis_set_p_type), DIMENSION(:), POINTER :: basis_set_list
123  TYPE(gto_basis_set_type), POINTER :: basis_set_a, basis_set_b
124  TYPE(local_potential_type), POINTER :: tnadd_potential
125  TYPE(neighbor_list_iterator_p_type), &
126  DIMENSION(:), POINTER :: ap_iterator, nl_iterator
127  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
128  POINTER :: sac_kin
129 
130  IF (calculate_forces) THEN
131  CALL timeset(routinen//"_forces", handle)
132  ELSE
133  CALL timeset(routinen, handle)
134  END IF
135 
136  NULLIFY (matrix_kg)
137  IF (ASSOCIATED(kg_env%tnadd_mat)) THEN
138  CALL dbcsr_deallocate_matrix_set(kg_env%tnadd_mat)
139  END IF
140  sac_kin => kg_env%sac_kin
141  atom_to_molecule => kg_env%atom_to_molecule
142 
143  nkind = SIZE(atomic_kind_set)
144  natom = SIZE(particle_set)
145 
146  CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
147 
148  IF (calculate_forces) THEN
149  nder = 1
150  IF (SIZE(matrix_p, 1) == 2) THEN
151  DO img = 1, SIZE(matrix_p, 2)
152  CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
153  alpha_scalar=1.0_dp, beta_scalar=1.0_dp)
154  CALL dbcsr_add(matrix_p(2, img)%matrix, matrix_p(1, img)%matrix, &
155  alpha_scalar=-2.0_dp, beta_scalar=1.0_dp)
156  END DO
157  END IF
158  ELSE
159  nder = 0
160  END IF
161 
162  maxder = ncoset(nder)
163 
164  CALL get_qs_kind_set(qs_kind_set, maxpol=maxpol, maxco=maxco, maxlgto=maxlgto, &
165  maxsgf=maxsgf, maxnset=maxnset)
166 
167  maxl = max(maxlgto, maxpol)
168  CALL init_orbital_pointers(maxl + nder + 1)
169 
170  ldsab = max(maxco, ncoset(maxpol), maxsgf, maxpol)
171  ldai = ncoset(maxl + nder + 1)
172 
173  ALLOCATE (basis_set_list(nkind))
174  DO ikind = 1, nkind
175  CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set_a)
176  IF (ASSOCIATED(basis_set_a)) THEN
177  basis_set_list(ikind)%gto_basis_set => basis_set_a
178  ELSE
179  NULLIFY (basis_set_list(ikind)%gto_basis_set)
180  END IF
181  END DO
182 
183  ! build the matrix
184  ALLOCATE (row_blk_sizes(natom), col_blk_sizes(natom))
185 
186  CALL get_particle_set(particle_set, qs_kind_set, nsgf=row_blk_sizes)
187  CALL get_particle_set(particle_set, qs_kind_set, nsgf=col_blk_sizes)
188 
189  CALL dbcsr_allocate_matrix_set(matrix_kg, 1)
190 
191  ALLOCATE (matrix_kg(1)%matrix)
192  CALL dbcsr_create(matrix=matrix_kg(1)%matrix, &
193  name="Nonadditive kinetic energy potential", &
194  dist=dbcsr_dist, matrix_type=dbcsr_type_symmetric, &
195  row_blk_size=row_blk_sizes, col_blk_size=col_blk_sizes, &
196  nze=0, reuse_arrays=.true.)
197  CALL cp_dbcsr_alloc_block_from_nbl(matrix_kg(1)%matrix, sab_orb)
198  CALL dbcsr_set(matrix_kg(1)%matrix, 0.0_dp)
199 
200  nthread = 1
201 !$ nthread = omp_get_max_threads()
202 
203  CALL neighbor_list_iterator_create(nl_iterator, sab_orb, nthread=nthread)
204  ! iterator for basis/potential list
205  CALL neighbor_list_iterator_create(ap_iterator, sac_kin, search=.true., nthread=nthread)
206 
207 !$OMP PARALLEL &
208 !$OMP DEFAULT (NONE) &
209 !$OMP SHARED (nl_iterator, ap_iterator, basis_set_list, calculate_forces, force, use_virial,&
210 !$OMP matrix_kg, matrix_p,virial, atomic_kind_set, qs_kind_set, particle_set,&
211 !$OMP sab_orb, sac_kin, nthread, ncoset, nkind,&
212 !$OMP atom_of_kind, ldsab, maxnset, maxder, &
213 !$OMP maxlgto, nder, maxco, atom_to_molecule) &
214 !$OMP PRIVATE (ikind, jkind, inode, iatom, jatom, rab, basis_set_a, basis_set_b, atom_a, &
215 !$OMP atom_b, first_sgfa, la_max, la_min, npgfa, nsgfa, sphi_a, &
216 !$OMP zeta, first_sgfb, lb_max, lb_min, npgfb, nsetb, rpgfb, set_radius_b, sphi_b, &
217 !$OMP zetb, last_iatom, last_jatom, new_atom_pair, dab, irow, icol, h_block, found, iset, ncoa, &
218 !$OMP sgfa, jset, ncob, sgfb, nsgfb, p_block, work, pab, hab, kkind, nseta, &
219 !$OMP radius, alpha, nct, imol, jmol, kmol,&
220 !$OMP npol, ngau, cval, ccval, rac, dac, rbc, dbc, &
221 !$OMP set_radius_a, rpgfa, force_a, force_b, ppl_fwork, mepos, &
222 !$OMP f0, katom, ppl_work, atom_c,&
223 !$OMP ldai,tnadd_potential)
224 
225  mepos = 0
226 !$ mepos = omp_get_thread_num()
227 
228  ALLOCATE (hab(ldsab, ldsab, maxnset, maxnset), work(ldsab, ldsab*maxder))
229  ldai = ncoset(2*maxlgto + 2*nder)
230  ALLOCATE (ppl_work(ldai, ldai, max(maxder, 2*maxlgto + 2*nder + 1)))
231  IF (calculate_forces) THEN
232  ALLOCATE (pab(maxco, maxco, maxnset, maxnset))
233  ldai = ncoset(maxlgto)
234  ALLOCATE (ppl_fwork(ldai, ldai, maxder))
235  END IF
236 
237  last_iatom = 0
238  last_jatom = 0
239  DO WHILE (neighbor_list_iterate(nl_iterator, mepos=mepos) == 0)
240 
241  CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, inode=inode, &
242  iatom=iatom, jatom=jatom, r=rab)
243 
244  basis_set_a => basis_set_list(ikind)%gto_basis_set
245  IF (.NOT. ASSOCIATED(basis_set_a)) cycle
246  basis_set_b => basis_set_list(jkind)%gto_basis_set
247  IF (.NOT. ASSOCIATED(basis_set_b)) cycle
248 
249  atom_a = atom_of_kind(iatom)
250  atom_b = atom_of_kind(jatom)
251  imol = atom_to_molecule(iatom)
252  jmol = atom_to_molecule(jatom)
253  cpassert(imol == jmol)
254 
255  ! basis ikind
256  first_sgfa => basis_set_a%first_sgf
257  la_max => basis_set_a%lmax
258  la_min => basis_set_a%lmin
259  npgfa => basis_set_a%npgf
260  nseta = basis_set_a%nset
261  nsgfa => basis_set_a%nsgf_set
262  rpgfa => basis_set_a%pgf_radius
263  set_radius_a => basis_set_a%set_radius
264  sphi_a => basis_set_a%sphi
265  zeta => basis_set_a%zet
266  ! basis jkind
267  first_sgfb => basis_set_b%first_sgf
268  lb_max => basis_set_b%lmax
269  lb_min => basis_set_b%lmin
270  npgfb => basis_set_b%npgf
271  nsetb = basis_set_b%nset
272  nsgfb => basis_set_b%nsgf_set
273  rpgfb => basis_set_b%pgf_radius
274  set_radius_b => basis_set_b%set_radius
275  sphi_b => basis_set_b%sphi
276  zetb => basis_set_b%zet
277 
278  dab = sqrt(sum(rab*rab))
279 
280  IF (iatom == last_iatom .AND. jatom == last_jatom) THEN
281  new_atom_pair = .false.
282  ELSE
283  new_atom_pair = .true.
284  last_jatom = jatom
285  last_iatom = iatom
286  END IF
287 
288  ! *** Use the symmetry of the first derivatives ***
289  IF (iatom == jatom) THEN
290  f0 = 1.0_dp
291  ELSE
292  f0 = 2.0_dp
293  END IF
294 
295  ! *** Create matrix blocks for a new matrix block column ***
296  IF (new_atom_pair) THEN
297  IF (iatom <= jatom) THEN
298  irow = iatom
299  icol = jatom
300  ELSE
301  irow = jatom
302  icol = iatom
303  END IF
304  NULLIFY (h_block)
305  CALL dbcsr_get_block_p(matrix_kg(1)%matrix, irow, icol, h_block, found)
306  IF (ASSOCIATED(h_block)) THEN
307  IF (calculate_forces) THEN
308  cpassert(SIZE(matrix_p, 2) == 1)
309  NULLIFY (p_block)
310  CALL dbcsr_get_block_p(matrix_p(1, 1)%matrix, irow, icol, p_block, found)
311  IF (ASSOCIATED(p_block)) THEN
312  DO iset = 1, nseta
313  ncoa = npgfa(iset)*ncoset(la_max(iset))
314  sgfa = first_sgfa(1, iset)
315  DO jset = 1, nsetb
316  ncob = npgfb(jset)*ncoset(lb_max(jset))
317  sgfb = first_sgfb(1, jset)
318 
319  ! *** Decontract density matrix block ***
320  IF (iatom <= jatom) THEN
321  CALL dgemm("N", "N", ncoa, nsgfb(jset), nsgfa(iset), &
322  1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
323  p_block(sgfa, sgfb), SIZE(p_block, 1), &
324  0.0_dp, work(1, 1), SIZE(work, 1))
325  ELSE
326  CALL dgemm("N", "T", ncoa, nsgfb(jset), nsgfa(iset), &
327  1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
328  p_block(sgfb, sgfa), SIZE(p_block, 1), &
329  0.0_dp, work(1, 1), SIZE(work, 1))
330  END IF
331 
332  CALL dgemm("N", "T", ncoa, ncob, nsgfb(jset), &
333  1.0_dp, work(1, 1), SIZE(work, 1), &
334  sphi_b(1, sgfb), SIZE(sphi_b, 1), &
335  0.0_dp, pab(1, 1, iset, jset), SIZE(pab, 1))
336  END DO
337  END DO
338  END IF
339  END IF
340  END IF
341  END IF
342 
343  hab = 0._dp
344 
345  ! loop over all kinds for nonadditive kinetic energy potential atoms
346  DO kkind = 1, nkind
347 
348  CALL get_qs_kind(qs_kind_set(kkind), tnadd_potential=tnadd_potential)
349  IF (.NOT. ASSOCIATED(tnadd_potential)) cycle
350  CALL get_potential(potential=tnadd_potential, &
351  alpha=alpha, cval=cval, ngau=ngau, npol=npol, radius=radius)
352  nct = npol
353  ALLOCATE (ccval(npol, ngau))
354  ccval(1:npol, 1:ngau) = transpose(cval(1:ngau, 1:npol))
355 
356  CALL nl_set_sub_iterator(ap_iterator, ikind, kkind, iatom, mepos)
357 
358  DO WHILE (nl_sub_iterate(ap_iterator, mepos) == 0)
359 
360  CALL get_iterator_info(ap_iterator, mepos, jatom=katom, r=rac)
361 
362  ! this potential only acts on other moleclules
363  kmol = atom_to_molecule(katom)
364  IF (kmol == imol) cycle
365 
366  dac = sqrt(sum(rac*rac))
367  rbc(:) = rac(:) - rab(:)
368  dbc = sqrt(sum(rbc*rbc))
369  IF ((maxval(set_radius_a(:)) + radius < dac) .OR. &
370  (maxval(set_radius_b(:)) + radius < dbc)) THEN
371  cycle
372  END IF
373 
374  DO iset = 1, nseta
375  IF (set_radius_a(iset) + radius < dac) cycle
376  ncoa = npgfa(iset)*ncoset(la_max(iset))
377  sgfa = first_sgfa(1, iset)
378  DO jset = 1, nsetb
379  IF (set_radius_b(jset) + radius < dbc) cycle
380  ncob = npgfb(jset)*ncoset(lb_max(jset))
381  sgfb = first_sgfb(1, jset)
382  IF (set_radius_a(iset) + set_radius_b(jset) < dab) cycle
383  ! *** Calculate the GTH pseudo potential forces ***
384  IF (calculate_forces) THEN
385 
386  CALL ppl_integral( &
387  la_max(iset), la_min(iset), npgfa(iset), &
388  rpgfa(:, iset), zeta(:, iset), &
389  lb_max(jset), lb_min(jset), npgfb(jset), &
390  rpgfb(:, jset), zetb(:, jset), &
391  ngau, alpha, nct, ccval, radius, &
392  rab, dab, rac, dac, rbc, dbc, &
393  hab(:, :, iset, jset), ppl_work, pab(:, :, iset, jset), &
394  force_a, force_b, ppl_fwork)
395  ! *** The derivatives w.r.t. atomic center c are ***
396  ! *** calculated using the translational invariance ***
397  ! *** of the first derivatives ***
398  atom_c = atom_of_kind(katom)
399 
400 !$OMP CRITICAL(force_critical)
401  force(ikind)%kinetic(1, atom_a) = force(ikind)%kinetic(1, atom_a) + f0*force_a(1)
402  force(ikind)%kinetic(2, atom_a) = force(ikind)%kinetic(2, atom_a) + f0*force_a(2)
403  force(ikind)%kinetic(3, atom_a) = force(ikind)%kinetic(3, atom_a) + f0*force_a(3)
404  force(kkind)%kinetic(1, atom_c) = force(kkind)%kinetic(1, atom_c) - f0*force_a(1)
405  force(kkind)%kinetic(2, atom_c) = force(kkind)%kinetic(2, atom_c) - f0*force_a(2)
406  force(kkind)%kinetic(3, atom_c) = force(kkind)%kinetic(3, atom_c) - f0*force_a(3)
407 
408  force(jkind)%kinetic(1, atom_b) = force(jkind)%kinetic(1, atom_b) + f0*force_b(1)
409  force(jkind)%kinetic(2, atom_b) = force(jkind)%kinetic(2, atom_b) + f0*force_b(2)
410  force(jkind)%kinetic(3, atom_b) = force(jkind)%kinetic(3, atom_b) + f0*force_b(3)
411  force(kkind)%kinetic(1, atom_c) = force(kkind)%kinetic(1, atom_c) - f0*force_b(1)
412  force(kkind)%kinetic(2, atom_c) = force(kkind)%kinetic(2, atom_c) - f0*force_b(2)
413  force(kkind)%kinetic(3, atom_c) = force(kkind)%kinetic(3, atom_c) - f0*force_b(3)
414 
415  IF (use_virial) THEN
416  CALL virial_pair_force(virial%pv_virial, f0, force_a, rac)
417  CALL virial_pair_force(virial%pv_virial, f0, force_b, rbc)
418  END IF
419 !$OMP END CRITICAL(force_critical)
420 
421  ELSE
422  CALL ppl_integral( &
423  la_max(iset), la_min(iset), npgfa(iset), &
424  rpgfa(:, iset), zeta(:, iset), &
425  lb_max(jset), lb_min(jset), npgfb(jset), &
426  rpgfb(:, jset), zetb(:, jset), &
427  ngau, alpha, nct, ccval, radius, &
428  rab, dab, rac, dac, rbc, dbc, hab(:, :, iset, jset), ppl_work)
429  END IF
430  END DO
431  END DO
432  END DO
433  DEALLOCATE (ccval)
434  END DO
435 
436  ! *** Contract integrals
437  DO iset = 1, nseta
438  ncoa = npgfa(iset)*ncoset(la_max(iset))
439  sgfa = first_sgfa(1, iset)
440  DO jset = 1, nsetb
441  ncob = npgfb(jset)*ncoset(lb_max(jset))
442  sgfb = first_sgfb(1, jset)
443 
444  CALL dgemm("N", "N", ncoa, nsgfb(jset), ncob, &
445  1.0_dp, hab(1, 1, iset, jset), SIZE(hab, 1), &
446  sphi_b(1, sgfb), SIZE(sphi_b, 1), &
447  0.0_dp, work(1, 1), SIZE(work, 1))
448 
449 !$OMP CRITICAL(h_block_critical)
450  IF (iatom <= jatom) THEN
451  CALL dgemm("T", "N", nsgfa(iset), nsgfb(jset), ncoa, &
452  1.0_dp, sphi_a(1, sgfa), SIZE(sphi_a, 1), &
453  work(1, 1), SIZE(work, 1), &
454  1.0_dp, h_block(sgfa, sgfb), SIZE(h_block, 1))
455  ELSE
456  CALL dgemm("T", "N", nsgfb(jset), nsgfa(iset), ncoa, &
457  1.0_dp, work(1, 1), SIZE(work, 1), &
458  sphi_a(1, sgfa), SIZE(sphi_a, 1), &
459  1.0_dp, h_block(sgfb, sgfa), SIZE(h_block, 1))
460  END IF
461 !$OMP END CRITICAL(h_block_critical)
462 
463  END DO
464  END DO
465  END DO
466 
467  DEALLOCATE (hab, work, ppl_work)
468 
469  IF (calculate_forces) THEN
470  DEALLOCATE (pab, ppl_fwork)
471  END IF
472 
473 !$OMP END PARALLEL
474 
475  CALL neighbor_list_iterator_release(ap_iterator)
476  CALL neighbor_list_iterator_release(nl_iterator)
477 
478  DO i = 1, SIZE(matrix_kg)
479  CALL dbcsr_finalize(matrix_kg(i)%matrix)
480  END DO
481  kg_env%tnadd_mat => matrix_kg
482 
483  DEALLOCATE (basis_set_list)
484 
485  IF (calculate_forces) THEN
486  ! *** If LSD, then recover alpha density and beta density ***
487  ! *** from the total density (1) and the spin density (2) ***
488  IF (SIZE(matrix_p, 1) == 2) THEN
489  DO img = 1, SIZE(matrix_p, 2)
490  CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
491  alpha_scalar=0.5_dp, beta_scalar=0.5_dp)
492  CALL dbcsr_add(matrix_p(2, img)%matrix, matrix_p(1, img)%matrix, &
493  alpha_scalar=-1.0_dp, beta_scalar=1.0_dp)
494  END DO
495  END IF
496  END IF
497 
498  CALL timestop(handle)
499 
500  END SUBROUTINE build_tnadd_mat
501 
502 !==========================================================================================================
503 
504 END MODULE kg_tnadd_mat
dgemm
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Definition: grid_cpu_task_list.c:195
cp_dbcsr_operations::dbcsr_allocate_matrix_set
Definition: cp_dbcsr_operations.F:81
cp_dbcsr_operations::dbcsr_deallocate_matrix_set
Definition: cp_dbcsr_operations.F:89
ai_overlap_ppl
Calculation of three-center overlap integrals over Cartesian Gaussian-type functions for the second t...
Definition: ai_overlap_ppl.F:27
ai_overlap_ppl::ppl_integral
subroutine, public ppl_integral(la_max_set, la_min_set, npgfa, rpgfa, zeta, lb_max_set, lb_min_set, npgfb, rpgfb, zetb, nexp_ppl, alpha_ppl, nct_ppl, cexp_ppl, rpgfc, rab, dab, rac, dac, rbc, dbc, vab, s, pab, force_a, force_b, fs, hab2, hab2_work, deltaR, iatom, jatom, katom)
Calculation of three-center overlap integrals <a|c|b> over Cartesian Gaussian functions for the local...
Definition: ai_overlap_ppl.F:103
atomic_kind_types
Define the atomic kind types and their sub types.
Definition: atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind_set
subroutine, public get_atomic_kind_set(atomic_kind_set, atom_of_kind, kind_of, natom_of_kind, maxatom, natom, nshell, fist_potential_present, shell_present, shell_adiabatic, shell_check_distance, damping_present)
Get attributes of an atomic kind set.
Definition: atomic_kind_types.F:223
basis_set_types
Definition: basis_set_types.F:15
cp_dbcsr_cp2k_link
Routines that link DBCSR and CP2K concepts together.
Definition: cp_dbcsr_cp2k_link.F:14
cp_dbcsr_cp2k_link::cp_dbcsr_alloc_block_from_nbl
subroutine, public cp_dbcsr_alloc_block_from_nbl(matrix, sab_orb, desymmetrize)
allocate the blocks of a dbcsr based on the neighbor list
Definition: cp_dbcsr_cp2k_link.F:445
cp_dbcsr_operations
DBCSR operations in CP2K.
Definition: cp_dbcsr_operations.F:18
external_potential_types
Definition of the atomic potential types.
Definition: external_potential_types.F:14
kg_environment_types
Types needed for a Kim-Gordon-like partitioning into molecular subunits.
Definition: kg_environment_types.F:15
kg_tnadd_mat
Calculation of the local potential contribution of the nonadditive kinetic energy <a|V(local)|b> = <a...
Definition: kg_tnadd_mat.F:13
kg_tnadd_mat::build_tnadd_mat
subroutine, public build_tnadd_mat(kg_env, matrix_p, force, virial, calculate_forces, use_virial, qs_kind_set, atomic_kind_set, particle_set, sab_orb, dbcsr_dist)
...
Definition: kg_tnadd_mat.F:85
kinds
Defines the basic variable types.
Definition: kinds.F:23
kinds::dp
integer, parameter, public dp
Definition: kinds.F:34
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition: orbital_pointers.F:15
orbital_pointers::init_orbital_pointers
subroutine, public init_orbital_pointers(maxl)
Initialize or update the orbital pointers.
Definition: orbital_pointers.F:253
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition: orbital_pointers.F:42
particle_methods
Define methods related to particle_type.
Definition: particle_methods.F:14
particle_methods::get_particle_set
subroutine, public get_particle_set(particle_set, qs_kind_set, first_sgf, last_sgf, nsgf, nmao, basis)
Get the components of a particle set.
Definition: particle_methods.F:98
particle_types
Define the data structure for the particle information.
Definition: particle_types.F:19
qs_force_types
Definition: qs_force_types.F:13
qs_kind_types
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_r3d_rs_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, U_of_dft_plus_u, J_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, J0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition: qs_kind_types.F:451
qs_kind_types::get_qs_kind_set
subroutine, public get_qs_kind_set(qs_kind_set, all_potential_present, tnadd_potential_present, gth_potential_present, sgp_potential_present, paw_atom_present, dft_plus_u_atom_present, maxcgf, maxsgf, maxco, maxco_proj, maxgtops, maxlgto, maxlprj, maxnset, maxsgf_set, ncgf, npgf, nset, nsgf, nshell, maxpol, maxlppl, maxlppnl, maxppnl, nelectron, maxder, max_ngrid_rad, max_sph_harm, maxg_iso_not0, lmax_rho0, basis_rcut, basis_type, total_zeff_corr)
Get attributes of an atomic kind set.
Definition: qs_kind_types.F:864
qs_neighbor_list_types
Define the neighbor list data types and the corresponding functionality.
Definition: qs_neighbor_list_types.F:17
qs_neighbor_list_types::nl_sub_iterate
integer function, public nl_sub_iterate(iterator_set, mepos)
...
Definition: qs_neighbor_list_types.F:440
qs_neighbor_list_types::neighbor_list_iterator_create
subroutine, public neighbor_list_iterator_create(iterator_set, nl, search, nthread)
Neighbor list iterator functions.
Definition: qs_neighbor_list_types.F:165
qs_neighbor_list_types::nl_set_sub_iterator
subroutine, public nl_set_sub_iterator(iterator_set, ikind, jkind, iatom, mepos)
...
Definition: qs_neighbor_list_types.F:288
qs_neighbor_list_types::neighbor_list_iterator_release
subroutine, public neighbor_list_iterator_release(iterator_set)
...
Definition: qs_neighbor_list_types.F:251
qs_neighbor_list_types::neighbor_list_iterate
integer function, public neighbor_list_iterate(iterator_set, mepos)
...
Definition: qs_neighbor_list_types.F:351
qs_neighbor_list_types::get_iterator_info
subroutine, public get_iterator_info(iterator_set, mepos, ikind, jkind, nkind, ilist, nlist, inode, nnode, iatom, jatom, r, cell)
...
Definition: qs_neighbor_list_types.F:549
virial_methods
Definition: virial_methods.F:12
virial_methods::virial_pair_force
pure subroutine, public virial_pair_force(pv_virial, f0, force, rab)
Computes the contribution to the stress tensor from two-body pair-wise forces.
Definition: virial_methods.F:139
virial_types
Definition: virial_types.F:13