Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2025 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : !> \brief Calculates electric field gradients
10 : !> H.M. Petrili, P.E. Blochl, P. Blaha, K. Schwarz, PRB 57, 14690 (1998)
11 : !> \par History
12 : !> 12.2007 Added checksum for interpolation regtest [rdeclerck]
13 : !> \author JGH (03-05-2006)
14 : ! **************************************************************************************************
15 : MODULE qs_electric_field_gradient
16 : USE atomic_kind_types, ONLY: atomic_kind_type,&
17 : get_atomic_kind
18 : USE basis_set_types, ONLY: get_gto_basis_set,&
19 : gto_basis_set_type
20 : USE cp_control_types, ONLY: dft_control_type
21 : USE cp_log_handling, ONLY: cp_get_default_logger,&
22 : cp_logger_type
23 : USE cp_output_handling, ONLY: cp_print_key_unit_nr
24 : USE eigenvalueproblems, ONLY: diagonalise
25 : USE input_section_types, ONLY: section_get_lval,&
26 : section_vals_get_subs_vals,&
27 : section_vals_type,&
28 : section_vals_val_get
29 : USE kinds, ONLY: dp
30 : USE mathconstants, ONLY: fac,&
31 : fourpi
32 : USE message_passing, ONLY: mp_para_env_type
33 : USE orbital_pointers, ONLY: indso,&
34 : nsoset
35 : USE particle_types, ONLY: particle_type
36 : USE paw_basis_types, ONLY: get_paw_basis_info
37 : USE physcon, ONLY: a_bohr,&
38 : e_charge,&
39 : joule
40 : USE pw_env_types, ONLY: pw_env_get,&
41 : pw_env_type
42 : USE pw_methods, ONLY: pw_dr2,&
43 : pw_integral_ab,&
44 : pw_smoothing,&
45 : pw_structure_factor,&
46 : pw_transfer
47 : USE pw_poisson_methods, ONLY: pw_poisson_solve
48 : USE pw_poisson_types, ONLY: pw_poisson_type
49 : USE pw_pool_types, ONLY: pw_pool_type
50 : USE pw_spline_utils, ONLY: &
51 : Eval_Interp_Spl3_pbc, Eval_d_Interp_Spl3_pbc, find_coeffs, pw_spline_do_precond, &
52 : pw_spline_precond_create, pw_spline_precond_release, pw_spline_precond_set_kind, &
53 : pw_spline_precond_type, spl3_pbc
54 : USE pw_types, ONLY: pw_c1d_gs_type,&
55 : pw_r3d_rs_type
56 : USE qs_environment_types, ONLY: get_qs_env,&
57 : qs_environment_type
58 : USE qs_gapw_densities, ONLY: prepare_gapw_den
59 : USE qs_harmonics_atom, ONLY: get_none0_cg_list,&
60 : harmonics_atom_type
61 : USE qs_kind_types, ONLY: get_qs_kind,&
62 : qs_kind_type
63 : USE qs_ks_methods, ONLY: calc_rho_tot_gspace
64 : USE qs_rho_atom_types, ONLY: rho_atom_type
65 : USE qs_rho_types, ONLY: qs_rho_type
66 : USE util, ONLY: get_limit,&
67 : sort
68 : #include "./base/base_uses.f90"
69 :
70 : IMPLICIT NONE
71 :
72 : PRIVATE
73 : PUBLIC :: qs_efg_calc
74 :
75 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'qs_electric_field_gradient'
76 :
77 : CONTAINS
78 :
79 : ! **************************************************************************************************
80 : !> \brief ...
81 : !> \param qs_env ...
82 : ! **************************************************************************************************
83 30 : SUBROUTINE qs_efg_calc(qs_env)
84 :
85 : TYPE(qs_environment_type), POINTER :: qs_env
86 :
87 : CHARACTER(LEN=*), PARAMETER :: routineN = 'qs_efg_calc'
88 :
89 : CHARACTER(LEN=2) :: element_symbol
90 : INTEGER :: aint_precond, handle, i, iat, iatom, ij, &
91 : ikind, j, max_iter, natom, natomkind, &
92 : nkind, nspins, precond_kind, unit_nr
93 30 : INTEGER, DIMENSION(:), POINTER :: atom_list
94 : LOGICAL :: efg_debug, efg_interpolation, gapw, &
95 : paw_atom, smoothing, success
96 : REAL(KIND=dp) :: chk_spl, ecut, efg_units, efg_val, &
97 : ehartree, eps_r, eps_x, f1, f2, sigma
98 60 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: efg_diagval, vh0
99 30 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: efg_pw, efg_tensor
100 : REAL(KIND=dp), DIMENSION(3) :: eigenvalues, ra
101 : REAL(KIND=dp), DIMENSION(3, 3) :: eigenvectors
102 30 : REAL(KIND=dp), DIMENSION(:), POINTER :: rvals
103 30 : TYPE(atomic_kind_type), DIMENSION(:), POINTER :: atomic_kind_set
104 : TYPE(cp_logger_type), POINTER :: logger
105 : TYPE(dft_control_type), POINTER :: dft_control
106 : TYPE(mp_para_env_type), POINTER :: para_env
107 30 : TYPE(particle_type), DIMENSION(:), POINTER :: particle_set
108 : TYPE(pw_c1d_gs_type) :: rho_tot_gspace, structure_factor, &
109 : v_hartree_gspace
110 210 : TYPE(pw_c1d_gs_type), DIMENSION(6) :: dvr2
111 : TYPE(pw_env_type), POINTER :: pw_env
112 : TYPE(pw_poisson_type), POINTER :: poisson_env
113 : TYPE(pw_pool_type), POINTER :: auxbas_pw_pool
114 : TYPE(pw_r3d_rs_type) :: dvr2rs
115 210 : TYPE(pw_r3d_rs_type), DIMENSION(6) :: dvspl
116 : TYPE(pw_spline_precond_type) :: precond
117 30 : TYPE(qs_kind_type), DIMENSION(:), POINTER :: qs_kind_set
118 : TYPE(qs_rho_type), POINTER :: rho
119 30 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
120 : TYPE(section_vals_type), POINTER :: dft_section, input, interp_section
121 :
122 30 : NULLIFY (atomic_kind_set, qs_kind_set, dft_control, para_env, particle_set, rho, &
123 30 : rho_atom_set, input, dft_section, interp_section)
124 :
125 30 : CALL timeset(routineN, handle)
126 :
127 30 : logger => cp_get_default_logger()
128 :
129 30 : chk_spl = 0.0_dp
130 30 : efg_units = Joule/a_bohr**2/e_charge*1.e-21_dp
131 30 : f1 = SQRT(15._dp/fourpi)
132 30 : f2 = SQRT(5._dp/fourpi)
133 :
134 : CALL get_qs_env(qs_env=qs_env, dft_control=dft_control, &
135 : rho=rho, qs_kind_set=qs_kind_set, &
136 : atomic_kind_set=atomic_kind_set, &
137 : rho_atom_set=rho_atom_set, pw_env=pw_env, &
138 : particle_set=particle_set, para_env=para_env, &
139 30 : input=input)
140 :
141 30 : dft_section => section_vals_get_subs_vals(input, "DFT")
142 :
143 : efg_interpolation = section_get_lval(section_vals=dft_section, &
144 30 : keyword_name="PRINT%ELECTRIC_FIELD_GRADIENT%INTERPOLATION")
145 : efg_debug = section_get_lval(section_vals=dft_section, &
146 30 : keyword_name="PRINT%ELECTRIC_FIELD_GRADIENT%DEBUG")
147 : CALL section_vals_val_get(dft_section, &
148 : "PRINT%ELECTRIC_FIELD_GRADIENT%GSPACE_SMOOTHING", &
149 30 : r_vals=rvals)
150 30 : ecut = rvals(1)
151 30 : sigma = rvals(2)
152 30 : IF (ecut == 0._dp .AND. sigma <= 0._dp) THEN
153 0 : smoothing = .FALSE.
154 0 : ecut = 1.e10_dp ! not used, just to have vars defined
155 0 : sigma = 1._dp ! not used, just to have vars defined
156 30 : ELSEIF (ecut == -1._dp .AND. sigma == -1._dp) THEN
157 30 : smoothing = .TRUE.
158 30 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool)
159 30 : CALL auxbas_pw_pool%create_pw(dvr2rs)
160 30 : ecut = 2._dp*dvr2rs%pw_grid%cutoff*0.875_dp
161 30 : sigma = 2._dp*dvr2rs%pw_grid%cutoff*0.125_dp
162 30 : CALL auxbas_pw_pool%give_back_pw(dvr2rs)
163 : ELSE
164 : smoothing = .TRUE.
165 : END IF
166 30 : CPASSERT(ecut > 0._dp)
167 30 : CPASSERT(sigma > 0._dp)
168 :
169 : unit_nr = cp_print_key_unit_nr(logger, dft_section, "PRINT%ELECTRIC_FIELD_GRADIENT", &
170 30 : extension=".efg", log_filename=.FALSE.)
171 :
172 30 : IF (unit_nr > 0) THEN
173 15 : WRITE (unit_nr, "(/,A,/)") " ELECTRIC FIELD GRADIENTS [10**21 V/m^2]"
174 15 : IF (efg_interpolation) THEN
175 : WRITE (unit_nr, "(T16,A)") &
176 1 : " Smooth potential contribution calculated by spline interpolation"
177 : ELSE
178 : WRITE (unit_nr, "(T12,A)") &
179 14 : " Smooth potential contribution calculated by plane wave interpolation"
180 : END IF
181 15 : IF (smoothing) THEN
182 : WRITE (unit_nr, "(T36,A)") &
183 15 : " G-Space potential smoothed by Fermi function"
184 15 : WRITE (unit_nr, "(T36,A,T71,F10.4)") " Cutoff (eV) ", ecut
185 15 : WRITE (unit_nr, "(T36,A,T71,F10.4)") " Width (eV) ", sigma
186 : END IF
187 15 : WRITE (unit_nr, *)
188 : END IF
189 :
190 30 : gapw = dft_control%qs_control%gapw
191 30 : nspins = dft_control%nspins
192 :
193 30 : natom = SIZE(particle_set, 1)
194 90 : ALLOCATE (efg_tensor(3, 3, natom))
195 1096 : efg_tensor = 0._dp
196 30 : IF (efg_debug) THEN
197 4 : ALLOCATE (efg_pw(3, 3, natom))
198 80 : efg_pw = 0._dp
199 : END IF
200 90 : ALLOCATE (efg_diagval(3, natom))
201 358 : efg_diagval = 0._dp
202 :
203 90 : ALLOCATE (vh0(1:natom, -2:2))
204 590 : vh0 = 0._dp
205 :
206 : !prepare calculation
207 : CALL pw_env_get(pw_env, auxbas_pw_pool=auxbas_pw_pool, &
208 30 : poisson_env=poisson_env)
209 30 : IF (gapw) CALL prepare_gapw_den(qs_env, do_rho0=.TRUE.)
210 :
211 : !calculate electrostatic potential
212 30 : CALL auxbas_pw_pool%create_pw(v_hartree_gspace)
213 30 : CALL auxbas_pw_pool%create_pw(rho_tot_gspace)
214 30 : CALL calc_rho_tot_gspace(rho_tot_gspace, qs_env, rho)
215 :
216 : CALL pw_poisson_solve(poisson_env, rho_tot_gspace, ehartree, &
217 30 : v_hartree_gspace)
218 30 : CALL auxbas_pw_pool%give_back_pw(rho_tot_gspace)
219 :
220 : ! smoothing of potential
221 30 : IF (smoothing) CALL pw_smoothing(v_hartree_gspace, ecut, sigma)
222 :
223 120 : DO i = 1, 3
224 300 : DO j = 1, i
225 180 : ij = (i*(i - 1))/2 + j
226 180 : CALL auxbas_pw_pool%create_pw(dvr2(ij))
227 270 : CALL pw_dr2(v_hartree_gspace, dvr2(ij), i, j)
228 : END DO
229 : END DO
230 :
231 30 : IF (.NOT. efg_interpolation) THEN
232 28 : CALL auxbas_pw_pool%create_pw(structure_factor)
233 : ELSE
234 :
235 : interp_section => section_vals_get_subs_vals(dft_section, &
236 2 : "PRINT%ELECTRIC_FIELD_GRADIENT%INTERPOLATOR")
237 : CALL section_vals_val_get(interp_section, "aint_precond", &
238 2 : i_val=aint_precond)
239 2 : CALL section_vals_val_get(interp_section, "precond", i_val=precond_kind)
240 2 : CALL section_vals_val_get(interp_section, "max_iter", i_val=max_iter)
241 2 : CALL section_vals_val_get(interp_section, "eps_r", r_val=eps_r)
242 2 : CALL section_vals_val_get(interp_section, "eps_x", r_val=eps_x)
243 :
244 2 : CALL auxbas_pw_pool%create_pw(dvr2rs)
245 14 : DO i = 1, 6
246 12 : CALL auxbas_pw_pool%create_pw(dvspl(i))
247 12 : CALL pw_transfer(dvr2(i), dvr2rs)
248 : ! calculate spline coefficients
249 : CALL pw_spline_precond_create(precond, precond_kind=aint_precond, &
250 12 : pool=auxbas_pw_pool, pbc=.TRUE., transpose=.FALSE.)
251 12 : CALL pw_spline_do_precond(precond, dvr2rs, dvspl(i))
252 12 : CALL pw_spline_precond_set_kind(precond, precond_kind)
253 : success = find_coeffs(values=dvr2rs, coeffs=dvspl(i), &
254 : linOp=spl3_pbc, preconditioner=precond, pool=auxbas_pw_pool, &
255 12 : eps_r=eps_r, eps_x=eps_x, max_iter=max_iter)
256 12 : CPASSERT(success)
257 12 : CALL pw_spline_precond_release(precond)
258 14 : CALL auxbas_pw_pool%give_back_pw(dvr2(i))
259 : END DO
260 2 : CALL auxbas_pw_pool%give_back_pw(dvr2rs)
261 : END IF
262 :
263 30 : nkind = SIZE(qs_kind_set)
264 :
265 86 : DO ikind = 1, nkind
266 56 : NULLIFY (atom_list)
267 56 : CALL get_atomic_kind(atomic_kind_set(ikind), atom_list=atom_list, natom=natomkind)
268 56 : CALL get_qs_kind(qs_kind_set(ikind), paw_atom=paw_atom)
269 138 : DO iat = 1, natomkind
270 82 : iatom = atom_list(iat)
271 328 : ra = particle_set(iatom)%r
272 82 : IF (efg_interpolation) THEN
273 24 : DO i = 1, 3
274 60 : DO j = 1, i
275 36 : ij = (i*(i - 1))/2 + j
276 36 : efg_val = Eval_Interp_Spl3_pbc(ra, dvspl(ij))
277 36 : efg_tensor(i, j, iatom) = -efg_val
278 36 : efg_tensor(j, i, iatom) = efg_tensor(i, j, iatom)
279 54 : IF (efg_debug) THEN
280 : chk_spl = chk_spl + efg_val + &
281 144 : SUM(Eval_d_Interp_Spl3_pbc(ra, dvspl(ij)))
282 : END IF
283 : END DO
284 : END DO
285 : ELSE
286 76 : CALL pw_structure_factor(structure_factor, ra)
287 304 : DO i = 1, 3
288 760 : DO j = 1, i
289 456 : ij = (i*(i - 1))/2 + j
290 456 : efg_tensor(i, j, iatom) = -pw_integral_ab(dvr2(ij), structure_factor)
291 684 : efg_tensor(j, i, iatom) = efg_tensor(i, j, iatom)
292 : END DO
293 : END DO
294 988 : efg_tensor(:, :, iatom) = efg_tensor(:, :, iatom)/structure_factor%pw_grid%vol
295 : END IF
296 138 : IF (efg_debug) THEN
297 78 : efg_pw(:, :, iatom) = efg_tensor(:, :, iatom)
298 : END IF
299 : END DO
300 :
301 56 : IF (paw_atom) THEN
302 : CALL vlimit_atom(para_env, vh0, rho_atom_set, qs_kind_set(ikind), &
303 0 : atom_list, natomkind, nspins)
304 0 : DO iat = 1, natomkind
305 0 : iatom = atom_list(iat)
306 : efg_tensor(1, 1, iatom) = efg_tensor(1, 1, iatom) &
307 0 : + f1*(vh0(iatom, 2)) - f2*(vh0(iatom, 0))
308 : efg_tensor(2, 2, iatom) = efg_tensor(2, 2, iatom) &
309 0 : - f1*(vh0(iatom, 2)) - f2*(vh0(iatom, 0))
310 0 : efg_tensor(3, 3, iatom) = efg_tensor(3, 3, iatom) + 2._dp*f2*(vh0(iatom, 0))
311 0 : efg_tensor(1, 2, iatom) = efg_tensor(1, 2, iatom) + f1*(vh0(iatom, -2))
312 0 : efg_tensor(2, 1, iatom) = efg_tensor(2, 1, iatom) + f1*(vh0(iatom, -2))
313 0 : efg_tensor(1, 3, iatom) = efg_tensor(1, 3, iatom) + f1*(vh0(iatom, 1))
314 0 : efg_tensor(3, 1, iatom) = efg_tensor(3, 1, iatom) + f1*(vh0(iatom, 1))
315 0 : efg_tensor(2, 3, iatom) = efg_tensor(2, 3, iatom) + f1*(vh0(iatom, -1))
316 0 : efg_tensor(3, 2, iatom) = efg_tensor(3, 2, iatom) + f1*(vh0(iatom, -1))
317 : END DO
318 : END IF
319 :
320 224 : DO iat = 1, natomkind
321 82 : iatom = atom_list(iat)
322 : CALL diagonalise(efg_tensor(:, :, iatom), 3, "U", &
323 82 : eigenvalues, eigenvectors)
324 138 : CALL efgsort(eigenvalues, efg_diagval(:, iatom))
325 : END DO
326 : END DO ! ikind
327 :
328 1096 : efg_tensor(:, :, :) = efg_tensor(:, :, :)*efg_units
329 358 : efg_diagval(:, :) = efg_diagval(:, :)*efg_units
330 :
331 30 : IF (efg_debug) THEN
332 80 : efg_pw(:, :, :) = efg_pw(:, :, :)*efg_units
333 8 : DO iatom = 1, natom
334 8 : IF (unit_nr > 0) THEN
335 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, &
336 3 : element_symbol=element_symbol)
337 : WRITE (UNIT=unit_nr, FMT="(T2,I5,T8,A,T12,A,T15,6F11.5)") &
338 30 : iatom, element_symbol, "PW", ((efg_pw(i, j, iatom), i=1, j), j=1, 3)
339 : WRITE (UNIT=unit_nr, FMT="(T12,A,T15,6F11.5)") &
340 30 : "AT", ((efg_tensor(i, j, iatom) - efg_pw(i, j, iatom), i=1, j), j=1, 3)
341 : END IF
342 : END DO
343 2 : IF (unit_nr > 0) THEN
344 1 : WRITE (UNIT=unit_nr, FMT=*)
345 : END IF
346 2 : IF (efg_interpolation) THEN
347 2 : IF (unit_nr > 0) THEN
348 1 : WRITE (UNIT=unit_nr, FMT="(T2,A,E24.16)") "CheckSum splines =", &
349 2 : chk_spl
350 1 : WRITE (UNIT=unit_nr, FMT=*)
351 : END IF
352 : END IF
353 : END IF
354 :
355 112 : DO iatom = 1, natom
356 112 : IF (unit_nr > 0) THEN
357 : CALL get_atomic_kind(particle_set(iatom)%atomic_kind, &
358 41 : element_symbol=element_symbol)
359 : WRITE (UNIT=unit_nr, FMT="(T2,I5,T8,A,T12,A,3(T39,3F14.7,/))") &
360 164 : iatom, element_symbol, "EFG Tensor", (efg_tensor(i, :, iatom), i=1, 3)
361 : WRITE (UNIT=unit_nr, FMT="(T12,A,T39,3F14.7)") &
362 41 : "EFG Tensor eigenvalues", efg_diagval(:, iatom)
363 41 : WRITE (UNIT=unit_nr, FMT="(T12,A,T67,F14.7)") "EFG Tensor anisotropy", &
364 82 : (efg_diagval(1, iatom) - efg_diagval(2, iatom))/efg_diagval(3, iatom)
365 41 : WRITE (UNIT=unit_nr, FMT=*)
366 : END IF
367 : END DO
368 :
369 30 : CALL auxbas_pw_pool%give_back_pw(v_hartree_gspace)
370 30 : IF (.NOT. efg_interpolation) THEN
371 28 : CALL auxbas_pw_pool%give_back_pw(structure_factor)
372 196 : DO i = 1, 6
373 196 : CALL auxbas_pw_pool%give_back_pw(dvr2(i))
374 : END DO
375 : ELSE
376 14 : DO i = 1, 6
377 14 : CALL auxbas_pw_pool%give_back_pw(dvspl(i))
378 : END DO
379 : END IF
380 :
381 30 : DEALLOCATE (efg_tensor)
382 30 : IF (efg_debug) THEN
383 2 : DEALLOCATE (efg_pw)
384 : END IF
385 :
386 30 : DEALLOCATE (vh0)
387 :
388 30 : CALL timestop(handle)
389 :
390 390 : END SUBROUTINE qs_efg_calc
391 :
392 : ! **************************************************************************************************
393 : !> \brief ...
394 : !> \param para_env ...
395 : !> \param vlimit ...
396 : !> \param rho_atom_set ...
397 : !> \param qs_kind ...
398 : !> \param atom_list ...
399 : !> \param natom ...
400 : !> \param nspins ...
401 : ! **************************************************************************************************
402 0 : SUBROUTINE vlimit_atom(para_env, vlimit, rho_atom_set, qs_kind, &
403 0 : atom_list, natom, nspins)
404 :
405 : ! calculate : Limit(r->0) V_hartree(r)/r^2
406 :
407 : TYPE(mp_para_env_type), POINTER :: para_env
408 : REAL(dp), DIMENSION(:, -2:), INTENT(inout) :: vlimit
409 : TYPE(rho_atom_type), DIMENSION(:), POINTER :: rho_atom_set
410 : TYPE(qs_kind_type), INTENT(IN) :: qs_kind
411 : INTEGER, DIMENSION(:), INTENT(IN) :: atom_list
412 : INTEGER, INTENT(IN) :: natom, nspins
413 :
414 : INTEGER :: i, i1, i2, iat, iatom, icg, ipgf1, ipgf2, iset1, iset2, iso, iso1, iso1_first, &
415 : iso1_last, iso2, iso2_first, iso2_last, l, l_iso, llmax, m1s, m2s, m_iso, max_iso_not0, &
416 : max_iso_not0_local, max_s_harm, maxl, maxso, mepos, n1s, n2s, nset, num_pe, size1, size2
417 : INTEGER, ALLOCATABLE, DIMENSION(:) :: cg_n_list
418 : INTEGER, ALLOCATABLE, DIMENSION(:, :, :) :: cg_list
419 : INTEGER, DIMENSION(2) :: bo
420 0 : INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf, o2nindex
421 : REAL(dp) :: zet12
422 0 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: CPC_sphere
423 : REAL(dp), DIMENSION(20) :: vgg
424 0 : REAL(dp), DIMENSION(:, :), POINTER :: coeff, zet
425 : REAL(dp), DIMENSION(:, :, :), POINTER :: my_CG
426 : TYPE(gto_basis_set_type), POINTER :: basis_1c
427 : TYPE(harmonics_atom_type), POINTER :: harmonics
428 :
429 0 : NULLIFY (basis_1c)
430 0 : NULLIFY (harmonics)
431 0 : NULLIFY (lmin, lmax, npgf, zet, my_CG, coeff)
432 :
433 : CALL get_qs_kind(qs_kind=qs_kind, basis_set=basis_1c, basis_type="GAPW_1C", &
434 0 : harmonics=harmonics)
435 :
436 : CALL get_gto_basis_set(gto_basis_set=basis_1c, lmax=lmax, lmin=lmin, &
437 : maxl=maxl, npgf=npgf, nset=nset, zet=zet, &
438 0 : maxso=maxso)
439 0 : CALL get_paw_basis_info(basis_1c, o2nindex=o2nindex)
440 :
441 0 : max_iso_not0 = harmonics%max_iso_not0
442 0 : max_s_harm = harmonics%max_s_harm
443 0 : llmax = harmonics%llmax
444 :
445 : ! Distribute the atoms of this kind
446 0 : num_pe = para_env%num_pe
447 0 : mepos = para_env%mepos
448 0 : bo = get_limit(natom, num_pe, mepos)
449 :
450 0 : my_CG => harmonics%my_CG
451 :
452 0 : ALLOCATE (CPC_sphere(nsoset(maxl), nsoset(maxl)))
453 0 : ALLOCATE (cg_list(2, nsoset(maxl)**2, max_s_harm), cg_n_list(max_s_harm))
454 :
455 0 : m1s = 0
456 0 : DO iset1 = 1, nset
457 : m2s = 0
458 0 : DO iset2 = 1, nset
459 : CALL get_none0_cg_list(my_CG, lmin(iset1), lmax(iset1), lmin(iset2), lmax(iset2), &
460 0 : max_s_harm, llmax, cg_list, cg_n_list, max_iso_not0_local)
461 0 : CPASSERT(max_iso_not0_local .LE. max_iso_not0)
462 :
463 0 : n1s = nsoset(lmax(iset1))
464 0 : DO ipgf1 = 1, npgf(iset1)
465 0 : iso1_first = nsoset(lmin(iset1) - 1) + 1 + n1s*(ipgf1 - 1) + m1s
466 0 : iso1_last = nsoset(lmax(iset1)) + n1s*(ipgf1 - 1) + m1s
467 0 : size1 = iso1_last - iso1_first + 1
468 0 : iso1_first = o2nindex(iso1_first)
469 0 : iso1_last = o2nindex(iso1_last)
470 0 : i1 = iso1_last - iso1_first + 1
471 0 : CPASSERT(size1 == i1)
472 0 : i1 = nsoset(lmin(iset1) - 1) + 1
473 :
474 0 : n2s = nsoset(lmax(iset2))
475 0 : DO ipgf2 = 1, npgf(iset2)
476 0 : iso2_first = nsoset(lmin(iset2) - 1) + 1 + n2s*(ipgf2 - 1) + m2s
477 0 : iso2_last = nsoset(lmax(iset2)) + n2s*(ipgf2 - 1) + m2s
478 0 : size2 = iso2_last - iso2_first + 1
479 0 : iso2_first = o2nindex(iso2_first)
480 0 : iso2_last = o2nindex(iso2_last)
481 0 : i2 = iso2_last - iso2_first + 1
482 0 : CPASSERT(size2 == i2)
483 0 : i2 = nsoset(lmin(iset2) - 1) + 1
484 :
485 0 : zet12 = zet(ipgf1, iset1) + zet(ipgf2, iset2)
486 :
487 0 : vgg = 0.0_dp
488 :
489 0 : DO iso = 1, max_iso_not0_local
490 0 : l_iso = indso(1, iso)
491 0 : IF (l_iso /= 2) CYCLE
492 0 : DO icg = 1, cg_n_list(iso)
493 0 : iso1 = cg_list(1, icg, iso)
494 0 : iso2 = cg_list(2, icg, iso)
495 0 : l = indso(1, iso1) + indso(1, iso2)
496 0 : IF (MOD(l, 2) == 0 .AND. l > 0) THEN
497 0 : vgg(l/2) = fourpi/10._dp*fac(l - 2)/zet12**(l/2)
498 : END IF
499 : END DO
500 : END DO
501 :
502 0 : DO iat = bo(1), bo(2)
503 0 : iatom = atom_list(iat)
504 :
505 0 : CPC_sphere = 0.0_dp
506 0 : DO i = 1, nspins
507 0 : coeff => rho_atom_set(iatom)%cpc_h(i)%r_coef
508 : CPC_sphere(i1:i1 + size1 - 1, i2:i2 + size2 - 1) = &
509 : CPC_sphere(i1:i1 + size1 - 1, i2:i2 + size2 - 1) + &
510 0 : coeff(iso1_first:iso1_last, iso2_first:iso2_last)
511 0 : coeff => rho_atom_set(iatom)%cpc_s(i)%r_coef
512 : CPC_sphere(i1:i1 + size1 - 1, i2:i2 + size2 - 1) = &
513 : CPC_sphere(i1:i1 + size1 - 1, i2:i2 + size2 - 1) - &
514 0 : coeff(iso1_first:iso1_last, iso2_first:iso2_last)
515 : END DO ! i
516 :
517 0 : DO iso = 1, max_iso_not0_local
518 0 : l_iso = indso(1, iso)
519 0 : m_iso = indso(1, iso)
520 0 : IF (l_iso /= 2) CYCLE
521 0 : DO icg = 1, cg_n_list(iso)
522 0 : iso1 = cg_list(1, icg, iso)
523 0 : iso2 = cg_list(2, icg, iso)
524 0 : l = indso(1, iso1) + indso(1, iso2)
525 0 : IF (MOD(l, 2) == 0 .AND. l > 0) THEN
526 : vlimit(iatom, m_iso) = vlimit(iatom, m_iso) + &
527 0 : vgg(l/2)*CPC_sphere(iso1, iso2)*my_CG(iso1, iso2, iso)
528 : END IF
529 : END DO ! icg
530 : END DO ! iso
531 :
532 : END DO ! iat
533 :
534 : END DO ! ipgf2
535 : END DO ! ipgf1
536 0 : m2s = m2s + maxso
537 : END DO ! iset2
538 0 : m1s = m1s + maxso
539 : END DO ! iset1
540 :
541 0 : CALL para_env%sum(vlimit)
542 :
543 0 : DEALLOCATE (o2nindex)
544 0 : DEALLOCATE (CPC_sphere)
545 0 : DEALLOCATE (cg_list, cg_n_list)
546 :
547 0 : END SUBROUTINE vlimit_atom
548 :
549 : ! **************************************************************************************************
550 : !> \brief ...
551 : !> \param ein ...
552 : !> \param eout ...
553 : ! **************************************************************************************************
554 82 : SUBROUTINE efgsort(ein, eout)
555 : REAL(dp), DIMENSION(3), INTENT(in) :: ein
556 : REAL(dp), DIMENSION(3), INTENT(inout) :: eout
557 :
558 : INTEGER :: i
559 : INTEGER, DIMENSION(3) :: ind
560 : REAL(dp), DIMENSION(3) :: eab
561 :
562 328 : eab = ABS(ein)
563 82 : CALL sort(eab, 3, ind)
564 328 : DO i = 1, 3
565 328 : eout(i) = ein(ind(i))
566 : END DO
567 :
568 82 : END SUBROUTINE efgsort
569 :
570 : END MODULE qs_electric_field_gradient
571 :
|
