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 : ! **************************************************************************************************
9 : !> \par History
10 : !> none
11 : !> \author JHU (9.2022)
12 : ! **************************************************************************************************
13 : MODULE gapw_1c_basis_set
14 :
15 : USE basis_set_types, ONLY: allocate_gto_basis_set,&
16 : combine_basis_sets,&
17 : copy_gto_basis_set,&
18 : create_primitive_basis_set,&
19 : deallocate_gto_basis_set,&
20 : get_gto_basis_set,&
21 : gto_basis_set_type
22 : USE kinds, ONLY: dp
23 : USE orbital_transformation_matrices, ONLY: init_spherical_harmonics
24 : #include "base/base_uses.f90"
25 :
26 : IMPLICIT NONE
27 :
28 : PRIVATE
29 :
30 : ! *** Global parameters (only in this module)
31 :
32 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'gapw_1c_basis_set'
33 :
34 : INTEGER, PARAMETER :: max_name_length = 60
35 :
36 : ! *** Public subroutines ***
37 :
38 : PUBLIC :: create_1c_basis
39 :
40 : CONTAINS
41 :
42 : ! **************************************************************************************************
43 : !> \brief create the one center basis from the orbital basis
44 : !> \param orb_basis ...
45 : !> \param soft_basis ...
46 : !> \param gapw_1c_basis ...
47 : !> \param basis_1c_level ...
48 : !> \version 1.0
49 : ! **************************************************************************************************
50 3034 : SUBROUTINE create_1c_basis(orb_basis, soft_basis, gapw_1c_basis, basis_1c_level)
51 :
52 : TYPE(gto_basis_set_type), POINTER :: orb_basis, soft_basis, gapw_1c_basis
53 : INTEGER, INTENT(IN) :: basis_1c_level
54 :
55 : INTEGER :: i, ipgf, iset, j, l, lbas, maxl, maxlo, &
56 : maxls, mp, n1, n2, nn, nseto, nsets
57 3034 : INTEGER, ALLOCATABLE, DIMENSION(:) :: nps
58 3034 : INTEGER, DIMENSION(:), POINTER :: lmaxo, lmaxs, lmino, lmins, npgfo, npgfs
59 : REAL(KIND=dp) :: fmin, fr1, fr2, fz, rr, xz, yz, zmall, &
60 : zms, zz
61 3034 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: z1, z2, zmaxs
62 3034 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: zeta, zexs
63 3034 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: zeto, zets
64 : TYPE(gto_basis_set_type), POINTER :: ext_basis, p_basis
65 :
66 0 : CPASSERT(.NOT. ASSOCIATED(gapw_1c_basis))
67 :
68 3034 : IF (basis_1c_level == 0) THEN
69 : ! we use the orbital basis set
70 2814 : CALL copy_gto_basis_set(orb_basis, gapw_1c_basis)
71 : ELSE
72 220 : CALL copy_gto_basis_set(orb_basis, gapw_1c_basis)
73 220 : NULLIFY (ext_basis)
74 220 : CALL allocate_gto_basis_set(ext_basis)
75 : ! get information on orbital basis and soft basis
76 : CALL get_gto_basis_set(gto_basis_set=orb_basis, maxl=maxlo, nset=nseto, &
77 220 : lmin=lmino, lmax=lmaxo, npgf=npgfo, zet=zeto)
78 : CALL get_gto_basis_set(gto_basis_set=soft_basis, maxl=maxls, nset=nsets, &
79 220 : lmin=lmins, lmax=lmaxs, npgf=npgfs, zet=zets)
80 : ! determine max soft exponent per l-qn
81 220 : maxl = MAX(maxls, maxlo)
82 1100 : ALLOCATE (zmaxs(0:maxl), nps(0:maxl))
83 220 : zmaxs = 0.0_dp
84 586 : DO iset = 1, nsets
85 1570 : zms = MAXVAL(zets(:, iset))
86 1116 : DO l = lmins(iset), lmaxs(iset)
87 896 : zmaxs(l) = MAX(zmaxs(l), zms)
88 : END DO
89 : END DO
90 732 : zmall = MAXVAL(zmaxs)
91 : ! in case of missing soft basis!
92 220 : zmall = MAX(zmall, 0.20_dp)
93 : ! create pools of exponents for each l-qn
94 220 : nps = 0
95 586 : DO iset = 1, nsets
96 1116 : DO l = lmins(iset), lmaxs(iset)
97 896 : nps(l) = nps(l) + npgfs(iset)
98 : END DO
99 : END DO
100 732 : mp = MAXVAL(nps)
101 880 : ALLOCATE (zexs(1:mp, 0:maxl))
102 220 : zexs = 0.0_dp
103 220 : nps = 0
104 586 : DO iset = 1, nsets
105 1528 : DO ipgf = 1, npgfs(iset)
106 2832 : DO l = lmins(iset), lmaxs(iset)
107 1524 : nps(l) = nps(l) + 1
108 2466 : zexs(nps(l), l) = zets(ipgf, iset)
109 : END DO
110 : END DO
111 : END DO
112 :
113 198 : SELECT CASE (basis_1c_level)
114 : CASE (1)
115 198 : lbas = maxl
116 198 : fr1 = 2.50_dp
117 198 : fr2 = 2.50_dp
118 : CASE (2)
119 6 : lbas = maxl
120 6 : fr1 = 2.00_dp
121 6 : fr2 = 2.50_dp
122 : CASE (3)
123 10 : lbas = maxl + 1
124 10 : fr1 = 1.75_dp
125 10 : fr2 = 2.50_dp
126 : CASE (4)
127 6 : lbas = maxl + 2
128 6 : fr1 = 1.50_dp
129 6 : fr2 = 2.50_dp
130 : CASE DEFAULT
131 220 : CPABORT("unknown case")
132 : END SELECT
133 220 : lbas = MIN(lbas, 5)
134 : !
135 220 : CALL init_spherical_harmonics(lbas, 0)
136 : !
137 220 : rr = LOG(zmall/0.05_dp)/LOG(fr1)
138 220 : n1 = INT(rr) + 1
139 220 : rr = LOG(zmall/0.05_dp)/LOG(fr2)
140 220 : n2 = INT(rr) + 1
141 1100 : ALLOCATE (z1(n1), z2(n2))
142 220 : z1 = 0.0_dp
143 220 : zz = zmall*SQRT(fr1)
144 1302 : DO i = 1, n1
145 1302 : z1(i) = zz/(fr1**(i - 1))
146 : END DO
147 220 : z2 = 0.0_dp
148 220 : zz = zmall
149 1238 : DO i = 1, n2
150 1238 : z2(i) = zz/(fr2**(i - 1))
151 : END DO
152 880 : ALLOCATE (zeta(MAX(n1, n2), lbas + 1))
153 220 : zeta = 0.0_dp
154 : !
155 220 : ext_basis%nset = lbas + 1
156 1100 : ALLOCATE (ext_basis%lmin(lbas + 1), ext_basis%lmax(lbas + 1))
157 440 : ALLOCATE (ext_basis%npgf(lbas + 1))
158 754 : DO l = 0, lbas
159 534 : ext_basis%lmin(l + 1) = l
160 534 : ext_basis%lmax(l + 1) = l
161 754 : IF (l <= maxl) THEN
162 : nn = 0
163 3052 : DO i = 1, n1
164 2540 : xz = z1(i)
165 : ! Each candidate is compared independently with the soft basis.
166 2540 : fmin = 10.0_dp
167 10242 : DO j = 1, nps(l)
168 7702 : yz = zexs(j, l)
169 7702 : fz = MAX(xz, yz)/MIN(xz, yz)
170 10242 : fmin = MIN(fmin, fz)
171 : END DO
172 3052 : IF (fmin > fr1**0.25_dp) THEN
173 1878 : nn = nn + 1
174 1878 : zeta(nn, l + 1) = xz
175 : END IF
176 : END DO
177 512 : CPASSERT(nn > 0)
178 512 : ext_basis%npgf(l + 1) = nn
179 : ELSE
180 22 : ext_basis%npgf(l + 1) = n2
181 124 : zeta(1:n2, l + 1) = z2(1:n2)
182 : END IF
183 : END DO
184 754 : nn = MAXVAL(ext_basis%npgf)
185 880 : ALLOCATE (ext_basis%zet(nn, lbas + 1))
186 754 : DO i = 1, lbas + 1
187 534 : nn = ext_basis%npgf(i)
188 2734 : ext_basis%zet(1:nn, i) = zeta(1:nn, i)
189 : END DO
190 220 : ext_basis%name = "extbas"
191 220 : ext_basis%kind_radius = orb_basis%kind_radius
192 220 : ext_basis%short_kind_radius = orb_basis%short_kind_radius
193 220 : ext_basis%norm_type = orb_basis%norm_type
194 :
195 220 : NULLIFY (p_basis)
196 220 : CALL create_primitive_basis_set(ext_basis, p_basis)
197 220 : CALL combine_basis_sets(gapw_1c_basis, p_basis)
198 :
199 220 : CALL deallocate_gto_basis_set(ext_basis)
200 220 : CALL deallocate_gto_basis_set(p_basis)
201 440 : DEALLOCATE (zmaxs, zexs, nps, z1, z2, zeta)
202 : END IF
203 :
204 3034 : END SUBROUTINE create_1c_basis
205 :
206 : END MODULE gapw_1c_basis_set
|
