(git:374b731)
Loading...
Searching...
No Matches
atom_sgp.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2024 CP2K developers group <https://cp2k.org> !
4! !
5! SPDX-License-Identifier: GPL-2.0-or-later !
6!--------------------------------------------------------------------------------------------------!
7
8! **************************************************************************************************
9MODULE 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
68CONTAINS
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 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 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 IF (PRESENT(orb_basis)) THEN
96 CALL set_kind_basis_atomic(basis, orb_basis, has_pp=.true., cp2k_norm=.true.)
97 ELSE
98 CALL init_atom_basis_default_pp(basis)
99 END IF
100
101 is_ecp = .false.
102 IF (PRESENT(ecp_pot)) is_ecp = .true.
103 is_upf = .false.
104 IF (PRESENT(upf_pot)) is_upf = .true.
105 cpassert(.NOT. (is_ecp .AND. is_upf))
106
107 ! upf has often very small grids, use a smooth cutoff function
108 IF (is_upf) THEN
109 n = SIZE(upf_pot%r)
110 ALLOCATE (cutpotu(n))
111 rcut = maxval(upf_pot%r)
112 CALL cutpot(cutpotu, upf_pot%r, rcut, 2.5_dp)
113 n = basis%grid%nr
114 ALLOCATE (cutpots(n))
115 CALL cutpot(cutpots, basis%grid%rad, rcut, 2.5_dp)
116 ELSE
117 n = basis%grid%nr
118 ALLOCATE (cutpots(n))
119 cutpots = 1.0_dp
120 END IF
121
122 ! generate the transformed potentials
123 IF (is_ecp) THEN
124 CALL ecp_sgp_constr(ecp_pot, sgp_pot, basis)
125 ELSEIF (is_upf) THEN
126 CALL upf_sgp_constr(upf_pot, sgp_pot, basis)
127 ELSE
128 cpabort("")
129 END IF
130
131 NULLIFY (core, hnl)
132 CALL create_opmat(core, basis%nbas)
133 CALL create_opmat(hnl, basis%nbas, 5)
134 NULLIFY (score, shnl)
135 CALL create_opmat(score, basis%nbas)
136 CALL create_opmat(shnl, basis%nbas, 5)
137 !
138 IF (is_ecp) THEN
139 CALL ecpints(hnl%op, basis, ecp_pot)
140 ELSEIF (is_upf) THEN
141 CALL upfints(core%op, hnl%op, basis, upf_pot, cutpotu, sgp_pot%ac_local)
142 ELSE
143 cpabort("")
144 END IF
145 !
146 CALL sgpints(score%op, shnl%op, basis, sgp_pot, cutpots)
147 !
148 error = 0.0_dp
149 IF (sgp_pot%has_local) THEN
150 n = min(3, ubound(core%op, 3))
151 error(1) = maxval(abs(core%op(:, :, 0:n) - score%op(:, :, 0:n)))
152 mloc = maxloc(abs(core%op(:, :, 0:n) - score%op(:, :, 0:n)))
153 error(4) = error(1)/max(abs(score%op(mloc(1), mloc(2), mloc(3))), 1.e-6_dp)
154 END IF
155 IF (sgp_pot%has_nonlocal) THEN
156 n = min(3, ubound(hnl%op, 3))
157 error(2) = maxval(abs(hnl%op(:, :, 0:n) - shnl%op(:, :, 0:n)))
158 mloc = maxloc(abs(hnl%op(:, :, 0:n) - shnl%op(:, :, 0:n)))
159 error(5) = error(1)/max(abs(hnl%op(mloc(1), mloc(2), mloc(3))), 1.e-6_dp)
160 END IF
161 IF (sgp_pot%has_nlcc) THEN
162 IF (is_upf) THEN
163 n = SIZE(upf_pot%r)
164 ALLOCATE (cgauss(n))
165 cgauss = 0.0_dp
166 DO i = 1, sgp_pot%n_nlcc
167 cgauss(:) = cgauss(:) + sgp_pot%c_nlcc(i)*exp(-sgp_pot%a_nlcc(i)*upf_pot%r(:)**2)
168 END DO
169 errcc = sum((cgauss(:) - upf_pot%rho_nlcc(:))**2*upf_pot%r(:)**2*upf_pot%rab(:))
170 errcc = sqrt(errcc/real(n, kind=dp))
171 DEALLOCATE (cgauss)
172 ELSE
173 cpabort("")
174 END IF
175 error(3) = errcc
176 END IF
177 !
178 IF (is_upf) THEN
179 DEALLOCATE (cutpotu)
180 DEALLOCATE (cutpots)
181 ELSE
182 DEALLOCATE (cutpots)
183 END IF
184 !
185 CALL release_opmat(score)
186 CALL release_opmat(shnl)
187 CALL release_opmat(core)
188 CALL release_opmat(hnl)
189
190 CALL release_atom_basis(basis)
191
192 END SUBROUTINE sgp_construction
193
194! **************************************************************************************************
195!> \brief ...
196!> \param atom_info ...
197!> \param input_section ...
198!> \param iw ...
199! **************************************************************************************************
200 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 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 CALL section_vals_get(input_section, explicit=explicit)
217 IF (.NOT. explicit) RETURN
218
219 IF (iw > 0) WRITE (iw, '(/," ",79("*"),/,T24,A,/," ",79("*"))') "SEPARABLE GAUSSIAN PSEUDOPOTENTIAL"
220
221 atom_ref => atom_info(1, 1)%atom
222
223 ppot_type = atom_ref%potential%ppot_type
224 SELECT CASE (ppot_type)
225 CASE (gth_pseudo)
226 IF (iw > 0) WRITE (iw, '(" GTH Pseudopotential is already in SGP form. ")')
227 do_transform = .false.
228 CASE (ecp_pseudo)
229 do_transform = .true.
230 is_ecp = .true.
231 is_upf = .false.
232 ecp_pot => atom_ref%potential%ecp_pot
233 CASE (upf_pseudo)
234 do_transform = .true.
235 is_ecp = .false.
236 is_upf = .true.
237 upf_pot => atom_ref%potential%upf_pot
238 CASE (no_pseudo)
239 IF (iw > 0) WRITE (iw, '(" No Pseudopotential available for transformation. ")')
240 do_transform = .false.
241 CASE DEFAULT
242 cpabort("")
243 END SELECT
244
245 ! generate the transformed potentials
246 IF (do_transform) THEN
247 IF (is_ecp) THEN
248 CALL ecp_sgp_constr(ecp_pot, sgp_pot, atom_ref%basis)
249 ELSEIF (is_upf) THEN
250 CALL upf_sgp_constr(upf_pot, sgp_pot, atom_ref%basis)
251 ELSE
252 cpabort("")
253 END IF
254 END IF
255
256 ! Check the result
257 IF (do_transform) THEN
258 NULLIFY (core, hnl)
259 CALL create_opmat(core, atom_ref%basis%nbas)
260 CALL create_opmat(hnl, atom_ref%basis%nbas, 5)
261 NULLIFY (score, shnl)
262 CALL create_opmat(score, atom_ref%basis%nbas)
263 CALL create_opmat(shnl, atom_ref%basis%nbas, 5)
264 !
265 ! upf has often very small grids, use a smooth cutoff function
266 IF (is_upf) THEN
267 n = SIZE(upf_pot%r)
268 ALLOCATE (cutpotu(n))
269 rcut = maxval(upf_pot%r)
270 CALL cutpot(cutpotu, upf_pot%r, rcut, 2.5_dp)
271 n = atom_ref%basis%grid%nr
272 ALLOCATE (cutpots(n))
273 CALL cutpot(cutpots, atom_ref%basis%grid%rad, rcut, 2.5_dp)
274 ELSE
275 n = atom_ref%basis%grid%nr
276 ALLOCATE (cutpots(n))
277 cutpots = 1.0_dp
278 END IF
279 !
280 IF (is_ecp) THEN
281 CALL ecpints(hnl%op, atom_ref%basis, ecp_pot)
282 ELSEIF (is_upf) THEN
283 CALL upfints(core%op, hnl%op, atom_ref%basis, upf_pot, cutpotu, sgp_pot%ac_local)
284 ELSE
285 cpabort("")
286 END IF
287 !
288 CALL sgpints(score%op, shnl%op, atom_ref%basis, sgp_pot, cutpots)
289 !
290 IF (sgp_pot%has_local) THEN
291 n = min(3, ubound(core%op, 3))
292 errcc = maxval(abs(core%op(:, :, 0:n) - score%op(:, :, 0:n)))
293 IF (iw > 0) THEN
294 WRITE (iw, '(" Local part of pseudopotential")')
295 WRITE (iw, '(" Number of basis functions ",T77,i4)') sgp_pot%n_local
296 WRITE (iw, '(" Max. abs. error of matrix elements ",T65,f16.8)') errcc
297 END IF
298 END IF
299 IF (sgp_pot%has_nonlocal) THEN
300 IF (iw > 0) THEN
301 WRITE (iw, '(" Nonlocal part of pseudopotential")')
302 WRITE (iw, '(" Max. l-quantum number",T77,i4)') sgp_pot%lmax
303 WRITE (iw, '(" Number of projector basis functions ",T77,i4)') sgp_pot%n_nonlocal
304 DO i = 0, sgp_pot%lmax
305 errcc = maxval(abs(hnl%op(:, :, i) - shnl%op(:, :, i)))
306 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 IF (sgp_pot%has_nlcc) THEN
311 IF (is_upf) THEN
312 n = SIZE(upf_pot%r)
313 ALLOCATE (cgauss(n))
314 cgauss = 0.0_dp
315 DO i = 1, sgp_pot%n_nlcc
316 cgauss(:) = cgauss(:) + sgp_pot%c_nlcc(i)*exp(-sgp_pot%a_nlcc(i)*upf_pot%r(:)**2)
317 END DO
318 errcc = sum((cgauss(:) - upf_pot%rho_nlcc(:))**2*upf_pot%r(:)**2*upf_pot%rab(:))
319 errcc = sqrt(errcc/real(n, kind=dp))
320 DEALLOCATE (cgauss)
321 ELSE
322 cpabort("")
323 END IF
324 IF (iw > 0) THEN
325 WRITE (iw, '(" Non-linear core correction: core density")')
326 WRITE (iw, '(" Number of basis functions ",T77,i4)') sgp_pot%n_nlcc
327 WRITE (iw, '(" RMS error of core density ",T69,f12.8)') errcc
328 END IF
329 END IF
330 !
331 IF (is_upf) THEN
332 DEALLOCATE (cutpotu)
333 DEALLOCATE (cutpots)
334 ELSE
335 DEALLOCATE (cutpots)
336 END IF
337 !
338 CALL release_opmat(score)
339 CALL release_opmat(shnl)
340 CALL release_opmat(core)
341 CALL release_opmat(hnl)
342 END IF
343
344 CALL atom_sgp_release(sgp_pot)
345
346 IF (iw > 0) WRITE (iw, '(" ",79("*"))')
347
348 END SUBROUTINE atom_sgp_construction
349
350! **************************************************************************************************
351!> \brief ...
352!> \param ecp_pot ...
353!> \param sgp_pot ...
354!> \param basis ...
355! **************************************************************************************************
356 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: al, cl, cpot, pgauss
365 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: cmat, qmat, score, sinv, smat, tmat
366 REAL(kind=dp), DIMENSION(:), POINTER :: rad
367
368 sgp_pot%has_nlcc = .false.
369 sgp_pot%n_nlcc = 0
370 sgp_pot%has_local = .false.
371 sgp_pot%n_local = 0
372
373 ! transform semilocal potential into a separable local form
374 sgp_pot%has_nonlocal = .true.
375 !
376 ! hardcoded number of projectors (equal for all l values)
377 nl = 8
378 !
379 sgp_pot%n_nonlocal = nl
380 sgp_pot%lmax = ecp_pot%lmax
381 ALLOCATE (sgp_pot%a_nonlocal(nl))
382 ALLOCATE (sgp_pot%h_nonlocal(nl, 0:ecp_pot%lmax))
383 ALLOCATE (sgp_pot%c_nonlocal(nl, nl, 0:ecp_pot%lmax))
384 !
385 ALLOCATE (al(nl), cl(nl))
386 ALLOCATE (smat(nl, nl), sinv(nl, nl))
387 ALLOCATE (tmat(nl, nl), cmat(nl, nl))
388 al = 0.0_dp
389 amx = maxval(ecp_pot%bpot)
390 cmx = amx/0.15_dp
391 cmx = cmx**(1._dp/real(nl - 1, dp))
392 cmx = 1._dp/cmx
393 DO ir = 1, nl
394 al(ir) = amx*cmx**(ir - 1)
395 END DO
396 !
397 sgp_pot%a_nonlocal(1:nl) = al(1:nl)
398 !
399 nr = basis%grid%nr
400 rad => basis%grid%rad
401 ALLOCATE (cpot(nr), pgauss(nr))
402 DO l = 0, ecp_pot%lmax
403 na = basis%nbas(l)
404 ALLOCATE (score(na, na), qmat(na, nl))
405 cpot = 0._dp
406 DO k = 1, ecp_pot%npot(l)
407 n = ecp_pot%nrpot(k, l)
408 alpha = ecp_pot%bpot(k, l)
409 cpot(:) = cpot + ecp_pot%apot(k, l)*rad**(n - 2)*exp(-alpha*rad**2)
410 END DO
411 DO i = 1, na
412 DO j = i, na
413 score(i, j) = integrate_grid(cpot, basis%bf(:, i, l), basis%bf(:, j, l), basis%grid)
414 score(j, i) = score(i, j)
415 END DO
416 END DO
417 ! overlap basis with projectors
418 DO i = 1, nl
419 pgauss(:) = exp(-al(i)*rad(:)**2)*rad(:)**l
420 DO ia = 1, na
421 qmat(ia, i) = integrate_grid(basis%bf(:, ia, l), pgauss(:), basis%grid)
422 END DO
423 END DO
424 qmat = sqrt(fourpi)*qmat
425 ! tmat = qmat * score * qmat
426 tmat(1:nl, 1:nl) = matmul(transpose(qmat(1:na, 1:nl)), matmul(score(1:na, 1:na), qmat(1:na, 1:nl)))
427 smat(1:nl, 1:nl) = matmul(transpose(qmat(1:na, 1:nl)), qmat(1:na, 1:nl))
428 CALL get_pseudo_inverse_diag(smat(1:nl, 1:nl), sinv(1:nl, 1:nl), 1.e-10_dp)
429 cmat(1:nl, 1:nl) = matmul(sinv(1:nl, 1:nl), matmul(tmat(1:nl, 1:nl), sinv(1:nl, 1:nl)))
430 cmat(1:nl, 1:nl) = (cmat(1:nl, 1:nl) + transpose(cmat(1:nl, 1:nl)))*0.5_dp
431 !
432 ! Residium
433 !
434 CALL diamat_all(cmat(1:nl, 1:nl), cl(1:nl), .true.)
435 !
436 sgp_pot%h_nonlocal(1:nl, l) = cl(1:nl)
437 sgp_pot%c_nonlocal(1:nl, 1:nl, l) = cmat(1:nl, 1:nl)
438 sgp_pot%is_nonlocal(l) = .true.
439 !
440 DEALLOCATE (score, qmat)
441 END DO
442 DEALLOCATE (cpot, pgauss)
443 DEALLOCATE (al, cl, smat, sinv, tmat, cmat)
444
445 END SUBROUTINE ecp_sgp_constr
446
447! **************************************************************************************************
448!> \brief ...
449!> \param upf_pot ...
450!> \param sgp_pot ...
451!> \param basis ...
452! **************************************************************************************************
453 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: al, ccharge, cgauss, cl, pgauss, pproa, &
466 pprob, tv, vgauss, vloc, ww
467 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 IF (upf_pot%is_ultrasoft .OR. upf_pot%is_paw .OR. upf_pot%is_coulomb) THEN
472 sgp_pot%has_nonlocal = .false.
473 sgp_pot%n_nonlocal = 0
474 sgp_pot%has_local = .false.
475 sgp_pot%n_local = 0
476 sgp_pot%has_nlcc = .false.
477 sgp_pot%n_nlcc = 0
478 RETURN
479 END IF
480
481 ! radial grid
482 nr = SIZE(upf_pot%r)
483 ! weights for integration
484 ALLOCATE (ww(nr))
485 ww(:) = upf_pot%r(:)**2*upf_pot%rab(:)
486
487 ! start with local potential
488 sgp_pot%has_local = .true.
489 ! fit local potential to Gaussian form
490 ALLOCATE (vloc(nr), vgauss(nr))
491 ! smearing of core charge
492 zval = upf_pot%zion
493 ! Try to find an optimal Gaussian charge distribution
494 CALL erffit(sgp_pot%ac_local, upf_pot%vlocal, upf_pot%r, zval)
495 sgp_pot%zval = zval
496 DO ir = 1, nr
497 IF (upf_pot%r(ir) < 1.e-12_dp) THEN
498 vgauss(ir) = -sqrt(2.0_dp)*zval/rootpi/sgp_pot%ac_local
499 ELSE
500 rc = upf_pot%r(ir)/sgp_pot%ac_local/sqrt(2.0_dp)
501 vgauss(ir) = -zval/upf_pot%r(ir)*erf(rc)
502 END IF
503 END DO
504 vloc(:) = upf_pot%vlocal(:) - vgauss(:)
505 !
506 CALL atom_basis_gridrep(basis, gbasis, upf_pot%r, upf_pot%rab)
507 !
508 nl = 12
509 ALLOCATE (al(nl), cl(nl))
510 ostate%nf = 0
511 ostate%nvar = 2
512 x(1) = 1.00_dp !starting point of geometric series
513 x(2) = 1.20_dp !factor of series
514 ostate%rhoend = 1.e-12_dp
515 ostate%rhobeg = 5.e-2_dp
516 ostate%maxfun = 1000
517 ostate%iprint = 1
518 ostate%unit = -1
519 ostate%state = 0
520 DO
521 IF (ostate%state == 2) THEN
522 DO ir = 1, nl
523 al(ir) = x(1)*x(2)**(ir - 1)
524 END DO
525 CALL pplocal_error(nl, al, cl, vloc, vgauss, gbasis, upf_pot%r, ww, 1, ostate%f)
526 END IF
527 IF (ostate%state == -1) EXIT
528 CALL powell_optimize(ostate%nvar, x, ostate)
529 END DO
530 ostate%state = 8
531 CALL powell_optimize(ostate%nvar, x, ostate)
532 DO ir = 1, nl
533 al(ir) = x(1)*x(2)**(ir - 1)
534 END DO
535 CALL pplocal_error(nl, al, cl, vloc, vgauss, gbasis, upf_pot%r, ww, 1, errloc)
536 !
537 ALLOCATE (sgp_pot%a_local(nl), sgp_pot%c_local(nl))
538 sgp_pot%n_local = nl
539 sgp_pot%a_local(1:nl) = al(1:nl)
540 sgp_pot%c_local(1:nl) = cl(1:nl)
541 DEALLOCATE (vloc, vgauss)
542 DEALLOCATE (al, cl)
543 CALL release_atom_basis(gbasis)
544 !
545 ptype = adjustl(trim(upf_pot%pseudo_type))
546 IF (ptype(1:2) == "NC" .OR. ptype(1:2) == "US") THEN
547 nl_trans = .false.
548 ELSE IF (ptype(1:2) == "SL") THEN
549 nl_trans = .true.
550 ELSE
551 cpabort("Pseudopotential type: ["//adjustl(trim(ptype))//"] not known")
552 END IF
553
554 ! purely local pseudopotentials
555 IF (upf_pot%l_max < 0) THEN
556 sgp_pot%n_nonlocal = 0
557 sgp_pot%lmax = -1
558 sgp_pot%has_nonlocal = .false.
559 ELSE
560 ! Non-local pseudopotential in Gaussian form
561 IF (nl_trans) THEN
562 sgp_pot%has_nonlocal = .true.
563 ! semi local pseudopotential
564 ! fit to nonlocal form
565 ! get basis representation on UPF grid
566 nl = 8
567 !
568 sgp_pot%n_nonlocal = nl
569 sgp_pot%lmax = upf_pot%l_max
570 ALLOCATE (sgp_pot%a_nonlocal(nl))
571 ALLOCATE (sgp_pot%h_nonlocal(nl, 0:upf_pot%l_max))
572 ALLOCATE (sgp_pot%c_nonlocal(nl, nl, 0:upf_pot%l_max))
573 !
574 ALLOCATE (al(nl), cl(nl))
575 ALLOCATE (smat(nl, nl), sinv(nl, nl))
576 ALLOCATE (tmat(nl, nl), cmat(nl, nl))
577 al = 0.0_dp
578 DO ir = 1, nl
579 al(ir) = 10.0_dp*0.60_dp**(ir - 1)
580 END DO
581 !
582 sgp_pot%a_nonlocal(1:nl) = al(1:nl)
583 !
584 CALL atom_basis_gridrep(basis, gbasis, upf_pot%r, upf_pot%rab)
585 ALLOCATE (pgauss(nr), vloc(nr))
586 DO la = 0, upf_pot%l_max
587 IF (la == upf_pot%l_local) cycle
588 sgp_pot%is_nonlocal(la) = .true.
589 na = gbasis%nbas(la)
590 ALLOCATE (score(na, na), qmat(na, nl))
591 ! Reference matrix
592 vloc(:) = upf_pot%vsemi(:, la + 1) - upf_pot%vlocal(:)
593 DO ia = 1, na
594 DO ib = ia, na
595 score(ia, ib) = sum(vloc(:)*gbasis%bf(:, ia, la)*gbasis%bf(:, ib, la)*ww(:))
596 score(ib, ia) = score(ia, ib)
597 END DO
598 END DO
599 ! overlap basis with projectors
600 DO ir = 1, nl
601 pgauss(:) = exp(-al(ir)*upf_pot%r(:)**2)*upf_pot%r(:)**la
602 eee = rootpi/(2._dp**(la + 2)*dfac(2*la + 1))/(2._dp*al(ir))**(la + 1.5_dp)
603 pgauss(:) = pgauss(:)/sqrt(eee)
604 DO ia = 1, na
605 qmat(ia, ir) = sum(gbasis%bf(:, ia, la)*pgauss(:)*ww)
606 END DO
607 END DO
608 ! tmat = qmat * score * qmat
609 tmat(1:nl, 1:nl) = matmul(transpose(qmat(1:na, 1:nl)), matmul(score(1:na, 1:na), qmat(1:na, 1:nl)))
610 smat(1:nl, 1:nl) = matmul(transpose(qmat(1:na, 1:nl)), qmat(1:na, 1:nl))
611 CALL get_pseudo_inverse_diag(smat(1:nl, 1:nl), sinv(1:nl, 1:nl), 1.e-10_dp)
612 cmat(1:nl, 1:nl) = matmul(sinv(1:nl, 1:nl), matmul(tmat(1:nl, 1:nl), sinv(1:nl, 1:nl)))
613 cmat(1:nl, 1:nl) = (cmat(1:nl, 1:nl) + transpose(cmat(1:nl, 1:nl)))*0.5_dp
614 CALL diamat_all(cmat(1:nl, 1:nl), cl(1:nl), .true.)
615 !
616 ! get back unnormalized Gaussians
617 DO ir = 1, nl
618 ei = rootpi/(2._dp**(la + 2)*dfac(2*la + 1))/(2._dp*al(ir))**(la + 1.5_dp)
619 cmat(ir, 1:nl) = cmat(ir, 1:nl)/sqrt(ei)
620 END DO
621 sgp_pot%h_nonlocal(1:nl, la) = cl(1:nl)
622 sgp_pot%c_nonlocal(1:nl, 1:nl, la) = cmat(1:nl, 1:nl)
623 sgp_pot%is_nonlocal(la) = .true.
624 DEALLOCATE (score, qmat)
625 END DO
626 ! SQRT(4PI)
627 sgp_pot%c_nonlocal = sgp_pot%c_nonlocal/sqrt(fourpi)
628 CALL release_atom_basis(gbasis)
629 DEALLOCATE (pgauss, vloc)
630 DEALLOCATE (al, cl, smat, sinv, tmat, cmat)
631 ELSE
632 sgp_pot%has_nonlocal = .true.
633 ! non local pseudopotential
634 ALLOCATE (pproa(nr), pprob(nr), pgauss(nr))
635 np = upf_pot%number_of_proj
636 nl = 8
637 ALLOCATE (al(nl), cl(nl))
638 ALLOCATE (smat(nl, nl), sinv(nl, nl))
639 ALLOCATE (tmat(nl, nl), cmat(nl, nl))
640 al = 0.0_dp
641 cl = 0.0_dp
642 DO ir = 1, nl
643 al(ir) = 10.0_dp*0.60_dp**(ir - 1)
644 END DO
645 !
646 sgp_pot%lmax = maxval(upf_pot%lbeta(:))
647 sgp_pot%n_nonlocal = nl
648 ALLOCATE (sgp_pot%a_nonlocal(nl))
649 ALLOCATE (sgp_pot%h_nonlocal(nl, 0:sgp_pot%lmax))
650 ALLOCATE (sgp_pot%c_nonlocal(nl, nl, 0:sgp_pot%lmax))
651 !
652 sgp_pot%a_nonlocal(1:nl) = al(1:nl)
653 !
654 CALL atom_basis_gridrep(basis, gbasis, upf_pot%r, upf_pot%rab)
655 DO la = 0, sgp_pot%lmax
656 sgp_pot%is_nonlocal(la) = .true.
657 na = gbasis%nbas(la)
658 ALLOCATE (score(na, na), qmat(na, nl))
659 ! Reference matrix
660 score = 0.0_dp
661 DO ipa = 1, np
662 lb = upf_pot%lbeta(ipa)
663 IF (la /= lb) cycle
664 pproa(:) = upf_pot%beta(:, ipa)
665 DO ipb = 1, np
666 lb = upf_pot%lbeta(ipb)
667 IF (la /= lb) cycle
668 pprob(:) = upf_pot%beta(:, ipb)
669 eee = upf_pot%dion(ipa, ipb)
670 DO ia = 1, na
671 cpa = sum(pproa(:)*gbasis%bf(:, ia, la)*ww(:))
672 DO ib = ia, na
673 cpb = sum(pprob(:)*gbasis%bf(:, ib, la)*ww(:))
674 score(ia, ib) = score(ia, ib) + cpa*eee*cpb
675 score(ib, ia) = score(ia, ib)
676 END DO
677 END DO
678 END DO
679 END DO
680 ! overlap basis with projectors
681 DO ir = 1, nl
682 pgauss(:) = exp(-al(ir)*upf_pot%r(:)**2)*upf_pot%r(:)**la
683 eee = rootpi/(2._dp**(la + 2)*dfac(2*la + 1))/(2._dp*al(ir))**(la + 1.5_dp)
684 pgauss(:) = pgauss(:)/sqrt(eee)
685 DO ia = 1, na
686 qmat(ia, ir) = sum(gbasis%bf(:, ia, la)*pgauss(:)*ww)
687 END DO
688 END DO
689 ! tmat = qmat * score * qmat
690 tmat(1:nl, 1:nl) = matmul(transpose(qmat(1:na, 1:nl)), matmul(score(1:na, 1:na), qmat(1:na, 1:nl)))
691 smat(1:nl, 1:nl) = matmul(transpose(qmat(1:na, 1:nl)), qmat(1:na, 1:nl))
692 CALL get_pseudo_inverse_diag(smat(1:nl, 1:nl), sinv(1:nl, 1:nl), 1.e-10_dp)
693 cmat(1:nl, 1:nl) = matmul(sinv(1:nl, 1:nl), matmul(tmat(1:nl, 1:nl), sinv(1:nl, 1:nl)))
694 cmat(1:nl, 1:nl) = (cmat(1:nl, 1:nl) + transpose(cmat(1:nl, 1:nl)))*0.5_dp
695 CALL diamat_all(cmat(1:nl, 1:nl), cl(1:nl), .true.)
696 !
697 ! get back unnormalized Gaussians
698 DO ir = 1, nl
699 ei = rootpi/(2._dp**(la + 2)*dfac(2*la + 1))/(2._dp*al(ir))**(la + 1.5_dp)
700 cmat(ir, 1:nl) = cmat(ir, 1:nl)/sqrt(ei)
701 END DO
702 sgp_pot%h_nonlocal(1:nl, la) = cl(1:nl)
703 sgp_pot%c_nonlocal(1:nl, 1:nl, la) = cmat(1:nl, 1:nl)
704 sgp_pot%is_nonlocal(la) = .true.
705 DEALLOCATE (score, qmat)
706 END DO
707 ! SQRT(4PI)
708 sgp_pot%c_nonlocal = sgp_pot%c_nonlocal/sqrt(fourpi)
709 CALL release_atom_basis(gbasis)
710 DEALLOCATE (pgauss, pproa, pprob)
711 DEALLOCATE (al, cl, smat, sinv, tmat, cmat)
712 END IF
713 END IF
714
715 IF (upf_pot%core_correction) THEN
716 sgp_pot%has_nlcc = .true.
717 ELSE
718 sgp_pot%has_nlcc = .false.
719 sgp_pot%n_nlcc = 0
720 END IF
721
722 ! fit core charge to Gaussian form
723 IF (sgp_pot%has_nlcc) THEN
724 ALLOCATE (ccharge(nr), cgauss(nr))
725 ccharge(:) = upf_pot%rho_nlcc(:)
726 nl = 8
727 ALLOCATE (al(nl), cl(nl), tv(nl))
728 ALLOCATE (smat(nl, nl), sinv(nl, nl))
729 al = 0.0_dp
730 cl = 0.0_dp
731 DO ir = 1, nl
732 al(ir) = 10.0_dp*0.6_dp**(ir - 1)
733 END DO
734 ! calculate integrals
735 smat = 0.0_dp
736 sinv = 0.0_dp
737 tv = 0.0_dp
738 CALL sg_overlap(smat(1:nl, 1:nl), 0, al(1:nl), al(1:nl))
739 DO ir = 1, nl
740 cgauss(:) = exp(-al(ir)*upf_pot%r(:)**2)
741 tv(ir) = sum(cgauss*ccharge*ww)
742 END DO
743 CALL get_pseudo_inverse_diag(smat(1:nl, 1:nl), sinv(1:nl, 1:nl), 1.e-10_dp)
744 cl(1:nl) = matmul(sinv(1:nl, 1:nl), tv(1:nl))
745 cgauss = 0.0_dp
746 DO ir = 1, nl
747 cgauss(:) = cgauss(:) + cl(ir)*exp(-al(ir)*upf_pot%r(:)**2)
748 END DO
749 errcc = sum((cgauss - ccharge)**2*ww)
750 ALLOCATE (sgp_pot%a_local(nl), sgp_pot%c_local(nl))
751 sgp_pot%n_nlcc = nl
752 sgp_pot%a_nlcc(1:nl) = al(1:nl)
753 sgp_pot%c_nlcc(1:nl) = cl(1:nl)
754 DEALLOCATE (ccharge, cgauss)
755 DEALLOCATE (al, cl, tv, smat, sinv)
756 END IF
757
758 DEALLOCATE (ww)
759
760 END SUBROUTINE upf_sgp_constr
761
762! **************************************************************************************************
763!> \brief ...
764!> \param sgp_pot ...
765! **************************************************************************************************
766 SUBROUTINE atom_sgp_release(sgp_pot)
767
768 TYPE(atom_sgp_potential_type) :: sgp_pot
769
770 IF (ASSOCIATED(sgp_pot%a_nonlocal)) DEALLOCATE (sgp_pot%a_nonlocal)
771 IF (ASSOCIATED(sgp_pot%h_nonlocal)) DEALLOCATE (sgp_pot%h_nonlocal)
772 IF (ASSOCIATED(sgp_pot%c_nonlocal)) DEALLOCATE (sgp_pot%c_nonlocal)
773
774 IF (ASSOCIATED(sgp_pot%a_local)) DEALLOCATE (sgp_pot%a_local)
775 IF (ASSOCIATED(sgp_pot%c_local)) DEALLOCATE (sgp_pot%c_local)
776
777 IF (ASSOCIATED(sgp_pot%a_nlcc)) DEALLOCATE (sgp_pot%a_nlcc)
778 IF (ASSOCIATED(sgp_pot%c_nlcc)) DEALLOCATE (sgp_pot%c_nlcc)
779
780 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: spot
802 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: spmat
803 TYPE(atom_basis_type) :: gbasis
804
805 ! get basis representation on UPF grid
806 CALL atom_basis_gridrep(basis, gbasis, upf_pot%r, upf_pot%rab)
807
808 ! local pseudopotential
809 core = 0._dp
810 n = SIZE(upf_pot%r)
811 ALLOCATE (spot(n))
812 spot(:) = upf_pot%vlocal(:)
813 zval = upf_pot%zion
814 DO i = 1, n
815 IF (upf_pot%r(i) < 1.e-12_dp) THEN
816 spot(i) = spot(i) + sqrt2*zval/rootpi/ac_local
817 ELSE
818 rc = upf_pot%r(i)/ac_local/sqrt2
819 spot(i) = spot(i) + zval/upf_pot%r(i)*erf(rc)
820 END IF
821 END DO
822 spot(:) = spot(:)*cutpotu(:)
823
824 CALL numpot_matrix(core, spot, gbasis, 0)
825 DEALLOCATE (spot)
826
827 hnl = 0._dp
828 ptype = adjustl(trim(upf_pot%pseudo_type))
829 IF (ptype(1:2) == "NC" .OR. ptype(1:2) == "US") THEN
830 ! non local pseudopotential
831 n = maxval(gbasis%nbas(:))
832 m = upf_pot%number_of_proj
833 ALLOCATE (spmat(n, m))
834 spmat = 0.0_dp
835 DO i = 1, m
836 la = upf_pot%lbeta(i)
837 DO j = 1, gbasis%nbas(la)
838 spmat(j, i) = integrate_grid(upf_pot%beta(:, i), gbasis%bf(:, j, la), gbasis%grid)
839 END DO
840 END DO
841 DO i = 1, m
842 la = upf_pot%lbeta(i)
843 DO j = 1, m
844 lb = upf_pot%lbeta(j)
845 IF (la == lb) THEN
846 DO k1 = 1, gbasis%nbas(la)
847 DO k2 = 1, gbasis%nbas(la)
848 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 DEALLOCATE (spmat)
855 ELSE IF (ptype(1:2) == "SL") THEN
856 ! semi local pseudopotential
857 DO la = 0, upf_pot%l_max
858 IF (la == upf_pot%l_local) cycle
859 m = SIZE(upf_pot%vsemi(:, la + 1))
860 ALLOCATE (spot(m))
861 spot(:) = upf_pot%vsemi(:, la + 1) - upf_pot%vlocal(:)
862 spot(:) = spot(:)*cutpotu(:)
863 n = basis%nbas(la)
864 DO i = 1, n
865 DO j = i, n
866 hnl(i, j, la) = hnl(i, j, la) + &
867 integrate_grid(spot(:), &
868 gbasis%bf(:, i, la), gbasis%bf(:, j, la), gbasis%grid)
869 hnl(j, i, la) = hnl(i, j, la)
870 END DO
871 END DO
872 DEALLOCATE (spot)
873 END DO
874 ELSE
875 cpabort("Pseudopotential type: ["//adjustl(trim(ptype))//"] not known")
876 END IF
877
878 ! release basis representation on UPF grid
879 CALL release_atom_basis(gbasis)
880
881 END SUBROUTINE upfints
882
883! **************************************************************************************************
884!> \brief ...
885!> \param hnl ...
886!> \param basis ...
887!> \param ecp_pot ...
888! **************************************************************************************************
889 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: cpot
897 REAL(kind=dp), DIMENSION(:), POINTER :: rad
898
899 rad => basis%grid%rad
900 m = basis%grid%nr
901 ALLOCATE (cpot(1:m))
902
903 ! non local pseudopotential
904 hnl = 0.0_dp
905 DO l = 0, ecp_pot%lmax
906 cpot = 0._dp
907 DO k = 1, ecp_pot%npot(l)
908 n = ecp_pot%nrpot(k, l)
909 alpha = ecp_pot%bpot(k, l)
910 cpot(:) = cpot(:) + ecp_pot%apot(k, l)*rad(:)**(n - 2)*exp(-alpha*rad(:)**2)
911 END DO
912 DO i = 1, basis%nbas(l)
913 DO j = i, basis%nbas(l)
914 hnl(i, j, l) = integrate_grid(cpot, basis%bf(:, i, l), basis%bf(:, j, l), basis%grid)
915 hnl(j, i, l) = hnl(i, j, l)
916 END DO
917 END DO
918 END DO
919 DEALLOCATE (cpot)
920
921 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: cpot, pgauss
940 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: qmat
941 REAL(kind=dp), DIMENSION(:), POINTER :: rad
942
943 rad => basis%grid%rad
944 m = basis%grid%nr
945
946 ! local pseudopotential
947 ALLOCATE (cpot(m))
948 IF (sgp_pot%has_local) THEN
949 zval = sgp_pot%zval
950 core = 0._dp
951 cpot = 0.0_dp
952 DO i = 1, sgp_pot%n_local
953 cpot(:) = cpot(:) + sgp_pot%c_local(i)*exp(-sgp_pot%a_local(i)*rad(:)**2)
954 END DO
955 cpot(:) = cpot(:)*cutpots(:)
956 CALL numpot_matrix(core, cpot, basis, 0)
957 END IF
958 DEALLOCATE (cpot)
959
960 ! nonlocal pseudopotential
961 IF (sgp_pot%has_nonlocal) THEN
962 hnl = 0.0_dp
963 ALLOCATE (pgauss(1:m))
964 n = sgp_pot%n_nonlocal
965 !
966 DO l = 0, sgp_pot%lmax
967 cpassert(l <= ubound(basis%nbas, 1))
968 IF (.NOT. sgp_pot%is_nonlocal(l)) cycle
969 ! overlap (a|p)
970 na = basis%nbas(l)
971 ALLOCATE (qmat(na, n))
972 DO i = 1, n
973 a = sgp_pot%a_nonlocal(i)
974 pgauss(:) = cutpots(:)*exp(-a*rad(:)**2)*rad(:)**l
975 DO ia = 1, na
976 qmat(ia, i) = integrate_grid(basis%bf(:, ia, l), pgauss(:), basis%grid)
977 END DO
978 END DO
979 qmat(1:na, 1:n) = sqrt(fourpi)*matmul(qmat(1:na, 1:n), sgp_pot%c_nonlocal(1:n, 1:n, l))
980 DO i = 1, na
981 DO j = i, na
982 DO ia = 1, n
983 hnl(i, j, l) = hnl(i, j, l) + qmat(i, ia)*qmat(j, ia)*sgp_pot%h_nonlocal(ia, l)
984 END DO
985 hnl(j, i, l) = hnl(i, j, l)
986 END DO
987 END DO
988 DEALLOCATE (qmat)
989 END DO
990 DEALLOCATE (pgauss)
991 END IF
992
993 END SUBROUTINE sgpints
994
995! **************************************************************************************************
996!> \brief ...
997!> \param ac ...
998!> \param vlocal ...
999!> \param r ...
1000!> \param z ...
1001! **************************************************************************************************
1002 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: epot, rval, vpot
1012
1013 m = SIZE(r)
1014 ALLOCATE (epot(m), vpot(m), rval(m))
1015 cpassert(SIZE(vlocal) == m)
1016 IF (r(1) > r(m)) THEN
1017 DO i = 1, m
1018 vpot(m - i + 1) = vlocal(i)
1019 rval(m - i + 1) = r(i)
1020 END DO
1021 ELSE
1022 vpot(1:m) = vlocal(1:m)
1023 rval(1:m) = r(1:m)
1024 END IF
1025 m1 = 1
1026 DO i = 1, m
1027 IF (rval(i) > rcut) THEN
1028 m1 = i
1029 EXIT
1030 END IF
1031 END DO
1032
1033 a1 = 0.2_dp
1034 a2 = 0.2_dp
1035 e2 = 1.e20_dp
1036 epot = 0.0_dp
1037 DO i = 0, 20
1038 an = a1 + i*0.025_dp
1039 rc = 1._dp/(an*sqrt(2.0_dp))
1040 DO j = m1, m
1041 epot(j) = vpot(j) + z/rval(j)*erf(rval(j)*rc)
1042 END DO
1043 en = sum(abs(epot(m1:m)*rval(m1:m)**2))
1044 IF (en < e2) THEN
1045 e2 = en
1046 a2 = an
1047 END IF
1048 END DO
1049 ac = a2
1050
1051 DEALLOCATE (epot, vpot, rval)
1052
1053 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 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 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: tv
1078 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: rmat, sinv, smat
1079 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :, :) :: gmat
1080
1081 cl = 0.0_dp
1082 IF (method == 1) THEN
1083 ALLOCATE (tv(nl), smat(nl, nl), sinv(nl, nl))
1084 DO ir = 1, nl
1085 vgauss(:) = exp(-al(ir)*rad(:)**2)
1086 DO ix = 1, nl
1087 smat(ir, ix) = sum(vgauss(:)*exp(-al(ix)*rad(:)**2)*ww(:))
1088 END DO
1089 tv(ir) = sum(vloc(:)*vgauss(:)*ww(:))
1090 END DO
1091 CALL get_pseudo_inverse_diag(smat(1:nl, 1:nl), sinv(1:nl, 1:nl), 1.e-12_dp)
1092 cl(1:nl) = matmul(sinv(1:nl, 1:nl), tv(1:nl))
1093 ELSE
1094 !
1095 ALLOCATE (tv(nl), smat(nl, nl), sinv(nl, nl))
1096 !
1097 smat = 0.0_dp
1098 tv = 0.0_dp
1099 DO la = 0, min(ubound(gbasis%nbas, 1), 3)
1100 na = gbasis%nbas(la)
1101 ALLOCATE (rmat(na, na), gmat(na, na, nl))
1102 rmat = 0.0_dp
1103 gmat = 0.0_dp
1104 DO ia = 1, na
1105 DO ib = ia, na
1106 rmat(ia, ib) = sum(gbasis%bf(:, ia, la)*gbasis%bf(:, ib, la)*vloc(:)*ww(:))
1107 rmat(ib, ia) = rmat(ia, ib)
1108 END DO
1109 END DO
1110 DO ir = 1, nl
1111 vgauss(:) = exp(-al(ir)*rad(:)**2)
1112 DO ia = 1, na
1113 DO ib = ia, na
1114 gmat(ia, ib, ir) = sum(gbasis%bf(:, ia, la)*gbasis%bf(:, ib, la)*vgauss(:)*ww(:))
1115 gmat(ib, ia, ir) = gmat(ia, ib, ir)
1116 END DO
1117 END DO
1118 END DO
1119 DO ir = 1, nl
1120 tv(ir) = tv(ir) + accurate_dot_product(rmat, gmat(:, :, ir))
1121 DO ix = ir, nl
1122 smat(ir, ix) = smat(ir, ix) + accurate_dot_product(gmat(:, :, ix), gmat(:, :, ir))
1123 smat(ix, ir) = smat(ir, ix)
1124 END DO
1125 END DO
1126 DEALLOCATE (rmat, gmat)
1127 END DO
1128 sinv = 0.0_dp
1129 CALL get_pseudo_inverse_diag(smat(1:nl, 1:nl), sinv(1:nl, 1:nl), 1.e-12_dp)
1130 cl(1:nl) = matmul(sinv(1:nl, 1:nl), tv(1:nl))
1131 END IF
1132 !
1133 vgauss = 0.0_dp
1134 DO ir = 1, nl
1135 vgauss(:) = vgauss(:) + cl(ir)*exp(-al(ir)*rad(:)**2)
1136 END DO
1137 errloc = sum((vgauss - vloc)**2*ww)
1138 !
1139 DEALLOCATE (tv, smat, sinv)
1140 !
1141 END SUBROUTINE pplocal_error
1142
1143! **************************************************************************************************
1144!> \brief ...
1145!> \param pot ...
1146!> \param r ...
1147!> \param rcut ...
1148!> \param rsmooth ...
1149! **************************************************************************************************
1150 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 n = SIZE(pot)
1159 cpassert(n <= SIZE(r))
1160
1161 pot(:) = 1.0_dp
1162 DO i = 1, n
1163 rab = r(i)
1164 IF (rab > rcut) THEN
1165 pot(i) = 0.0_dp
1166 ELSE IF (rab > rcut - rsmooth) THEN
1167 rx = rab - (rcut - rsmooth)
1168 x = rx/rsmooth
1169 pot(i) = -6._dp*x**5 + 15._dp*x**4 - 10._dp*x**3 + 1._dp
1170 END IF
1171 END DO
1172
1173 END SUBROUTINE cutpot
1174
1175END MODULE atom_sgp
atom_utils::integrate_grid
Definition atom_utils.F:83
kahan_sum::accurate_dot_product
Definition kahan_sum.F:52
ai_onecenter
Definition ai_onecenter.F:28
ai_onecenter::sg_overlap
subroutine, public sg_overlap(smat, l, pa, pb)
...
Definition ai_onecenter.F:65
atom_set_basis
Definition atom_set_basis.F:9
atom_set_basis::set_kind_basis_atomic
subroutine, public set_kind_basis_atomic(basis, orb_basis_set, has_pp, agrid, cp2k_norm)
...
Definition atom_set_basis.F:58
atom_sgp
Definition atom_sgp.F:9
atom_sgp::atom_sgp_release
subroutine, public atom_sgp_release(sgp_pot)
...
Definition atom_sgp.F:767
atom_sgp::sgp_construction
subroutine, public sgp_construction(sgp_pot, ecp_pot, upf_pot, orb_basis, error)
...
Definition atom_sgp.F:79
atom_sgp::atom_sgp_construction
subroutine, public atom_sgp_construction(atom_info, input_section, iw)
...
Definition atom_sgp.F:201
atom_types
Define the atom type and its sub types.
Definition atom_types.F:15
atom_types::release_opmat
subroutine, public release_opmat(opmat)
...
Definition atom_types.F:1368
atom_types::release_atom_basis
subroutine, public release_atom_basis(basis)
...
Definition atom_types.F:910
atom_types::init_atom_basis_default_pp
subroutine, public init_atom_basis_default_pp(basis)
...
Definition atom_types.F:707
atom_types::create_opmat
subroutine, public create_opmat(opmat, n, lmax)
...
Definition atom_types.F:1340
atom_types::atom_basis_gridrep
subroutine, public atom_basis_gridrep(basis, gbasis, r, rab)
...
Definition atom_types.F:778
atom_upf
Routines that process Quantum Espresso UPF files.
Definition atom_upf.F:14
atom_utils
Some basic routines for atomic calculations.
Definition atom_utils.F:15
atom_utils::numpot_matrix
subroutine, public numpot_matrix(imat, cpot, basis, derivatives)
Calculate a potential matrix by integrating the potential on an atomic radial grid.
Definition atom_utils.F:1840
basis_set_types
Definition basis_set_types.F:15
input_constants
collects all constants needed in input so that they can be used without circular dependencies
Definition input_constants.F:17
input_constants::gth_pseudo
integer, parameter, public gth_pseudo
Definition input_constants.F:174
input_constants::ecp_pseudo
integer, parameter, public ecp_pseudo
Definition input_constants.F:174
input_constants::upf_pseudo
integer, parameter, public upf_pseudo
Definition input_constants.F:174
input_constants::no_pseudo
integer, parameter, public no_pseudo
Definition input_constants.F:174
input_constants::gaussian
integer, parameter, public gaussian
Definition input_constants.F:153
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
input_section_types::section_vals_get
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
Definition input_section_types.F:697
kahan_sum
sums arrays of real/complex numbers with much reduced round-off as compared to a naive implementation...
Definition kahan_sum.F:29
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::dfac
real(kind=dp), dimension(-1:2 *maxfac+1), parameter, public dfac
Definition mathconstants.F:61
mathconstants::rootpi
real(kind=dp), parameter, public rootpi
Definition mathconstants.F:114
mathconstants::fourpi
real(kind=dp), parameter, public fourpi
Definition mathconstants.F:113
mathconstants::sqrt2
real(kind=dp), parameter, public sqrt2
Definition mathconstants.F:127
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
mathlib::get_pseudo_inverse_diag
subroutine, public get_pseudo_inverse_diag(a, a_pinverse, rskip)
returns the pseudoinverse of a real, symmetric and positive definite matrix using diagonalization.
Definition mathlib.F:1030
mathlib::diamat_all
subroutine, public diamat_all(a, eigval, dac)
Diagonalize the symmetric n by n matrix a using the LAPACK library. Only the upper triangle of matrix...
Definition mathlib.F:372
powell
Definition powell.F:9
powell::powell_optimize
subroutine, public powell_optimize(n, x, optstate)
...
Definition powell.F:55
atom_sgp::atom_sgp_potential_type
Definition atom_sgp.F:40
atom_types::atom_basis_type
Provides all information about a basis set.
Definition atom_types.F:78
atom_types::atom_ecppot_type
Definition atom_types.F:130
atom_types::atom_p_type
Definition atom_types.F:320
atom_types::atom_type
Provides all information about an atomic kind.
Definition atom_types.F:290
atom_types::opmat_type
Operator matrices.
Definition atom_types.F:248
atom_upf::atom_upfpot_type
Definition atom_upf.F:32
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:71
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:126
powell::opt_state_type
Definition powell.F:18