Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief routines that build the Kohn-Sham matrix for the LRIGPW
10 : !> and xc parts
11 : !> \par History
12 : !> 09.2013 created [Dorothea Golze]
13 : !> \author Dorothea Golze
14 : ! **************************************************************************************************
15 : MODULE lri_ks_methods
16 : USE atomic_kind_types, ONLY: atomic_kind_type,&
17 : get_atomic_kind_set
18 : USE cp_dbcsr_api, ONLY: dbcsr_add,&
19 : dbcsr_finalize,&
20 : dbcsr_get_block_p,&
21 : dbcsr_p_type,&
22 : dbcsr_reserve_blocks,&
23 : dbcsr_type
24 : USE kinds, ONLY: dp
25 : USE lri_compression, ONLY: lri_decomp_i
26 : USE lri_environment_types, ONLY: lri_environment_type,&
27 : lri_int_type,&
28 : lri_kind_type
29 : USE qs_neighbor_list_types, ONLY: get_iterator_info,&
30 : neighbor_list_iterate,&
31 : neighbor_list_iterator_create,&
32 : neighbor_list_iterator_p_type,&
33 : neighbor_list_iterator_release,&
34 : neighbor_list_set_p_type
35 : USE qs_o3c_methods, ONLY: contract3_o3c
36 : USE qs_o3c_types, ONLY: get_o3c_vec,&
37 : o3c_vec_create,&
38 : o3c_vec_release,&
39 : o3c_vec_type
40 : USE ri_environment_methods, ONLY: ri_metric_solver
41 :
42 : !$ USE OMP_LIB, ONLY: omp_get_max_threads, omp_get_thread_num
43 : #include "./base/base_uses.f90"
44 :
45 : IMPLICIT NONE
46 :
47 : PRIVATE
48 :
49 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'lri_ks_methods'
50 :
51 : PUBLIC :: calculate_lri_ks_matrix, calculate_ri_ks_matrix
52 :
53 : CONTAINS
54 :
55 : !*****************************************************************************
56 : !> \brief update of LRIGPW KS matrix
57 : !> \param lri_env ...
58 : !> \param lri_v_int integrals of potential * ri basis set
59 : !> \param h_matrix KS matrix, on entry containing the core hamiltonian
60 : !> \param atomic_kind_set ...
61 : !> \param cell_to_index ...
62 : !> \note including this in lri_environment_methods?
63 : ! **************************************************************************************************
64 724 : SUBROUTINE calculate_lri_ks_matrix(lri_env, lri_v_int, h_matrix, atomic_kind_set, cell_to_index)
65 :
66 : TYPE(lri_environment_type), POINTER :: lri_env
67 : TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_v_int
68 : TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: h_matrix
69 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
70 : INTEGER, DIMENSION(:, :, :), OPTIONAL, POINTER :: cell_to_index
71 :
72 : CHARACTER(*), PARAMETER :: routineN = 'calculate_lri_ks_matrix'
73 :
74 : INTEGER :: atom_a, atom_b, col, handle, i, iac, iatom, ic, ikind, ilist, jatom, jkind, &
75 : jneighbor, mepos, nba, nbb, nfa, nfb, nkind, nlist, nm, nn, nthread, row
76 724 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind
77 : INTEGER, DIMENSION(3) :: cell
78 : LOGICAL :: found, trans, use_cell_mapping
79 : REAL(KIND=dp) :: dab, fw, isn, isna, isnb, rab(3), &
80 : threshold
81 724 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: vi, via, vib
82 724 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: hf_work, hs_work, int3, wab, wbb
83 724 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: h_block
84 : TYPE(dbcsr_type), POINTER :: hmat
85 : TYPE(lri_int_type), POINTER :: lrii
86 : TYPE(neighbor_list_iterator_p_type), &
87 724 : DIMENSION(:), POINTER :: nl_iterator
88 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
89 724 : POINTER :: soo_list
90 :
91 724 : CALL timeset(routineN, handle)
92 724 : NULLIFY (h_block, lrii, nl_iterator, soo_list)
93 :
94 724 : threshold = lri_env%eps_o3_int
95 :
96 724 : use_cell_mapping = (SIZE(h_matrix, 1) > 1)
97 724 : IF (use_cell_mapping) THEN
98 6 : CPASSERT(PRESENT(cell_to_index))
99 : END IF
100 :
101 724 : IF (ASSOCIATED(lri_env%soo_list)) THEN
102 724 : soo_list => lri_env%soo_list
103 :
104 724 : nkind = lri_env%lri_ints%nkind
105 : nthread = 1
106 724 : !$ nthread = omp_get_max_threads()
107 :
108 724 : CALL get_atomic_kind_set(atomic_kind_set=atomic_kind_set, atom_of_kind=atom_of_kind)
109 724 : CALL neighbor_list_iterator_create(nl_iterator, soo_list, nthread=nthread)
110 : !$OMP PARALLEL DEFAULT(NONE)&
111 : !$OMP SHARED (nthread,nl_iterator,nkind,atom_of_kind,threshold,lri_env,lri_v_int,&
112 : !$OMP h_matrix,use_cell_mapping,cell_to_index)&
113 : !$OMP PRIVATE (mepos,ikind,jkind,iatom,jatom,nlist,ilist,jneighbor,rab,iac,dab,lrii,&
114 : !$OMP nfa,nfb,nba,nbb,nn,hs_work,hf_work,h_block,row,col,trans,found,wab,wbb,&
115 724 : !$OMP atom_a,atom_b,isn,nm,vi,isna,isnb,via,vib,fw,int3,cell,ic,hmat)
116 :
117 : mepos = 0
118 : !$ mepos = omp_get_thread_num()
119 :
120 : DO WHILE (neighbor_list_iterate(nl_iterator, mepos) == 0)
121 : CALL get_iterator_info(nl_iterator, mepos=mepos, ikind=ikind, jkind=jkind, iatom=iatom, &
122 : jatom=jatom, nlist=nlist, ilist=ilist, inode=jneighbor, &
123 : r=rab, cell=cell)
124 :
125 : iac = ikind + nkind*(jkind - 1)
126 : dab = SQRT(SUM(rab*rab))
127 :
128 : IF (.NOT. ASSOCIATED(lri_env%lri_ints%lri_atom(iac)%lri_node)) CYCLE
129 : IF (lri_env%exact_1c_terms) THEN
130 : IF (iatom == jatom .AND. dab < lri_env%delta) CYCLE
131 : END IF
132 :
133 : lrii => lri_env%lri_ints%lri_atom(iac)%lri_node(ilist)%lri_int(jneighbor)
134 :
135 : nfa = lrii%nfa
136 : nfb = lrii%nfb
137 : nba = lrii%nba
138 : nbb = lrii%nbb
139 : nn = nfa + nfb
140 :
141 : atom_a = atom_of_kind(iatom)
142 : atom_b = atom_of_kind(jatom)
143 :
144 : IF (use_cell_mapping) THEN
145 : ic = cell_to_index(cell(1), cell(2), cell(3))
146 : CPASSERT(ic > 0)
147 : ELSE
148 : ic = 1
149 : END IF
150 : hmat => h_matrix(ic)%matrix
151 :
152 : ALLOCATE (int3(nba, nbb))
153 : IF (lrii%lrisr) THEN
154 : ALLOCATE (hs_work(nba, nbb))
155 : IF (iatom == jatom .AND. dab < lri_env%delta) THEN
156 : nm = nfa
157 : ALLOCATE (vi(nfa))
158 : vi(1:nfa) = lri_v_int(ikind)%v_int(atom_a, 1:nfa)
159 : ELSE
160 : nm = nn
161 : ALLOCATE (vi(nn))
162 : vi(1:nfa) = lri_v_int(ikind)%v_int(atom_a, 1:nfa)
163 : vi(nfa + 1:nn) = lri_v_int(jkind)%v_int(atom_b, 1:nfb)
164 : END IF
165 : isn = SUM(lrii%sn(1:nm)*vi(1:nm))/lrii%nsn
166 : vi(1:nm) = MATMUL(lrii%sinv(1:nm, 1:nm), vi(1:nm)) - isn*lrii%sn(1:nm)
167 : hs_work(1:nba, 1:nbb) = isn*lrii%soo(1:nba, 1:nbb)
168 : IF (iatom == jatom .AND. dab < lri_env%delta) THEN
169 : DO i = 1, nfa
170 : CALL lri_decomp_i(int3, lrii%cabai, i)
171 : hs_work(1:nba, 1:nbb) = hs_work(1:nba, 1:nbb) + vi(i)*int3(1:nba, 1:nbb)
172 : END DO
173 : ELSE
174 : DO i = 1, nfa
175 : CALL lri_decomp_i(int3, lrii%cabai, i)
176 : hs_work(1:nba, 1:nbb) = hs_work(1:nba, 1:nbb) + vi(i)*int3(1:nba, 1:nbb)
177 : END DO
178 : DO i = 1, nfb
179 : CALL lri_decomp_i(int3, lrii%cabbi, i)
180 : hs_work(1:nba, 1:nbb) = hs_work(1:nba, 1:nbb) + vi(nfa + i)*int3(1:nba, 1:nbb)
181 : END DO
182 : END IF
183 : DEALLOCATE (vi)
184 : END IF
185 :
186 : IF (lrii%lriff) THEN
187 : ALLOCATE (hf_work(nba, nbb), wab(nba, nbb), wbb(nba, nbb))
188 : wab(1:nba, 1:nbb) = lri_env%wmat(ikind, jkind)%mat(1:nba, 1:nbb)
189 : wbb(1:nba, 1:nbb) = 1.0_dp - lri_env%wmat(ikind, jkind)%mat(1:nba, 1:nbb)
190 : !
191 : ALLOCATE (via(nfa), vib(nfb))
192 : via(1:nfa) = lri_v_int(ikind)%v_int(atom_a, 1:nfa)
193 : vib(1:nfb) = lri_v_int(jkind)%v_int(atom_b, 1:nfb)
194 : !
195 : isna = SUM(lrii%sna(1:nfa)*via(1:nfa))/lrii%nsna
196 : isnb = SUM(lrii%snb(1:nfb)*vib(1:nfb))/lrii%nsnb
197 : via(1:nfa) = MATMUL(lrii%asinv(1:nfa, 1:nfa), via(1:nfa)) - isna*lrii%sna(1:nfa)
198 : vib(1:nfb) = MATMUL(lrii%bsinv(1:nfb, 1:nfb), vib(1:nfb)) - isnb*lrii%snb(1:nfb)
199 : !
200 : hf_work(1:nba, 1:nbb) = (isna*wab(1:nba, 1:nbb) + isnb*wbb(1:nba, 1:nbb))*lrii%soo(1:nba, 1:nbb)
201 : !
202 : DO i = 1, nfa
203 : IF (lrii%abascr(i) > threshold) THEN
204 : CALL lri_decomp_i(int3, lrii%cabai, i)
205 : hf_work(1:nba, 1:nbb) = hf_work(1:nba, 1:nbb) + &
206 : via(i)*int3(1:nba, 1:nbb)*wab(1:nba, 1:nbb)
207 : END IF
208 : END DO
209 : DO i = 1, nfb
210 : IF (lrii%abbscr(i) > threshold) THEN
211 : CALL lri_decomp_i(int3, lrii%cabbi, i)
212 : hf_work(1:nba, 1:nbb) = hf_work(1:nba, 1:nbb) + &
213 : vib(i)*int3(1:nba, 1:nbb)*wbb(1:nba, 1:nbb)
214 : END IF
215 : END DO
216 : !
217 : DEALLOCATE (via, vib, wab, wbb)
218 : END IF
219 : DEALLOCATE (int3)
220 :
221 : ! add h_work to core hamiltonian
222 : IF (iatom <= jatom) THEN
223 : row = iatom
224 : col = jatom
225 : trans = .FALSE.
226 : ELSE
227 : row = jatom
228 : col = iatom
229 : trans = .TRUE.
230 : END IF
231 : !$OMP CRITICAL(addhamiltonian)
232 : NULLIFY (h_block)
233 : CALL dbcsr_get_block_p(hmat, row, col, h_block, found)
234 : IF (.NOT. ASSOCIATED(h_block)) THEN
235 : CALL dbcsr_reserve_blocks(hmat, rows=[row], cols=[col])
236 : CALL dbcsr_get_block_p(hmat, row, col, h_block, found)
237 : END IF
238 : IF (lrii%lrisr) THEN
239 : fw = lrii%wsr
240 : IF (trans) THEN
241 : h_block(1:nbb, 1:nba) = h_block(1:nbb, 1:nba) + fw*TRANSPOSE(hs_work(1:nba, 1:nbb))
242 : ELSE
243 : h_block(1:nba, 1:nbb) = h_block(1:nba, 1:nbb) + fw*hs_work(1:nba, 1:nbb)
244 : END IF
245 : END IF
246 : IF (lrii%lriff) THEN
247 : fw = lrii%wff
248 : IF (trans) THEN
249 : h_block(1:nbb, 1:nba) = h_block(1:nbb, 1:nba) + fw*TRANSPOSE(hf_work(1:nba, 1:nbb))
250 : ELSE
251 : h_block(1:nba, 1:nbb) = h_block(1:nba, 1:nbb) + fw*hf_work(1:nba, 1:nbb)
252 : END IF
253 : END IF
254 : !$OMP END CRITICAL(addhamiltonian)
255 :
256 : IF (lrii%lrisr) DEALLOCATE (hs_work)
257 : IF (lrii%lriff) DEALLOCATE (hf_work)
258 : END DO
259 : !$OMP END PARALLEL
260 :
261 3590 : DO ic = 1, SIZE(h_matrix, 1)
262 3590 : CALL dbcsr_finalize(h_matrix(ic)%matrix)
263 : END DO
264 :
265 724 : CALL neighbor_list_iterator_release(nl_iterator)
266 :
267 : END IF
268 :
269 724 : CALL timestop(handle)
270 :
271 1448 : END SUBROUTINE calculate_lri_ks_matrix
272 :
273 : !*****************************************************************************
274 : !> \brief update of RIGPW KS matrix
275 : !> \param lri_env ...
276 : !> \param lri_v_int integrals of potential * ri basis set
277 : !> \param h_matrix KS matrix, on entry containing the core hamiltonian
278 : !> \param s_matrix overlap matrix
279 : !> \param atomic_kind_set ...
280 : !> \param ispin ...
281 : !> \note including this in lri_environment_methods?
282 : ! **************************************************************************************************
283 0 : SUBROUTINE calculate_ri_ks_matrix(lri_env, lri_v_int, h_matrix, s_matrix, &
284 : atomic_kind_set, ispin)
285 :
286 : TYPE(lri_environment_type), POINTER :: lri_env
287 : TYPE(lri_kind_type), DIMENSION(:), POINTER :: lri_v_int
288 : TYPE(dbcsr_type), POINTER :: h_matrix, s_matrix
289 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
290 : INTEGER, INTENT(IN) :: ispin
291 :
292 : CHARACTER(*), PARAMETER :: routineN = 'calculate_ri_ks_matrix'
293 :
294 : INTEGER :: atom_a, handle, i1, i2, iatom, ikind, n, &
295 : natom, nbas
296 0 : INTEGER, ALLOCATABLE, DIMENSION(:) :: atom_of_kind, kind_of, nsize
297 : INTEGER, DIMENSION(:, :), POINTER :: bas_ptr
298 : REAL(KIND=dp) :: fscal, ftrm1n
299 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: fout, fvec
300 0 : REAL(KIND=dp), DIMENSION(:), POINTER :: v
301 0 : TYPE(o3c_vec_type), DIMENSION(:), POINTER :: o3c_vec
302 :
303 0 : CALL timeset(routineN, handle)
304 :
305 0 : bas_ptr => lri_env%ri_fit%bas_ptr
306 0 : natom = SIZE(bas_ptr, 2)
307 0 : nbas = bas_ptr(2, natom)
308 0 : ALLOCATE (fvec(nbas), fout(nbas))
309 0 : CALL get_atomic_kind_set(atomic_kind_set, atom_of_kind=atom_of_kind, kind_of=kind_of)
310 0 : DO iatom = 1, natom
311 0 : ikind = kind_of(iatom)
312 0 : atom_a = atom_of_kind(iatom)
313 0 : i1 = bas_ptr(1, iatom)
314 0 : i2 = bas_ptr(2, iatom)
315 0 : n = i2 - i1 + 1
316 0 : fvec(i1:i2) = lri_v_int(ikind)%v_int(atom_a, 1:n)
317 : END DO
318 : ! f(T) * R^(-1)*n
319 0 : ftrm1n = SUM(fvec(:)*lri_env%ri_fit%rm1n(:))
320 0 : lri_env%ri_fit%ftrm1n(ispin) = ftrm1n
321 0 : fscal = ftrm1n/lri_env%ri_fit%ntrm1n
322 : ! renormalize fvec -> fvec - fscal * n
323 0 : fvec(:) = fvec(:) - fscal*lri_env%ri_fit%nvec(:)
324 : ! solve Rx=f'
325 : CALL ri_metric_solver(mat=lri_env%ri_smat(1)%matrix, &
326 : vecr=fvec(:), &
327 : vecx=fout(:), &
328 : matp=lri_env%ri_sinv(1)%matrix, &
329 : solver=lri_env%ri_sinv_app, &
330 0 : ptr=bas_ptr)
331 0 : lri_env%ri_fit%fout(:, ispin) = fout(:)
332 :
333 : ! add overlap matrix contribution
334 0 : CALL dbcsr_add(h_matrix, s_matrix, 1.0_dp, fscal)
335 :
336 : ! create a o3c_vec from fout
337 0 : ALLOCATE (nsize(natom), o3c_vec(natom))
338 0 : DO iatom = 1, natom
339 0 : i1 = bas_ptr(1, iatom)
340 0 : i2 = bas_ptr(2, iatom)
341 0 : n = i2 - i1 + 1
342 0 : nsize(iatom) = n
343 : END DO
344 0 : CALL o3c_vec_create(o3c_vec, nsize)
345 0 : DEALLOCATE (nsize)
346 0 : DO iatom = 1, natom
347 0 : i1 = bas_ptr(1, iatom)
348 0 : i2 = bas_ptr(2, iatom)
349 0 : n = i2 - i1 + 1
350 0 : CALL get_o3c_vec(o3c_vec, iatom, v)
351 0 : v(1:n) = fout(i1:i2)
352 : END DO
353 : ! add <T.f'>
354 0 : CALL contract3_o3c(lri_env%o3c, o3c_vec, h_matrix)
355 : !
356 0 : CALL o3c_vec_release(o3c_vec)
357 0 : DEALLOCATE (o3c_vec, fvec, fout)
358 :
359 0 : CALL timestop(handle)
360 :
361 0 : END SUBROUTINE calculate_ri_ks_matrix
362 :
363 : END MODULE lri_ks_methods
|
