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 Routines for the construction of the coefficients
10 : !> for the expansion of the atomic
11 : !> densities rho1_hard and rho1_soft in terms of primitive spherical gaussians.
12 : !> \par History
13 : !> 05-2004 created
14 : !> \author MI
15 : ! **************************************************************************************************
16 : MODULE qs_oce_methods
17 :
18 : USE ai_overlap, ONLY: overlap
19 : USE ao_util, ONLY: exp_radius
20 : USE basis_set_types, ONLY: get_gto_basis_set,&
21 : gto_basis_set_type
22 : USE block_p_types, ONLY: block_p_type
23 : USE kinds, ONLY: dp
24 : USE orbital_pointers, ONLY: init_orbital_pointers,&
25 : nco,&
26 : ncoset,&
27 : nso
28 : USE orbital_transformation_matrices, ONLY: orbtramat
29 : USE particle_types, ONLY: particle_type
30 : USE paw_basis_types, ONLY: get_paw_basis_info
31 : USE paw_proj_set_types, ONLY: get_paw_proj_set,&
32 : paw_proj_set_type
33 : USE qs_kind_types, ONLY: get_qs_kind,&
34 : get_qs_kind_set,&
35 : qs_kind_type
36 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
37 : USE sap_kind_types, ONLY: clist_type,&
38 : sap_int_type,&
39 : sap_sort
40 : #include "./base/base_uses.f90"
41 :
42 : IMPLICIT NONE
43 :
44 : PRIVATE
45 :
46 : ! Global parameters (only in this module)
47 :
48 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_oce_methods'
49 :
50 : ! Public subroutines
51 :
52 : PUBLIC :: build_oce_matrices, proj_blk, prj_scatter, prj_gather
53 :
54 : CONTAINS
55 :
56 : ! **************************************************************************************************
57 : !> \brief ...
58 : !> \param oces ...
59 : !> \param atom_ka ...
60 : !> \param atom_kb ...
61 : !> \param rab ...
62 : !> \param nder ...
63 : !> \param sgf_list ...
64 : !> \param nsgf_cnt ...
65 : !> \param sgf_soft_only ...
66 : !> \param eps_fit ...
67 : ! **************************************************************************************************
68 161450 : SUBROUTINE build_oce_block(oces, atom_ka, atom_kb, rab, nder, sgf_list, nsgf_cnt, sgf_soft_only, &
69 : eps_fit)
70 :
71 : TYPE(block_p_type), DIMENSION(:), POINTER :: oces
72 : TYPE(qs_kind_type), POINTER :: atom_ka, atom_kb
73 : REAL(KIND=dp), DIMENSION(3) :: rab
74 : INTEGER, INTENT(IN) :: nder
75 : INTEGER, DIMENSION(:), INTENT(OUT) :: sgf_list
76 : INTEGER, INTENT(OUT) :: nsgf_cnt
77 : LOGICAL, INTENT(OUT) :: sgf_soft_only
78 : REAL(KIND=dp), INTENT(IN) :: eps_fit
79 :
80 : INTEGER :: first_col, ic, ider, ig1, igau, ip, ipgf, is, isgfb, isgfb_cnt, isp, jc, jset, &
81 : lds, lm, lpoint, lprj, lsgfb, lsgfb_cnt, lshell, m, m1, maxcob, maxder, maxlb, maxlprj, &
82 : maxnprja, maxsoa, msab, n, ncob, np_car, np_sph, nsatbas, nseta, nsetb, nsoatot, &
83 : ntotsgfb, sgf_hard_only
84 161450 : INTEGER, DIMENSION(:), POINTER :: fp_cara, fp_spha, lb_max, lb_min, npgfb, &
85 161450 : nprjla, nsgfb
86 161450 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfb
87 : LOGICAL :: calculate_forces, paw_atom_a, paw_atom_b
88 : REAL(KIND=dp) :: dab, hard_radius_a, hard_radius_b, &
89 : radius, rcprja
90 161450 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ovs, spa_sb
91 161450 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: s
92 161450 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_b, zisomina, zisominb
93 161450 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: csprj, rpgfb, rzetprja, spa_tmp, sphi_b, &
94 161450 : zetb, zetprja
95 : TYPE(gto_basis_set_type), POINTER :: basis_1c_a, orb_basis_b
96 : TYPE(paw_proj_set_type), POINTER :: paw_proj_a, paw_proj_b
97 :
98 161450 : NULLIFY (basis_1c_a, paw_proj_a)
99 : CALL get_qs_kind(qs_kind=atom_ka, basis_set=basis_1c_a, basis_type="GAPW_1C", &
100 : paw_proj_set=paw_proj_a, paw_atom=paw_atom_a, &
101 161450 : hard_radius=hard_radius_a)
102 :
103 161450 : NULLIFY (orb_basis_b, paw_proj_b)
104 : CALL get_qs_kind(qs_kind=atom_kb, basis_set=orb_basis_b, &
105 : paw_proj_set=paw_proj_b, paw_atom=paw_atom_b, &
106 161450 : hard_radius=hard_radius_b)
107 :
108 161450 : IF (.NOT. paw_atom_a) RETURN
109 :
110 161450 : NULLIFY (nprjla, fp_cara, fp_spha, rzetprja, zetprja)
111 : CALL get_paw_proj_set(paw_proj_set=paw_proj_a, csprj=csprj, maxl=maxlprj, &
112 : nprj=nprjla, ncgauprj=np_car, nsgauprj=np_sph, nsatbas=nsatbas, rcprj=rcprja, &
113 : first_prj=fp_cara, first_prjs=fp_spha, &
114 161450 : zisomin=zisomina, rzetprj=rzetprja, zetprj=zetprja)
115 :
116 161450 : NULLIFY (first_sgfb, lb_max, lb_min, npgfb, nsgfb, rpgfb, sphi_b, set_radius_b, zetb, zisominb)
117 : CALL get_gto_basis_set(gto_basis_set=orb_basis_b, nset=nsetb, nsgf=ntotsgfb, &
118 : set_radius=set_radius_b, lmax=lb_max, lmin=lb_min, &
119 : npgf=npgfb, nsgf_set=nsgfb, pgf_radius=rpgfb, &
120 : sphi=sphi_b, zet=zetb, first_sgf=first_sgfb, &
121 161450 : maxco=maxcob, maxl=maxlb)
122 :
123 161450 : CALL get_gto_basis_set(gto_basis_set=basis_1c_a, nset=nseta, maxso=maxsoa)
124 :
125 : ! Add the block ab
126 645800 : dab = SQRT(SUM(rab*rab))
127 :
128 161450 : maxder = ncoset(nder)
129 161450 : nsoatot = maxsoa*nseta
130 161450 : maxnprja = SIZE(zetprja, 1)
131 :
132 : calculate_forces = .FALSE.
133 : IF (nder > 0) THEN
134 161450 : calculate_forces = .TRUE.
135 : END IF
136 :
137 161450 : lm = MAX(maxlb, maxlprj)
138 161450 : lds = ncoset(lm + nder + 1)
139 161450 : msab = MAX(maxnprja*ncoset(maxlprj), maxcob)
140 :
141 807250 : ALLOCATE (s(lds, lds, ncoset(nder + 1)))
142 645800 : ALLOCATE (spa_sb(np_car, ntotsgfb))
143 645800 : ALLOCATE (spa_tmp(msab, msab*maxder))
144 645800 : ALLOCATE (ovs(np_sph, maxcob*nsetb*maxder))
145 :
146 161450 : m1 = 0
147 161450 : nsgf_cnt = 0
148 161450 : isgfb_cnt = 1
149 161450 : sgf_hard_only = 0
150 514474 : DO jset = 1, nsetb
151 : !
152 : ! Set the contribution list
153 353024 : IF (hard_radius_a + set_radius_b(jset) >= dab) THEN
154 192900 : isgfb = first_sgfb(1, jset)
155 192900 : lsgfb = isgfb - 1 + nsgfb(jset)
156 1066996 : DO jc = isgfb, lsgfb
157 874096 : nsgf_cnt = nsgf_cnt + 1
158 1066996 : sgf_list(nsgf_cnt) = jc
159 : END DO
160 :
161 : ! check if this function is hard
162 1086466 : radius = exp_radius(lb_max(jset), MAXVAL(zetb(1:npgfb(jset), jset)), eps_fit, 1.0_dp)
163 192900 : IF (radius <= hard_radius_b) sgf_hard_only = sgf_hard_only + 1
164 :
165 : ! Integral between proj of iatom and primitives of jatom
166 : ! Calculate the primitives overlap
167 1676788512 : spa_tmp = 0.0_dp
168 579751842 : ovs = 0.0_dp
169 861188636 : s = 0.0_dp
170 192900 : ncob = npgfb(jset)*ncoset(lb_max(jset))
171 192900 : isgfb = first_sgfb(1, jset)
172 192900 : lsgfb = isgfb - 1 + nsgfb(jset)
173 :
174 192900 : lsgfb_cnt = isgfb_cnt - 1 + nsgfb(jset)
175 :
176 672568 : DO lprj = 0, maxlprj
177 : CALL overlap(lprj, lprj, nprjla(lprj), &
178 : rzetprja(:, lprj), zetprja(:, lprj), &
179 : lb_max(jset), lb_min(jset), npgfb(jset), &
180 : rpgfb(:, jset), zetb(:, jset), &
181 : -rab, dab, spa_tmp, &
182 1918672 : nder, .TRUE., s, lds)
183 1728038 : DO ider = 1, maxder
184 1055470 : is = (ider - 1)*SIZE(spa_tmp, 1)
185 1055470 : isp = (ider - 1)*maxcob*nsetb
186 7567580 : DO ipgf = 1, nprjla(lprj)
187 6032442 : lpoint = ncoset(lprj - 1) + 1 + (ipgf - 1)*ncoset(lprj)
188 6032442 : m = fp_spha(lprj) + (ipgf - 1)*nso(lprj)
189 34555106 : DO ip = 1, npgfb(jset)
190 27467194 : ic = (ip - 1)*ncoset(lb_max(jset))
191 27467194 : igau = isp + ic + m1 + ncoset(lb_min(jset) - 1) + 1
192 27467194 : ig1 = is + ic + ncoset(lb_min(jset) - 1) + 1
193 27467194 : n = ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1)
194 : ovs(m:m + nso(lprj) - 1, igau:igau + n - 1) = &
195 27467194 : MATMUL(orbtramat(lprj)%slm(1:nso(lprj), 1:nco(lprj)), &
196 1627042886 : spa_tmp(lpoint:lpoint + nco(lprj) - 1, ig1:ig1 + n - 1))
197 : END DO
198 : END DO
199 : END DO
200 : END DO
201 :
202 192900 : IF (paw_atom_b) THEN
203 181696 : CALL get_paw_proj_set(paw_proj_set=paw_proj_b, zisomin=zisominb)
204 853486 : DO ipgf = 1, npgfb(jset)
205 2040778 : DO lshell = lb_min(jset), lb_max(jset)
206 1859082 : IF (zetb(ipgf, jset) >= zisominb(lshell)) THEN
207 21278 : n = ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1)
208 21278 : igau = n*(ipgf - 1) + ncoset(lshell - 1)
209 64222 : DO ider = 1, maxder
210 42944 : is = maxcob*(ider - 1)
211 42944 : isp = (ider - 1)*maxcob*nsetb
212 1688338 : ovs(:, igau + 1 + isp + m1:igau + nco(lshell) + isp + m1) = 0.0_dp
213 : END DO
214 : END IF
215 : END DO
216 : END DO
217 : END IF
218 :
219 : ! Contraction step (integrals and derivatives)
220 615588 : DO ider = 1, maxder
221 422688 : first_col = (ider - 1)*maxcob*nsetb + 1 + m1
222 : ! CALL dgemm("N", "N", np_sph, nsgfb(jset), ncob, &
223 : ! 1.0_dp, ovs(1, first_col), SIZE(ovs, 1), &
224 : ! sphi_b(1, isgfb), SIZE(sphi_b, 1), &
225 : ! 0.0_dp, spa_sb(1, isgfb), SIZE(spa_sb, 1))
226 : spa_sb(1:np_sph, isgfb:isgfb + nsgfb(jset) - 1) = &
227 435978 : MATMUL(ovs(1:np_sph, first_col:first_col + ncob - 1), &
228 1244546842 : sphi_b(1:ncob, isgfb:isgfb + nsgfb(jset) - 1))
229 :
230 : ! CALL dgemm("T", "N", nsatbas, nsgfb(jset), np_sph, &
231 : ! 1.0_dp, csprj(1, 1), SIZE(csprj, 1), &
232 : ! spa_sb(1, isgfb), SIZE(spa_sb, 1), &
233 : ! 1.0_dp, oces(ider)%block(1, isgfb_cnt), SIZE(oces(ider)%block, 1))
234 : oces(ider)%block(1:nsatbas, isgfb_cnt:isgfb_cnt + nsgfb(jset) - 1) = &
235 : oces(ider)%block(1:nsatbas, isgfb_cnt:isgfb_cnt + nsgfb(jset) - 1) + &
236 1567912 : MATMUL(TRANSPOSE(csprj(1:np_sph, 1:nsatbas)), &
237 1759212240 : spa_sb(1:np_sph, isgfb:isgfb + nsgfb(jset) - 1))
238 : END DO
239 192900 : isgfb_cnt = isgfb_cnt + nsgfb(jset)
240 : END IF ! radius
241 514474 : m1 = m1 + maxcob
242 : END DO !jset
243 :
244 : ! Check if the screened functions are all soft
245 161450 : sgf_soft_only = .FALSE.
246 161450 : IF (sgf_hard_only == 0) sgf_soft_only = .TRUE.
247 :
248 161450 : DEALLOCATE (s, spa_sb, spa_tmp, ovs)
249 :
250 322900 : END SUBROUTINE build_oce_block
251 :
252 : ! **************************************************************************************************
253 : !> \brief ...
254 : !> \param oceh ...
255 : !> \param oces ...
256 : !> \param atom_ka ...
257 : !> \param sgf_list ...
258 : !> \param nsgf_cnt ...
259 : ! **************************************************************************************************
260 10252 : SUBROUTINE build_oce_block_local(oceh, oces, atom_ka, sgf_list, nsgf_cnt)
261 :
262 : TYPE(block_p_type), DIMENSION(:), POINTER :: oceh, oces
263 : TYPE(qs_kind_type), POINTER :: atom_ka
264 : INTEGER, DIMENSION(:), INTENT(OUT) :: sgf_list
265 : INTEGER, INTENT(OUT) :: nsgf_cnt
266 :
267 : INTEGER :: i, iset, isgfa, j, jc, lsgfa, maxlprj, &
268 : maxso1a, n, nsatbas, nset1a, nseta, &
269 : nsgfa
270 10252 : INTEGER, DIMENSION(:), POINTER :: n2oindex, nsgf_seta
271 10252 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
272 : LOGICAL :: paw_atom_a
273 10252 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: prjloc_h, prjloc_s
274 10252 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: local_oce_h, local_oce_s
275 : TYPE(gto_basis_set_type), POINTER :: basis_1c_a, orb_basis_a
276 : TYPE(paw_proj_set_type), POINTER :: paw_proj_a
277 :
278 10252 : NULLIFY (orb_basis_a, basis_1c_a, paw_proj_a)
279 : CALL get_qs_kind(qs_kind=atom_ka, basis_set=orb_basis_a, &
280 10252 : paw_proj_set=paw_proj_a, paw_atom=paw_atom_a)
281 :
282 10252 : IF (.NOT. paw_atom_a) RETURN
283 :
284 10252 : CALL get_paw_proj_set(paw_proj_set=paw_proj_a, maxl=maxlprj)
285 10252 : CALL get_qs_kind(qs_kind=atom_ka, basis_set=basis_1c_a, basis_type="GAPW_1C")
286 10252 : CALL get_gto_basis_set(gto_basis_set=basis_1c_a, nset=nset1a, maxso=maxso1a)
287 10252 : NULLIFY (n2oindex)
288 10252 : CALL get_paw_basis_info(basis_1c_a, n2oindex=n2oindex, nsatbas=nsatbas)
289 :
290 : CALL get_gto_basis_set(gto_basis_set=orb_basis_a, first_sgf=first_sgfa, &
291 10252 : nsgf=nsgfa, nsgf_set=nsgf_seta, nset=nseta)
292 :
293 10252 : NULLIFY (local_oce_h, local_oce_s)
294 : CALL get_paw_proj_set(paw_proj_set=paw_proj_a, &
295 : local_oce_sphi_h=local_oce_h, &
296 10252 : local_oce_sphi_s=local_oce_s)
297 :
298 61512 : ALLOCATE (prjloc_h(nset1a*maxso1a, nsgfa), prjloc_s(nset1a*maxso1a, nsgfa))
299 5301108 : prjloc_h = 0._dp
300 5301108 : prjloc_s = 0._dp
301 :
302 10252 : nsgf_cnt = 0
303 35538 : DO iset = 1, nseta
304 25286 : isgfa = first_sgfa(1, iset)
305 25286 : lsgfa = isgfa - 1 + nsgf_seta(iset)
306 111280 : DO jc = isgfa, lsgfa
307 85994 : nsgf_cnt = nsgf_cnt + 1
308 111280 : sgf_list(nsgf_cnt) = jc
309 : END DO
310 : ! this asumes that the first sets are the same for basis_1c/orb_basis!
311 25286 : n = maxso1a*(iset - 1)
312 1589144 : prjloc_h(n + 1:n + maxso1a, isgfa:lsgfa) = local_oce_h(1:maxso1a, isgfa:lsgfa)
313 1599396 : prjloc_s(n + 1:n + maxso1a, isgfa:lsgfa) = local_oce_s(1:maxso1a, isgfa:lsgfa)
314 : END DO
315 :
316 96246 : DO i = 1, nsgfa
317 2301836 : DO j = 1, nsatbas
318 2205590 : jc = n2oindex(j)
319 2205590 : oceh(1)%block(j, i) = prjloc_h(jc, i)
320 2291584 : oces(1)%block(j, i) = prjloc_s(jc, i)
321 : END DO
322 : END DO
323 :
324 10252 : DEALLOCATE (prjloc_h, prjloc_s)
325 10252 : DEALLOCATE (n2oindex)
326 :
327 30756 : END SUBROUTINE build_oce_block_local
328 :
329 : ! **************************************************************************************************
330 : !> \brief ...
331 : !> \param oceh ...
332 : !> \param oces ...
333 : !> \param atom_ka ...
334 : !> \param sgf_list ...
335 : !> \param nsgf_cnt ...
336 : !> \param eps_fit ...
337 : ! **************************************************************************************************
338 0 : SUBROUTINE build_oce_block_1c(oceh, oces, atom_ka, sgf_list, nsgf_cnt, eps_fit)
339 :
340 : TYPE(block_p_type), DIMENSION(:), POINTER :: oceh, oces
341 : TYPE(qs_kind_type), POINTER :: atom_ka
342 : INTEGER, DIMENSION(:), INTENT(OUT) :: sgf_list
343 : INTEGER, INTENT(OUT) :: nsgf_cnt
344 : REAL(KIND=dp), INTENT(IN) :: eps_fit
345 :
346 : INTEGER :: first_col, ic, ig1, igau, ip, ipgf, isgfb, isgfb_cnt, jc, jset, lds, lm, lpoint, &
347 : lprj, lsgfb, lsgfb_cnt, lshell, m, m1, maxcob, maxlb, maxlprj, maxnprja, maxsoa, msab, n, &
348 : ncob, np_car, np_sph, nsatbas, nseta, nsetb, nsoatot, ntotsgfb
349 0 : INTEGER, DIMENSION(:), POINTER :: fp_cara, fp_spha, lb_max, lb_min, npgfb, &
350 0 : nprjla, nsgfb
351 0 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfb
352 : LOGICAL :: paw_atom_a
353 : REAL(KIND=dp) :: radius, rc, rcprja
354 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ovh, ovs, spa_sb
355 0 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: s
356 0 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_b, zisominb
357 0 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: chprj, csprj, rpgfb, rzetprja, spa_tmp, &
358 0 : sphi_b, zetb, zetprja
359 : TYPE(gto_basis_set_type), POINTER :: basis_1c_a, orb_basis_b
360 : TYPE(paw_proj_set_type), POINTER :: paw_proj_a
361 :
362 0 : NULLIFY (orb_basis_b, basis_1c_a, paw_proj_a)
363 0 : CALL get_qs_kind(qs_kind=atom_ka, paw_atom=paw_atom_a)
364 :
365 0 : IF (.NOT. paw_atom_a) RETURN
366 :
367 0 : CALL get_qs_kind(qs_kind=atom_ka, paw_proj_set=paw_proj_a)
368 0 : NULLIFY (nprjla, fp_cara, fp_spha, rzetprja, zetprja)
369 : CALL get_paw_proj_set(paw_proj_set=paw_proj_a, csprj=csprj, chprj=chprj, maxl=maxlprj, &
370 : nprj=nprjla, ncgauprj=np_car, nsgauprj=np_sph, nsatbas=nsatbas, rcprj=rcprja, &
371 : first_prj=fp_cara, first_prjs=fp_spha, &
372 0 : rzetprj=rzetprja, zetprj=zetprja)
373 :
374 0 : CALL get_qs_kind(qs_kind=atom_ka, basis_set=orb_basis_b)
375 0 : NULLIFY (first_sgfb, lb_max, lb_min, npgfb, nsgfb, rpgfb, sphi_b, set_radius_b, zetb, zisominb)
376 : CALL get_gto_basis_set(gto_basis_set=orb_basis_b, nset=nsetb, nsgf=ntotsgfb, &
377 : set_radius=set_radius_b, lmax=lb_max, lmin=lb_min, &
378 : npgf=npgfb, nsgf_set=nsgfb, pgf_radius=rpgfb, &
379 : sphi=sphi_b, zet=zetb, first_sgf=first_sgfb, &
380 0 : maxco=maxcob, maxl=maxlb)
381 :
382 0 : CALL get_qs_kind(qs_kind=atom_ka, hard_radius=rc)
383 :
384 0 : CALL get_qs_kind(qs_kind=atom_ka, basis_set=basis_1c_a, basis_type="GAPW_1C")
385 0 : CALL get_gto_basis_set(gto_basis_set=basis_1c_a, nset=nseta, maxso=maxsoa)
386 :
387 : ! Add the block ab
388 0 : nsoatot = maxsoa*nseta
389 0 : maxnprja = SIZE(zetprja, 1)
390 :
391 0 : lm = MAX(maxlb, maxlprj)
392 0 : lds = ncoset(lm + 1)
393 0 : msab = MAX(maxnprja*ncoset(maxlprj), maxcob)
394 :
395 0 : ALLOCATE (s(lds, lds, 1))
396 0 : ALLOCATE (spa_sb(np_car, ntotsgfb))
397 0 : ALLOCATE (spa_tmp(msab, msab))
398 0 : ALLOCATE (ovs(np_sph, maxcob*nsetb))
399 0 : ALLOCATE (ovh(np_sph, maxcob*nsetb))
400 :
401 0 : m1 = 0
402 0 : nsgf_cnt = 0
403 0 : isgfb_cnt = 1
404 0 : DO jset = 1, nsetb
405 : !
406 : ! Set the contribution list
407 0 : isgfb = first_sgfb(1, jset)
408 0 : lsgfb = isgfb - 1 + nsgfb(jset)
409 0 : DO jc = isgfb, lsgfb
410 0 : nsgf_cnt = nsgf_cnt + 1
411 0 : sgf_list(nsgf_cnt) = jc
412 : END DO
413 :
414 : ! Integral between proj of iatom and primitives of iatom
415 : ! Calculate the primitives overlap
416 0 : spa_tmp = 0.0_dp
417 0 : ovs = 0.0_dp
418 0 : ovh = 0.0_dp
419 0 : s = 0.0_dp
420 0 : ncob = npgfb(jset)*ncoset(lb_max(jset))
421 0 : isgfb = first_sgfb(1, jset)
422 0 : lsgfb = isgfb - 1 + nsgfb(jset)
423 :
424 0 : lsgfb_cnt = isgfb_cnt - 1 + nsgfb(jset)
425 :
426 0 : DO lprj = 0, maxlprj
427 : CALL overlap(lprj, lprj, nprjla(lprj), &
428 : rzetprja(:, lprj), zetprja(:, lprj), &
429 : lb_max(jset), lb_min(jset), npgfb(jset), &
430 : rpgfb(:, jset), zetb(:, jset), &
431 : -[0._dp, 0._dp, 0._dp], 0.0_dp, spa_tmp, &
432 0 : 0, .TRUE., s, lds)
433 0 : DO ipgf = 1, nprjla(lprj)
434 0 : lpoint = ncoset(lprj - 1) + 1 + (ipgf - 1)*ncoset(lprj)
435 0 : m = fp_spha(lprj) + (ipgf - 1)*nso(lprj)
436 0 : DO ip = 1, npgfb(jset)
437 0 : ic = (ip - 1)*ncoset(lb_max(jset))
438 0 : igau = ic + m1 + ncoset(lb_min(jset) - 1) + 1
439 0 : ig1 = ic + ncoset(lb_min(jset) - 1) + 1
440 0 : n = ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1)
441 : ovs(m:m + nso(lprj) - 1, igau:igau + n - 1) = &
442 0 : MATMUL(orbtramat(lprj)%slm(1:nso(lprj), 1:nco(lprj)), &
443 0 : spa_tmp(lpoint:lpoint + nco(lprj) - 1, ig1:ig1 + n - 1))
444 : END DO
445 : END DO
446 : END DO
447 :
448 0 : ovh(:, :) = ovs(:, :)
449 :
450 0 : CALL get_paw_proj_set(paw_proj_set=paw_proj_a, zisomin=zisominb)
451 0 : DO ipgf = 1, npgfb(jset)
452 0 : DO lshell = lb_min(jset), lb_max(jset)
453 0 : radius = exp_radius(lshell, zetb(ipgf, jset), eps_fit, 1.0_dp)
454 0 : IF (radius < rc) THEN
455 0 : n = ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1)
456 0 : igau = n*(ipgf - 1) + ncoset(lshell - 1)
457 0 : ovs(:, igau + 1 + m1:igau + nco(lshell) + m1) = 0.0_dp
458 : END IF
459 : END DO
460 : END DO
461 :
462 : ! Contraction step (integrals and derivatives)
463 0 : first_col = 1 + m1
464 : spa_sb(1:np_sph, isgfb:isgfb + nsgfb(jset) - 1) = &
465 0 : MATMUL(ovs(1:np_sph, first_col:first_col + ncob - 1), &
466 0 : sphi_b(1:ncob, isgfb:isgfb + nsgfb(jset) - 1))
467 :
468 : oces(1)%block(1:nsatbas, isgfb_cnt:isgfb_cnt + nsgfb(jset) - 1) = &
469 : oces(1)%block(1:nsatbas, isgfb_cnt:isgfb_cnt + nsgfb(jset) - 1) + &
470 0 : MATMUL(TRANSPOSE(csprj(1:np_sph, 1:nsatbas)), &
471 0 : spa_sb(1:np_sph, isgfb:isgfb + nsgfb(jset) - 1))
472 :
473 : spa_sb(1:np_sph, isgfb:isgfb + nsgfb(jset) - 1) = &
474 0 : MATMUL(ovh(1:np_sph, first_col:first_col + ncob - 1), &
475 0 : sphi_b(1:ncob, isgfb:isgfb + nsgfb(jset) - 1))
476 :
477 : oceh(1)%block(1:nsatbas, isgfb_cnt:isgfb_cnt + nsgfb(jset) - 1) = &
478 : oceh(1)%block(1:nsatbas, isgfb_cnt:isgfb_cnt + nsgfb(jset) - 1) + &
479 0 : MATMUL(TRANSPOSE(chprj(1:np_sph, 1:nsatbas)), &
480 0 : spa_sb(1:np_sph, isgfb:isgfb + nsgfb(jset) - 1))
481 :
482 0 : isgfb_cnt = isgfb_cnt + nsgfb(jset)
483 0 : m1 = m1 + maxcob
484 : END DO !jset
485 :
486 0 : DEALLOCATE (s, spa_sb, spa_tmp, ovs)
487 :
488 0 : END SUBROUTINE build_oce_block_1c
489 :
490 : ! **************************************************************************************************
491 : !> \brief Set up the sparse matrix for the coefficients of one center expansions
492 : !> This routine uses the same logic as the nonlocal pseudopotential
493 : !> \param intac TYPE that holds the integrals (a=basis; c=projector)
494 : !> \param calculate_forces ...
495 : !> \param nder ...
496 : !> \param qs_kind_set ...
497 : !> \param particle_set ...
498 : !> \param sap_oce ...
499 : !> \param eps_fit ...
500 : !> \par History
501 : !> 02.2009 created
502 : !> \author jgh
503 : ! **************************************************************************************************
504 3578 : SUBROUTINE build_oce_matrices(intac, calculate_forces, nder, &
505 : qs_kind_set, particle_set, sap_oce, eps_fit)
506 :
507 : TYPE(sap_int_type), DIMENSION(:), POINTER :: intac
508 : LOGICAL, INTENT(IN) :: calculate_forces
509 : INTEGER :: nder
510 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
511 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
512 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
513 : POINTER :: sap_oce
514 : REAL(KIND=dp), INTENT(IN) :: eps_fit
515 :
516 : CHARACTER(LEN=*), PARAMETER :: routineN = 'build_oce_matrices'
517 :
518 : INTEGER :: atom_a, atom_b, handle, i, iac, ikind, ilist, jkind, jneighbor, ldai, ldsab, &
519 : maxco, maxder, maxl, maxlgto, maxlprj, maxprj, maxsgf, maxsoa, maxsob, mlprj, natom, &
520 : ncoa_sum, nkind, nlist, nneighbor, nsatbas, nseta, nsetb, nsgf_cnt, nsgfa, nsobtot, slot
521 3578 : INTEGER, ALLOCATABLE, DIMENSION(:) :: sgf_list
522 : INTEGER, DIMENSION(3) :: cell_b
523 3578 : INTEGER, DIMENSION(:), POINTER :: fp_car, fp_sph, la_max, la_min, npgfa, &
524 3578 : nprjla, nsgf_seta
525 3578 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
526 : LOGICAL :: local, paw_atom_b, sgf_soft_only
527 3578 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: has_1c_basis, has_orb_basis, paw_kind
528 : REAL(KIND=dp) :: dab, rcprj
529 3578 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: sab, work
530 3578 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: ai_work
531 : REAL(KIND=dp), DIMENSION(3) :: rab
532 3578 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a
533 3578 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rzetprj, sphi_a, zeta, zetb
534 3578 : TYPE(block_p_type), DIMENSION(:), POINTER :: oceh, oces
535 : TYPE(clist_type), POINTER :: clist
536 : TYPE(gto_basis_set_type), POINTER :: orb_basis_paw, orb_basis_set
537 : TYPE(paw_proj_set_type), POINTER :: paw_proj_b
538 : TYPE(qs_kind_type), POINTER :: at_a, at_b, qs_kind
539 :
540 7156 : IF (calculate_forces) THEN
541 1262 : CALL timeset(routineN//"_forces", handle)
542 : ELSE
543 2316 : CALL timeset(routineN, handle)
544 : END IF
545 :
546 3578 : IF (ASSOCIATED(sap_oce)) THEN
547 :
548 3578 : nkind = SIZE(qs_kind_set)
549 3578 : natom = SIZE(particle_set)
550 :
551 3578 : maxder = ncoset(nder)
552 :
553 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
554 : maxco=maxco, &
555 : maxlgto=maxlgto, &
556 : maxlprj=maxlprj, &
557 : maxco_proj=maxprj, &
558 3578 : maxsgf=maxsgf)
559 :
560 3578 : maxl = MAX(maxlgto, maxlprj)
561 3578 : CALL init_orbital_pointers(maxl + nder + 1)
562 :
563 18050 : DO i = 1, nkind*nkind
564 14472 : NULLIFY (intac(i)%alist, intac(i)%asort, intac(i)%aindex)
565 18050 : intac(i)%nalist = 0
566 : END DO
567 :
568 17890 : ALLOCATE (has_orb_basis(nkind), paw_kind(nkind), has_1c_basis(nkind))
569 10498 : has_orb_basis(:) = .FALSE.
570 10498 : paw_kind(:) = .FALSE.
571 10498 : has_1c_basis(:) = .FALSE.
572 10498 : DO ikind = 1, nkind
573 6920 : NULLIFY (orb_basis_set)
574 6920 : CALL get_qs_kind(qs_kind=qs_kind_set(ikind), basis_set=orb_basis_set)
575 6920 : has_orb_basis(ikind) = ASSOCIATED(orb_basis_set)
576 6920 : NULLIFY (paw_proj_b)
577 6920 : CALL get_qs_kind(qs_kind=qs_kind_set(ikind), paw_proj_set=paw_proj_b, paw_atom=paw_kind(ikind))
578 6920 : NULLIFY (orb_basis_paw)
579 6920 : CALL get_qs_kind(qs_kind=qs_kind_set(ikind), basis_set=orb_basis_paw, basis_type="GAPW_1C")
580 10498 : has_1c_basis(ikind) = ASSOCIATED(orb_basis_paw)
581 : END DO
582 :
583 : ! Allocate memory list
584 175280 : DO slot = 1, sap_oce(1)%nl_size
585 171702 : ikind = sap_oce(1)%nlist_task(slot)%ikind
586 171702 : jkind = sap_oce(1)%nlist_task(slot)%jkind
587 171702 : atom_a = sap_oce(1)%nlist_task(slot)%iatom
588 171702 : nlist = sap_oce(1)%nlist_task(slot)%nlist
589 171702 : ilist = sap_oce(1)%nlist_task(slot)%ilist
590 171702 : nneighbor = sap_oce(1)%nlist_task(slot)%nnode
591 :
592 171702 : iac = ikind + nkind*(jkind - 1)
593 :
594 171702 : IF (.NOT. has_orb_basis(ikind)) CYCLE
595 171702 : IF (.NOT. paw_kind(jkind)) CYCLE
596 171702 : IF (.NOT. has_1c_basis(jkind)) CYCLE
597 171702 : IF (.NOT. ASSOCIATED(intac(iac)%alist)) THEN
598 13212 : intac(iac)%a_kind = ikind
599 13212 : intac(iac)%p_kind = jkind
600 13212 : intac(iac)%nalist = nlist
601 60380 : ALLOCATE (intac(iac)%alist(nlist))
602 33956 : DO i = 1, nlist
603 20744 : NULLIFY (intac(iac)%alist(i)%clist)
604 20744 : intac(iac)%alist(i)%aatom = 0
605 33956 : intac(iac)%alist(i)%nclist = 0
606 : END DO
607 : END IF
608 175280 : IF (.NOT. ASSOCIATED(intac(iac)%alist(ilist)%clist)) THEN
609 20744 : intac(iac)%alist(ilist)%aatom = atom_a
610 20744 : intac(iac)%alist(ilist)%nclist = nneighbor
611 379142 : ALLOCATE (intac(iac)%alist(ilist)%clist(nneighbor))
612 : END IF
613 : END DO
614 :
615 3578 : ldsab = MAX(maxco, ncoset(maxlprj), maxsgf, maxprj)
616 3578 : ldai = ncoset(maxl + nder + 1)
617 :
618 : !calculate the overlap integrals <a|p>
619 : !$OMP PARALLEL DEFAULT(NONE) &
620 : !$OMP SHARED (intac, ldsab, ldai, nder, nkind, maxder, ncoset, sap_oce, qs_kind_set, eps_fit, &
621 : !$OMP has_orb_basis, paw_kind, has_1c_basis, maxsgf) &
622 : !$OMP PRIVATE (sab, work, ai_work, oceh, oces, slot, ikind, jkind, atom_a, atom_b, ilist, jneighbor, rab, cell_b, &
623 : !$OMP iac, dab, qs_kind, orb_basis_set, first_sgfa, la_max, la_min, ncoa_sum, maxsoa, npgfa, nseta, &
624 : !$OMP nsgfa, nsgf_seta, rpgfa, set_radius_a, sphi_a, zeta, paw_proj_b, paw_atom_b, orb_basis_paw, &
625 : !$OMP maxsob, nsetb, mlprj, nprjla, nsatbas, rcprj, fp_car, fp_sph, rzetprj, zetb, nsobtot, clist, &
626 3578 : !$OMP sgf_list, at_a, at_b, local, i, sgf_soft_only, nsgf_cnt)
627 :
628 : ALLOCATE (sab(ldsab, ldsab*maxder), work(ldsab, ldsab*maxder))
629 : sab = 0.0_dp
630 : ALLOCATE (ai_work(ldai, ldai, ncoset(nder + 1)))
631 : ai_work = 0.0_dp
632 : ALLOCATE (oceh(maxder), oces(maxder))
633 : ALLOCATE (sgf_list(maxsgf))
634 :
635 : !$OMP DO SCHEDULE(GUIDED)
636 : DO slot = 1, sap_oce(1)%nl_size
637 : ikind = sap_oce(1)%nlist_task(slot)%ikind
638 : jkind = sap_oce(1)%nlist_task(slot)%jkind
639 : atom_a = sap_oce(1)%nlist_task(slot)%iatom
640 : atom_b = sap_oce(1)%nlist_task(slot)%jatom
641 : ilist = sap_oce(1)%nlist_task(slot)%ilist
642 : jneighbor = sap_oce(1)%nlist_task(slot)%inode
643 : rab(1:3) = sap_oce(1)%nlist_task(slot)%r(1:3)
644 : cell_b(1:3) = sap_oce(1)%nlist_task(slot)%cell(1:3)
645 :
646 : iac = ikind + nkind*(jkind - 1)
647 : dab = SQRT(SUM(rab*rab))
648 :
649 : IF (.NOT. has_orb_basis(ikind)) CYCLE
650 : IF (.NOT. paw_kind(jkind)) CYCLE
651 : IF (.NOT. has_1c_basis(jkind)) CYCLE
652 :
653 : qs_kind => qs_kind_set(ikind)
654 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis_set)
655 :
656 : IF (.NOT. ASSOCIATED(orb_basis_set)) CYCLE
657 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
658 : first_sgf=first_sgfa, &
659 : lmax=la_max, &
660 : lmin=la_min, &
661 : nco_sum=ncoa_sum, &
662 : maxso=maxsoa, &
663 : npgf=npgfa, &
664 : nset=nseta, &
665 : nsgf=nsgfa, &
666 : nsgf_set=nsgf_seta, &
667 : pgf_radius=rpgfa, &
668 : set_radius=set_radius_a, &
669 : sphi=sphi_a, &
670 : zet=zeta)
671 :
672 : qs_kind => qs_kind_set(jkind)
673 :
674 : NULLIFY (paw_proj_b)
675 : CALL get_qs_kind(qs_kind=qs_kind, paw_proj_set=paw_proj_b, paw_atom=paw_atom_b)
676 : IF (.NOT. paw_atom_b) CYCLE
677 :
678 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis_paw, basis_type="GAPW_1C")
679 : IF (.NOT. ASSOCIATED(orb_basis_paw)) CYCLE
680 : CALL get_gto_basis_set(gto_basis_set=orb_basis_paw, maxso=maxsob, nset=nsetb)
681 :
682 : CALL get_paw_proj_set(paw_proj_set=paw_proj_b, &
683 : maxl=mlprj, &
684 : nprj=nprjla, &
685 : nsatbas=nsatbas, &
686 : rcprj=rcprj, &
687 : first_prj=fp_car, &
688 : first_prjs=fp_sph, &
689 : rzetprj=rzetprj, &
690 : zetprj=zetb)
691 :
692 : nsobtot = nsatbas
693 :
694 : clist => intac(iac)%alist(ilist)%clist(jneighbor)
695 : clist%catom = atom_b
696 : clist%cell = cell_b
697 : clist%rac = rab
698 : clist%nsgf_cnt = 0
699 : clist%maxac = 0.0_dp
700 : clist%maxach = 0.0_dp
701 : NULLIFY (clist%acint, clist%achint, clist%sgf_list)
702 :
703 : at_a => qs_kind_set(jkind)
704 : at_b => qs_kind_set(ikind)
705 :
706 : local = (atom_a == atom_b .AND. ALL(cell_b == 0))
707 :
708 : IF (local) THEN
709 : DO i = 1, maxder
710 : ALLOCATE (oceh(i)%block(nsobtot, nsgfa), oces(i)%block(nsobtot, nsgfa))
711 : oceh(i)%block = 0._dp
712 : oces(i)%block = 0._dp
713 : END DO
714 : CALL build_oce_block_local(oceh, oces, at_a, sgf_list, nsgf_cnt)
715 : clist%nsgf_cnt = nsgf_cnt
716 : clist%sgf_soft_only = .FALSE.
717 : IF (nsgf_cnt > 0) THEN
718 : ALLOCATE (clist%acint(nsgf_cnt, nsobtot, maxder), clist%sgf_list(nsgf_cnt))
719 : clist%acint(:, :, :) = 0._dp
720 : clist%sgf_list(:) = HUGE(0)
721 : CPASSERT(nsgf_cnt == nsgfa)
722 : ! *** Special case: A=B
723 : ALLOCATE (clist%achint(nsgfa, nsobtot, maxder))
724 : clist%achint = 0._dp
725 : clist%acint(1:nsgfa, 1:nsobtot, 1) = TRANSPOSE(oces(1)%block(1:nsobtot, 1:nsgfa))
726 : clist%achint(1:nsgfa, 1:nsobtot, 1) = TRANSPOSE(oceh(1)%block(1:nsobtot, 1:nsgfa))
727 : clist%maxac = MAXVAL(ABS(clist%acint(:, :, 1)))
728 : clist%maxach = 0._dp
729 : clist%sgf_list(1:nsgf_cnt) = sgf_list(1:nsgf_cnt)
730 : END IF
731 : DO i = 1, maxder
732 : DEALLOCATE (oceh(i)%block, oces(i)%block)
733 : END DO
734 : ELSE
735 : DO i = 1, maxder
736 : ALLOCATE (oces(i)%block(nsobtot, nsgfa))
737 : oces(i)%block = 0._dp
738 : END DO
739 : CALL build_oce_block(oces, at_a, at_b, rab, nder, sgf_list, nsgf_cnt, sgf_soft_only, eps_fit)
740 : clist%nsgf_cnt = nsgf_cnt
741 : clist%sgf_soft_only = sgf_soft_only
742 : IF (nsgf_cnt > 0) THEN
743 : ALLOCATE (clist%acint(nsgf_cnt, nsobtot, maxder), clist%sgf_list(nsgf_cnt))
744 : clist%acint(:, :, :) = 0._dp
745 : clist%sgf_list(:) = HUGE(0)
746 : DO i = 1, maxder
747 : clist%acint(1:nsgf_cnt, 1:nsobtot, i) = TRANSPOSE(oces(i)%block(1:nsobtot, 1:nsgf_cnt))
748 : END DO
749 : clist%maxac = MAXVAL(ABS(clist%acint(:, :, 1)))
750 : clist%maxach = 0._dp
751 : clist%sgf_list(1:nsgf_cnt) = sgf_list(1:nsgf_cnt)
752 : END IF
753 : DO i = 1, maxder
754 : DEALLOCATE (oces(i)%block)
755 : END DO
756 : END IF
757 :
758 : END DO
759 :
760 : DEALLOCATE (sab, work, ai_work)
761 : DEALLOCATE (sgf_list)
762 : DEALLOCATE (oceh, oces)
763 : !$OMP END PARALLEL
764 :
765 : ! Setup sort index
766 3578 : CALL sap_sort(intac)
767 :
768 7156 : DEALLOCATE (has_orb_basis, paw_kind, has_1c_basis)
769 :
770 : END IF
771 :
772 3578 : CALL timestop(handle)
773 :
774 7156 : END SUBROUTINE build_oce_matrices
775 :
776 : ! **************************************************************************************************
777 : !> \brief Project a matrix block onto the local atomic functions.
778 : !>
779 : !> \param h_a ...
780 : !> \param s_a ...
781 : !> \param na ...
782 : !> \param h_b ...
783 : !> \param s_b ...
784 : !> \param nb ...
785 : !> \param blk ...
786 : !> \param ldb ...
787 : !> \param proj_h ...
788 : !> \param proj_s ...
789 : !> \param nso ...
790 : !> \param len1 ...
791 : !> \param len2 ...
792 : !> \param fac ...
793 : !> \param distab ...
794 : !> \param work1 ...
795 : !> \param work2 ...
796 : !> \par History
797 : !> 02.2009 created
798 : !> 09.2016 use automatic arrays [M Tucker]
799 : !> 05.2026 move work arrays to callers
800 : ! **************************************************************************************************
801 8063318 : SUBROUTINE proj_blk(h_a, s_a, na, h_b, s_b, nb, blk, ldb, proj_h, proj_s, nso, len1, len2, fac, distab, &
802 8063318 : work1, work2)
803 :
804 : INTEGER, INTENT(IN) :: na
805 : REAL(KIND=dp), INTENT(IN) :: s_a(na, *), h_a(na, *)
806 : INTEGER, INTENT(IN) :: nb
807 : REAL(KIND=dp), INTENT(IN) :: s_b(nb, *), h_b(nb, *)
808 : INTEGER, INTENT(IN) :: ldb
809 : REAL(KIND=dp), INTENT(IN) :: blk(ldb, *)
810 : INTEGER, INTENT(IN) :: nso
811 : REAL(KIND=dp), INTENT(INOUT) :: proj_s(nso, *), proj_h(nso, *)
812 : INTEGER, INTENT(IN) :: len1, len2
813 : REAL(KIND=dp), INTENT(IN) :: fac
814 : LOGICAL, INTENT(IN) :: distab
815 : REAL(KIND=dp), INTENT(INOUT) :: work1(len1), work2(len2)
816 :
817 8063318 : IF (na == 0 .OR. nb == 0 .OR. nso == 0) RETURN
818 :
819 : ! Handle special cases
820 8063318 : IF (na == 1 .AND. nb == 1) THEN
821 : ! hard
822 282312 : CALL dger(nso, nso, fac*blk(1, 1), h_a(1, 1), 1, h_b(1, 1), 1, proj_h(1, 1), nso)
823 : ! soft
824 282312 : CALL dger(nso, nso, fac*blk(1, 1), s_a(1, 1), 1, s_b(1, 1), 1, proj_s(1, 1), nso)
825 : ELSE
826 7781006 : IF (distab) THEN
827 : ! hard
828 7292033 : CALL dgemm('N', 'N', na, nso, nb, fac, blk(1, 1), ldb, h_b(1, 1), nb, 0.0_dp, work1(1), na)
829 7292033 : CALL dgemm('T', 'N', nso, nso, na, 1.0_dp, h_a(1, 1), na, work1(1), na, 0.0_dp, work2(1), nso)
830 7292033 : CALL daxpy(nso*nso, 1.0_dp, work2(1), 1, proj_h(1, 1), 1)
831 : ! soft
832 7292033 : CALL daxpy(nso*nso, 1.0_dp, work2(1), 1, proj_s(1, 1), 1)
833 : ELSE
834 : ! hard
835 488973 : CALL dgemm('N', 'N', na, nso, nb, fac, blk(1, 1), ldb, h_b(1, 1), nb, 0.0_dp, work1(1), na)
836 488973 : CALL dgemm('T', 'N', nso, nso, na, 1.0_dp, h_a(1, 1), na, work1(1), na, 1.0_dp, proj_h(1, 1), nso)
837 : ! soft
838 488973 : CALL dgemm('N', 'N', na, nso, nb, fac, blk(1, 1), ldb, s_b(1, 1), nb, 0.0_dp, work1(1), na)
839 488973 : CALL dgemm('T', 'N', nso, nso, na, 1.0_dp, s_a(1, 1), na, work1(1), na, 1.0_dp, proj_s(1, 1), nso)
840 : END IF
841 : END IF
842 :
843 : END SUBROUTINE proj_blk
844 :
845 : ! **************************************************************************************************
846 : !> \brief ...
847 : !> \param ain matrix in old indexing
848 : !> \param aout matrix in new compressed indexing
849 : !> \param atom ...
850 : ! **************************************************************************************************
851 111824 : SUBROUTINE prj_gather(ain, aout, atom)
852 :
853 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: ain
854 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: aout
855 : TYPE(qs_kind_type), INTENT(IN) :: atom
856 :
857 : INTEGER :: i, ip, j, jp, nbas
858 111824 : INTEGER, DIMENSION(:), POINTER :: n2oindex
859 : TYPE(gto_basis_set_type), POINTER :: basis_1c
860 :
861 111824 : NULLIFY (basis_1c)
862 111824 : CALL get_qs_kind(qs_kind=atom, basis_set=basis_1c, basis_type="GAPW_1C")
863 111824 : NULLIFY (n2oindex)
864 111824 : CALL get_paw_basis_info(basis_1c, n2oindex=n2oindex, nsatbas=nbas)
865 :
866 1953902 : DO i = 1, nbas
867 1842078 : ip = n2oindex(i)
868 53076208 : DO j = 1, nbas
869 51122306 : jp = n2oindex(j)
870 52964384 : aout(j, i) = ain(jp, ip)
871 : END DO
872 : END DO
873 :
874 111824 : DEALLOCATE (n2oindex)
875 :
876 111824 : END SUBROUTINE prj_gather
877 :
878 : ! **************************************************************************************************
879 : !> \brief ...
880 : !> \param ain matrix in new compressed indexing
881 : !> \param aout matrix in old indexing (addup)
882 : !> \param atom ...
883 : ! **************************************************************************************************
884 183064 : SUBROUTINE prj_scatter(ain, aout, atom)
885 :
886 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: ain
887 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: aout
888 : TYPE(qs_kind_type), INTENT(IN) :: atom
889 :
890 : INTEGER :: i, ip, j, jp, nbas
891 183064 : INTEGER, DIMENSION(:), POINTER :: n2oindex
892 : TYPE(gto_basis_set_type), POINTER :: basis_1c
893 :
894 183064 : NULLIFY (basis_1c)
895 183064 : CALL get_qs_kind(qs_kind=atom, basis_set=basis_1c, basis_type="GAPW_1C")
896 183064 : NULLIFY (n2oindex)
897 183064 : CALL get_paw_basis_info(basis_1c, n2oindex=n2oindex, nsatbas=nbas)
898 :
899 3294580 : DO i = 1, nbas
900 3111516 : ip = n2oindex(i)
901 93691008 : DO j = 1, nbas
902 90396428 : jp = n2oindex(j)
903 93507944 : aout(jp, ip) = aout(jp, ip) + ain(j, i)
904 : END DO
905 : END DO
906 :
907 183064 : DEALLOCATE (n2oindex)
908 :
909 183064 : END SUBROUTINE prj_scatter
910 :
911 422688 : END MODULE qs_oce_methods
|
