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 2914 : 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 2914 : INTEGER, ALLOCATABLE, DIMENSION(:) :: nps
58 2914 : 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 2914 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: z1, z2, zmaxs
62 2914 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: zeta, zexs
63 2914 : 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 2914 : IF (basis_1c_level == 0) THEN
69 : ! we use the orbital basis set
70 2702 : CALL copy_gto_basis_set(orb_basis, gapw_1c_basis)
71 : ELSE
72 212 : CALL copy_gto_basis_set(orb_basis, gapw_1c_basis)
73 212 : NULLIFY (ext_basis)
74 212 : 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 212 : lmin=lmino, lmax=lmaxo, npgf=npgfo, zet=zeto)
78 : CALL get_gto_basis_set(gto_basis_set=soft_basis, maxl=maxls, nset=nsets, &
79 212 : lmin=lmins, lmax=lmaxs, npgf=npgfs, zet=zets)
80 : ! determine max soft exponent per l-qn
81 212 : maxl = MAX(maxls, maxlo)
82 1060 : ALLOCATE (zmaxs(0:maxl), nps(0:maxl))
83 212 : zmaxs = 0.0_dp
84 570 : DO iset = 1, nsets
85 1546 : zms = MAXVAL(zets(:, iset))
86 1092 : DO l = lmins(iset), lmaxs(iset)
87 880 : zmaxs(l) = MAX(zmaxs(l), zms)
88 : END DO
89 : END DO
90 716 : zmall = MAXVAL(zmaxs)
91 : ! in case of missing soft basis!
92 212 : zmall = MAX(zmall, 0.20_dp)
93 : ! create pools of exponents for each l-qn
94 212 : nps = 0
95 570 : DO iset = 1, nsets
96 1092 : DO l = lmins(iset), lmaxs(iset)
97 880 : nps(l) = nps(l) + npgfs(iset)
98 : END DO
99 : END DO
100 716 : mp = MAXVAL(nps)
101 848 : ALLOCATE (zexs(1:mp, 0:maxl))
102 212 : zexs = 0.0_dp
103 212 : nps = 0
104 570 : DO iset = 1, nsets
105 1496 : DO ipgf = 1, npgfs(iset)
106 2792 : DO l = lmins(iset), lmaxs(iset)
107 1508 : nps(l) = nps(l) + 1
108 2434 : zexs(nps(l), l) = zets(ipgf, iset)
109 : END DO
110 : END DO
111 : END DO
112 :
113 196 : SELECT CASE (basis_1c_level)
114 : CASE (1)
115 196 : lbas = maxl
116 196 : fr1 = 2.50_dp
117 196 : fr2 = 2.50_dp
118 : CASE (2)
119 4 : lbas = maxl
120 4 : fr1 = 2.00_dp
121 4 : fr2 = 2.50_dp
122 : CASE (3)
123 8 : lbas = maxl + 1
124 8 : fr1 = 1.75_dp
125 8 : fr2 = 2.50_dp
126 : CASE (4)
127 4 : lbas = maxl + 2
128 4 : fr1 = 1.50_dp
129 4 : fr2 = 2.50_dp
130 : CASE DEFAULT
131 212 : CPABORT("unknown case")
132 : END SELECT
133 212 : lbas = MIN(lbas, 5)
134 : !
135 212 : CALL init_spherical_harmonics(lbas, 0)
136 : !
137 212 : rr = LOG(zmall/0.05_dp)/LOG(fr1)
138 212 : n1 = INT(rr) + 1
139 212 : rr = LOG(zmall/0.05_dp)/LOG(fr2)
140 212 : n2 = INT(rr) + 1
141 1060 : ALLOCATE (z1(n1), z2(n2))
142 212 : z1 = 0.0_dp
143 212 : zz = zmall*SQRT(fr1)
144 1248 : DO i = 1, n1
145 1248 : z1(i) = zz/(fr1**(i - 1))
146 : END DO
147 212 : z2 = 0.0_dp
148 212 : zz = zmall
149 1198 : DO i = 1, n2
150 1198 : z2(i) = zz/(fr2**(i - 1))
151 : END DO
152 848 : ALLOCATE (zeta(MAX(n1, n2), lbas + 1))
153 212 : zeta = 0.0_dp
154 : !
155 212 : ext_basis%nset = lbas + 1
156 1060 : ALLOCATE (ext_basis%lmin(lbas + 1), ext_basis%lmax(lbas + 1))
157 424 : ALLOCATE (ext_basis%npgf(lbas + 1))
158 732 : DO l = 0, lbas
159 520 : ext_basis%lmin(l + 1) = l
160 520 : ext_basis%lmax(l + 1) = l
161 732 : IF (l <= maxl) THEN
162 : fmin = 10.0_dp
163 : nn = 0
164 2998 : DO i = 1, n1
165 2494 : xz = z1(i)
166 10104 : DO j = 1, nps(l)
167 7610 : yz = zexs(j, l)
168 7610 : fz = MAX(xz, yz)/MIN(xz, yz)
169 10104 : fmin = MIN(fmin, fz)
170 : END DO
171 2998 : IF (fmin > fr1**0.25) THEN
172 1458 : nn = nn + 1
173 1458 : zeta(nn, l + 1) = xz
174 : END IF
175 : END DO
176 504 : CPASSERT(nn > 0)
177 504 : ext_basis%npgf(l + 1) = nn
178 : ELSE
179 16 : ext_basis%npgf(l + 1) = n2
180 94 : zeta(1:n2, l + 1) = z2(1:n2)
181 : END IF
182 : END DO
183 732 : nn = MAXVAL(ext_basis%npgf)
184 848 : ALLOCATE (ext_basis%zet(nn, lbas + 1))
185 732 : DO i = 1, lbas + 1
186 520 : nn = ext_basis%npgf(i)
187 2268 : ext_basis%zet(1:nn, i) = zeta(1:nn, i)
188 : END DO
189 212 : ext_basis%name = "extbas"
190 212 : ext_basis%kind_radius = orb_basis%kind_radius
191 212 : ext_basis%short_kind_radius = orb_basis%short_kind_radius
192 212 : ext_basis%norm_type = orb_basis%norm_type
193 :
194 212 : NULLIFY (p_basis)
195 212 : CALL create_primitive_basis_set(ext_basis, p_basis)
196 212 : CALL combine_basis_sets(gapw_1c_basis, p_basis)
197 :
198 212 : CALL deallocate_gto_basis_set(ext_basis)
199 212 : CALL deallocate_gto_basis_set(p_basis)
200 424 : DEALLOCATE (zmaxs, zexs, nps, z1, z2, zeta)
201 : END IF
202 :
203 2914 : END SUBROUTINE create_1c_basis
204 :
205 : END MODULE gapw_1c_basis_set
|
