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 : !> \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_api, ONLY: dbcsr_add,&
20 : dbcsr_create,&
21 : dbcsr_distribution_type,&
22 : dbcsr_finalize,&
23 : dbcsr_get_block_p,&
24 : dbcsr_p_type,&
25 : dbcsr_set,&
26 : dbcsr_type_symmetric
27 : USE cp_dbcsr_cp2k_link, ONLY: cp_dbcsr_alloc_block_from_nbl
28 : USE cp_dbcsr_operations, ONLY: dbcsr_allocate_matrix_set,&
29 : dbcsr_deallocate_matrix_set
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 42 : 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 42 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
107 42 : INTEGER, DIMENSION(:), POINTER :: atom_to_molecule, col_blk_sizes, la_max, &
108 42 : la_min, lb_max, lb_min, npgfa, npgfb, &
109 42 : nsgfa, nsgfb, row_blk_sizes
110 42 : 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 42 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: work
115 42 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: ppl_fwork, ppl_work
116 42 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: hab, pab
117 : REAL(KIND=dp), DIMENSION(3) :: force_a, force_b, rab, rac, rbc
118 42 : REAL(KIND=dp), DIMENSION(:), POINTER :: alpha, set_radius_a, set_radius_b
119 42 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: ccval, cval, h_block, p_block, rpgfa, &
120 42 : rpgfb, sphi_a, sphi_b, zeta, zetb
121 42 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_kg
122 42 : 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 42 : DIMENSION(:), POINTER :: ap_iterator, nl_iterator
127 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
128 42 : POINTER :: sac_kin
129 :
130 84 : IF (calculate_forces) THEN
131 20 : CALL timeset(routineN//"_forces", handle)
132 : ELSE
133 22 : CALL timeset(routineN, handle)
134 : END IF
135 :
136 42 : NULLIFY (matrix_kg)
137 42 : IF (ASSOCIATED(kg_env%tnadd_mat)) THEN
138 30 : CALL dbcsr_deallocate_matrix_set(kg_env%tnadd_mat)
139 : END IF
140 42 : sac_kin => kg_env%sac_kin
141 42 : atom_to_molecule => kg_env%atom_to_molecule
142 :
143 42 : nkind = SIZE(atomic_kind_set)
144 42 : natom = SIZE(particle_set)
145 :
146 42 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
147 :
148 42 : IF (calculate_forces) THEN
149 20 : nder = 1
150 20 : IF (SIZE(matrix_p, 1) == 2) THEN
151 0 : DO img = 1, SIZE(matrix_p, 2)
152 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
153 0 : 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 0 : 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 42 : maxder = ncoset(nder)
163 :
164 : CALL get_qs_kind_set(qs_kind_set, maxpol=maxpol, maxco=maxco, maxlgto=maxlgto, &
165 42 : maxsgf=maxsgf, maxnset=maxnset)
166 :
167 42 : maxl = MAX(maxlgto, maxpol)
168 42 : CALL init_orbital_pointers(maxl + nder + 1)
169 :
170 42 : ldsab = MAX(maxco, ncoset(maxpol), maxsgf, maxpol)
171 42 : ldai = ncoset(maxl + nder + 1)
172 :
173 192 : ALLOCATE (basis_set_list(nkind))
174 108 : DO ikind = 1, nkind
175 66 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_set_a)
176 108 : IF (ASSOCIATED(basis_set_a)) THEN
177 66 : basis_set_list(ikind)%gto_basis_set => basis_set_a
178 : ELSE
179 0 : NULLIFY (basis_set_list(ikind)%gto_basis_set)
180 : END IF
181 : END DO
182 :
183 : ! build the matrix
184 168 : ALLOCATE (row_blk_sizes(natom), col_blk_sizes(natom))
185 :
186 42 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=row_blk_sizes)
187 42 : CALL get_particle_set(particle_set, qs_kind_set, nsgf=col_blk_sizes)
188 :
189 42 : CALL dbcsr_allocate_matrix_set(matrix_kg, 1)
190 :
191 42 : 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 42 : nze=0, reuse_arrays=.TRUE.)
197 42 : CALL cp_dbcsr_alloc_block_from_nbl(matrix_kg(1)%matrix, sab_orb)
198 42 : CALL dbcsr_set(matrix_kg(1)%matrix, 0.0_dp)
199 :
200 : nthread = 1
201 42 : !$ nthread = omp_get_max_threads()
202 :
203 42 : CALL neighbor_list_iterator_create(nl_iterator, sab_orb, nthread=nthread)
204 : ! iterator for basis/potential list
205 42 : 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 42 : !$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 42 : CALL neighbor_list_iterator_release(ap_iterator)
476 42 : CALL neighbor_list_iterator_release(nl_iterator)
477 :
478 84 : DO i = 1, SIZE(matrix_kg)
479 84 : CALL dbcsr_finalize(matrix_kg(i)%matrix)
480 : END DO
481 42 : kg_env%tnadd_mat => matrix_kg
482 :
483 42 : DEALLOCATE (basis_set_list)
484 :
485 42 : 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 20 : IF (SIZE(matrix_p, 1) == 2) THEN
489 0 : DO img = 1, SIZE(matrix_p, 2)
490 : CALL dbcsr_add(matrix_p(1, img)%matrix, matrix_p(2, img)%matrix, &
491 0 : 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 0 : alpha_scalar=-1.0_dp, beta_scalar=1.0_dp)
494 : END DO
495 : END IF
496 : END IF
497 :
498 42 : CALL timestop(handle)
499 :
500 126 : END SUBROUTINE build_tnadd_mat
501 :
502 : !==========================================================================================================
503 :
504 : END MODULE kg_tnadd_mat
|
