(git:58e3e09)
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 ! **************************************************************************************************
9 MODULE atom_sgp
10  USE ai_onecenter, ONLY: sg_overlap
12  USE atom_types, ONLY: &
13  atom_basis_gridrep, atom_basis_type, atom_ecppot_type, atom_p_type, atom_type, &
15  USE atom_upf, ONLY: atom_upfpot_type
16  USE atom_utils, ONLY: integrate_grid,&
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,&
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,&
34  USE powell, ONLY: opt_state_type,&
36 #include "./base/base_uses.f90"
37 
38  IMPLICIT NONE
39 
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
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  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 
1175 END MODULE atom_sgp
subroutine, public sg_overlap(smat, l, pa, pb)
...
Definition: ai_onecenter.F:64
subroutine, public set_kind_basis_atomic(basis, orb_basis_set, has_pp, agrid, cp2k_norm)
...
subroutine, public atom_sgp_release(sgp_pot)
...
Definition: atom_sgp.F:767
subroutine, public sgp_construction(sgp_pot, ecp_pot, upf_pot, orb_basis, error)
...
Definition: atom_sgp.F:79
subroutine, public atom_sgp_construction(atom_info, input_section, iw)
...
Definition: atom_sgp.F:201
Define the atom type and its sub types.
Definition: atom_types.F:15
subroutine, public release_opmat(opmat)
...
Definition: atom_types.F:1368
subroutine, public release_atom_basis(basis)
...
Definition: atom_types.F:910
subroutine, public init_atom_basis_default_pp(basis)
...
Definition: atom_types.F:707
subroutine, public create_opmat(opmat, n, lmax)
...
Definition: atom_types.F:1340
subroutine, public atom_basis_gridrep(basis, gbasis, r, rab)
...
Definition: atom_types.F:778
Routines that process Quantum Espresso UPF files.
Definition: atom_upf.F:14
Some basic routines for atomic calculations.
Definition: atom_utils.F:15
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
collects all constants needed in input so that they can be used without circular dependencies
integer, parameter, public gth_pseudo
integer, parameter, public ecp_pseudo
integer, parameter, public upf_pseudo
integer, parameter, public no_pseudo
integer, parameter, public gaussian
objects that represent the structure of input sections and the data contained in an input section
subroutine, public section_vals_get(section_vals, ref_count, n_repetition, n_subs_vals_rep, section, explicit)
returns various attributes about the section_vals
sums arrays of real/complex numbers with much reduced round-off as compared to a naive implementation...
Definition: kahan_sum.F:29
Defines the basic variable types.
Definition: kinds.F:23
integer, parameter, public dp
Definition: kinds.F:34
Definition of mathematical constants and functions.
Definition: mathconstants.F:16
real(kind=dp), dimension(-1:2 *maxfac+1), parameter, public dfac
Definition: mathconstants.F:61
real(kind=dp), parameter, public rootpi
real(kind=dp), parameter, public fourpi
real(kind=dp), parameter, public sqrt2
Collection of simple mathematical functions and subroutines.
Definition: mathlib.F:15
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:1036
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:376
Definition: powell.F:9
subroutine, public powell_optimize(n, x, optstate)
...
Definition: powell.F:55