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 MI (08.01.2004)
12 : ! **************************************************************************************************
13 : MODULE paw_proj_set_types
14 :
15 : USE ao_util, ONLY: exp_radius,&
16 : gauss_exponent
17 : USE basis_set_types, ONLY: get_gto_basis_set,&
18 : gto_basis_set_type
19 : USE cp_control_types, ONLY: qs_control_type
20 : USE cp_log_handling, ONLY: cp_get_default_logger,&
21 : cp_logger_type
22 : USE cp_output_handling, ONLY: cp_print_key_finished_output,&
23 : cp_print_key_unit_nr
24 : USE input_section_types, ONLY: section_vals_type
25 : USE kinds, ONLY: default_string_length,&
26 : dp
27 : USE mathconstants, ONLY: dfac,&
28 : rootpi
29 : USE mathlib, ONLY: invert_matrix
30 : USE memory_utilities, ONLY: reallocate
31 : USE orbital_pointers, ONLY: indco,&
32 : indso,&
33 : nco,&
34 : ncoset,&
35 : nso,&
36 : nsoset
37 : USE orbital_transformation_matrices, ONLY: orbtramat
38 : #include "./base/base_uses.f90"
39 :
40 : IMPLICIT NONE
41 :
42 : PRIVATE
43 :
44 : ! Global parameters (only in this module)
45 :
46 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'paw_proj_set_types'
47 :
48 : INTEGER, PARAMETER :: max_name_length = 60
49 :
50 : ! Define the projector types
51 :
52 : TYPE paw_proj_set_type
53 : INTEGER :: maxl = -1, ncgauprj = -1, nsgauprj = -1
54 : INTEGER, DIMENSION(:), POINTER :: nprj => NULL() ! 0:maxl
55 : INTEGER, DIMENSION(:), POINTER :: lx => NULL(), ly => NULL(), lz => NULL() ! ncgauprj
56 : INTEGER, DIMENSION(:), POINTER :: ll => NULL(), m => NULL() ! nsgauprj
57 : INTEGER, DIMENSION(:), POINTER :: first_prj => NULL(), last_prj => NULL() ! 0:maxl
58 : INTEGER, DIMENSION(:), POINTER :: first_prjs => NULL() ! 0:maxl
59 : REAL(KIND=dp) :: rcprj = 0.0_dp
60 : REAL(KIND=dp), DIMENSION(:), POINTER :: zisomin => NULL()
61 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: zetprj => NULL() ! maxnprj,0:maxl
62 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: rzetprj => NULL() ! maxnprj,0:maxl
63 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: csprj => NULL() ! ncgauprj, np_so
64 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: chprj => NULL() ! ncgauprj, np_so
65 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: local_oce_sphi_h => NULL(), local_oce_sphi_s => NULL() ! maxco,nsgf
66 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: sphi_h => NULL(), sphi_s => NULL()
67 : LOGICAL, DIMENSION(:, :), POINTER :: isoprj => NULL() ! maxnprj,0:maxl
68 : INTEGER :: nsatbas = -1
69 : INTEGER :: nsotot = -1
70 : INTEGER, DIMENSION(:), POINTER :: o2nindex => NULL() ! maxso*nset
71 : INTEGER, DIMENSION(:), POINTER :: n2oindex => NULL() ! maxso*nset
72 :
73 : END TYPE paw_proj_set_type
74 :
75 : ! Public subroutines
76 :
77 : PUBLIC :: allocate_paw_proj_set, &
78 : deallocate_paw_proj_set, &
79 : get_paw_proj_set, &
80 : projectors, &
81 : set_paw_proj_set
82 :
83 : ! Public data types
84 :
85 : PUBLIC :: paw_proj_set_type
86 :
87 : CONTAINS
88 :
89 : ! **************************************************************************************************
90 : !> \brief Allocate projector type for GAPW
91 : !> \param paw_proj_set ...
92 : !> \version 1.0
93 : ! **************************************************************************************************
94 1618 : SUBROUTINE allocate_paw_proj_set(paw_proj_set)
95 :
96 : TYPE(paw_proj_set_type), POINTER :: paw_proj_set
97 :
98 1618 : IF (ASSOCIATED(paw_proj_set)) CALL deallocate_paw_proj_set(paw_proj_set)
99 :
100 1618 : ALLOCATE (paw_proj_set)
101 :
102 : NULLIFY (paw_proj_set%nprj)
103 : NULLIFY (paw_proj_set%lx)
104 : NULLIFY (paw_proj_set%ly)
105 : NULLIFY (paw_proj_set%lz)
106 : NULLIFY (paw_proj_set%ll)
107 : NULLIFY (paw_proj_set%m)
108 : NULLIFY (paw_proj_set%first_prj)
109 : NULLIFY (paw_proj_set%last_prj)
110 : NULLIFY (paw_proj_set%first_prjs)
111 :
112 : NULLIFY (paw_proj_set%zisomin)
113 : NULLIFY (paw_proj_set%zetprj)
114 : NULLIFY (paw_proj_set%csprj)
115 : NULLIFY (paw_proj_set%chprj)
116 : NULLIFY (paw_proj_set%local_oce_sphi_h)
117 : NULLIFY (paw_proj_set%local_oce_sphi_s)
118 : NULLIFY (paw_proj_set%sphi_h)
119 : NULLIFY (paw_proj_set%sphi_s)
120 : NULLIFY (paw_proj_set%rzetprj)
121 :
122 : NULLIFY (paw_proj_set%isoprj)
123 :
124 : NULLIFY (paw_proj_set%o2nindex)
125 : NULLIFY (paw_proj_set%n2oindex)
126 :
127 1618 : END SUBROUTINE allocate_paw_proj_set
128 :
129 : ! **************************************************************************************************
130 : !> \brief Deallocate a projector-type set data set.
131 : !> \param paw_proj_set ...
132 : !> \version 1.0
133 : ! **************************************************************************************************
134 1618 : SUBROUTINE deallocate_paw_proj_set(paw_proj_set)
135 : TYPE(paw_proj_set_type), POINTER :: paw_proj_set
136 :
137 1618 : IF (ASSOCIATED(paw_proj_set)) THEN
138 :
139 1618 : IF (ASSOCIATED(paw_proj_set%zisomin)) DEALLOCATE (paw_proj_set%zisomin)
140 1618 : IF (ASSOCIATED(paw_proj_set%nprj)) DEALLOCATE (paw_proj_set%nprj)
141 1618 : IF (ASSOCIATED(paw_proj_set%lx)) DEALLOCATE (paw_proj_set%lx)
142 1618 : IF (ASSOCIATED(paw_proj_set%ly)) DEALLOCATE (paw_proj_set%ly)
143 1618 : IF (ASSOCIATED(paw_proj_set%lz)) DEALLOCATE (paw_proj_set%lz)
144 1618 : IF (ASSOCIATED(paw_proj_set%ll)) DEALLOCATE (paw_proj_set%ll)
145 1618 : IF (ASSOCIATED(paw_proj_set%m)) DEALLOCATE (paw_proj_set%m)
146 1618 : IF (ASSOCIATED(paw_proj_set%first_prj)) DEALLOCATE (paw_proj_set%first_prj)
147 1618 : IF (ASSOCIATED(paw_proj_set%last_prj)) DEALLOCATE (paw_proj_set%last_prj)
148 1618 : IF (ASSOCIATED(paw_proj_set%first_prjs)) DEALLOCATE (paw_proj_set%first_prjs)
149 1618 : IF (ASSOCIATED(paw_proj_set%zetprj)) DEALLOCATE (paw_proj_set%zetprj)
150 1618 : IF (ASSOCIATED(paw_proj_set%csprj)) DEALLOCATE (paw_proj_set%csprj)
151 1618 : IF (ASSOCIATED(paw_proj_set%chprj)) DEALLOCATE (paw_proj_set%chprj)
152 1618 : IF (ASSOCIATED(paw_proj_set%local_oce_sphi_h)) DEALLOCATE (paw_proj_set%local_oce_sphi_h)
153 1618 : IF (ASSOCIATED(paw_proj_set%local_oce_sphi_s)) DEALLOCATE (paw_proj_set%local_oce_sphi_s)
154 1618 : IF (ASSOCIATED(paw_proj_set%sphi_h)) DEALLOCATE (paw_proj_set%sphi_h)
155 1618 : IF (ASSOCIATED(paw_proj_set%sphi_s)) DEALLOCATE (paw_proj_set%sphi_s)
156 1618 : IF (ASSOCIATED(paw_proj_set%isoprj)) DEALLOCATE (paw_proj_set%isoprj)
157 1618 : IF (ASSOCIATED(paw_proj_set%rzetprj)) DEALLOCATE (paw_proj_set%rzetprj)
158 1618 : IF (ASSOCIATED(paw_proj_set%o2nindex)) DEALLOCATE (paw_proj_set%o2nindex)
159 1618 : IF (ASSOCIATED(paw_proj_set%n2oindex)) DEALLOCATE (paw_proj_set%n2oindex)
160 :
161 1618 : DEALLOCATE (paw_proj_set)
162 :
163 : END IF
164 :
165 1618 : END SUBROUTINE deallocate_paw_proj_set
166 :
167 : ! **************************************************************************************************
168 : !> \brief Initialize the projector-type set data set.
169 : !> \param paw_proj ...
170 : !> \param basis_1c Basis set used for the one-center expansions
171 : !> \param orb_basis Orbital basis set
172 : !> \param rc ...
173 : !> \param qs_control ...
174 : !> \param max_rad_local_type ...
175 : !> \param force_env_section ...
176 : !> \version 1.0
177 : ! **************************************************************************************************
178 1618 : SUBROUTINE projectors(paw_proj, basis_1c, orb_basis, rc, qs_control, max_rad_local_type, &
179 : force_env_section)
180 :
181 : TYPE(paw_proj_set_type), POINTER :: paw_proj
182 : TYPE(gto_basis_set_type), POINTER :: basis_1c, orb_basis
183 : REAL(KIND=dp) :: rc
184 : TYPE(qs_control_type), INTENT(IN) :: qs_control
185 : REAL(KIND=dp), INTENT(IN) :: max_rad_local_type
186 : TYPE(section_vals_type), POINTER :: force_env_section
187 :
188 : REAL(KIND=dp) :: eps_fit, eps_iso, eps_orb, eps_svd, &
189 : max_rad_local
190 :
191 1618 : eps_fit = qs_control%gapw_control%eps_fit
192 1618 : eps_iso = qs_control%gapw_control%eps_iso
193 1618 : eps_svd = qs_control%gapw_control%eps_svd
194 1618 : max_rad_local = qs_control%gapw_control%max_rad_local
195 1618 : IF (max_rad_local_type .LT. max_rad_local) THEN
196 1618 : max_rad_local = max_rad_local_type
197 : END IF
198 1618 : eps_orb = qs_control%eps_pgf_orb
199 :
200 : CALL build_projector(paw_proj, basis_1c, orb_basis, eps_fit, eps_iso, eps_svd, &
201 1618 : rc, eps_orb, max_rad_local, force_env_section)
202 :
203 1618 : END SUBROUTINE projectors
204 :
205 : ! **************************************************************************************************
206 : !> \brief initialize the projector-type set data set.
207 : !> \param paw_proj ...
208 : !> \param basis_1c Basis set used for the one-center expansions
209 : !> \param orb_basis Orbital basis set
210 : !> \param eps_fit ...
211 : !> \param eps_iso ...
212 : !> \param eps_svd ...
213 : !> \param rc ...
214 : !> \param eps_orb ...
215 : !> \param max_rad_local To eliminate very smooth functions from the 1c basis
216 : !> \param force_env_section ...
217 : !> \version 1.0
218 : ! **************************************************************************************************
219 1618 : SUBROUTINE build_projector(paw_proj, basis_1c, orb_basis, eps_fit, eps_iso, eps_svd, &
220 : rc, eps_orb, max_rad_local, force_env_section)
221 :
222 : TYPE(paw_proj_set_type), POINTER :: paw_proj
223 : TYPE(gto_basis_set_type), POINTER :: basis_1c, orb_basis
224 : REAL(KIND=dp), INTENT(IN) :: eps_fit, eps_iso, eps_svd, rc, eps_orb, &
225 : max_rad_local
226 : TYPE(section_vals_type), POINTER :: force_env_section
227 :
228 : CHARACTER(LEN=default_string_length) :: bsname
229 : INTEGER :: ic, ico, icomax, icomin, il, info, ip, ipgf, ipp, iprjfirst, iprjs, is, iset, &
230 : isgf, isgfmax, isgfmin, ishell, iso, iso_pgf, iso_set, isomin, jp, k, lshell, lwork, lx, &
231 : ly, lz, maxco, maxl, maxnprj, maxpgf, maxso, mp, ms, n, ncgauprj, ncgf, ncgfo, nisop, np, &
232 : npgfg, ns, nset, nseta, nsgauprj, nsgf, nsgfo, nsox, output_unit
233 1618 : INTEGER, ALLOCATABLE, DIMENSION(:) :: IWORK
234 1618 : INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf, nshell
235 1618 : INTEGER, DIMENSION(:, :), POINTER :: first_cgf, first_sgf, l, last_cgf, &
236 1618 : last_sgf
237 1618 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: isoprj
238 : REAL(KIND=dp) :: expzet, my_error, prefac, radius, x, &
239 : zetmin, zetval
240 1618 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: S, Work_dgesdd
241 1618 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: U, VT
242 1618 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius, zet, zetp
243 1618 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cprj_h, cprj_s, gcc, gcch, smat, sphi, &
244 1618 : work, zetb
245 1618 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: gcca, set_radius2
246 1618 : REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: gcchprj, gccprj
247 : TYPE(cp_logger_type), POINTER :: logger
248 :
249 1618 : NULLIFY (logger)
250 3236 : logger => cp_get_default_logger()
251 :
252 1618 : NULLIFY (first_cgf, first_sgf, last_cgf, last_sgf, gcc, l, set_radius, set_radius2)
253 1618 : NULLIFY (sphi, lmax, lmin, npgf, nshell, zet, zetb, zetp, smat, work, gcca)
254 :
255 1618 : CPASSERT(ASSOCIATED(paw_proj))
256 1618 : CPASSERT(ASSOCIATED(orb_basis))
257 1618 : CPASSERT(ASSOCIATED(basis_1c))
258 :
259 : CALL get_gto_basis_set(gto_basis_set=basis_1c, name=bsname, &
260 : ncgf=ncgf, nset=nset, nsgf=nsgf, &
261 : lmax=lmax, lmin=lmin, npgf=npgf, &
262 : nshell=nshell, sphi=sphi, &
263 : first_cgf=first_cgf, first_sgf=first_sgf, &
264 : l=l, last_cgf=last_cgf, last_sgf=last_sgf, &
265 : maxco=maxco, maxso=maxso, maxl=maxl, maxpgf=maxpgf, &
266 1618 : zet=zetb, gcc=gcca)
267 :
268 1618 : paw_proj%maxl = maxl
269 1618 : CPASSERT(.NOT. ASSOCIATED(paw_proj%zisomin))
270 1618 : CPASSERT(.NOT. ASSOCIATED(paw_proj%nprj))
271 :
272 4854 : ALLOCATE (paw_proj%zisomin(0:maxl))
273 5462 : paw_proj%zisomin(0:maxl) = 0.0_dp
274 4854 : ALLOCATE (paw_proj%nprj(0:maxl))
275 5462 : paw_proj%nprj(0:maxl) = 0
276 :
277 : output_unit = cp_print_key_unit_nr(logger, force_env_section, &
278 1618 : "DFT%PRINT%GAPW%PROJECTORS", extension=".Log")
279 :
280 1618 : IF (output_unit > 0) THEN
281 : WRITE (UNIT=output_unit, FMT="(/,T2,A)") &
282 0 : "Projectors for the basis functions of "//TRIM(bsname)
283 : END IF
284 :
285 4854 : ALLOCATE (set_radius(nset))
286 6724 : set_radius = 0.0_dp
287 6724 : DO iset = 1, nset
288 14972 : DO is = 1, nshell(iset)
289 : set_radius(iset) = MAX(set_radius(iset), &
290 : exp_radius(l(is, iset), zetb(npgf(iset), iset), &
291 : eps_orb, gcca(npgf(iset), is, iset), &
292 13354 : rlow=set_radius(iset)))
293 : END DO ! is
294 : END DO ! iset
295 :
296 8090 : ALLOCATE (set_radius2(maxpgf, 0:maxl, nset))
297 88672 : set_radius2 = 0.0_dp
298 6724 : DO iset = 1, nset
299 13188 : DO lshell = lmin(iset), lmax(iset)
300 28008 : DO ip = 1, npgf(iset)
301 : set_radius2(ip, lshell, iset) = &
302 22902 : exp_radius(lshell, zetb(ip, iset), eps_orb, 1.0_dp)
303 : END DO
304 : END DO ! is
305 : END DO ! iset
306 :
307 1618 : maxnprj = 0
308 5462 : DO lshell = 0, maxl ! lshell
309 3844 : np = 0
310 17064 : DO iset = 1, nset
311 17064 : IF (lshell >= lmin(iset) .AND. lshell <= lmax(iset)) THEN
312 22902 : DO ip = 1, npgf(iset)
313 22902 : IF (set_radius2(ip, lshell, iset) < max_rad_local) THEN
314 16304 : np = np + 1
315 : END IF
316 : END DO
317 : END IF
318 : END DO
319 3844 : maxnprj = MAX(maxnprj, np)
320 3844 : paw_proj%nprj(lshell) = np
321 5462 : IF (np < 1) THEN
322 0 : CPABORT("No Projector for lshell found")
323 : END IF
324 : END DO ! lshell
325 :
326 : ! Allocate exponents and coefficients
327 6472 : ALLOCATE (paw_proj%zetprj(maxnprj, 0:maxl))
328 30946 : paw_proj%zetprj(1:maxnprj, 0:maxl) = 0.0_dp
329 4854 : ALLOCATE (paw_proj%rzetprj(maxnprj, 0:maxl))
330 30946 : paw_proj%rzetprj(1:maxnprj, 0:maxl) = 0.0_dp
331 6472 : ALLOCATE (paw_proj%isoprj(maxnprj, 0:maxl))
332 30946 : paw_proj%isoprj = .FALSE.
333 9708 : ALLOCATE (gccprj(maxnprj, maxpgf, 0:maxl, nset))
334 711562 : gccprj = 0.0_dp
335 8090 : ALLOCATE (gcchprj(maxnprj, maxpgf, 0:maxl, nset))
336 711562 : gcchprj = 0.0_dp
337 :
338 1618 : NULLIFY (zet, zetp, gcc, smat, work)
339 : ! Generate the projetor basis for each ang. mom. q.n.
340 5462 : DO lshell = 0, maxl ! lshell
341 :
342 3844 : np = paw_proj%nprj(lshell)
343 :
344 11532 : ALLOCATE (isoprj(np))
345 20148 : isoprj = .FALSE.
346 :
347 49972 : ALLOCATE (zet(np), zetp(np), gcc(np, np), gcch(np, np), smat(np, np), work(np, np))
348 :
349 20148 : zet(:) = 0.0_dp
350 20148 : zetp(:) = 0.0_dp
351 123528 : gcc(:, :) = 0.0_dp
352 123528 : gcch(:, :) = 0.0_dp
353 123528 : smat(:, :) = 0.0_dp
354 123528 : work(:, :) = 0.0_dp
355 :
356 3844 : npgfg = 0
357 : ! Collect all the exponent which contribute to lshell
358 17064 : DO iset = 1, nset ! iset
359 17064 : IF (lshell >= lmin(iset) .AND. lshell <= lmax(iset)) THEN
360 22902 : DO ip = 1, npgf(iset)
361 22902 : IF (set_radius2(ip, lshell, iset) < max_rad_local) THEN
362 16304 : npgfg = npgfg + 1
363 16304 : zet(npgfg) = zetb(ip, iset)
364 : END IF
365 : END DO
366 : END IF
367 : END DO ! iset
368 :
369 : ! *** Smallest exp. due to eps_iso: concerned as an isolated projector ***
370 3844 : paw_proj%zisomin(lshell) = gauss_exponent(lshell, rc, eps_iso, 1.0_dp)
371 :
372 : ! maybe order the exponents here?
373 : ! zet(1) > zet(2) ...
374 : !
375 3844 : nisop = 0
376 20148 : DO ip = 1, np
377 : ! Check for equal exponents
378 59842 : DO ipp = 1, ip - 1
379 59842 : IF (zet(ip) == zet(ipp)) THEN
380 : CALL cp_abort(__LOCATION__, &
381 : "Linear dependency in the construction of the GAPW projectors:"// &
382 : " different sets of the BASIS SET contain identical exponents"// &
383 0 : " for the same l quantum numbers")
384 : END IF
385 : END DO
386 :
387 20148 : IF (zet(ip) >= paw_proj%zisomin(lshell)) THEN
388 3414 : isoprj(ip) = .TRUE.
389 3414 : nisop = nisop + 1
390 : ELSE
391 12890 : isoprj(ip) = .FALSE.
392 : END IF
393 : END DO
394 :
395 : ! Smallest exp. due to eps_fit: where to start geometric progression
396 3844 : zetmin = gauss_exponent(lshell, rc, eps_fit, 1.0_dp)
397 :
398 : ! Generate the projectors by the geometric progression
399 3844 : IF (np - nisop - 1 > 2) THEN
400 2196 : x = (80.0_dp/zetmin)**(1.0_dp/REAL(np - nisop - 1, dp))
401 : ELSE
402 : x = 2.0_dp
403 : END IF
404 2196 : IF (x > 2.0_dp) x = 2.0_dp
405 :
406 3844 : zetval = zetmin
407 20148 : DO ip = np, 1, -1
408 20148 : IF (.NOT. isoprj(ip)) THEN
409 12890 : zetp(ip) = zetval
410 12890 : zetval = x*zetval
411 : END IF
412 : END DO
413 :
414 3844 : nisop = 0
415 20148 : DO ip = np, 1, -1
416 20148 : IF (isoprj(ip)) THEN
417 3414 : zetp(ip) = zetval
418 3414 : zetval = x*zetval
419 3414 : nisop = nisop + 1
420 : END IF
421 : END DO
422 :
423 : ! *** Build the overlap matrix: <projector|primitive> ***
424 3844 : prefac = 0.5_dp**(lshell + 2)*rootpi*dfac(2*lshell + 1)
425 3844 : expzet = REAL(lshell, dp) + 1.5_dp
426 :
427 20148 : DO ip = 1, np
428 20148 : IF (isoprj(ip)) THEN
429 36726 : DO jp = 1, np
430 36726 : IF (isoprj(jp)) THEN
431 17750 : smat(ip, jp) = prefac/(zetp(ip) + zet(jp))**expzet
432 : END IF
433 : END DO
434 : ELSE
435 82958 : DO jp = 1, np
436 82958 : IF (.NOT. isoprj(jp)) THEN
437 54506 : smat(ip, jp) = prefac/(zetp(ip) + zet(jp))**expzet
438 : END IF
439 : END DO
440 : END IF
441 : END DO
442 :
443 : ! Compute inverse of the transpose
444 3844 : IF (eps_svd .EQ. 0.0_dp) THEN
445 222 : CALL invert_matrix(smat, gcc, my_error, "T")
446 : ELSE
447 229854 : work = TRANSPOSE(smat)
448 : ! Workspace query
449 28976 : ALLOCATE (iwork(8*np), S(np), U(np, np), VT(np, np), work_dgesdd(1))
450 3622 : lwork = -1
451 3622 : CALL DGESDD('S', np, np, work, np, S, U, np, vt, np, work_dgesdd, lwork, iwork, info)
452 3622 : lwork = INT(work_dgesdd(1))
453 10866 : DEALLOCATE (work_dgesdd); ALLOCATE (work_dgesdd(lwork))
454 3622 : CALL DGESDD('S', np, np, work, np, S, U, np, vt, np, work_dgesdd, lwork, iwork, info)
455 : ! Construct the inverse
456 19026 : DO k = 1, np
457 : ! invert SV
458 15404 : IF (S(k) < eps_svd) THEN
459 1084 : S(k) = 0.0_dp
460 : ELSE
461 14320 : S(k) = 1.0_dp/S(k)
462 : END IF
463 116738 : VT(k, :) = VT(k, :)*S(k)
464 : END DO
465 3622 : CALL DGEMM('T', 'T', np, np, np, 1.0_dp, VT, np, U, np, 0.0_dp, gcc, np)
466 3622 : DEALLOCATE (iwork, S, U, VT, work_dgesdd)
467 : END IF
468 :
469 : ! Set the coefficient of the isolated projectors to 0
470 243212 : gcch(:, :) = gcc(:, :)
471 20148 : DO ip = 1, np
472 20148 : IF (isoprj(ip)) THEN
473 36726 : gcc(:, ip) = 0.0_dp
474 36726 : gcc(ip, :) = 0.0_dp
475 : END IF
476 : END DO
477 :
478 : ! Transfer data from local to global variables
479 :
480 36452 : paw_proj%zetprj(1:np, lshell) = zetp(1:np)
481 20148 : paw_proj%isoprj(1:np, lshell) = isoprj(1:np)
482 :
483 3844 : npgfg = 0
484 17064 : DO iset = 1, nset ! iset
485 17064 : IF (lshell >= lmin(iset) .AND. lshell <= lmax(iset)) THEN
486 22902 : DO ip = 1, npgf(iset)
487 22902 : IF (set_radius2(ip, lshell, iset) < max_rad_local) THEN
488 16304 : npgfg = npgfg + 1
489 223064 : gccprj(1:np, ip, lshell, iset) = gcc(1:np, npgfg)
490 223064 : gcchprj(1:np, ip, lshell, iset) = gcch(1:np, npgfg)
491 : ELSE
492 980 : gccprj(1:np, ip, lshell, iset) = 0.0_dp
493 980 : gcchprj(1:np, ip, lshell, iset) = 0.0_dp
494 : END IF
495 : END DO
496 : END IF
497 : END DO ! iset
498 :
499 : ! Print exponents and coefficients of the projectors
500 3844 : IF (output_unit > 0) THEN
501 : WRITE (UNIT=output_unit, FMT="(/,/,T2,A,I2)") &
502 0 : "Built projector for angular momentum quantum number l= ", lshell
503 : WRITE (UNIT=output_unit, FMT="(T2,A,I2)") &
504 0 : "Number of isolated projectors = ", nisop
505 0 : DO iset = 1, nset ! iset
506 0 : IF (lshell >= lmin(iset) .AND. lshell <= lmax(iset)) THEN
507 : WRITE (UNIT=output_unit, FMT="(/,T2,A,I5,/,/,T4,A9,(T13,4f15.6))") &
508 0 : "Set ", iset, "exp prj: ", &
509 0 : (paw_proj%zetprj(ip, lshell), ip=1, np)
510 0 : DO jp = 1, npgf(iset)
511 : WRITE (UNIT=output_unit, FMT="(/,T4,A9,F15.6,/,T4,A9,(t13,4E15.6))") &
512 0 : "exp gto: ", zetb(jp, iset), &
513 0 : "coeff.: ", (gccprj(ip, jp, lshell, iset), ip=1, np)
514 : END DO
515 : END IF
516 : END DO ! iset
517 : END IF
518 :
519 : ! Release the working storage for the current value lshell
520 3844 : DEALLOCATE (isoprj)
521 5462 : DEALLOCATE (gcc, gcch, zet, zetp, smat, work)
522 :
523 : END DO ! lshell
524 : CALL cp_print_key_finished_output(output_unit, logger, force_env_section, &
525 1618 : "DFT%PRINT%GAPW%PROJECTORS")
526 :
527 : ! Release the working storage for the current value lshell
528 1618 : DEALLOCATE (set_radius)
529 1618 : DEALLOCATE (set_radius2)
530 :
531 : ! Count primitives basis functions for the projectors
532 1618 : paw_proj%ncgauprj = 0
533 1618 : paw_proj%nsgauprj = 0
534 5462 : DO lshell = 0, maxl
535 3844 : paw_proj%ncgauprj = paw_proj%ncgauprj + nco(lshell)*paw_proj%nprj(lshell)
536 5462 : paw_proj%nsgauprj = paw_proj%nsgauprj + nso(lshell)*paw_proj%nprj(lshell)
537 : END DO
538 :
539 1618 : ncgauprj = paw_proj%ncgauprj
540 1618 : nsgauprj = paw_proj%nsgauprj
541 1618 : CALL reallocate(paw_proj%lx, 1, ncgauprj)
542 1618 : CALL reallocate(paw_proj%ly, 1, ncgauprj)
543 1618 : CALL reallocate(paw_proj%lz, 1, ncgauprj)
544 1618 : CALL reallocate(paw_proj%first_prj, 0, maxl)
545 1618 : CALL reallocate(paw_proj%last_prj, 0, maxl)
546 1618 : CALL reallocate(paw_proj%ll, 1, nsgauprj)
547 1618 : CALL reallocate(paw_proj%m, 1, nsgauprj)
548 1618 : CALL reallocate(paw_proj%first_prjs, 0, maxl)
549 :
550 6472 : ALLOCATE (cprj_s(1:nsgauprj, 1:maxso*nset))
551 4854 : ALLOCATE (cprj_h(1:nsgauprj, 1:maxso*nset))
552 2892622 : cprj_s = 0.0_dp
553 2892622 : cprj_h = 0.0_dp
554 :
555 1618 : ncgauprj = 0
556 1618 : nsgauprj = 0
557 5462 : DO lshell = 0, maxl
558 3844 : np = paw_proj%nprj(lshell)
559 3844 : paw_proj%first_prj(lshell) = ncgauprj + 1
560 3844 : paw_proj%first_prjs(lshell) = nsgauprj + 1
561 3844 : paw_proj%last_prj(lshell) = ncgauprj + nco(lshell)*np
562 21766 : DO ip = 1, np
563 49454 : DO ico = ncoset(lshell - 1) + 1, ncoset(lshell)
564 33150 : ncgauprj = ncgauprj + 1
565 33150 : paw_proj%lx(ncgauprj) = indco(1, ico)
566 33150 : paw_proj%ly(ncgauprj) = indco(2, ico)
567 49454 : paw_proj%lz(ncgauprj) = indco(3, ico)
568 : END DO ! ico
569 51820 : DO iso = nsoset(lshell - 1) + 1, nsoset(lshell)
570 31672 : nsgauprj = nsgauprj + 1
571 31672 : paw_proj%ll(nsgauprj) = indso(1, iso)
572 47976 : paw_proj%m(nsgauprj) = indso(2, iso)
573 : END DO
574 : END DO ! ip
575 : END DO ! lshell
576 :
577 1618 : ms = 0
578 6724 : DO iset = 1, nset
579 5106 : ns = nsoset(lmax(iset))
580 11570 : DO lshell = lmin(iset), lmax(iset)
581 6464 : iprjfirst = paw_proj%first_prjs(lshell)
582 6464 : np = paw_proj%nprj(lshell)
583 28008 : DO ipgf = 1, npgf(iset)
584 127128 : DO ip = 1, np
585 291602 : DO il = 1, nso(lshell)
586 170938 : iprjs = iprjfirst - 1 + il + (ip - 1)*nso(lshell)
587 170938 : iso = nsoset(lshell - 1) + 1 + (lshell + paw_proj%m(iprjs))
588 :
589 170938 : iso = iso + (ipgf - 1)*ns + ms
590 170938 : cprj_s(iprjs, iso) = gccprj(ip, ipgf, lshell, iset)
591 275164 : cprj_h(iprjs, iso) = gcchprj(ip, ipgf, lshell, iset)
592 : END DO ! iprjs
593 : END DO ! ip
594 : END DO ! ipgf
595 : END DO ! lshell
596 6724 : ms = ms + maxso
597 : END DO ! iset
598 :
599 : ! Local coefficients for the one center expansions : oce
600 : ! the coefficients are calculated for the full and soft expansions
601 : CALL get_gto_basis_set(gto_basis_set=orb_basis, &
602 1618 : nset=nseta, ncgf=ncgfo, nsgf=nsgfo)
603 :
604 6472 : ALLOCATE (paw_proj%local_oce_sphi_h(maxco, nsgfo))
605 439448 : paw_proj%local_oce_sphi_h = 0.0_dp
606 4854 : ALLOCATE (paw_proj%sphi_h(maxco, nsgfo))
607 439448 : paw_proj%sphi_h = 0.0_dp
608 :
609 4854 : ALLOCATE (paw_proj%local_oce_sphi_s(maxco, nsgfo))
610 439448 : paw_proj%local_oce_sphi_s = 0.0_dp
611 4854 : ALLOCATE (paw_proj%sphi_s(maxco, nsgfo))
612 439448 : paw_proj%sphi_s = 0.0_dp
613 :
614 : ! only use first nset of orb basis local projection!
615 6648 : DO iset = 1, nseta
616 5030 : n = ncoset(lmax(iset))
617 18982 : DO ipgf = 1, npgf(iset)
618 43288 : DO ishell = 1, nshell(iset)
619 25924 : lshell = l(ishell, iset)
620 25924 : icomin = ncoset(lshell - 1) + 1 + n*(ipgf - 1)
621 25924 : icomax = ncoset(lshell) + n*(ipgf - 1)
622 25924 : isgfmin = first_sgf(ishell, iset)
623 25924 : isgfmax = last_sgf(ishell, iset)
624 : radius = exp_radius(lshell, basis_1c%zet(ipgf, iset), &
625 25924 : eps_fit, 1.0_dp)
626 90050 : DO isgf = isgfmin, isgfmax
627 : paw_proj%sphi_h(icomin:icomax, isgf) = &
628 394848 : sphi(icomin:icomax, isgf)
629 77716 : IF (radius < rc) THEN
630 78268 : paw_proj%sphi_s(icomin:icomax, isgf) = 0.0_dp
631 : ELSE
632 : paw_proj%sphi_s(icomin:icomax, isgf) = &
633 256906 : sphi(icomin:icomax, isgf)
634 : END IF
635 : END DO
636 : END DO ! ishell
637 : END DO ! ipgf
638 : END DO ! iset
639 :
640 : ! only use first nset of orb basis local projection!
641 6648 : DO iset = 1, nseta
642 5030 : n = ncoset(lmax(iset))
643 5030 : ns = nsoset(lmax(iset))
644 18982 : DO ipgf = 1, npgf(iset)
645 43288 : DO ishell = 1, nshell(iset)
646 25924 : lshell = l(ishell, iset)
647 25924 : icomin = ncoset(lshell - 1) + 1 + n*(ipgf - 1)
648 25924 : icomax = ncoset(lshell) + n*(ipgf - 1)
649 25924 : isgfmin = first_sgf(ishell, iset)
650 25924 : isgfmax = last_sgf(ishell, iset)
651 25924 : isomin = nsoset(lshell - 1) + 1 + ns*(ipgf - 1)
652 90050 : DO is = 1, nso(lshell)
653 51792 : iso = isomin + is - 1
654 249244 : DO ic = 1, nco(lshell)
655 171528 : ico = icomin + ic - 1
656 171528 : lx = indco(1, ic + ncoset(lshell - 1))
657 171528 : ly = indco(2, ic + ncoset(lshell - 1))
658 171528 : lz = indco(3, ic + ncoset(lshell - 1))
659 934644 : DO isgf = isgfmin, isgfmax
660 : paw_proj%local_oce_sphi_h(iso, isgf) = &
661 : paw_proj%local_oce_sphi_h(iso, isgf) + &
662 711324 : orbtramat(lshell)%slm_inv(is, ic)*paw_proj%sphi_h(ico, isgf)
663 : paw_proj%local_oce_sphi_s(iso, isgf) = &
664 : paw_proj%local_oce_sphi_s(iso, isgf) + &
665 882852 : orbtramat(lshell)%slm_inv(is, ic)*paw_proj%sphi_s(ico, isgf)
666 : END DO ! isgf
667 : END DO ! ic
668 : END DO ! is
669 : END DO ! ishell
670 : END DO ! ipgf
671 : END DO ! iset
672 :
673 : ! Index transformation OLD-NEW
674 4854 : ALLOCATE (paw_proj%o2nindex(maxso*nset))
675 3236 : ALLOCATE (paw_proj%n2oindex(maxso*nset))
676 65830 : paw_proj%o2nindex = 0
677 65830 : paw_proj%n2oindex = 0
678 1618 : ico = 1
679 6724 : DO iset = 1, nset
680 5106 : iso_set = (iset - 1)*maxso + 1
681 5106 : nsox = nsoset(lmax(iset))
682 19332 : DO ipgf = 1, npgf(iset)
683 12608 : iso_pgf = iso_set + (ipgf - 1)*nsox
684 12608 : iso = iso_pgf + nsoset(lmin(iset) - 1)
685 34152 : DO lx = lmin(iset), lmax(iset)
686 61008 : DO k = 1, nso(lx)
687 31962 : paw_proj%n2oindex(ico) = iso
688 31962 : paw_proj%o2nindex(iso) = ico
689 31962 : iso = iso + 1
690 48400 : ico = ico + 1
691 : END DO
692 : END DO
693 : END DO
694 : END DO
695 1618 : mp = ico - 1
696 1618 : paw_proj%nsatbas = mp
697 1618 : paw_proj%nsotot = nset*maxso
698 6472 : ALLOCATE (paw_proj%csprj(nsgauprj, mp))
699 1118356 : paw_proj%csprj = 0.0_dp
700 33580 : DO k = 1, mp
701 31962 : ico = paw_proj%n2oindex(k)
702 2203132 : paw_proj%csprj(:, k) = cprj_s(:, ico)
703 : END DO
704 4854 : ALLOCATE (paw_proj%chprj(nsgauprj, mp))
705 1118356 : paw_proj%chprj = 0.0_dp
706 33580 : DO k = 1, mp
707 31962 : ico = paw_proj%n2oindex(k)
708 2203132 : paw_proj%chprj(:, k) = cprj_h(:, ico)
709 : END DO
710 1618 : DEALLOCATE (cprj_s, cprj_h, gcchprj, gccprj)
711 :
712 6472 : END SUBROUTINE build_projector
713 :
714 : ! **************************************************************************************************
715 : !> \brief Get informations about a paw projectors set.
716 : !> \param paw_proj_set ...
717 : !> \param csprj ...
718 : !> \param chprj ...
719 : !> \param first_prj ...
720 : !> \param first_prjs ...
721 : !> \param last_prj ...
722 : !> \param local_oce_sphi_h ...
723 : !> \param local_oce_sphi_s ...
724 : !> \param maxl ...
725 : !> \param ncgauprj ...
726 : !> \param nsgauprj ...
727 : !> \param nsatbas ...
728 : !> \param nsotot ...
729 : !> \param nprj ...
730 : !> \param o2nindex ...
731 : !> \param n2oindex ...
732 : !> \param rcprj ...
733 : !> \param rzetprj ...
734 : !> \param zisomin ...
735 : !> \param zetprj ...
736 : !> \version 1.0
737 : ! **************************************************************************************************
738 241630 : SUBROUTINE get_paw_proj_set(paw_proj_set, csprj, chprj, &
739 : first_prj, first_prjs, last_prj, &
740 : local_oce_sphi_h, local_oce_sphi_s, &
741 : maxl, ncgauprj, nsgauprj, nsatbas, nsotot, nprj, &
742 : o2nindex, n2oindex, &
743 : rcprj, rzetprj, zisomin, zetprj)
744 :
745 : TYPE(paw_proj_set_type), POINTER :: paw_proj_set
746 : REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: csprj, chprj
747 : INTEGER, DIMENSION(:), OPTIONAL, POINTER :: first_prj, first_prjs, last_prj
748 : REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: local_oce_sphi_h, local_oce_sphi_s
749 : INTEGER, INTENT(OUT), OPTIONAL :: maxl, ncgauprj, nsgauprj, nsatbas, nsotot
750 : INTEGER, DIMENSION(:), OPTIONAL, POINTER :: nprj, o2nindex, n2oindex
751 : REAL(KIND=dp), INTENT(OUT), OPTIONAL :: rcprj
752 : REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: rzetprj
753 : REAL(KIND=dp), DIMENSION(:), OPTIONAL, POINTER :: zisomin
754 : REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: zetprj
755 :
756 241630 : IF (ASSOCIATED(paw_proj_set)) THEN
757 241630 : IF (PRESENT(csprj)) csprj => paw_proj_set%csprj
758 241630 : IF (PRESENT(chprj)) chprj => paw_proj_set%chprj
759 241630 : IF (PRESENT(local_oce_sphi_h)) local_oce_sphi_h => paw_proj_set%local_oce_sphi_h
760 241630 : IF (PRESENT(local_oce_sphi_s)) local_oce_sphi_s => paw_proj_set%local_oce_sphi_s
761 241630 : IF (PRESENT(first_prj)) first_prj => paw_proj_set%first_prj
762 241630 : IF (PRESENT(last_prj)) last_prj => paw_proj_set%last_prj
763 241630 : IF (PRESENT(first_prjs)) first_prjs => paw_proj_set%first_prjs
764 241630 : IF (PRESENT(maxl)) maxl = paw_proj_set%maxl
765 241630 : IF (PRESENT(ncgauprj)) ncgauprj = paw_proj_set%ncgauprj
766 241630 : IF (PRESENT(nsgauprj)) nsgauprj = paw_proj_set%nsgauprj
767 241630 : IF (PRESENT(nsatbas)) nsatbas = paw_proj_set%nsatbas
768 241630 : IF (PRESENT(nsotot)) nsotot = paw_proj_set%nsotot
769 241630 : IF (PRESENT(nprj)) nprj => paw_proj_set%nprj
770 241630 : IF (PRESENT(rcprj)) rcprj = paw_proj_set%rcprj
771 241630 : IF (PRESENT(rzetprj)) rzetprj => paw_proj_set%rzetprj
772 241630 : IF (PRESENT(zisomin)) zisomin => paw_proj_set%zisomin
773 241630 : IF (PRESENT(zetprj)) zetprj => paw_proj_set%zetprj
774 241630 : IF (PRESENT(o2nindex)) o2nindex => paw_proj_set%o2nindex
775 241630 : IF (PRESENT(n2oindex)) n2oindex => paw_proj_set%n2oindex
776 : ELSE
777 0 : CPABORT("The pointer paw_proj_set is not associated")
778 : END IF
779 :
780 241630 : END SUBROUTINE get_paw_proj_set
781 :
782 : ! **************************************************************************************************
783 : !> \brief Set informations about a paw projectors set.
784 : !> \param paw_proj_set ...
785 : !> \param rzetprj ...
786 : !> \param rcprj ...
787 : !> \version 1.0
788 : ! **************************************************************************************************
789 1682 : SUBROUTINE set_paw_proj_set(paw_proj_set, rzetprj, rcprj)
790 :
791 : TYPE(paw_proj_set_type), POINTER :: paw_proj_set
792 : REAL(KIND=dp), DIMENSION(:, :), OPTIONAL, POINTER :: rzetprj
793 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: rcprj
794 :
795 1682 : IF (ASSOCIATED(paw_proj_set)) THEN
796 63556 : IF (PRESENT(rzetprj)) paw_proj_set%rzetprj(:, 0:) = rzetprj(:, 0:)
797 1682 : IF (PRESENT(rcprj)) paw_proj_set%rcprj = rcprj
798 : ELSE
799 0 : CPABORT("The pointer paw_proj_set is not associated")
800 : END IF
801 :
802 1682 : END SUBROUTINE set_paw_proj_set
803 :
804 0 : END MODULE paw_proj_set_types
|
