Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Calculates Nabla V_KS (local part if PSP) on the different grids
10 : !> \par History
11 : !> created 06-2007 [RD]
12 : !> \author RD
13 : ! **************************************************************************************************
14 : MODULE qs_linres_epr_nablavks
15 : USE atomic_kind_types, ONLY: atomic_kind_type,&
16 : get_atomic_kind
17 : USE cell_types, ONLY: cell_type
18 : USE cp_control_types, ONLY: dft_control_type
19 : USE cp_dbcsr_api, ONLY: dbcsr_p_type
20 : USE cp_log_handling, ONLY: cp_get_default_logger,&
21 : cp_logger_type
22 : USE cp_output_handling, ONLY: cp_p_file,&
23 : cp_print_key_finished_output,&
24 : cp_print_key_should_output,&
25 : cp_print_key_unit_nr
26 : USE cp_realspace_grid_cube, ONLY: cp_pw_to_cube
27 : USE external_potential_types, ONLY: all_potential_type,&
28 : get_potential,&
29 : gth_potential_type,&
30 : sgp_potential_type
31 : USE hartree_local_methods, ONLY: calculate_Vh_1center
32 : USE input_section_types, ONLY: section_get_ivals,&
33 : section_vals_get_subs_vals,&
34 : section_vals_type,&
35 : section_vals_val_get
36 : USE kinds, ONLY: default_string_length,&
37 : dp
38 : USE mathconstants, ONLY: rootpi,&
39 : twopi
40 : USE message_passing, ONLY: mp_para_env_type
41 : USE particle_list_types, ONLY: particle_list_type
42 : USE particle_types, ONLY: particle_type
43 : USE pw_env_types, ONLY: pw_env_get,&
44 : pw_env_type
45 : USE pw_methods, ONLY: pw_axpy,&
46 : pw_copy,&
47 : pw_derive,&
48 : pw_transfer,&
49 : pw_zero
50 : USE pw_poisson_methods, ONLY: pw_poisson_solve
51 : USE pw_poisson_types, ONLY: pw_poisson_type
52 : USE pw_pool_types, ONLY: pw_pool_type
53 : USE pw_types, ONLY: pw_c1d_gs_type,&
54 : pw_r3d_rs_type
55 : USE qs_cneo_types, ONLY: cneo_potential_type
56 : USE qs_environment_types, ONLY: get_qs_env,&
57 : qs_environment_type
58 : USE qs_gapw_densities, ONLY: prepare_gapw_den
59 : USE qs_grid_atom, ONLY: grid_atom_type
60 : USE qs_harmonics_atom, ONLY: harmonics_atom_type
61 : USE qs_kind_types, ONLY: get_qs_kind,&
62 : qs_kind_type
63 : USE qs_ks_methods, ONLY: calc_rho_tot_gspace
64 : USE qs_ks_types, ONLY: qs_ks_env_type
65 : USE qs_linres_types, ONLY: epr_env_type,&
66 : get_epr_env,&
67 : nablavks_atom_type
68 : ! R0
69 : USE qs_local_rho_types, ONLY: rhoz_type
70 : USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
71 : USE qs_oce_types, ONLY: oce_matrix_type
72 : USE qs_rho0_types, ONLY: rho0_atom_type
73 : USE qs_rho_atom_methods, ONLY: calculate_rho_atom_coeff
74 : USE qs_rho_atom_types, ONLY: get_rho_atom,&
75 : rho_atom_coeff,&
76 : rho_atom_type
77 : ! R0
78 : USE qs_rho_types, ONLY: qs_rho_get,&
79 : qs_rho_p_type,&
80 : qs_rho_type
81 : USE qs_subsys_types, ONLY: qs_subsys_get,&
82 : qs_subsys_type
83 : USE qs_vxc, ONLY: qs_vxc_create
84 : USE qs_vxc_atom, ONLY: calculate_vxc_atom_epr
85 : USE util, ONLY: get_limit
86 : #include "./base/base_uses.f90"
87 :
88 : IMPLICIT NONE
89 :
90 : PRIVATE
91 : PUBLIC :: epr_nablavks
92 :
93 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_linres_epr_nablavks'
94 :
95 : CONTAINS
96 :
97 : ! **************************************************************************************************
98 : !> \brief Evaluates Nabla V_KS on the grids
99 : !> \param epr_env ...
100 : !> \param qs_env ...
101 : !> \par History
102 : !> 06.2006 created [RD]
103 : !> \author RD
104 : ! **************************************************************************************************
105 14 : SUBROUTINE epr_nablavks(epr_env, qs_env)
106 :
107 : TYPE(epr_env_type) :: epr_env
108 : TYPE(qs_environment_type), POINTER :: qs_env
109 :
110 : CHARACTER(LEN=default_string_length) :: ext, filename
111 : COMPLEX(KIND=dp) :: gtemp
112 : INTEGER :: bo_atom(2), ia, iat, iatom, idir, iexp, &
113 : ig, ikind, ir, iso, ispin, ix, iy, iz, &
114 : natom, nexp_ppl, nkind, nspins, &
115 : output_unit, unit_nr
116 : INTEGER, DIMENSION(2, 3) :: bo
117 14 : INTEGER, DIMENSION(:), POINTER :: atom_list
118 : LOGICAL :: gapw, gapw_xc, gth_gspace, ionode, &
119 : make_soft, mpi_io, paw_atom
120 : REAL(KIND=dp) :: alpha, alpha_core, arg, charge, ehartree, exc, exc1, exp_rap, &
121 : gapw_max_alpha, hard_radius, hard_value, soft_value, sqrt_alpha, sqrt_rap
122 14 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: vh1_rad_h, vh1_rad_s
123 : REAL(KIND=dp), DIMENSION(3) :: rap, ratom, roffset, rpoint
124 14 : REAL(KIND=dp), DIMENSION(:), POINTER :: cexp_ppl, rho_rad_z
125 14 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rho_rad_0
126 : TYPE(all_potential_type), POINTER :: all_potential
127 14 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
128 : TYPE(cell_type), POINTER :: cell
129 : TYPE(cneo_potential_type), POINTER :: cneo_potential
130 : TYPE(cp_logger_type), POINTER :: logger
131 14 : TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: rho_ao
132 : TYPE(dft_control_type), POINTER :: dft_control
133 : TYPE(grid_atom_type), POINTER :: grid_atom
134 : TYPE(gth_potential_type), POINTER :: gth_potential
135 : TYPE(harmonics_atom_type), POINTER :: harmonics
136 : TYPE(mp_para_env_type), POINTER :: para_env
137 14 : TYPE(nablavks_atom_type), DIMENSION(:), POINTER :: nablavks_atom_set
138 : TYPE(neighbor_list_set_p_type), DIMENSION(:), &
139 14 : POINTER :: sab
140 : TYPE(oce_matrix_type), POINTER :: oce
141 : TYPE(particle_list_type), POINTER :: particles
142 14 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
143 : TYPE(pw_c1d_gs_type) :: rho_tot_gspace, v_coulomb_gspace, &
144 : v_coulomb_gtemp, v_hartree_gspace, &
145 : v_hartree_gtemp, v_xc_gtemp
146 : TYPE(pw_env_type), POINTER :: pw_env
147 : TYPE(pw_poisson_type), POINTER :: poisson_env
148 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
149 : TYPE(pw_r3d_rs_type) :: v_coulomb_rtemp, v_hartree_rtemp, &
150 : v_xc_rtemp, wf_r
151 14 : TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho1_r, rho2_r, rho_r, v_rspace_new, &
152 14 : v_tau_rspace
153 : TYPE(pw_r3d_rs_type), POINTER :: pwx, pwy, pwz
154 14 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
155 : TYPE(qs_ks_env_type), POINTER :: ks_env
156 14 : TYPE(qs_rho_p_type), DIMENSION(:, :), POINTER :: nablavks_set
157 : TYPE(qs_rho_type), POINTER :: rho, rho_xc
158 : TYPE(qs_subsys_type), POINTER :: subsys
159 14 : TYPE(rho0_atom_type), DIMENSION(:), POINTER :: rho0_atom_set
160 14 : TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: rho_rad_h, rho_rad_s
161 14 : TYPE(rho_atom_coeff), DIMENSION(:, :), POINTER :: nablavks_vec_rad_h, nablavks_vec_rad_s
162 14 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
163 : TYPE(rho_atom_type), POINTER :: rho_atom
164 14 : TYPE(rhoz_type), DIMENSION(:), POINTER :: rhoz_set
165 : TYPE(section_vals_type), POINTER :: g_section, input, lr_section, xc_section
166 : TYPE(sgp_potential_type), POINTER :: sgp_potential
167 :
168 : ! R0
169 : ! R0
170 : ! R0
171 : ! R0
172 : ! R0
173 : ! R0
174 :
175 14 : NULLIFY (auxbas_pw_pool)
176 14 : NULLIFY (cell)
177 14 : NULLIFY (dft_control)
178 14 : NULLIFY (g_section)
179 14 : NULLIFY (logger)
180 14 : NULLIFY (lr_section)
181 14 : NULLIFY (nablavks_set)
182 14 : NULLIFY (nablavks_atom_set) ! R0
183 14 : NULLIFY (nablavks_vec_rad_h) ! R0
184 14 : NULLIFY (nablavks_vec_rad_s) ! R0
185 14 : NULLIFY (para_env)
186 14 : NULLIFY (particle_set)
187 14 : NULLIFY (particles)
188 14 : NULLIFY (poisson_env)
189 14 : NULLIFY (pw_env)
190 14 : NULLIFY (pwx) ! R0
191 14 : NULLIFY (pwy) ! R0
192 14 : NULLIFY (pwz) ! R0
193 14 : NULLIFY (rho)
194 14 : NULLIFY (rho_xc)
195 14 : NULLIFY (rho0_atom_set)
196 14 : NULLIFY (rho_atom_set)
197 14 : NULLIFY (rhoz_set)
198 14 : NULLIFY (subsys)
199 14 : NULLIFY (v_rspace_new)
200 14 : NULLIFY (v_tau_rspace)
201 14 : NULLIFY (xc_section)
202 14 : NULLIFY (input)
203 14 : NULLIFY (ks_env)
204 14 : NULLIFY (rho_r, rho_ao, rho1_r, rho2_r)
205 14 : NULLIFY (oce, qs_kind_set, sab)
206 :
207 28 : logger => cp_get_default_logger()
208 14 : lr_section => section_vals_get_subs_vals(qs_env%input, "PROPERTIES%LINRES")
209 : ionode = logger%para_env%is_source()
210 :
211 : output_unit = cp_print_key_unit_nr(logger, lr_section, "PRINT%PROGRAM_RUN_INFO", &
212 14 : extension=".linresLog")
213 :
214 : ! -------------------------------------
215 : ! Read settings
216 : ! -------------------------------------
217 :
218 : g_section => section_vals_get_subs_vals(lr_section, &
219 14 : "EPR%PRINT%G_TENSOR")
220 :
221 14 : CALL section_vals_val_get(g_section, "gapw_max_alpha", r_val=gapw_max_alpha)
222 :
223 14 : gth_gspace = .TRUE.
224 :
225 : ! -------------------------------------
226 : ! Get nablavks arrays
227 : ! -------------------------------------
228 :
229 : CALL get_epr_env(epr_env, nablavks_set=nablavks_set, & ! R0
230 14 : nablavks_atom_set=nablavks_atom_set) ! R0
231 : ! R0
232 :
233 42 : DO ispin = 1, SIZE(nablavks_set, 2)
234 126 : DO idir = 1, SIZE(nablavks_set, 1)
235 84 : CALL qs_rho_get(nablavks_set(idir, ispin)%rho, rho_r=rho_r)
236 112 : CALL pw_zero(rho_r(1))
237 : END DO
238 : END DO
239 :
240 14 : CALL qs_rho_get(nablavks_set(1, 1)%rho, rho_r=rho_r)
241 14 : pwx => rho_r(1)
242 14 : CALL qs_rho_get(nablavks_set(2, 1)%rho, rho_r=rho_r)
243 14 : pwy => rho_r(1)
244 14 : CALL qs_rho_get(nablavks_set(3, 1)%rho, rho_r=rho_r)
245 14 : pwz => rho_r(1)
246 :
247 56 : roffset = -REAL(MODULO(pwx%pw_grid%npts, 2), dp)*pwx%pw_grid%dr/2.0_dp
248 :
249 : ! -------------------------------------
250 : ! Get grids / atom info
251 : ! -------------------------------------
252 :
253 : CALL get_qs_env(qs_env=qs_env, &
254 : atomic_kind_set=atomic_kind_set, &
255 : qs_kind_set=qs_kind_set, &
256 : input=input, &
257 : cell=cell, &
258 : dft_control=dft_control, &
259 : para_env=para_env, &
260 : particle_set=particle_set, &
261 : pw_env=pw_env, &
262 : rho=rho, &
263 : rho_xc=rho_xc, &
264 : rho_atom_set=rho_atom_set, &
265 : rho0_atom_set=rho0_atom_set, &
266 : rhoz_set=rhoz_set, &
267 : subsys=subsys, &
268 : ks_env=ks_env, &
269 14 : oce=oce, sab_orb=sab)
270 :
271 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
272 14 : poisson_env=poisson_env)
273 :
274 14 : CALL qs_subsys_get(subsys, particles=particles)
275 :
276 14 : gapw = dft_control%qs_control%gapw
277 14 : gapw_xc = dft_control%qs_control%gapw_xc
278 14 : nkind = SIZE(atomic_kind_set)
279 14 : nspins = dft_control%nspins
280 :
281 : ! -------------------------------------
282 : ! Add Hartree potential
283 : ! -------------------------------------
284 :
285 14 : CALL auxbas_pw_pool%create_pw(v_hartree_gspace)
286 14 : CALL auxbas_pw_pool%create_pw(v_hartree_gtemp)
287 14 : CALL auxbas_pw_pool%create_pw(v_hartree_rtemp)
288 14 : CALL auxbas_pw_pool%create_pw(rho_tot_gspace)
289 :
290 14 : IF (gapw) THEN
291 : ! need to rebuild the coeff !
292 10 : CALL qs_rho_get(rho, rho_ao_kp=rho_ao)
293 10 : CALL calculate_rho_atom_coeff(qs_env, rho_ao, rho_atom_set, qs_kind_set, oce, sab, para_env)
294 :
295 10 : CALL prepare_gapw_den(qs_env, do_rho0=.TRUE.)
296 : END IF
297 :
298 14 : CALL pw_zero(rho_tot_gspace)
299 :
300 : CALL calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho, &
301 14 : skip_nuclear_density=.NOT. gapw)
302 :
303 : CALL pw_poisson_solve(poisson_env, rho_tot_gspace, ehartree, &
304 14 : v_hartree_gspace)
305 :
306 : ! -------------------------------------
307 : ! Atomic grids part
308 : ! -------------------------------------
309 :
310 14 : IF (gapw) THEN
311 :
312 30 : DO ikind = 1, nkind ! loop over atom types
313 :
314 20 : NULLIFY (atom_list)
315 20 : NULLIFY (grid_atom)
316 20 : NULLIFY (harmonics)
317 : NULLIFY (rho_rad_z)
318 :
319 20 : rho_rad_z => rhoz_set(ikind)%r_coef
320 :
321 20 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
322 : CALL get_qs_kind(qs_kind_set(ikind), &
323 : grid_atom=grid_atom, &
324 : harmonics=harmonics, &
325 : hard_radius=hard_radius, &
326 : paw_atom=paw_atom, &
327 : zeff=charge, &
328 20 : alpha_core_charge=alpha_core)
329 :
330 30 : IF (paw_atom) THEN
331 :
332 80 : ALLOCATE (vh1_rad_h(grid_atom%nr, harmonics%max_iso_not0))
333 60 : ALLOCATE (vh1_rad_s(grid_atom%nr, harmonics%max_iso_not0))
334 :
335 : ! DO iat = 1, natom ! natom = # atoms for ikind
336 : !
337 : ! iatom = atom_list(iat)
338 : ! ratom = particle_set(iatom)%r
339 : !
340 : ! DO ig = v_hartree_gspace%pw_grid%first_gne0,v_hartree_gspace%pw_grid%ngpts_cut_local
341 : !
342 : ! gtemp = fourpi * charge / cell%deth * &
343 : ! EXP ( - v_hartree_gspace%pw_grid%gsq(ig) / (4.0_dp * alpha_core) ) &
344 : ! / v_hartree_gspace%pw_grid%gsq(ig)
345 : !
346 : ! arg = DOT_PRODUCT(v_hartree_gspace%pw_grid%g(:,ig),ratom)
347 : !
348 : ! gtemp = gtemp * CMPLX(COS(arg),-SIN(arg),KIND=dp)
349 : !
350 : ! v_hartree_gspace%array(ig) = v_hartree_gspace%array(ig) + gtemp
351 : ! END DO
352 : ! IF ( v_hartree_gspace%pw_grid%have_g0 ) v_hartree_gspace%array(1) = 0.0_dp
353 : !
354 : ! END DO
355 :
356 20 : bo_atom = get_limit(natom, para_env%num_pe, para_env%mepos)
357 :
358 35 : DO iat = bo_atom(1), bo_atom(2) ! natomkind = # atoms for ikind
359 :
360 15 : iatom = atom_list(iat)
361 60 : ratom = particle_set(iatom)%r
362 :
363 15 : nablavks_vec_rad_h => nablavks_atom_set(iatom)%nablavks_vec_rad_h
364 15 : nablavks_vec_rad_s => nablavks_atom_set(iatom)%nablavks_vec_rad_s
365 :
366 45 : DO ispin = 1, nspins
367 135 : DO idir = 1, 3
368 229590 : nablavks_vec_rad_h(idir, ispin)%r_coef(:, :) = 0.0_dp
369 229620 : nablavks_vec_rad_s(idir, ispin)%r_coef(:, :) = 0.0_dp
370 : END DO ! idir
371 : END DO ! ispin
372 :
373 15 : rho_atom => rho_atom_set(iatom)
374 15 : NULLIFY (rho_rad_h, rho_rad_s, rho_rad_0)
375 : CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=rho_rad_h, &
376 15 : rho_rad_s=rho_rad_s)
377 15 : rho_rad_0 => rho0_atom_set(iatom)%rho0_rad_h%r_coef
378 9348 : vh1_rad_h = 0.0_dp
379 9348 : vh1_rad_s = 0.0_dp
380 :
381 15 : CALL calculate_Vh_1center(vh1_rad_h, vh1_rad_s, rho_rad_h, rho_rad_s, rho_rad_0, rho_rad_z, grid_atom)
382 :
383 750 : DO ir = 2, grid_atom%nr
384 :
385 735 : IF (grid_atom%rad(ir) >= hard_radius) CYCLE
386 :
387 21435 : DO ia = 1, grid_atom%ng_sphere
388 :
389 : ! hard part
390 :
391 84000 : DO idir = 1, 3
392 63000 : hard_value = 0.0_dp
393 831600 : DO iso = 1, harmonics%max_iso_not0
394 : hard_value = hard_value + &
395 : vh1_rad_h(ir, iso)*harmonics%dslm_dxyz(idir, ia, iso) + &
396 : harmonics%slm(ia, iso)* &
397 : (vh1_rad_h(ir - 1, iso) - vh1_rad_h(ir, iso))/ &
398 : (grid_atom%rad(ir - 1) - grid_atom%rad(ir))* &
399 831600 : (harmonics%a(idir, ia))
400 : END DO
401 84000 : nablavks_vec_rad_h(idir, 1)%r_coef(ir, ia) = hard_value
402 : END DO
403 :
404 : ! soft part
405 :
406 84735 : DO idir = 1, 3
407 63000 : soft_value = 0.0_dp
408 831600 : DO iso = 1, harmonics%max_iso_not0
409 : soft_value = soft_value + &
410 : vh1_rad_s(ir, iso)*harmonics%dslm_dxyz(idir, ia, iso) + &
411 : harmonics%slm(ia, iso)* &
412 : (vh1_rad_s(ir - 1, iso) - vh1_rad_s(ir, iso))/ &
413 : (grid_atom%rad(ir - 1) - grid_atom%rad(ir))* &
414 831600 : (harmonics%a(idir, ia))
415 : END DO
416 84000 : nablavks_vec_rad_s(idir, 1)%r_coef(ir, ia) = soft_value
417 : END DO
418 :
419 : END DO ! ia
420 :
421 : END DO ! ir
422 :
423 80 : DO idir = 1, 3
424 114795 : nablavks_vec_rad_h(idir, 2)%r_coef(:, :) = nablavks_vec_rad_h(idir, 1)%r_coef(:, :)
425 114810 : nablavks_vec_rad_s(idir, 2)%r_coef(:, :) = nablavks_vec_rad_s(idir, 1)%r_coef(:, :)
426 : END DO
427 :
428 : END DO ! iat
429 :
430 20 : DEALLOCATE (vh1_rad_h)
431 20 : DEALLOCATE (vh1_rad_s)
432 :
433 : END IF ! paw_atom
434 :
435 : END DO ! ikind
436 :
437 : END IF ! gapw
438 :
439 14 : CALL pw_copy(v_hartree_gspace, v_hartree_gtemp)
440 14 : CALL pw_derive(v_hartree_gtemp, [1, 0, 0])
441 14 : CALL pw_transfer(v_hartree_gtemp, v_hartree_rtemp)
442 14 : CALL pw_copy(v_hartree_rtemp, pwx)
443 :
444 14 : CALL pw_copy(v_hartree_gspace, v_hartree_gtemp)
445 14 : CALL pw_derive(v_hartree_gtemp, [0, 1, 0])
446 14 : CALL pw_transfer(v_hartree_gtemp, v_hartree_rtemp)
447 14 : CALL pw_copy(v_hartree_rtemp, pwy)
448 :
449 14 : CALL pw_copy(v_hartree_gspace, v_hartree_gtemp)
450 14 : CALL pw_derive(v_hartree_gtemp, [0, 0, 1])
451 14 : CALL pw_transfer(v_hartree_gtemp, v_hartree_rtemp)
452 14 : CALL pw_copy(v_hartree_rtemp, pwz)
453 :
454 14 : CALL auxbas_pw_pool%give_back_pw(v_hartree_gspace)
455 14 : CALL auxbas_pw_pool%give_back_pw(v_hartree_gtemp)
456 14 : CALL auxbas_pw_pool%give_back_pw(v_hartree_rtemp)
457 14 : CALL auxbas_pw_pool%give_back_pw(rho_tot_gspace)
458 :
459 : ! -------------------------------------
460 : ! Add Coulomb potential
461 : ! -------------------------------------
462 :
463 42 : DO ikind = 1, nkind ! loop over atom types
464 :
465 28 : NULLIFY (atom_list)
466 28 : NULLIFY (grid_atom)
467 28 : NULLIFY (harmonics)
468 :
469 28 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
470 : CALL get_qs_kind(qs_kind_set(ikind), &
471 : grid_atom=grid_atom, &
472 : harmonics=harmonics, &
473 : hard_radius=hard_radius, &
474 : gth_potential=gth_potential, &
475 : sgp_potential=sgp_potential, &
476 : all_potential=all_potential, &
477 : cneo_potential=cneo_potential, &
478 28 : paw_atom=paw_atom)
479 :
480 42 : IF (ASSOCIATED(gth_potential)) THEN
481 :
482 12 : NULLIFY (cexp_ppl)
483 :
484 : CALL get_potential(potential=gth_potential, &
485 : zeff=charge, &
486 : alpha_ppl=alpha, &
487 : nexp_ppl=nexp_ppl, &
488 12 : cexp_ppl=cexp_ppl)
489 :
490 12 : sqrt_alpha = SQRT(alpha)
491 :
492 12 : IF (gapw .AND. paw_atom .AND. alpha > gapw_max_alpha) THEN
493 : make_soft = .TRUE.
494 : ELSE
495 8 : make_soft = .FALSE.
496 : END IF
497 :
498 : ! -------------------------------------
499 : ! PW grid part
500 : ! -------------------------------------
501 :
502 : IF (gth_gspace) THEN
503 :
504 12 : CALL auxbas_pw_pool%create_pw(v_coulomb_gspace)
505 12 : CALL auxbas_pw_pool%create_pw(v_coulomb_gtemp)
506 12 : CALL auxbas_pw_pool%create_pw(v_coulomb_rtemp)
507 :
508 12 : CALL pw_zero(v_coulomb_gspace)
509 :
510 30 : DO iat = 1, natom ! natom = # atoms for ikind
511 :
512 18 : iatom = atom_list(iat)
513 72 : ratom = particle_set(iatom)%r
514 :
515 820134 : DO ig = v_coulomb_gspace%pw_grid%first_gne0, v_coulomb_gspace%pw_grid%ngpts_cut_local
516 820116 : gtemp = 0.0_dp
517 : ! gtemp = - fourpi * charge / cell%deth * &
518 : ! EXP ( - v_coulomb_gspace%pw_grid%gsq(ig) / (4.0_dp * alpha) ) &
519 : ! / v_coulomb_gspace%pw_grid%gsq(ig)
520 :
521 820116 : IF (.NOT. make_soft) THEN
522 :
523 0 : SELECT CASE (nexp_ppl)
524 : CASE (1)
525 : gtemp = gtemp + &
526 : (twopi)**(1.5_dp)/(cell%deth*(2.0_dp*alpha)**(1.5_dp))* &
527 : EXP(-v_coulomb_gspace%pw_grid%gsq(ig)/(4.0_dp*alpha))*( &
528 : ! C1
529 : +cexp_ppl(1) &
530 0 : )
531 : CASE (2)
532 : gtemp = gtemp + &
533 : (twopi)**(1.5_dp)/(cell%deth*(2.0_dp*alpha)**(1.5_dp))* &
534 : EXP(-v_coulomb_gspace%pw_grid%gsq(ig)/(4.0_dp*alpha))*( &
535 : ! C1
536 : +cexp_ppl(1) &
537 : ! C2
538 : + cexp_ppl(2)/(2.0_dp*alpha)* &
539 : (3.0_dp - v_coulomb_gspace%pw_grid%gsq(ig)/(2.0_dp*alpha)) &
540 546744 : )
541 : CASE (3)
542 : gtemp = gtemp + &
543 : (twopi)**(1.5_dp)/(cell%deth*(2.0_dp*alpha)**(1.5_dp))* &
544 : EXP(-v_coulomb_gspace%pw_grid%gsq(ig)/(4.0_dp*alpha))*( &
545 : ! C1
546 : +cexp_ppl(1) &
547 : ! C2
548 : + cexp_ppl(2)/(2.0_dp*alpha)* &
549 : (3.0_dp - v_coulomb_gspace%pw_grid%gsq(ig)/(2.0_dp*alpha)) &
550 : ! C3
551 : + cexp_ppl(3)/(2.0_dp*alpha)**2* &
552 : (15.0_dp - 10.0_dp*v_coulomb_gspace%pw_grid%gsq(ig)/(2.0_dp*alpha) &
553 : + (v_coulomb_gspace%pw_grid%gsq(ig)/(2.0_dp*alpha))**2) &
554 0 : )
555 : CASE (4)
556 : gtemp = gtemp + &
557 : (twopi)**(1.5_dp)/(cell%deth*(2.0_dp*alpha)**(1.5_dp))* &
558 : EXP(-v_coulomb_gspace%pw_grid%gsq(ig)/(4.0_dp*alpha))*( &
559 : ! C1
560 : +cexp_ppl(1) &
561 : ! C2
562 : + cexp_ppl(2)/(2.0_dp*alpha)* &
563 : (3.0_dp - v_coulomb_gspace%pw_grid%gsq(ig)/(2.0_dp*alpha)) &
564 : ! C3
565 : + cexp_ppl(3)/(2.0_dp*alpha)**2* &
566 : (15.0_dp - 10.0_dp*v_coulomb_gspace%pw_grid%gsq(ig)/(2.0_dp*alpha) &
567 : + (v_coulomb_gspace%pw_grid%gsq(ig)/(2.0_dp*alpha))**2) &
568 : ! C4
569 : + cexp_ppl(4)/(2.0_dp*alpha)**3* &
570 : (105.0_dp - 105.0_dp*v_coulomb_gspace%pw_grid%gsq(ig)/(2.0_dp*alpha) &
571 : + 21.0_dp*(v_coulomb_gspace%pw_grid%gsq(ig)/(2.0_dp*alpha))**2 &
572 : - (v_coulomb_gspace%pw_grid%gsq(ig)/(2.0_dp*alpha))**3) &
573 546744 : )
574 : END SELECT
575 :
576 : END IF
577 :
578 3280464 : arg = DOT_PRODUCT(v_coulomb_gspace%pw_grid%g(:, ig), ratom)
579 :
580 820116 : gtemp = gtemp*CMPLX(COS(arg), -SIN(arg), KIND=dp)
581 820134 : v_coulomb_gspace%array(ig) = v_coulomb_gspace%array(ig) + gtemp
582 : END DO
583 30 : IF (v_coulomb_gspace%pw_grid%have_g0) v_coulomb_gspace%array(1) = 0.0_dp
584 :
585 : END DO
586 :
587 12 : CALL pw_copy(v_coulomb_gspace, v_coulomb_gtemp)
588 12 : CALL pw_derive(v_coulomb_gtemp, [1, 0, 0])
589 12 : CALL pw_transfer(v_coulomb_gtemp, v_coulomb_rtemp)
590 12 : CALL pw_axpy(v_coulomb_rtemp, pwx)
591 :
592 12 : CALL pw_copy(v_coulomb_gspace, v_coulomb_gtemp)
593 12 : CALL pw_derive(v_coulomb_gtemp, [0, 1, 0])
594 12 : CALL pw_transfer(v_coulomb_gtemp, v_coulomb_rtemp)
595 12 : CALL pw_axpy(v_coulomb_rtemp, pwy)
596 :
597 12 : CALL pw_copy(v_coulomb_gspace, v_coulomb_gtemp)
598 12 : CALL pw_derive(v_coulomb_gtemp, [0, 0, 1])
599 12 : CALL pw_transfer(v_coulomb_gtemp, v_coulomb_rtemp)
600 12 : CALL pw_axpy(v_coulomb_rtemp, pwz)
601 :
602 12 : CALL auxbas_pw_pool%give_back_pw(v_coulomb_gspace)
603 12 : CALL auxbas_pw_pool%give_back_pw(v_coulomb_gtemp)
604 12 : CALL auxbas_pw_pool%give_back_pw(v_coulomb_rtemp)
605 : ELSE
606 :
607 : ! Attic of the atomic parallellisation
608 : !
609 : ! bo(2)
610 : ! bo = get_limit(natom, para_env%num_pe, para_env%mepos)
611 : ! DO iat = bo(1),bo(2) ! natom = # atoms for ikind
612 : ! DO ix = lbound(pwx%array,1), ubound(pwx%array,1)
613 : ! DO iy = lbound(pwx%array,2), ubound(pwx%array,2)
614 : ! DO iz = lbound(pwx%array,3), ubound(pwx%array,3)
615 :
616 : bo = pwx%pw_grid%bounds_local
617 :
618 : DO iat = 1, natom ! natom = # atoms for ikind
619 :
620 : iatom = atom_list(iat)
621 : ratom = particle_set(iatom)%r
622 :
623 : DO ix = bo(1, 1), bo(2, 1)
624 : DO iy = bo(1, 2), bo(2, 2)
625 : DO iz = bo(1, 3), bo(2, 3)
626 : rpoint = [REAL(ix, dp)*pwx%pw_grid%dr(1), &
627 : REAL(iy, dp)*pwx%pw_grid%dr(2), &
628 : REAL(iz, dp)*pwx%pw_grid%dr(3)]
629 : rpoint = rpoint + roffset
630 : rap = rpoint - ratom
631 : rap(1) = MODULO(rap(1), cell%hmat(1, 1)) - cell%hmat(1, 1)/2._dp
632 : rap(2) = MODULO(rap(2), cell%hmat(2, 2)) - cell%hmat(2, 2)/2._dp
633 : rap(3) = MODULO(rap(3), cell%hmat(3, 3)) - cell%hmat(3, 3)/2._dp
634 : sqrt_rap = SQRT(DOT_PRODUCT(rap, rap))
635 : exp_rap = EXP(-alpha*sqrt_rap**2)
636 : sqrt_rap = MAX(sqrt_rap, 1.e-10_dp)
637 : ! d_x
638 :
639 : pwx%array(ix, iy, iz) = pwx%array(ix, iy, iz) + charge*( &
640 : -2.0_dp*sqrt_alpha*EXP(-sqrt_rap**2*sqrt_alpha**2)*rap(1) &
641 : /(rootpi*sqrt_rap**2) &
642 : + erf(sqrt_rap*sqrt_alpha)*rap(1) &
643 : /sqrt_rap**3)
644 :
645 : ! d_y
646 :
647 : pwy%array(ix, iy, iz) = pwy%array(ix, iy, iz) + charge*( &
648 : -2.0_dp*sqrt_alpha*EXP(-sqrt_rap**2*sqrt_alpha**2)*rap(2) &
649 : /(rootpi*sqrt_rap**2) &
650 : + erf(sqrt_rap*sqrt_alpha)*rap(2) &
651 : /sqrt_rap**3)
652 :
653 : ! d_z
654 :
655 : pwz%array(ix, iy, iz) = pwz%array(ix, iy, iz) + charge*( &
656 : -2.0_dp*sqrt_alpha*EXP(-sqrt_rap**2*sqrt_alpha**2)*rap(3) &
657 : /(rootpi*sqrt_rap**2) &
658 : + erf(sqrt_rap*sqrt_alpha)*rap(3) &
659 : /sqrt_rap**3)
660 :
661 : IF (make_soft) CYCLE
662 :
663 : ! d_x
664 :
665 : DO iexp = 1, nexp_ppl
666 : pwx%array(ix, iy, iz) = pwx%array(ix, iy, iz) + ( &
667 : -2.0_dp*alpha*rap(1)*exp_rap* &
668 : cexp_ppl(iexp)*(sqrt_rap**2)**(iexp - 1))
669 : IF (iexp > 1) THEN
670 : pwx%array(ix, iy, iz) = pwx%array(ix, iy, iz) + ( &
671 : 2.0_dp*exp_rap*cexp_ppl(iexp)* &
672 : (sqrt_rap**2)**(iexp - 2)*REAL(iexp - 1, dp)*rap(1))
673 : END IF
674 : END DO
675 :
676 : ! d_y
677 :
678 : DO iexp = 1, nexp_ppl
679 : pwy%array(ix, iy, iz) = pwy%array(ix, iy, iz) + ( &
680 : -2.0_dp*alpha*rap(2)*exp_rap* &
681 : cexp_ppl(iexp)*(sqrt_rap**2)**(iexp - 1))
682 : IF (iexp > 1) THEN
683 : pwy%array(ix, iy, iz) = pwy%array(ix, iy, iz) + ( &
684 : 2.0_dp*exp_rap*cexp_ppl(iexp)* &
685 : (sqrt_rap**2)**(iexp - 2)*REAL(iexp - 1, dp)*rap(2))
686 : END IF
687 : END DO
688 :
689 : ! d_z
690 :
691 : DO iexp = 1, nexp_ppl
692 : pwz%array(ix, iy, iz) = pwz%array(ix, iy, iz) + ( &
693 : -2.0_dp*alpha*rap(3)*exp_rap* &
694 : cexp_ppl(iexp)*(sqrt_rap**2)**(iexp - 1))
695 : IF (iexp > 1) THEN
696 : pwz%array(ix, iy, iz) = pwz%array(ix, iy, iz) + ( &
697 : 2.0_dp*exp_rap*cexp_ppl(iexp)* &
698 : (sqrt_rap**2)**(iexp - 2)*REAL(iexp - 1, dp)*rap(3))
699 : END IF
700 : END DO
701 :
702 : END DO ! iz
703 : END DO ! iy
704 : END DO ! ix
705 : END DO ! iat
706 : END IF ! gth_gspace
707 :
708 : ! -------------------------------------
709 : ! Atomic grids part
710 : ! -------------------------------------
711 :
712 12 : IF (gapw .AND. paw_atom) THEN
713 :
714 4 : bo_atom = get_limit(natom, para_env%num_pe, para_env%mepos)
715 :
716 7 : DO iat = bo_atom(1), bo_atom(2) ! natom = # atoms for ikind
717 :
718 3 : iatom = atom_list(iat)
719 :
720 3 : nablavks_vec_rad_h => nablavks_atom_set(iatom)%nablavks_vec_rad_h
721 3 : nablavks_vec_rad_s => nablavks_atom_set(iatom)%nablavks_vec_rad_s
722 :
723 153 : DO ir = 1, grid_atom%nr
724 :
725 150 : IF (grid_atom%rad(ir) >= hard_radius) CYCLE
726 :
727 84 : exp_rap = EXP(-alpha*grid_atom%rad(ir)**2)
728 :
729 4287 : DO ia = 1, grid_atom%ng_sphere
730 :
731 16950 : DO idir = 1, 3
732 12600 : hard_value = 0.0_dp
733 : hard_value = charge*( &
734 : -2.0_dp*sqrt_alpha*EXP(-grid_atom%rad(ir)**2*sqrt_alpha**2) &
735 : *grid_atom%rad(ir)*harmonics%a(idir, ia) &
736 : /(rootpi*grid_atom%rad(ir)**2) &
737 : + erf(grid_atom%rad(ir)*sqrt_alpha) &
738 : *grid_atom%rad(ir)*harmonics%a(idir, ia) &
739 12600 : /grid_atom%rad(ir)**3)
740 12600 : soft_value = hard_value
741 37800 : DO iexp = 1, nexp_ppl
742 : hard_value = hard_value + ( &
743 : -2.0_dp*alpha*grid_atom%rad(ir)*harmonics%a(idir, ia) &
744 25200 : *exp_rap*cexp_ppl(iexp)*(grid_atom%rad(ir)**2)**(iexp - 1))
745 37800 : IF (iexp > 1) THEN
746 : hard_value = hard_value + ( &
747 : 2.0_dp*exp_rap*cexp_ppl(iexp) &
748 : *(grid_atom%rad(ir)**2)**(iexp - 2)*REAL(iexp - 1, dp) &
749 12600 : *grid_atom%rad(ir)*harmonics%a(idir, ia))
750 : END IF
751 : END DO
752 : nablavks_vec_rad_h(idir, 1)%r_coef(ir, ia) = &
753 12600 : nablavks_vec_rad_h(idir, 1)%r_coef(ir, ia) + hard_value
754 16800 : IF (make_soft) THEN
755 : nablavks_vec_rad_s(idir, 1)%r_coef(ir, ia) = &
756 12600 : nablavks_vec_rad_s(idir, 1)%r_coef(ir, ia) + soft_value
757 : ELSE
758 : nablavks_vec_rad_s(idir, 1)%r_coef(ir, ia) = &
759 0 : nablavks_vec_rad_s(idir, 1)%r_coef(ir, ia) + hard_value
760 : END IF
761 : END DO
762 :
763 : END DO ! ia
764 : END DO ! ir
765 :
766 10 : DO ispin = 2, nspins
767 15 : DO idir = 1, 3
768 22959 : nablavks_vec_rad_h(idir, ispin)%r_coef(:, :) = nablavks_vec_rad_h(idir, 1)%r_coef(:, :)
769 22962 : nablavks_vec_rad_s(idir, ispin)%r_coef(:, :) = nablavks_vec_rad_s(idir, 1)%r_coef(:, :)
770 : END DO
771 : END DO
772 :
773 : END DO
774 :
775 : END IF
776 :
777 16 : ELSE IF (ASSOCIATED(sgp_potential)) THEN
778 :
779 0 : CPABORT("EPR with SGP potentials is not implemented")
780 :
781 16 : ELSE IF (ASSOCIATED(cneo_potential)) THEN
782 :
783 0 : CPABORT("EPR with CNEO potentials is not implemented")
784 :
785 16 : ELSE IF (ASSOCIATED(all_potential)) THEN
786 :
787 : CALL get_potential(potential=all_potential, &
788 : alpha_core_charge=alpha, &
789 16 : zeff=charge)
790 :
791 16 : sqrt_alpha = SQRT(alpha)
792 :
793 : ! -------------------------------------
794 : ! Atomic grids part
795 : ! -------------------------------------
796 :
797 16 : bo_atom = get_limit(natom, para_env%num_pe, para_env%mepos)
798 :
799 28 : DO iat = bo_atom(1), bo_atom(2) ! natom = # atoms for ikind
800 :
801 12 : iatom = atom_list(iat)
802 :
803 12 : nablavks_vec_rad_h => nablavks_atom_set(iatom)%nablavks_vec_rad_h
804 :
805 612 : DO ir = 1, grid_atom%nr
806 :
807 600 : IF (grid_atom%rad(ir) >= hard_radius) CYCLE
808 :
809 17148 : DO ia = 1, grid_atom%ng_sphere
810 :
811 67800 : DO idir = 1, 3
812 50400 : hard_value = 0.0_dp
813 : hard_value = charge*( &
814 : 2.0_dp*sqrt_alpha*EXP(-grid_atom%rad(ir)**2*sqrt_alpha**2) &
815 : *grid_atom%rad(ir)*harmonics%a(idir, ia) &
816 : /(rootpi*grid_atom%rad(ir)**2) &
817 : + erfc(grid_atom%rad(ir)*sqrt_alpha) &
818 : *grid_atom%rad(ir)*harmonics%a(idir, ia) &
819 50400 : /grid_atom%rad(ir)**3)
820 : nablavks_vec_rad_h(idir, 1)%r_coef(ir, ia) = &
821 67200 : nablavks_vec_rad_h(idir, 1)%r_coef(ir, ia) + hard_value
822 : END DO
823 :
824 : END DO ! ia
825 : END DO ! ir
826 :
827 40 : DO ispin = 2, nspins
828 60 : DO idir = 1, 3
829 91848 : nablavks_vec_rad_h(idir, ispin)%r_coef(:, :) = nablavks_vec_rad_h(idir, 1)%r_coef(:, :)
830 : END DO
831 : END DO
832 :
833 : END DO
834 :
835 : ELSE
836 : CYCLE
837 : END IF
838 :
839 : END DO
840 :
841 56 : DO idir = 1, 3
842 42 : CALL qs_rho_get(nablavks_set(idir, 1)%rho, rho_r=rho1_r)
843 42 : CALL qs_rho_get(nablavks_set(idir, 2)%rho, rho_r=rho2_r)
844 56 : CALL pw_copy(rho1_r(1), rho2_r(1))
845 : END DO
846 :
847 : ! -------------------------------------
848 : ! Add V_xc potential
849 : ! -------------------------------------
850 :
851 14 : CALL auxbas_pw_pool%create_pw(v_xc_gtemp)
852 14 : CALL auxbas_pw_pool%create_pw(v_xc_rtemp)
853 :
854 14 : xc_section => section_vals_get_subs_vals(input, "PROPERTIES%LINRES%EPR%PRINT%G_TENSOR%XC")
855 : ! a possible vdW section in xc_section will be ignored
856 :
857 14 : IF (gapw_xc) THEN
858 : CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho_xc, xc_section=xc_section, &
859 0 : vxc_rho=v_rspace_new, vxc_tau=v_tau_rspace, exc=exc, just_energy=.FALSE.)
860 : ELSE
861 : CALL qs_vxc_create(ks_env=ks_env, rho_struct=rho, xc_section=xc_section, &
862 14 : vxc_rho=v_rspace_new, vxc_tau=v_tau_rspace, exc=exc, just_energy=.FALSE.)
863 : END IF
864 :
865 14 : IF (ASSOCIATED(v_rspace_new)) THEN
866 :
867 42 : DO ispin = 1, nspins
868 :
869 28 : CALL pw_transfer(v_rspace_new(ispin), v_xc_gtemp)
870 28 : CALL pw_derive(v_xc_gtemp, [1, 0, 0])
871 28 : CALL pw_transfer(v_xc_gtemp, v_xc_rtemp)
872 28 : CALL qs_rho_get(nablavks_set(1, ispin)%rho, rho_r=rho_r)
873 28 : CALL pw_axpy(v_xc_rtemp, rho_r(1))
874 :
875 28 : CALL pw_transfer(v_rspace_new(ispin), v_xc_gtemp)
876 28 : CALL pw_derive(v_xc_gtemp, [0, 1, 0])
877 28 : CALL pw_transfer(v_xc_gtemp, v_xc_rtemp)
878 28 : CALL qs_rho_get(nablavks_set(2, ispin)%rho, rho_r=rho_r)
879 28 : CALL pw_axpy(v_xc_rtemp, rho_r(1))
880 :
881 28 : CALL pw_transfer(v_rspace_new(ispin), v_xc_gtemp)
882 28 : CALL pw_derive(v_xc_gtemp, [0, 0, 1])
883 28 : CALL pw_transfer(v_xc_gtemp, v_xc_rtemp)
884 28 : CALL qs_rho_get(nablavks_set(3, ispin)%rho, rho_r=rho_r)
885 28 : CALL pw_axpy(v_xc_rtemp, rho_r(1))
886 :
887 42 : CALL auxbas_pw_pool%give_back_pw(v_rspace_new(ispin))
888 :
889 : END DO
890 :
891 14 : DEALLOCATE (v_rspace_new)
892 : END IF
893 :
894 14 : IF (ASSOCIATED(v_tau_rspace)) THEN
895 :
896 0 : DO ispin = 1, nspins
897 :
898 0 : CALL pw_transfer(v_tau_rspace(ispin), v_xc_gtemp)
899 0 : CALL pw_derive(v_xc_gtemp, [1, 0, 0])
900 0 : CALL pw_transfer(v_xc_gtemp, v_xc_rtemp)
901 0 : CALL qs_rho_get(nablavks_set(1, ispin)%rho, rho_r=rho_r)
902 0 : CALL pw_axpy(v_xc_rtemp, rho_r(1))
903 :
904 0 : CALL pw_transfer(v_tau_rspace(ispin), v_xc_gtemp)
905 0 : CALL pw_derive(v_xc_gtemp, [0, 1, 0])
906 0 : CALL pw_transfer(v_xc_gtemp, v_xc_rtemp)
907 0 : CALL qs_rho_get(nablavks_set(2, ispin)%rho, rho_r=rho_r)
908 0 : CALL pw_axpy(v_xc_rtemp, rho_r(1))
909 :
910 0 : CALL pw_transfer(v_tau_rspace(ispin), v_xc_gtemp)
911 0 : CALL pw_derive(v_xc_gtemp, [0, 0, 1])
912 0 : CALL pw_transfer(v_xc_gtemp, v_xc_rtemp)
913 0 : CALL qs_rho_get(nablavks_set(3, ispin)%rho, rho_r=rho_r)
914 0 : CALL pw_axpy(v_xc_rtemp, rho_r(1))
915 :
916 0 : CALL auxbas_pw_pool%give_back_pw(v_tau_rspace(ispin))
917 :
918 : END DO
919 :
920 0 : DEALLOCATE (v_tau_rspace)
921 : END IF
922 :
923 14 : CALL auxbas_pw_pool%give_back_pw(v_xc_gtemp)
924 14 : CALL auxbas_pw_pool%give_back_pw(v_xc_rtemp)
925 :
926 14 : IF (gapw .OR. gapw_xc) THEN
927 10 : CALL calculate_vxc_atom_epr(qs_env=qs_env, exc1=exc1, gradient_atom_set=nablavks_atom_set)
928 : END IF
929 :
930 : ! -------------------------------------
931 : ! Write Nabla V_KS (local) to cubes
932 : ! -------------------------------------
933 :
934 14 : IF (BTEST(cp_print_key_should_output(logger%iter_info, lr_section, &
935 : "EPR%PRINT%NABLAVKS_CUBES"), cp_p_file)) THEN
936 0 : CALL auxbas_pw_pool%create_pw(wf_r)
937 0 : DO idir = 1, 3
938 0 : CALL pw_zero(wf_r)
939 0 : CALL qs_rho_get(nablavks_set(idir, 1)%rho, rho_r=rho_r)
940 0 : CALL pw_copy(rho_r(1), wf_r) ! RA
941 0 : filename = "nablavks"
942 0 : mpi_io = .TRUE.
943 0 : WRITE (ext, '(a2,I1,a5)') "_d", idir, ".cube"
944 : unit_nr = cp_print_key_unit_nr(logger, lr_section, "EPR%PRINT%NABLAVKS_CUBES", &
945 : extension=TRIM(ext), middle_name=TRIM(filename), &
946 : log_filename=.FALSE., file_position="REWIND", &
947 0 : mpi_io=mpi_io)
948 : CALL cp_pw_to_cube(wf_r, unit_nr, "NABLA V_KS ", &
949 : particles=particles, &
950 : stride=section_get_ivals(lr_section, &
951 : "EPR%PRINT%NABLAVKS_CUBES%STRIDE"), &
952 0 : mpi_io=mpi_io)
953 : CALL cp_print_key_finished_output(unit_nr, logger, lr_section, &
954 0 : "EPR%PRINT%NABLAVKS_CUBES", mpi_io=mpi_io)
955 : END DO
956 0 : CALL auxbas_pw_pool%give_back_pw(wf_r)
957 : END IF
958 :
959 : CALL cp_print_key_finished_output(output_unit, logger, lr_section, &
960 14 : "PRINT%PROGRAM_RUN_INFO")
961 :
962 28 : END SUBROUTINE epr_nablavks
963 :
964 : END MODULE qs_linres_epr_nablavks
965 :
|
