Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 CP2K developers group <https://cp2k.org> !
4 : ! !
5 : ! SPDX-License-Identifier: GPL-2.0-or-later !
6 : !--------------------------------------------------------------------------------------------------!
7 :
8 2 : PROGRAM gapw_unittest
9 : USE basis_set_types, ONLY: allocate_gto_basis_set,&
10 2 : create_primitive_basis_set,&
11 : deallocate_gto_basis_set,&
12 : gto_basis_set_type
13 : USE cp_log_handling, ONLY: cp_logger_get_default_io_unit
14 : USE gapw_1c_basis_set, ONLY: create_1c_basis
15 : USE hartree_local_methods, ONLY: gapw_background_potential
16 : USE kinds, ONLY: dp
17 : USE libcp2k, ONLY: cp2k_finalize,&
18 : cp2k_init
19 : USE mathconstants, ONLY: fourpi,&
20 : pi
21 : USE orbital_pointers, ONLY: init_orbital_pointers
22 : USE orbital_transformation_matrices, ONLY: init_spherical_harmonics
23 :
24 : IMPLICIT NONE
25 :
26 : INTEGER, PARAMETER :: nr = 1601, ng = 32
27 : REAL(dp), PARAMETER :: alpha = 1.4_dp, beta = 0.55_dp, box = 12.0_dp, dr = 0.00625_dp
28 : REAL(dp), PARAMETER :: charges(2) = [-1.0_dp, 0.3_dp]
29 : REAL(dp) :: rad2(nr), weights(nr), delta_rho(nr), potential(nr), shifted(nr)
30 : REAL(dp) :: q, g2, periodic_difference, local_difference, correction, error, h, eplus, eminus
31 : INTEGER :: i, ix, iy, iz, io_unit
32 :
33 2 : CALL cp2k_init()
34 2 : io_unit = cp_logger_get_default_io_unit()
35 2 : CALL check_1c_basis()
36 3204 : DO i = 1, nr
37 3204 : rad2(i) = ((i - 1)*dr)**2
38 : END DO
39 3204 : weights = 2.0_dp*dr/3.0_dp
40 1602 : weights(2:nr - 1:2) = 4.0_dp*dr/3.0_dp
41 2 : weights(1) = dr/3.0_dp
42 2 : weights(nr) = dr/3.0_dp
43 3204 : weights = weights*rad2
44 :
45 6 : DO i = 1, SIZE(charges)
46 4 : q = charges(i)
47 : ! Electron-minus-nuclear l=0 density. Both normalized Gaussians carry q,
48 : ! so their difference has zero monopole and a nonzero second moment.
49 : delta_rho = q*SQRT(fourpi)*((alpha/pi)**1.5_dp*EXP(-alpha*rad2) - &
50 6408 : (beta/pi)**1.5_dp*EXP(-beta*rad2))
51 4 : potential = gapw_background_potential(-q, box**3, rad2)
52 6408 : correction = SUM(potential*delta_rho*weights)
53 :
54 : ! Independent reference: exact reciprocal-space Coulomb energies of the
55 : ! two periodic Gaussians with G=0 omitted, minus their free-space energies.
56 : periodic_difference = 0.0_dp
57 264 : DO iz = -ng, ng
58 17164 : DO iy = -ng, ng
59 1115660 : DO ix = -ng, ng
60 1098500 : IF (ix == 0 .AND. iy == 0 .AND. iz == 0) CYCLE
61 1098496 : g2 = (2.0_dp*pi/box)**2*REAL(ix*ix + iy*iy + iz*iz, dp)
62 : periodic_difference = periodic_difference + &
63 1115400 : 2.0_dp*pi*q*q/box**3/g2*(EXP(-g2/(2.0_dp*alpha)) - EXP(-g2/(2.0_dp*beta)))
64 : END DO
65 : END DO
66 : END DO
67 4 : local_difference = q*q*(SQRT(alpha) - SQRT(beta))/SQRT(2.0_dp*pi)
68 4 : error = ABS(periodic_difference - local_difference - correction)
69 4 : IF (io_unit > 0) WRITE (io_unit, '(A,F8.3,ES16.7)') "GAPW periodic Gaussian background:", q, error
70 4 : IF (error > 1.e-12_dp) ERROR STOP "GAPW background disagrees with periodic Gaussian electrostatics"
71 :
72 : ! The explicit volume derivative contributes E_background to each
73 : ! diagonal component of -dE/dstrain (local radial functions held fixed).
74 4 : h = 1.e-5_dp
75 4 : shifted = gapw_background_potential(-q, box**3*(1.0_dp + h), rad2)
76 6408 : eplus = SUM(shifted*delta_rho*weights)
77 4 : shifted = gapw_background_potential(-q, box**3*(1.0_dp - h), rad2)
78 6408 : eminus = SUM(shifted*delta_rho*weights)
79 6 : IF (ABS(-(eplus - eminus)/(2.0_dp*h) - correction) > 1.e-11_dp) THEN
80 0 : ERROR STOP "GAPW background volume derivative failed"
81 : END IF
82 : END DO
83 2 : potential = gapw_background_potential(0.0_dp, box**3, rad2)
84 3204 : IF (ANY(potential /= 0.0_dp)) ERROR STOP "Neutral cell has a background potential"
85 2 : CALL cp2k_finalize()
86 : CONTAINS
87 :
88 : ! **************************************************************************************************
89 : !> \brief Check that one-centre basis extension retains diffuse exponents without duplicates.
90 : ! **************************************************************************************************
91 2 : SUBROUTINE check_1c_basis()
92 : INTEGER, PARAMETER :: expected(4) = [3, 4, 5, 7]
93 : REAL(KIND=dp), PARAMETER :: spacing(4) = [2.5_dp, 2.0_dp, 1.75_dp, 1.5_dp]
94 :
95 : INTEGER :: level
96 : REAL(KIND=dp) :: diffuse, ratio
97 : TYPE(gto_basis_set_type), POINTER :: extended, orbitals, seed
98 :
99 2 : NULLIFY (extended, orbitals, seed)
100 2 : CALL init_orbital_pointers(5)
101 2 : CALL init_spherical_harmonics(5, 0)
102 10 : DO level = 1, 4
103 8 : ratio = SPACING(level)
104 8 : CALL allocate_gto_basis_set(seed)
105 8 : seed%name = 'two-s-primitives'
106 8 : seed%norm_type = 2
107 8 : seed%nset = 1
108 8 : ALLOCATE (seed%lmin(1), seed%lmax(1), seed%npgf(1), seed%zet(2, 1))
109 16 : seed%lmin = 0
110 16 : seed%lmax = 0
111 16 : seed%npgf = 2
112 : ! The second extension candidate duplicates the second soft exponent.
113 : ! More diffuse, well-separated candidates must still be admitted.
114 24 : seed%zet(:, 1) = [1.0_dp, 1.0_dp/SQRT(ratio)]
115 8 : CALL create_primitive_basis_set(seed, orbitals)
116 16 : orbitals%set_radius = 0.0_dp
117 8 : CALL create_1c_basis(orbitals, orbitals, extended, level)
118 8 : diffuse = 1.0_dp/(ratio*SQRT(ratio))
119 8 : IF (io_unit > 0) WRITE (io_unit, '(A,I2,A,I3)') &
120 4 : 'GAPW extension level ', level, ': added s primitives ', extended%npgf(2)
121 8 : IF (extended%npgf(2) /= expected(level)) ERROR STOP 'GAPW extension lost admissible diffuse exponents'
122 46 : IF (MINVAL(ABS(extended%zet(:extended%npgf(2), 2) - diffuse)) > 1.e-12_dp) THEN
123 0 : ERROR STOP 'GAPW extension is missing the first admissible diffuse exponent'
124 : END IF
125 46 : IF (MINVAL(ABS(extended%zet(:extended%npgf(2), 2) - seed%zet(2, 1))) < 1.e-12_dp) THEN
126 0 : ERROR STOP 'GAPW extension duplicated an existing soft exponent'
127 : END IF
128 8 : CALL deallocate_gto_basis_set(extended)
129 8 : CALL deallocate_gto_basis_set(orbitals)
130 10 : CALL deallocate_gto_basis_set(seed)
131 : END DO
132 2 : END SUBROUTINE check_1c_basis
133 : END PROGRAM gapw_unittest
|
