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 GAPW reciprocal-space reconstruction and its discrete adjoint.
10 : !> \author Matthias Krack (original density reconstruction)
11 : ! **************************************************************************************************
12 : MODULE gapw_gspace_reconstruction
13 : USE atomic_kind_types, ONLY: atomic_kind_type,&
14 : get_atomic_kind
15 : USE basis_set_types, ONLY: get_gto_basis_set,&
16 : gto_basis_set_type
17 : USE cell_types, ONLY: cell_type,&
18 : pbc
19 : USE cp_control_types, ONLY: dft_control_type
20 : USE cp_log_handling, ONLY: cp_get_default_logger,&
21 : cp_logger_get_default_io_unit,&
22 : cp_logger_type
23 : USE kinds, ONLY: dp
24 : USE mathconstants, ONLY: pi
25 : USE memory_utilities, ONLY: reallocate
26 : USE message_passing, ONLY: mp_para_env_type
27 : USE orbital_pointers, ONLY: indco,&
28 : nco,&
29 : ncoset,&
30 : nso,&
31 : nsoset
32 : USE orbital_transformation_matrices, ONLY: orbtramat
33 : USE particle_types, ONLY: particle_type
34 : USE paw_basis_types, ONLY: get_paw_basis_info
35 : USE pw_grid_types, ONLY: HALFSPACE
36 : USE pw_grids, ONLY: get_pw_grid_info
37 : USE pw_methods, ONLY: pw_axpy,&
38 : pw_integrate_function,&
39 : pw_scale,&
40 : pw_transfer,&
41 : pw_zero
42 : USE pw_pool_types, ONLY: pw_pool_type
43 : USE pw_types, ONLY: pw_c1d_gs_type,&
44 : pw_r3d_rs_type
45 : USE qs_environment_types, ONLY: get_qs_env,&
46 : qs_environment_type
47 : USE qs_kind_types, ONLY: get_qs_kind,&
48 : qs_kind_type
49 : USE qs_rho_atom_types, ONLY: get_rho_atom,&
50 : rho_atom_coeff,&
51 : rho_atom_type
52 : USE qs_rho_types, ONLY: qs_rho_get,&
53 : qs_rho_type
54 : #include "./base/base_uses.f90"
55 :
56 : IMPLICIT NONE
57 :
58 : PRIVATE
59 :
60 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gapw_gspace_reconstruction'
61 :
62 : PUBLIC :: calculate_rhotot_elec_gspace, &
63 : integrate_gapw_composite_vxc_gspace
64 :
65 : CONTAINS
66 :
67 : ! **************************************************************************************************
68 : !> \brief The total electronic density in reciprocal space (g-space) is
69 : !> calculated.
70 : !> \param qs_env ...
71 : !> \param auxbas_pw_pool ...
72 : !> \param rhotot_elec_gspace ...
73 : !> \param q_max ...
74 : !> \param rho_hard ...
75 : !> \param rho_soft ...
76 : !> \param fsign ...
77 : !> \param compute_tau reconstruct the kinetic-energy density instead of rho
78 : !> \param rho_source density structure providing the smooth rho or tau field
79 : !> \param allow_nonorthorhombic permit a general cell for non-XRD reconstruction
80 : !> \date 14.03.2008 (splitted from the routine xray_diffraction_spectrum)
81 : !> \author Matthias Krack
82 : !> \note This code assumes that the g-vectors are ordered (in gsq and %cc)
83 : ! **************************************************************************************************
84 708 : SUBROUTINE calculate_rhotot_elec_gspace(qs_env, auxbas_pw_pool, &
85 : rhotot_elec_gspace, q_max, rho_hard, &
86 : rho_soft, fsign, compute_tau, rho_source, &
87 : allow_nonorthorhombic)
88 :
89 : TYPE(qs_environment_type), POINTER :: qs_env
90 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
91 : TYPE(pw_c1d_gs_type), INTENT(INOUT) :: rhotot_elec_gspace
92 : REAL(KIND=dp), INTENT(IN) :: q_max
93 : REAL(KIND=dp), INTENT(OUT) :: rho_hard, rho_soft
94 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: fsign
95 : LOGICAL, INTENT(IN), OPTIONAL :: compute_tau
96 : TYPE(qs_rho_type), OPTIONAL, POINTER :: rho_source
97 : LOGICAL, INTENT(IN), OPTIONAL :: allow_nonorthorhombic
98 :
99 : CHARACTER(LEN=*), PARAMETER :: routineN = 'calculate_rhotot_elec_gspace'
100 :
101 : INTEGER :: atom, handle, iatom, ico, ico1_pgf, ico1_set, ikind, ipgf, iset, iso, iso1_pgf, &
102 : iso1_set, ison, ispin, jco, jco1_pgf, jco1_set, jpgf, jset, jso, jso1_pgf, jso1_set, &
103 : json, la, lb, maxco, maxso, na, natom, nb, ncoa, ncob, ncotot, nkind, nsatbas, nset, &
104 : nsoa, nsob, nsotot, nspin
105 236 : INTEGER, DIMENSION(:), POINTER :: atom_list, lmax, lmin, npgf, o2nindex
106 : LOGICAL :: my_compute_tau, orthorhombic, paw_atom, &
107 : permit_nonorthorhombic
108 : REAL(KIND=dp) :: alpha, eps_rho_gspace, rho_total, scale, &
109 : volume
110 : REAL(KIND=dp), DIMENSION(3) :: ra
111 236 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: delta_cpc, pab, work, zet
112 236 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
113 : TYPE(cell_type), POINTER :: cell
114 : TYPE(dft_control_type), POINTER :: dft_control
115 : TYPE(gto_basis_set_type), POINTER :: basis_1c_set
116 236 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
117 : TYPE(pw_c1d_gs_type) :: rho_elec_gspace
118 236 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, source_r, tau_r
119 236 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
120 : TYPE(qs_rho_type), POINTER :: rho
121 236 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: cpc_h, cpc_s
122 236 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
123 : TYPE(rho_atom_type), POINTER :: rho_atom
124 :
125 0 : CPASSERT(ASSOCIATED(qs_env))
126 236 : CPASSERT(ASSOCIATED(auxbas_pw_pool))
127 :
128 236 : CALL timeset(routineN, handle)
129 :
130 236 : NULLIFY (atom_list)
131 236 : NULLIFY (atomic_kind_set)
132 236 : NULLIFY (qs_kind_set)
133 236 : NULLIFY (cell)
134 236 : NULLIFY (cpc_h)
135 236 : NULLIFY (cpc_s)
136 236 : NULLIFY (delta_cpc)
137 236 : NULLIFY (dft_control)
138 236 : NULLIFY (lmax)
139 236 : NULLIFY (lmin)
140 236 : NULLIFY (npgf)
141 236 : NULLIFY (basis_1c_set)
142 236 : NULLIFY (pab)
143 236 : NULLIFY (particle_set)
144 236 : NULLIFY (rho, rho_r, source_r, tau_r)
145 236 : NULLIFY (rho_atom)
146 236 : NULLIFY (rho_atom_set)
147 236 : NULLIFY (work)
148 236 : NULLIFY (zet)
149 :
150 : CALL get_qs_env(qs_env=qs_env, &
151 : atomic_kind_set=atomic_kind_set, &
152 : qs_kind_set=qs_kind_set, &
153 : cell=cell, &
154 : dft_control=dft_control, &
155 : particle_set=particle_set, &
156 236 : rho_atom_set=rho_atom_set)
157 :
158 236 : IF (PRESENT(rho_source)) THEN
159 192 : rho => rho_source
160 : ELSE
161 44 : CALL get_qs_env(qs_env=qs_env, rho=rho)
162 : END IF
163 :
164 236 : CALL qs_rho_get(rho, rho_r=rho_r, tau_r=tau_r)
165 236 : my_compute_tau = .FALSE.
166 236 : IF (PRESENT(compute_tau)) my_compute_tau = compute_tau
167 236 : IF (my_compute_tau) THEN
168 96 : CPASSERT(ASSOCIATED(tau_r))
169 96 : source_r => tau_r
170 : ELSE
171 140 : source_r => rho_r
172 : END IF
173 236 : eps_rho_gspace = dft_control%qs_control%eps_rho_gspace
174 236 : nkind = SIZE(atomic_kind_set)
175 236 : nspin = dft_control%nspins
176 :
177 : ! Load the soft contribution of the electronic density
178 :
179 236 : CALL auxbas_pw_pool%create_pw(pw=rho_elec_gspace)
180 :
181 236 : CALL pw_zero(rhotot_elec_gspace)
182 :
183 484 : DO ispin = 1, nspin
184 248 : CALL pw_zero(rho_elec_gspace)
185 248 : CALL pw_transfer(source_r(ispin), rho_elec_gspace)
186 248 : IF (PRESENT(fsign) .AND. (ispin == 2)) THEN
187 6 : alpha = fsign
188 : ELSE
189 242 : alpha = 1.0_dp
190 : END IF
191 484 : CALL pw_axpy(rho_elec_gspace, rhotot_elec_gspace, alpha=alpha)
192 : END DO
193 :
194 : ! Release the auxiliary PW grid for the calculation of the soft
195 : ! contribution
196 :
197 236 : CALL auxbas_pw_pool%give_back_pw(rho_elec_gspace)
198 :
199 236 : rho_soft = pw_integrate_function(rhotot_elec_gspace, isign=-1)
200 :
201 : CALL get_pw_grid_info(pw_grid=rhotot_elec_gspace%pw_grid, vol=volume, &
202 236 : orthorhombic=orthorhombic)
203 236 : permit_nonorthorhombic = .FALSE.
204 236 : IF (PRESENT(allow_nonorthorhombic)) permit_nonorthorhombic = allow_nonorthorhombic
205 236 : IF (.NOT. orthorhombic .AND. .NOT. permit_nonorthorhombic) THEN
206 : CALL cp_abort(__LOCATION__, &
207 0 : "The calculation of XRD spectra for non-orthorhombic cells is not implemented")
208 : END IF
209 :
210 236 : CALL pw_scale(rhotot_elec_gspace, volume)
211 :
212 : ! Add the hard contribution of the electronic density
213 :
214 : ! Each process has to loop over all PAW atoms, since the g-space grid
215 : ! is already distributed over all processes
216 :
217 552 : DO ikind = 1, nkind
218 :
219 : CALL get_atomic_kind(atomic_kind_set(ikind), &
220 : atom_list=atom_list, &
221 316 : natom=natom)
222 :
223 : CALL get_qs_kind(qs_kind_set(ikind), &
224 : basis_set=basis_1c_set, &
225 : basis_type="GAPW_1C", &
226 316 : paw_atom=paw_atom)
227 :
228 316 : IF (.NOT. paw_atom) CYCLE ! no PAW atom: nothing to do
229 :
230 264 : CALL get_paw_basis_info(basis_1c_set, o2nindex=o2nindex, nsatbas=nsatbas)
231 :
232 : CALL get_gto_basis_set(gto_basis_set=basis_1c_set, &
233 : lmax=lmax, &
234 : lmin=lmin, &
235 : maxco=maxco, &
236 : maxso=maxso, &
237 : npgf=npgf, &
238 : nset=nset, &
239 264 : zet=zet)
240 :
241 264 : ncotot = maxco*nset
242 264 : nsotot = maxso*nset
243 264 : CALL reallocate(delta_cpc, 1, nsatbas, 1, nsatbas)
244 264 : CALL reallocate(pab, 1, ncotot, 1, ncotot)
245 264 : CALL reallocate(work, 1, maxso, 1, maxco)
246 :
247 680 : DO iatom = 1, natom
248 :
249 416 : atom = atom_list(iatom)
250 416 : rho_atom => rho_atom_set(atom)
251 :
252 : CALL get_rho_atom(rho_atom=rho_atom, &
253 : cpc_h=cpc_h, &
254 416 : cpc_s=cpc_s)
255 :
256 416 : ra(:) = pbc(particle_set(atom)%r, cell)
257 :
258 48048 : delta_cpc = 0.0_dp
259 :
260 852 : DO ispin = 1, nspin
261 436 : IF (PRESENT(fsign) .AND. (ispin == 2)) THEN
262 10 : alpha = fsign
263 : ELSE
264 426 : alpha = 1.0_dp
265 : END IF
266 57844 : delta_cpc = delta_cpc + alpha*(cpc_h(ispin)%r_coef - cpc_s(ispin)%r_coef)
267 : END DO
268 :
269 416 : scale = 1.0_dp
270 :
271 1248 : DO iset = 1, nset
272 568 : ico1_set = (iset - 1)*maxco + 1
273 568 : iso1_set = (iset - 1)*maxso + 1
274 568 : ncoa = ncoset(lmax(iset))
275 568 : nsoa = nsoset(lmax(iset))
276 1968 : DO jset = 1, nset
277 984 : jco1_set = (jset - 1)*maxco + 1
278 984 : jso1_set = (jset - 1)*maxso + 1
279 984 : ncob = ncoset(lmax(jset))
280 984 : nsob = nsoset(lmax(jset))
281 4644 : DO ipgf = 1, npgf(iset)
282 3092 : ico1_pgf = ico1_set + (ipgf - 1)*ncoa
283 3092 : iso1_pgf = iso1_set + (ipgf - 1)*nsoa
284 14280 : DO jpgf = 1, npgf(jset)
285 10204 : jco1_pgf = jco1_set + (jpgf - 1)*ncob
286 10204 : jso1_pgf = jso1_set + (jpgf - 1)*nsob
287 10204 : ico = ico1_pgf + ncoset(lmin(iset) - 1)
288 10204 : iso = iso1_pgf + nsoset(lmin(iset) - 1)
289 :
290 : ! Transformation spherical to Cartesian
291 :
292 22352 : DO la = lmin(iset), lmax(iset)
293 12148 : jco = jco1_pgf + ncoset(lmin(jset) - 1)
294 12148 : jso = jso1_pgf + nsoset(lmin(jset) - 1)
295 27336 : DO lb = lmin(jset), lmax(jset)
296 15188 : ison = o2nindex(iso)
297 15188 : json = o2nindex(jso)
298 : CALL dgemm("N", "N", nso(la), nco(lb), nso(lb), 1.0_dp, &
299 : delta_cpc(ison:ison + nso(la) - 1, json), SIZE(delta_cpc, 1), &
300 : orbtramat(lb)%slm, nso(lb), 0.0_dp, work, &
301 15188 : maxso)
302 : CALL dgemm("T", "N", nco(la), nco(lb), nso(la), 1.0_dp, &
303 : orbtramat(la)%slm, nso(la), work, maxso, &
304 15188 : 0.0_dp, pab(ico:ico + nco(la) - 1, jco), SIZE(pab, 1))
305 15188 : jco = jco + nco(lb)
306 27336 : jso = jso + nso(lb)
307 : END DO ! next lb
308 12148 : ico = ico + nco(la)
309 22352 : iso = iso + nso(la)
310 : END DO ! next la
311 :
312 : ! Collocate current product of primitive Cartesian functions
313 :
314 10204 : na = ico1_pgf - 1
315 10204 : nb = jco1_pgf - 1
316 :
317 : CALL collocate_pgf_product_gspace( &
318 : la_max=lmax(iset), &
319 : zeta=zet(ipgf, iset), &
320 : la_min=lmin(iset), &
321 : lb_max=lmax(jset), &
322 : zetb=zet(jpgf, jset), &
323 : lb_min=lmin(jset), &
324 : ra=ra, &
325 : rab=[0.0_dp, 0.0_dp, 0.0_dp], &
326 : rab2=0.0_dp, &
327 : scale=scale, &
328 : pab=pab, &
329 : na=na, &
330 : nb=nb, &
331 : eps_rho_gspace=eps_rho_gspace, &
332 : gsq_max=q_max*q_max, &
333 : pw=rhotot_elec_gspace, &
334 13296 : compute_tau=my_compute_tau)
335 :
336 : END DO ! next primitive Gaussian function "jpgf"
337 : END DO ! next primitive Gaussian function "ipgf"
338 : END DO ! next shell set "jset"
339 : END DO ! next shell set "iset"
340 : END DO ! next atom "iatom" of atomic kind "ikind"
341 1080 : DEALLOCATE (o2nindex)
342 : END DO ! next atomic kind "ikind"
343 :
344 236 : rho_total = pw_integrate_function(rhotot_elec_gspace, isign=-1)/volume
345 :
346 236 : rho_hard = rho_total - rho_soft
347 :
348 : ! Release work storage
349 :
350 236 : IF (ASSOCIATED(delta_cpc)) THEN
351 208 : DEALLOCATE (delta_cpc)
352 : END IF
353 :
354 236 : IF (ASSOCIATED(work)) THEN
355 208 : DEALLOCATE (work)
356 : END IF
357 :
358 236 : IF (ASSOCIATED(pab)) THEN
359 208 : DEALLOCATE (pab)
360 : END IF
361 :
362 236 : CALL timestop(handle)
363 :
364 236 : END SUBROUTINE calculate_rhotot_elec_gspace
365 :
366 : ! **************************************************************************************************
367 : !> \brief Apply the discrete adjoint of the common-grid GAPW hard-minus-soft reconstruction.
368 : !> \param qs_env ...
369 : !> \param auxbas_pw_pool ...
370 : !> \param vxc_rho derivative of the XC energy with respect to rho
371 : !> \param vxc_tau derivative of the XC energy with respect to tau
372 : !> \param atom_force explicit force from translating the one-center reconstruction
373 : !> \param strain_virial derivative of the one-center reconstruction with respect to strain
374 : !> \param one_center_contraction contraction of VXC with the hard-minus-soft fields
375 : !> \param one_center_rho_contraction rho contribution to one_center_contraction
376 : !> \param one_center_tau_contraction tau contribution to one_center_contraction
377 : ! **************************************************************************************************
378 96 : SUBROUTINE integrate_gapw_composite_vxc_gspace(qs_env, auxbas_pw_pool, &
379 96 : vxc_rho, vxc_tau, atom_force, strain_virial, &
380 : one_center_contraction, one_center_rho_contraction, &
381 : one_center_tau_contraction)
382 :
383 : TYPE(qs_environment_type), POINTER :: qs_env
384 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
385 : TYPE(pw_r3d_rs_type), DIMENSION(:), INTENT(IN) :: vxc_rho, vxc_tau
386 : REAL(KIND=dp), DIMENSION(:, :), INTENT(INOUT), &
387 : OPTIONAL :: atom_force, strain_virial
388 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: one_center_contraction, &
389 : one_center_rho_contraction, &
390 : one_center_tau_contraction
391 :
392 : CHARACTER(LEN=*), PARAMETER :: routineN = 'integrate_gapw_composite_vxc_gspace'
393 :
394 : INTEGER :: atom, handle, iatom, ico, ico1_pgf, ico1_set, ikind, ipgf, iset, iso, iso1_pgf, &
395 : iso1_set, ison, ispin, jco, jco1_pgf, jco1_set, jpgf, jset, jso, jso1_pgf, jso1_set, &
396 : json, la, lb, maxco, maxso, na, natom, nb, ncoa, ncob, ncotot, nkind, nsatbas, nset, &
397 : nsoa, nsob, nsotot, nspin, output_unit
398 96 : INTEGER, DIMENSION(:), POINTER :: atom_list, lmax, lmin, npgf, o2nindex
399 : LOGICAL :: adjoint_consistent, paw_atom
400 : REAL(KIND=dp) :: adjoint_cart, adjoint_rho_cart, &
401 : adjoint_spherical, adjoint_tau_cart, &
402 : eps_rho_gspace, q_max
403 : REAL(KIND=dp), DIMENSION(3) :: center_force, ra
404 : REAL(KIND=dp), DIMENSION(3, 3) :: atom_strain_virial
405 96 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: delta_cpc, hab, hab_rho, hab_tau, pab, &
406 96 : v_spherical, work, zet
407 96 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
408 : TYPE(cell_type), POINTER :: cell
409 : TYPE(cp_logger_type), POINTER :: logger
410 : TYPE(dft_control_type), POINTER :: dft_control
411 : TYPE(gto_basis_set_type), POINTER :: basis_1c_set
412 : TYPE(mp_para_env_type), POINTER :: para_env
413 96 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
414 96 : TYPE(pw_c1d_gs_type), ALLOCATABLE, DIMENSION(:) :: vrho_g, vtau_g
415 96 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
416 96 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: cpc_h, cpc_s, int_h, int_s
417 96 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
418 : TYPE(rho_atom_type), POINTER :: rho_atom
419 :
420 96 : NULLIFY (logger)
421 192 : logger => cp_get_default_logger()
422 96 : output_unit = cp_logger_get_default_io_unit(logger)
423 :
424 96 : CALL timeset(routineN, handle)
425 :
426 96 : NULLIFY (atom_list, atomic_kind_set, basis_1c_set, cell, cpc_h, cpc_s, delta_cpc, &
427 96 : dft_control, hab, hab_rho, hab_tau, int_h, int_s, lmax, lmin, npgf, o2nindex, &
428 96 : pab, para_env, &
429 96 : particle_set, qs_kind_set, rho_atom, rho_atom_set, work, zet)
430 96 : NULLIFY (v_spherical)
431 :
432 96 : CPASSERT(ASSOCIATED(qs_env))
433 96 : CPASSERT(ASSOCIATED(auxbas_pw_pool))
434 96 : CPASSERT(SIZE(vxc_rho) == SIZE(vxc_tau))
435 96 : IF (PRESENT(one_center_contraction)) one_center_contraction = 0.0_dp
436 96 : IF (PRESENT(one_center_rho_contraction)) one_center_rho_contraction = 0.0_dp
437 96 : IF (PRESENT(one_center_tau_contraction)) one_center_tau_contraction = 0.0_dp
438 96 : IF (PRESENT(atom_force)) THEN
439 20 : CPASSERT(SIZE(atom_force, 1) == 3)
440 : END IF
441 96 : IF (PRESENT(strain_virial)) THEN
442 20 : CPASSERT(SIZE(strain_virial, 1) == 3)
443 20 : CPASSERT(SIZE(strain_virial, 2) == 3)
444 : END IF
445 :
446 : CALL get_qs_env(qs_env=qs_env, &
447 : atomic_kind_set=atomic_kind_set, &
448 : qs_kind_set=qs_kind_set, &
449 : cell=cell, &
450 : dft_control=dft_control, &
451 : para_env=para_env, &
452 : particle_set=particle_set, &
453 96 : rho_atom_set=rho_atom_set)
454 :
455 96 : eps_rho_gspace = dft_control%qs_control%eps_rho_gspace
456 96 : nkind = SIZE(atomic_kind_set)
457 96 : nspin = SIZE(vxc_rho)
458 :
459 576 : ALLOCATE (vrho_g(nspin), vtau_g(nspin))
460 192 : DO ispin = 1, nspin
461 96 : CALL auxbas_pw_pool%create_pw(vrho_g(ispin))
462 96 : CALL auxbas_pw_pool%create_pw(vtau_g(ispin))
463 96 : CALL pw_transfer(vxc_rho(ispin), vrho_g(ispin))
464 192 : CALL pw_transfer(vxc_tau(ispin), vtau_g(ispin))
465 : END DO
466 255056 : q_max = SQRT(MAXVAL(vrho_g(1)%pw_grid%gsq))
467 :
468 216 : DO ikind = 1, nkind
469 120 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
470 : CALL get_qs_kind(qs_kind_set(ikind), basis_set=basis_1c_set, &
471 120 : basis_type="GAPW_1C", paw_atom=paw_atom)
472 120 : IF (.NOT. paw_atom) CYCLE
473 :
474 120 : CALL get_paw_basis_info(basis_1c_set, o2nindex=o2nindex, nsatbas=nsatbas)
475 : CALL get_gto_basis_set(gto_basis_set=basis_1c_set, &
476 : lmax=lmax, lmin=lmin, maxco=maxco, maxso=maxso, &
477 120 : npgf=npgf, nset=nset, zet=zet)
478 120 : ncotot = maxco*nset
479 120 : nsotot = maxso*nset
480 120 : CALL reallocate(hab, 1, ncotot, 1, ncotot)
481 120 : CALL reallocate(hab_rho, 1, ncotot, 1, ncotot)
482 120 : CALL reallocate(hab_tau, 1, ncotot, 1, ncotot)
483 120 : CALL reallocate(pab, 1, ncotot, 1, ncotot)
484 120 : CALL reallocate(delta_cpc, 1, nsatbas, 1, nsatbas)
485 120 : CALL reallocate(v_spherical, 1, nsotot, 1, nsotot)
486 120 : CALL reallocate(work, 1, maxso, 1, maxco)
487 38576 : pab = 0.0_dp
488 :
489 312 : DO iatom = 1, natom
490 192 : atom = atom_list(iatom)
491 192 : rho_atom => rho_atom_set(atom)
492 : CALL get_rho_atom(rho_atom=rho_atom, cpc_h=cpc_h, cpc_s=cpc_s, &
493 192 : ga_Vlocal_gb_h=int_h, ga_Vlocal_gb_s=int_s)
494 192 : ra(:) = pbc(particle_set(atom)%r, cell)
495 :
496 504 : DO ispin = 1, nspin
497 40504 : hab = 0.0_dp
498 40504 : hab_rho = 0.0_dp
499 40504 : hab_tau = 0.0_dp
500 40504 : pab = 0.0_dp
501 192 : center_force = 0.0_dp
502 192 : atom_strain_virial = 0.0_dp
503 17736 : delta_cpc = cpc_h(ispin)%r_coef - cpc_s(ispin)%r_coef
504 432 : DO iset = 1, nset
505 240 : ico1_set = (iset - 1)*maxco + 1
506 240 : iso1_set = (iset - 1)*maxso + 1
507 240 : ncoa = ncoset(lmax(iset))
508 240 : nsoa = nsoset(lmax(iset))
509 776 : DO jset = 1, nset
510 344 : jco1_set = (jset - 1)*maxco + 1
511 344 : jso1_set = (jset - 1)*maxso + 1
512 344 : ncob = ncoset(lmax(jset))
513 344 : nsob = nsoset(lmax(jset))
514 1724 : DO ipgf = 1, npgf(iset)
515 1140 : ico1_pgf = ico1_set + (ipgf - 1)*ncoa
516 1140 : iso1_pgf = iso1_set + (ipgf - 1)*nsoa
517 5432 : DO jpgf = 1, npgf(jset)
518 3948 : jco1_pgf = jco1_set + (jpgf - 1)*ncob
519 3948 : jso1_pgf = jso1_set + (jpgf - 1)*nsob
520 3948 : ico = ico1_pgf + ncoset(lmin(iset) - 1)
521 3948 : iso = iso1_pgf + nsoset(lmin(iset) - 1)
522 8516 : DO la = lmin(iset), lmax(iset)
523 4568 : jco = jco1_pgf + ncoset(lmin(jset) - 1)
524 4568 : jso = jso1_pgf + nsoset(lmin(jset) - 1)
525 10176 : DO lb = lmin(jset), lmax(jset)
526 5608 : ison = o2nindex(iso)
527 5608 : json = o2nindex(jso)
528 : CALL dgemm("N", "N", nso(la), nco(lb), nso(lb), 1.0_dp, &
529 : delta_cpc(ison:ison + nso(la) - 1, json), &
530 : SIZE(delta_cpc, 1), orbtramat(lb)%slm, nso(lb), &
531 5608 : 0.0_dp, work, maxso)
532 : CALL dgemm("T", "N", nco(la), nco(lb), nso(la), 1.0_dp, &
533 : orbtramat(la)%slm, nso(la), work, maxso, 0.0_dp, &
534 5608 : pab(ico:ico + nco(la) - 1, jco), SIZE(pab, 1))
535 5608 : jco = jco + nco(lb)
536 10176 : jso = jso + nso(lb)
537 : END DO
538 4568 : ico = ico + nco(la)
539 8516 : iso = iso + nso(la)
540 : END DO
541 3948 : na = ico1_pgf - 1
542 3948 : nb = jco1_pgf - 1
543 5088 : IF (PRESENT(strain_virial)) THEN
544 : CALL collocate_pgf_product_gspace( &
545 : la_max=lmax(iset), zeta=zet(ipgf, iset), la_min=lmin(iset), &
546 : lb_max=lmax(jset), zetb=zet(jpgf, jset), lb_min=lmin(jset), &
547 : ra=ra, rab=[0.0_dp, 0.0_dp, 0.0_dp], rab2=0.0_dp, scale=1.0_dp, &
548 : pab=pab, na=na, nb=nb, eps_rho_gspace=eps_rho_gspace, &
549 : gsq_max=q_max*q_max, pw=vrho_g(ispin), &
550 : potential_gspace=vrho_g(ispin), hab=hab_rho, &
551 1078 : center_force=center_force, strain_virial=atom_strain_virial)
552 : CALL collocate_pgf_product_gspace( &
553 : la_max=lmax(iset), zeta=zet(ipgf, iset), la_min=lmin(iset), &
554 : lb_max=lmax(jset), zetb=zet(jpgf, jset), lb_min=lmin(jset), &
555 : ra=ra, rab=[0.0_dp, 0.0_dp, 0.0_dp], rab2=0.0_dp, scale=1.0_dp, &
556 : pab=pab, na=na, nb=nb, eps_rho_gspace=eps_rho_gspace, &
557 : gsq_max=q_max*q_max, pw=vtau_g(ispin), compute_tau=.TRUE., &
558 : potential_gspace=vtau_g(ispin), hab=hab_tau, &
559 1078 : center_force=center_force, strain_virial=atom_strain_virial)
560 2870 : ELSE IF (PRESENT(atom_force)) THEN
561 : CALL collocate_pgf_product_gspace( &
562 : la_max=lmax(iset), zeta=zet(ipgf, iset), la_min=lmin(iset), &
563 : lb_max=lmax(jset), zetb=zet(jpgf, jset), lb_min=lmin(jset), &
564 : ra=ra, rab=[0.0_dp, 0.0_dp, 0.0_dp], rab2=0.0_dp, scale=1.0_dp, &
565 : pab=pab, na=na, nb=nb, eps_rho_gspace=eps_rho_gspace, &
566 : gsq_max=q_max*q_max, pw=vrho_g(ispin), &
567 : potential_gspace=vrho_g(ispin), hab=hab_rho, &
568 0 : center_force=center_force)
569 : CALL collocate_pgf_product_gspace( &
570 : la_max=lmax(iset), zeta=zet(ipgf, iset), la_min=lmin(iset), &
571 : lb_max=lmax(jset), zetb=zet(jpgf, jset), lb_min=lmin(jset), &
572 : ra=ra, rab=[0.0_dp, 0.0_dp, 0.0_dp], rab2=0.0_dp, scale=1.0_dp, &
573 : pab=pab, na=na, nb=nb, eps_rho_gspace=eps_rho_gspace, &
574 : gsq_max=q_max*q_max, pw=vtau_g(ispin), compute_tau=.TRUE., &
575 : potential_gspace=vtau_g(ispin), hab=hab_tau, &
576 0 : center_force=center_force)
577 : ELSE
578 : CALL collocate_pgf_product_gspace( &
579 : la_max=lmax(iset), zeta=zet(ipgf, iset), la_min=lmin(iset), &
580 : lb_max=lmax(jset), zetb=zet(jpgf, jset), lb_min=lmin(jset), &
581 : ra=ra, rab=[0.0_dp, 0.0_dp, 0.0_dp], rab2=0.0_dp, scale=1.0_dp, &
582 : pab=pab, na=na, nb=nb, eps_rho_gspace=eps_rho_gspace, &
583 : gsq_max=q_max*q_max, pw=vrho_g(ispin), &
584 2870 : potential_gspace=vrho_g(ispin), hab=hab_rho)
585 : CALL collocate_pgf_product_gspace( &
586 : la_max=lmax(iset), zeta=zet(ipgf, iset), la_min=lmin(iset), &
587 : lb_max=lmax(jset), zetb=zet(jpgf, jset), lb_min=lmin(jset), &
588 : ra=ra, rab=[0.0_dp, 0.0_dp, 0.0_dp], rab2=0.0_dp, scale=1.0_dp, &
589 : pab=pab, na=na, nb=nb, eps_rho_gspace=eps_rho_gspace, &
590 : gsq_max=q_max*q_max, pw=vtau_g(ispin), compute_tau=.TRUE., &
591 2870 : potential_gspace=vtau_g(ispin), hab=hab_tau)
592 : END IF
593 : END DO
594 : END DO
595 : END DO
596 : END DO
597 80816 : CALL para_env%sum(hab_rho)
598 80816 : CALL para_env%sum(hab_tau)
599 : ! Atom forces and virials remain rank-local here. Quickstep performs the
600 : ! corresponding global reductions after all force components are assembled.
601 80816 : hab = hab_rho + hab_tau
602 40504 : adjoint_rho_cart = SUM(pab*hab_rho)
603 40504 : adjoint_tau_cart = SUM(pab*hab_tau)
604 192 : adjoint_cart = adjoint_rho_cart + adjoint_tau_cart
605 192 : IF (PRESENT(one_center_contraction)) THEN
606 192 : one_center_contraction = one_center_contraction + adjoint_cart
607 : END IF
608 192 : IF (PRESENT(one_center_rho_contraction)) THEN
609 192 : one_center_rho_contraction = one_center_rho_contraction + adjoint_rho_cart
610 : END IF
611 192 : IF (PRESENT(one_center_tau_contraction)) THEN
612 192 : one_center_tau_contraction = one_center_tau_contraction + adjoint_tau_cart
613 : END IF
614 :
615 192 : IF (ASSOCIATED(int_h) .AND. ASSOCIATED(int_s)) THEN
616 96 : IF (ASSOCIATED(int_h(ispin)%r_coef) .AND. ASSOCIATED(int_s(ispin)%r_coef)) THEN
617 20252 : v_spherical = 0.0_dp
618 216 : DO iset = 1, nset
619 120 : ico1_set = (iset - 1)*maxco + 1
620 120 : iso1_set = (iset - 1)*maxso + 1
621 120 : ncoa = ncoset(lmax(iset))
622 120 : nsoa = nsoset(lmax(iset))
623 388 : DO jset = 1, nset
624 172 : jco1_set = (jset - 1)*maxco + 1
625 172 : jso1_set = (jset - 1)*maxso + 1
626 172 : ncob = ncoset(lmax(jset))
627 172 : nsob = nsoset(lmax(jset))
628 862 : DO ipgf = 1, npgf(iset)
629 570 : ico1_pgf = ico1_set + (ipgf - 1)*ncoa
630 570 : iso1_pgf = iso1_set + (ipgf - 1)*nsoa
631 2716 : DO jpgf = 1, npgf(jset)
632 1974 : jco1_pgf = jco1_set + (jpgf - 1)*ncob
633 1974 : jso1_pgf = jso1_set + (jpgf - 1)*nsob
634 1974 : ico = ico1_pgf + ncoset(lmin(iset) - 1)
635 1974 : iso = iso1_pgf + nsoset(lmin(iset) - 1)
636 4828 : DO la = lmin(iset), lmax(iset)
637 2284 : jco = jco1_pgf + ncoset(lmin(jset) - 1)
638 2284 : jso = jso1_pgf + nsoset(lmin(jset) - 1)
639 5088 : DO lb = lmin(jset), lmax(jset)
640 : CALL dgemm("N", "N", nso(la), nco(lb), nco(la), 1.0_dp, &
641 : orbtramat(la)%slm, nso(la), &
642 : hab(ico:ico + nco(la) - 1, jco), &
643 2804 : SIZE(hab, 1), 0.0_dp, work, maxso)
644 : CALL dgemm("N", "T", nso(la), nso(lb), nco(lb), 1.0_dp, &
645 : work, maxso, orbtramat(lb)%slm, nso(lb), 1.0_dp, &
646 2804 : v_spherical(iso:iso + nso(la) - 1, jso), nsotot)
647 2804 : jco = jco + nco(lb)
648 5088 : jso = jso + nso(lb)
649 : END DO
650 2284 : ico = ico + nco(la)
651 4258 : iso = iso + nso(la)
652 : END DO
653 : END DO
654 : END DO
655 : END DO
656 : END DO
657 40408 : int_h(ispin)%r_coef = int_h(ispin)%r_coef + v_spherical
658 40408 : int_s(ispin)%r_coef = int_s(ispin)%r_coef + v_spherical
659 96 : adjoint_spherical = 0.0_dp
660 916 : DO jso = 1, nsotot
661 820 : json = o2nindex(jso)
662 820 : IF (json == 0) CYCLE
663 13068 : DO iso = 1, nsotot
664 12332 : ison = o2nindex(iso)
665 12332 : IF (ison == 0) CYCLE
666 : adjoint_spherical = adjoint_spherical + &
667 13152 : delta_cpc(ison, json)*v_spherical(iso, jso)
668 : END DO
669 : END DO
670 : adjoint_consistent = ABS(adjoint_cart - adjoint_spherical) <= &
671 96 : 1.0E-9_dp*MAX(1.0_dp, ABS(adjoint_cart))
672 96 : IF (.NOT. adjoint_consistent) THEN
673 : WRITE (UNIT=output_unit, FMT="(A,3I6,3ES24.15)") &
674 0 : " SKALA_GPW| Composite adjoint mismatch", ikind, atom, ispin, &
675 0 : adjoint_cart, adjoint_spherical, adjoint_cart - adjoint_spherical
676 : END IF
677 96 : CPASSERT(adjoint_consistent)
678 : END IF
679 : END IF
680 312 : IF (PRESENT(atom_force)) atom_force(:, atom) = atom_force(:, atom) + center_force
681 864 : IF (PRESENT(strain_virial)) strain_virial = strain_virial + atom_strain_virial
682 : END DO
683 : END DO
684 :
685 456 : DEALLOCATE (delta_cpc, hab, hab_rho, hab_tau, pab, v_spherical, work, o2nindex)
686 : END DO
687 :
688 192 : DO ispin = 1, nspin
689 96 : CALL auxbas_pw_pool%give_back_pw(vrho_g(ispin))
690 192 : CALL auxbas_pw_pool%give_back_pw(vtau_g(ispin))
691 : END DO
692 96 : DEALLOCATE (vrho_g, vtau_g)
693 :
694 96 : CALL timestop(handle)
695 :
696 192 : END SUBROUTINE integrate_gapw_composite_vxc_gspace
697 :
698 : ! **************************************************************************************************
699 : !> \brief low level collocation of primitive gaussian functions in g-space
700 : !> \param la_max ...
701 : !> \param zeta ...
702 : !> \param la_min ...
703 : !> \param lb_max ...
704 : !> \param zetb ...
705 : !> \param lb_min ...
706 : !> \param ra ...
707 : !> \param rab ...
708 : !> \param rab2 ...
709 : !> \param scale ...
710 : !> \param pab ...
711 : !> \param na ...
712 : !> \param nb ...
713 : !> \param eps_rho_gspace ...
714 : !> \param gsq_max ...
715 : !> \param pw ...
716 : !> \param compute_tau apply the symmetric kinetic-density DADB transformation
717 : !> \param potential_gspace integrate each primitive product against this potential
718 : !> \param hab Cartesian primitive-product potential integrals
719 : !> \param center_force force from translating the collocated product density
720 : !> \param strain_virial derivative of the product density with respect to strain
721 : ! **************************************************************************************************
722 18100 : SUBROUTINE collocate_pgf_product_gspace(la_max, zeta, la_min, &
723 : lb_max, zetb, lb_min, &
724 : ra, rab, rab2, scale, pab, na, nb, &
725 : eps_rho_gspace, gsq_max, pw, compute_tau, &
726 18100 : potential_gspace, hab, center_force, strain_virial)
727 :
728 : INTEGER, INTENT(IN) :: la_max
729 : REAL(dp), INTENT(IN) :: zeta
730 : INTEGER, INTENT(IN) :: la_min, lb_max
731 : REAL(dp), INTENT(IN) :: zetb
732 : INTEGER, INTENT(IN) :: lb_min
733 : REAL(dp), DIMENSION(3), INTENT(IN) :: ra, rab
734 : REAL(dp), INTENT(IN) :: rab2, scale
735 : REAL(dp), DIMENSION(:, :), POINTER :: pab
736 : INTEGER, INTENT(IN) :: na, nb
737 : REAL(dp), INTENT(IN) :: eps_rho_gspace, gsq_max
738 : TYPE(pw_c1d_gs_type), INTENT(IN) :: pw
739 : LOGICAL, INTENT(IN), OPTIONAL :: compute_tau
740 : TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: potential_gspace
741 : REAL(dp), DIMENSION(:, :), INTENT(INOUT), OPTIONAL :: hab
742 : REAL(dp), DIMENSION(3), INTENT(INOUT), OPTIONAL :: center_force
743 : REAL(dp), DIMENSION(3, 3), INTENT(INOUT), OPTIONAL :: strain_virial
744 :
745 : CHARACTER(LEN=*), PARAMETER :: routineN = 'collocate_pgf_product_gspace'
746 :
747 : COMPLEX(dp) :: plane_product, rag, rbg, shifted_product, &
748 : tau_product
749 : COMPLEX(dp), DIMENSION(3) :: plane_gradient
750 18100 : COMPLEX(dp), DIMENSION(:, :, :), ALLOCATABLE :: cubeaxis
751 : INTEGER :: ax, ay, az, bx, by, bz, handle, i, ico, &
752 : idir, jco, jdir, la, la_rec_max, lb, &
753 : lb_grid, lb_rec_max, ub_grid
754 : LOGICAL :: integrate_potential, integrate_strain, &
755 : my_compute_tau
756 : REAL(dp) :: f, fa, fb, gcoord, pij, potential_term, &
757 : prefactor, reciprocal_weight, rzetp, &
758 : twozetp, zetp
759 : REAL(dp), DIMENSION(3) :: fap, fbp, force_product, gvec, rap, rbp, rp
760 :
761 18100 : CALL timeset(routineN, handle)
762 :
763 18100 : my_compute_tau = .FALSE.
764 18100 : IF (PRESENT(compute_tau)) my_compute_tau = compute_tau
765 18100 : integrate_potential = PRESENT(potential_gspace) .OR. PRESENT(hab)
766 18100 : integrate_strain = PRESENT(strain_virial)
767 18100 : CPASSERT(PRESENT(potential_gspace) .EQV. PRESENT(hab))
768 18100 : CPASSERT(.NOT. integrate_strain .OR. integrate_potential)
769 18100 : IF (integrate_potential) THEN
770 7896 : CPASSERT(ASSOCIATED(pw%pw_grid, potential_gspace%pw_grid))
771 : END IF
772 :
773 18100 : la_rec_max = la_max
774 18100 : lb_rec_max = lb_max
775 18100 : IF (my_compute_tau) THEN
776 7896 : la_rec_max = la_rec_max + 1
777 7896 : lb_rec_max = lb_rec_max + 1
778 : END IF
779 18100 : IF (integrate_strain) la_rec_max = la_rec_max + 1
780 :
781 18100 : zetp = zeta + zetb
782 18100 : rzetp = 1.0_dp/zetp
783 18100 : f = zetb*rzetp
784 72400 : rap(:) = f*rab(:)
785 72400 : rbp(:) = rap(:) - rab(:)
786 72400 : rp(:) = ra(:) + rap(:)
787 18100 : twozetp = 2.0_dp*zetp
788 72400 : fap(:) = twozetp*rap(:)
789 72400 : fbp(:) = twozetp*rbp(:)
790 18100 : prefactor = scale*SQRT((pi*rzetp)**3)*EXP(-zeta*f*rab2)
791 :
792 72400 : ALLOCATE (cubeaxis(3, 0:la_rec_max, 0:lb_rec_max))
793 18100 : lb_grid = LBOUND(pw%array, 1)
794 18100 : ub_grid = UBOUND(pw%array, 1)
795 :
796 216263668 : DO i = lb_grid, ub_grid
797 216245568 : IF (pw%pw_grid%gsq(i) > gsq_max) CYCLE
798 773838704 : gvec = pw%pw_grid%g(:, i)
799 193459676 : CALL build_cubeaxis(gvec)
800 193459676 : reciprocal_weight = 1.0_dp
801 193459676 : IF (pw%pw_grid%grid_span == HALFSPACE .AND. pw%pw_grid%gsq(i) > 0.0_dp) THEN
802 0 : reciprocal_weight = 2.0_dp
803 : END IF
804 :
805 637085998 : DO ico = ncoset(la_min - 1) + 1, ncoset(la_max)
806 443608222 : ax = indco(1, ico)
807 443608222 : ay = indco(2, ico)
808 443608222 : az = indco(3, ico)
809 :
810 1789499748 : DO jco = ncoset(lb_min - 1) + 1, ncoset(lb_max)
811 1129645958 : bx = indco(1, jco)
812 1129645958 : by = indco(2, jco)
813 1129645958 : bz = indco(3, jco)
814 1129645958 : pij = prefactor*pab(na + ico, nb + jco)
815 1129645958 : IF (.NOT. integrate_potential .AND. ABS(pij) < eps_rho_gspace) CYCLE
816 :
817 484421456 : IF (.NOT. my_compute_tau) THEN
818 418105414 : plane_product = cubeaxis(1, ax, bx)*cubeaxis(2, ay, by)*cubeaxis(3, az, bz)
819 : ELSE
820 : ! Equivalent to GRID_FUNC_DADB: 0.5 * grad(phi_a).grad(phi_b).
821 : tau_product = 2.0_dp*zeta*zetb* &
822 : cubeaxis(1, ax + 1, bx + 1)* &
823 66316042 : cubeaxis(2, ay, by)*cubeaxis(3, az, bz)
824 66316042 : IF (ax > 0) THEN
825 : tau_product = tau_product - REAL(ax, dp)*zetb* &
826 : cubeaxis(1, ax - 1, bx + 1)* &
827 16506453 : cubeaxis(2, ay, by)*cubeaxis(3, az, bz)
828 : END IF
829 66316042 : IF (bx > 0) THEN
830 : tau_product = tau_product - zeta*REAL(bx, dp)* &
831 : cubeaxis(1, ax + 1, bx - 1)* &
832 16506453 : cubeaxis(2, ay, by)*cubeaxis(3, az, bz)
833 : END IF
834 66316042 : IF (ax > 0 .AND. bx > 0) THEN
835 : tau_product = tau_product + 0.5_dp*REAL(ax*bx, dp)* &
836 : cubeaxis(1, ax - 1, bx - 1)* &
837 4952040 : cubeaxis(2, ay, by)*cubeaxis(3, az, bz)
838 : END IF
839 :
840 : tau_product = tau_product + 2.0_dp*zeta*zetb* &
841 : cubeaxis(1, ax, bx)* &
842 66316042 : cubeaxis(2, ay + 1, by + 1)*cubeaxis(3, az, bz)
843 66316042 : IF (ay > 0) THEN
844 : tau_product = tau_product - REAL(ay, dp)*zetb* &
845 : cubeaxis(1, ax, bx)* &
846 16520276 : cubeaxis(2, ay - 1, by + 1)*cubeaxis(3, az, bz)
847 : END IF
848 66316042 : IF (by > 0) THEN
849 : tau_product = tau_product - zeta*REAL(by, dp)* &
850 : cubeaxis(1, ax, bx)* &
851 16520276 : cubeaxis(2, ay + 1, by - 1)*cubeaxis(3, az, bz)
852 : END IF
853 66316042 : IF (ay > 0 .AND. by > 0) THEN
854 : tau_product = tau_product + 0.5_dp*REAL(ay*by, dp)* &
855 : cubeaxis(1, ax, bx)* &
856 4952040 : cubeaxis(2, ay - 1, by - 1)*cubeaxis(3, az, bz)
857 : END IF
858 :
859 : tau_product = tau_product + 2.0_dp*zeta*zetb* &
860 : cubeaxis(1, ax, bx)*cubeaxis(2, ay, by)* &
861 66316042 : cubeaxis(3, az + 1, bz + 1)
862 66316042 : IF (az > 0) THEN
863 : tau_product = tau_product - REAL(az, dp)*zetb* &
864 : cubeaxis(1, ax, bx)*cubeaxis(2, ay, by)* &
865 16572764 : cubeaxis(3, az - 1, bz + 1)
866 : END IF
867 66316042 : IF (bz > 0) THEN
868 : tau_product = tau_product - zeta*REAL(bz, dp)* &
869 : cubeaxis(1, ax, bx)*cubeaxis(2, ay, by)* &
870 16572764 : cubeaxis(3, az + 1, bz - 1)
871 : END IF
872 66316042 : IF (az > 0 .AND. bz > 0) THEN
873 : tau_product = tau_product + 0.5_dp*REAL(az*bz, dp)* &
874 : cubeaxis(1, ax, bx)*cubeaxis(2, ay, by)* &
875 4952040 : cubeaxis(3, az - 1, bz - 1)
876 : END IF
877 : plane_product = tau_product
878 : END IF
879 :
880 484421456 : IF (integrate_strain) THEN
881 58366088 : IF (.NOT. my_compute_tau) THEN
882 : shifted_product = cubeaxis(1, ax + 1, bx)* &
883 29183044 : cubeaxis(2, ay, by)*cubeaxis(3, az, bz)
884 : plane_gradient(1) = CMPLX(0.0_dp, -1.0_dp, KIND=dp)* &
885 29183044 : (shifted_product + ra(1)*plane_product)
886 : shifted_product = cubeaxis(1, ax, bx)* &
887 29183044 : cubeaxis(2, ay + 1, by)*cubeaxis(3, az, bz)
888 : plane_gradient(2) = CMPLX(0.0_dp, -1.0_dp, KIND=dp)* &
889 29183044 : (shifted_product + ra(2)*plane_product)
890 : shifted_product = cubeaxis(1, ax, bx)* &
891 29183044 : cubeaxis(2, ay, by)*cubeaxis(3, az + 1, bz)
892 : plane_gradient(3) = CMPLX(0.0_dp, -1.0_dp, KIND=dp)* &
893 29183044 : (shifted_product + ra(3)*plane_product)
894 : ELSE
895 116732176 : DO idir = 1, 3
896 87549132 : shifted_product = tau_plane_product_shifted(idir)
897 : plane_gradient(idir) = CMPLX(0.0_dp, -1.0_dp, KIND=dp)* &
898 116732176 : (shifted_product + ra(idir)*plane_product)
899 : END DO
900 : END IF
901 : END IF
902 :
903 928029678 : IF (integrate_potential) THEN
904 : potential_term = reciprocal_weight*prefactor* &
905 122767888 : REAL(CONJG(potential_gspace%array(i))*plane_product, KIND=dp)
906 122767888 : hab(na + ico, nb + jco) = hab(na + ico, nb + jco) + potential_term
907 122767888 : IF (PRESENT(center_force) .OR. integrate_strain) THEN
908 : force_product = reciprocal_weight*pij* &
909 : REAL(CONJG(potential_gspace%array(i))* &
910 233464352 : CMPLX(0.0_dp, -gvec, KIND=dp)*plane_product, KIND=dp)
911 233464352 : IF (PRESENT(center_force)) center_force = center_force + force_product
912 : END IF
913 122767888 : IF (integrate_strain) THEN
914 233464352 : DO idir = 1, 3
915 700393056 : DO jdir = 1, 3
916 : strain_virial(idir, jdir) = strain_virial(idir, jdir) - &
917 : reciprocal_weight*pij*gvec(idir)* &
918 : REAL(CONJG(potential_gspace%array(i))*plane_gradient(jdir), KIND=dp) + &
919 700393056 : force_product(idir)*ra(jdir)
920 : END DO
921 : strain_virial(idir, idir) = strain_virial(idir, idir) - &
922 233464352 : pab(na + ico, nb + jco)*potential_term
923 : END DO
924 : END IF
925 : ELSE
926 361653568 : pw%array(i) = pw%array(i) + pij*plane_product
927 : END IF
928 : END DO
929 : END DO
930 : END DO
931 :
932 18100 : DEALLOCATE (cubeaxis)
933 36200 : CALL timestop(handle)
934 :
935 : CONTAINS
936 :
937 : ! **************************************************************************************************
938 : !> \brief ...
939 : !> \param gvector ...
940 : ! **************************************************************************************************
941 193459676 : SUBROUTINE build_cubeaxis(gvector)
942 : REAL(dp), DIMENSION(3), INTENT(IN) :: gvector
943 :
944 193459676 : cubeaxis = CMPLX(0.0_dp, 0.0_dp, KIND=dp)
945 773838704 : DO idir = 1, 3
946 580379028 : gcoord = gvector(idir)
947 : cubeaxis(idir, 0, 0) = &
948 580379028 : EXP(CMPLX(-0.25_dp*rzetp*gcoord*gcoord, -rp(idir)*gcoord, KIND=dp))
949 :
950 580379028 : IF (la_rec_max > 0) THEN
951 280967070 : rag = CMPLX(fap(idir), -gcoord, KIND=dp)
952 280967070 : cubeaxis(idir, 1, 0) = rag*cubeaxis(idir, 0, 0)
953 361962780 : DO la = 2, la_rec_max
954 80995710 : fa = REAL(la - 1, dp)*twozetp
955 : cubeaxis(idir, la, 0) = rag*cubeaxis(idir, la - 1, 0) + &
956 361962780 : fa*cubeaxis(idir, la - 2, 0)
957 : END DO
958 280967070 : IF (lb_rec_max > 0) THEN
959 170411154 : rbg = CMPLX(fbp(idir), -gcoord, KIND=dp)
960 170411154 : fa = twozetp
961 170411154 : cubeaxis(idir, 0, 1) = rbg*cubeaxis(idir, 0, 0)
962 : cubeaxis(idir, 1, 1) = rbg*cubeaxis(idir, 1, 0) + &
963 170411154 : fa*cubeaxis(idir, 0, 0)
964 213351516 : DO lb = 2, lb_rec_max
965 42940362 : fb = REAL(lb - 1, dp)*twozetp
966 : cubeaxis(idir, 0, lb) = rbg*cubeaxis(idir, 0, lb - 1) + &
967 42940362 : fb*cubeaxis(idir, 0, lb - 2)
968 : cubeaxis(idir, 1, lb) = rbg*cubeaxis(idir, 1, lb - 1) + &
969 : fb*cubeaxis(idir, 1, lb - 2) + &
970 213351516 : fa*cubeaxis(idir, 0, lb - 1)
971 : END DO
972 230326356 : DO la = 2, la_rec_max
973 59915202 : fa = REAL(la, dp)*twozetp
974 : cubeaxis(idir, la, 1) = rbg*cubeaxis(idir, la, 0) + &
975 59915202 : fa*cubeaxis(idir, la - 1, 0)
976 263130714 : DO lb = 2, lb_rec_max
977 32804358 : fb = REAL(lb - 1, dp)*twozetp
978 : cubeaxis(idir, la, lb) = rbg*cubeaxis(idir, la, lb - 1) + &
979 : fb*cubeaxis(idir, la, lb - 2) + &
980 92719560 : fa*cubeaxis(idir, la - 1, lb - 1)
981 : END DO
982 : END DO
983 : END IF
984 299411958 : ELSE IF (lb_rec_max > 0) THEN
985 103235724 : rbg = CMPLX(fbp(idir), -gcoord, KIND=dp)
986 103235724 : cubeaxis(idir, 0, 1) = rbg*cubeaxis(idir, 0, 0)
987 123232824 : DO lb = 2, lb_rec_max
988 19997100 : fb = REAL(lb - 1, dp)*twozetp
989 : cubeaxis(idir, 0, lb) = rbg*cubeaxis(idir, 0, lb - 1) + &
990 123232824 : fb*cubeaxis(idir, 0, lb - 2)
991 : END DO
992 : END IF
993 :
994 1716180512 : DO la = 0, la_rec_max
995 3107718060 : DO lb = 0, lb_rec_max
996 1584997224 : IF (la + lb == 0) CYCLE
997 1004618196 : fa = (1.0_dp/twozetp)**(la + lb)
998 2527339032 : cubeaxis(idir, la, lb) = fa*cubeaxis(idir, la, lb)
999 : END DO
1000 : END DO
1001 : END DO
1002 193459676 : END SUBROUTINE build_cubeaxis
1003 :
1004 : ! **************************************************************************************************
1005 : !> \brief ...
1006 : !> \param shift_axis ...
1007 : !> \return ...
1008 : ! **************************************************************************************************
1009 87549132 : FUNCTION tau_plane_product_shifted(shift_axis) RESULT(value)
1010 : INTEGER, INTENT(IN) :: shift_axis
1011 : COMPLEX(dp) :: value
1012 :
1013 : INTEGER :: sx, sy, sz
1014 :
1015 87549132 : sx = 0
1016 87549132 : sy = 0
1017 87549132 : sz = 0
1018 87549132 : IF (shift_axis == 1) sx = 1
1019 87549132 : IF (shift_axis == 2) sy = 1
1020 87549132 : IF (shift_axis == 3) sz = 1
1021 :
1022 : value = 2.0_dp*zeta*zetb* &
1023 : cubeaxis(1, ax + 1 + sx, bx + 1)* &
1024 87549132 : cubeaxis(2, ay + sy, by)*cubeaxis(3, az + sz, bz)
1025 87549132 : IF (ax > 0) THEN
1026 : value = value - REAL(ax, dp)*zetb* &
1027 : cubeaxis(1, ax - 1 + sx, bx + 1)* &
1028 24120108 : cubeaxis(2, ay + sy, by)*cubeaxis(3, az + sz, bz)
1029 : END IF
1030 87549132 : IF (bx > 0) THEN
1031 : value = value - zeta*REAL(bx, dp)* &
1032 : cubeaxis(1, ax + 1 + sx, bx - 1)* &
1033 24120108 : cubeaxis(2, ay + sy, by)*cubeaxis(3, az + sz, bz)
1034 : END IF
1035 87549132 : IF (ax > 0 .AND. bx > 0) THEN
1036 : value = value + 0.5_dp*REAL(ax*bx, dp)* &
1037 : cubeaxis(1, ax - 1 + sx, bx - 1)* &
1038 7165620 : cubeaxis(2, ay + sy, by)*cubeaxis(3, az + sz, bz)
1039 : END IF
1040 :
1041 : value = value + 2.0_dp*zeta*zetb* &
1042 : cubeaxis(1, ax + sx, bx)* &
1043 87549132 : cubeaxis(2, ay + 1 + sy, by + 1)*cubeaxis(3, az + sz, bz)
1044 87549132 : IF (ay > 0) THEN
1045 : value = value - REAL(ay, dp)*zetb* &
1046 : cubeaxis(1, ax + sx, bx)* &
1047 24120108 : cubeaxis(2, ay - 1 + sy, by + 1)*cubeaxis(3, az + sz, bz)
1048 : END IF
1049 87549132 : IF (by > 0) THEN
1050 : value = value - zeta*REAL(by, dp)* &
1051 : cubeaxis(1, ax + sx, bx)* &
1052 24120108 : cubeaxis(2, ay + 1 + sy, by - 1)*cubeaxis(3, az + sz, bz)
1053 : END IF
1054 87549132 : IF (ay > 0 .AND. by > 0) THEN
1055 : value = value + 0.5_dp*REAL(ay*by, dp)* &
1056 : cubeaxis(1, ax + sx, bx)* &
1057 7165620 : cubeaxis(2, ay - 1 + sy, by - 1)*cubeaxis(3, az + sz, bz)
1058 : END IF
1059 :
1060 : value = value + 2.0_dp*zeta*zetb* &
1061 : cubeaxis(1, ax + sx, bx)*cubeaxis(2, ay + sy, by)* &
1062 87549132 : cubeaxis(3, az + 1 + sz, bz + 1)
1063 87549132 : IF (az > 0) THEN
1064 : value = value - REAL(az, dp)*zetb* &
1065 : cubeaxis(1, ax + sx, bx)*cubeaxis(2, ay + sy, by)* &
1066 24120108 : cubeaxis(3, az - 1 + sz, bz + 1)
1067 : END IF
1068 87549132 : IF (bz > 0) THEN
1069 : value = value - zeta*REAL(bz, dp)* &
1070 : cubeaxis(1, ax + sx, bx)*cubeaxis(2, ay + sy, by)* &
1071 24120108 : cubeaxis(3, az + 1 + sz, bz - 1)
1072 : END IF
1073 87549132 : IF (az > 0 .AND. bz > 0) THEN
1074 : value = value + 0.5_dp*REAL(az*bz, dp)* &
1075 : cubeaxis(1, ax + sx, bx)*cubeaxis(2, ay + sy, by)* &
1076 7165620 : cubeaxis(3, az - 1 + sz, bz - 1)
1077 : END IF
1078 87549132 : END FUNCTION tau_plane_product_shifted
1079 :
1080 : END SUBROUTINE collocate_pgf_product_gspace
1081 :
1082 : ! **************************************************************************************************
1083 : !> \brief Buffered reference retained while validating the streaming implementation.
1084 : !> \param la_max ...
1085 : !> \param zeta ...
1086 : !> \param la_min ...
1087 : !> \param lb_max ...
1088 : !> \param zetb ...
1089 : !> \param lb_min ...
1090 : !> \param ra ...
1091 : !> \param rab ...
1092 : !> \param rab2 ...
1093 : !> \param scale ...
1094 : !> \param pab ...
1095 : !> \param na ...
1096 : !> \param nb ...
1097 : !> \param eps_rho_gspace ...
1098 : !> \param gsq_max ...
1099 : !> \param pw ...
1100 : !> \param compute_tau ...
1101 : !> \param potential_gspace ...
1102 : !> \param hab ...
1103 : !> \param center_force ...
1104 : !> \param strain_virial ...
1105 : ! **************************************************************************************************
1106 0 : SUBROUTINE collocate_pgf_product_gspace_buffered(la_max, zeta, la_min, &
1107 : lb_max, zetb, lb_min, &
1108 : ra, rab, rab2, scale, pab, na, nb, &
1109 : eps_rho_gspace, gsq_max, pw, compute_tau, &
1110 0 : potential_gspace, hab, center_force, strain_virial)
1111 :
1112 : ! NOTE: this routine is much slower than the real-space version of collocate_pgf_product
1113 :
1114 : INTEGER, INTENT(IN) :: la_max
1115 : REAL(dp), INTENT(IN) :: zeta
1116 : INTEGER, INTENT(IN) :: la_min, lb_max
1117 : REAL(dp), INTENT(IN) :: zetb
1118 : INTEGER, INTENT(IN) :: lb_min
1119 : REAL(dp), DIMENSION(3), INTENT(IN) :: ra, rab
1120 : REAL(dp), INTENT(IN) :: rab2, scale
1121 : REAL(dp), DIMENSION(:, :), POINTER :: pab
1122 : INTEGER, INTENT(IN) :: na, nb
1123 : REAL(dp), INTENT(IN) :: eps_rho_gspace, gsq_max
1124 : TYPE(pw_c1d_gs_type), INTENT(IN) :: pw
1125 : LOGICAL, INTENT(IN), OPTIONAL :: compute_tau
1126 : TYPE(pw_c1d_gs_type), INTENT(IN), OPTIONAL :: potential_gspace
1127 : REAL(dp), DIMENSION(:, :), INTENT(INOUT), OPTIONAL :: hab
1128 : REAL(dp), DIMENSION(3), INTENT(INOUT), OPTIONAL :: center_force
1129 : REAL(dp), DIMENSION(3, 3), INTENT(INOUT), OPTIONAL :: strain_virial
1130 :
1131 : CHARACTER(LEN=*), PARAMETER :: routineN = 'collocate_pgf_product_gspace'
1132 :
1133 : COMPLEX(dp) :: plane_product, rag, rbg, shifted_product, &
1134 : tau_product
1135 : COMPLEX(dp), DIMENSION(3) :: plane_gradient
1136 0 : COMPLEX(dp), DIMENSION(:, :, :, :), POINTER :: cubeaxis
1137 : INTEGER :: ax, ay, az, bx, by, bz, handle, i, ico, &
1138 : idir, jco, jdir, la, la_rec_max, lb, &
1139 : lb_grid, lb_rec_max, ub_grid
1140 : LOGICAL :: integrate_potential, integrate_strain, my_compute_tau
1141 : REAL(dp) :: f, fa, fb, g0_contribution, gcoord, pij, &
1142 : potential_product, prefactor, rzetp, &
1143 : twozetp, zetp
1144 : REAL(dp), DIMENSION(3) :: force_product, gvec
1145 : REAL(dp), DIMENSION(3) :: fap, fbp, rap, rbp, rp
1146 : REAL(dp), DIMENSION(3, 3) :: strain_product
1147 :
1148 0 : CALL timeset(routineN, handle)
1149 :
1150 0 : my_compute_tau = .FALSE.
1151 0 : IF (PRESENT(compute_tau)) my_compute_tau = compute_tau
1152 0 : integrate_potential = PRESENT(potential_gspace) .OR. PRESENT(hab)
1153 0 : integrate_strain = PRESENT(strain_virial)
1154 0 : CPASSERT(PRESENT(potential_gspace) .EQV. PRESENT(hab))
1155 0 : CPASSERT(.NOT. integrate_strain .OR. integrate_potential)
1156 0 : IF (integrate_potential) THEN
1157 0 : CPASSERT(ASSOCIATED(pw%pw_grid, potential_gspace%pw_grid))
1158 : END IF
1159 0 : la_rec_max = la_max
1160 0 : lb_rec_max = lb_max
1161 0 : IF (my_compute_tau) THEN
1162 0 : la_rec_max = la_rec_max + 1
1163 0 : lb_rec_max = lb_rec_max + 1
1164 : END IF
1165 0 : IF (integrate_strain) la_rec_max = la_rec_max + 1
1166 :
1167 0 : zetp = zeta + zetb
1168 0 : rzetp = 1.0_dp/zetp
1169 0 : f = zetb*rzetp
1170 0 : rap(:) = f*rab(:)
1171 0 : rbp(:) = rap(:) - rab(:)
1172 0 : rp(:) = ra(:) + rap(:)
1173 0 : twozetp = 2.0_dp*zetp
1174 0 : fap(:) = twozetp*rap(:)
1175 0 : fbp(:) = twozetp*rbp(:)
1176 :
1177 0 : prefactor = scale*SQRT((pi*rzetp)**3)*EXP(-zeta*f*rab2)
1178 0 : lb_grid = LBOUND(pw%array, 1)
1179 0 : ub_grid = UBOUND(pw%array, 1)
1180 0 : NULLIFY (cubeaxis)
1181 0 : CALL reallocate(cubeaxis, lb_grid, ub_grid, 1, 3, 0, la_rec_max, 0, lb_rec_max)
1182 :
1183 0 : DO idir = 1, 3
1184 0 : DO i = lb_grid, ub_grid
1185 0 : gcoord = pw%pw_grid%g(idir, i)
1186 : cubeaxis(i, idir, 0, 0) = &
1187 0 : EXP(CMPLX(-0.25_dp*rzetp*gcoord*gcoord, -rp(idir)*gcoord, KIND=dp))
1188 : END DO
1189 :
1190 0 : IF (la_rec_max > 0) THEN
1191 0 : DO i = lb_grid, ub_grid
1192 0 : gcoord = pw%pw_grid%g(idir, i)
1193 0 : rag = CMPLX(fap(idir), -gcoord, KIND=dp)
1194 0 : cubeaxis(i, idir, 1, 0) = rag*cubeaxis(i, idir, 0, 0)
1195 : END DO
1196 0 : DO la = 2, la_rec_max
1197 0 : fa = REAL(la - 1, dp)*twozetp
1198 0 : DO i = lb_grid, ub_grid
1199 0 : gcoord = pw%pw_grid%g(idir, i)
1200 0 : rag = CMPLX(fap(idir), -gcoord, KIND=dp)
1201 : cubeaxis(i, idir, la, 0) = rag*cubeaxis(i, idir, la - 1, 0) + &
1202 0 : fa*cubeaxis(i, idir, la - 2, 0)
1203 : END DO
1204 : END DO
1205 0 : IF (lb_rec_max > 0) THEN
1206 0 : fa = twozetp
1207 0 : DO i = lb_grid, ub_grid
1208 0 : gcoord = pw%pw_grid%g(idir, i)
1209 0 : rbg = CMPLX(fbp(idir), -gcoord, KIND=dp)
1210 0 : cubeaxis(i, idir, 0, 1) = rbg*cubeaxis(i, idir, 0, 0)
1211 : cubeaxis(i, idir, 1, 1) = rbg*cubeaxis(i, idir, 1, 0) + &
1212 0 : fa*cubeaxis(i, idir, 0, 0)
1213 : END DO
1214 0 : DO lb = 2, lb_rec_max
1215 0 : fb = REAL(lb - 1, dp)*twozetp
1216 0 : DO i = lb_grid, ub_grid
1217 0 : gcoord = pw%pw_grid%g(idir, i)
1218 0 : rbg = CMPLX(fbp(idir), -gcoord, KIND=dp)
1219 : cubeaxis(i, idir, 0, lb) = rbg*cubeaxis(i, idir, 0, lb - 1) + &
1220 0 : fb*cubeaxis(i, idir, 0, lb - 2)
1221 : cubeaxis(i, idir, 1, lb) = rbg*cubeaxis(i, idir, 1, lb - 1) + &
1222 : fb*cubeaxis(i, idir, 1, lb - 2) + &
1223 0 : fa*cubeaxis(i, idir, 0, lb - 1)
1224 : END DO
1225 : END DO
1226 0 : DO la = 2, la_rec_max
1227 0 : fa = REAL(la, dp)*twozetp
1228 0 : DO i = lb_grid, ub_grid
1229 0 : gcoord = pw%pw_grid%g(idir, i)
1230 0 : rbg = CMPLX(fbp(idir), -gcoord, KIND=dp)
1231 : cubeaxis(i, idir, la, 1) = rbg*cubeaxis(i, idir, la, 0) + &
1232 0 : fa*cubeaxis(i, idir, la - 1, 0)
1233 : END DO
1234 0 : DO lb = 2, lb_rec_max
1235 0 : fb = REAL(lb - 1, dp)*twozetp
1236 0 : DO i = lb_grid, ub_grid
1237 0 : gcoord = pw%pw_grid%g(idir, i)
1238 0 : rbg = CMPLX(fbp(idir), -gcoord, KIND=dp)
1239 : cubeaxis(i, idir, la, lb) = rbg*cubeaxis(i, idir, la, lb - 1) + &
1240 : fb*cubeaxis(i, idir, la, lb - 2) + &
1241 0 : fa*cubeaxis(i, idir, la - 1, lb - 1)
1242 : END DO
1243 : END DO
1244 : END DO
1245 : END IF
1246 : ELSE
1247 0 : IF (lb_rec_max > 0) THEN
1248 0 : DO i = lb_grid, ub_grid
1249 0 : gcoord = pw%pw_grid%g(idir, i)
1250 0 : rbg = CMPLX(fbp(idir), -gcoord, KIND=dp)
1251 0 : cubeaxis(i, idir, 0, 1) = rbg*cubeaxis(i, idir, 0, 0)
1252 : END DO
1253 0 : DO lb = 2, lb_rec_max
1254 0 : fb = REAL(lb - 1, dp)*twozetp
1255 0 : DO i = lb_grid, ub_grid
1256 0 : gcoord = pw%pw_grid%g(idir, i)
1257 0 : rbg = CMPLX(fbp(idir), -gcoord, KIND=dp)
1258 : cubeaxis(i, idir, 0, lb) = rbg*cubeaxis(i, idir, 0, lb - 1) + &
1259 0 : fb*cubeaxis(i, idir, 0, lb - 2)
1260 : END DO
1261 : END DO
1262 : END IF
1263 : END IF
1264 : END DO
1265 :
1266 0 : DO la = 0, la_rec_max
1267 0 : DO lb = 0, lb_rec_max
1268 0 : IF (la + lb == 0) CYCLE
1269 0 : fa = (1.0_dp/twozetp)**(la + lb)
1270 0 : DO idir = 1, 3
1271 0 : DO i = lb_grid, ub_grid
1272 0 : cubeaxis(i, idir, la, lb) = fa*cubeaxis(i, idir, la, lb)
1273 : END DO
1274 : END DO
1275 : END DO
1276 : END DO
1277 :
1278 : ! Add the current primitive Gaussian function product to grid
1279 :
1280 0 : DO ico = ncoset(la_min - 1) + 1, ncoset(la_max)
1281 :
1282 0 : ax = indco(1, ico)
1283 0 : ay = indco(2, ico)
1284 0 : az = indco(3, ico)
1285 :
1286 0 : DO jco = ncoset(lb_min - 1) + 1, ncoset(lb_max)
1287 :
1288 0 : IF (integrate_potential) THEN
1289 0 : g0_contribution = 0.0_dp
1290 0 : potential_product = 0.0_dp
1291 0 : force_product = 0.0_dp
1292 0 : strain_product = 0.0_dp
1293 0 : IF (PRESENT(center_force) .OR. integrate_strain) THEN
1294 0 : pij = prefactor*pab(na + ico, nb + jco)
1295 : END IF
1296 : ELSE
1297 0 : pij = prefactor*pab(na + ico, nb + jco)
1298 0 : IF (ABS(pij) < eps_rho_gspace) CYCLE
1299 : END IF
1300 :
1301 0 : bx = indco(1, jco)
1302 0 : by = indco(2, jco)
1303 0 : bz = indco(3, jco)
1304 :
1305 0 : DO i = lb_grid, ub_grid
1306 0 : IF (pw%pw_grid%gsq(i) > gsq_max) CYCLE
1307 0 : IF (.NOT. my_compute_tau) THEN
1308 : plane_product = cubeaxis(i, 1, ax, bx)* &
1309 : cubeaxis(i, 2, ay, by)* &
1310 0 : cubeaxis(i, 3, az, bz)
1311 : ELSE
1312 : ! Equivalent to GRID_FUNC_DADB: 0.5 * grad(phi_a).grad(phi_b).
1313 : tau_product = 2.0_dp*zeta*zetb* &
1314 : cubeaxis(i, 1, ax + 1, bx + 1)* &
1315 0 : cubeaxis(i, 2, ay, by)*cubeaxis(i, 3, az, bz)
1316 0 : IF (ax > 0) THEN
1317 : tau_product = tau_product - REAL(ax, dp)*zetb* &
1318 : cubeaxis(i, 1, ax - 1, bx + 1)* &
1319 0 : cubeaxis(i, 2, ay, by)*cubeaxis(i, 3, az, bz)
1320 : END IF
1321 0 : IF (bx > 0) THEN
1322 : tau_product = tau_product - zeta*REAL(bx, dp)* &
1323 : cubeaxis(i, 1, ax + 1, bx - 1)* &
1324 0 : cubeaxis(i, 2, ay, by)*cubeaxis(i, 3, az, bz)
1325 : END IF
1326 0 : IF (ax > 0 .AND. bx > 0) THEN
1327 : tau_product = tau_product + 0.5_dp*REAL(ax*bx, dp)* &
1328 : cubeaxis(i, 1, ax - 1, bx - 1)* &
1329 0 : cubeaxis(i, 2, ay, by)*cubeaxis(i, 3, az, bz)
1330 : END IF
1331 :
1332 : tau_product = tau_product + 2.0_dp*zeta*zetb* &
1333 : cubeaxis(i, 1, ax, bx)* &
1334 0 : cubeaxis(i, 2, ay + 1, by + 1)*cubeaxis(i, 3, az, bz)
1335 0 : IF (ay > 0) THEN
1336 : tau_product = tau_product - REAL(ay, dp)*zetb* &
1337 : cubeaxis(i, 1, ax, bx)* &
1338 0 : cubeaxis(i, 2, ay - 1, by + 1)*cubeaxis(i, 3, az, bz)
1339 : END IF
1340 0 : IF (by > 0) THEN
1341 : tau_product = tau_product - zeta*REAL(by, dp)* &
1342 : cubeaxis(i, 1, ax, bx)* &
1343 0 : cubeaxis(i, 2, ay + 1, by - 1)*cubeaxis(i, 3, az, bz)
1344 : END IF
1345 0 : IF (ay > 0 .AND. by > 0) THEN
1346 : tau_product = tau_product + 0.5_dp*REAL(ay*by, dp)* &
1347 : cubeaxis(i, 1, ax, bx)* &
1348 0 : cubeaxis(i, 2, ay - 1, by - 1)*cubeaxis(i, 3, az, bz)
1349 : END IF
1350 :
1351 : tau_product = tau_product + 2.0_dp*zeta*zetb* &
1352 : cubeaxis(i, 1, ax, bx)*cubeaxis(i, 2, ay, by)* &
1353 0 : cubeaxis(i, 3, az + 1, bz + 1)
1354 0 : IF (az > 0) THEN
1355 : tau_product = tau_product - REAL(az, dp)*zetb* &
1356 : cubeaxis(i, 1, ax, bx)*cubeaxis(i, 2, ay, by)* &
1357 0 : cubeaxis(i, 3, az - 1, bz + 1)
1358 : END IF
1359 0 : IF (bz > 0) THEN
1360 : tau_product = tau_product - zeta*REAL(bz, dp)* &
1361 : cubeaxis(i, 1, ax, bx)*cubeaxis(i, 2, ay, by)* &
1362 0 : cubeaxis(i, 3, az + 1, bz - 1)
1363 : END IF
1364 0 : IF (az > 0 .AND. bz > 0) THEN
1365 : tau_product = tau_product + 0.5_dp*REAL(az*bz, dp)* &
1366 : cubeaxis(i, 1, ax, bx)*cubeaxis(i, 2, ay, by)* &
1367 0 : cubeaxis(i, 3, az - 1, bz - 1)
1368 : END IF
1369 : plane_product = tau_product
1370 : END IF
1371 0 : IF (integrate_strain) THEN
1372 0 : IF (.NOT. my_compute_tau) THEN
1373 : shifted_product = cubeaxis(i, 1, ax + 1, bx)* &
1374 0 : cubeaxis(i, 2, ay, by)*cubeaxis(i, 3, az, bz)
1375 : plane_gradient(1) = CMPLX(0.0_dp, -1.0_dp, KIND=dp)* &
1376 0 : (shifted_product + ra(1)*plane_product)
1377 : shifted_product = cubeaxis(i, 1, ax, bx)* &
1378 0 : cubeaxis(i, 2, ay + 1, by)*cubeaxis(i, 3, az, bz)
1379 : plane_gradient(2) = CMPLX(0.0_dp, -1.0_dp, KIND=dp)* &
1380 0 : (shifted_product + ra(2)*plane_product)
1381 : shifted_product = cubeaxis(i, 1, ax, bx)* &
1382 0 : cubeaxis(i, 2, ay, by)*cubeaxis(i, 3, az + 1, bz)
1383 : plane_gradient(3) = CMPLX(0.0_dp, -1.0_dp, KIND=dp)* &
1384 0 : (shifted_product + ra(3)*plane_product)
1385 : ELSE
1386 0 : DO idir = 1, 3
1387 0 : shifted_product = tau_plane_product_shifted(idir)
1388 : plane_gradient(idir) = CMPLX(0.0_dp, -1.0_dp, KIND=dp)* &
1389 0 : (shifted_product + ra(idir)*plane_product)
1390 : END DO
1391 : END IF
1392 : END IF
1393 0 : IF (integrate_potential) THEN
1394 : potential_product = potential_product + prefactor* &
1395 0 : REAL(CONJG(potential_gspace%array(i))*plane_product, KIND=dp)
1396 0 : IF (PRESENT(center_force) .OR. integrate_strain) THEN
1397 0 : gvec = pw%pw_grid%g(:, i)
1398 : force_product = force_product + pij* &
1399 : REAL(CONJG(potential_gspace%array(i))* &
1400 0 : CMPLX(0.0_dp, -gvec, KIND=dp)*plane_product, KIND=dp)
1401 : END IF
1402 0 : IF (integrate_strain) THEN
1403 0 : DO idir = 1, 3
1404 0 : DO jdir = 1, 3
1405 : strain_product(idir, jdir) = strain_product(idir, jdir) - &
1406 : pij*gvec(idir)*REAL(CONJG(potential_gspace%array(i))* &
1407 0 : plane_gradient(jdir), KIND=dp)
1408 : END DO
1409 : END DO
1410 : END IF
1411 0 : IF (pw%pw_grid%have_g0 .AND. i == 1) THEN
1412 : g0_contribution = prefactor* &
1413 0 : REAL(CONJG(potential_gspace%array(i))*plane_product, KIND=dp)
1414 : END IF
1415 : ELSE
1416 0 : pw%array(i) = pw%array(i) + pij*plane_product
1417 : END IF
1418 : END DO
1419 :
1420 0 : IF (integrate_potential) THEN
1421 0 : IF (pw%pw_grid%grid_span == HALFSPACE) THEN
1422 0 : potential_product = 2.0_dp*potential_product - g0_contribution
1423 0 : force_product = 2.0_dp*force_product
1424 0 : strain_product = 2.0_dp*strain_product
1425 : END IF
1426 0 : hab(na + ico, nb + jco) = hab(na + ico, nb + jco) + potential_product
1427 0 : IF (PRESENT(center_force)) center_force = center_force + force_product
1428 0 : IF (integrate_strain) THEN
1429 0 : DO idir = 1, 3
1430 0 : DO jdir = 1, 3
1431 : strain_product(idir, jdir) = strain_product(idir, jdir) + &
1432 0 : force_product(idir)*ra(jdir)
1433 : END DO
1434 : strain_product(idir, idir) = strain_product(idir, idir) - &
1435 0 : pab(na + ico, nb + jco)*potential_product
1436 : END DO
1437 0 : strain_virial = strain_virial + strain_product
1438 : END IF
1439 : END IF
1440 :
1441 : END DO
1442 :
1443 : END DO
1444 :
1445 0 : DEALLOCATE (cubeaxis)
1446 :
1447 0 : CALL timestop(handle)
1448 :
1449 : CONTAINS
1450 :
1451 : ! **************************************************************************************************
1452 : !> \brief ...
1453 : !> \param shift_axis ...
1454 : !> \return ...
1455 : ! **************************************************************************************************
1456 0 : FUNCTION tau_plane_product_shifted(shift_axis) RESULT(value)
1457 : INTEGER, INTENT(IN) :: shift_axis
1458 : COMPLEX(dp) :: value
1459 :
1460 : INTEGER :: sx, sy, sz
1461 :
1462 0 : sx = 0
1463 0 : sy = 0
1464 0 : sz = 0
1465 0 : IF (shift_axis == 1) sx = 1
1466 0 : IF (shift_axis == 2) sy = 1
1467 0 : IF (shift_axis == 3) sz = 1
1468 :
1469 : value = 2.0_dp*zeta*zetb* &
1470 : cubeaxis(i, 1, ax + 1 + sx, bx + 1)* &
1471 0 : cubeaxis(i, 2, ay + sy, by)*cubeaxis(i, 3, az + sz, bz)
1472 0 : IF (ax > 0) THEN
1473 : value = value - REAL(ax, dp)*zetb* &
1474 : cubeaxis(i, 1, ax - 1 + sx, bx + 1)* &
1475 0 : cubeaxis(i, 2, ay + sy, by)*cubeaxis(i, 3, az + sz, bz)
1476 : END IF
1477 0 : IF (bx > 0) THEN
1478 : value = value - zeta*REAL(bx, dp)* &
1479 : cubeaxis(i, 1, ax + 1 + sx, bx - 1)* &
1480 0 : cubeaxis(i, 2, ay + sy, by)*cubeaxis(i, 3, az + sz, bz)
1481 : END IF
1482 0 : IF (ax > 0 .AND. bx > 0) THEN
1483 : value = value + 0.5_dp*REAL(ax*bx, dp)* &
1484 : cubeaxis(i, 1, ax - 1 + sx, bx - 1)* &
1485 0 : cubeaxis(i, 2, ay + sy, by)*cubeaxis(i, 3, az + sz, bz)
1486 : END IF
1487 :
1488 : value = value + 2.0_dp*zeta*zetb* &
1489 : cubeaxis(i, 1, ax + sx, bx)* &
1490 0 : cubeaxis(i, 2, ay + 1 + sy, by + 1)*cubeaxis(i, 3, az + sz, bz)
1491 0 : IF (ay > 0) THEN
1492 : value = value - REAL(ay, dp)*zetb* &
1493 : cubeaxis(i, 1, ax + sx, bx)* &
1494 0 : cubeaxis(i, 2, ay - 1 + sy, by + 1)*cubeaxis(i, 3, az + sz, bz)
1495 : END IF
1496 0 : IF (by > 0) THEN
1497 : value = value - zeta*REAL(by, dp)* &
1498 : cubeaxis(i, 1, ax + sx, bx)* &
1499 0 : cubeaxis(i, 2, ay + 1 + sy, by - 1)*cubeaxis(i, 3, az + sz, bz)
1500 : END IF
1501 0 : IF (ay > 0 .AND. by > 0) THEN
1502 : value = value + 0.5_dp*REAL(ay*by, dp)* &
1503 : cubeaxis(i, 1, ax + sx, bx)* &
1504 0 : cubeaxis(i, 2, ay - 1 + sy, by - 1)*cubeaxis(i, 3, az + sz, bz)
1505 : END IF
1506 :
1507 : value = value + 2.0_dp*zeta*zetb* &
1508 : cubeaxis(i, 1, ax + sx, bx)*cubeaxis(i, 2, ay + sy, by)* &
1509 0 : cubeaxis(i, 3, az + 1 + sz, bz + 1)
1510 0 : IF (az > 0) THEN
1511 : value = value - REAL(az, dp)*zetb* &
1512 : cubeaxis(i, 1, ax + sx, bx)*cubeaxis(i, 2, ay + sy, by)* &
1513 0 : cubeaxis(i, 3, az - 1 + sz, bz + 1)
1514 : END IF
1515 0 : IF (bz > 0) THEN
1516 : value = value - zeta*REAL(bz, dp)* &
1517 : cubeaxis(i, 1, ax + sx, bx)*cubeaxis(i, 2, ay + sy, by)* &
1518 0 : cubeaxis(i, 3, az + 1 + sz, bz - 1)
1519 : END IF
1520 0 : IF (az > 0 .AND. bz > 0) THEN
1521 : value = value + 0.5_dp*REAL(az*bz, dp)* &
1522 : cubeaxis(i, 1, ax + sx, bx)*cubeaxis(i, 2, ay + sy, by)* &
1523 0 : cubeaxis(i, 3, az - 1 + sz, bz - 1)
1524 : END IF
1525 :
1526 0 : END FUNCTION tau_plane_product_shifted
1527 :
1528 : END SUBROUTINE collocate_pgf_product_gspace_buffered
1529 :
1530 : END MODULE gapw_gspace_reconstruction
|
