Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Analytic translations of finite or periodized Gaussian AO functions.
10 : !> \author Thomas D. Kuehne, tkuehne@cp2k.org
11 : ! **************************************************************************************************
12 : MODULE qs_ao_translations
13 : USE ai_contraction, ONLY: block_add,&
14 : contraction
15 : USE ai_overlap, ONLY: overlap_ab
16 : USE basis_set_types, ONLY: gto_basis_set_p_type,&
17 : gto_basis_set_type
18 : USE cell_types, ONLY: cell_type
19 : USE cp_dbcsr_api, ONLY: dbcsr_create,&
20 : dbcsr_distribution_get,&
21 : dbcsr_distribution_type,&
22 : dbcsr_finalize,&
23 : dbcsr_get_info,&
24 : dbcsr_p_type,&
25 : dbcsr_put_block,&
26 : dbcsr_type_no_symmetry
27 : USE kinds, ONLY: dp
28 : USE message_passing, ONLY: mp_para_env_type
29 : USE orbital_pointers, ONLY: ncoset
30 : USE particle_types, ONLY: particle_type
31 : USE qs_environment_types, ONLY: get_qs_env,&
32 : qs_environment_type
33 : USE qs_integral_utils, ONLY: basis_set_list_setup,&
34 : get_memory_usage
35 : USE qs_kind_types, ONLY: qs_kind_type
36 : #include "./base/base_uses.f90"
37 :
38 : IMPLICIT NONE
39 : PRIVATE
40 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_ao_translations'
41 : PUBLIC :: ao_translation_matrices
42 : CONTAINS
43 : ! **************************************************************************************************
44 : !> \brief Overlaps with translated ket Gaussians, U(a) psi(r) = psi(r+a).
45 : !> \param qs_env Quickstep environment supplying primitive AOs and their distribution
46 : !> \param shifts Cartesian translation vectors in internal bohr
47 : !> \param translations Newly allocated nonsymmetric primitive-size matrices, released by caller
48 : !> \param periodic_grid If present, sum images of the grid-expanded periodic cell; otherwise finite
49 : !> \note Arbitrary shifts use their own overlap screening, not the SCF neighbour list.
50 : ! **************************************************************************************************
51 16 : SUBROUTINE ao_translation_matrices(qs_env, shifts, translations, periodic_grid)
52 : TYPE(qs_environment_type), POINTER :: qs_env
53 : REAL(dp), INTENT(IN) :: shifts(:, :)
54 : TYPE(dbcsr_p_type), POINTER :: translations(:)
55 : INTEGER, INTENT(IN), OPTIONAL :: periodic_grid(3)
56 :
57 : INTEGER :: axis, ia, ib, iset, it, jset, ka, kb, &
58 : ld, lower(3), na, nb, nx, ny, nz, sa, &
59 : sb, upper(3)
60 16 : INTEGER, POINTER :: col_dist(:), pgrid(:, :), row_dist(:)
61 : LOGICAL :: do_kpoints
62 : REAL(dp) :: distance, h(3, 3), hinv(3, 3), &
63 : offset(3), rab(3), radius, reach(3), &
64 : scaled(3)
65 16 : REAL(dp), ALLOCATABLE :: block(:, :), integrals(:, :), work(:, :)
66 : TYPE(cell_type), POINTER :: cell
67 : TYPE(dbcsr_distribution_type) :: distribution
68 16 : TYPE(dbcsr_p_type), POINTER :: matrix_s(:, :)
69 16 : TYPE(gto_basis_set_p_type), ALLOCATABLE :: bases(:)
70 : TYPE(gto_basis_set_type), POINTER :: a, b
71 : TYPE(mp_para_env_type), POINTER :: para_env
72 16 : TYPE(particle_type), POINTER :: particles(:)
73 16 : TYPE(qs_kind_type), POINTER :: kinds(:)
74 :
75 : CALL get_qs_env(qs_env, cell=cell, do_kpoints=do_kpoints, matrix_s_kp=matrix_s, &
76 16 : para_env=para_env, particle_set=particles, qs_kind_set=kinds)
77 16 : h = 0.0_dp
78 16 : hinv = 0.0_dp
79 16 : reach = 0.0_dp
80 16 : IF (PRESENT(periodic_grid)) THEN
81 56 : CPASSERT(ALL(periodic_grid >= 1))
82 56 : CPASSERT(ALL(periodic_grid == 1 .OR. cell%perd /= 0))
83 56 : DO axis = 1, 3
84 168 : h(:, axis) = cell%hmat(:, axis)*periodic_grid(axis)
85 168 : hinv(axis, :) = cell%h_inv(axis, :)/periodic_grid(axis)
86 182 : IF (cell%perd(axis) /= 0) reach(axis) = SQRT(SUM(hinv(axis, :)**2))
87 : END DO
88 8 : ELSE IF (ANY(cell%perd /= 0) .OR. do_kpoints) THEN
89 0 : CPABORT('Finite AO translations require isolated Gamma systems.')
90 : END IF
91 16 : CPASSERT(SIZE(shifts, 1) == 3)
92 64 : ALLOCATE (bases(SIZE(kinds)))
93 16 : CALL basis_set_list_setup(bases, 'ORB', kinds)
94 16 : ld = get_memory_usage(kinds, 'ORB')
95 146 : ALLOCATE (integrals(ld, ld), work(ld, ld), translations(SIZE(shifts, 2)))
96 16 : CALL dbcsr_get_info(matrix_s(1, 1)%matrix, distribution=distribution)
97 16 : CALL dbcsr_distribution_get(distribution, pgrid=pgrid, row_dist=row_dist, col_dist=col_dist)
98 34 : DO it = 1, SIZE(shifts, 2)
99 18 : ALLOCATE (translations(it)%matrix)
100 : CALL dbcsr_create(translations(it)%matrix, template=matrix_s(1, 1)%matrix, &
101 18 : name='AO TRANSLATION', matrix_type=dbcsr_type_no_symmetry)
102 44 : DO ia = 1, SIZE(particles)
103 26 : ka = particles(ia)%atomic_kind%kind_number
104 26 : a => bases(ka)%gto_basis_set
105 26 : IF (.NOT. ASSOCIATED(a)) CYCLE
106 90 : DO ib = 1, SIZE(particles)
107 46 : IF (pgrid(row_dist(ia), col_dist(ib)) /= para_env%mepos) CYCLE
108 23 : kb = particles(ib)%atomic_kind%kind_number
109 23 : b => bases(kb)%gto_basis_set
110 23 : IF (.NOT. ASSOCIATED(b)) CYCLE
111 92 : offset = particles(ib)%r - particles(ia)%r - shifts(:, it)
112 23 : lower = 0
113 23 : upper = 0
114 23 : IF (PRESENT(periodic_grid)) THEN
115 : ! Reciprocal-plane distances bound every image, including skew cells.
116 195 : scaled = MATMUL(hinv, offset)
117 75 : radius = MAXVAL(a%set_radius) + MAXVAL(b%set_radius)
118 60 : DO axis = 1, 3
119 45 : IF (cell%perd(axis) == 0) CYCLE
120 45 : lower(axis) = CEILING(-scaled(axis) - radius*reach(axis))
121 60 : upper(axis) = FLOOR(-scaled(axis) + radius*reach(axis))
122 : END DO
123 : END IF
124 92 : ALLOCATE (block(a%nsgf, b%nsgf))
125 23 : block = 0.0_dp
126 106 : DO nz = lower(3), upper(3)
127 489 : DO ny = lower(2), upper(2)
128 1099 : DO nx = lower(1), upper(1)
129 12027 : rab = offset + MATMUL(h, REAL([nx, ny, nz], dp))
130 2532 : distance = SQRT(SUM(rab**2))
131 2274 : DO iset = 1, a%nset
132 1258 : sa = a%first_sgf(1, iset)
133 1258 : na = a%npgf(iset)*(ncoset(a%lmax(iset)) - ncoset(a%lmin(iset) - 1))
134 4399 : DO jset = 1, b%nset
135 2508 : IF (a%set_radius(iset) + b%set_radius(jset) < distance) CYCLE
136 1540 : sb = b%first_sgf(1, jset)
137 1540 : nb = b%npgf(jset)*(ncoset(b%lmax(jset)) - ncoset(b%lmin(jset) - 1))
138 : CALL overlap_ab(a%lmax(iset), a%lmin(iset), a%npgf(iset), &
139 : a%pgf_radius(:, iset), a%zet(:, iset), &
140 : b%lmax(jset), b%lmin(jset), b%npgf(jset), &
141 1540 : b%pgf_radius(:, jset), b%zet(:, jset), rab, sab=integrals)
142 : CALL contraction(integrals, work, ca=a%scon(:, sa:), na=na, ma=a%nsgf_set(iset), &
143 : cb=b%scon(:, sb:), nb=nb, mb=b%nsgf_set(jset), &
144 1540 : fscale=1.0_dp, trans=.FALSE.)
145 : CALL block_add('IN', work, a%nsgf_set(iset), b%nsgf_set(jset), &
146 3766 : block, sa, sb, trans=.FALSE.)
147 : END DO
148 : END DO
149 : END DO
150 : END DO
151 : END DO
152 989 : IF (MAXVAL(ABS(block)) > 0.0_dp) CALL dbcsr_put_block(translations(it)%matrix, ia, ib, block)
153 72 : DEALLOCATE (block)
154 : END DO
155 : END DO
156 34 : CALL dbcsr_finalize(translations(it)%matrix)
157 : END DO
158 48 : END SUBROUTINE ao_translation_matrices
159 : END MODULE qs_ao_translations
|
