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 Main setup file for RI-RS grids {r_l}.
10 : !> \par History
11 : !> 09.2026 created
12 : ! **************************************************************************************************
13 : MODULE gw_ri_rs_grid_setup_main
14 : USE gw_ri_rs_grid_from_file, ONLY: read_ri_rs_grid_from_file
15 : USE gw_ri_rs_grid_optimization, ONLY: optimize_ri_rs_grid
16 : USE kinds, ONLY: dp,&
17 : int_8
18 : USE particle_types, ONLY: particle_type
19 : USE post_scf_bandstructure_types, ONLY: post_scf_bandstructure_type
20 : USE util, ONLY: sort
21 : #include "./base/base_uses.f90"
22 :
23 : IMPLICIT NONE
24 : PRIVATE
25 :
26 : PUBLIC :: setup_ri_rs_grid
27 :
28 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gw_ri_rs_grid_setup_main'
29 :
30 : CONTAINS
31 :
32 : ! **************************************************************************************************
33 : !> \brief Get RI-RS grid points {r_l}, either by on-the-fly optimization or
34 : !> reading pretabulated atomic grids
35 : !> \param bs_env Band-structure environment containing GW RI-RS parameters.
36 : !> \param grid_points x,y,z RI-RS grid coordinates, size (3, ngrid)
37 : ! **************************************************************************************************
38 48 : SUBROUTINE setup_ri_rs_grid(bs_env, grid_points)
39 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
40 : REAL(KIND=dp), ALLOCATABLE, INTENT(OUT) :: grid_points(:, :)
41 :
42 48 : IF (bs_env%ri_rs%grid_opt%enabled) THEN
43 :
44 : ! Initialize Lebedev grids for every element from CP2K routines and then
45 : ! optimize the coordinates of these grid points to minimize the RIRS error
46 : !
47 : ! [(μν|P) - Σ_l ϕ_μ(r_l) ϕ_ν(r_l) Z_lP^(A)]²
48 : !
49 10 : CALL optimize_ri_rs_grid(bs_env)
50 :
51 : ELSE
52 :
53 : ! Read pretabulated atom-relative RI-RS grids from the data files
54 38 : CALL read_ri_rs_grid_from_file(bs_env)
55 :
56 : END IF
57 :
58 : ! Move atom-relative grids to their atomic centres and collect the global grid.
59 48 : CALL assemble_ri_rs_grid(bs_env, grid_points)
60 :
61 48 : END SUBROUTINE setup_ri_rs_grid
62 :
63 : ! **************************************************************************************************
64 : !> \brief Move atom-relative RI-RS grids to their atomic centres and assemble the global grid.
65 : !> \param bs_env ...
66 : !> \param ri_rs_grid_points ...
67 : ! **************************************************************************************************
68 48 : SUBROUTINE assemble_ri_rs_grid(bs_env, ri_rs_grid_points)
69 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
70 : REAL(KIND=dp), ALLOCATABLE, INTENT(OUT) :: ri_rs_grid_points(:, :)
71 :
72 : CHARACTER(LEN=*), PARAMETER :: routineN = 'assemble_ri_rs_grid'
73 :
74 : INTEGER :: handle
75 :
76 48 : CALL timeset(routineN, handle)
77 :
78 48 : CPASSERT(ALLOCATED(bs_env%ri_rs%atomic_grids))
79 48 : CALL assemble_ri_rs_grid_points(bs_env, ri_rs_grid_points)
80 48 : CALL release_atomic_grids(bs_env)
81 :
82 48 : CALL timestop(handle)
83 :
84 48 : END SUBROUTINE assemble_ri_rs_grid
85 :
86 : ! **************************************************************************************************
87 : !> \brief Assemble the global RI-RS grid in spatial atom order from atom-relative grids.
88 : !> \param bs_env ...
89 : !> \param ri_rs_grid_points ...
90 : ! **************************************************************************************************
91 48 : SUBROUTINE assemble_ri_rs_grid_points(bs_env, ri_rs_grid_points)
92 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
93 : REAL(KIND=dp), ALLOCATABLE, INTENT(OUT) :: ri_rs_grid_points(:, :)
94 :
95 : INTEGER :: atom_grid_end, atom_grid_start, iatom, &
96 : ilayout, natom
97 48 : INTEGER, ALLOCATABLE :: atom_grid_offsets(:), atom_order(:)
98 : REAL(KIND=dp) :: atom_center(3)
99 :
100 48 : natom = bs_env%n_atom
101 48 : CALL spatial_atom_order(bs_env%ri_rs%particle_set, atom_order)
102 :
103 240 : ALLOCATE (bs_env%ri_rs%grid_atom_boundaries(natom + 1), atom_grid_offsets(natom))
104 48 : bs_env%ri_rs%n_grid_points = 0
105 166 : DO ilayout = 1, natom
106 118 : iatom = atom_order(ilayout)
107 118 : atom_grid_offsets(iatom) = bs_env%ri_rs%n_grid_points + 1
108 118 : bs_env%ri_rs%grid_atom_boundaries(ilayout) = bs_env%ri_rs%n_grid_points + 1
109 : bs_env%ri_rs%n_grid_points = bs_env%ri_rs%n_grid_points + &
110 166 : bs_env%ri_rs%atomic_grids(iatom)%npts
111 : END DO
112 48 : bs_env%ri_rs%grid_atom_boundaries(natom + 1) = bs_env%ri_rs%n_grid_points + 1
113 :
114 48 : IF (bs_env%unit_nr > 0) THEN
115 : WRITE (bs_env%unit_nr, FMT="(T2,A,T69,I12)") &
116 24 : 'Total grid points used for RI-RS:', bs_env%ri_rs%n_grid_points
117 24 : WRITE (bs_env%unit_nr, "(A)") ' '
118 : END IF
119 :
120 144 : ALLOCATE (ri_rs_grid_points(3, bs_env%ri_rs%n_grid_points))
121 : !$OMP PARALLEL DO DEFAULT(NONE) &
122 : !$OMP SHARED(ri_rs_grid_points, atom_grid_offsets, bs_env, natom) &
123 : !$OMP PRIVATE(iatom, atom_center, atom_grid_start, atom_grid_end) &
124 48 : !$OMP SCHEDULE(DYNAMIC, 1)
125 : DO iatom = 1, natom
126 : atom_center(:) = bs_env%ri_rs%particle_set(iatom)%r(:)
127 : atom_grid_start = atom_grid_offsets(iatom)
128 : atom_grid_end = atom_grid_start + bs_env%ri_rs%atomic_grids(iatom)%npts - 1
129 :
130 : ri_rs_grid_points(1, atom_grid_start:atom_grid_end) = &
131 : bs_env%ri_rs%atomic_grids(iatom)%raw_points(1, :) + atom_center(1)
132 : ri_rs_grid_points(2, atom_grid_start:atom_grid_end) = &
133 : bs_env%ri_rs%atomic_grids(iatom)%raw_points(2, :) + atom_center(2)
134 : ri_rs_grid_points(3, atom_grid_start:atom_grid_end) = &
135 : bs_env%ri_rs%atomic_grids(iatom)%raw_points(3, :) + atom_center(3)
136 : END DO
137 : !$OMP END PARALLEL DO
138 48 : END SUBROUTINE assemble_ri_rs_grid_points
139 :
140 : ! **************************************************************************************************
141 : !> \brief Release the atom-relative RI-RS grids after assembling the global molecular grid.
142 : !> \param bs_env ...
143 : ! **************************************************************************************************
144 48 : SUBROUTINE release_atomic_grids(bs_env)
145 : TYPE(post_scf_bandstructure_type), POINTER :: bs_env
146 :
147 : INTEGER :: iatom
148 :
149 166 : DO iatom = 1, SIZE(bs_env%ri_rs%atomic_grids)
150 166 : DEALLOCATE (bs_env%ri_rs%atomic_grids(iatom)%raw_points)
151 : END DO
152 166 : DEALLOCATE (bs_env%ri_rs%atomic_grids)
153 48 : END SUBROUTINE release_atomic_grids
154 :
155 : ! **************************************************************************************************
156 : !> \brief Order atoms by a three-dimensional Morton code so consecutive grid runs remain local.
157 : !> \param particle_set ...
158 : !> \param order ...
159 : ! **************************************************************************************************
160 48 : SUBROUTINE spatial_atom_order(particle_set, order)
161 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
162 : INTEGER, ALLOCATABLE, DIMENSION(:), INTENT(OUT) :: order
163 :
164 : CHARACTER(LEN=*), PARAMETER :: routineN = 'spatial_atom_order'
165 : INTEGER, PARAMETER :: nbits = 21
166 :
167 : INTEGER :: handle, iatom, idimension, natom
168 : INTEGER(KIND=int_8) :: cmax, integer_coordinate(3), m1, m2, m3
169 48 : INTEGER(KIND=int_8), ALLOCATABLE :: morton_code(:)
170 : REAL(KIND=dp) :: hi(3), lo(3), span(3)
171 :
172 48 : CALL timeset(routineN, handle)
173 :
174 48 : natom = SIZE(particle_set)
175 240 : ALLOCATE (order(natom), morton_code(natom))
176 : cmax = ISHFT(1_int_8, nbits) - 1_int_8
177 :
178 192 : lo(:) = HUGE(1.0_dp)
179 192 : hi(:) = -HUGE(1.0_dp)
180 166 : DO iatom = 1, natom
181 520 : DO idimension = 1, 3
182 354 : lo(idimension) = MIN(lo(idimension), particle_set(iatom)%r(idimension))
183 472 : hi(idimension) = MAX(hi(idimension), particle_set(iatom)%r(idimension))
184 : END DO
185 : END DO
186 192 : span(:) = hi(:) - lo(:)
187 192 : DO idimension = 1, 3
188 192 : IF (span(idimension) <= 0.0_dp) span(idimension) = 1.0_dp
189 : END DO
190 :
191 166 : DO iatom = 1, natom
192 472 : DO idimension = 1, 3
193 : integer_coordinate(idimension) = &
194 : INT(((particle_set(iatom)%r(idimension) - lo(idimension))/span(idimension))* &
195 354 : REAL(cmax, dp), int_8)
196 : integer_coordinate(idimension) = &
197 472 : MIN(cmax, MAX(0_int_8, integer_coordinate(idimension)))
198 : END DO
199 118 : CALL morton_split3(integer_coordinate(1), m1)
200 118 : CALL morton_split3(integer_coordinate(2), m2)
201 118 : CALL morton_split3(integer_coordinate(3), m3)
202 166 : morton_code(iatom) = IOR(IOR(m1, ISHFT(m2, 1)), ISHFT(m3, 2))
203 : END DO
204 :
205 48 : CALL sort(morton_code, natom, order)
206 48 : DEALLOCATE (morton_code)
207 :
208 48 : CALL timestop(handle)
209 :
210 48 : END SUBROUTINE spatial_atom_order
211 :
212 : ! **************************************************************************************************
213 : !> \brief Spread the low 21 bits of an integer over every third bit of a Morton code.
214 : !> \param input_integer ...
215 : !> \param spread_integer ...
216 : ! **************************************************************************************************
217 354 : SUBROUTINE morton_split3(input_integer, spread_integer)
218 : INTEGER(KIND=int_8), INTENT(IN) :: input_integer
219 : INTEGER(KIND=int_8), INTENT(OUT) :: spread_integer
220 :
221 354 : spread_integer = IAND(input_integer, INT(z'1FFFFF', int_8))
222 : spread_integer = IAND(IOR(spread_integer, ISHFT(spread_integer, 32)), &
223 354 : INT(z'1F00000000FFFF', int_8))
224 : spread_integer = IAND(IOR(spread_integer, ISHFT(spread_integer, 16)), &
225 354 : INT(z'1F0000FF0000FF', int_8))
226 : spread_integer = IAND(IOR(spread_integer, ISHFT(spread_integer, 8)), &
227 354 : INT(z'100F00F00F00F00F', int_8))
228 : spread_integer = IAND(IOR(spread_integer, ISHFT(spread_integer, 4)), &
229 354 : INT(z'10C30C30C30C30C3', int_8))
230 : spread_integer = IAND(IOR(spread_integer, ISHFT(spread_integer, 2)), &
231 354 : INT(z'1249249249249249', int_8))
232 354 : END SUBROUTINE morton_split3
233 :
234 : END MODULE gw_ri_rs_grid_setup_main
|
