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