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 2476 : 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 2476 : INTEGER, ALLOCATABLE, DIMENSION(:) :: nps
58 2476 : 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 2476 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: z1, z2, zmaxs
62 2476 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: zeta, zexs
63 2476 : 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 2476 : IF (basis_1c_level == 0) THEN
69 : ! we use the orbital basis set
70 2360 : CALL copy_gto_basis_set(orb_basis, gapw_1c_basis)
71 : ELSE
72 116 : CALL copy_gto_basis_set(orb_basis, gapw_1c_basis)
73 116 : NULLIFY (ext_basis)
74 116 : 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 116 : lmin=lmino, lmax=lmaxo, npgf=npgfo, zet=zeto)
78 : CALL get_gto_basis_set(gto_basis_set=soft_basis, maxl=maxls, nset=nsets, &
79 116 : lmin=lmins, lmax=lmaxs, npgf=npgfs, zet=zets)
80 : ! determine max soft exponent per l-qn
81 116 : maxl = MAX(maxls, maxlo)
82 580 : ALLOCATE (zmaxs(0:maxl), nps(0:maxl))
83 404 : zmaxs = 0.0_dp
84 322 : DO iset = 1, nsets
85 1072 : zms = MAXVAL(zets(:, iset))
86 628 : DO l = lmins(iset), lmaxs(iset)
87 512 : zmaxs(l) = MAX(zmaxs(l), zms)
88 : END DO
89 : END DO
90 520 : zmall = MAXVAL(zmaxs)
91 : ! in case of missing soft basis!
92 116 : zmall = MAX(zmall, 0.20_dp)
93 : ! create pools of exponents for each l-qn
94 404 : nps = 0
95 322 : DO iset = 1, nsets
96 628 : DO l = lmins(iset), lmaxs(iset)
97 512 : nps(l) = nps(l) + npgfs(iset)
98 : END DO
99 : END DO
100 404 : mp = MAXVAL(nps)
101 464 : ALLOCATE (zexs(1:mp, 0:maxl))
102 1384 : zexs = 0.0_dp
103 404 : nps = 0
104 322 : DO iset = 1, nsets
105 820 : DO ipgf = 1, npgfs(iset)
106 1526 : DO l = lmins(iset), lmaxs(iset)
107 822 : nps(l) = nps(l) + 1
108 1320 : zexs(nps(l), l) = zets(ipgf, iset)
109 : END DO
110 : END DO
111 : END DO
112 :
113 100 : SELECT CASE (basis_1c_level)
114 : CASE (1)
115 100 : lbas = maxl
116 100 : fr1 = 2.50_dp
117 100 : 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 116 : CPABORT("unknown case")
132 : END SELECT
133 116 : lbas = MIN(lbas, 5)
134 : !
135 116 : CALL init_spherical_harmonics(lbas, 0)
136 : !
137 116 : rr = LOG(zmall/0.05_dp)/LOG(fr1)
138 116 : n1 = INT(rr) + 1
139 116 : rr = LOG(zmall/0.05_dp)/LOG(fr2)
140 116 : n2 = INT(rr) + 1
141 580 : ALLOCATE (z1(n1), z2(n2))
142 698 : z1 = 0.0_dp
143 116 : zz = zmall*SQRT(fr1)
144 698 : DO i = 1, n1
145 698 : z1(i) = zz/(fr1**(i - 1))
146 : END DO
147 648 : z2 = 0.0_dp
148 648 : zz = zmall
149 648 : DO i = 1, n2
150 648 : z2(i) = zz/(fr2**(i - 1))
151 : END DO
152 464 : ALLOCATE (zeta(MAX(n1, n2), lbas + 1))
153 2034 : zeta = 0.0_dp
154 : !
155 116 : ext_basis%nset = lbas + 1
156 464 : ALLOCATE (ext_basis%lmin(lbas + 1), ext_basis%lmax(lbas + 1))
157 232 : ALLOCATE (ext_basis%npgf(lbas + 1))
158 420 : DO l = 0, lbas
159 304 : ext_basis%lmin(l + 1) = l
160 304 : ext_basis%lmax(l + 1) = l
161 420 : IF (l <= maxl) THEN
162 : fmin = 10.0_dp
163 : nn = 0
164 1754 : DO i = 1, n1
165 1466 : xz = z1(i)
166 5742 : DO j = 1, nps(l)
167 4276 : yz = zexs(j, l)
168 4276 : fz = MAX(xz, yz)/MIN(xz, yz)
169 5742 : fmin = MIN(fmin, fz)
170 : END DO
171 1754 : IF (fmin > fr1**0.25) THEN
172 796 : nn = nn + 1
173 796 : zeta(nn, l + 1) = xz
174 : END IF
175 : END DO
176 288 : CPASSERT(nn > 0)
177 288 : 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 420 : nn = MAXVAL(ext_basis%npgf)
184 464 : ALLOCATE (ext_basis%zet(nn, lbas + 1))
185 420 : DO i = 1, lbas + 1
186 304 : nn = ext_basis%npgf(i)
187 1294 : ext_basis%zet(1:nn, i) = zeta(1:nn, i)
188 : END DO
189 116 : ext_basis%name = "extbas"
190 116 : ext_basis%kind_radius = orb_basis%kind_radius
191 116 : ext_basis%short_kind_radius = orb_basis%short_kind_radius
192 116 : ext_basis%norm_type = orb_basis%norm_type
193 :
194 116 : NULLIFY (p_basis)
195 116 : CALL create_primitive_basis_set(ext_basis, p_basis)
196 116 : CALL combine_basis_sets(gapw_1c_basis, p_basis)
197 :
198 116 : CALL deallocate_gto_basis_set(ext_basis)
199 116 : CALL deallocate_gto_basis_set(p_basis)
200 232 : DEALLOCATE (zmaxs, zexs, nps, z1, z2, zeta)
201 : END IF
202 :
203 2476 : END SUBROUTINE create_1c_basis
204 :
205 : END MODULE gapw_1c_basis_set
|
