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