Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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 66076 : 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 66076 : INTEGER, DIMENSION(:), POINTER :: fp_cara, fp_spha, lb_max, lb_min, npgfb, &
85 66076 : nprjla, nsgfb
86 66076 : 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 66076 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: ovs, spa_sb
91 66076 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: s
92 66076 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_b, zisomina, zisominb
93 66076 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: csprj, rpgfb, rzetprja, spa_tmp, sphi_b, &
94 66076 : 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 66076 : 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 66076 : hard_radius=hard_radius_a)
102 :
103 66076 : 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 66076 : hard_radius=hard_radius_b)
107 :
108 66076 : IF (.NOT. paw_atom_a) RETURN
109 :
110 66076 : 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 66076 : zisomin=zisomina, rzetprj=rzetprja, zetprj=zetprja)
115 :
116 66076 : 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 66076 : maxco=maxcob, maxl=maxlb)
122 :
123 66076 : CALL get_gto_basis_set(gto_basis_set=basis_1c_a, nset=nseta, maxso=maxsoa)
124 :
125 : ! Add the block ab
126 264304 : dab = SQRT(SUM(rab*rab))
127 :
128 66076 : maxder = ncoset(nder)
129 66076 : nsoatot = maxsoa*nseta
130 66076 : maxnprja = SIZE(zetprja, 1)
131 :
132 : calculate_forces = .FALSE.
133 : IF (nder > 0) THEN
134 66076 : calculate_forces = .TRUE.
135 : END IF
136 :
137 66076 : lm = MAX(maxlb, maxlprj)
138 66076 : lds = ncoset(lm + nder + 1)
139 66076 : msab = MAX(maxnprja*ncoset(maxlprj), maxcob)
140 :
141 330380 : ALLOCATE (s(lds, lds, ncoset(nder + 1)))
142 264304 : ALLOCATE (spa_sb(np_car, ntotsgfb))
143 264304 : ALLOCATE (spa_tmp(msab, msab*maxder))
144 264304 : ALLOCATE (ovs(np_sph, maxcob*nsetb*maxder))
145 :
146 66076 : m1 = 0
147 66076 : nsgf_cnt = 0
148 66076 : isgfb_cnt = 1
149 66076 : sgf_hard_only = 0
150 232736 : DO jset = 1, nsetb
151 : !
152 : ! Set the contribution list
153 166660 : IF (hard_radius_a + set_radius_b(jset) >= dab) THEN
154 85498 : isgfb = first_sgfb(1, jset)
155 85498 : lsgfb = isgfb - 1 + nsgfb(jset)
156 460058 : DO jc = isgfb, lsgfb
157 374560 : nsgf_cnt = nsgf_cnt + 1
158 460058 : sgf_list(nsgf_cnt) = jc
159 : END DO
160 :
161 : ! check if this function is hard
162 403906 : radius = exp_radius(lb_max(jset), MAXVAL(zetb(1:npgfb(jset), jset)), eps_fit, 1.0_dp)
163 85498 : IF (radius .LE. 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 510803194 : spa_tmp = 0.0_dp
168 150547990 : ovs = 0.0_dp
169 338160686 : s = 0.0_dp
170 85498 : ncob = npgfb(jset)*ncoset(lb_max(jset))
171 85498 : isgfb = first_sgfb(1, jset)
172 85498 : lsgfb = isgfb - 1 + nsgfb(jset)
173 :
174 85498 : lsgfb_cnt = isgfb_cnt - 1 + nsgfb(jset)
175 :
176 295010 : 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 838048 : nder, .TRUE., s, lds)
183 708120 : DO ider = 1, maxder
184 413110 : is = (ider - 1)*SIZE(spa_tmp, 1)
185 413110 : isp = (ider - 1)*maxcob*nsetb
186 2243948 : DO ipgf = 1, nprjla(lprj)
187 1621326 : lpoint = ncoset(lprj - 1) + 1 + (ipgf - 1)*ncoset(lprj)
188 1621326 : m = fp_spha(lprj) + (ipgf - 1)*nso(lprj)
189 6556582 : DO ip = 1, npgfb(jset)
190 4522146 : ic = (ip - 1)*ncoset(lb_max(jset))
191 4522146 : igau = isp + ic + m1 + ncoset(lb_min(jset) - 1) + 1
192 4522146 : ig1 = is + ic + ncoset(lb_min(jset) - 1) + 1
193 4522146 : n = ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1)
194 : ovs(m:m + nso(lprj) - 1, igau:igau + n - 1) = &
195 4522146 : MATMUL(orbtramat(lprj)%slm(1:nso(lprj), 1:nco(lprj)), &
196 238348958 : 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 85498 : IF (paw_atom_b) THEN
203 75614 : CALL get_paw_proj_set(paw_proj_set=paw_proj_b, zisomin=zisominb)
204 281910 : DO ipgf = 1, npgfb(jset)
205 605292 : DO lshell = lb_min(jset), lb_max(jset)
206 529678 : IF (zetb(ipgf, jset) >= zisominb(lshell)) THEN
207 14206 : n = ncoset(lb_max(jset)) - ncoset(lb_min(jset) - 1)
208 14206 : igau = n*(ipgf - 1) + ncoset(lshell - 1)
209 41510 : DO ider = 1, maxder
210 27304 : is = maxcob*(ider - 1)
211 27304 : isp = (ider - 1)*maxcob*nsetb
212 751950 : 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 253994 : DO ider = 1, maxder
221 168496 : 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 173382 : MATMUL(ovs(1:np_sph, first_col:first_col + ncob - 1), &
228 207547452 : 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 669624 : MATMUL(TRANSPOSE(csprj(1:np_sph, 1:nsatbas)), &
237 452211742 : spa_sb(1:np_sph, isgfb:isgfb + nsgfb(jset) - 1))
238 : END DO
239 85498 : isgfb_cnt = isgfb_cnt + nsgfb(jset)
240 : END IF ! radius
241 232736 : m1 = m1 + maxcob
242 : END DO !jset
243 :
244 : ! Check if the screened functions are all soft
245 66076 : sgf_soft_only = .FALSE.
246 66076 : IF (sgf_hard_only .EQ. 0) sgf_soft_only = .TRUE.
247 :
248 66076 : DEALLOCATE (s, spa_sb, spa_tmp, ovs)
249 :
250 132152 : 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 6848 : 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 6848 : INTEGER, DIMENSION(:), POINTER :: n2oindex, nsgf_seta
271 6848 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
272 : LOGICAL :: paw_atom_a
273 6848 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: prjloc_h, prjloc_s
274 6848 : 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 6848 : 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 6848 : paw_proj_set=paw_proj_a, paw_atom=paw_atom_a)
281 :
282 6848 : IF (.NOT. paw_atom_a) RETURN
283 :
284 6848 : CALL get_paw_proj_set(paw_proj_set=paw_proj_a, maxl=maxlprj)
285 6848 : CALL get_qs_kind(qs_kind=atom_ka, basis_set=basis_1c_a, basis_type="GAPW_1C")
286 6848 : CALL get_gto_basis_set(gto_basis_set=basis_1c_a, nset=nset1a, maxso=maxso1a)
287 6848 : NULLIFY (n2oindex)
288 6848 : 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 6848 : nsgf=nsgfa, nsgf_set=nsgf_seta, nset=nseta)
292 :
293 6848 : 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 6848 : local_oce_sphi_s=local_oce_s)
297 :
298 41088 : ALLOCATE (prjloc_h(nset1a*maxso1a, nsgfa), prjloc_s(nset1a*maxso1a, nsgfa))
299 3318440 : prjloc_h = 0._dp
300 3318440 : prjloc_s = 0._dp
301 :
302 6848 : nsgf_cnt = 0
303 24294 : DO iset = 1, nseta
304 17446 : isgfa = first_sgfa(1, iset)
305 17446 : lsgfa = isgfa - 1 + nsgf_seta(iset)
306 74574 : DO jc = isgfa, lsgfa
307 57128 : nsgf_cnt = nsgf_cnt + 1
308 74574 : sgf_list(nsgf_cnt) = jc
309 : END DO
310 : ! this asumes that the first sets are the same for basis_1c/orb_basis!
311 17446 : n = maxso1a*(iset - 1)
312 981972 : prjloc_h(n + 1:n + maxso1a, isgfa:lsgfa) = local_oce_h(1:maxso1a, isgfa:lsgfa)
313 988820 : prjloc_s(n + 1:n + maxso1a, isgfa:lsgfa) = local_oce_s(1:maxso1a, isgfa:lsgfa)
314 : END DO
315 :
316 63976 : DO i = 1, nsgfa
317 1489978 : DO j = 1, nsatbas
318 1426002 : jc = n2oindex(j)
319 1426002 : oceh(1)%block(j, i) = prjloc_h(jc, i)
320 1483130 : oces(1)%block(j, i) = prjloc_s(jc, i)
321 : END DO
322 : END DO
323 :
324 6848 : DEALLOCATE (prjloc_h, prjloc_s)
325 6848 : DEALLOCATE (n2oindex)
326 :
327 20544 : 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 2352 : 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 2352 : INTEGER, ALLOCATABLE, DIMENSION(:) :: sgf_list
522 : INTEGER, DIMENSION(3) :: cell_b
523 2352 : INTEGER, DIMENSION(:), POINTER :: fp_car, fp_sph, la_max, la_min, npgfa, &
524 2352 : nprjla, nsgf_seta
525 2352 : INTEGER, DIMENSION(:, :), POINTER :: first_sgfa
526 : LOGICAL :: local, paw_atom_b, sgf_soft_only
527 : REAL(KIND=dp) :: dab, rcprj
528 2352 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: sab, work
529 2352 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: ai_work
530 : REAL(KIND=dp), DIMENSION(3) :: rab
531 2352 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius_a
532 2352 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rpgfa, rzetprj, sphi_a, zeta, zetb
533 2352 : TYPE(block_p_type), DIMENSION(:), POINTER :: oceh, oces
534 : TYPE(clist_type), POINTER :: clist
535 : TYPE(gto_basis_set_type), POINTER :: orb_basis_paw, orb_basis_set
536 : TYPE(paw_proj_set_type), POINTER :: paw_proj_b
537 : TYPE(qs_kind_type), POINTER :: at_a, at_b, qs_kind
538 :
539 4704 : IF (calculate_forces) THEN
540 786 : CALL timeset(routineN//"_forces", handle)
541 : ELSE
542 1566 : CALL timeset(routineN, handle)
543 : END IF
544 :
545 2352 : IF (ASSOCIATED(sap_oce)) THEN
546 :
547 2352 : nkind = SIZE(qs_kind_set)
548 2352 : natom = SIZE(particle_set)
549 :
550 2352 : maxder = ncoset(nder)
551 :
552 : CALL get_qs_kind_set(qs_kind_set=qs_kind_set, &
553 : maxco=maxco, &
554 : maxlgto=maxlgto, &
555 : maxlprj=maxlprj, &
556 : maxco_proj=maxprj, &
557 2352 : maxsgf=maxsgf)
558 :
559 2352 : maxl = MAX(maxlgto, maxlprj)
560 2352 : CALL init_orbital_pointers(maxl + nder + 1)
561 :
562 12164 : DO i = 1, nkind*nkind
563 9812 : NULLIFY (intac(i)%alist, intac(i)%asort, intac(i)%aindex)
564 12164 : intac(i)%nalist = 0
565 : END DO
566 :
567 : ! Allocate memory list
568 75276 : DO slot = 1, sap_oce(1)%nl_size
569 72924 : ikind = sap_oce(1)%nlist_task(slot)%ikind
570 72924 : jkind = sap_oce(1)%nlist_task(slot)%jkind
571 72924 : atom_a = sap_oce(1)%nlist_task(slot)%iatom
572 72924 : nlist = sap_oce(1)%nlist_task(slot)%nlist
573 72924 : ilist = sap_oce(1)%nlist_task(slot)%ilist
574 72924 : nneighbor = sap_oce(1)%nlist_task(slot)%nnode
575 :
576 72924 : iac = ikind + nkind*(jkind - 1)
577 :
578 72924 : qs_kind => qs_kind_set(ikind)
579 72924 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis_set)
580 72924 : IF (.NOT. ASSOCIATED(orb_basis_set)) CYCLE
581 72924 : qs_kind => qs_kind_set(jkind)
582 72924 : NULLIFY (paw_proj_b)
583 72924 : CALL get_qs_kind(qs_kind=qs_kind, paw_proj_set=paw_proj_b, paw_atom=paw_atom_b)
584 72924 : IF (.NOT. paw_atom_b) CYCLE
585 72924 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis_paw, basis_type="GAPW_1C")
586 72924 : IF (.NOT. ASSOCIATED(orb_basis_paw)) CYCLE
587 72924 : IF (.NOT. ASSOCIATED(intac(iac)%alist)) THEN
588 8892 : intac(iac)%a_kind = ikind
589 8892 : intac(iac)%p_kind = jkind
590 8892 : intac(iac)%nalist = nlist
591 40884 : ALLOCATE (intac(iac)%alist(nlist))
592 23100 : DO i = 1, nlist
593 14208 : NULLIFY (intac(iac)%alist(i)%clist)
594 14208 : intac(iac)%alist(i)%aatom = 0
595 23100 : intac(iac)%alist(i)%nclist = 0
596 : END DO
597 : END IF
598 148200 : IF (.NOT. ASSOCIATED(intac(iac)%alist(ilist)%clist)) THEN
599 14208 : intac(iac)%alist(ilist)%aatom = atom_a
600 14208 : intac(iac)%alist(ilist)%nclist = nneighbor
601 215004 : ALLOCATE (intac(iac)%alist(ilist)%clist(nneighbor))
602 : END IF
603 : END DO
604 :
605 2352 : ldsab = MAX(maxco, ncoset(maxlprj), maxsgf, maxprj)
606 2352 : ldai = ncoset(maxl + nder + 1)
607 :
608 : !calculate the overlap integrals <a|p>
609 : !$OMP PARALLEL DEFAULT(NONE) &
610 : !$OMP SHARED (intac, ldsab, ldai, nder, nkind, maxder, ncoset, sap_oce, qs_kind_set, eps_fit) &
611 : !$OMP PRIVATE (sab, work, ai_work, oceh, oces, slot, ikind, jkind, atom_a, atom_b, ilist, jneighbor, rab, cell_b, &
612 : !$OMP iac, dab, qs_kind, orb_basis_set, first_sgfa, la_max, la_min, ncoa_sum, maxsoa, npgfa, nseta, &
613 : !$OMP nsgfa, nsgf_seta, rpgfa, set_radius_a, sphi_a, zeta, paw_proj_b, paw_atom_b, orb_basis_paw, &
614 : !$OMP maxsob, nsetb, mlprj, nprjla, nsatbas, rcprj, fp_car, fp_sph, rzetprj, zetb, nsobtot, clist, &
615 2352 : !$OMP sgf_list, at_a, at_b, local, i, sgf_soft_only, nsgf_cnt)
616 :
617 : ALLOCATE (sab(ldsab, ldsab*maxder), work(ldsab, ldsab*maxder))
618 : sab = 0.0_dp
619 : ALLOCATE (ai_work(ldai, ldai, ncoset(nder + 1)))
620 : ai_work = 0.0_dp
621 : ALLOCATE (oceh(maxder), oces(maxder))
622 :
623 : !$OMP DO SCHEDULE(GUIDED)
624 : DO slot = 1, sap_oce(1)%nl_size
625 : ikind = sap_oce(1)%nlist_task(slot)%ikind
626 : jkind = sap_oce(1)%nlist_task(slot)%jkind
627 : atom_a = sap_oce(1)%nlist_task(slot)%iatom
628 : atom_b = sap_oce(1)%nlist_task(slot)%jatom
629 : ilist = sap_oce(1)%nlist_task(slot)%ilist
630 : jneighbor = sap_oce(1)%nlist_task(slot)%inode
631 : rab(1:3) = sap_oce(1)%nlist_task(slot)%r(1:3)
632 : cell_b(1:3) = sap_oce(1)%nlist_task(slot)%cell(1:3)
633 :
634 : iac = ikind + nkind*(jkind - 1)
635 : dab = SQRT(SUM(rab*rab))
636 :
637 : qs_kind => qs_kind_set(ikind)
638 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis_set)
639 :
640 : IF (.NOT. ASSOCIATED(orb_basis_set)) CYCLE
641 : CALL get_gto_basis_set(gto_basis_set=orb_basis_set, &
642 : first_sgf=first_sgfa, &
643 : lmax=la_max, &
644 : lmin=la_min, &
645 : nco_sum=ncoa_sum, &
646 : maxso=maxsoa, &
647 : npgf=npgfa, &
648 : nset=nseta, &
649 : nsgf=nsgfa, &
650 : nsgf_set=nsgf_seta, &
651 : pgf_radius=rpgfa, &
652 : set_radius=set_radius_a, &
653 : sphi=sphi_a, &
654 : zet=zeta)
655 :
656 : qs_kind => qs_kind_set(jkind)
657 :
658 : NULLIFY (paw_proj_b)
659 : CALL get_qs_kind(qs_kind=qs_kind, paw_proj_set=paw_proj_b, paw_atom=paw_atom_b)
660 : IF (.NOT. paw_atom_b) CYCLE
661 :
662 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=orb_basis_paw, basis_type="GAPW_1C")
663 : IF (.NOT. ASSOCIATED(orb_basis_paw)) CYCLE
664 : CALL get_gto_basis_set(gto_basis_set=orb_basis_paw, maxso=maxsob, nset=nsetb)
665 :
666 : CALL get_paw_proj_set(paw_proj_set=paw_proj_b, &
667 : maxl=mlprj, &
668 : nprj=nprjla, &
669 : nsatbas=nsatbas, &
670 : rcprj=rcprj, &
671 : first_prj=fp_car, &
672 : first_prjs=fp_sph, &
673 : rzetprj=rzetprj, &
674 : zetprj=zetb)
675 :
676 : nsobtot = nsatbas
677 :
678 : clist => intac(iac)%alist(ilist)%clist(jneighbor)
679 : clist%catom = atom_b
680 : clist%cell = cell_b
681 : clist%rac = rab
682 : clist%nsgf_cnt = 0
683 : clist%maxac = 0.0_dp
684 : clist%maxach = 0.0_dp
685 : NULLIFY (clist%acint, clist%achint, clist%sgf_list)
686 :
687 : ALLOCATE (sgf_list(nsgfa))
688 :
689 : at_a => qs_kind_set(jkind)
690 : at_b => qs_kind_set(ikind)
691 :
692 : local = (atom_a == atom_b .AND. ALL(cell_b == 0))
693 :
694 : IF (local) THEN
695 : DO i = 1, maxder
696 : ALLOCATE (oceh(i)%block(nsobtot, nsgfa), oces(i)%block(nsobtot, nsgfa))
697 : oceh(i)%block = 0._dp
698 : oces(i)%block = 0._dp
699 : END DO
700 : CALL build_oce_block_local(oceh, oces, at_a, sgf_list, nsgf_cnt)
701 : clist%nsgf_cnt = nsgf_cnt
702 : clist%sgf_soft_only = .FALSE.
703 : IF (nsgf_cnt > 0) THEN
704 : ALLOCATE (clist%acint(nsgf_cnt, nsobtot, maxder), clist%sgf_list(nsgf_cnt))
705 : clist%acint(:, :, :) = 0._dp
706 : clist%sgf_list(:) = HUGE(0)
707 : CPASSERT(nsgf_cnt == nsgfa)
708 : ! *** Special case: A=B
709 : ALLOCATE (clist%achint(nsgfa, nsobtot, maxder))
710 : clist%achint = 0._dp
711 : clist%acint(1:nsgfa, 1:nsobtot, 1) = TRANSPOSE(oces(1)%block(1:nsobtot, 1:nsgfa))
712 : clist%achint(1:nsgfa, 1:nsobtot, 1) = TRANSPOSE(oceh(1)%block(1:nsobtot, 1:nsgfa))
713 : clist%maxac = MAXVAL(ABS(clist%acint(:, :, 1)))
714 : clist%maxach = 0._dp
715 : clist%sgf_list(1:nsgf_cnt) = sgf_list(1:nsgf_cnt)
716 : END IF
717 : DO i = 1, maxder
718 : DEALLOCATE (oceh(i)%block, oces(i)%block)
719 : END DO
720 : ELSE
721 : DO i = 1, maxder
722 : ALLOCATE (oces(i)%block(nsobtot, nsgfa))
723 : oces(i)%block = 0._dp
724 : END DO
725 : CALL build_oce_block(oces, at_a, at_b, rab, nder, sgf_list, nsgf_cnt, sgf_soft_only, eps_fit)
726 : clist%nsgf_cnt = nsgf_cnt
727 : clist%sgf_soft_only = sgf_soft_only
728 : IF (nsgf_cnt > 0) THEN
729 : ALLOCATE (clist%acint(nsgf_cnt, nsobtot, maxder), clist%sgf_list(nsgf_cnt))
730 : clist%acint(:, :, :) = 0._dp
731 : clist%sgf_list(:) = HUGE(0)
732 : DO i = 1, maxder
733 : clist%acint(1:nsgf_cnt, 1:nsobtot, i) = TRANSPOSE(oces(i)%block(1:nsobtot, 1:nsgf_cnt))
734 : END DO
735 : clist%maxac = MAXVAL(ABS(clist%acint(:, :, 1)))
736 : clist%maxach = 0._dp
737 : clist%sgf_list(1:nsgf_cnt) = sgf_list(1:nsgf_cnt)
738 : END IF
739 : DO i = 1, maxder
740 : DEALLOCATE (oces(i)%block)
741 : END DO
742 : END IF
743 :
744 : DEALLOCATE (sgf_list)
745 :
746 : END DO
747 :
748 : DEALLOCATE (sab, work, ai_work)
749 : DEALLOCATE (oceh, oces)
750 : !$OMP END PARALLEL
751 :
752 : ! Setup sort index
753 2352 : CALL sap_sort(intac)
754 :
755 : END IF
756 :
757 2352 : CALL timestop(handle)
758 :
759 4704 : END SUBROUTINE build_oce_matrices
760 :
761 : ! **************************************************************************************************
762 : !> \brief Project a matrix block onto the local atomic functions.
763 : !>
764 : !> \param h_a ...
765 : !> \param s_a ...
766 : !> \param na ...
767 : !> \param h_b ...
768 : !> \param s_b ...
769 : !> \param nb ...
770 : !> \param blk ...
771 : !> \param ldb ...
772 : !> \param proj_h ...
773 : !> \param proj_s ...
774 : !> \param nso ...
775 : !> \param len1 ...
776 : !> \param len2 ...
777 : !> \param fac ...
778 : !> \param distab ...
779 : !> \par History
780 : !> 02.2009 created
781 : !> 09.2016 use automatic arrays [M Tucker]
782 : ! **************************************************************************************************
783 4489514 : SUBROUTINE proj_blk(h_a, s_a, na, h_b, s_b, nb, blk, ldb, proj_h, proj_s, nso, len1, len2, fac, distab)
784 :
785 : INTEGER, INTENT(IN) :: na
786 : REAL(KIND=dp), INTENT(IN) :: s_a(na, *), h_a(na, *)
787 : INTEGER, INTENT(IN) :: nb
788 : REAL(KIND=dp), INTENT(IN) :: s_b(nb, *), h_b(nb, *)
789 : INTEGER, INTENT(IN) :: ldb
790 : REAL(KIND=dp), INTENT(IN) :: blk(ldb, *)
791 : INTEGER, INTENT(IN) :: nso
792 : REAL(KIND=dp), INTENT(INOUT) :: proj_s(nso, *), proj_h(nso, *)
793 : INTEGER, INTENT(IN) :: len1, len2
794 : REAL(KIND=dp), INTENT(IN) :: fac
795 : LOGICAL, INTENT(IN) :: distab
796 :
797 0 : REAL(KIND=dp) :: buf1(len1), buf2(len2)
798 :
799 4489514 : IF (na .EQ. 0 .OR. nb .EQ. 0 .OR. nso .EQ. 0) RETURN
800 :
801 : ! Handle special cases
802 4489514 : IF (na .EQ. 1 .AND. nb .EQ. 1) THEN
803 : ! hard
804 225229 : CALL dger(nso, nso, fac*blk(1, 1), h_a(1, 1), 1, h_b(1, 1), 1, proj_h(1, 1), nso)
805 : ! soft
806 225229 : CALL dger(nso, nso, fac*blk(1, 1), s_a(1, 1), 1, s_b(1, 1), 1, proj_s(1, 1), nso)
807 : ELSE
808 4264285 : IF (distab) THEN
809 : ! hard
810 3947652 : CALL dgemm('N', 'N', na, nso, nb, fac, blk(1, 1), ldb, h_b(1, 1), nb, 0.0_dp, buf1(1), na)
811 3947652 : CALL dgemm('T', 'N', nso, nso, na, 1.0_dp, h_a(1, 1), na, buf1(1), na, 0.0_dp, buf2(1), nso)
812 3947652 : CALL daxpy(nso*nso, 1.0_dp, buf2(1), 1, proj_h(1, 1), 1)
813 : ! soft
814 3947652 : CALL daxpy(nso*nso, 1.0_dp, buf2(1), 1, proj_s(1, 1), 1)
815 : ELSE
816 : ! hard
817 316633 : CALL dgemm('N', 'N', na, nso, nb, fac, blk(1, 1), ldb, h_b(1, 1), nb, 0.0_dp, buf1(1), na)
818 316633 : CALL dgemm('T', 'N', nso, nso, na, 1.0_dp, h_a(1, 1), na, buf1(1), na, 1.0_dp, proj_h(1, 1), nso)
819 : ! soft
820 316633 : CALL dgemm('N', 'N', na, nso, nb, fac, blk(1, 1), ldb, s_b(1, 1), nb, 0.0_dp, buf1(1), na)
821 316633 : CALL dgemm('T', 'N', nso, nso, na, 1.0_dp, s_a(1, 1), na, buf1(1), na, 1.0_dp, proj_s(1, 1), nso)
822 : END IF
823 : END IF
824 :
825 4489514 : END SUBROUTINE proj_blk
826 :
827 : ! **************************************************************************************************
828 : !> \brief ...
829 : !> \param ain matrix in old indexing
830 : !> \param aout matrix in new compressed indexing
831 : !> \param atom ...
832 : ! **************************************************************************************************
833 71344 : SUBROUTINE prj_gather(ain, aout, atom)
834 :
835 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: ain
836 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: aout
837 : TYPE(qs_kind_type), INTENT(IN) :: atom
838 :
839 : INTEGER :: i, ip, j, jp, nbas
840 71344 : INTEGER, DIMENSION(:), POINTER :: n2oindex
841 : TYPE(gto_basis_set_type), POINTER :: basis_1c
842 :
843 71344 : NULLIFY (basis_1c)
844 71344 : CALL get_qs_kind(qs_kind=atom, basis_set=basis_1c, basis_type="GAPW_1C")
845 71344 : NULLIFY (n2oindex)
846 71344 : CALL get_paw_basis_info(basis_1c, n2oindex=n2oindex, nsatbas=nbas)
847 :
848 1141034 : DO i = 1, nbas
849 1069690 : ip = n2oindex(i)
850 28134428 : DO j = 1, nbas
851 26993394 : jp = n2oindex(j)
852 28063084 : aout(j, i) = ain(jp, ip)
853 : END DO
854 : END DO
855 :
856 71344 : DEALLOCATE (n2oindex)
857 :
858 71344 : END SUBROUTINE prj_gather
859 :
860 : ! **************************************************************************************************
861 : !> \brief ...
862 : !> \param ain matrix in new compressed indexing
863 : !> \param aout matrix in old indexing (addup)
864 : !> \param atom ...
865 : ! **************************************************************************************************
866 115332 : SUBROUTINE prj_scatter(ain, aout, atom)
867 :
868 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: ain
869 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT) :: aout
870 : TYPE(qs_kind_type), INTENT(IN) :: atom
871 :
872 : INTEGER :: i, ip, j, jp, nbas
873 115332 : INTEGER, DIMENSION(:), POINTER :: n2oindex
874 : TYPE(gto_basis_set_type), POINTER :: basis_1c
875 :
876 115332 : NULLIFY (basis_1c)
877 115332 : CALL get_qs_kind(qs_kind=atom, basis_set=basis_1c, basis_type="GAPW_1C")
878 115332 : NULLIFY (n2oindex)
879 115332 : CALL get_paw_basis_info(basis_1c, n2oindex=n2oindex, nsatbas=nbas)
880 :
881 1904472 : DO i = 1, nbas
882 1789140 : ip = n2oindex(i)
883 50014068 : DO j = 1, nbas
884 48109596 : jp = n2oindex(j)
885 49898736 : aout(jp, ip) = aout(jp, ip) + ain(j, i)
886 : END DO
887 : END DO
888 :
889 115332 : DEALLOCATE (n2oindex)
890 :
891 115332 : END SUBROUTINE prj_scatter
892 :
893 168496 : END MODULE qs_oce_methods
|
