(git:58e3e09)
xc_pot_saop.F
Go to the documentation of this file.
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 Calculate the saop potential
10 ! **************************************************************************************************
12  USE atomic_kind_types, ONLY: atomic_kind_type,&
14  USE basis_set_types, ONLY: gto_basis_set_type
15  USE cp_array_utils, ONLY: cp_1d_r_p_type
16  USE cp_control_types, ONLY: dft_control_type
20  USE cp_fm_types, ONLY: cp_fm_create,&
23  cp_fm_p_type,&
24  cp_fm_release,&
27  cp_fm_type
28  USE dbcsr_api, ONLY: dbcsr_copy,&
29  dbcsr_deallocate_matrix,&
30  dbcsr_p_type,&
31  dbcsr_set
32  USE input_constants, ONLY: do_method_gapw,&
33  oe_gllb,&
34  oe_lb,&
35  oe_saop,&
37  USE input_section_types, ONLY: &
41  USE kinds, ONLY: dp
42  USE mathconstants, ONLY: pi
43  USE message_passing, ONLY: mp_para_env_type
44  USE pw_env_types, ONLY: pw_env_get,&
45  pw_env_type
46  USE pw_methods, ONLY: pw_axpy,&
47  pw_copy,&
48  pw_scale,&
49  pw_zero
50  USE pw_pool_types, ONLY: pw_pool_type
51  USE pw_types, ONLY: pw_c1d_gs_type,&
52  pw_r3d_rs_type
54  USE qs_environment_types, ONLY: get_qs_env,&
55  qs_environment_type
57  USE qs_grid_atom, ONLY: grid_atom_type
58  USE qs_harmonics_atom, ONLY: harmonics_atom_type
59  USE qs_integrate_potential, ONLY: integrate_v_rspace
60  USE qs_kind_types, ONLY: get_qs_kind,&
61  qs_kind_type
62  USE qs_ks_atom, ONLY: update_ks_atom
63  USE qs_ks_types, ONLY: qs_ks_env_type
66  local_rho_type
67  USE qs_mo_types, ONLY: get_mo_set,&
68  mo_set_type
69  USE qs_neighbor_list_types, ONLY: neighbor_list_set_p_type
70  USE qs_oce_types, ONLY: oce_matrix_type
73  USE qs_rho_atom_types, ONLY: get_rho_atom,&
74  rho_atom_coeff,&
75  rho_atom_type
76  USE qs_rho_types, ONLY: qs_rho_get,&
77  qs_rho_type
78  USE qs_vxc_atom, ONLY: calc_rho_angular,&
80  USE util, ONLY: get_limit
81  USE virial_types, ONLY: virial_type
82  USE xc, ONLY: xc_vxc_pw_create1
83  USE xc_atom, ONLY: fill_rho_set,&
84  vxc_of_r_new,&
86  USE xc_derivative_set_types, ONLY: xc_derivative_set_type,&
92  xc_derivative_type
95  xc_rho_cflags_type
98  xc_rho_set_type,&
100  USE xc_xbecke88, ONLY: xb88_lda_info,&
102 #include "./base/base_uses.f90"
103 
104  IMPLICIT NONE
105 
106  PRIVATE
107 
108  PUBLIC :: add_saop_pot
109 
110  CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'xc_pot_saop'
111 
112  ! should be eliminated
113  REAL(KIND=dp), PARAMETER :: alpha = 1.19_dp, beta = 0.01_dp, k_rho = 0.42_dp
114  REAL(KIND=dp), PARAMETER :: kappa = 0.804_dp, mu = 0.21951_dp, &
115  beta_ec = 0.066725_dp, gamma_saop = 0.031091_dp
116 
117 CONTAINS
118 
119 ! **************************************************************************************************
120 !> \brief ...
121 !> \param ks_matrix ...
122 !> \param qs_env ...
123 !> \param oe_corr ...
124 ! **************************************************************************************************
125  SUBROUTINE add_saop_pot(ks_matrix, qs_env, oe_corr)
126 
127  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: ks_matrix
128  TYPE(qs_environment_type), POINTER :: qs_env
129  INTEGER, INTENT(IN) :: oe_corr
130 
131  INTEGER :: dft_method_id, homo, i, ispin, j, k, &
132  nspins, orb, xc_deriv_method_id, &
133  xc_rho_smooth_id
134  INTEGER, DIMENSION(2) :: ncol, nrow
135  INTEGER, DIMENSION(2, 3) :: bo
136  LOGICAL :: compute_virial, gapw, lsd
137  REAL(kind=dp) :: density_cut, efac, gradient_cut, &
138  tau_cut, we_gllb, we_lb, xc_energy
139  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: coeff_col
140  REAL(kind=dp), DIMENSION(3, 3) :: virial_xc_tmp
141  REAL(kind=dp), DIMENSION(:), POINTER :: mo_eigenvalues
142  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: e_uniform
143  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: single_mo_coeff
144  TYPE(cp_fm_type), POINTER :: mo_coeff
145  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: orbital_density_matrix, rho_struct_ao
146  TYPE(mo_set_type), DIMENSION(:), POINTER :: molecular_orbitals
147  TYPE(pw_c1d_gs_type) :: orbital_g
148  TYPE(pw_c1d_gs_type), DIMENSION(:), POINTER :: rho_g
149  TYPE(pw_env_type), POINTER :: pw_env
150  TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
151  TYPE(pw_r3d_rs_type) :: orbital
152  TYPE(pw_r3d_rs_type), ALLOCATABLE, DIMENSION(:) :: vxc_gllb, vxc_saop
153  TYPE(pw_r3d_rs_type), DIMENSION(:), POINTER :: rho_r, rho_struct_r, tau, vxc_lb, &
154  vxc_tau, vxc_tmp
155  TYPE(qs_ks_env_type), POINTER :: ks_env
156  TYPE(qs_rho_type), POINTER :: rho_struct
157  TYPE(section_vals_type), POINTER :: input, xc_fun_section_orig, &
158  xc_fun_section_tmp, xc_section_orig, &
159  xc_section_tmp
160  TYPE(virial_type), POINTER :: virial
161  TYPE(xc_derivative_set_type) :: deriv_set
162  TYPE(xc_derivative_type), POINTER :: deriv
163  TYPE(xc_rho_cflags_type) :: needs
164  TYPE(xc_rho_set_type) :: rho_set
165 
166  NULLIFY (ks_env, pw_env, auxbas_pw_pool, input)
167  NULLIFY (rho_g, rho_r, tau, rho_struct, e_uniform)
168  NULLIFY (vxc_lb, vxc_tmp, vxc_tau)
169  NULLIFY (mo_eigenvalues, deriv, rho_struct_r, rho_struct_ao)
170  NULLIFY (orbital_density_matrix, xc_section_tmp, xc_fun_section_tmp)
171 
172  CALL get_qs_env(qs_env, &
173  ks_env=ks_env, &
174  rho=rho_struct, &
175  pw_env=pw_env, &
176  input=input, &
177  virial=virial, &
178  mos=molecular_orbitals)
179  compute_virial = virial%pv_calculate .AND. (.NOT. virial%pv_numer)
180  CALL section_vals_val_get(input, "DFT%QS%METHOD", i_val=dft_method_id)
181  gapw = (dft_method_id == do_method_gapw)
182 
183  xc_section_orig => section_vals_get_subs_vals(input, "DFT%XC")
184  CALL section_vals_retain(xc_section_orig)
185  CALL section_vals_duplicate(xc_section_orig, xc_section_tmp)
186 
187  CALL section_vals_val_get(xc_section_orig, "DENSITY_CUTOFF", &
188  r_val=density_cut)
189  CALL section_vals_val_get(xc_section_orig, "GRADIENT_CUTOFF", &
190  r_val=gradient_cut)
191  CALL section_vals_val_get(xc_section_orig, "TAU_CUTOFF", &
192  r_val=tau_cut)
193 
194  CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
195 
196  CALL section_vals_val_get(input, "DFT%LSD", l_val=lsd)
197  IF (lsd) THEN
198  nspins = 2
199  ELSE
200  nspins = 1
201  END IF
202 
203  ALLOCATE (single_mo_coeff(nspins))
204  CALL dbcsr_allocate_matrix_set(orbital_density_matrix, nspins)
205  CALL qs_rho_get(rho_struct, rho_r=rho_struct_r, rho_ao=rho_struct_ao)
206  rho_r => rho_struct_r
207  DO ispin = 1, nspins
208  ALLOCATE (orbital_density_matrix(ispin)%matrix)
209  CALL dbcsr_copy(orbital_density_matrix(ispin)%matrix, &
210  rho_struct_ao(ispin)%matrix, "orbital density")
211  END DO
212  bo = rho_r(1)%pw_grid%bounds_local
213 
214  !---------------------------!
215  ! create the density needed !
216  !---------------------------!
217  CALL xc_rho_set_create(rho_set, bo, &
218  density_cut, &
219  gradient_cut, &
220  tau_cut)
221  CALL xc_rho_cflags_setall(needs, .false.)
222  IF (lsd) THEN
223  CALL xb88_lsd_info(needs=needs)
224  needs%norm_drho = .true.
225  ELSE
226  CALL xb88_lda_info(needs=needs)
227  END IF
228  CALL section_vals_val_get(xc_section_orig, "XC_GRID%XC_DERIV", &
229  i_val=xc_deriv_method_id)
230  CALL section_vals_val_get(xc_section_orig, "XC_GRID%XC_SMOOTH_RHO", &
231  i_val=xc_rho_smooth_id)
232  CALL xc_rho_set_update(rho_set, rho_r, rho_g, tau, needs, &
233  xc_deriv_method_id, &
234  xc_rho_smooth_id, &
235  auxbas_pw_pool)
236 
237  !----------------------------------------!
238  ! Construct the LB94 potential in vxc_LB !
239  !----------------------------------------!
240  xc_fun_section_orig => section_vals_get_subs_vals(xc_section_orig, &
241  "XC_FUNCTIONAL")
242  CALL section_vals_create(xc_fun_section_tmp, xc_fun_section_orig%section)
243  CALL section_vals_val_set(xc_fun_section_tmp, "_SECTION_PARAMETERS_", &
244  i_val=xc_funct_no_shortcut)
245  CALL section_vals_val_set(xc_fun_section_tmp, "XALPHA%_SECTION_PARAMETERS_", &
246  l_val=.true.)
247  CALL section_vals_set_subs_vals(xc_section_tmp, "XC_FUNCTIONAL", &
248  xc_fun_section_tmp)
249 
250  cpassert(.NOT. compute_virial)
251 ! CALL xc_vxc_pw_create(vxc_tmp, vxc_tau, xc_energy, rho_r, rho_g, tau, &
252 ! xc_section_tmp, auxbas_pw_pool, &
253 ! compute_virial=.FALSE., virial_xc=virial_xc_tmp)
254  CALL xc_vxc_pw_create1(vxc_tmp, vxc_tau, rho_r, rho_g, tau, xc_energy, &
255  xc_section_tmp, auxbas_pw_pool, &
256  compute_virial=.false., virial_xc=virial_xc_tmp)
257 
258  CALL section_vals_val_set(xc_fun_section_tmp, "XALPHA%_SECTION_PARAMETERS_", &
259  l_val=.false.)
260  CALL section_vals_val_set(xc_fun_section_tmp, "PZ81%_SECTION_PARAMETERS_", &
261  l_val=.true.)
262 
263  cpassert(.NOT. compute_virial)
264 ! CALL xc_vxc_pw_create(vxc_LB, vxc_tau, xc_energy, rho_r, rho_g, tau, &
265 ! xc_section_tmp, auxbas_pw_pool, &
266 ! compute_virial=.FALSE., virial_xc=virial_xc_tmp)
267  CALL xc_vxc_pw_create1(vxc_lb, vxc_tau, rho_r, rho_g, tau, xc_energy, &
268  xc_section_tmp, auxbas_pw_pool, &
269  compute_virial=.false., virial_xc=virial_xc_tmp)
270 
271  DO ispin = 1, nspins
272  CALL pw_axpy(vxc_tmp(ispin), vxc_lb(ispin), alpha)
273  END DO
274 
275  DO ispin = 1, nspins
276  CALL add_lb_pot(vxc_tmp(ispin)%array, rho_set, lsd, ispin)
277  CALL pw_axpy(vxc_tmp(ispin), vxc_lb(ispin), -1.0_dp)
278  END DO
279 
280  !-----------------------------------------------------------------------------------!
281  ! Construct 2 times PBE one particle density from the PZ correlation energy density !
282  !-----------------------------------------------------------------------------------!
283  CALL xc_dset_create(deriv_set, local_bounds=bo)
284  CALL xc_functionals_eval(xc_fun_section_tmp, &
285  lsd=lsd, &
286  rho_set=rho_set, &
287  deriv_set=deriv_set, &
288  deriv_order=0)
289 
290  deriv => xc_dset_get_derivative(deriv_set, [INTEGER::])
291  CALL xc_derivative_get(deriv, deriv_data=e_uniform)
292 
293  ALLOCATE (vxc_gllb(nspins))
294  DO ispin = 1, nspins
295  CALL auxbas_pw_pool%create_pw(vxc_gllb(ispin))
296  END DO
297 
298  DO ispin = 1, nspins
299  CALL calc_2excpbe(vxc_gllb(ispin)%array, rho_set, e_uniform, lsd)
300  END DO
301 
302  CALL xc_dset_release(deriv_set)
303 
304  CALL auxbas_pw_pool%create_pw(orbital)
305  CALL auxbas_pw_pool%create_pw(orbital_g)
306 
307  DO ispin = 1, nspins
308 
309  CALL get_mo_set(molecular_orbitals(ispin), &
310  mo_coeff=mo_coeff, &
311  eigenvalues=mo_eigenvalues, &
312  homo=homo)
313  CALL cp_fm_create(single_mo_coeff(ispin), &
314  mo_coeff%matrix_struct, &
315  "orbital density matrix")
316 
317  CALL cp_fm_get_info(single_mo_coeff(ispin), &
318  nrow_global=nrow(ispin), ncol_global=ncol(ispin))
319  ALLOCATE (coeff_col(nrow(ispin), 1))
320 
321  CALL pw_zero(vxc_tmp(ispin))
322 
323  DO orb = 1, homo - 1
324 
325  efac = k_rho*sqrt(mo_eigenvalues(homo) - mo_eigenvalues(orb))
326  IF (.NOT. lsd) efac = 2.0_dp*efac
327 
328  CALL cp_fm_set_all(single_mo_coeff(ispin), 0.0_dp)
329  CALL cp_fm_get_submatrix(mo_coeff, coeff_col, &
330  1, orb, nrow(ispin), 1)
331  CALL cp_fm_set_submatrix(single_mo_coeff(ispin), coeff_col, &
332  1, orb)
333  CALL dbcsr_set(orbital_density_matrix(ispin)%matrix, 0.0_dp)
334  CALL cp_dbcsr_plus_fm_fm_t(orbital_density_matrix(ispin)%matrix, &
335  matrix_v=single_mo_coeff(ispin), &
336  ncol=ncol(ispin), &
337  alpha=1.0_dp)
338  CALL pw_zero(orbital)
339  CALL pw_zero(orbital_g)
340  CALL calculate_rho_elec(matrix_p=orbital_density_matrix(ispin)%matrix, &
341  rho=orbital, rho_gspace=orbital_g, &
342  ks_env=ks_env)
343 
344  CALL pw_axpy(orbital, vxc_tmp(ispin), efac)
345 
346  END DO
347  DEALLOCATE (coeff_col)
348 
349  DO k = bo(1, 3), bo(2, 3)
350  DO j = bo(1, 2), bo(2, 2)
351  DO i = bo(1, 1), bo(2, 1)
352  IF (rho_r(ispin)%array(i, j, k) > density_cut) THEN
353  vxc_tmp(ispin)%array(i, j, k) = vxc_tmp(ispin)%array(i, j, k)/ &
354  rho_r(ispin)%array(i, j, k)
355  ELSE
356  vxc_tmp(ispin)%array(i, j, k) = 0.0_dp
357  END IF
358  END DO
359  END DO
360  END DO
361 
362  CALL pw_axpy(vxc_tmp(ispin), vxc_gllb(ispin), 1.0_dp)
363 
364  END DO
365 
366  !---------------!
367  ! Assemble SAOP !
368  !---------------!
369  ALLOCATE (vxc_saop(nspins))
370 
371  DO ispin = 1, nspins
372 
373  CALL get_mo_set(molecular_orbitals(ispin), &
374  mo_coeff=mo_coeff, &
375  eigenvalues=mo_eigenvalues, &
376  homo=homo)
377  CALL auxbas_pw_pool%create_pw(vxc_saop(ispin))
378  CALL pw_zero(vxc_saop(ispin))
379 
380  ALLOCATE (coeff_col(nrow(ispin), 1))
381 
382  DO orb = 1, homo
383 
384  we_lb = exp(-2.0_dp*(mo_eigenvalues(homo) - mo_eigenvalues(orb))**2)
385  we_gllb = 1.0_dp - we_lb
386  IF (.NOT. lsd) THEN
387  we_lb = 2.0_dp*we_lb
388  we_gllb = 2.0_dp*we_gllb
389  END IF
390 
391  vxc_tmp(ispin)%array = we_lb*vxc_lb(ispin)%array + &
392  we_gllb*vxc_gllb(ispin)%array
393 
394  CALL cp_fm_set_all(single_mo_coeff(ispin), 0.0_dp)
395  CALL cp_fm_get_submatrix(mo_coeff, coeff_col, &
396  1, orb, nrow(ispin), 1)
397  CALL cp_fm_set_submatrix(single_mo_coeff(ispin), coeff_col, &
398  1, orb)
399  CALL dbcsr_set(orbital_density_matrix(ispin)%matrix, 0.0_dp)
400  CALL cp_dbcsr_plus_fm_fm_t(orbital_density_matrix(ispin)%matrix, &
401  matrix_v=single_mo_coeff(ispin), &
402  ncol=ncol(ispin), &
403  alpha=1.0_dp)
404  CALL pw_zero(orbital)
405  CALL pw_zero(orbital_g)
406  CALL calculate_rho_elec(matrix_p=orbital_density_matrix(ispin)%matrix, &
407  rho=orbital, rho_gspace=orbital_g, &
408  ks_env=ks_env)
409 
410  vxc_saop(ispin)%array = vxc_saop(ispin)%array + &
411  orbital%array*vxc_tmp(ispin)%array
412 
413  END DO
414 
415  CALL dbcsr_deallocate_matrix(orbital_density_matrix(ispin)%matrix)
416 
417  DEALLOCATE (coeff_col)
418 
419  DO k = bo(1, 3), bo(2, 3)
420  DO j = bo(1, 2), bo(2, 2)
421  DO i = bo(1, 1), bo(2, 1)
422  IF (rho_r(ispin)%array(i, j, k) > density_cut) THEN
423  vxc_saop(ispin)%array(i, j, k) = vxc_saop(ispin)%array(i, j, k)/ &
424  rho_r(ispin)%array(i, j, k)
425  ELSE
426  vxc_saop(ispin)%array(i, j, k) = 0.0_dp
427  END IF
428  END DO
429  END DO
430  END DO
431 
432  END DO
433 
434  CALL cp_fm_release(single_mo_coeff)
435 
436  CALL xc_rho_set_release(rho_set, auxbas_pw_pool)
437  CALL auxbas_pw_pool%give_back_pw(orbital)
438  CALL auxbas_pw_pool%give_back_pw(orbital_g)
439 
440  !--------------------!
441  ! Do the integration !
442  !--------------------!
443  DO ispin = 1, nspins
444 
445  IF (oe_corr == oe_lb) THEN
446  CALL pw_copy(vxc_lb(ispin), vxc_saop(ispin))
447  ELSE IF (oe_corr == oe_gllb) THEN
448  CALL pw_copy(vxc_gllb(ispin), vxc_saop(ispin))
449  END IF
450  CALL pw_scale(vxc_saop(ispin), vxc_saop(ispin)%pw_grid%dvol)
451 
452  CALL integrate_v_rspace(v_rspace=vxc_saop(ispin), pmat=rho_struct_ao(ispin), &
453  hmat=ks_matrix(ispin), qs_env=qs_env, &
454  calculate_forces=.false., &
455  gapw=gapw)
456 
457  END DO
458 
459  DO ispin = 1, nspins
460  CALL auxbas_pw_pool%give_back_pw(vxc_saop(ispin))
461  CALL auxbas_pw_pool%give_back_pw(vxc_gllb(ispin))
462  CALL vxc_lb(ispin)%release()
463  CALL vxc_tmp(ispin)%release()
464  END DO
465  DEALLOCATE (vxc_gllb, vxc_lb, vxc_tmp, orbital_density_matrix)
466 
467  DEALLOCATE (vxc_saop)
468 
469  CALL section_vals_release(xc_fun_section_tmp)
470  CALL section_vals_release(xc_section_tmp)
471  CALL section_vals_release(xc_section_orig)
472 
473  !-----------------------!
474  ! Call the GAPW routine !
475  !-----------------------!
476  IF (gapw) THEN
477  CALL gapw_add_atomic_saop_pot(qs_env, oe_corr)
478  END IF
479 
480  END SUBROUTINE add_saop_pot
481 
482 ! **************************************************************************************************
483 !> \brief ...
484 !> \param qs_env ...
485 !> \param oe_corr ...
486 ! **************************************************************************************************
487  SUBROUTINE gapw_add_atomic_saop_pot(qs_env, oe_corr)
488 
489  TYPE(qs_environment_type), POINTER :: qs_env
490  INTEGER, INTENT(IN) :: oe_corr
491 
492  INTEGER :: ia, iat, iatom, ikind, ir, ispin, na, &
493  natom, nr, ns, nspins, orb
494  INTEGER, DIMENSION(2) :: bo, homo, ncol, nrow
495  INTEGER, DIMENSION(2, 3) :: bounds
496  INTEGER, DIMENSION(:), POINTER :: atom_list
497  LOGICAL :: lsd, paw_atom
498  REAL(dp), DIMENSION(:, :, :), POINTER :: tau
499  REAL(kind=dp) :: density_cut, efac, exc, gradient_cut, &
500  tau_cut, we_gllb, we_lb
501  REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: coeff_col
502  REAL(kind=dp), DIMENSION(:, :), POINTER :: weight
503  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: dummy, e_uniform, rho_h, rho_s, vtau, &
504  vxc_gllb_h, vxc_gllb_s, vxc_lb_h, vxc_lb_s, vxc_saop_h, vxc_saop_s, vxc_tmp_h, vxc_tmp_s
505  REAL(kind=dp), DIMENSION(:, :, :, :), POINTER :: drho_h, drho_s, vxg
506  TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
507  TYPE(cp_1d_r_p_type), DIMENSION(:), POINTER :: mo_eigenvalues
508  TYPE(cp_fm_p_type), ALLOCATABLE, DIMENSION(:) :: mo_coeff
509  TYPE(cp_fm_type), ALLOCATABLE, DIMENSION(:) :: single_mo_coeff
510  TYPE(dbcsr_p_type), DIMENSION(:), POINTER :: matrix_ks, orbital_density_matrix, &
511  rho_struct_ao
512  TYPE(dbcsr_p_type), DIMENSION(:, :), POINTER :: ksmat, psmat
513  TYPE(dft_control_type), POINTER :: dft_control
514  TYPE(grid_atom_type), POINTER :: atomic_grid, grid_atom
515  TYPE(gto_basis_set_type), POINTER :: orb_basis
516  TYPE(harmonics_atom_type), POINTER :: harmonics
517  TYPE(local_rho_type), POINTER :: local_rho_set
518  TYPE(mo_set_type), DIMENSION(:), POINTER :: molecular_orbitals
519  TYPE(mp_para_env_type), POINTER :: para_env
520  TYPE(neighbor_list_set_p_type), DIMENSION(:), &
521  POINTER :: sab
522  TYPE(oce_matrix_type), POINTER :: oce
523  TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
524  TYPE(qs_rho_type), POINTER :: rho_structure
525  TYPE(rho_atom_coeff), DIMENSION(:), POINTER :: dr_h, dr_s, int_hh, int_ss, r_h, r_s
526  TYPE(rho_atom_coeff), DIMENSION(:, :), POINTER :: r_h_d, r_s_d
527  TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
528  TYPE(rho_atom_type), POINTER :: rho_atom
529  TYPE(section_vals_type), POINTER :: input, xc_fun_section_orig, &
530  xc_fun_section_tmp, xc_section_orig, &
531  xc_section_tmp
532  TYPE(xc_derivative_set_type) :: deriv_set
533  TYPE(xc_derivative_type), POINTER :: deriv
534  TYPE(xc_rho_cflags_type) :: needs, needs_orbs
535  TYPE(xc_rho_set_type) :: orb_rho_set_h, orb_rho_set_s, rho_set_h, &
536  rho_set_s
537 
538  NULLIFY (weight, rho_h, rho_s, vxc_lb_h, vxc_lb_s, vxc_gllb_h, vxc_gllb_s, &
539  vxc_tmp_h, vxc_tmp_s, vtau, dummy, e_uniform, drho_h, drho_s, vxg, atom_list, &
540  atomic_kind_set, qs_kind_set, deriv, atomic_grid, rho_struct_ao, &
541  harmonics, molecular_orbitals, rho_structure, r_h, r_s, dr_h, dr_s, &
542  r_h_d, r_s_d, rho_atom_set, rho_atom, para_env, &
543  mo_eigenvalues, local_rho_set, matrix_ks, &
544  orbital_density_matrix, vxc_saop_h, vxc_saop_s)
545 
546  ! tau is needed for fill_rho_set, but should never be used
547  NULLIFY (tau)
548  NULLIFY (dft_control, oce, sab)
549 
550  CALL get_qs_env(qs_env, input=input, &
551  rho=rho_structure, &
552  mos=molecular_orbitals, &
553  atomic_kind_set=atomic_kind_set, &
554  qs_kind_set=qs_kind_set, &
555  rho_atom_set=rho_atom_set, &
556  matrix_ks=matrix_ks, &
557  dft_control=dft_control, &
558  para_env=para_env, &
559  oce=oce, sab_orb=sab)
560 
561  CALL qs_rho_get(rho_structure, rho_ao=rho_struct_ao)
562 
563  xc_section_orig => section_vals_get_subs_vals(input, "DFT%XC")
564  CALL section_vals_retain(xc_section_orig)
565  CALL section_vals_duplicate(xc_section_orig, xc_section_tmp)
566 
567  ! [SC] the following code can be traced back to SVN rev. 4296 (git:f97138b) that
568  ! has removed the component 'nspins' from the derived type 'dft_control_type'.
569  ! Is it worth to remove the code below in favour of 'dft_control%nspins'
570  ! since its reintroduction? Note that in case of ROKS calculations,
571  ! 'lsd == .FALSE.' but 'dft_control%nspins == 2'.
572  CALL section_vals_val_get(input, "DFT%LSD", l_val=lsd)
573  IF (lsd) THEN
574  nspins = 2
575  ELSE
576  nspins = 1
577  END IF
578 
579  CALL section_vals_val_get(xc_section_orig, "DENSITY_CUTOFF", &
580  r_val=density_cut)
581  CALL section_vals_val_get(xc_section_orig, "GRADIENT_CUTOFF", &
582  r_val=gradient_cut)
583  CALL section_vals_val_get(xc_section_orig, "TAU_CUTOFF", &
584  r_val=tau_cut)
585 
586  ! remap pointer
587  ns = SIZE(rho_struct_ao)
588  psmat(1:ns, 1:1) => rho_struct_ao(1:ns)
589  CALL calculate_rho_atom_coeff(qs_env, psmat, rho_atom_set, qs_kind_set, oce, sab, para_env)
590  CALL prepare_gapw_den(qs_env)
591 
592  ALLOCATE (mo_coeff(nspins), single_mo_coeff(nspins), mo_eigenvalues(nspins))
593 
594  CALL dbcsr_allocate_matrix_set(orbital_density_matrix, nspins)
595 
596  DO ispin = 1, nspins
597  CALL get_mo_set(molecular_orbitals(ispin), &
598  mo_coeff=mo_coeff(ispin)%matrix, &
599  eigenvalues=mo_eigenvalues(ispin)%array, &
600  homo=homo(ispin))
601  CALL cp_fm_create(single_mo_coeff(ispin), &
602  mo_coeff(ispin)%matrix%matrix_struct, &
603  "orbital density matrix")
604  CALL cp_fm_get_info(single_mo_coeff(ispin), &
605  nrow_global=nrow(ispin), ncol_global=ncol(ispin))
606  ALLOCATE (orbital_density_matrix(ispin)%matrix)
607  CALL dbcsr_copy(orbital_density_matrix(ispin)%matrix, &
608  rho_struct_ao(ispin)%matrix, &
609  "orbital density")
610  END DO
611  CALL local_rho_set_create(local_rho_set)
612  CALL allocate_rho_atom_internals(local_rho_set%rho_atom_set, atomic_kind_set, &
613  qs_kind_set, dft_control, para_env)
614 
615  DO ikind = 1, SIZE(atomic_kind_set)
616  CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natom)
617 
618  CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom, &
619  harmonics=harmonics, grid_atom=atomic_grid)
620  IF (.NOT. paw_atom) cycle
621 
622  nr = atomic_grid%nr
623  na = atomic_grid%ng_sphere
624  bounds(1:2, 1:3) = 1
625  bounds(2, 1) = na
626  bounds(2, 2) = nr
627 
628  CALL xc_dset_create(deriv_set, local_bounds=bounds)
629 
630  CALL xc_rho_set_create(rho_set_h, bounds, density_cut, &
631  gradient_cut, tau_cut)
632  CALL xc_rho_set_create(rho_set_s, bounds, density_cut, &
633  gradient_cut, tau_cut)
634  CALL xc_rho_set_create(orb_rho_set_h, bounds, density_cut, &
635  gradient_cut, tau_cut)
636  CALL xc_rho_set_create(orb_rho_set_s, bounds, density_cut, &
637  gradient_cut, tau_cut)
638 
639  CALL xc_rho_cflags_setall(needs, .false.)
640  IF (lsd) THEN
641  CALL xb88_lsd_info(needs=needs)
642  needs%norm_drho = .true.
643  ELSE
644  CALL xb88_lda_info(needs=needs)
645  END IF
646  CALL xc_rho_set_atom_update(rho_set_h, needs, nspins, bounds)
647  CALL xc_rho_set_atom_update(rho_set_s, needs, nspins, bounds)
648  CALL xc_rho_cflags_setall(needs_orbs, .false.)
649  needs_orbs%rho = .true.
650  IF (lsd) needs_orbs%rho_spin = .true.
651  CALL xc_rho_set_atom_update(orb_rho_set_h, needs, nspins, bounds)
652  CALL xc_rho_set_atom_update(orb_rho_set_s, needs, nspins, bounds)
653 
654  ALLOCATE (rho_h(1:na, 1:nr, 1:nspins), rho_s(1:na, 1:nr, 1:nspins))
655  ALLOCATE (weight(1:na, 1:nr))
656  ALLOCATE (vxc_lb_h(1:na, 1:nr, 1:nspins), vxc_lb_s(1:na, 1:nr, 1:nspins))
657  ALLOCATE (vxc_gllb_h(1:na, 1:nr, 1:nspins), vxc_gllb_s(1:na, 1:nr, 1:nspins))
658  ALLOCATE (vxc_tmp_h(1:na, 1:nr, 1:nspins), vxc_tmp_s(1:na, 1:nr, 1:nspins))
659  ALLOCATE (vxc_saop_h(1:na, 1:nr, 1:nspins), vxc_saop_s(1:na, 1:nr, 1:nspins))
660  ALLOCATE (drho_h(1:4, 1:na, 1:nr, 1:nspins), drho_s(1:4, 1:na, 1:nr, 1:nspins))
661 
662  ! Distribute the atoms of this kind
663  bo = get_limit(natom, para_env%num_pe, para_env%mepos)
664 
665  DO iat = 1, natom !bo(1),bo(2)
666  iatom = atom_list(iat)
667 
668  rho_atom => rho_atom_set(iatom)
669  NULLIFY (r_h, r_s, dr_h, dr_s, r_h_d, r_s_d)
670  CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=r_h, &
671  rho_rad_s=r_s, drho_rad_h=dr_h, &
672  drho_rad_s=dr_s, rho_rad_h_d=r_h_d, &
673  rho_rad_s_d=r_s_d)
674  rho_h = 0.0_dp
675  rho_s = 0.0_dp
676  drho_h = 0.0_dp
677  drho_s = 0.0_dp
678  DO ir = 1, nr
679  CALL calc_rho_angular(atomic_grid, harmonics, nspins, .true., &
680  ir, r_h, r_s, rho_h, rho_s, &
681  dr_h, dr_s, r_h_d, r_s_d, drho_h, drho_s)
682  END DO
683  DO ir = 1, nr
684  CALL fill_rho_set(rho_set_h, lsd, nspins, needs, rho_h, drho_h, tau, na, ir)
685  CALL fill_rho_set(rho_set_s, lsd, nspins, needs, rho_s, drho_s, tau, na, ir)
686  END DO
687  DO ir = 1, nr
688  DO ia = 1, na
689  weight(ia, ir) = atomic_grid%wr(ir)*atomic_grid%wa(ia)
690  END DO
691  END DO
692 
693  !-----------------------------!
694  ! 1. Slater exchange for LB94 !
695  !-----------------------------!
696  xc_fun_section_orig => section_vals_get_subs_vals(xc_section_orig, &
697  "XC_FUNCTIONAL")
698  CALL section_vals_create(xc_fun_section_tmp, xc_fun_section_orig%section)
699  CALL section_vals_val_set(xc_fun_section_tmp, "_SECTION_PARAMETERS_", &
700  i_val=xc_funct_no_shortcut)
701  CALL section_vals_val_set(xc_fun_section_tmp, "XALPHA%_SECTION_PARAMETERS_", &
702  l_val=.true.)
703  CALL section_vals_set_subs_vals(xc_section_tmp, "XC_FUNCTIONAL", &
704  xc_fun_section_tmp)
705 
706  !---------------------!
707  ! Both: hard and soft !
708  !---------------------!
709  CALL xc_dset_zero_all(deriv_set)
710  CALL vxc_of_r_new(xc_fun_section_tmp, rho_set_h, deriv_set, 1, needs, &
711  weight, lsd, na, nr, exc, vxc_tmp_h, vxg, vtau)
712  CALL xc_dset_zero_all(deriv_set)
713  CALL vxc_of_r_new(xc_fun_section_tmp, rho_set_s, deriv_set, 1, needs, &
714  weight, lsd, na, nr, exc, vxc_tmp_s, vxg, vtau)
715 
716  !-------------------------------------------!
717  ! 2. PZ correlation for LB94 and ec_uniform !
718  !-------------------------------------------!
719  CALL section_vals_val_set(xc_fun_section_tmp, "XALPHA%_SECTION_PARAMETERS_", &
720  l_val=.false.)
721  CALL section_vals_val_set(xc_fun_section_tmp, "PZ81%_SECTION_PARAMETERS_", &
722  l_val=.true.)
723 
724  !------!
725  ! Hard !
726  !------!
727  CALL xc_dset_zero_all(deriv_set)
728  CALL vxc_of_r_new(xc_fun_section_tmp, rho_set_h, deriv_set, 1, needs, &
729  weight, lsd, na, nr, exc, vxc_lb_h, vxg, vtau)
730  vxc_lb_h = vxc_lb_h + alpha*vxc_tmp_h
731  DO ispin = 1, nspins
732  dummy => vxc_tmp_h(:, :, ispin:ispin)
733  CALL add_lb_pot(dummy, rho_set_h, lsd, ispin)
734  vxc_lb_h(:, :, ispin) = vxc_lb_h(:, :, ispin) - weight(:, :)*vxc_tmp_h(:, :, ispin)
735  END DO
736  NULLIFY (dummy)
737 
738  vxc_gllb_h = 0.0_dp
739  deriv => xc_dset_get_derivative(deriv_set, [INTEGER::])
740  cpassert(ASSOCIATED(deriv))
741  CALL xc_derivative_get(deriv, deriv_data=e_uniform)
742  DO ispin = 1, nspins
743  dummy => vxc_gllb_h(:, :, ispin:ispin)
744  CALL calc_2excpbe(dummy, rho_set_h, e_uniform, lsd)
745  vxc_gllb_h(:, :, ispin) = vxc_gllb_h(:, :, ispin)*weight(:, :)
746  END DO
747  NULLIFY (deriv, dummy, e_uniform)
748 
749  !------!
750  ! Soft !
751  !------!
752  CALL xc_dset_zero_all(deriv_set)
753  CALL vxc_of_r_new(xc_fun_section_tmp, rho_set_s, deriv_set, 1, needs, &
754  weight, lsd, na, nr, exc, vxc_lb_s, vxg, vtau)
755 
756  vxc_lb_s = vxc_lb_s + alpha*vxc_tmp_s
757  DO ispin = 1, nspins
758  dummy => vxc_tmp_s(:, :, ispin:ispin)
759  CALL add_lb_pot(dummy, rho_set_s, lsd, ispin)
760  vxc_lb_s(:, :, ispin) = vxc_lb_s(:, :, ispin) - weight(:, :)*vxc_tmp_s(:, :, ispin)
761  END DO
762  NULLIFY (dummy)
763 
764  vxc_gllb_s = 0.0_dp
765  deriv => xc_dset_get_derivative(deriv_set, [INTEGER::])
766  cpassert(ASSOCIATED(deriv))
767  CALL xc_derivative_get(deriv, deriv_data=e_uniform)
768  DO ispin = 1, nspins
769  dummy => vxc_gllb_s(:, :, ispin:ispin)
770  CALL calc_2excpbe(dummy, rho_set_s, e_uniform, lsd)
771  vxc_gllb_s(:, :, ispin) = vxc_gllb_s(:, :, ispin)*weight(:, :)
772  END DO
773  NULLIFY (deriv, dummy, e_uniform)
774 
775  !------------------!
776  ! Now the orbitals !
777  !------------------!
778  vxc_tmp_h = 0.0_dp; vxc_tmp_s = 0.0_dp
779 
780  DO ispin = 1, nspins
781 
782  DO orb = 1, homo(ispin) - 1
783 
784  ALLOCATE (coeff_col(nrow(ispin), 1))
785 
786  efac = k_rho*sqrt(mo_eigenvalues(ispin)%array(homo(ispin)) - &
787  mo_eigenvalues(ispin)%array(orb))
788  IF (.NOT. lsd) efac = 2.0_dp*efac
789 
790  CALL cp_fm_set_all(single_mo_coeff(ispin), 0.0_dp)
791  CALL cp_fm_get_submatrix(mo_coeff(ispin)%matrix, coeff_col, &
792  1, orb, nrow(ispin), 1)
793  CALL cp_fm_set_submatrix(single_mo_coeff(ispin), coeff_col, &
794  1, orb)
795  CALL dbcsr_set(orbital_density_matrix(ispin)%matrix, 0.0_dp)
796  CALL cp_dbcsr_plus_fm_fm_t(orbital_density_matrix(ispin)%matrix, &
797  matrix_v=single_mo_coeff(ispin), &
798  ncol=ncol(ispin), &
799  alpha=1.0_dp)
800 
801  DEALLOCATE (coeff_col)
802 
803  ! This calculates the CPC and density on the grids for every atom even though
804  ! we need it only for iatom at the moment. It seems that to circumvent this,
805  ! the routines must be adapted to calculate just iatom
806  ! remap pointer
807  ns = SIZE(orbital_density_matrix)
808  psmat(1:ns, 1:1) => orbital_density_matrix(1:ns)
809  CALL calculate_rho_atom_coeff(qs_env, psmat, local_rho_set%rho_atom_set, qs_kind_set, oce, sab, para_env)
810  CALL prepare_gapw_den(qs_env, local_rho_set, .false.)
811 
812  rho_atom => local_rho_set%rho_atom_set(iatom)
813  NULLIFY (r_h, r_s, dr_h, dr_s, r_h_d, r_s_d)
814  CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=r_h, rho_rad_s=r_s)
815  rho_h = 0.0_dp
816  rho_s = 0.0_dp
817  drho_h = 0.0_dp
818  drho_s = 0.0_dp
819  DO ir = 1, nr
820  CALL calc_rho_angular(atomic_grid, harmonics, nspins, .false., &
821  ir, r_h, r_s, rho_h, rho_s, &
822  dr_h, dr_s, r_h_d, r_s_d, drho_h, drho_s)
823  END DO
824  DO ir = 1, nr
825  CALL fill_rho_set(orb_rho_set_h, lsd, nspins, needs_orbs, rho_h, drho_h, tau, na, ir)
826  CALL fill_rho_set(orb_rho_set_s, lsd, nspins, needs_orbs, rho_s, drho_s, tau, na, ir)
827  END DO
828 
829  IF (lsd) THEN
830  IF (ispin == 1) THEN
831  vxc_tmp_h(:, :, 1) = vxc_tmp_h(:, :, 1) + efac*orb_rho_set_h%rhoa(:, :, 1)
832  vxc_tmp_s(:, :, 1) = vxc_tmp_s(:, :, 1) + efac*orb_rho_set_s%rhoa(:, :, 1)
833  ELSE
834  vxc_tmp_h(:, :, 2) = vxc_tmp_h(:, :, 2) + efac*orb_rho_set_h%rhob(:, :, 1)
835  vxc_tmp_s(:, :, 2) = vxc_tmp_s(:, :, 2) + efac*orb_rho_set_s%rhob(:, :, 1)
836  END IF
837  ELSE
838  vxc_tmp_h(:, :, 1) = vxc_tmp_h(:, :, 1) + efac*orb_rho_set_h%rho(:, :, 1)
839  vxc_tmp_s(:, :, 1) = vxc_tmp_s(:, :, 1) + efac*orb_rho_set_s%rho(:, :, 1)
840  END IF
841 
842  END DO ! orb
843 
844  END DO ! ispin
845 
846  IF (lsd) THEN
847  DO ir = 1, nr
848  DO ia = 1, na
849  IF (rho_set_h%rhoa(ia, ir, 1) > rho_set_h%rho_cutoff) &
850  vxc_gllb_h(ia, ir, 1) = vxc_gllb_h(ia, ir, 1) + &
851  weight(ia, ir)*vxc_tmp_h(ia, ir, 1)/rho_set_h%rhoa(ia, ir, 1)
852  IF (rho_set_h%rhob(ia, ir, 1) > rho_set_h%rho_cutoff) &
853  vxc_gllb_h(ia, ir, 2) = vxc_gllb_h(ia, ir, 2) + &
854  weight(ia, ir)*vxc_tmp_h(ia, ir, 2)/rho_set_h%rhob(ia, ir, 1)
855  IF (rho_set_s%rhoa(ia, ir, 1) > rho_set_s%rho_cutoff) &
856  vxc_gllb_s(ia, ir, 1) = vxc_gllb_s(ia, ir, 1) + &
857  weight(ia, ir)*vxc_tmp_s(ia, ir, 1)/rho_set_s%rhoa(ia, ir, 1)
858  IF (rho_set_s%rhob(ia, ir, 1) > rho_set_s%rho_cutoff) &
859  vxc_gllb_s(ia, ir, 2) = vxc_gllb_s(ia, ir, 2) + &
860  weight(ia, ir)*vxc_tmp_s(ia, ir, 2)/rho_set_s%rhob(ia, ir, 1)
861  END DO
862  END DO
863  ELSE
864  DO ir = 1, nr
865  DO ia = 1, na
866  IF (rho_set_h%rho(ia, ir, 1) > rho_set_h%rho_cutoff) &
867  vxc_gllb_h(ia, ir, 1) = vxc_gllb_h(ia, ir, 1) + &
868  weight(ia, ir)*vxc_tmp_h(ia, ir, 1)/rho_set_h%rho(ia, ir, 1)
869  IF (rho_set_s%rho(ia, ir, 1) > rho_set_s%rho_cutoff) &
870  vxc_gllb_s(ia, ir, 1) = vxc_gllb_s(ia, ir, 1) + &
871  weight(ia, ir)*vxc_tmp_s(ia, ir, 1)/rho_set_s%rho(ia, ir, 1)
872  END DO
873  END DO
874  END IF
875 
876  vxc_saop_h = 0.0_dp; vxc_saop_s = 0.0_dp
877 
878  DO ispin = 1, nspins
879 
880  DO orb = 1, homo(ispin)
881 
882  ALLOCATE (coeff_col(nrow(ispin), 1))
883 
884  we_lb = exp(-2.0_dp*(mo_eigenvalues(ispin)%array(homo(ispin)) - &
885  mo_eigenvalues(ispin)%array(orb))**2)
886  we_gllb = 1.0_dp - we_lb
887  IF (.NOT. lsd) THEN
888  we_lb = 2.0_dp*we_lb
889  we_gllb = 2.0_dp*we_gllb
890  END IF
891 
892  vxc_tmp_h(:, :, ispin) = we_lb*vxc_lb_h(:, :, ispin) + &
893  we_gllb*vxc_gllb_h(:, :, ispin)
894  vxc_tmp_s(:, :, ispin) = we_lb*vxc_lb_s(:, :, ispin) + &
895  we_gllb*vxc_gllb_s(:, :, ispin)
896 
897  CALL cp_fm_set_all(single_mo_coeff(ispin), 0.0_dp)
898  CALL cp_fm_get_submatrix(mo_coeff(ispin)%matrix, coeff_col, &
899  1, orb, nrow(ispin), 1)
900  CALL cp_fm_set_submatrix(single_mo_coeff(ispin), coeff_col, &
901  1, orb)
902  CALL dbcsr_set(orbital_density_matrix(ispin)%matrix, 0.0_dp)
903  CALL cp_dbcsr_plus_fm_fm_t(orbital_density_matrix(ispin)%matrix, &
904  matrix_v=single_mo_coeff(ispin), &
905  ncol=ncol(ispin), &
906  alpha=1.0_dp)
907 
908  DEALLOCATE (coeff_col)
909 
910  ! This calculates the CPC and density on the grids for every atom even though
911  ! we need it only for iatom at the moment. It seems that to circumvent this,
912  ! the routines must be adapted to calculate just iatom
913  ! remap pointer
914  ns = SIZE(orbital_density_matrix)
915  psmat(1:ns, 1:1) => orbital_density_matrix(1:ns)
916  CALL calculate_rho_atom_coeff(qs_env, psmat, local_rho_set%rho_atom_set, qs_kind_set, oce, sab, para_env)
917  CALL prepare_gapw_den(qs_env, local_rho_set, .false.)
918 
919  rho_atom => local_rho_set%rho_atom_set(iatom)
920  NULLIFY (r_h, r_s, dr_h, dr_s, r_h_d, r_s_d)
921  CALL get_rho_atom(rho_atom=rho_atom, rho_rad_h=r_h, rho_rad_s=r_s)
922  rho_h = 0.0_dp
923  rho_s = 0.0_dp
924  drho_h = 0.0_dp
925  drho_s = 0.0_dp
926  DO ir = 1, nr
927  CALL calc_rho_angular(atomic_grid, harmonics, nspins, .false., &
928  ir, r_h, r_s, rho_h, rho_s, &
929  dr_h, dr_s, r_h_d, r_s_d, drho_h, drho_s)
930  END DO
931  DO ir = 1, nr
932  CALL fill_rho_set(orb_rho_set_h, lsd, nspins, needs_orbs, rho_h, drho_h, tau, na, ir)
933  CALL fill_rho_set(orb_rho_set_s, lsd, nspins, needs_orbs, rho_s, drho_s, tau, na, ir)
934  END DO
935 
936  IF (lsd) THEN
937  IF (ispin == 1) THEN
938  vxc_saop_h(:, :, 1) = vxc_saop_h(:, :, 1) + vxc_tmp_h(:, :, 1)*orb_rho_set_h%rhoa(:, :, 1)
939  vxc_saop_s(:, :, 1) = vxc_saop_s(:, :, 1) + vxc_tmp_s(:, :, 1)*orb_rho_set_s%rhoa(:, :, 1)
940  ELSE
941  vxc_saop_h(:, :, 2) = vxc_saop_h(:, :, 2) + vxc_tmp_h(:, :, 2)*orb_rho_set_h%rhob(:, :, 1)
942  vxc_saop_s(:, :, 2) = vxc_saop_s(:, :, 2) + vxc_tmp_s(:, :, 2)*orb_rho_set_s%rhob(:, :, 1)
943  END IF
944  ELSE
945  vxc_saop_h(:, :, 1) = vxc_saop_h(:, :, 1) + vxc_tmp_h(:, :, 1)*orb_rho_set_h%rho(:, :, 1)
946  vxc_saop_s(:, :, 1) = vxc_saop_s(:, :, 1) + vxc_tmp_s(:, :, 1)*orb_rho_set_s%rho(:, :, 1)
947  END IF
948 
949  END DO ! orb
950 
951  END DO ! ispin
952 
953  IF (lsd) THEN
954  DO ir = 1, nr
955  DO ia = 1, na
956  IF (rho_set_h%rhoa(ia, ir, 1) > rho_set_h%rho_cutoff) THEN
957  vxc_saop_h(ia, ir, 1) = vxc_saop_h(ia, ir, 1)/rho_set_h%rhoa(ia, ir, 1)
958  ELSE
959  vxc_saop_h(ia, ir, 1) = 0.0_dp
960  END IF
961  IF (rho_set_h%rhob(ia, ir, 1) > rho_set_h%rho_cutoff) THEN
962  vxc_saop_h(ia, ir, 2) = vxc_saop_h(ia, ir, 2)/rho_set_h%rhob(ia, ir, 1)
963  ELSE
964  vxc_saop_h(ia, ir, 2) = 0.0_dp
965  END IF
966  IF (rho_set_s%rhoa(ia, ir, 1) > rho_set_s%rho_cutoff) THEN
967  vxc_saop_s(ia, ir, 1) = vxc_saop_s(ia, ir, 1)/rho_set_s%rhoa(ia, ir, 1)
968  ELSE
969  vxc_saop_s(ia, ir, 1) = 0.0_dp
970  END IF
971  IF (rho_set_s%rhob(ia, ir, 1) > rho_set_s%rho_cutoff) THEN
972  vxc_saop_s(ia, ir, 2) = vxc_saop_s(ia, ir, 2)/rho_set_s%rhob(ia, ir, 1)
973  ELSE
974  vxc_saop_s(ia, ir, 2) = 0.0_dp
975  END IF
976  END DO
977  END DO
978  ELSE
979  DO ir = 1, nr
980  DO ia = 1, na
981  IF (rho_set_h%rho(ia, ir, 1) > rho_set_h%rho_cutoff) THEN
982  vxc_saop_h(ia, ir, 1) = vxc_saop_h(ia, ir, 1)/rho_set_h%rho(ia, ir, 1)
983  ELSE
984  vxc_saop_h(ia, ir, 1) = 0.0_dp
985  END IF
986  IF (rho_set_s%rho(ia, ir, 1) > rho_set_s%rho_cutoff) THEN
987  vxc_saop_s(ia, ir, 1) = vxc_saop_s(ia, ir, 1)/rho_set_s%rho(ia, ir, 1)
988  ELSE
989  vxc_saop_s(ia, ir, 1) = 0.0_dp
990  END IF
991  END DO
992  END DO
993  END IF
994 
995  rho_atom => rho_atom_set(iatom)
996  CALL get_rho_atom(rho_atom=rho_atom, ga_vlocal_gb_h=int_hh, ga_vlocal_gb_s=int_ss)
997  CALL get_qs_kind(qs_kind_set(ikind), basis_set=orb_basis, &
998  harmonics=harmonics, grid_atom=grid_atom)
999  SELECT CASE (oe_corr)
1000  CASE (oe_lb)
1001  CALL gavxcgb_nogc(vxc_lb_h, vxc_lb_s, int_hh, int_ss, grid_atom, orb_basis, harmonics, nspins)
1002  CASE (oe_gllb)
1003  CALL gavxcgb_nogc(vxc_gllb_h, vxc_gllb_s, int_hh, int_ss, grid_atom, orb_basis, harmonics, nspins)
1004  CASE (oe_saop)
1005  CALL gavxcgb_nogc(vxc_saop_h, vxc_saop_s, int_hh, int_ss, grid_atom, orb_basis, harmonics, nspins)
1006  CASE default
1007  cpabort("Unknown correction!")
1008  END SELECT
1009 
1010  END DO
1011 
1012  DEALLOCATE (rho_h, rho_s, weight)
1013  DEALLOCATE (vxc_lb_h, vxc_lb_s)
1014  DEALLOCATE (vxc_gllb_h, vxc_gllb_s)
1015  DEALLOCATE (vxc_tmp_h, vxc_tmp_s)
1016  DEALLOCATE (vxc_saop_h, vxc_saop_s)
1017  DEALLOCATE (drho_h, drho_s)
1018 
1019  CALL xc_dset_release(deriv_set)
1020  CALL xc_rho_set_release(rho_set_h)
1021  CALL xc_rho_set_release(rho_set_s)
1022  CALL xc_rho_set_release(orb_rho_set_h)
1023  CALL xc_rho_set_release(orb_rho_set_s)
1024 
1025  END DO
1026 
1027  ! remap pointer
1028  ns = SIZE(matrix_ks)
1029  ksmat(1:ns, 1:1) => matrix_ks(1:ns)
1030  ns = SIZE(rho_struct_ao)
1031  psmat(1:ns, 1:1) => rho_struct_ao(1:ns)
1032 
1033  CALL update_ks_atom(qs_env, ksmat, psmat, forces=.false.)
1034 
1035  !---------!
1036  ! Cleanup !
1037  !---------!
1038  CALL section_vals_release(xc_fun_section_tmp)
1039  CALL section_vals_release(xc_section_tmp)
1040  CALL section_vals_release(xc_section_orig)
1041 
1042  CALL local_rho_set_release(local_rho_set)
1043  CALL cp_fm_release(single_mo_coeff)
1044  DEALLOCATE (mo_coeff, mo_eigenvalues)
1045  CALL dbcsr_deallocate_matrix_set(orbital_density_matrix)
1046 
1047  END SUBROUTINE gapw_add_atomic_saop_pot
1048 
1049 ! **************************************************************************************************
1050 !> \brief ...
1051 !> \param pot ...
1052 !> \param rho_set ...
1053 !> \param lsd ...
1054 !> \param spin ...
1055 ! **************************************************************************************************
1056  SUBROUTINE add_lb_pot(pot, rho_set, lsd, spin)
1057 
1058  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: pot
1059  TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
1060  LOGICAL, INTENT(IN) :: lsd
1061  INTEGER, INTENT(IN) :: spin
1062 
1063  REAL(kind=dp), PARAMETER :: ob3 = 1.0_dp/3.0_dp
1064 
1065  INTEGER :: i, j, k
1066  INTEGER, DIMENSION(2, 3) :: bo
1067  REAL(kind=dp) :: n, n_13, x, x2
1068 
1069  bo = rho_set%local_bounds
1070 
1071  DO k = bo(1, 3), bo(2, 3)
1072  DO j = bo(1, 2), bo(2, 2)
1073  DO i = bo(1, 1), bo(2, 1)
1074  IF (.NOT. lsd) THEN
1075  IF (rho_set%rho(i, j, k) > rho_set%rho_cutoff) THEN
1076  n = rho_set%rho(i, j, k)/2.0_dp
1077  n_13 = n**ob3
1078  x = (rho_set%norm_drho(i, j, k)/2.0_dp)/(n*n_13)
1079  x2 = x*x
1080  pot(i, j, k) = beta*x2*n_13/(1.0_dp + 3.0_dp*beta*x*log(x + sqrt(x2 + 1.0_dp)))
1081  END IF
1082  ELSE
1083  IF (spin == 1) THEN
1084  IF (rho_set%rhoa(i, j, k) > rho_set%rho_cutoff) THEN
1085  n_13 = rho_set%rhoa_1_3(i, j, k)
1086  x = rho_set%norm_drhoa(i, j, k)/(rho_set%rhoa(i, j, k)*n_13)
1087  x2 = x*x
1088  pot(i, j, k) = beta*x2*n_13/(1.0_dp + 3.0_dp*beta*x*log(sqrt(x2 + 1.0_dp) + x))
1089  END IF
1090  ELSE IF (spin == 2) THEN
1091  IF (rho_set%rhob(i, j, k) > rho_set%rho_cutoff) THEN
1092  n_13 = rho_set%rhob_1_3(i, j, k)
1093  x = rho_set%norm_drhob(i, j, k)/(rho_set%rhob(i, j, k)*n_13)
1094  x2 = x*x
1095  pot(i, j, k) = beta*x2*n_13/(1.0_dp + 3.0_dp*beta*x*log(sqrt(x2 + 1.0_dp) + x))
1096  END IF
1097  END IF
1098  END IF
1099  END DO
1100  END DO
1101  END DO
1102 
1103  END SUBROUTINE add_lb_pot
1104 
1105 ! **************************************************************************************************
1106 !> \brief ...
1107 !> \param pot ...
1108 !> \param rho_set ...
1109 !> \param e_uniform ...
1110 !> \param lsd ...
1111 ! **************************************************************************************************
1112  SUBROUTINE calc_2excpbe(pot, rho_set, e_uniform, lsd)
1113 
1114  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: pot
1115  TYPE(xc_rho_set_type), INTENT(IN) :: rho_set
1116  REAL(kind=dp), DIMENSION(:, :, :), POINTER :: e_uniform
1117  LOGICAL, INTENT(IN) :: lsd
1118 
1119  INTEGER :: i, j, k
1120  INTEGER, DIMENSION(2, 3) :: bo
1121  REAL(kind=dp) :: e_unif, rho
1122 
1123  bo = rho_set%local_bounds
1124 
1125  DO k = bo(1, 3), bo(2, 3)
1126  DO j = bo(1, 2), bo(2, 2)
1127  DO i = bo(1, 1), bo(2, 1)
1128  IF (.NOT. lsd) THEN
1129  IF (rho_set%rho(i, j, k) > rho_set%rho_cutoff) THEN
1130  e_unif = e_uniform(i, j, k)/rho_set%rho(i, j, k)
1131  ELSE
1132  e_unif = 0.0_dp
1133  END IF
1134  pot(i, j, k) = &
1135  2.0_dp* &
1136  calc_ecpbe_r(rho_set%rho(i, j, k), rho_set%norm_drho(i, j, k), &
1137  e_unif, rho_set%rho_cutoff, rho_set%drho_cutoff) + &
1138  2.0_dp* &
1139  calc_expbe_r(rho_set%rho(i, j, k), rho_set%norm_drho(i, j, k), &
1140  rho_set%rho_cutoff, rho_set%drho_cutoff)
1141  ELSE
1142  rho = rho_set%rhoa(i, j, k) + rho_set%rhob(i, j, k)
1143  IF (rho > rho_set%rho_cutoff) THEN
1144  e_unif = e_uniform(i, j, k)/rho
1145  ELSE
1146  e_unif = 0.0_dp
1147  END IF
1148  pot(i, j, k) = &
1149  2.0_dp* &
1150  calc_ecpbe_u(rho_set%rhoa(i, j, k), rho_set%rhob(i, j, k), rho_set%norm_drho(i, j, k), &
1151  e_unif, &
1152  rho_set%rho_cutoff, rho_set%drho_cutoff) + &
1153  2.0_dp* &
1154  calc_expbe_u(rho_set%rhoa(i, j, k), rho_set%rhob(i, j, k), rho_set%norm_drho(i, j, k), &
1155  rho_set%rho_cutoff, rho_set%drho_cutoff)
1156  END IF
1157  END DO
1158  END DO
1159  END DO
1160 
1161  END SUBROUTINE calc_2excpbe
1162 
1163 ! **************************************************************************************************
1164 !> \brief ...
1165 !> \param ra ...
1166 !> \param rb ...
1167 !> \param ngr ...
1168 !> \param ec_unif ...
1169 !> \param rc ...
1170 !> \param ngrc ...
1171 !> \return ...
1172 ! **************************************************************************************************
1173  FUNCTION calc_ecpbe_u(ra, rb, ngr, ec_unif, rc, ngrc) RESULT(res)
1174 
1175  REAL(kind=dp), INTENT(in) :: ra, rb, ngr, ec_unif, rc, ngrc
1176  REAL(kind=dp) :: res
1177 
1178  REAL(kind=dp), PARAMETER :: ob3 = 1.0_dp/3.0_dp, tb3 = 2.0_dp/3.0_dp
1179 
1180  REAL(kind=dp) :: a, at2, h, kf, kl, ks, phi, phi3, r, t2, &
1181  zeta
1182 
1183  r = ra + rb
1184  h = 0.0_dp
1185  IF (r > rc .AND. ngr > ngrc) THEN
1186  zeta = (ra - rb)/r
1187  IF (zeta > 1.0_dp) zeta = 1.0_dp ! machine precision problem
1188  IF (zeta < -1.0_dp) zeta = -1.0_dp ! machine precision problem
1189  phi = ((1.0_dp + zeta)**tb3 + (1.0_dp - zeta)**tb3)/2.0_dp
1190  phi3 = phi*phi*phi
1191  kf = (3.0_dp*r*pi*pi)**ob3
1192  ks = sqrt(4.0_dp*kf/pi)
1193  t2 = (ngr/(2.0_dp*phi*ks*r))**2
1194  a = beta_ec/gamma_saop/(exp(-ec_unif/(gamma_saop*phi3)) - 1.0_dp)
1195  at2 = a*t2
1196  kl = (1.0_dp + at2)/(1.0_dp + at2 + at2*at2)
1197  h = gamma_saop*log(1.0_dp + beta_ec/gamma_saop*t2*kl)
1198  END IF
1199  res = ec_unif + h
1200 
1201  END FUNCTION calc_ecpbe_u
1202 
1203 ! **************************************************************************************************
1204 !> \brief ...
1205 !> \param r ...
1206 !> \param ngr ...
1207 !> \param ec_unif ...
1208 !> \param rc ...
1209 !> \param ngrc ...
1210 !> \return ...
1211 ! **************************************************************************************************
1212  FUNCTION calc_ecpbe_r(r, ngr, ec_unif, rc, ngrc) RESULT(res)
1213 
1214  REAL(kind=dp), INTENT(in) :: r, ngr, ec_unif, rc, ngrc
1215  REAL(kind=dp) :: res
1216 
1217  REAL(kind=dp) :: a, at2, h, kf, kl, ks, t2
1218 
1219  h = 0.0_dp
1220  IF (r > rc .AND. ngr > ngrc) THEN
1221  kf = (3.0_dp*r*pi*pi)**(1.0_dp/3.0_dp)
1222  ks = sqrt(4.0_dp*kf/pi)
1223  t2 = (ngr/(2.0_dp*ks*r))**2
1224  a = beta_ec/gamma_saop/(exp(-ec_unif/gamma_saop) - 1.0_dp)
1225  at2 = a*t2
1226  kl = (1.0_dp + at2)/(1.0_dp + at2 + at2*at2)
1227  h = gamma_saop*log(1.0_dp + beta_ec/gamma_saop*t2*kl)
1228  END IF
1229  res = ec_unif + h
1230 
1231  END FUNCTION calc_ecpbe_r
1232 
1233 ! **************************************************************************************************
1234 !> \brief ...
1235 !> \param ra ...
1236 !> \param rb ...
1237 !> \param ngr ...
1238 !> \param rc ...
1239 !> \param ngrc ...
1240 !> \return ...
1241 ! **************************************************************************************************
1242  FUNCTION calc_expbe_u(ra, rb, ngr, rc, ngrc) RESULT(res)
1243 
1244  REAL(kind=dp), INTENT(in) :: ra, rb, ngr, rc, ngrc
1245  REAL(kind=dp) :: res
1246 
1247  REAL(kind=dp) :: r
1248 
1249  r = ra + rb
1250  res = calc_expbe_r(r, ngr, rc, ngrc)
1251 
1252  END FUNCTION calc_expbe_u
1253 
1254 ! **************************************************************************************************
1255 !> \brief ...
1256 !> \param r ...
1257 !> \param ngr ...
1258 !> \param rc ...
1259 !> \param ngrc ...
1260 !> \return ...
1261 ! **************************************************************************************************
1262  FUNCTION calc_expbe_r(r, ngr, rc, ngrc) RESULT(res)
1263 
1264  REAL(kind=dp), INTENT(in) :: r, ngr, rc, ngrc
1265  REAL(kind=dp) :: res
1266 
1267  REAL(kind=dp) :: ex_unif, fx, kf, s
1268 
1269  IF (r > rc) THEN
1270  kf = (3.0_dp*r*pi*pi)**(1.0_dp/3.0_dp)
1271  ex_unif = -3.0_dp*kf/(4.0_dp*pi)
1272  fx = 1.0_dp
1273  IF (ngr > ngrc) THEN
1274  s = ngr/(2.0_dp*kf*r)
1275  fx = fx + kappa - kappa/(1.0_dp + mu*s*s/kappa)
1276  END IF
1277  res = ex_unif*fx
1278  ELSE
1279  res = 0.0_dp
1280  END IF
1281 
1282  END FUNCTION calc_expbe_r
1283 
1284 END MODULE xc_pot_saop
Define the atomic kind types and their sub types.
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.
various utilities that regard array of different kinds: output, allocation,... maybe it is not a good...
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
DBCSR operations in CP2K.
subroutine, public cp_dbcsr_plus_fm_fm_t(sparse_matrix, matrix_v, matrix_g, ncol, alpha, keep_sparsity, symmetry_mode)
performs the multiplication sparse_matrix+dense_mat*dens_mat^T if matrix_g is not explicitly given,...
represent a full matrix distributed on many processors
Definition: cp_fm_types.F:15
subroutine, public cp_fm_get_info(matrix, name, nrow_global, ncol_global, nrow_block, ncol_block, nrow_local, ncol_local, row_indices, col_indices, local_data, context, nrow_locals, ncol_locals, matrix_struct, para_env)
returns all kind of information about the full matrix
Definition: cp_fm_types.F:1016
subroutine, public cp_fm_set_submatrix(fm, new_values, start_row, start_col, n_rows, n_cols, alpha, beta, transpose)
sets a submatrix of a full matrix fm(start_row:start_row+n_rows,start_col:start_col+n_cols) = alpha*o...
Definition: cp_fm_types.F:768
subroutine, public cp_fm_set_all(matrix, alpha, beta)
set all elements of a matrix to the same value, and optionally the diagonal to a different one
Definition: cp_fm_types.F:535
subroutine, public cp_fm_get_submatrix(fm, target_m, start_row, start_col, n_rows, n_cols, transpose)
gets a submatrix of a full matrix op(target_m)(1:n_rows,1:n_cols) =fm(start_row:start_row+n_rows,...
Definition: cp_fm_types.F:901
subroutine, public cp_fm_create(matrix, matrix_struct, name, use_sp)
creates a new full matrix with the given structure
Definition: cp_fm_types.F:167
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public oe_saop
integer, parameter, public xc_funct_no_shortcut
integer, parameter, public do_method_gapw
integer, parameter, public oe_lb
integer, parameter, public oe_gllb
objects that represent the structure of input sections and the data contained in an input section
recursive subroutine, public section_vals_create(section_vals, section)
creates a object where to store the values of a section
subroutine, public section_vals_val_set(section_vals, keyword_name, i_rep_section, i_rep_val, val, l_val, i_val, r_val, c_val, l_vals_ptr, i_vals_ptr, r_vals_ptr, c_vals_ptr)
sets the requested value
subroutine, public section_vals_retain(section_vals)
retains the given section values (see doc/ReferenceCounting.html)
recursive type(section_vals_type) function, pointer, public section_vals_get_subs_vals(section_vals, subsection_name, i_rep_section, can_return_null)
returns the values of the requested subsection
subroutine, public section_vals_set_subs_vals(section_vals, subsection_name, new_section_vals, i_rep_section)
replaces of the requested subsection with the one given
subroutine, public section_vals_duplicate(section_vals_in, section_vals_out, i_rep_start, i_rep_end)
creates a deep copy from section_vals_in to section_vals_out
subroutine, public section_vals_val_get(section_vals, keyword_name, i_rep_section, i_rep_val, n_rep_val, val, l_val, i_val, r_val, c_val, l_vals, i_vals, r_vals, c_vals, explicit)
returns the requested value
recursive subroutine, public section_vals_release(section_vals)
releases the given object
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
real(kind=dp), parameter, public pi
Interface to the message passing library MPI.
container for various plainwaves related things
Definition: pw_env_types.F:14
subroutine, public pw_env_get(pw_env, pw_pools, cube_info, gridlevel_info, auxbas_pw_pool, auxbas_grid, auxbas_rs_desc, auxbas_rs_grid, rs_descs, rs_grids, xc_pw_pool, vdw_pw_pool, poisson_env, interp_section)
returns the various attributes of the pw env
Definition: pw_env_types.F:113
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
Calculate the plane wave density by collocating the primitive Gaussian functions (pgf).
subroutine, public calculate_rho_elec(matrix_p, matrix_p_kp, rho, rho_gspace, total_rho, ks_env, soft_valid, compute_tau, compute_grad, basis_type, der_type, idir, task_list_external, pw_env_external)
computes the density corresponding to a given density matrix on the grid
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, 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_nonbond, sab_almo, sab_kp, sab_kp_nosym, 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, 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, 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, ecoul_1c, rho0_s_rs, rho0_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, 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, rhs)
Get the QUICKSTEP environment.
subroutine, public prepare_gapw_den(qs_env, local_rho_set, do_rho0, kind_set_external)
...
Integrate single or product functions over a potential on a RS grid.
Define the quickstep kind type and their sub types.
Definition: qs_kind_types.F:23
subroutine, public get_qs_kind(qs_kind, basis_set, basis_type, ncgf, nsgf, all_potential, tnadd_potential, gth_potential, sgp_potential, upf_potential, se_parameter, dftb_parameter, xtb_parameter, dftb3_param, 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_r3d_rs_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, 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, init_u_ramping_each_scf, reltmat, ghost, floating, name, element_symbol, pao_basis_size, pao_potentials, pao_descriptors, nelec)
Get attributes of an atomic kind.
routines that build the Kohn-Sham matrix contributions coming from local atomic densities
Definition: qs_ks_atom.F:12
subroutine, public update_ks_atom(qs_env, ksmat, pmat, forces, tddft, rho_atom_external, kind_set_external, oce_external, sab_external, kscale, kintegral, kforce, fscale)
The correction to the KS matrix due to the GAPW local terms to the hartree and XC contributions is he...
Definition: qs_ks_atom.F:110
subroutine, public local_rho_set_create(local_rho_set)
...
subroutine, public local_rho_set_release(local_rho_set)
...
Definition and initialisation of the mo data type.
Definition: qs_mo_types.F:22
subroutine, public get_mo_set(mo_set, maxocc, homo, lfomo, nao, nelectron, n_el_f, nmo, eigenvalues, occupation_numbers, mo_coeff, mo_coeff_b, uniform_occupation, kTS, mu, flexible_electron_count)
Get the components of a MO set data structure.
Definition: qs_mo_types.F:397
Define the neighbor list data types and the corresponding functionality.
subroutine, public allocate_rho_atom_internals(rho_atom_set, atomic_kind_set, qs_kind_set, dft_control, para_env)
...
subroutine, public calculate_rho_atom_coeff(qs_env, rho_ao, rho_atom_set, qs_kind_set, oce, sab, para_env)
...
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)
...
superstucture that hold various representations of the density and keeps track of which ones are vali...
Definition: qs_rho_types.F:18
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
routines that build the integrals of the Vxc potential calculated for the atomic density in the basis...
Definition: qs_vxc_atom.F:12
subroutine, public calc_rho_angular(grid_atom, harmonics, nspins, grad_func, ir, r_h, r_s, rho_h, rho_s, dr_h, dr_s, r_h_d, r_s_d, drho_h, drho_s)
...
Definition: qs_vxc_atom.F:1537
subroutine, public gavxcgb_nogc(vxc_h, vxc_s, int_hh, int_ss, grid_atom, basis_1c, harmonics, nspins)
...
Definition: qs_vxc_atom.F:1954
All kind of helpful little routines.
Definition: util.F:14
pure integer function, dimension(2), public get_limit(m, n, me)
divide m entries into n parts, return size of part me
Definition: util.F:333
Definition: xc_atom.F:9
subroutine, public xc_rho_set_atom_update(rho_set, needs, nspins, bo)
...
Definition: xc_atom.F:432
subroutine, public fill_rho_set(rho_set, lsd, nspins, needs, rho, drho, tau, na, ir)
...
Definition: xc_atom.F:592
subroutine, public vxc_of_r_new(xc_fun_section, rho_set, deriv_set, deriv_order, needs, w, lsd, na, nr, exc, vxc, vxg, vtau, energy_only, epr_xc, adiabatic_rescale_factor)
...
Definition: xc_atom.F:65
represent a group ofunctional derivatives
subroutine, public xc_dset_zero_all(deriv_set)
...
type(xc_derivative_type) function, pointer, public xc_dset_get_derivative(derivative_set, description, allocate_deriv)
returns the requested xc_derivative
subroutine, public xc_dset_release(derivative_set)
releases a derivative set
subroutine, public xc_dset_create(derivative_set, pw_pool, local_bounds)
creates a derivative set object
Provides types for the management of the xc-functionals and their derivatives.
subroutine, public xc_derivative_get(deriv, split_desc, order, deriv_data, accept_null_data)
returns various information on the given derivative
subroutine, public xc_functionals_eval(functionals, lsd, rho_set, deriv_set, deriv_order)
...
Calculate the saop potential.
Definition: xc_pot_saop.F:11
subroutine, public add_saop_pot(ks_matrix, qs_env, oe_corr)
...
Definition: xc_pot_saop.F:126
contains the structure
elemental subroutine, public xc_rho_cflags_setall(cflags, value)
sets all the flags to the given value
contains the structure
subroutine, public xc_rho_set_update(rho_set, rho_r, rho_g, tau, needs, xc_deriv_method_id, xc_rho_smooth_id, pw_pool)
updates the given rho set with the density given by rho_r (and rho_g). The rho set will contain the c...
subroutine, public xc_rho_set_create(rho_set, local_bounds, rho_cutoff, drho_cutoff, tau_cutoff)
allocates and does (minimal) initialization of a rho_set
subroutine, public xc_rho_set_release(rho_set, pw_pool)
releases the given rho_set
calculates the Becke 88 exchange functional
Definition: xc_xbecke88.F:14
subroutine, public xb88_lsd_info(reference, shortform, needs, max_deriv)
return various information on the functional
Definition: xc_xbecke88.F:90
subroutine, public xb88_lda_info(reference, shortform, needs, max_deriv)
return various information on the functional
Definition: xc_xbecke88.F:59
Exchange and Correlation functional calculations.
Definition: xc.F:17
subroutine, public xc_vxc_pw_create1(vxc_rho, vxc_tau, rho_r, rho_g, tau, exc, xc_section, pw_pool, compute_virial, virial_xc, exc_r)
Exchange and Correlation functional calculations. depending on the selected functional_routine calls ...
Definition: xc.F:150
Orbital angular momentum.
Definition: grid_common.h:125