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