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