(git:98357aa)
Loading...
Searching...
No Matches
gapw_gspace_reconstruction.F
Go to the documentation of this file.
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! **************************************************************************************************
12MODULE 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
65CONTAINS
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 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 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 REAL(kind=dp), DIMENSION(:, :), POINTER :: delta_cpc, pab, work, zet
112 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 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
117 TYPE(pw_c1d_gs_type) :: rho_elec_gspace
118 TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, source_r, tau_r
119 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
120 TYPE(qs_rho_type), POINTER :: rho
121 TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: cpc_h, cpc_s
122 TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
123 TYPE(rho_atom_type), POINTER :: rho_atom
124
125 cpassert(ASSOCIATED(qs_env))
126 cpassert(ASSOCIATED(auxbas_pw_pool))
127
128 CALL timeset(routinen, handle)
129
130 NULLIFY (atom_list)
131 NULLIFY (atomic_kind_set)
132 NULLIFY (qs_kind_set)
133 NULLIFY (cell)
134 NULLIFY (cpc_h)
135 NULLIFY (cpc_s)
136 NULLIFY (delta_cpc)
137 NULLIFY (dft_control)
138 NULLIFY (lmax)
139 NULLIFY (lmin)
140 NULLIFY (npgf)
141 NULLIFY (basis_1c_set)
142 NULLIFY (pab)
143 NULLIFY (particle_set)
144 NULLIFY (rho, rho_r, source_r, tau_r)
145 NULLIFY (rho_atom)
146 NULLIFY (rho_atom_set)
147 NULLIFY (work)
148 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 rho_atom_set=rho_atom_set)
157
158 IF (PRESENT(rho_source)) THEN
159 rho => rho_source
160 ELSE
161 CALL get_qs_env(qs_env=qs_env, rho=rho)
162 END IF
163
164 CALL qs_rho_get(rho, rho_r=rho_r, tau_r=tau_r)
165 my_compute_tau = .false.
166 IF (PRESENT(compute_tau)) my_compute_tau = compute_tau
167 IF (my_compute_tau) THEN
168 cpassert(ASSOCIATED(tau_r))
169 source_r => tau_r
170 ELSE
171 source_r => rho_r
172 END IF
173 eps_rho_gspace = dft_control%qs_control%eps_rho_gspace
174 nkind = SIZE(atomic_kind_set)
175 nspin = dft_control%nspins
176
177 ! Load the soft contribution of the electronic density
178
179 CALL auxbas_pw_pool%create_pw(pw=rho_elec_gspace)
180
181 CALL pw_zero(rhotot_elec_gspace)
182
183 DO ispin = 1, nspin
184 CALL pw_zero(rho_elec_gspace)
185 CALL pw_transfer(source_r(ispin), rho_elec_gspace)
186 IF (PRESENT(fsign) .AND. (ispin == 2)) THEN
187 alpha = fsign
188 ELSE
189 alpha = 1.0_dp
190 END IF
191 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 CALL auxbas_pw_pool%give_back_pw(rho_elec_gspace)
198
199 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 orthorhombic=orthorhombic)
203 permit_nonorthorhombic = .false.
204 IF (PRESENT(allow_nonorthorhombic)) permit_nonorthorhombic = allow_nonorthorhombic
205 IF (.NOT. orthorhombic .AND. .NOT. permit_nonorthorhombic) THEN
206 CALL cp_abort(__location__, &
207 "The calculation of XRD spectra for non-orthorhombic cells is not implemented")
208 END IF
209
210 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 DO ikind = 1, nkind
218
219 CALL get_atomic_kind(atomic_kind_set(ikind), &
220 atom_list=atom_list, &
221 natom=natom)
222
223 CALL get_qs_kind(qs_kind_set(ikind), &
224 basis_set=basis_1c_set, &
225 basis_type="GAPW_1C", &
226 paw_atom=paw_atom)
227
228 IF (.NOT. paw_atom) cycle ! no PAW atom: nothing to do
229
230 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 zet=zet)
240
241 ncotot = maxco*nset
242 nsotot = maxso*nset
243 CALL reallocate(delta_cpc, 1, nsatbas, 1, nsatbas)
244 CALL reallocate(pab, 1, ncotot, 1, ncotot)
245 CALL reallocate(work, 1, maxso, 1, maxco)
246
247 DO iatom = 1, natom
248
249 atom = atom_list(iatom)
250 rho_atom => rho_atom_set(atom)
251
252 CALL get_rho_atom(rho_atom=rho_atom, &
253 cpc_h=cpc_h, &
254 cpc_s=cpc_s)
255
256 ra(:) = pbc(particle_set(atom)%r, cell)
257
258 delta_cpc = 0.0_dp
259
260 DO ispin = 1, nspin
261 IF (PRESENT(fsign) .AND. (ispin == 2)) THEN
262 alpha = fsign
263 ELSE
264 alpha = 1.0_dp
265 END IF
266 delta_cpc = delta_cpc + alpha*(cpc_h(ispin)%r_coef - cpc_s(ispin)%r_coef)
267 END DO
268
269 scale = 1.0_dp
270
271 DO iset = 1, nset
272 ico1_set = (iset - 1)*maxco + 1
273 iso1_set = (iset - 1)*maxso + 1
274 ncoa = ncoset(lmax(iset))
275 nsoa = nsoset(lmax(iset))
276 DO jset = 1, nset
277 jco1_set = (jset - 1)*maxco + 1
278 jso1_set = (jset - 1)*maxso + 1
279 ncob = ncoset(lmax(jset))
280 nsob = nsoset(lmax(jset))
281 DO ipgf = 1, npgf(iset)
282 ico1_pgf = ico1_set + (ipgf - 1)*ncoa
283 iso1_pgf = iso1_set + (ipgf - 1)*nsoa
284 DO jpgf = 1, npgf(jset)
285 jco1_pgf = jco1_set + (jpgf - 1)*ncob
286 jso1_pgf = jso1_set + (jpgf - 1)*nsob
287 ico = ico1_pgf + ncoset(lmin(iset) - 1)
288 iso = iso1_pgf + nsoset(lmin(iset) - 1)
289
290 ! Transformation spherical to Cartesian
291
292 DO la = lmin(iset), lmax(iset)
293 jco = jco1_pgf + ncoset(lmin(jset) - 1)
294 jso = jso1_pgf + nsoset(lmin(jset) - 1)
295 DO lb = lmin(jset), lmax(jset)
296 ison = o2nindex(iso)
297 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 maxso)
302 CALL dgemm("T", "N", nco(la), nco(lb), nso(la), 1.0_dp, &
303 orbtramat(la)%slm, nso(la), work, maxso, &
304 0.0_dp, pab(ico:ico + nco(la) - 1, jco), SIZE(pab, 1))
305 jco = jco + nco(lb)
306 jso = jso + nso(lb)
307 END DO ! next lb
308 ico = ico + nco(la)
309 iso = iso + nso(la)
310 END DO ! next la
311
312 ! Collocate current product of primitive Cartesian functions
313
314 na = ico1_pgf - 1
315 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 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 DEALLOCATE (o2nindex)
342 END DO ! next atomic kind "ikind"
343
344 rho_total = pw_integrate_function(rhotot_elec_gspace, isign=-1)/volume
345
346 rho_hard = rho_total - rho_soft
347
348 ! Release work storage
349
350 IF (ASSOCIATED(delta_cpc)) THEN
351 DEALLOCATE (delta_cpc)
352 END IF
353
354 IF (ASSOCIATED(work)) THEN
355 DEALLOCATE (work)
356 END IF
357
358 IF (ASSOCIATED(pab)) THEN
359 DEALLOCATE (pab)
360 END IF
361
362 CALL timestop(handle)
363
364 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 SUBROUTINE integrate_gapw_composite_vxc_gspace(qs_env, auxbas_pw_pool, &
379 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 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 REAL(kind=dp), DIMENSION(:, :), POINTER :: delta_cpc, hab, hab_rho, hab_tau, pab, &
406 v_spherical, work, zet
407 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 TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
414 TYPE(pw_c1d_gs_type), ALLOCATABLE, DIMENSION(:) :: vrho_g, vtau_g
415 TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
416 TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: cpc_h, cpc_s, int_h, int_s
417 TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
418 TYPE(rho_atom_type), POINTER :: rho_atom
419
420 NULLIFY (logger)
421 logger => cp_get_default_logger()
422 output_unit = cp_logger_get_default_io_unit(logger)
423
424 CALL timeset(routinen, handle)
425
426 NULLIFY (atom_list, atomic_kind_set, basis_1c_set, cell, cpc_h, cpc_s, delta_cpc, &
427 dft_control, hab, hab_rho, hab_tau, int_h, int_s, lmax, lmin, npgf, o2nindex, &
428 pab, para_env, &
429 particle_set, qs_kind_set, rho_atom, rho_atom_set, work, zet)
430 NULLIFY (v_spherical)
431
432 cpassert(ASSOCIATED(qs_env))
433 cpassert(ASSOCIATED(auxbas_pw_pool))
434 cpassert(SIZE(vxc_rho) == SIZE(vxc_tau))
435 IF (PRESENT(one_center_contraction)) one_center_contraction = 0.0_dp
436 IF (PRESENT(one_center_rho_contraction)) one_center_rho_contraction = 0.0_dp
437 IF (PRESENT(one_center_tau_contraction)) one_center_tau_contraction = 0.0_dp
438 IF (PRESENT(atom_force)) THEN
439 cpassert(SIZE(atom_force, 1) == 3)
440 END IF
441 IF (PRESENT(strain_virial)) THEN
442 cpassert(SIZE(strain_virial, 1) == 3)
443 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 rho_atom_set=rho_atom_set)
454
455 eps_rho_gspace = dft_control%qs_control%eps_rho_gspace
456 nkind = SIZE(atomic_kind_set)
457 nspin = SIZE(vxc_rho)
458
459 ALLOCATE (vrho_g(nspin), vtau_g(nspin))
460 DO ispin = 1, nspin
461 CALL auxbas_pw_pool%create_pw(vrho_g(ispin))
462 CALL auxbas_pw_pool%create_pw(vtau_g(ispin))
463 CALL pw_transfer(vxc_rho(ispin), vrho_g(ispin))
464 CALL pw_transfer(vxc_tau(ispin), vtau_g(ispin))
465 END DO
466 q_max = sqrt(maxval(vrho_g(1)%pw_grid%gsq))
467
468 DO ikind = 1, nkind
469 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 basis_type="GAPW_1C", paw_atom=paw_atom)
472 IF (.NOT. paw_atom) cycle
473
474 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 npgf=npgf, nset=nset, zet=zet)
478 ncotot = maxco*nset
479 nsotot = maxso*nset
480 CALL reallocate(hab, 1, ncotot, 1, ncotot)
481 CALL reallocate(hab_rho, 1, ncotot, 1, ncotot)
482 CALL reallocate(hab_tau, 1, ncotot, 1, ncotot)
483 CALL reallocate(pab, 1, ncotot, 1, ncotot)
484 CALL reallocate(delta_cpc, 1, nsatbas, 1, nsatbas)
485 CALL reallocate(v_spherical, 1, nsotot, 1, nsotot)
486 CALL reallocate(work, 1, maxso, 1, maxco)
487 pab = 0.0_dp
488
489 DO iatom = 1, natom
490 atom = atom_list(iatom)
491 rho_atom => rho_atom_set(atom)
492 CALL get_rho_atom(rho_atom=rho_atom, cpc_h=cpc_h, cpc_s=cpc_s, &
493 ga_vlocal_gb_h=int_h, ga_vlocal_gb_s=int_s)
494 ra(:) = pbc(particle_set(atom)%r, cell)
495
496 DO ispin = 1, nspin
497 hab = 0.0_dp
498 hab_rho = 0.0_dp
499 hab_tau = 0.0_dp
500 pab = 0.0_dp
501 center_force = 0.0_dp
502 atom_strain_virial = 0.0_dp
503 delta_cpc = cpc_h(ispin)%r_coef - cpc_s(ispin)%r_coef
504 DO iset = 1, nset
505 ico1_set = (iset - 1)*maxco + 1
506 iso1_set = (iset - 1)*maxso + 1
507 ncoa = ncoset(lmax(iset))
508 nsoa = nsoset(lmax(iset))
509 DO jset = 1, nset
510 jco1_set = (jset - 1)*maxco + 1
511 jso1_set = (jset - 1)*maxso + 1
512 ncob = ncoset(lmax(jset))
513 nsob = nsoset(lmax(jset))
514 DO ipgf = 1, npgf(iset)
515 ico1_pgf = ico1_set + (ipgf - 1)*ncoa
516 iso1_pgf = iso1_set + (ipgf - 1)*nsoa
517 DO jpgf = 1, npgf(jset)
518 jco1_pgf = jco1_set + (jpgf - 1)*ncob
519 jso1_pgf = jso1_set + (jpgf - 1)*nsob
520 ico = ico1_pgf + ncoset(lmin(iset) - 1)
521 iso = iso1_pgf + nsoset(lmin(iset) - 1)
522 DO la = lmin(iset), lmax(iset)
523 jco = jco1_pgf + ncoset(lmin(jset) - 1)
524 jso = jso1_pgf + nsoset(lmin(jset) - 1)
525 DO lb = lmin(jset), lmax(jset)
526 ison = o2nindex(iso)
527 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 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 pab(ico:ico + nco(la) - 1, jco), SIZE(pab, 1))
535 jco = jco + nco(lb)
536 jso = jso + nso(lb)
537 END DO
538 ico = ico + nco(la)
539 iso = iso + nso(la)
540 END DO
541 na = ico1_pgf - 1
542 nb = jco1_pgf - 1
543 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 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 center_force=center_force, strain_virial=atom_strain_virial)
560 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 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 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 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 potential_gspace=vtau_g(ispin), hab=hab_tau)
592 END IF
593 END DO
594 END DO
595 END DO
596 END DO
597 CALL para_env%sum(hab_rho)
598 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 hab = hab_rho + hab_tau
602 adjoint_rho_cart = sum(pab*hab_rho)
603 adjoint_tau_cart = sum(pab*hab_tau)
604 adjoint_cart = adjoint_rho_cart + adjoint_tau_cart
605 IF (PRESENT(one_center_contraction)) THEN
606 one_center_contraction = one_center_contraction + adjoint_cart
607 END IF
608 IF (PRESENT(one_center_rho_contraction)) THEN
609 one_center_rho_contraction = one_center_rho_contraction + adjoint_rho_cart
610 END IF
611 IF (PRESENT(one_center_tau_contraction)) THEN
612 one_center_tau_contraction = one_center_tau_contraction + adjoint_tau_cart
613 END IF
614
615 IF (ASSOCIATED(int_h) .AND. ASSOCIATED(int_s)) THEN
616 IF (ASSOCIATED(int_h(ispin)%r_coef) .AND. ASSOCIATED(int_s(ispin)%r_coef)) THEN
617 v_spherical = 0.0_dp
618 DO iset = 1, nset
619 ico1_set = (iset - 1)*maxco + 1
620 iso1_set = (iset - 1)*maxso + 1
621 ncoa = ncoset(lmax(iset))
622 nsoa = nsoset(lmax(iset))
623 DO jset = 1, nset
624 jco1_set = (jset - 1)*maxco + 1
625 jso1_set = (jset - 1)*maxso + 1
626 ncob = ncoset(lmax(jset))
627 nsob = nsoset(lmax(jset))
628 DO ipgf = 1, npgf(iset)
629 ico1_pgf = ico1_set + (ipgf - 1)*ncoa
630 iso1_pgf = iso1_set + (ipgf - 1)*nsoa
631 DO jpgf = 1, npgf(jset)
632 jco1_pgf = jco1_set + (jpgf - 1)*ncob
633 jso1_pgf = jso1_set + (jpgf - 1)*nsob
634 ico = ico1_pgf + ncoset(lmin(iset) - 1)
635 iso = iso1_pgf + nsoset(lmin(iset) - 1)
636 DO la = lmin(iset), lmax(iset)
637 jco = jco1_pgf + ncoset(lmin(jset) - 1)
638 jso = jso1_pgf + nsoset(lmin(jset) - 1)
639 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 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 v_spherical(iso:iso + nso(la) - 1, jso), nsotot)
647 jco = jco + nco(lb)
648 jso = jso + nso(lb)
649 END DO
650 ico = ico + nco(la)
651 iso = iso + nso(la)
652 END DO
653 END DO
654 END DO
655 END DO
656 END DO
657 int_h(ispin)%r_coef = int_h(ispin)%r_coef + v_spherical
658 int_s(ispin)%r_coef = int_s(ispin)%r_coef + v_spherical
659 adjoint_spherical = 0.0_dp
660 DO jso = 1, nsotot
661 json = o2nindex(jso)
662 IF (json == 0) cycle
663 DO iso = 1, nsotot
664 ison = o2nindex(iso)
665 IF (ison == 0) cycle
666 adjoint_spherical = adjoint_spherical + &
667 delta_cpc(ison, json)*v_spherical(iso, jso)
668 END DO
669 END DO
670 adjoint_consistent = abs(adjoint_cart - adjoint_spherical) <= &
671 1.0e-9_dp*max(1.0_dp, abs(adjoint_cart))
672 IF (.NOT. adjoint_consistent) THEN
673 WRITE (unit=output_unit, fmt="(A,3I6,3ES24.15)") &
674 " SKALA_GPW| Composite adjoint mismatch", ikind, atom, ispin, &
675 adjoint_cart, adjoint_spherical, adjoint_cart - adjoint_spherical
676 END IF
677 cpassert(adjoint_consistent)
678 END IF
679 END IF
680 IF (PRESENT(atom_force)) atom_force(:, atom) = atom_force(:, atom) + center_force
681 IF (PRESENT(strain_virial)) strain_virial = strain_virial + atom_strain_virial
682 END DO
683 END DO
684
685 DEALLOCATE (delta_cpc, hab, hab_rho, hab_tau, pab, v_spherical, work, o2nindex)
686 END DO
687
688 DO ispin = 1, nspin
689 CALL auxbas_pw_pool%give_back_pw(vrho_g(ispin))
690 CALL auxbas_pw_pool%give_back_pw(vtau_g(ispin))
691 END DO
692 DEALLOCATE (vrho_g, vtau_g)
693
694 CALL timestop(handle)
695
696 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 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 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 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 CALL timeset(routinen, handle)
762
763 my_compute_tau = .false.
764 IF (PRESENT(compute_tau)) my_compute_tau = compute_tau
765 integrate_potential = PRESENT(potential_gspace) .OR. PRESENT(hab)
766 integrate_strain = PRESENT(strain_virial)
767 cpassert(PRESENT(potential_gspace) .EQV. PRESENT(hab))
768 cpassert(.NOT. integrate_strain .OR. integrate_potential)
769 IF (integrate_potential) THEN
770 cpassert(ASSOCIATED(pw%pw_grid, potential_gspace%pw_grid))
771 END IF
772
773 la_rec_max = la_max
774 lb_rec_max = lb_max
775 IF (my_compute_tau) THEN
776 la_rec_max = la_rec_max + 1
777 lb_rec_max = lb_rec_max + 1
778 END IF
779 IF (integrate_strain) la_rec_max = la_rec_max + 1
780
781 zetp = zeta + zetb
782 rzetp = 1.0_dp/zetp
783 f = zetb*rzetp
784 rap(:) = f*rab(:)
785 rbp(:) = rap(:) - rab(:)
786 rp(:) = ra(:) + rap(:)
787 twozetp = 2.0_dp*zetp
788 fap(:) = twozetp*rap(:)
789 fbp(:) = twozetp*rbp(:)
790 prefactor = scale*sqrt((pi*rzetp)**3)*exp(-zeta*f*rab2)
791
792 ALLOCATE (cubeaxis(3, 0:la_rec_max, 0:lb_rec_max))
793 lb_grid = lbound(pw%array, 1)
794 ub_grid = ubound(pw%array, 1)
795
796 DO i = lb_grid, ub_grid
797 IF (pw%pw_grid%gsq(i) > gsq_max) cycle
798 gvec = pw%pw_grid%g(:, i)
799 CALL build_cubeaxis(gvec)
800 reciprocal_weight = 1.0_dp
801 IF (pw%pw_grid%grid_span == halfspace .AND. pw%pw_grid%gsq(i) > 0.0_dp) THEN
802 reciprocal_weight = 2.0_dp
803 END IF
804
805 DO ico = ncoset(la_min - 1) + 1, ncoset(la_max)
806 ax = indco(1, ico)
807 ay = indco(2, ico)
808 az = indco(3, ico)
809
810 DO jco = ncoset(lb_min - 1) + 1, ncoset(lb_max)
811 bx = indco(1, jco)
812 by = indco(2, jco)
813 bz = indco(3, jco)
814 pij = prefactor*pab(na + ico, nb + jco)
815 IF (.NOT. integrate_potential .AND. abs(pij) < eps_rho_gspace) cycle
816
817 IF (.NOT. my_compute_tau) THEN
818 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 cubeaxis(2, ay, by)*cubeaxis(3, az, bz)
824 IF (ax > 0) THEN
825 tau_product = tau_product - real(ax, dp)*zetb* &
826 cubeaxis(1, ax - 1, bx + 1)* &
827 cubeaxis(2, ay, by)*cubeaxis(3, az, bz)
828 END IF
829 IF (bx > 0) THEN
830 tau_product = tau_product - zeta*real(bx, dp)* &
831 cubeaxis(1, ax + 1, bx - 1)* &
832 cubeaxis(2, ay, by)*cubeaxis(3, az, bz)
833 END IF
834 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 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 cubeaxis(2, ay + 1, by + 1)*cubeaxis(3, az, bz)
843 IF (ay > 0) THEN
844 tau_product = tau_product - real(ay, dp)*zetb* &
845 cubeaxis(1, ax, bx)* &
846 cubeaxis(2, ay - 1, by + 1)*cubeaxis(3, az, bz)
847 END IF
848 IF (by > 0) THEN
849 tau_product = tau_product - zeta*real(by, dp)* &
850 cubeaxis(1, ax, bx)* &
851 cubeaxis(2, ay + 1, by - 1)*cubeaxis(3, az, bz)
852 END IF
853 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 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 cubeaxis(3, az + 1, bz + 1)
862 IF (az > 0) THEN
863 tau_product = tau_product - real(az, dp)*zetb* &
864 cubeaxis(1, ax, bx)*cubeaxis(2, ay, by)* &
865 cubeaxis(3, az - 1, bz + 1)
866 END IF
867 IF (bz > 0) THEN
868 tau_product = tau_product - zeta*real(bz, dp)* &
869 cubeaxis(1, ax, bx)*cubeaxis(2, ay, by)* &
870 cubeaxis(3, az + 1, bz - 1)
871 END IF
872 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 cubeaxis(3, az - 1, bz - 1)
876 END IF
877 plane_product = tau_product
878 END IF
879
880 IF (integrate_strain) THEN
881 IF (.NOT. my_compute_tau) THEN
882 shifted_product = cubeaxis(1, ax + 1, bx)* &
883 cubeaxis(2, ay, by)*cubeaxis(3, az, bz)
884 plane_gradient(1) = cmplx(0.0_dp, -1.0_dp, kind=dp)* &
885 (shifted_product + ra(1)*plane_product)
886 shifted_product = cubeaxis(1, ax, bx)* &
887 cubeaxis(2, ay + 1, by)*cubeaxis(3, az, bz)
888 plane_gradient(2) = cmplx(0.0_dp, -1.0_dp, kind=dp)* &
889 (shifted_product + ra(2)*plane_product)
890 shifted_product = cubeaxis(1, ax, bx)* &
891 cubeaxis(2, ay, by)*cubeaxis(3, az + 1, bz)
892 plane_gradient(3) = cmplx(0.0_dp, -1.0_dp, kind=dp)* &
893 (shifted_product + ra(3)*plane_product)
894 ELSE
895 DO idir = 1, 3
896 shifted_product = tau_plane_product_shifted(idir)
897 plane_gradient(idir) = cmplx(0.0_dp, -1.0_dp, kind=dp)* &
898 (shifted_product + ra(idir)*plane_product)
899 END DO
900 END IF
901 END IF
902
903 IF (integrate_potential) THEN
904 potential_term = reciprocal_weight*prefactor* &
905 REAL(conjg(potential_gspace%array(i))*plane_product, kind=dp)
906 hab(na + ico, nb + jco) = hab(na + ico, nb + jco) + potential_term
907 IF (PRESENT(center_force) .OR. integrate_strain) THEN
908 force_product = reciprocal_weight*pij* &
909 REAL(conjg(potential_gspace%array(i))* &
910 cmplx(0.0_dp, -gvec, kind=dp)*plane_product, kind=dp)
911 IF (PRESENT(center_force)) center_force = center_force + force_product
912 END IF
913 IF (integrate_strain) THEN
914 DO idir = 1, 3
915 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 force_product(idir)*ra(jdir)
920 END DO
921 strain_virial(idir, idir) = strain_virial(idir, idir) - &
922 pab(na + ico, nb + jco)*potential_term
923 END DO
924 END IF
925 ELSE
926 pw%array(i) = pw%array(i) + pij*plane_product
927 END IF
928 END DO
929 END DO
930 END DO
931
932 DEALLOCATE (cubeaxis)
933 CALL timestop(handle)
934
935 CONTAINS
936
937! **************************************************************************************************
938!> \brief ...
939!> \param gvector ...
940! **************************************************************************************************
941 SUBROUTINE build_cubeaxis(gvector)
942 REAL(dp), DIMENSION(3), INTENT(IN) :: gvector
943
944 cubeaxis = cmplx(0.0_dp, 0.0_dp, kind=dp)
945 DO idir = 1, 3
946 gcoord = gvector(idir)
947 cubeaxis(idir, 0, 0) = &
948 exp(cmplx(-0.25_dp*rzetp*gcoord*gcoord, -rp(idir)*gcoord, kind=dp))
949
950 IF (la_rec_max > 0) THEN
951 rag = cmplx(fap(idir), -gcoord, kind=dp)
952 cubeaxis(idir, 1, 0) = rag*cubeaxis(idir, 0, 0)
953 DO la = 2, la_rec_max
954 fa = real(la - 1, dp)*twozetp
955 cubeaxis(idir, la, 0) = rag*cubeaxis(idir, la - 1, 0) + &
956 fa*cubeaxis(idir, la - 2, 0)
957 END DO
958 IF (lb_rec_max > 0) THEN
959 rbg = cmplx(fbp(idir), -gcoord, kind=dp)
960 fa = twozetp
961 cubeaxis(idir, 0, 1) = rbg*cubeaxis(idir, 0, 0)
962 cubeaxis(idir, 1, 1) = rbg*cubeaxis(idir, 1, 0) + &
963 fa*cubeaxis(idir, 0, 0)
964 DO lb = 2, lb_rec_max
965 fb = real(lb - 1, dp)*twozetp
966 cubeaxis(idir, 0, lb) = rbg*cubeaxis(idir, 0, lb - 1) + &
967 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 fa*cubeaxis(idir, 0, lb - 1)
971 END DO
972 DO la = 2, la_rec_max
973 fa = real(la, dp)*twozetp
974 cubeaxis(idir, la, 1) = rbg*cubeaxis(idir, la, 0) + &
975 fa*cubeaxis(idir, la - 1, 0)
976 DO lb = 2, lb_rec_max
977 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 fa*cubeaxis(idir, la - 1, lb - 1)
981 END DO
982 END DO
983 END IF
984 ELSE IF (lb_rec_max > 0) THEN
985 rbg = cmplx(fbp(idir), -gcoord, kind=dp)
986 cubeaxis(idir, 0, 1) = rbg*cubeaxis(idir, 0, 0)
987 DO lb = 2, lb_rec_max
988 fb = real(lb - 1, dp)*twozetp
989 cubeaxis(idir, 0, lb) = rbg*cubeaxis(idir, 0, lb - 1) + &
990 fb*cubeaxis(idir, 0, lb - 2)
991 END DO
992 END IF
993
994 DO la = 0, la_rec_max
995 DO lb = 0, lb_rec_max
996 IF (la + lb == 0) cycle
997 fa = (1.0_dp/twozetp)**(la + lb)
998 cubeaxis(idir, la, lb) = fa*cubeaxis(idir, la, lb)
999 END DO
1000 END DO
1001 END DO
1002 END SUBROUTINE build_cubeaxis
1003
1004! **************************************************************************************************
1005!> \brief ...
1006!> \param shift_axis ...
1007!> \return ...
1008! **************************************************************************************************
1009 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 sx = 0
1016 sy = 0
1017 sz = 0
1018 IF (shift_axis == 1) sx = 1
1019 IF (shift_axis == 2) sy = 1
1020 IF (shift_axis == 3) sz = 1
1021
1022 value = 2.0_dp*zeta*zetb* &
1023 cubeaxis(1, ax + 1 + sx, bx + 1)* &
1024 cubeaxis(2, ay + sy, by)*cubeaxis(3, az + sz, bz)
1025 IF (ax > 0) THEN
1026 value = value - real(ax, dp)*zetb* &
1027 cubeaxis(1, ax - 1 + sx, bx + 1)* &
1028 cubeaxis(2, ay + sy, by)*cubeaxis(3, az + sz, bz)
1029 END IF
1030 IF (bx > 0) THEN
1031 value = value - zeta*real(bx, dp)* &
1032 cubeaxis(1, ax + 1 + sx, bx - 1)* &
1033 cubeaxis(2, ay + sy, by)*cubeaxis(3, az + sz, bz)
1034 END IF
1035 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 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 cubeaxis(2, ay + 1 + sy, by + 1)*cubeaxis(3, az + sz, bz)
1044 IF (ay > 0) THEN
1045 value = value - real(ay, dp)*zetb* &
1046 cubeaxis(1, ax + sx, bx)* &
1047 cubeaxis(2, ay - 1 + sy, by + 1)*cubeaxis(3, az + sz, bz)
1048 END IF
1049 IF (by > 0) THEN
1050 value = value - zeta*real(by, dp)* &
1051 cubeaxis(1, ax + sx, bx)* &
1052 cubeaxis(2, ay + 1 + sy, by - 1)*cubeaxis(3, az + sz, bz)
1053 END IF
1054 IF (ay > 0 .AND. by > 0) THEN
1055 value = value + 0.5_dp*real(ay*by, dp)* &
1056 cubeaxis(1, ax + sx, bx)* &
1057 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 cubeaxis(3, az + 1 + sz, bz + 1)
1063 IF (az > 0) THEN
1064 value = value - real(az, dp)*zetb* &
1065 cubeaxis(1, ax + sx, bx)*cubeaxis(2, ay + sy, by)* &
1066 cubeaxis(3, az - 1 + sz, bz + 1)
1067 END IF
1068 IF (bz > 0) THEN
1069 value = value - zeta*real(bz, dp)* &
1070 cubeaxis(1, ax + sx, bx)*cubeaxis(2, ay + sy, by)* &
1071 cubeaxis(3, az + 1 + sz, bz - 1)
1072 END IF
1073 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 cubeaxis(3, az - 1 + sz, bz - 1)
1077 END IF
1078 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 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 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 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 CALL timeset(routinen, handle)
1149
1150 my_compute_tau = .false.
1151 IF (PRESENT(compute_tau)) my_compute_tau = compute_tau
1152 integrate_potential = PRESENT(potential_gspace) .OR. PRESENT(hab)
1153 integrate_strain = PRESENT(strain_virial)
1154 cpassert(PRESENT(potential_gspace) .EQV. PRESENT(hab))
1155 cpassert(.NOT. integrate_strain .OR. integrate_potential)
1156 IF (integrate_potential) THEN
1157 cpassert(ASSOCIATED(pw%pw_grid, potential_gspace%pw_grid))
1158 END IF
1159 la_rec_max = la_max
1160 lb_rec_max = lb_max
1161 IF (my_compute_tau) THEN
1162 la_rec_max = la_rec_max + 1
1163 lb_rec_max = lb_rec_max + 1
1164 END IF
1165 IF (integrate_strain) la_rec_max = la_rec_max + 1
1166
1167 zetp = zeta + zetb
1168 rzetp = 1.0_dp/zetp
1169 f = zetb*rzetp
1170 rap(:) = f*rab(:)
1171 rbp(:) = rap(:) - rab(:)
1172 rp(:) = ra(:) + rap(:)
1173 twozetp = 2.0_dp*zetp
1174 fap(:) = twozetp*rap(:)
1175 fbp(:) = twozetp*rbp(:)
1176
1177 prefactor = scale*sqrt((pi*rzetp)**3)*exp(-zeta*f*rab2)
1178 lb_grid = lbound(pw%array, 1)
1179 ub_grid = ubound(pw%array, 1)
1180 NULLIFY (cubeaxis)
1181 CALL reallocate(cubeaxis, lb_grid, ub_grid, 1, 3, 0, la_rec_max, 0, lb_rec_max)
1182
1183 DO idir = 1, 3
1184 DO i = lb_grid, ub_grid
1185 gcoord = pw%pw_grid%g(idir, i)
1186 cubeaxis(i, idir, 0, 0) = &
1187 exp(cmplx(-0.25_dp*rzetp*gcoord*gcoord, -rp(idir)*gcoord, kind=dp))
1188 END DO
1189
1190 IF (la_rec_max > 0) THEN
1191 DO i = lb_grid, ub_grid
1192 gcoord = pw%pw_grid%g(idir, i)
1193 rag = cmplx(fap(idir), -gcoord, kind=dp)
1194 cubeaxis(i, idir, 1, 0) = rag*cubeaxis(i, idir, 0, 0)
1195 END DO
1196 DO la = 2, la_rec_max
1197 fa = real(la - 1, dp)*twozetp
1198 DO i = lb_grid, ub_grid
1199 gcoord = pw%pw_grid%g(idir, i)
1200 rag = cmplx(fap(idir), -gcoord, kind=dp)
1201 cubeaxis(i, idir, la, 0) = rag*cubeaxis(i, idir, la - 1, 0) + &
1202 fa*cubeaxis(i, idir, la - 2, 0)
1203 END DO
1204 END DO
1205 IF (lb_rec_max > 0) THEN
1206 fa = twozetp
1207 DO i = lb_grid, ub_grid
1208 gcoord = pw%pw_grid%g(idir, i)
1209 rbg = cmplx(fbp(idir), -gcoord, kind=dp)
1210 cubeaxis(i, idir, 0, 1) = rbg*cubeaxis(i, idir, 0, 0)
1211 cubeaxis(i, idir, 1, 1) = rbg*cubeaxis(i, idir, 1, 0) + &
1212 fa*cubeaxis(i, idir, 0, 0)
1213 END DO
1214 DO lb = 2, lb_rec_max
1215 fb = real(lb - 1, dp)*twozetp
1216 DO i = lb_grid, ub_grid
1217 gcoord = pw%pw_grid%g(idir, i)
1218 rbg = cmplx(fbp(idir), -gcoord, kind=dp)
1219 cubeaxis(i, idir, 0, lb) = rbg*cubeaxis(i, idir, 0, lb - 1) + &
1220 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 fa*cubeaxis(i, idir, 0, lb - 1)
1224 END DO
1225 END DO
1226 DO la = 2, la_rec_max
1227 fa = real(la, dp)*twozetp
1228 DO i = lb_grid, ub_grid
1229 gcoord = pw%pw_grid%g(idir, i)
1230 rbg = cmplx(fbp(idir), -gcoord, kind=dp)
1231 cubeaxis(i, idir, la, 1) = rbg*cubeaxis(i, idir, la, 0) + &
1232 fa*cubeaxis(i, idir, la - 1, 0)
1233 END DO
1234 DO lb = 2, lb_rec_max
1235 fb = real(lb - 1, dp)*twozetp
1236 DO i = lb_grid, ub_grid
1237 gcoord = pw%pw_grid%g(idir, i)
1238 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 fa*cubeaxis(i, idir, la - 1, lb - 1)
1242 END DO
1243 END DO
1244 END DO
1245 END IF
1246 ELSE
1247 IF (lb_rec_max > 0) THEN
1248 DO i = lb_grid, ub_grid
1249 gcoord = pw%pw_grid%g(idir, i)
1250 rbg = cmplx(fbp(idir), -gcoord, kind=dp)
1251 cubeaxis(i, idir, 0, 1) = rbg*cubeaxis(i, idir, 0, 0)
1252 END DO
1253 DO lb = 2, lb_rec_max
1254 fb = real(lb - 1, dp)*twozetp
1255 DO i = lb_grid, ub_grid
1256 gcoord = pw%pw_grid%g(idir, i)
1257 rbg = cmplx(fbp(idir), -gcoord, kind=dp)
1258 cubeaxis(i, idir, 0, lb) = rbg*cubeaxis(i, idir, 0, lb - 1) + &
1259 fb*cubeaxis(i, idir, 0, lb - 2)
1260 END DO
1261 END DO
1262 END IF
1263 END IF
1264 END DO
1265
1266 DO la = 0, la_rec_max
1267 DO lb = 0, lb_rec_max
1268 IF (la + lb == 0) cycle
1269 fa = (1.0_dp/twozetp)**(la + lb)
1270 DO idir = 1, 3
1271 DO i = lb_grid, ub_grid
1272 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 DO ico = ncoset(la_min - 1) + 1, ncoset(la_max)
1281
1282 ax = indco(1, ico)
1283 ay = indco(2, ico)
1284 az = indco(3, ico)
1285
1286 DO jco = ncoset(lb_min - 1) + 1, ncoset(lb_max)
1287
1288 IF (integrate_potential) THEN
1289 g0_contribution = 0.0_dp
1290 potential_product = 0.0_dp
1291 force_product = 0.0_dp
1292 strain_product = 0.0_dp
1293 IF (PRESENT(center_force) .OR. integrate_strain) THEN
1294 pij = prefactor*pab(na + ico, nb + jco)
1295 END IF
1296 ELSE
1297 pij = prefactor*pab(na + ico, nb + jco)
1298 IF (abs(pij) < eps_rho_gspace) cycle
1299 END IF
1300
1301 bx = indco(1, jco)
1302 by = indco(2, jco)
1303 bz = indco(3, jco)
1304
1305 DO i = lb_grid, ub_grid
1306 IF (pw%pw_grid%gsq(i) > gsq_max) cycle
1307 IF (.NOT. my_compute_tau) THEN
1308 plane_product = cubeaxis(i, 1, ax, bx)* &
1309 cubeaxis(i, 2, ay, by)* &
1310 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 cubeaxis(i, 2, ay, by)*cubeaxis(i, 3, az, bz)
1316 IF (ax > 0) THEN
1317 tau_product = tau_product - real(ax, dp)*zetb* &
1318 cubeaxis(i, 1, ax - 1, bx + 1)* &
1319 cubeaxis(i, 2, ay, by)*cubeaxis(i, 3, az, bz)
1320 END IF
1321 IF (bx > 0) THEN
1322 tau_product = tau_product - zeta*real(bx, dp)* &
1323 cubeaxis(i, 1, ax + 1, bx - 1)* &
1324 cubeaxis(i, 2, ay, by)*cubeaxis(i, 3, az, bz)
1325 END IF
1326 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 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 cubeaxis(i, 2, ay + 1, by + 1)*cubeaxis(i, 3, az, bz)
1335 IF (ay > 0) THEN
1336 tau_product = tau_product - real(ay, dp)*zetb* &
1337 cubeaxis(i, 1, ax, bx)* &
1338 cubeaxis(i, 2, ay - 1, by + 1)*cubeaxis(i, 3, az, bz)
1339 END IF
1340 IF (by > 0) THEN
1341 tau_product = tau_product - zeta*real(by, dp)* &
1342 cubeaxis(i, 1, ax, bx)* &
1343 cubeaxis(i, 2, ay + 1, by - 1)*cubeaxis(i, 3, az, bz)
1344 END IF
1345 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 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 cubeaxis(i, 3, az + 1, bz + 1)
1354 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 cubeaxis(i, 3, az - 1, bz + 1)
1358 END IF
1359 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 cubeaxis(i, 3, az + 1, bz - 1)
1363 END IF
1364 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 cubeaxis(i, 3, az - 1, bz - 1)
1368 END IF
1369 plane_product = tau_product
1370 END IF
1371 IF (integrate_strain) THEN
1372 IF (.NOT. my_compute_tau) THEN
1373 shifted_product = cubeaxis(i, 1, ax + 1, bx)* &
1374 cubeaxis(i, 2, ay, by)*cubeaxis(i, 3, az, bz)
1375 plane_gradient(1) = cmplx(0.0_dp, -1.0_dp, kind=dp)* &
1376 (shifted_product + ra(1)*plane_product)
1377 shifted_product = cubeaxis(i, 1, ax, bx)* &
1378 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 (shifted_product + ra(2)*plane_product)
1381 shifted_product = cubeaxis(i, 1, ax, bx)* &
1382 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 (shifted_product + ra(3)*plane_product)
1385 ELSE
1386 DO idir = 1, 3
1387 shifted_product = tau_plane_product_shifted(idir)
1388 plane_gradient(idir) = cmplx(0.0_dp, -1.0_dp, kind=dp)* &
1389 (shifted_product + ra(idir)*plane_product)
1390 END DO
1391 END IF
1392 END IF
1393 IF (integrate_potential) THEN
1394 potential_product = potential_product + prefactor* &
1395 REAL(conjg(potential_gspace%array(i))*plane_product, kind=dp)
1396 IF (PRESENT(center_force) .OR. integrate_strain) THEN
1397 gvec = pw%pw_grid%g(:, i)
1398 force_product = force_product + pij* &
1399 REAL(conjg(potential_gspace%array(i))* &
1400 cmplx(0.0_dp, -gvec, kind=dp)*plane_product, kind=dp)
1401 END IF
1402 IF (integrate_strain) THEN
1403 DO idir = 1, 3
1404 DO jdir = 1, 3
1405 strain_product(idir, jdir) = strain_product(idir, jdir) - &
1406 pij*gvec(idir)*real(conjg(potential_gspace%array(i))* &
1407 plane_gradient(jdir), kind=dp)
1408 END DO
1409 END DO
1410 END IF
1411 IF (pw%pw_grid%have_g0 .AND. i == 1) THEN
1412 g0_contribution = prefactor* &
1413 REAL(conjg(potential_gspace%array(i))*plane_product, kind=dp)
1414 END IF
1415 ELSE
1416 pw%array(i) = pw%array(i) + pij*plane_product
1417 END IF
1418 END DO
1419
1420 IF (integrate_potential) THEN
1421 IF (pw%pw_grid%grid_span == halfspace) THEN
1422 potential_product = 2.0_dp*potential_product - g0_contribution
1423 force_product = 2.0_dp*force_product
1424 strain_product = 2.0_dp*strain_product
1425 END IF
1426 hab(na + ico, nb + jco) = hab(na + ico, nb + jco) + potential_product
1427 IF (PRESENT(center_force)) center_force = center_force + force_product
1428 IF (integrate_strain) THEN
1429 DO idir = 1, 3
1430 DO jdir = 1, 3
1431 strain_product(idir, jdir) = strain_product(idir, jdir) + &
1432 force_product(idir)*ra(jdir)
1433 END DO
1434 strain_product(idir, idir) = strain_product(idir, idir) - &
1435 pab(na + ico, nb + jco)*potential_product
1436 END DO
1437 strain_virial = strain_virial + strain_product
1438 END IF
1439 END IF
1440
1441 END DO
1442
1443 END DO
1444
1445 DEALLOCATE (cubeaxis)
1446
1447 CALL timestop(handle)
1448
1449 CONTAINS
1450
1451! **************************************************************************************************
1452!> \brief ...
1453!> \param shift_axis ...
1454!> \return ...
1455! **************************************************************************************************
1456 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 sx = 0
1463 sy = 0
1464 sz = 0
1465 IF (shift_axis == 1) sx = 1
1466 IF (shift_axis == 2) sy = 1
1467 IF (shift_axis == 3) sz = 1
1468
1469 value = 2.0_dp*zeta*zetb* &
1470 cubeaxis(i, 1, ax + 1 + sx, bx + 1)* &
1471 cubeaxis(i, 2, ay + sy, by)*cubeaxis(i, 3, az + sz, bz)
1472 IF (ax > 0) THEN
1473 value = value - real(ax, dp)*zetb* &
1474 cubeaxis(i, 1, ax - 1 + sx, bx + 1)* &
1475 cubeaxis(i, 2, ay + sy, by)*cubeaxis(i, 3, az + sz, bz)
1476 END IF
1477 IF (bx > 0) THEN
1478 value = value - zeta*real(bx, dp)* &
1479 cubeaxis(i, 1, ax + 1 + sx, bx - 1)* &
1480 cubeaxis(i, 2, ay + sy, by)*cubeaxis(i, 3, az + sz, bz)
1481 END IF
1482 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 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 cubeaxis(i, 2, ay + 1 + sy, by + 1)*cubeaxis(i, 3, az + sz, bz)
1491 IF (ay > 0) THEN
1492 value = value - real(ay, dp)*zetb* &
1493 cubeaxis(i, 1, ax + sx, bx)* &
1494 cubeaxis(i, 2, ay - 1 + sy, by + 1)*cubeaxis(i, 3, az + sz, bz)
1495 END IF
1496 IF (by > 0) THEN
1497 value = value - zeta*real(by, dp)* &
1498 cubeaxis(i, 1, ax + sx, bx)* &
1499 cubeaxis(i, 2, ay + 1 + sy, by - 1)*cubeaxis(i, 3, az + sz, bz)
1500 END IF
1501 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 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 cubeaxis(i, 3, az + 1 + sz, bz + 1)
1510 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 cubeaxis(i, 3, az - 1 + sz, bz + 1)
1514 END IF
1515 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 cubeaxis(i, 3, az + 1 + sz, bz - 1)
1519 END IF
1520 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 cubeaxis(i, 3, az - 1 + sz, bz - 1)
1524 END IF
1525
1526 END FUNCTION tau_plane_product_shifted
1527
1528 END SUBROUTINE collocate_pgf_product_gspace_buffered
1529
1530END MODULE gapw_gspace_reconstruction
dgemm
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Definition grid_cpu_task_list.c:225
cell_types::pbc
Definition cell_types.F:104
memory_utilities::reallocate
Definition memory_utilities.F:27
pw_methods::pw_axpy
Definition pw_methods.F:165
pw_methods::pw_integrate_function
Definition pw_methods.F:116
pw_methods::pw_scale
Definition pw_methods.F:94
pw_methods::pw_transfer
Definition pw_methods.F:304
pw_methods::pw_zero
Definition pw_methods.F:83
atom
Definition atom.F:9
atomic_kind_types
Define the atomic kind types and their sub types.
Definition atomic_kind_types.F:23
atomic_kind_types::get_atomic_kind
subroutine, public get_atomic_kind(atomic_kind, fist_potential, element_symbol, name, mass, kind_number, natom, atom_list, rcov, rvdw, z, qeff, apol, cpol, mm_radius, shell, shell_active, damping)
Get attributes of an atomic kind.
Definition atomic_kind_types.F:138
basis_set_types
Definition basis_set_types.F:15
basis_set_types::get_gto_basis_set
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius, npgf_seg_sum, ccon)
...
Definition basis_set_types.F:643
cell_types
Handles all functions related to the CELL.
Definition cell_types.F:15
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_logger_get_default_io_unit
integer function, public cp_logger_get_default_io_unit(logger)
returns the unit nr for the ionode (-1 on all other processors) skips as well checks if the procs cal...
Definition cp_log_handling.F:523
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
gapw_gspace_reconstruction
GAPW reciprocal-space reconstruction and its discrete adjoint.
Definition gapw_gspace_reconstruction.F:12
gapw_gspace_reconstruction::integrate_gapw_composite_vxc_gspace
subroutine, public integrate_gapw_composite_vxc_gspace(qs_env, auxbas_pw_pool, vxc_rho, vxc_tau, atom_force, strain_virial, one_center_contraction, one_center_rho_contraction, one_center_tau_contraction)
Apply the discrete adjoint of the common-grid GAPW hard-minus-soft reconstruction.
Definition gapw_gspace_reconstruction.F:382
gapw_gspace_reconstruction::calculate_rhotot_elec_gspace
subroutine, public calculate_rhotot_elec_gspace(qs_env, auxbas_pw_pool, rhotot_elec_gspace, q_max, rho_hard, rho_soft, fsign, compute_tau, rho_source, allow_nonorthorhombic)
The total electronic density in reciprocal space (g-space) is calculated.
Definition gapw_gspace_reconstruction.F:88
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::pi
real(kind=dp), parameter, public pi
Definition mathconstants.F:110
memory_utilities
Utility routines for the memory handling.
Definition memory_utilities.F:14
message_passing
Interface to the message passing library MPI.
Definition message_passing.F:23
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::nco
integer, dimension(:), allocatable, public nco
Definition orbital_pointers.F:42
orbital_pointers::nsoset
integer, dimension(:), allocatable, public nsoset
Definition orbital_pointers.F:42
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition orbital_pointers.F:42
orbital_pointers::indco
integer, dimension(:, :), allocatable, public indco
Definition orbital_pointers.F:43
orbital_pointers::nso
integer, dimension(:), allocatable, public nso
Definition orbital_pointers.F:42
orbital_transformation_matrices
Calculation of the spherical harmonics and the corresponding orbital transformation matrices.
Definition orbital_transformation_matrices.F:17
orbital_transformation_matrices::orbtramat
type(orbtramat_type), dimension(:), pointer, public orbtramat
Definition orbital_transformation_matrices.F:47
particle_types
Define the data structure for the particle information.
Definition particle_types.F:19
paw_basis_types
Definition paw_basis_types.F:13
paw_basis_types::get_paw_basis_info
subroutine, public get_paw_basis_info(basis_1c, o2nindex, n2oindex, nsatbas)
Return some info on the PAW basis derived from a GTO basis set.
Definition paw_basis_types.F:40
pw_grid_types
Definition pw_grid_types.F:13
pw_grid_types::halfspace
integer, parameter, public halfspace
Definition pw_grid_types.F:26
pw_grids
This module defines the grid data type and some basic operations on it.
Definition pw_grids.F:36
pw_grids::get_pw_grid_info
subroutine, public get_pw_grid_info(pw_grid, id_nr, mode, vol, dvol, npts, ngpts, ngpts_cut, dr, cutoff, orthorhombic, gvectors, gsquare)
Access to information stored in the pw_grid_type.
Definition pw_grids.F:185
pw_methods
Definition pw_methods.F:25
pw_pool_types
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:24
pw_types
Definition pw_types.F:24
qs_environment_types
Definition qs_environment_types.F:14
qs_environment_types::get_qs_env
subroutine, public get_qs_env(qs_env, atomic_kind_set, qs_kind_set, cell, super_cell, cell_ref, use_ref_cell, kpoints, dft_control, mos, sab_orb, sab_all, qmmm, qmmm_periodic, mimic, sac_ae, sac_ppl, sac_lri, sap_ppnl, sab_vdw, sab_scp, sap_oce, sab_lrc, sab_se, sab_xtbe, sab_tbe, sab_core, sab_xb, sab_xtb_pp, sab_xtb_nonbond, sab_almo, sab_kp, sab_kp_nosym, sab_cneo, particle_set, energy, force, matrix_h, matrix_h_im, matrix_ks, matrix_ks_im, matrix_vxc, run_rtp, rtp, matrix_h_kp, matrix_h_im_kp, matrix_ks_kp, matrix_ks_im_kp, matrix_vxc_kp, kinetic_kp, matrix_s_kp, matrix_w_kp, matrix_s_ri_aux_kp, matrix_s, matrix_s_ri_aux, matrix_w, matrix_p_mp2, matrix_p_mp2_admm, matrix_vhxc, rho, rho_xc, pw_env, ewald_env, ewald_pw, active_space, mpools, input, para_env, blacs_env, scf_control, rel_control, kinetic, qs_charges, vppl, xcint_weights, rho_core, rho_nlcc, rho_nlcc_g, ks_env, ks_qmmm_env, wf_history, scf_env, local_particles, local_molecules, distribution_2d, dbcsr_dist, molecule_kind_set, molecule_set, subsys, cp_subsys, oce, local_rho_set, rho_atom_set, task_list, task_list_soft, rho0_atom_set, rho0_mpole, rhoz_set, rhoz_cneo_set, ecoul_1c, rho0_s_rs, rho0_s_gs, rhoz_cneo_s_rs, rhoz_cneo_s_gs, do_kpoints, has_unit_metric, requires_mo_derivs, mo_derivs, mo_loc_history, nkind, natom, nelectron_total, nelectron_spin, efield, neighbor_list_id, linres_control, xas_env, virial, cp_ddapc_env, cp_ddapc_ewald, outer_scf_history, outer_scf_ihistory, x_data, et_coupling, dftb_potential, results, se_taper, se_store_int_env, se_nddo_mpole, se_nonbond_env, admm_env, lri_env, lri_density, exstate_env, ec_env, harris_env, dispersion_env, gcp_env, vee, rho_external, external_vxc, mask, mp2_env, bs_env, kg_env, wanniercentres, atprop, ls_scf_env, do_transport, transport_env, v_hartree_rspace, s_mstruct_changed, rho_changed, potential_changed, forces_up_to_date, mscfg_env, almo_scf_env, gradient_history, variable_history, embed_pot, spin_embed_pot, polar_env, mos_last_converged, eeq, rhs, do_rixs, tb_tblite)
Get the QUICKSTEP environment.
Definition qs_environment_types.F:534
qs_kind_types
Define the quickstep kind type and their sub types.
Definition qs_kind_types.F:23
qs_kind_types::get_qs_kind
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, cneo_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, zatom, zeff, elec_conf, mao, lmax_dftb, alpha_core_charge, ccore_charge, core_charge, core_charge_radius, paw_proj_set, paw_atom, hard_radius, hard0_radius, max_rad_local, covalent_radius, vdw_radius, gpw_type_forced, harmonics, max_iso_not0, max_s_harm, grid_atom, ngrid_ang, ngrid_rad, lmax_rho0, dft_plus_u_atom, l_of_dft_plus_u, n_of_dft_plus_u, u_minus_j, hund_j, u_of_dft_plus_u, j_of_dft_plus_u, alpha_of_dft_plus_u, beta_of_dft_plus_u, j0_of_dft_plus_u, occupation_of_dft_plus_u, dispersion, bs_occupation, magnetization, no_optimize, addel, laddel, naddel, orbitals, max_scf, eps_scf, smear, u_ramping, u_minus_j_target, eps_u_ramping, proj_shell_charge, lr_atom, do_mtlr, u_j_loop, ao_coef, init_u_ramping_each_scf, reltmat, ghost, monovalent, floating, name, element_symbol, pao_basis_size, pao_model_file, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
Definition qs_kind_types.F:491
qs_rho_atom_types
Definition qs_rho_atom_types.F:8
qs_rho_atom_types::get_rho_atom
subroutine, public get_rho_atom(rho_atom, cpc_h, cpc_s, rho_rad_h, rho_rad_s, drho_rad_h, drho_rad_s, vrho_rad_h, vrho_rad_s, rho_rad_h_d, rho_rad_s_d, ga_vlocal_gb_h, ga_vlocal_gb_s, int_scr_h, int_scr_s)
...
Definition qs_rho_atom_types.F:195
qs_rho_types
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition qs_rho_types.F:18
qs_rho_types::qs_rho_get
subroutine, public qs_rho_get(rho_struct, rho_ao, rho_ao_im, rho_ao_kp, rho_ao_im_kp, rho_r, drho_r, rho_g, drho_g, tau_r, tau_g, rho_r_valid, drho_r_valid, rho_g_valid, drho_g_valid, tau_r_valid, tau_g_valid, tot_rho_r, tot_rho_g, rho_r_sccs, soft_valid, complex_rho_ao)
returns info about the density described by this object. If some representation is not available an e...
Definition qs_rho_types.F:229
atomic_kind_types::atomic_kind_type
Provides all information about an atomic kind.
Definition atomic_kind_types.F:45
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
cell_types::cell_type
Type defining parameters related to the simulation cell.
Definition cell_types.F:60
cp_control_types::dft_control_type
Definition cp_control_types.F:757
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
message_passing::mp_para_env_type
stores all the informations relevant to an mpi environment
Definition message_passing.F:750
particle_types::particle_type
Definition particle_types.F:35
pw_pool_types::pw_pool_type
Manages a pool of grids (to be used for example as tmp objects), but can also be used to instantiate ...
Definition pw_pool_types.F:83
pw_types::pw_c1d_gs_type
Definition pw_types.F:127
pw_types::pw_r3d_rs_type
Definition pw_types.F:62
qs_environment_types::qs_environment_type
Definition qs_environment_types.F:222
qs_kind_types::qs_kind_type
Provides all information about a quickstep kind.
Definition qs_kind_types.F:190
qs_rho_atom_types::rho_atom_coeff
Definition qs_rho_atom_types.F:19
qs_rho_atom_types::rho_atom_type
Definition qs_rho_atom_types.F:23
qs_rho_types::qs_rho_type
keeps the density in various representations, keeping track of which ones are valid.
Definition qs_rho_types.F:62