Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 : ! **************************************************************************************************
9 : MODULE atom_sgp
10 : USE atom_set_basis, ONLY: set_kind_basis_atomic
11 : USE atom_types, ONLY: &
12 : atom_basis_gridrep, atom_basis_type, atom_ecppot_type, atom_p_type, atom_type, &
13 : create_opmat, init_atom_basis_default_pp, opmat_type, release_atom_basis, release_opmat
14 : USE atom_utils, ONLY: integrate_grid,&
15 : numpot_matrix
16 : USE basis_set_types, ONLY: gto_basis_set_type
17 : USE input_constants, ONLY: ecp_pseudo,&
18 : gaussian,&
19 : gth_pseudo,&
20 : no_pseudo,&
21 : upf_pseudo
22 : USE input_section_types, ONLY: section_vals_get,&
23 : section_vals_type
24 : USE kahan_sum, ONLY: accurate_dot_product
25 : USE kinds, ONLY: dp
26 : USE mathconstants, ONLY: dfac,&
27 : fourpi,&
28 : rootpi,&
29 : sqrt2
30 : USE mathlib, ONLY: diamat_all,&
31 : get_pseudo_inverse_diag
32 : USE powell, ONLY: opt_state_type,&
33 : powell_optimize
34 : USE upf_types, ONLY: atom_upfpot_type
35 : #include "../base/base_uses.f90"
36 :
37 : IMPLICIT NONE
38 :
39 : TYPE atom_sgp_potential_type
40 : LOGICAL :: has_nonlocal = .FALSE.
41 : INTEGER :: n_nonlocal = 0
42 : INTEGER :: lmax = 0
43 : LOGICAL, DIMENSION(0:5) :: is_nonlocal = .FALSE.
44 : REAL(KIND=dp), DIMENSION(:), POINTER :: a_nonlocal => Null()
45 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: h_nonlocal => Null()
46 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: c_nonlocal => Null()
47 : LOGICAL :: has_local = .FALSE.
48 : INTEGER :: n_local = 0
49 : REAL(KIND=dp) :: zval = 0.0_dp
50 : REAL(KIND=dp) :: ac_local = 0.0_dp
51 : REAL(KIND=dp), DIMENSION(:), POINTER :: a_local => Null()
52 : REAL(KIND=dp), DIMENSION(:), POINTER :: c_local => Null()
53 : LOGICAL :: has_nlcc = .FALSE.
54 : INTEGER :: n_nlcc = 0
55 : REAL(KIND=dp), DIMENSION(:), POINTER :: a_nlcc => Null()
56 : REAL(KIND=dp), DIMENSION(:), POINTER :: c_nlcc => Null()
57 : END TYPE atom_sgp_potential_type
58 :
59 : PRIVATE
60 : PUBLIC :: atom_sgp_potential_type, atom_sgp_release
61 : PUBLIC :: atom_sgp_construction, sgp_construction
62 : PUBLIC :: fit_core_density_gaussians
63 :
64 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atom_sgp'
65 :
66 : ! **************************************************************************************************
67 :
68 : CONTAINS
69 :
70 : ! **************************************************************************************************
71 : !> \brief ...
72 : !> \param sgp_pot ...
73 : !> \param ecp_pot ...
74 : !> \param upf_pot ...
75 : !> \param orb_basis ...
76 : !> \param error ...
77 : ! **************************************************************************************************
78 24 : SUBROUTINE sgp_construction(sgp_pot, ecp_pot, upf_pot, orb_basis, error)
79 :
80 : TYPE(atom_sgp_potential_type) :: sgp_pot
81 : TYPE(atom_ecppot_type), OPTIONAL :: ecp_pot
82 : TYPE(atom_upfpot_type), OPTIONAL :: upf_pot
83 : TYPE(gto_basis_set_type), OPTIONAL, POINTER :: orb_basis
84 : REAL(KIND=dp), DIMENSION(6) :: error
85 :
86 : INTEGER :: i, n
87 : INTEGER, DIMENSION(3) :: mloc
88 : LOGICAL :: is_ecp, is_upf
89 : REAL(KIND=dp) :: errcc, rcut
90 24 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cgauss, cutpots, cutpotu
91 : TYPE(atom_basis_type) :: basis
92 : TYPE(opmat_type), POINTER :: core, hnl, score, shnl
93 :
94 : ! define basis
95 24 : IF (PRESENT(orb_basis)) THEN
96 0 : CALL set_kind_basis_atomic(basis, orb_basis, has_pp=.TRUE., cp2k_norm=.TRUE.)
97 : ELSE
98 24 : CALL init_atom_basis_default_pp(basis)
99 : END IF
100 :
101 24 : is_ecp = .FALSE.
102 24 : IF (PRESENT(ecp_pot)) is_ecp = .TRUE.
103 24 : is_upf = .FALSE.
104 24 : IF (PRESENT(upf_pot)) is_upf = .TRUE.
105 24 : CPASSERT(.NOT. (is_ecp .AND. is_upf))
106 :
107 : ! upf has often very small grids, use a smooth cutoff function
108 24 : IF (is_upf) THEN
109 24 : n = SIZE(upf_pot%r)
110 72 : ALLOCATE (cutpotu(n))
111 24274 : rcut = MAXVAL(upf_pot%r)
112 24 : CALL cutpot(cutpotu, upf_pot%r, rcut, 2.5_dp)
113 24 : n = basis%grid%nr
114 72 : ALLOCATE (cutpots(n))
115 24 : CALL cutpot(cutpots, basis%grid%rad, rcut, 2.5_dp)
116 : ELSE
117 0 : n = basis%grid%nr
118 0 : ALLOCATE (cutpots(n))
119 0 : cutpots = 1.0_dp
120 : END IF
121 :
122 : ! generate the transformed potentials
123 24 : IF (is_ecp) THEN
124 0 : CALL ecp_sgp_constr(ecp_pot, sgp_pot, basis)
125 24 : ELSE IF (is_upf) THEN
126 24 : CALL upf_sgp_constr(upf_pot, sgp_pot, basis)
127 : ELSE
128 0 : CPABORT("Either ecp_pot or upf_pot is needed for sgp_construction")
129 : END IF
130 :
131 24 : NULLIFY (core, hnl)
132 24 : CALL create_opmat(core, basis%nbas)
133 24 : CALL create_opmat(hnl, basis%nbas, 5)
134 24 : NULLIFY (score, shnl)
135 24 : CALL create_opmat(score, basis%nbas)
136 24 : CALL create_opmat(shnl, basis%nbas, 5)
137 : !
138 24 : IF (is_ecp) THEN
139 0 : CALL ecpints(hnl%op, basis, ecp_pot)
140 24 : ELSE IF (is_upf) THEN
141 24 : CALL upfints(core%op, hnl%op, basis, upf_pot, cutpotu, sgp_pot%ac_local)
142 : ELSE
143 0 : CPABORT("Either ecp_pot or upf_pot is needed for sgp_construction")
144 : END IF
145 : !
146 24 : CALL sgpints(score%op, shnl%op, basis, sgp_pot, cutpots)
147 : !
148 24 : error = 0.0_dp
149 24 : IF (sgp_pot%has_local) THEN
150 24 : n = MIN(3, UBOUND(core%op, 3))
151 40440 : error(1) = MAXVAL(ABS(core%op(:, :, 0:n) - score%op(:, :, 0:n)))
152 40512 : mloc = MAXLOC(ABS(core%op(:, :, 0:n) - score%op(:, :, 0:n)))
153 : ! MAXLOC returns positions starting at one, whereas the angular index starts at zero.
154 24 : mloc(3) = mloc(3) - 1
155 24 : error(4) = error(1)/MAX(ABS(core%op(mloc(1), mloc(2), mloc(3))), 1.E-6_dp)
156 : END IF
157 24 : IF (sgp_pot%has_nonlocal) THEN
158 12 : n = MIN(3, UBOUND(hnl%op, 3))
159 20220 : error(2) = MAXVAL(ABS(hnl%op(:, :, 0:n) - shnl%op(:, :, 0:n)))
160 20256 : mloc = MAXLOC(ABS(hnl%op(:, :, 0:n) - shnl%op(:, :, 0:n)))
161 12 : mloc(3) = mloc(3) - 1
162 12 : error(5) = error(2)/MAX(ABS(hnl%op(mloc(1), mloc(2), mloc(3))), 1.E-6_dp)
163 : END IF
164 24 : IF (sgp_pot%has_nlcc) THEN
165 6 : IF (is_upf) THEN
166 6 : n = SIZE(upf_pot%r)
167 18 : ALLOCATE (cgauss(n))
168 6 : cgauss = 0.0_dp
169 294 : DO i = 1, sgp_pot%n_nlcc
170 315942 : cgauss(:) = cgauss(:) + sgp_pot%c_nlcc(i)*EXP(-sgp_pot%a_nlcc(i)*upf_pot%r(:)**2)
171 : END DO
172 6582 : errcc = SUM((cgauss(:) - upf_pot%rho_nlcc(:))**2*upf_pot%r(:)**2*upf_pot%rab(:))
173 6 : errcc = SQRT(errcc/REAL(n, KIND=dp))
174 6 : DEALLOCATE (cgauss)
175 : ELSE
176 0 : CPABORT("sgp_pot%has_nlcc requires is_upf")
177 : END IF
178 6 : error(3) = errcc
179 : END IF
180 : !
181 24 : IF (is_upf) THEN
182 24 : DEALLOCATE (cutpotu)
183 24 : DEALLOCATE (cutpots)
184 : ELSE
185 0 : DEALLOCATE (cutpots)
186 : END IF
187 : !
188 24 : CALL release_opmat(score)
189 24 : CALL release_opmat(shnl)
190 24 : CALL release_opmat(core)
191 24 : CALL release_opmat(hnl)
192 :
193 24 : CALL release_atom_basis(basis)
194 :
195 480 : END SUBROUTINE sgp_construction
196 :
197 : ! **************************************************************************************************
198 : !> \brief ...
199 : !> \param atom_info ...
200 : !> \param input_section ...
201 : !> \param iw ...
202 : ! **************************************************************************************************
203 5 : SUBROUTINE atom_sgp_construction(atom_info, input_section, iw)
204 :
205 : TYPE(atom_p_type), DIMENSION(:, :), POINTER :: atom_info
206 : TYPE(section_vals_type), POINTER :: input_section
207 : INTEGER, INTENT(IN) :: iw
208 :
209 : INTEGER :: i, n, ppot_type
210 : LOGICAL :: do_transform, explicit, is_ecp, is_upf
211 : REAL(KIND=dp) :: errcc, rcut
212 5 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cgauss, cutpots, cutpotu
213 : TYPE(atom_ecppot_type), POINTER :: ecp_pot
214 : TYPE(atom_sgp_potential_type) :: sgp_pot
215 : TYPE(atom_type), POINTER :: atom_ref
216 : TYPE(atom_upfpot_type), POINTER :: upf_pot
217 : TYPE(opmat_type), POINTER :: core, hnl, score, shnl
218 :
219 5 : CALL section_vals_get(input_section, explicit=explicit)
220 5 : IF (.NOT. explicit) RETURN
221 :
222 5 : IF (iw > 0) WRITE (iw, '(/," ",79("*"),/,T24,A,/," ",79("*"))') "SEPARABLE GAUSSIAN PSEUDOPOTENTIAL"
223 :
224 5 : atom_ref => atom_info(1, 1)%atom
225 :
226 5 : ppot_type = atom_ref%potential%ppot_type
227 0 : SELECT CASE (ppot_type)
228 : CASE (gth_pseudo)
229 0 : IF (iw > 0) WRITE (iw, '(" GTH Pseudopotential is already in SGP form. ")')
230 3 : do_transform = .FALSE.
231 : CASE (ecp_pseudo)
232 3 : do_transform = .TRUE.
233 3 : is_ecp = .TRUE.
234 3 : is_upf = .FALSE.
235 3 : ecp_pot => atom_ref%potential%ecp_pot
236 : CASE (upf_pseudo)
237 2 : do_transform = .TRUE.
238 2 : is_ecp = .FALSE.
239 2 : is_upf = .TRUE.
240 2 : upf_pot => atom_ref%potential%upf_pot
241 : CASE (no_pseudo)
242 0 : IF (iw > 0) WRITE (iw, '(" No Pseudopotential available for transformation. ")')
243 0 : do_transform = .FALSE.
244 : CASE DEFAULT
245 5 : CPABORT("Unknown ppot_type")
246 : END SELECT
247 :
248 : ! generate the transformed potentials
249 0 : IF (do_transform) THEN
250 5 : IF (is_ecp) THEN
251 3 : CALL ecp_sgp_constr(ecp_pot, sgp_pot, atom_ref%basis)
252 2 : ELSE IF (is_upf) THEN
253 2 : CALL upf_sgp_constr(upf_pot, sgp_pot, atom_ref%basis)
254 : ELSE
255 0 : CPABORT("Either ecp_pseudo or upf_pseudo is needed for atom_sgp_construction")
256 : END IF
257 : END IF
258 :
259 : ! Check the result
260 5 : IF (do_transform) THEN
261 5 : NULLIFY (core, hnl)
262 5 : CALL create_opmat(core, atom_ref%basis%nbas)
263 5 : CALL create_opmat(hnl, atom_ref%basis%nbas, 5)
264 5 : NULLIFY (score, shnl)
265 5 : CALL create_opmat(score, atom_ref%basis%nbas)
266 5 : CALL create_opmat(shnl, atom_ref%basis%nbas, 5)
267 : !
268 : ! upf has often very small grids, use a smooth cutoff function
269 5 : IF (is_upf) THEN
270 2 : n = SIZE(upf_pot%r)
271 6 : ALLOCATE (cutpotu(n))
272 1536 : rcut = MAXVAL(upf_pot%r)
273 2 : CALL cutpot(cutpotu, upf_pot%r, rcut, 2.5_dp)
274 2 : n = atom_ref%basis%grid%nr
275 6 : ALLOCATE (cutpots(n))
276 2 : CALL cutpot(cutpots, atom_ref%basis%grid%rad, rcut, 2.5_dp)
277 : ELSE
278 3 : n = atom_ref%basis%grid%nr
279 9 : ALLOCATE (cutpots(n))
280 1203 : cutpots = 1.0_dp
281 : END IF
282 : !
283 5 : IF (is_ecp) THEN
284 3 : CALL ecpints(hnl%op, atom_ref%basis, ecp_pot)
285 2 : ELSE IF (is_upf) THEN
286 2 : CALL upfints(core%op, hnl%op, atom_ref%basis, upf_pot, cutpotu, sgp_pot%ac_local)
287 : ELSE
288 0 : CPABORT("Either ecp_pseudo or upf_pseudo is needed for atom_sgp_construction")
289 : END IF
290 : !
291 5 : CALL sgpints(score%op, shnl%op, atom_ref%basis, sgp_pot, cutpots)
292 : !
293 5 : IF (sgp_pot%has_local) THEN
294 2 : n = MIN(3, UBOUND(core%op, 3))
295 2186 : errcc = MAXVAL(ABS(core%op(:, :, 0:n) - score%op(:, :, 0:n)))
296 2 : IF (iw > 0) THEN
297 2 : WRITE (iw, '(" Local part of pseudopotential")')
298 2 : WRITE (iw, '(" Number of basis functions ",T77,i4)') sgp_pot%n_local
299 2 : WRITE (iw, '(" Max. abs. error of matrix elements ",T65,f16.8)') errcc
300 : END IF
301 : END IF
302 5 : IF (sgp_pot%has_nonlocal) THEN
303 5 : IF (iw > 0) THEN
304 5 : WRITE (iw, '(" Nonlocal part of pseudopotential")')
305 5 : WRITE (iw, '(" Max. l-quantum number",T77,i4)') sgp_pot%lmax
306 5 : WRITE (iw, '(" Number of projector basis functions ",T77,i4)') sgp_pot%n_nonlocal
307 21 : DO i = 0, sgp_pot%lmax
308 4368 : errcc = MAXVAL(ABS(hnl%op(:, :, i) - shnl%op(:, :, i)))
309 21 : WRITE (iw, '(" Max. abs. error of matrix elements: l=",i2,T69,f12.8)') i, errcc
310 : END DO
311 : END IF
312 : END IF
313 5 : IF (sgp_pot%has_nlcc) THEN
314 0 : IF (is_upf) THEN
315 0 : n = SIZE(upf_pot%r)
316 0 : ALLOCATE (cgauss(n))
317 0 : cgauss = 0.0_dp
318 0 : DO i = 1, sgp_pot%n_nlcc
319 0 : cgauss(:) = cgauss(:) + sgp_pot%c_nlcc(i)*EXP(-sgp_pot%a_nlcc(i)*upf_pot%r(:)**2)
320 : END DO
321 0 : errcc = SUM((cgauss(:) - upf_pot%rho_nlcc(:))**2*upf_pot%r(:)**2*upf_pot%rab(:))
322 0 : errcc = SQRT(errcc/REAL(n, KIND=dp))
323 0 : DEALLOCATE (cgauss)
324 : ELSE
325 0 : CPABORT("sgp_pot%has_nlcc requires is_upf")
326 : END IF
327 0 : IF (iw > 0) THEN
328 0 : WRITE (iw, '(" Non-linear core correction: core density")')
329 0 : WRITE (iw, '(" Number of basis functions ",T77,i4)') sgp_pot%n_nlcc
330 0 : WRITE (iw, '(" RMS error of core density ",T69,f12.8)') errcc
331 : END IF
332 : END IF
333 : !
334 5 : IF (is_upf) THEN
335 2 : DEALLOCATE (cutpotu)
336 2 : DEALLOCATE (cutpots)
337 : ELSE
338 3 : DEALLOCATE (cutpots)
339 : END IF
340 : !
341 5 : CALL release_opmat(score)
342 5 : CALL release_opmat(shnl)
343 5 : CALL release_opmat(core)
344 5 : CALL release_opmat(hnl)
345 : END IF
346 :
347 5 : CALL atom_sgp_release(sgp_pot)
348 :
349 5 : IF (iw > 0) WRITE (iw, '(" ",79("*"))')
350 :
351 40 : END SUBROUTINE atom_sgp_construction
352 :
353 : ! **************************************************************************************************
354 : !> \brief ...
355 : !> \param ecp_pot ...
356 : !> \param sgp_pot ...
357 : !> \param basis ...
358 : ! **************************************************************************************************
359 3 : SUBROUTINE ecp_sgp_constr(ecp_pot, sgp_pot, basis)
360 :
361 : TYPE(atom_ecppot_type) :: ecp_pot
362 : TYPE(atom_sgp_potential_type) :: sgp_pot
363 : TYPE(atom_basis_type) :: basis
364 :
365 : INTEGER :: i, ia, ir, j, k, l, n, na, nl, nr
366 : REAL(KIND=dp) :: alpha, amx, cmx
367 3 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: al, cl, cpot, pgauss
368 3 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: cmat, qmat, score, sinv, smat, tmat
369 3 : REAL(KIND=dp), DIMENSION(:), POINTER :: rad
370 :
371 3 : sgp_pot%has_nlcc = .FALSE.
372 3 : sgp_pot%n_nlcc = 0
373 3 : sgp_pot%has_local = .FALSE.
374 3 : sgp_pot%n_local = 0
375 :
376 : ! transform semilocal potential into a separable local form
377 3 : sgp_pot%has_nonlocal = .TRUE.
378 : !
379 : ! hardcoded number of projectors (equal for all l values)
380 3 : nl = 8
381 : !
382 3 : sgp_pot%n_nonlocal = nl
383 3 : sgp_pot%lmax = ecp_pot%lmax
384 0 : ALLOCATE (sgp_pot%a_nonlocal(nl))
385 9 : ALLOCATE (sgp_pot%h_nonlocal(nl, 0:ecp_pot%lmax))
386 9 : ALLOCATE (sgp_pot%c_nonlocal(nl, nl, 0:ecp_pot%lmax))
387 : !
388 3 : ALLOCATE (al(nl), cl(nl))
389 3 : ALLOCATE (smat(nl, nl), sinv(nl, nl))
390 3 : ALLOCATE (tmat(nl, nl), cmat(nl, nl))
391 3 : al = 0.0_dp
392 531 : amx = MAXVAL(ecp_pot%bpot)
393 3 : cmx = amx/0.15_dp
394 3 : cmx = cmx**(1._dp/REAL(nl - 1, dp))
395 3 : cmx = 1._dp/cmx
396 27 : DO ir = 1, nl
397 27 : al(ir) = amx*cmx**(ir - 1)
398 : END DO
399 : !
400 27 : sgp_pot%a_nonlocal(1:nl) = al(1:nl)
401 : !
402 3 : nr = basis%grid%nr
403 3 : rad => basis%grid%rad
404 12 : ALLOCATE (cpot(nr), pgauss(nr))
405 14 : DO l = 0, ecp_pot%lmax
406 11 : na = basis%nbas(l)
407 66 : ALLOCATE (score(na, na), qmat(na, nl))
408 11 : cpot = 0._dp
409 26 : DO k = 1, ecp_pot%npot(l)
410 15 : n = ecp_pot%nrpot(k, l)
411 15 : alpha = ecp_pot%bpot(k, l)
412 6026 : cpot(:) = cpot + ecp_pot%apot(k, l)*rad**(n - 2)*EXP(-alpha*rad**2)
413 : END DO
414 187 : DO i = 1, na
415 1683 : DO j = i, na
416 1496 : score(i, j) = integrate_grid(cpot, basis%bf(:, i, l), basis%bf(:, j, l), basis%grid)
417 1672 : score(j, i) = score(i, j)
418 : END DO
419 : END DO
420 : ! overlap basis with projectors
421 99 : DO i = 1, nl
422 35288 : pgauss(:) = EXP(-al(i)*rad(:)**2)*rad(:)**l
423 1507 : DO ia = 1, na
424 1496 : qmat(ia, i) = integrate_grid(basis%bf(:, ia, l), pgauss(:), basis%grid)
425 : END DO
426 : END DO
427 1507 : qmat = SQRT(fourpi)*qmat
428 : ! tmat = qmat * score * qmat
429 70708 : tmat(1:nl, 1:nl) = MATMUL(TRANSPOSE(qmat(1:na, 1:nl)), MATMUL(score(1:na, 1:na), qmat(1:na, 1:nl)))
430 12859 : smat(1:nl, 1:nl) = MATMUL(TRANSPOSE(qmat(1:na, 1:nl)), qmat(1:na, 1:nl))
431 11 : CALL get_pseudo_inverse_diag(smat(1:nl, 1:nl), sinv(1:nl, 1:nl), 1.e-10_dp)
432 24915 : cmat(1:nl, 1:nl) = MATMUL(sinv(1:nl, 1:nl), MATMUL(tmat(1:nl, 1:nl), sinv(1:nl, 1:nl)))
433 1595 : cmat(1:nl, 1:nl) = (cmat(1:nl, 1:nl) + TRANSPOSE(cmat(1:nl, 1:nl)))*0.5_dp
434 : !
435 : ! Residium
436 : !
437 11 : CALL diamat_all(cmat(1:nl, 1:nl), cl(1:nl), .TRUE.)
438 : !
439 99 : sgp_pot%h_nonlocal(1:nl, l) = cl(1:nl)
440 803 : sgp_pot%c_nonlocal(1:nl, 1:nl, l) = cmat(1:nl, 1:nl)
441 11 : sgp_pot%is_nonlocal(l) = .TRUE.
442 : !
443 14 : DEALLOCATE (score, qmat)
444 : END DO
445 3 : DEALLOCATE (cpot, pgauss)
446 3 : DEALLOCATE (al, cl, smat, sinv, tmat, cmat)
447 :
448 6 : END SUBROUTINE ecp_sgp_constr
449 :
450 : ! **************************************************************************************************
451 : !> \brief ...
452 : !> \param upf_pot ...
453 : !> \param sgp_pot ...
454 : !> \param basis ...
455 : ! **************************************************************************************************
456 26 : SUBROUTINE upf_sgp_constr(upf_pot, sgp_pot, basis)
457 :
458 : TYPE(atom_upfpot_type) :: upf_pot
459 : TYPE(atom_sgp_potential_type) :: sgp_pot
460 : TYPE(atom_basis_type) :: basis
461 :
462 : CHARACTER(len=4) :: ptype
463 : INTEGER :: ia, ib, ipa, ipb, ir, la, lb, na, nl, &
464 : np, nr
465 : LOGICAL :: nl_trans
466 : REAL(KIND=dp) :: cpa, cpb, eee, ei, errloc, rc, x(2), zval
467 52 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: al, cl, pgauss, pproa, pprob, vgauss, &
468 26 : vloc, ww
469 26 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: cmat, qmat, score, sinv, smat, tmat
470 : TYPE(atom_basis_type) :: gbasis
471 : TYPE(opt_state_type) :: ostate
472 :
473 26 : IF (upf_pot%is_ultrasoft .OR. upf_pot%is_paw .OR. upf_pot%is_coulomb) THEN
474 0 : sgp_pot%has_nonlocal = .FALSE.
475 0 : sgp_pot%n_nonlocal = 0
476 0 : sgp_pot%has_local = .FALSE.
477 0 : sgp_pot%n_local = 0
478 0 : sgp_pot%has_nlcc = .FALSE.
479 0 : sgp_pot%n_nlcc = 0
480 0 : RETURN
481 : END IF
482 :
483 : ! radial grid
484 26 : nr = SIZE(upf_pot%r)
485 : ! weights for integration
486 78 : ALLOCATE (ww(nr))
487 25810 : ww(:) = upf_pot%r(:)**2*upf_pot%rab(:)
488 :
489 : ! start with local potential
490 26 : sgp_pot%has_local = .TRUE.
491 : ! fit local potential to Gaussian form
492 78 : ALLOCATE (vloc(nr), vgauss(nr))
493 : ! smearing of core charge
494 26 : zval = upf_pot%zion
495 : ! Try to find an optimal Gaussian charge distribution
496 26 : CALL erffit(sgp_pot%ac_local, upf_pot%vlocal, upf_pot%r, zval)
497 26 : sgp_pot%zval = zval
498 25810 : DO ir = 1, nr
499 25810 : IF (upf_pot%r(ir) < 1.e-12_dp) THEN
500 6 : vgauss(ir) = -SQRT(2.0_dp)*zval/rootpi/sgp_pot%ac_local
501 : ELSE
502 25778 : rc = upf_pot%r(ir)/sgp_pot%ac_local/SQRT(2.0_dp)
503 25778 : vgauss(ir) = -zval/upf_pot%r(ir)*erf(rc)
504 : END IF
505 : END DO
506 25810 : vloc(:) = upf_pot%vlocal(:) - vgauss(:)
507 : !
508 26 : CALL atom_basis_gridrep(basis, gbasis, upf_pot%r, upf_pot%rab)
509 : !
510 26 : nl = 12
511 26 : ALLOCATE (al(nl), cl(nl))
512 26 : ostate%nf = 0
513 26 : ostate%nvar = 2
514 26 : x(1) = 1.00_dp !starting point of geometric series
515 26 : x(2) = 1.20_dp !factor of series
516 26 : ostate%rhoend = 1.e-12_dp
517 26 : ostate%rhobeg = 5.e-2_dp
518 26 : ostate%maxfun = 1000
519 26 : ostate%iprint = 1
520 26 : ostate%unit = -1
521 26 : ostate%state = 0
522 2273 : DO
523 2299 : IF (ostate%state == 2) THEN
524 28873 : DO ir = 1, nl
525 28873 : al(ir) = x(1)*x(2)**(ir - 1)
526 : END DO
527 2221 : CALL pplocal_error(nl, al, cl, vloc, vgauss, gbasis, upf_pot%r, ww, 1, ostate%f)
528 : END IF
529 2299 : IF (ostate%state == -1) EXIT
530 2273 : CALL powell_optimize(ostate%nvar, x, ostate)
531 : END DO
532 26 : ostate%state = 8
533 26 : CALL powell_optimize(ostate%nvar, x, ostate)
534 338 : DO ir = 1, nl
535 338 : al(ir) = x(1)*x(2)**(ir - 1)
536 : END DO
537 26 : CALL pplocal_error(nl, al, cl, vloc, vgauss, gbasis, upf_pot%r, ww, 1, errloc)
538 : !
539 26 : ALLOCATE (sgp_pot%a_local(nl), sgp_pot%c_local(nl))
540 26 : sgp_pot%n_local = nl
541 338 : sgp_pot%a_local(1:nl) = al(1:nl)
542 338 : sgp_pot%c_local(1:nl) = cl(1:nl)
543 26 : DEALLOCATE (vloc, vgauss)
544 26 : DEALLOCATE (al, cl)
545 26 : CALL release_atom_basis(gbasis)
546 : !
547 26 : ptype = ADJUSTL(TRIM(upf_pot%pseudo_type))
548 26 : IF (ptype(1:2) == "NC" .OR. ptype(1:2) == "US") THEN
549 : nl_trans = .FALSE.
550 1 : ELSE IF (ptype(1:2) == "SL") THEN
551 : nl_trans = .TRUE.
552 : ELSE
553 0 : CPABORT("Pseudopotential type: ["//ADJUSTL(TRIM(ptype))//"] not known")
554 : END IF
555 :
556 : ! purely local pseudopotentials
557 26 : IF (upf_pot%l_max < 0) THEN
558 12 : sgp_pot%n_nonlocal = 0
559 12 : sgp_pot%lmax = -1
560 12 : sgp_pot%has_nonlocal = .FALSE.
561 : ELSE
562 : ! Non-local pseudopotential in Gaussian form
563 14 : IF (nl_trans) THEN
564 1 : sgp_pot%has_nonlocal = .TRUE.
565 : ! semi local pseudopotential
566 : ! fit to nonlocal form
567 : ! get basis representation on UPF grid
568 1 : nl = 8
569 : !
570 1 : sgp_pot%n_nonlocal = nl
571 1 : sgp_pot%lmax = upf_pot%l_max
572 1 : ALLOCATE (sgp_pot%a_nonlocal(nl))
573 3 : ALLOCATE (sgp_pot%h_nonlocal(nl, 0:upf_pot%l_max))
574 3 : ALLOCATE (sgp_pot%c_nonlocal(nl, nl, 0:upf_pot%l_max))
575 : !
576 1 : ALLOCATE (al(nl), cl(nl))
577 1 : ALLOCATE (smat(nl, nl), sinv(nl, nl))
578 1 : ALLOCATE (tmat(nl, nl), cmat(nl, nl))
579 1 : al = 0.0_dp
580 9 : DO ir = 1, nl
581 9 : al(ir) = 10.0_dp*0.60_dp**(ir - 1)
582 : END DO
583 : !
584 9 : sgp_pot%a_nonlocal(1:nl) = al(1:nl)
585 : !
586 1 : CALL atom_basis_gridrep(basis, gbasis, upf_pot%r, upf_pot%rab)
587 3 : ALLOCATE (pgauss(nr), vloc(nr))
588 5 : DO la = 0, upf_pot%l_max
589 4 : IF (la == upf_pot%l_local) CYCLE
590 3 : sgp_pot%is_nonlocal(la) = .TRUE.
591 3 : na = gbasis%nbas(la)
592 21 : ALLOCATE (score(na, na), qmat(na, nl))
593 : ! Reference matrix
594 1386 : vloc(:) = upf_pot%vsemi(:, la + 1) - upf_pot%vlocal(:)
595 51 : DO ia = 1, na
596 459 : DO ib = ia, na
597 188496 : score(ia, ib) = SUM(vloc(:)*gbasis%bf(:, ia, la)*gbasis%bf(:, ib, la)*ww(:))
598 456 : score(ib, ia) = score(ia, ib)
599 : END DO
600 : END DO
601 : ! overlap basis with projectors
602 27 : DO ir = 1, nl
603 11088 : pgauss(:) = EXP(-al(ir)*upf_pot%r(:)**2)*upf_pot%r(:)**la
604 24 : eee = rootpi/(2._dp**(la + 2)*dfac(2*la + 1))/(2._dp*al(ir))**(la + 1.5_dp)
605 11088 : pgauss(:) = pgauss(:)/SQRT(eee)
606 411 : DO ia = 1, na
607 177432 : qmat(ia, ir) = SUM(gbasis%bf(:, ia, la)*pgauss(:)*ww)
608 : END DO
609 : END DO
610 : ! tmat = qmat * score * qmat
611 19284 : tmat(1:nl, 1:nl) = MATMUL(TRANSPOSE(qmat(1:na, 1:nl)), MATMUL(score(1:na, 1:na), qmat(1:na, 1:nl)))
612 3507 : smat(1:nl, 1:nl) = MATMUL(TRANSPOSE(qmat(1:na, 1:nl)), qmat(1:na, 1:nl))
613 3 : CALL get_pseudo_inverse_diag(smat(1:nl, 1:nl), sinv(1:nl, 1:nl), 1.e-10_dp)
614 6795 : cmat(1:nl, 1:nl) = MATMUL(sinv(1:nl, 1:nl), MATMUL(tmat(1:nl, 1:nl), sinv(1:nl, 1:nl)))
615 435 : cmat(1:nl, 1:nl) = (cmat(1:nl, 1:nl) + TRANSPOSE(cmat(1:nl, 1:nl)))*0.5_dp
616 3 : CALL diamat_all(cmat(1:nl, 1:nl), cl(1:nl), .TRUE.)
617 : !
618 : ! get back unnormalized Gaussians
619 27 : DO ir = 1, nl
620 24 : ei = rootpi/(2._dp**(la + 2)*dfac(2*la + 1))/(2._dp*al(ir))**(la + 1.5_dp)
621 219 : cmat(ir, 1:nl) = cmat(ir, 1:nl)/SQRT(ei)
622 : END DO
623 27 : sgp_pot%h_nonlocal(1:nl, la) = cl(1:nl)
624 219 : sgp_pot%c_nonlocal(1:nl, 1:nl, la) = cmat(1:nl, 1:nl)
625 3 : sgp_pot%is_nonlocal(la) = .TRUE.
626 4 : DEALLOCATE (score, qmat)
627 : END DO
628 : ! SQRT(4PI)
629 293 : sgp_pot%c_nonlocal = sgp_pot%c_nonlocal/SQRT(fourpi)
630 1 : CALL release_atom_basis(gbasis)
631 1 : DEALLOCATE (pgauss, vloc)
632 1 : DEALLOCATE (al, cl, smat, sinv, tmat, cmat)
633 : ELSE
634 13 : sgp_pot%has_nonlocal = .TRUE.
635 : ! non local pseudopotential
636 52 : ALLOCATE (pproa(nr), pprob(nr), pgauss(nr))
637 13 : np = upf_pot%number_of_proj
638 13 : nl = 8
639 13 : ALLOCATE (al(nl), cl(nl))
640 13 : ALLOCATE (smat(nl, nl), sinv(nl, nl))
641 13 : ALLOCATE (tmat(nl, nl), cmat(nl, nl))
642 13 : al = 0.0_dp
643 13 : cl = 0.0_dp
644 117 : DO ir = 1, nl
645 117 : al(ir) = 10.0_dp*0.60_dp**(ir - 1)
646 : END DO
647 : !
648 56 : sgp_pot%lmax = MAXVAL(upf_pot%lbeta(:))
649 13 : sgp_pot%n_nonlocal = nl
650 13 : ALLOCATE (sgp_pot%a_nonlocal(nl))
651 39 : ALLOCATE (sgp_pot%h_nonlocal(nl, 0:sgp_pot%lmax))
652 39 : ALLOCATE (sgp_pot%c_nonlocal(nl, nl, 0:sgp_pot%lmax))
653 : !
654 117 : sgp_pot%a_nonlocal(1:nl) = al(1:nl)
655 : !
656 13 : CALL atom_basis_gridrep(basis, gbasis, upf_pot%r, upf_pot%rab)
657 38 : DO la = 0, sgp_pot%lmax
658 25 : sgp_pot%is_nonlocal(la) = .TRUE.
659 25 : na = gbasis%nbas(la)
660 175 : ALLOCATE (score(na, na), qmat(na, nl))
661 : ! Reference matrix
662 25 : score = 0.0_dp
663 140 : DO ipa = 1, np
664 115 : lb = upf_pot%lbeta(ipa)
665 115 : IF (la /= lb) CYCLE
666 47098 : pproa(:) = upf_pot%beta(:, ipa)
667 291 : DO ipb = 1, np
668 223 : lb = upf_pot%lbeta(ipb)
669 223 : IF (la /= lb) CYCLE
670 86590 : pprob(:) = upf_pot%beta(:, ipb)
671 79 : eee = upf_pot%dion(ipa, ipb)
672 1770 : DO ia = 1, na
673 1727504 : cpa = SUM(pproa(:)*gbasis%bf(:, ia, la)*ww(:))
674 18315 : DO ib = ia, na
675 18104424 : cpb = SUM(pprob(:)*gbasis%bf(:, ib, la)*ww(:))
676 16516 : score(ia, ib) = score(ia, ib) + cpa*eee*cpb
677 18092 : score(ib, ia) = score(ia, ib)
678 : END DO
679 : END DO
680 : END DO
681 : END DO
682 : ! overlap basis with projectors
683 225 : DO ir = 1, nl
684 218816 : pgauss(:) = EXP(-al(ir)*upf_pot%r(:)**2)*upf_pot%r(:)**la
685 200 : eee = rootpi/(2._dp**(la + 2)*dfac(2*la + 1))/(2._dp*al(ir))**(la + 1.5_dp)
686 218816 : pgauss(:) = pgauss(:)/SQRT(eee)
687 4193 : DO ia = 1, na
688 4342152 : qmat(ia, ir) = SUM(gbasis%bf(:, ia, la)*pgauss(:)*ww)
689 : END DO
690 : END DO
691 : ! tmat = qmat * score * qmat
692 229052 : tmat(1:nl, 1:nl) = MATMUL(TRANSPOSE(qmat(1:na, 1:nl)), MATMUL(score(1:na, 1:na), qmat(1:na, 1:nl)))
693 35369 : smat(1:nl, 1:nl) = MATMUL(TRANSPOSE(qmat(1:na, 1:nl)), qmat(1:na, 1:nl))
694 25 : CALL get_pseudo_inverse_diag(smat(1:nl, 1:nl), sinv(1:nl, 1:nl), 1.e-10_dp)
695 56625 : cmat(1:nl, 1:nl) = MATMUL(sinv(1:nl, 1:nl), MATMUL(tmat(1:nl, 1:nl), sinv(1:nl, 1:nl)))
696 3625 : cmat(1:nl, 1:nl) = (cmat(1:nl, 1:nl) + TRANSPOSE(cmat(1:nl, 1:nl)))*0.5_dp
697 25 : CALL diamat_all(cmat(1:nl, 1:nl), cl(1:nl), .TRUE.)
698 : !
699 : ! get back unnormalized Gaussians
700 225 : DO ir = 1, nl
701 200 : ei = rootpi/(2._dp**(la + 2)*dfac(2*la + 1))/(2._dp*al(ir))**(la + 1.5_dp)
702 1825 : cmat(ir, 1:nl) = cmat(ir, 1:nl)/SQRT(ei)
703 : END DO
704 225 : sgp_pot%h_nonlocal(1:nl, la) = cl(1:nl)
705 1825 : sgp_pot%c_nonlocal(1:nl, 1:nl, la) = cmat(1:nl, 1:nl)
706 25 : sgp_pot%is_nonlocal(la) = .TRUE.
707 38 : DEALLOCATE (score, qmat)
708 : END DO
709 : ! SQRT(4PI)
710 1838 : sgp_pot%c_nonlocal = sgp_pot%c_nonlocal/SQRT(fourpi)
711 13 : CALL release_atom_basis(gbasis)
712 13 : DEALLOCATE (pgauss, pproa, pprob)
713 13 : DEALLOCATE (al, cl, smat, sinv, tmat, cmat)
714 : END IF
715 : END IF
716 :
717 26 : IF (upf_pot%core_correction) THEN
718 6 : sgp_pot%has_nlcc = .TRUE.
719 : ELSE
720 20 : sgp_pot%has_nlcc = .FALSE.
721 20 : sgp_pot%n_nlcc = 0
722 : END IF
723 :
724 : ! fit core charge to Gaussian form
725 26 : IF (sgp_pot%has_nlcc) THEN
726 6 : CALL fit_core_density_gaussians(upf_pot%r, upf_pot%rab, upf_pot%rho_nlcc, al, cl)
727 6 : sgp_pot%n_nlcc = SIZE(al)
728 306 : ALLOCATE (sgp_pot%a_nlcc, SOURCE=al)
729 306 : ALLOCATE (sgp_pot%c_nlcc, SOURCE=cl)
730 6 : DEALLOCATE (al, cl)
731 : END IF
732 :
733 26 : DEALLOCATE (ww)
734 :
735 520 : END SUBROUTINE upf_sgp_constr
736 :
737 : ! **************************************************************************************************
738 : !> \brief Fit a tabulated core density without forming ill-conditioned normal equations.
739 : !> Column-normalized, radial least squares use a direct SVD. A broad Gaussian
740 : !> sequence resolves both the core and its smooth tail. Zero-density samples
741 : !> beyond the tabulated mesh prevent an unconstrained diffuse continuation.
742 : !> \param r Radial grid in bohr.
743 : !> \param rab Radial mesh Jacobian dr/dj.
744 : !> \param rho True spherical density (without the 4*pi*r**2 prefactor).
745 : !> \param alpha Gaussian exponents in inverse square bohr.
746 : !> \param coeff Coefficients multiplying exp(-alpha*r**2).
747 : ! **************************************************************************************************
748 14 : SUBROUTINE fit_core_density_gaussians(r, rab, rho, alpha, coeff)
749 : REAL(dp), INTENT(IN) :: r(:), rab(:), rho(:)
750 : REAL(dp), ALLOCATABLE, INTENT(OUT) :: alpha(:), coeff(:)
751 :
752 : INTEGER, PARAMETER :: ng = 48, ntail = 200
753 : REAL(dp), PARAMETER :: amax = 300.0_dp, amin = 0.03_dp, &
754 : rcond = 1.e-12_dp
755 :
756 : INTEGER :: i, info, lwork, m, nr, rank
757 : REAL(dp) :: dr, extent, norms(ng), singular(ng), &
758 : work_query(1)
759 14 : REAL(dp), ALLOCATABLE :: matrix(:, :), rhs(:, :), rr(:), sw(:), &
760 14 : work(:)
761 :
762 14 : nr = SIZE(r)
763 14 : CPASSERT(nr >= 2 .AND. SIZE(rab) == nr .AND. SIZE(rho) == nr)
764 38796 : CPASSERT(ALL(r >= 0.0_dp) .AND. ALL(rab > 0.0_dp))
765 19384 : CPASSERT(ALL(r(2:) > r(:nr - 1)))
766 14 : m = nr + ntail
767 98 : ALLOCATE (alpha(ng), coeff(ng), matrix(m, ng), rhs(m, 1), rr(m), sw(m))
768 19398 : rr(:nr) = r
769 19398 : sw(:nr) = r*SQRT(rab)
770 : ! At 50 bohr even the most diffuse fitting Gaussian is below exp(-75).
771 14 : extent = MAX(r(nr) + 1.0_dp, 50.0_dp)
772 14 : dr = (extent - r(nr))/REAL(ntail, dp)
773 2814 : DO i = 1, ntail
774 2800 : rr(nr + i) = r(nr) + REAL(i, dp)*dr
775 2814 : sw(nr + i) = rr(nr + i)*SQRT(dr)
776 : END DO
777 14 : sw(m) = sw(m)/SQRT(2.0_dp)
778 14 : rhs = 0.0_dp
779 19398 : rhs(:nr, 1) = rho*sw(:nr)
780 686 : DO i = 1, ng
781 672 : alpha(i) = amin*EXP(LOG(amax/amin)*REAL(i - 1, dp)/REAL(ng - 1, dp))
782 1065504 : matrix(:, i) = EXP(-alpha(i)*rr**2)*sw
783 1065504 : norms(i) = SQRT(SUM(matrix(:, i)**2))
784 672 : CPASSERT(norms(i) > 0.0_dp)
785 1065518 : matrix(:, i) = matrix(:, i)/norms(i)
786 : END DO
787 14 : CALL dgelss(m, ng, 1, matrix, m, rhs, m, singular, rcond, rank, work_query, -1, info)
788 14 : CPASSERT(info == 0)
789 14 : lwork = MAX(1, CEILING(work_query(1)))
790 42 : ALLOCATE (work(lwork))
791 14 : CALL dgelss(m, ng, 1, matrix, m, rhs, m, singular, rcond, rank, work, lwork, info)
792 14 : IF (info /= 0) CPABORT("SVD failed while fitting the UPF core density")
793 686 : coeff(:) = rhs(:ng, 1)/norms
794 14 : END SUBROUTINE fit_core_density_gaussians
795 :
796 : ! **************************************************************************************************
797 : !> \brief ...
798 : !> \param sgp_pot ...
799 : ! **************************************************************************************************
800 115 : SUBROUTINE atom_sgp_release(sgp_pot)
801 :
802 : TYPE(atom_sgp_potential_type) :: sgp_pot
803 :
804 115 : IF (ASSOCIATED(sgp_pot%a_nonlocal)) DEALLOCATE (sgp_pot%a_nonlocal)
805 115 : IF (ASSOCIATED(sgp_pot%h_nonlocal)) DEALLOCATE (sgp_pot%h_nonlocal)
806 115 : IF (ASSOCIATED(sgp_pot%c_nonlocal)) DEALLOCATE (sgp_pot%c_nonlocal)
807 :
808 115 : IF (ASSOCIATED(sgp_pot%a_local)) DEALLOCATE (sgp_pot%a_local)
809 115 : IF (ASSOCIATED(sgp_pot%c_local)) DEALLOCATE (sgp_pot%c_local)
810 :
811 115 : IF (ASSOCIATED(sgp_pot%a_nlcc)) DEALLOCATE (sgp_pot%a_nlcc)
812 115 : IF (ASSOCIATED(sgp_pot%c_nlcc)) DEALLOCATE (sgp_pot%c_nlcc)
813 :
814 805 : sgp_pot = atom_sgp_potential_type()
815 :
816 115 : END SUBROUTINE atom_sgp_release
817 :
818 : ! **************************************************************************************************
819 : !> \brief ...
820 : !> \param core ...
821 : !> \param hnl ...
822 : !> \param basis ...
823 : !> \param upf_pot ...
824 : !> \param cutpotu ...
825 : !> \param ac_local ...
826 : ! **************************************************************************************************
827 26 : SUBROUTINE upfints(core, hnl, basis, upf_pot, cutpotu, ac_local)
828 : REAL(KIND=dp), DIMENSION(:, :, 0:) :: core, hnl
829 : TYPE(atom_basis_type), INTENT(INOUT) :: basis
830 : TYPE(atom_upfpot_type) :: upf_pot
831 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: cutpotu
832 : REAL(KIND=dp), INTENT(IN) :: ac_local
833 :
834 : CHARACTER(len=4) :: ptype
835 : INTEGER :: i, j, k1, k2, la, lb, m, n
836 : REAL(KIND=dp) :: rc, zval
837 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: spot
838 26 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: spmat
839 : TYPE(atom_basis_type) :: gbasis
840 :
841 : ! get basis representation on UPF grid
842 26 : CALL atom_basis_gridrep(basis, gbasis, upf_pot%r, upf_pot%rab)
843 :
844 : ! local pseudopotential
845 63926 : core = 0._dp
846 26 : n = SIZE(upf_pot%r)
847 78 : ALLOCATE (spot(n))
848 25810 : spot(:) = upf_pot%vlocal(:)
849 26 : zval = upf_pot%zion
850 25810 : DO i = 1, n
851 25810 : IF (upf_pot%r(i) < 1.e-12_dp) THEN
852 6 : spot(i) = spot(i) + sqrt2*zval/rootpi/ac_local
853 : ELSE
854 25778 : rc = upf_pot%r(i)/ac_local/sqrt2
855 25778 : spot(i) = spot(i) + zval/upf_pot%r(i)*erf(rc)
856 : END IF
857 : END DO
858 25810 : spot(:) = spot(:)*cutpotu(:)
859 :
860 26 : CALL numpot_matrix(core, spot, gbasis, 0)
861 26 : DEALLOCATE (spot)
862 :
863 63926 : hnl = 0._dp
864 26 : ptype = ADJUSTL(TRIM(upf_pot%pseudo_type))
865 26 : IF (ptype(1:2) == "NC" .OR. ptype(1:2) == "US") THEN
866 : ! non local pseudopotential
867 175 : n = MAXVAL(gbasis%nbas(:))
868 25 : m = upf_pot%number_of_proj
869 88 : ALLOCATE (spmat(n, m))
870 25 : spmat = 0.0_dp
871 68 : DO i = 1, m
872 43 : la = upf_pot%lbeta(i)
873 924 : DO j = 1, gbasis%nbas(la)
874 899 : spmat(j, i) = integrate_grid(upf_pot%beta(:, i), gbasis%bf(:, j, la), gbasis%grid)
875 : END DO
876 : END DO
877 68 : DO i = 1, m
878 43 : la = upf_pot%lbeta(i)
879 291 : DO j = 1, m
880 223 : lb = upf_pot%lbeta(j)
881 266 : IF (la == lb) THEN
882 1655 : DO k1 = 1, gbasis%nbas(la)
883 33111 : DO k2 = 1, gbasis%nbas(la)
884 33032 : hnl(k1, k2, la) = hnl(k1, k2, la) + spmat(k1, i)*upf_pot%dion(i, j)*spmat(k2, j)
885 : END DO
886 : END DO
887 : END IF
888 : END DO
889 : END DO
890 25 : DEALLOCATE (spmat)
891 1 : ELSE IF (ptype(1:2) == "SL") THEN
892 : ! semi local pseudopotential
893 5 : DO la = 0, upf_pot%l_max
894 4 : IF (la == upf_pot%l_local) CYCLE
895 3 : m = SIZE(upf_pot%vsemi(:, la + 1))
896 9 : ALLOCATE (spot(m))
897 1386 : spot(:) = upf_pot%vsemi(:, la + 1) - upf_pot%vlocal(:)
898 1386 : spot(:) = spot(:)*cutpotu(:)
899 3 : n = basis%nbas(la)
900 51 : DO i = 1, n
901 459 : DO j = i, n
902 : hnl(i, j, la) = hnl(i, j, la) + &
903 : integrate_grid(spot(:), &
904 408 : gbasis%bf(:, i, la), gbasis%bf(:, j, la), gbasis%grid)
905 456 : hnl(j, i, la) = hnl(i, j, la)
906 : END DO
907 : END DO
908 5 : DEALLOCATE (spot)
909 : END DO
910 : ELSE
911 0 : CPABORT("Pseudopotential type: ["//ADJUSTL(TRIM(ptype))//"] not known")
912 : END IF
913 :
914 : ! release basis representation on UPF grid
915 26 : CALL release_atom_basis(gbasis)
916 :
917 520 : END SUBROUTINE upfints
918 :
919 : ! **************************************************************************************************
920 : !> \brief ...
921 : !> \param hnl ...
922 : !> \param basis ...
923 : !> \param ecp_pot ...
924 : ! **************************************************************************************************
925 3 : SUBROUTINE ecpints(hnl, basis, ecp_pot)
926 : REAL(KIND=dp), DIMENSION(:, :, 0:) :: hnl
927 : TYPE(atom_basis_type), INTENT(INOUT) :: basis
928 : TYPE(atom_ecppot_type) :: ecp_pot
929 :
930 : INTEGER :: i, j, k, l, m, n
931 : REAL(KIND=dp) :: alpha
932 3 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cpot
933 3 : REAL(KIND=dp), DIMENSION(:), POINTER :: rad
934 :
935 3 : rad => basis%grid%rad
936 3 : m = basis%grid%nr
937 9 : ALLOCATE (cpot(1:m))
938 :
939 : ! non local pseudopotential
940 4917 : hnl = 0.0_dp
941 14 : DO l = 0, ecp_pot%lmax
942 11 : cpot = 0._dp
943 26 : DO k = 1, ecp_pot%npot(l)
944 15 : n = ecp_pot%nrpot(k, l)
945 15 : alpha = ecp_pot%bpot(k, l)
946 6026 : cpot(:) = cpot(:) + ecp_pot%apot(k, l)*rad(:)**(n - 2)*EXP(-alpha*rad(:)**2)
947 : END DO
948 190 : DO i = 1, basis%nbas(l)
949 1683 : DO j = i, basis%nbas(l)
950 1496 : hnl(i, j, l) = integrate_grid(cpot, basis%bf(:, i, l), basis%bf(:, j, l), basis%grid)
951 1672 : hnl(j, i, l) = hnl(i, j, l)
952 : END DO
953 : END DO
954 : END DO
955 3 : DEALLOCATE (cpot)
956 :
957 3 : END SUBROUTINE ecpints
958 :
959 : ! **************************************************************************************************
960 : !> \brief ...
961 : !> \param core ...
962 : !> \param hnl ...
963 : !> \param basis ...
964 : !> \param sgp_pot ...
965 : !> \param cutpots ...
966 : ! **************************************************************************************************
967 29 : SUBROUTINE sgpints(core, hnl, basis, sgp_pot, cutpots)
968 : REAL(KIND=dp), DIMENSION(:, :, 0:) :: core, hnl
969 : TYPE(atom_basis_type), INTENT(INOUT) :: basis
970 : TYPE(atom_sgp_potential_type) :: sgp_pot
971 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: cutpots
972 :
973 : INTEGER :: i, ia, j, l, m, n, na
974 : REAL(KIND=dp) :: a, zval
975 29 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cpot, pgauss
976 29 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: qmat
977 29 : REAL(KIND=dp), DIMENSION(:), POINTER :: rad
978 :
979 29 : rad => basis%grid%rad
980 29 : m = basis%grid%nr
981 :
982 : ! local pseudopotential
983 87 : ALLOCATE (cpot(m))
984 29 : IF (sgp_pot%has_local) THEN
985 : zval = sgp_pot%zval
986 63926 : core = 0._dp
987 26 : cpot = 0.0_dp
988 338 : DO i = 1, sgp_pot%n_local
989 153938 : cpot(:) = cpot(:) + sgp_pot%c_local(i)*EXP(-sgp_pot%a_local(i)*rad(:)**2)
990 : END DO
991 12826 : cpot(:) = cpot(:)*cutpots(:)
992 26 : CALL numpot_matrix(core, cpot, basis, 0)
993 : END IF
994 29 : DEALLOCATE (cpot)
995 :
996 : ! nonlocal pseudopotential
997 29 : IF (sgp_pot%has_nonlocal) THEN
998 38519 : hnl = 0.0_dp
999 51 : ALLOCATE (pgauss(1:m))
1000 17 : n = sgp_pot%n_nonlocal
1001 : !
1002 57 : DO l = 0, sgp_pot%lmax
1003 40 : CPASSERT(l <= UBOUND(basis%nbas, 1))
1004 40 : IF (.NOT. sgp_pot%is_nonlocal(l)) CYCLE
1005 : ! overlap (a|p)
1006 39 : na = basis%nbas(l)
1007 156 : ALLOCATE (qmat(na, n))
1008 351 : DO i = 1, n
1009 312 : a = sgp_pot%a_nonlocal(i)
1010 144312 : pgauss(:) = cutpots(:)*EXP(-a*rad(:)**2)*rad(:)**l
1011 6111 : DO ia = 1, na
1012 6072 : qmat(ia, i) = integrate_grid(basis%bf(:, ia, l), pgauss(:), basis%grid)
1013 : END DO
1014 : END DO
1015 55116 : qmat(1:na, 1:n) = SQRT(fourpi)*MATMUL(qmat(1:na, 1:n), sgp_pot%c_nonlocal(1:n, 1:n, l))
1016 759 : DO i = 1, na
1017 7839 : DO j = i, na
1018 63720 : DO ia = 1, n
1019 63720 : hnl(i, j, l) = hnl(i, j, l) + qmat(i, ia)*qmat(j, ia)*sgp_pot%h_nonlocal(ia, l)
1020 : END DO
1021 7800 : hnl(j, i, l) = hnl(i, j, l)
1022 : END DO
1023 : END DO
1024 56 : DEALLOCATE (qmat)
1025 : END DO
1026 17 : DEALLOCATE (pgauss)
1027 : END IF
1028 :
1029 58 : END SUBROUTINE sgpints
1030 :
1031 : ! **************************************************************************************************
1032 : !> \brief ...
1033 : !> \param ac ...
1034 : !> \param vlocal ...
1035 : !> \param r ...
1036 : !> \param z ...
1037 : ! **************************************************************************************************
1038 26 : SUBROUTINE erffit(ac, vlocal, r, z)
1039 : REAL(KIND=dp), INTENT(OUT) :: ac
1040 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: vlocal, r
1041 : REAL(KIND=dp), INTENT(IN) :: z
1042 :
1043 : REAL(KIND=dp), PARAMETER :: rcut = 1.4_dp
1044 :
1045 : INTEGER :: i, j, m, m1
1046 : REAL(KIND=dp) :: a1, a2, an, e2, en, rc
1047 26 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: epot, rval, vpot
1048 :
1049 26 : m = SIZE(r)
1050 130 : ALLOCATE (epot(m), vpot(m), rval(m))
1051 26 : CPASSERT(SIZE(vlocal) == m)
1052 26 : IF (r(1) > r(m)) THEN
1053 0 : DO i = 1, m
1054 0 : vpot(m - i + 1) = vlocal(i)
1055 0 : rval(m - i + 1) = r(i)
1056 : END DO
1057 : ELSE
1058 25810 : vpot(1:m) = vlocal(1:m)
1059 25810 : rval(1:m) = r(1:m)
1060 : END IF
1061 26 : m1 = 1
1062 13421 : DO i = 1, m
1063 13421 : IF (rval(i) > rcut) THEN
1064 : m1 = i
1065 : EXIT
1066 : END IF
1067 : END DO
1068 :
1069 26 : a1 = 0.2_dp
1070 26 : a2 = 0.2_dp
1071 26 : e2 = 1.e20_dp
1072 26 : epot = 0.0_dp
1073 572 : DO i = 0, 20
1074 546 : an = a1 + i*0.025_dp
1075 546 : rc = 1._dp/(an*SQRT(2.0_dp))
1076 260715 : DO j = m1, m
1077 260715 : epot(j) = vpot(j) + z/rval(j)*erf(rval(j)*rc)
1078 : END DO
1079 260715 : en = SUM(ABS(epot(m1:m)*rval(m1:m)**2))
1080 572 : IF (en < e2) THEN
1081 145 : e2 = en
1082 145 : a2 = an
1083 : END IF
1084 : END DO
1085 26 : ac = a2
1086 :
1087 26 : DEALLOCATE (epot, vpot, rval)
1088 :
1089 26 : END SUBROUTINE erffit
1090 :
1091 : ! **************************************************************************************************
1092 : !> \brief ...
1093 : !> \param nl ...
1094 : !> \param al ...
1095 : !> \param cl ...
1096 : !> \param vloc ...
1097 : !> \param vgauss ...
1098 : !> \param gbasis ...
1099 : !> \param rad ...
1100 : !> \param ww ...
1101 : !> \param method ...
1102 : !> \param errloc ...
1103 : ! **************************************************************************************************
1104 2247 : SUBROUTINE pplocal_error(nl, al, cl, vloc, vgauss, gbasis, rad, ww, method, errloc)
1105 : INTEGER :: nl
1106 : REAL(KIND=dp), DIMENSION(:) :: al, cl, vloc, vgauss
1107 : TYPE(atom_basis_type) :: gbasis
1108 : REAL(KIND=dp), DIMENSION(:) :: rad, ww
1109 : INTEGER, INTENT(IN) :: method
1110 : REAL(KIND=dp) :: errloc
1111 :
1112 : INTEGER :: ia, ib, ir, ix, la, na
1113 2247 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: tv
1114 2247 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: rmat, sinv, smat
1115 2247 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: gmat
1116 :
1117 29211 : cl = 0.0_dp
1118 2247 : IF (method == 1) THEN
1119 17976 : ALLOCATE (tv(nl), smat(nl, nl), sinv(nl, nl))
1120 29211 : DO ir = 1, nl
1121 26446968 : vgauss(:) = EXP(-al(ir)*rad(:)**2)
1122 350532 : DO ix = 1, nl
1123 317390580 : smat(ir, ix) = SUM(vgauss(:)*EXP(-al(ix)*rad(:)**2)*ww(:))
1124 : END DO
1125 26449215 : tv(ir) = SUM(vloc(:)*vgauss(:)*ww(:))
1126 : END DO
1127 2247 : CALL get_pseudo_inverse_diag(smat(1:nl, 1:nl), sinv(1:nl, 1:nl), 1.e-12_dp)
1128 379743 : cl(1:nl) = MATMUL(sinv(1:nl, 1:nl), tv(1:nl))
1129 : ELSE
1130 : !
1131 0 : ALLOCATE (tv(nl), smat(nl, nl), sinv(nl, nl))
1132 : !
1133 0 : smat = 0.0_dp
1134 0 : tv = 0.0_dp
1135 0 : DO la = 0, MIN(UBOUND(gbasis%nbas, 1), 3)
1136 0 : na = gbasis%nbas(la)
1137 0 : ALLOCATE (rmat(na, na), gmat(na, na, nl))
1138 0 : rmat = 0.0_dp
1139 0 : gmat = 0.0_dp
1140 0 : DO ia = 1, na
1141 0 : DO ib = ia, na
1142 0 : rmat(ia, ib) = SUM(gbasis%bf(:, ia, la)*gbasis%bf(:, ib, la)*vloc(:)*ww(:))
1143 0 : rmat(ib, ia) = rmat(ia, ib)
1144 : END DO
1145 : END DO
1146 0 : DO ir = 1, nl
1147 0 : vgauss(:) = EXP(-al(ir)*rad(:)**2)
1148 0 : DO ia = 1, na
1149 0 : DO ib = ia, na
1150 0 : gmat(ia, ib, ir) = SUM(gbasis%bf(:, ia, la)*gbasis%bf(:, ib, la)*vgauss(:)*ww(:))
1151 0 : gmat(ib, ia, ir) = gmat(ia, ib, ir)
1152 : END DO
1153 : END DO
1154 : END DO
1155 0 : DO ir = 1, nl
1156 0 : tv(ir) = tv(ir) + accurate_dot_product(rmat, gmat(:, :, ir))
1157 0 : DO ix = ir, nl
1158 0 : smat(ir, ix) = smat(ir, ix) + accurate_dot_product(gmat(:, :, ix), gmat(:, :, ir))
1159 0 : smat(ix, ir) = smat(ir, ix)
1160 : END DO
1161 : END DO
1162 0 : DEALLOCATE (rmat, gmat)
1163 : END DO
1164 0 : sinv = 0.0_dp
1165 0 : CALL get_pseudo_inverse_diag(smat(1:nl, 1:nl), sinv(1:nl, 1:nl), 1.e-12_dp)
1166 0 : cl(1:nl) = MATMUL(sinv(1:nl, 1:nl), tv(1:nl))
1167 : END IF
1168 : !
1169 2203914 : vgauss = 0.0_dp
1170 29211 : DO ir = 1, nl
1171 26449215 : vgauss(:) = vgauss(:) + cl(ir)*EXP(-al(ir)*rad(:)**2)
1172 : END DO
1173 2203914 : errloc = SUM((vgauss - vloc)**2*ww)
1174 : !
1175 2247 : DEALLOCATE (tv, smat, sinv)
1176 : !
1177 2247 : END SUBROUTINE pplocal_error
1178 :
1179 : ! **************************************************************************************************
1180 : !> \brief ...
1181 : !> \param pot ...
1182 : !> \param r ...
1183 : !> \param rcut ...
1184 : !> \param rsmooth ...
1185 : ! **************************************************************************************************
1186 52 : SUBROUTINE cutpot(pot, r, rcut, rsmooth)
1187 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: pot
1188 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: r
1189 : REAL(KIND=dp), INTENT(IN) :: rcut, rsmooth
1190 :
1191 : INTEGER :: i, n
1192 : REAL(KIND=dp) :: rab, rx, x
1193 :
1194 52 : n = SIZE(pot)
1195 52 : CPASSERT(n <= SIZE(r))
1196 :
1197 38636 : pot(:) = 1.0_dp
1198 38636 : DO i = 1, n
1199 38584 : rab = r(i)
1200 38636 : IF (rab > rcut) THEN
1201 300 : pot(i) = 0.0_dp
1202 38284 : ELSE IF (rab > rcut - rsmooth) THEN
1203 1685 : rx = rab - (rcut - rsmooth)
1204 1685 : x = rx/rsmooth
1205 1685 : pot(i) = -6._dp*x**5 + 15._dp*x**4 - 10._dp*x**3 + 1._dp
1206 : END IF
1207 : END DO
1208 :
1209 52 : END SUBROUTINE cutpot
1210 :
1211 117 : END MODULE atom_sgp
|
