Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2024 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 1932 : 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 1932 : INTEGER, ALLOCATABLE, DIMENSION(:) :: nps
58 1932 : 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 1932 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: z1, z2, zmaxs
62 1932 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: zeta, zexs
63 1932 : 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 1932 : IF (basis_1c_level == 0) THEN
69 : ! we use the orbital basis set
70 1890 : CALL copy_gto_basis_set(orb_basis, gapw_1c_basis)
71 : ELSE
72 42 : CALL copy_gto_basis_set(orb_basis, gapw_1c_basis)
73 42 : NULLIFY (ext_basis)
74 42 : 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 42 : lmin=lmino, lmax=lmaxo, npgf=npgfo, zet=zeto)
78 : CALL get_gto_basis_set(gto_basis_set=soft_basis, maxl=maxls, nset=nsets, &
79 42 : lmin=lmins, lmax=lmaxs, npgf=npgfs, zet=zets)
80 : ! determine max soft exponent per l-qn
81 42 : maxl = MAX(maxls, maxlo)
82 210 : ALLOCATE (zmaxs(0:maxl), nps(0:maxl))
83 144 : zmaxs = 0.0_dp
84 120 : DO iset = 1, nsets
85 450 : zms = MAXVAL(zets(:, iset))
86 222 : DO l = lmins(iset), lmaxs(iset)
87 180 : zmaxs(l) = MAX(zmaxs(l), zms)
88 : END DO
89 : END DO
90 186 : zmall = MAXVAL(zmaxs)
91 : ! in case of missing soft basis!
92 42 : zmall = MAX(zmall, 0.20_dp)
93 : ! create pools of exponents for each l-qn
94 144 : nps = 0
95 120 : DO iset = 1, nsets
96 222 : DO l = lmins(iset), lmaxs(iset)
97 180 : nps(l) = nps(l) + npgfs(iset)
98 : END DO
99 : END DO
100 144 : mp = MAXVAL(nps)
101 168 : ALLOCATE (zexs(1:mp, 0:maxl))
102 532 : zexs = 0.0_dp
103 144 : nps = 0
104 120 : DO iset = 1, nsets
105 328 : DO ipgf = 1, npgfs(iset)
106 588 : DO l = lmins(iset), lmaxs(iset)
107 302 : nps(l) = nps(l) + 1
108 510 : zexs(nps(l), l) = zets(ipgf, iset)
109 : END DO
110 : END DO
111 : END DO
112 :
113 26 : SELECT CASE (basis_1c_level)
114 : CASE (1)
115 26 : lbas = maxl
116 26 : fr1 = 2.50_dp
117 26 : 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 42 : CPABORT("unknown case")
132 : END SELECT
133 42 : lbas = MIN(lbas, 5)
134 : !
135 42 : CALL init_spherical_harmonics(lbas, 0)
136 : !
137 42 : rr = LOG(zmall/0.05_dp)/LOG(fr1)
138 42 : n1 = INT(rr) + 1
139 42 : rr = LOG(zmall/0.05_dp)/LOG(fr2)
140 42 : n2 = INT(rr) + 1
141 210 : ALLOCATE (z1(n1), z2(n2))
142 286 : z1 = 0.0_dp
143 42 : zz = zmall*SQRT(fr1)
144 286 : DO i = 1, n1
145 286 : z1(i) = zz/(fr1**(i - 1))
146 : END DO
147 236 : z2 = 0.0_dp
148 236 : zz = zmall
149 236 : DO i = 1, n2
150 236 : z2(i) = zz/(fr2**(i - 1))
151 : END DO
152 168 : ALLOCATE (zeta(MAX(n1, n2), lbas + 1))
153 910 : zeta = 0.0_dp
154 : !
155 42 : ext_basis%nset = lbas + 1
156 168 : ALLOCATE (ext_basis%lmin(lbas + 1), ext_basis%lmax(lbas + 1))
157 84 : ALLOCATE (ext_basis%npgf(lbas + 1))
158 160 : DO l = 0, lbas
159 118 : ext_basis%lmin(l + 1) = l
160 118 : ext_basis%lmax(l + 1) = l
161 160 : IF (l <= maxl) THEN
162 : fmin = 10.0_dp
163 : nn = 0
164 704 : DO i = 1, n1
165 602 : xz = z1(i)
166 2448 : DO j = 1, nps(l)
167 1846 : yz = zexs(j, l)
168 1846 : fz = MAX(xz, yz)/MIN(xz, yz)
169 2448 : fmin = MIN(fmin, fz)
170 : END DO
171 704 : IF (fmin > fr1**0.25) THEN
172 342 : nn = nn + 1
173 342 : zeta(nn, l + 1) = xz
174 : END IF
175 : END DO
176 102 : CPASSERT(nn > 0)
177 102 : 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 160 : nn = MAXVAL(ext_basis%npgf)
184 168 : ALLOCATE (ext_basis%zet(nn, lbas + 1))
185 160 : DO i = 1, lbas + 1
186 118 : nn = ext_basis%npgf(i)
187 580 : ext_basis%zet(1:nn, i) = zeta(1:nn, i)
188 : END DO
189 42 : ext_basis%name = "extbas"
190 42 : ext_basis%kind_radius = orb_basis%kind_radius
191 42 : ext_basis%short_kind_radius = orb_basis%short_kind_radius
192 42 : ext_basis%norm_type = orb_basis%norm_type
193 :
194 42 : NULLIFY (p_basis)
195 42 : CALL create_primitive_basis_set(ext_basis, p_basis)
196 42 : CALL combine_basis_sets(gapw_1c_basis, p_basis)
197 :
198 42 : CALL deallocate_gto_basis_set(ext_basis)
199 42 : CALL deallocate_gto_basis_set(p_basis)
200 84 : DEALLOCATE (zmaxs, zexs, nps, z1, z2, zeta)
201 : END IF
202 :
203 1932 : END SUBROUTINE create_1c_basis
204 :
205 : END MODULE gapw_1c_basis_set
|
