(git:374b731)
Loading...
Searching...
No Matches
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! **************************************************************************************************
13MODULE 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
65CONTAINS
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
504END 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
external_potential_types::get_potential
Definition external_potential_types.F:276
qs_neighbor_list_types::nl_sub_iterate
Definition qs_neighbor_list_types.F:96
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::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
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
basis_set_types::gto_basis_set_p_type
Definition basis_set_types.F:93
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:71
external_potential_types::local_potential_type
Definition external_potential_types.F:98
kg_environment_types::kg_environment_type
Contains all the info needed for KG runs...
Definition kg_environment_types.F:102
particle_types::particle_type
Definition particle_types.F:35
qs_force_types::qs_force_type
Definition qs_force_types.F:28
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:171
qs_neighbor_list_types::neighbor_list_iterator_p_type
Definition qs_neighbor_list_types.F:117
qs_neighbor_list_types::neighbor_list_set_p_type
Definition qs_neighbor_list_types.F:67
virial_types::virial_type
Definition virial_types.F:26