54#include "./base/base_uses.f90"
62 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'qs_kinetic'
68 INTEGER,
DIMENSION(1:56),
SAVE :: ndod = (/0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 1, 1, 1, &
69 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, &
70 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, &
71 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1/)
99 basis_type, sab_nl, calculate_forces, matrix_p, matrixkp_p, &
100 eps_filter, nderivative)
106 POINTER :: matrixkp_t
107 CHARACTER(LEN=*),
INTENT(IN),
OPTIONAL :: matrix_name
108 CHARACTER(LEN=*),
INTENT(IN) :: basis_type
111 LOGICAL,
INTENT(IN),
OPTIONAL :: calculate_forces
112 TYPE(
dbcsr_type),
OPTIONAL,
POINTER :: matrix_p
114 POINTER :: matrixkp_p
115 REAL(kind=
dp),
INTENT(IN),
OPTIONAL :: eps_filter
116 INTEGER,
INTENT(IN),
OPTIONAL :: nderivative
122 CALL build_kinetic_matrix_low(ks_env, matrix_t, matrixkp_t, matrix_name, basis_type, &
123 sab_nl, calculate_forces, matrix_p, matrixkp_p, eps_filter, natom, &
144 SUBROUTINE build_kinetic_matrix_low(ks_env, matrix_t, matrixkp_t, matrix_name, basis_type, &
145 sab_nl, calculate_forces, matrix_p, matrixkp_p, eps_filter, natom, &
152 POINTER :: matrixkp_t
153 CHARACTER(LEN=*),
INTENT(IN),
OPTIONAL :: matrix_name
154 CHARACTER(LEN=*),
INTENT(IN) :: basis_type
157 LOGICAL,
INTENT(IN),
OPTIONAL :: calculate_forces
158 TYPE(
dbcsr_type),
OPTIONAL,
POINTER :: matrix_p
160 POINTER :: matrixkp_p
161 REAL(kind=
dp),
INTENT(IN),
OPTIONAL :: eps_filter
162 INTEGER,
INTENT(IN) :: natom
163 INTEGER,
INTENT(IN),
OPTIONAL :: nderivative
165 CHARACTER(len=*),
PARAMETER :: routinen =
'build_kinetic_matrix_low'
167 INTEGER :: atom_a, handle, i, iatom, ic, icol, ikind, irow, iset, jatom, jkind, jset, ldsab, &
168 maxder, ncoa, ncob, nder, nimg, nkind, nseta, nsetb, sgfa, sgfb, slot
169 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: atom_of_kind, kind_of
170 INTEGER,
DIMENSION(3) :: cell
171 INTEGER,
DIMENSION(:),
POINTER :: la_max, la_min, lb_max, lb_min, npgfa, &
173 INTEGER,
DIMENSION(:, :),
POINTER :: first_sgfa, first_sgfb
174 INTEGER,
DIMENSION(:, :, :),
POINTER :: cell_to_index
175 LOGICAL :: do_forces, do_symmetric, dokp, found, &
176 trans, use_cell_mapping, use_virial
177 REAL(kind=
dp) :: f, f0, ff, rab2, tab
178 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: pab, qab
179 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: kab
180 REAL(kind=
dp),
DIMENSION(3) :: force_a, rab
181 REAL(kind=
dp),
DIMENSION(3, 3) :: pv_thread
182 REAL(kind=
dp),
DIMENSION(3, natom) :: force_thread
183 REAL(kind=
dp),
DIMENSION(:),
POINTER :: set_radius_a, set_radius_b
184 REAL(kind=
dp),
DIMENSION(:, :),
POINTER :: p_block, rpgfa, rpgfb, scon_a, scon_b, &
186 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: dab
188 TYPE(
block_p_type),
ALLOCATABLE,
DIMENSION(:) :: k_block
194 TYPE(
qs_kind_type),
DIMENSION(:),
POINTER :: qs_kind_set
205 CALL timeset(routinen, handle)
208 IF (
PRESENT(matrix_t))
THEN
210 use_cell_mapping = .false.
211 ELSEIF (
PRESENT(matrixkp_t))
THEN
213 CALL get_ks_env(ks_env=ks_env, kpoints=kpoints)
215 use_cell_mapping = (
SIZE(cell_to_index) > 1)
220 NULLIFY (atomic_kind_set, qs_kind_set, p_block, dft_control)
222 dft_control=dft_control, &
223 atomic_kind_set=atomic_kind_set, &
224 qs_kind_set=qs_kind_set)
226 nimg = dft_control%nimages
227 nkind =
SIZE(atomic_kind_set)
230 IF (
PRESENT(calculate_forces)) do_forces = calculate_forces
233 IF (
PRESENT(nderivative)) nder = nderivative
237 cpassert(
SIZE(sab_nl) > 0)
241 ALLOCATE (basis_set_list(nkind))
246 CALL create_sab_matrix(ks_env, matrixkp_t, matrix_name, basis_set_list, basis_set_list, &
247 sab_nl, do_symmetric)
250 CALL create_sab_matrix(ks_env, matrix_t, matrix_name, basis_set_list, basis_set_list, &
251 sab_nl, do_symmetric)
257 CALL get_ks_env(ks_env=ks_env, force=force, virial=virial)
258 use_virial = virial%pv_availability .AND. (.NOT. virial%pv_numer)
261 cpassert(
PRESENT(matrixkp_p))
263 cpassert(
PRESENT(matrix_p))
267 force_thread = 0.0_dp
294 ALLOCATE (kab(ldsab, ldsab, maxder), qab(ldsab, ldsab))
296 ALLOCATE (dab(ldsab, ldsab, 3), pab(ldsab, ldsab))
299 ALLOCATE (k_block(maxder))
301 NULLIFY (k_block(i)%block)
305 DO slot = 1, sab_nl(1)%nl_size
307 ikind = sab_nl(1)%nlist_task(slot)%ikind
308 jkind = sab_nl(1)%nlist_task(slot)%jkind
309 iatom = sab_nl(1)%nlist_task(slot)%iatom
310 jatom = sab_nl(1)%nlist_task(slot)%jatom
311 cell(:) = sab_nl(1)%nlist_task(slot)%cell(:)
312 rab(1:3) = sab_nl(1)%nlist_task(slot)%r(1:3)
314 basis_set_a => basis_set_list(ikind)%gto_basis_set
315 IF (.NOT.
ASSOCIATED(basis_set_a)) cycle
316 basis_set_b => basis_set_list(jkind)%gto_basis_set
317 IF (.NOT.
ASSOCIATED(basis_set_b)) cycle
323 first_sgfa => basis_set_a%first_sgf
324 la_max => basis_set_a%lmax
325 la_min => basis_set_a%lmin
326 npgfa => basis_set_a%npgf
327 nseta = basis_set_a%nset
328 nsgfa => basis_set_a%nsgf_set
329 rpgfa => basis_set_a%pgf_radius
330 set_radius_a => basis_set_a%set_radius
331 scon_a => basis_set_a%scon
332 zeta => basis_set_a%zet
334 first_sgfb => basis_set_b%first_sgf
335 lb_max => basis_set_b%lmax
336 lb_min => basis_set_b%lmin
337 npgfb => basis_set_b%npgf
338 nsetb = basis_set_b%nset
339 nsgfb => basis_set_b%nsgf_set
340 rpgfb => basis_set_b%pgf_radius
341 set_radius_b => basis_set_b%set_radius
342 scon_b => basis_set_b%scon
343 zetb => basis_set_b%zet
345 IF (use_cell_mapping)
THEN
346 ic = cell_to_index(cell(1), cell(2), cell(3))
352 IF (do_symmetric)
THEN
353 IF (iatom <= jatom)
THEN
361 IF (iatom == jatom) f0 = 1.0_dp
371 row=irow, col=icol, block=k_block(1)%block, found=found)
375 NULLIFY (k_block(i)%block)
377 row=irow, col=icol, block=k_block(i)%block, found=found)
386 row=irow, col=icol, block=p_block, found=found)
390 block=p_block, found=found)
395 rab2 = rab(1)*rab(1) + rab(2)*rab(2) + rab(3)*rab(3)
397 trans = do_symmetric .AND. (iatom > jatom)
401 ncoa = npgfa(iset)*(
ncoset(la_max(iset)) -
ncoset(la_min(iset) - 1))
402 sgfa = first_sgfa(1, iset)
406 IF (set_radius_a(iset) + set_radius_b(jset) < tab) cycle
411 ncob = npgfb(jset)*(
ncoset(lb_max(jset)) -
ncoset(lb_min(jset) - 1))
412 sgfb = first_sgfb(1, jset)
414 IF (do_forces .AND.
ASSOCIATED(p_block) .AND. ((iatom /= jatom) .OR. use_virial))
THEN
417 CALL block_add(
"OUT", kab(:, :, 1), nsgfa(iset), nsgfb(jset), p_block, sgfa, sgfb, trans=trans)
418 CALL decontraction(kab(:, :, 1), pab, scon_a(:, sgfa:), ncoa, nsgfa(iset), scon_b(:, sgfb:), ncob, nsgfb(jset), &
421 CALL kinetic(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
422 lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
423 rab, kab(:, :, 1), dab)
425 force_thread(:, iatom) = force_thread(:, iatom) + ff*force_a(:)
426 force_thread(:, jatom) = force_thread(:, jatom) - ff*force_a(:)
433 CALL kinetic(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
434 lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
435 rab, kab=kab(:, :, 1))
436 ELSE IF (nder == 1)
THEN
437 CALL kinetic(la_max(iset), la_min(iset), npgfa(iset), rpgfa(:, iset), zeta(:, iset), &
438 lb_max(jset), lb_min(jset), npgfb(jset), rpgfb(:, jset), zetb(:, jset), &
439 rab, kab=kab(:, :, 1), dab=kab(:, :, 2:4))
444 IF (ndod(i) == 1 .AND. trans) f = -1.0_dp
446 CALL contraction(kab(:, :, i), qab, ca=scon_a(:, sgfa:), na=ncoa, ma=nsgfa(iset), &
447 cb=scon_b(:, sgfb:), nb=ncob, mb=nsgfb(jset), fscale=f, &
450 CALL block_add(
"IN", qab, nsgfa(iset), nsgfb(jset), k_block(i)%block, sgfa, sgfb, trans=trans)
457 DEALLOCATE (kab, qab)
458 IF (do_forces)
DEALLOCATE (pab, dab)
477 atom_a = atom_of_kind(iatom)
478 ikind = kind_of(iatom)
479 force(ikind)%kinetic(:, atom_a) = force(ikind)%kinetic(:, atom_a) + force_thread(:, iatom)
483 IF (do_forces .AND. use_virial)
THEN
484 virial%pv_ekinetic = virial%pv_ekinetic + pv_thread
485 virial%pv_virial = virial%pv_virial + pv_thread
491 IF (
PRESENT(eps_filter))
THEN
497 IF (
PRESENT(eps_filter))
THEN
502 DEALLOCATE (basis_set_list)
504 CALL timestop(handle)
506 END SUBROUTINE build_kinetic_matrix_low
Set of routines to: Contract integrals over primitive Gaussians Decontract (density) matrices Trace m...
Calculation of the kinetic energy integrals over Cartesian Gaussian-type functions.
subroutine, public kinetic(la_max, la_min, npgfa, rpgfa, zeta, lb_max, lb_min, npgfb, rpgfb, zetb, rab, kab, dab)
Calculation of the two-center kinetic energy integrals [a|T|b] over Cartesian Gaussian-type functions...
Define the atomic kind types and their sub types.
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.
collect pointers to a block of reals
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
subroutine, public dbcsr_get_block_p(matrix, row, col, block, found, row_size, col_size)
...
subroutine, public dbcsr_filter(matrix, eps)
...
subroutine, public dbcsr_finalize(matrix)
...
DBCSR operations in CP2K.
Defines the basic variable types.
integer, parameter, public int_8
integer, parameter, public dp
Types and basic routines needed for a kpoint calculation.
subroutine, public get_kpoint_info(kpoint, kp_scheme, nkp_grid, kp_shift, symmetry, verbose, full_grid, use_real_wfn, eps_geo, parallel_group_size, kp_range, nkp, xkp, wkp, para_env, blacs_env_all, para_env_kp, para_env_inter_kp, blacs_env, kp_env, kp_aux_env, mpools, iogrp, nkp_groups, kp_dist, cell_to_index, index_to_cell, sab_nl, sab_nl_nosym)
Retrieve information from a kpoint environment.
Provides Cartesian and spherical orbital pointers and indices.
integer, dimension(:), allocatable, public ncoset
Some utility functions for the calculation of integrals.
subroutine, public basis_set_list_setup(basis_set_list, basis_type, qs_kind_set)
Set up an easy accessible list of the basis sets for all kinds.
Define the quickstep kind type and their sub types.
Calculation of kinetic energy matrix and forces.
subroutine, public build_kinetic_matrix(ks_env, matrix_t, matrixkp_t, matrix_name, basis_type, sab_nl, calculate_forces, matrix_p, matrixkp_p, eps_filter, nderivative)
Calculation of the kinetic energy matrix over Cartesian Gaussian functions.
subroutine, public get_ks_env(ks_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, complex_ks, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, kinetic, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_ks_im_kp, rho, rho_xc, vppl, rho_core, rho_nlcc, rho_nlcc_g, vee, neighbor_list_id, sab_orb, sab_all, sac_ae, sac_ppl, sac_lri, sap_ppnl, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_vdw, sab_scp, sab_almo, sab_kp, sab_kp_nosym, task_list, task_list_soft, kpoints, do_kpoints, atomic_kind_set, qs_kind_set, cell, cell_ref, use_ref_cell, particle_set, energy, force, local_particles, local_molecules, molecule_kind_set, molecule_set, subsys, cp_subsys, virial, results, atprop, nkind, natom, dft_control, dbcsr_dist, distribution_2d, pw_env, para_env, blacs_env, nelectron_total, nelectron_spin)
...
Define the neighbor list data types and the corresponding functionality.
subroutine, public get_neighbor_list_set_p(neighbor_list_sets, nlist, symmetric)
Return the components of the first neighbor list set.
Calculation of overlap matrix, its derivatives and forces.
pure subroutine, public virial_pair_force(pv_virial, f0, force, rab)
Computes the contribution to the stress tensor from two-body pair-wise forces.
Provides all information about an atomic kind.
Contains information about kpoints.
Provides all information about a quickstep kind.
calculation environment to calculate the ks matrix, holds all the needed vars. assumes that the core ...