(git:33f85d8)
Loading...
Searching...
No Matches
paw_proj_set_types.F
Go to the documentation of this file.
1!--------------------------------------------------------------------------------------------------!
2! CP2K: A general program to perform molecular dynamics simulations !
3! Copyright 2000-2025 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! **************************************************************************************************
13MODULE 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
87CONTAINS
88
89! **************************************************************************************************
90!> \brief Allocate projector type for GAPW
91!> \param paw_proj_set ...
92!> \version 1.0
93! **************************************************************************************************
94 SUBROUTINE allocate_paw_proj_set(paw_proj_set)
95
96 TYPE(paw_proj_set_type), POINTER :: paw_proj_set
97
98 IF (ASSOCIATED(paw_proj_set)) CALL deallocate_paw_proj_set(paw_proj_set)
99
100 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 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 SUBROUTINE deallocate_paw_proj_set(paw_proj_set)
135 TYPE(paw_proj_set_type), POINTER :: paw_proj_set
136
137 IF (ASSOCIATED(paw_proj_set)) THEN
138
139 IF (ASSOCIATED(paw_proj_set%zisomin)) DEALLOCATE (paw_proj_set%zisomin)
140 IF (ASSOCIATED(paw_proj_set%nprj)) DEALLOCATE (paw_proj_set%nprj)
141 IF (ASSOCIATED(paw_proj_set%lx)) DEALLOCATE (paw_proj_set%lx)
142 IF (ASSOCIATED(paw_proj_set%ly)) DEALLOCATE (paw_proj_set%ly)
143 IF (ASSOCIATED(paw_proj_set%lz)) DEALLOCATE (paw_proj_set%lz)
144 IF (ASSOCIATED(paw_proj_set%ll)) DEALLOCATE (paw_proj_set%ll)
145 IF (ASSOCIATED(paw_proj_set%m)) DEALLOCATE (paw_proj_set%m)
146 IF (ASSOCIATED(paw_proj_set%first_prj)) DEALLOCATE (paw_proj_set%first_prj)
147 IF (ASSOCIATED(paw_proj_set%last_prj)) DEALLOCATE (paw_proj_set%last_prj)
148 IF (ASSOCIATED(paw_proj_set%first_prjs)) DEALLOCATE (paw_proj_set%first_prjs)
149 IF (ASSOCIATED(paw_proj_set%zetprj)) DEALLOCATE (paw_proj_set%zetprj)
150 IF (ASSOCIATED(paw_proj_set%csprj)) DEALLOCATE (paw_proj_set%csprj)
151 IF (ASSOCIATED(paw_proj_set%chprj)) DEALLOCATE (paw_proj_set%chprj)
152 IF (ASSOCIATED(paw_proj_set%local_oce_sphi_h)) DEALLOCATE (paw_proj_set%local_oce_sphi_h)
153 IF (ASSOCIATED(paw_proj_set%local_oce_sphi_s)) DEALLOCATE (paw_proj_set%local_oce_sphi_s)
154 IF (ASSOCIATED(paw_proj_set%sphi_h)) DEALLOCATE (paw_proj_set%sphi_h)
155 IF (ASSOCIATED(paw_proj_set%sphi_s)) DEALLOCATE (paw_proj_set%sphi_s)
156 IF (ASSOCIATED(paw_proj_set%isoprj)) DEALLOCATE (paw_proj_set%isoprj)
157 IF (ASSOCIATED(paw_proj_set%rzetprj)) DEALLOCATE (paw_proj_set%rzetprj)
158 IF (ASSOCIATED(paw_proj_set%o2nindex)) DEALLOCATE (paw_proj_set%o2nindex)
159 IF (ASSOCIATED(paw_proj_set%n2oindex)) DEALLOCATE (paw_proj_set%n2oindex)
160
161 DEALLOCATE (paw_proj_set)
162
163 END IF
164
165 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 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 eps_fit = qs_control%gapw_control%eps_fit
192 eps_iso = qs_control%gapw_control%eps_iso
193 eps_svd = qs_control%gapw_control%eps_svd
194 max_rad_local = qs_control%gapw_control%max_rad_local
195 IF (max_rad_local_type .LT. max_rad_local) THEN
196 max_rad_local = max_rad_local_type
197 END IF
198 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 rc, eps_orb, max_rad_local, force_env_section)
202
203 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 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 INTEGER, ALLOCATABLE, DIMENSION(:) :: iwork
234 INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf, nshell
235 INTEGER, DIMENSION(:, :), POINTER :: first_cgf, first_sgf, l, last_cgf, &
236 last_sgf
237 LOGICAL, ALLOCATABLE, DIMENSION(:) :: isoprj
238 REAL(kind=dp) :: expzet, my_error, prefac, radius, x, &
239 zetmin, zetval
240 REAL(kind=dp), ALLOCATABLE, DIMENSION(:) :: s, work_dgesdd
241 REAL(kind=dp), ALLOCATABLE, DIMENSION(:, :) :: u, vt
242 REAL(kind=dp), DIMENSION(:), POINTER :: set_radius, zet, zetp
243 REAL(kind=dp), DIMENSION(:, :), POINTER :: cprj_h, cprj_s, gcc, gcch, smat, sphi, &
244 work, zetb
245 REAL(kind=dp), DIMENSION(:, :, :), POINTER :: gcca, set_radius2
246 REAL(kind=dp), DIMENSION(:, :, :, :), POINTER :: gcchprj, gccprj
247 TYPE(cp_logger_type), POINTER :: logger
248
249 NULLIFY (logger)
250 logger => cp_get_default_logger()
251
252 NULLIFY (first_cgf, first_sgf, last_cgf, last_sgf, gcc, l, set_radius, set_radius2)
253 NULLIFY (sphi, lmax, lmin, npgf, nshell, zet, zetb, zetp, smat, work, gcca)
254
255 cpassert(ASSOCIATED(paw_proj))
256 cpassert(ASSOCIATED(orb_basis))
257 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 zet=zetb, gcc=gcca)
267
268 paw_proj%maxl = maxl
269 cpassert(.NOT. ASSOCIATED(paw_proj%zisomin))
270 cpassert(.NOT. ASSOCIATED(paw_proj%nprj))
271
272 ALLOCATE (paw_proj%zisomin(0:maxl))
273 paw_proj%zisomin(0:maxl) = 0.0_dp
274 ALLOCATE (paw_proj%nprj(0:maxl))
275 paw_proj%nprj(0:maxl) = 0
276
277 output_unit = cp_print_key_unit_nr(logger, force_env_section, &
278 "DFT%PRINT%GAPW%PROJECTORS", extension=".Log")
279
280 IF (output_unit > 0) THEN
281 WRITE (unit=output_unit, fmt="(/,T2,A)") &
282 "Projectors for the basis functions of "//trim(bsname)
283 END IF
284
285 ALLOCATE (set_radius(nset))
286 set_radius = 0.0_dp
287 DO iset = 1, nset
288 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 rlow=set_radius(iset)))
293 END DO ! is
294 END DO ! iset
295
296 ALLOCATE (set_radius2(maxpgf, 0:maxl, nset))
297 set_radius2 = 0.0_dp
298 DO iset = 1, nset
299 DO lshell = lmin(iset), lmax(iset)
300 DO ip = 1, npgf(iset)
301 set_radius2(ip, lshell, iset) = &
302 exp_radius(lshell, zetb(ip, iset), eps_orb, 1.0_dp)
303 END DO
304 END DO ! is
305 END DO ! iset
306
307 maxnprj = 0
308 DO lshell = 0, maxl ! lshell
309 np = 0
310 DO iset = 1, nset
311 IF (lshell >= lmin(iset) .AND. lshell <= lmax(iset)) THEN
312 DO ip = 1, npgf(iset)
313 IF (set_radius2(ip, lshell, iset) < max_rad_local) THEN
314 np = np + 1
315 END IF
316 END DO
317 END IF
318 END DO
319 maxnprj = max(maxnprj, np)
320 paw_proj%nprj(lshell) = np
321 IF (np < 1) THEN
322 cpabort("No Projector for lshell found")
323 END IF
324 END DO ! lshell
325
326 ! Allocate exponents and coefficients
327 ALLOCATE (paw_proj%zetprj(maxnprj, 0:maxl))
328 paw_proj%zetprj(1:maxnprj, 0:maxl) = 0.0_dp
329 ALLOCATE (paw_proj%rzetprj(maxnprj, 0:maxl))
330 paw_proj%rzetprj(1:maxnprj, 0:maxl) = 0.0_dp
331 ALLOCATE (paw_proj%isoprj(maxnprj, 0:maxl))
332 paw_proj%isoprj = .false.
333 ALLOCATE (gccprj(maxnprj, maxpgf, 0:maxl, nset))
334 gccprj = 0.0_dp
335 ALLOCATE (gcchprj(maxnprj, maxpgf, 0:maxl, nset))
336 gcchprj = 0.0_dp
337
338 NULLIFY (zet, zetp, gcc, smat, work)
339 ! Generate the projetor basis for each ang. mom. q.n.
340 DO lshell = 0, maxl ! lshell
341
342 np = paw_proj%nprj(lshell)
343
344 ALLOCATE (isoprj(np))
345 isoprj = .false.
346
347 ALLOCATE (zet(np), zetp(np), gcc(np, np), gcch(np, np), smat(np, np), work(np, np))
348
349 zet(:) = 0.0_dp
350 zetp(:) = 0.0_dp
351 gcc(:, :) = 0.0_dp
352 gcch(:, :) = 0.0_dp
353 smat(:, :) = 0.0_dp
354 work(:, :) = 0.0_dp
355
356 npgfg = 0
357 ! Collect all the exponent which contribute to lshell
358 DO iset = 1, nset ! iset
359 IF (lshell >= lmin(iset) .AND. lshell <= lmax(iset)) THEN
360 DO ip = 1, npgf(iset)
361 IF (set_radius2(ip, lshell, iset) < max_rad_local) THEN
362 npgfg = npgfg + 1
363 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 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 nisop = 0
376 DO ip = 1, np
377 ! Check for equal exponents
378 DO ipp = 1, ip - 1
379 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 " for the same l quantum numbers")
384 END IF
385 END DO
386
387 IF (zet(ip) >= paw_proj%zisomin(lshell)) THEN
388 isoprj(ip) = .true.
389 nisop = nisop + 1
390 ELSE
391 isoprj(ip) = .false.
392 END IF
393 END DO
394
395 ! Smallest exp. due to eps_fit: where to start geometric progression
396 zetmin = gauss_exponent(lshell, rc, eps_fit, 1.0_dp)
397
398 ! Generate the projectors by the geometric progression
399 IF (np - nisop - 1 > 2) THEN
400 x = (80.0_dp/zetmin)**(1.0_dp/real(np - nisop - 1, dp))
401 ELSE
402 x = 2.0_dp
403 END IF
404 IF (x > 2.0_dp) x = 2.0_dp
405
406 zetval = zetmin
407 DO ip = np, 1, -1
408 IF (.NOT. isoprj(ip)) THEN
409 zetp(ip) = zetval
410 zetval = x*zetval
411 END IF
412 END DO
413
414 nisop = 0
415 DO ip = np, 1, -1
416 IF (isoprj(ip)) THEN
417 zetp(ip) = zetval
418 zetval = x*zetval
419 nisop = nisop + 1
420 END IF
421 END DO
422
423 ! *** Build the overlap matrix: <projector|primitive> ***
424 prefac = 0.5_dp**(lshell + 2)*rootpi*dfac(2*lshell + 1)
425 expzet = real(lshell, dp) + 1.5_dp
426
427 DO ip = 1, np
428 IF (isoprj(ip)) THEN
429 DO jp = 1, np
430 IF (isoprj(jp)) THEN
431 smat(ip, jp) = prefac/(zetp(ip) + zet(jp))**expzet
432 END IF
433 END DO
434 ELSE
435 DO jp = 1, np
436 IF (.NOT. isoprj(jp)) THEN
437 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 IF (eps_svd .EQ. 0.0_dp) THEN
445 CALL invert_matrix(smat, gcc, my_error, "T")
446 ELSE
447 work = transpose(smat)
448 ! Workspace query
449 ALLOCATE (iwork(8*np), s(np), u(np, np), vt(np, np), work_dgesdd(1))
450 lwork = -1
451 CALL dgesdd('S', np, np, work, np, s, u, np, vt, np, work_dgesdd, lwork, iwork, info)
452 lwork = int(work_dgesdd(1))
453 DEALLOCATE (work_dgesdd); ALLOCATE (work_dgesdd(lwork))
454 CALL dgesdd('S', np, np, work, np, s, u, np, vt, np, work_dgesdd, lwork, iwork, info)
455 ! Construct the inverse
456 DO k = 1, np
457 ! invert SV
458 IF (s(k) < eps_svd) THEN
459 s(k) = 0.0_dp
460 ELSE
461 s(k) = 1.0_dp/s(k)
462 END IF
463 vt(k, :) = vt(k, :)*s(k)
464 END DO
465 CALL dgemm('T', 'T', np, np, np, 1.0_dp, vt, np, u, np, 0.0_dp, gcc, np)
466 DEALLOCATE (iwork, s, u, vt, work_dgesdd)
467 END IF
468
469 ! Set the coefficient of the isolated projectors to 0
470 gcch(:, :) = gcc(:, :)
471 DO ip = 1, np
472 IF (isoprj(ip)) THEN
473 gcc(:, ip) = 0.0_dp
474 gcc(ip, :) = 0.0_dp
475 END IF
476 END DO
477
478 ! Transfer data from local to global variables
479
480 paw_proj%zetprj(1:np, lshell) = zetp(1:np)
481 paw_proj%isoprj(1:np, lshell) = isoprj(1:np)
482
483 npgfg = 0
484 DO iset = 1, nset ! iset
485 IF (lshell >= lmin(iset) .AND. lshell <= lmax(iset)) THEN
486 DO ip = 1, npgf(iset)
487 IF (set_radius2(ip, lshell, iset) < max_rad_local) THEN
488 npgfg = npgfg + 1
489 gccprj(1:np, ip, lshell, iset) = gcc(1:np, npgfg)
490 gcchprj(1:np, ip, lshell, iset) = gcch(1:np, npgfg)
491 ELSE
492 gccprj(1:np, ip, lshell, iset) = 0.0_dp
493 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 IF (output_unit > 0) THEN
501 WRITE (unit=output_unit, fmt="(/,/,T2,A,I2)") &
502 "Built projector for angular momentum quantum number l= ", lshell
503 WRITE (unit=output_unit, fmt="(T2,A,I2)") &
504 "Number of isolated projectors = ", nisop
505 DO iset = 1, nset ! iset
506 IF (lshell >= lmin(iset) .AND. lshell <= lmax(iset)) THEN
507 WRITE (unit=output_unit, fmt="(/,T2,A,I5,/,/,T4,A9,(T13,4f15.6))") &
508 "Set ", iset, "exp prj: ", &
509 (paw_proj%zetprj(ip, lshell), ip=1, np)
510 DO jp = 1, npgf(iset)
511 WRITE (unit=output_unit, fmt="(/,T4,A9,F15.6,/,T4,A9,(t13,4E15.6))") &
512 "exp gto: ", zetb(jp, iset), &
513 "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 DEALLOCATE (isoprj)
521 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 "DFT%PRINT%GAPW%PROJECTORS")
526
527 ! Release the working storage for the current value lshell
528 DEALLOCATE (set_radius)
529 DEALLOCATE (set_radius2)
530
531 ! Count primitives basis functions for the projectors
532 paw_proj%ncgauprj = 0
533 paw_proj%nsgauprj = 0
534 DO lshell = 0, maxl
535 paw_proj%ncgauprj = paw_proj%ncgauprj + nco(lshell)*paw_proj%nprj(lshell)
536 paw_proj%nsgauprj = paw_proj%nsgauprj + nso(lshell)*paw_proj%nprj(lshell)
537 END DO
538
539 ncgauprj = paw_proj%ncgauprj
540 nsgauprj = paw_proj%nsgauprj
541 CALL reallocate(paw_proj%lx, 1, ncgauprj)
542 CALL reallocate(paw_proj%ly, 1, ncgauprj)
543 CALL reallocate(paw_proj%lz, 1, ncgauprj)
544 CALL reallocate(paw_proj%first_prj, 0, maxl)
545 CALL reallocate(paw_proj%last_prj, 0, maxl)
546 CALL reallocate(paw_proj%ll, 1, nsgauprj)
547 CALL reallocate(paw_proj%m, 1, nsgauprj)
548 CALL reallocate(paw_proj%first_prjs, 0, maxl)
549
550 ALLOCATE (cprj_s(1:nsgauprj, 1:maxso*nset))
551 ALLOCATE (cprj_h(1:nsgauprj, 1:maxso*nset))
552 cprj_s = 0.0_dp
553 cprj_h = 0.0_dp
554
555 ncgauprj = 0
556 nsgauprj = 0
557 DO lshell = 0, maxl
558 np = paw_proj%nprj(lshell)
559 paw_proj%first_prj(lshell) = ncgauprj + 1
560 paw_proj%first_prjs(lshell) = nsgauprj + 1
561 paw_proj%last_prj(lshell) = ncgauprj + nco(lshell)*np
562 DO ip = 1, np
563 DO ico = ncoset(lshell - 1) + 1, ncoset(lshell)
564 ncgauprj = ncgauprj + 1
565 paw_proj%lx(ncgauprj) = indco(1, ico)
566 paw_proj%ly(ncgauprj) = indco(2, ico)
567 paw_proj%lz(ncgauprj) = indco(3, ico)
568 END DO ! ico
569 DO iso = nsoset(lshell - 1) + 1, nsoset(lshell)
570 nsgauprj = nsgauprj + 1
571 paw_proj%ll(nsgauprj) = indso(1, iso)
572 paw_proj%m(nsgauprj) = indso(2, iso)
573 END DO
574 END DO ! ip
575 END DO ! lshell
576
577 ms = 0
578 DO iset = 1, nset
579 ns = nsoset(lmax(iset))
580 DO lshell = lmin(iset), lmax(iset)
581 iprjfirst = paw_proj%first_prjs(lshell)
582 np = paw_proj%nprj(lshell)
583 DO ipgf = 1, npgf(iset)
584 DO ip = 1, np
585 DO il = 1, nso(lshell)
586 iprjs = iprjfirst - 1 + il + (ip - 1)*nso(lshell)
587 iso = nsoset(lshell - 1) + 1 + (lshell + paw_proj%m(iprjs))
588
589 iso = iso + (ipgf - 1)*ns + ms
590 cprj_s(iprjs, iso) = gccprj(ip, ipgf, lshell, iset)
591 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 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 nset=nseta, ncgf=ncgfo, nsgf=nsgfo)
603
604 ALLOCATE (paw_proj%local_oce_sphi_h(maxco, nsgfo))
605 paw_proj%local_oce_sphi_h = 0.0_dp
606 ALLOCATE (paw_proj%sphi_h(maxco, nsgfo))
607 paw_proj%sphi_h = 0.0_dp
608
609 ALLOCATE (paw_proj%local_oce_sphi_s(maxco, nsgfo))
610 paw_proj%local_oce_sphi_s = 0.0_dp
611 ALLOCATE (paw_proj%sphi_s(maxco, nsgfo))
612 paw_proj%sphi_s = 0.0_dp
613
614 ! only use first nset of orb basis local projection!
615 DO iset = 1, nseta
616 n = ncoset(lmax(iset))
617 DO ipgf = 1, npgf(iset)
618 DO ishell = 1, nshell(iset)
619 lshell = l(ishell, iset)
620 icomin = ncoset(lshell - 1) + 1 + n*(ipgf - 1)
621 icomax = ncoset(lshell) + n*(ipgf - 1)
622 isgfmin = first_sgf(ishell, iset)
623 isgfmax = last_sgf(ishell, iset)
624 radius = exp_radius(lshell, basis_1c%zet(ipgf, iset), &
625 eps_fit, 1.0_dp)
626 DO isgf = isgfmin, isgfmax
627 paw_proj%sphi_h(icomin:icomax, isgf) = &
628 sphi(icomin:icomax, isgf)
629 IF (radius < rc) THEN
630 paw_proj%sphi_s(icomin:icomax, isgf) = 0.0_dp
631 ELSE
632 paw_proj%sphi_s(icomin:icomax, isgf) = &
633 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 DO iset = 1, nseta
642 n = ncoset(lmax(iset))
643 ns = nsoset(lmax(iset))
644 DO ipgf = 1, npgf(iset)
645 DO ishell = 1, nshell(iset)
646 lshell = l(ishell, iset)
647 icomin = ncoset(lshell - 1) + 1 + n*(ipgf - 1)
648 icomax = ncoset(lshell) + n*(ipgf - 1)
649 isgfmin = first_sgf(ishell, iset)
650 isgfmax = last_sgf(ishell, iset)
651 isomin = nsoset(lshell - 1) + 1 + ns*(ipgf - 1)
652 DO is = 1, nso(lshell)
653 iso = isomin + is - 1
654 DO ic = 1, nco(lshell)
655 ico = icomin + ic - 1
656 lx = indco(1, ic + ncoset(lshell - 1))
657 ly = indco(2, ic + ncoset(lshell - 1))
658 lz = indco(3, ic + ncoset(lshell - 1))
659 DO isgf = isgfmin, isgfmax
660 paw_proj%local_oce_sphi_h(iso, isgf) = &
661 paw_proj%local_oce_sphi_h(iso, isgf) + &
662 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 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 ALLOCATE (paw_proj%o2nindex(maxso*nset))
675 ALLOCATE (paw_proj%n2oindex(maxso*nset))
676 paw_proj%o2nindex = 0
677 paw_proj%n2oindex = 0
678 ico = 1
679 DO iset = 1, nset
680 iso_set = (iset - 1)*maxso + 1
681 nsox = nsoset(lmax(iset))
682 DO ipgf = 1, npgf(iset)
683 iso_pgf = iso_set + (ipgf - 1)*nsox
684 iso = iso_pgf + nsoset(lmin(iset) - 1)
685 DO lx = lmin(iset), lmax(iset)
686 DO k = 1, nso(lx)
687 paw_proj%n2oindex(ico) = iso
688 paw_proj%o2nindex(iso) = ico
689 iso = iso + 1
690 ico = ico + 1
691 END DO
692 END DO
693 END DO
694 END DO
695 mp = ico - 1
696 paw_proj%nsatbas = mp
697 paw_proj%nsotot = nset*maxso
698 ALLOCATE (paw_proj%csprj(nsgauprj, mp))
699 paw_proj%csprj = 0.0_dp
700 DO k = 1, mp
701 ico = paw_proj%n2oindex(k)
702 paw_proj%csprj(:, k) = cprj_s(:, ico)
703 END DO
704 ALLOCATE (paw_proj%chprj(nsgauprj, mp))
705 paw_proj%chprj = 0.0_dp
706 DO k = 1, mp
707 ico = paw_proj%n2oindex(k)
708 paw_proj%chprj(:, k) = cprj_h(:, ico)
709 END DO
710 DEALLOCATE (cprj_s, cprj_h, gcchprj, gccprj)
711
712 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 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 IF (ASSOCIATED(paw_proj_set)) THEN
757 IF (PRESENT(csprj)) csprj => paw_proj_set%csprj
758 IF (PRESENT(chprj)) chprj => paw_proj_set%chprj
759 IF (PRESENT(local_oce_sphi_h)) local_oce_sphi_h => paw_proj_set%local_oce_sphi_h
760 IF (PRESENT(local_oce_sphi_s)) local_oce_sphi_s => paw_proj_set%local_oce_sphi_s
761 IF (PRESENT(first_prj)) first_prj => paw_proj_set%first_prj
762 IF (PRESENT(last_prj)) last_prj => paw_proj_set%last_prj
763 IF (PRESENT(first_prjs)) first_prjs => paw_proj_set%first_prjs
764 IF (PRESENT(maxl)) maxl = paw_proj_set%maxl
765 IF (PRESENT(ncgauprj)) ncgauprj = paw_proj_set%ncgauprj
766 IF (PRESENT(nsgauprj)) nsgauprj = paw_proj_set%nsgauprj
767 IF (PRESENT(nsatbas)) nsatbas = paw_proj_set%nsatbas
768 IF (PRESENT(nsotot)) nsotot = paw_proj_set%nsotot
769 IF (PRESENT(nprj)) nprj => paw_proj_set%nprj
770 IF (PRESENT(rcprj)) rcprj = paw_proj_set%rcprj
771 IF (PRESENT(rzetprj)) rzetprj => paw_proj_set%rzetprj
772 IF (PRESENT(zisomin)) zisomin => paw_proj_set%zisomin
773 IF (PRESENT(zetprj)) zetprj => paw_proj_set%zetprj
774 IF (PRESENT(o2nindex)) o2nindex => paw_proj_set%o2nindex
775 IF (PRESENT(n2oindex)) n2oindex => paw_proj_set%n2oindex
776 ELSE
777 cpabort("The pointer paw_proj_set is not associated")
778 END IF
779
780 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 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 IF (ASSOCIATED(paw_proj_set)) THEN
796 IF (PRESENT(rzetprj)) paw_proj_set%rzetprj(:, 0:) = rzetprj(:, 0:)
797 IF (PRESENT(rcprj)) paw_proj_set%rcprj = rcprj
798 ELSE
799 cpabort("The pointer paw_proj_set is not associated")
800 END IF
801
802 END SUBROUTINE set_paw_proj_set
803
804END MODULE paw_proj_set_types
dgemm
static void dgemm(const char transa, const char transb, const int m, const int n, const int k, const double alpha, const double *a, const int lda, const double *b, const int ldb, const double beta, double *c, const int ldc)
Convenient wrapper to hide Fortran nature of dgemm_, swapping a and b.
Definition grid_cpu_task_list.c:206
mathlib::invert_matrix
Definition mathlib.F:70
memory_utilities::reallocate
Definition memory_utilities.F:27
ao_util
All kind of helpful little routines.
Definition ao_util.F:14
ao_util::exp_radius
real(kind=dp) function, public exp_radius(l, alpha, threshold, prefactor, epsabs, epsrel, rlow)
The radius of a primitive Gaussian function for a given threshold is calculated. g(r) = prefactor*r**...
Definition ao_util.F:96
ao_util::gauss_exponent
real(kind=dp) function, public gauss_exponent(l, radius, threshold, prefactor)
The exponent of a primitive Gaussian function for a given radius and threshold is calculated.
Definition ao_util.F:60
basis_set_types
Definition basis_set_types.F:15
basis_set_types::get_gto_basis_set
subroutine, public get_gto_basis_set(gto_basis_set, name, aliases, norm_type, kind_radius, ncgf, nset, nsgf, cgf_symbol, sgf_symbol, norm_cgf, set_radius, lmax, lmin, lx, ly, lz, m, ncgf_set, npgf, nsgf_set, nshell, cphi, pgf_radius, sphi, scon, zet, first_cgf, first_sgf, l, last_cgf, last_sgf, n, gcc, maxco, maxl, maxpgf, maxsgf_set, maxshell, maxso, nco_sum, npgf_sum, nshell_sum, maxder, short_kind_radius, npgf_seg_sum)
...
Definition basis_set_types.F:624
cp_control_types
Defines control structures, which contain the parameters and the settings for the DFT-based calculati...
Definition cp_control_types.F:12
cp_log_handling
various routines to log and control the output. The idea is that decisions about where to log should ...
Definition cp_log_handling.F:41
cp_log_handling::cp_get_default_logger
type(cp_logger_type) function, pointer, public cp_get_default_logger()
returns the default logger
Definition cp_log_handling.F:234
cp_output_handling
routines to handle the output, The idea is to remove the decision of wheter to output and what to out...
Definition cp_output_handling.F:25
cp_output_handling::cp_print_key_unit_nr
integer function, public cp_print_key_unit_nr(logger, basis_section, print_key_path, extension, middle_name, local, log_filename, ignore_should_output, file_form, file_position, file_action, file_status, do_backup, on_file, is_new_file, mpi_io, fout)
...
Definition cp_output_handling.F:853
cp_output_handling::cp_print_key_finished_output
subroutine, public cp_print_key_finished_output(unit_nr, logger, basis_section, print_key_path, local, ignore_should_output, on_file, mpi_io)
should be called after you finish working with a unit obtained with cp_print_key_unit_nr,...
Definition cp_output_handling.F:1063
input_section_types
objects that represent the structure of input sections and the data contained in an input section
Definition input_section_types.F:15
kinds
Defines the basic variable types.
Definition kinds.F:23
kinds::dp
integer, parameter, public dp
Definition kinds.F:34
kinds::default_string_length
integer, parameter, public default_string_length
Definition kinds.F:57
mathconstants
Definition of mathematical constants and functions.
Definition mathconstants.F:16
mathconstants::dfac
real(kind=dp), dimension(-1:2 *maxfac+1), parameter, public dfac
Definition mathconstants.F:61
mathconstants::rootpi
real(kind=dp), parameter, public rootpi
Definition mathconstants.F:114
mathlib
Collection of simple mathematical functions and subroutines.
Definition mathlib.F:15
memory_utilities
Utility routines for the memory handling.
Definition memory_utilities.F:14
orbital_pointers
Provides Cartesian and spherical orbital pointers and indices.
Definition orbital_pointers.F:15
orbital_pointers::nco
integer, dimension(:), allocatable, public nco
Definition orbital_pointers.F:42
orbital_pointers::nsoset
integer, dimension(:), allocatable, public nsoset
Definition orbital_pointers.F:42
orbital_pointers::indso
integer, dimension(:, :), allocatable, public indso
Definition orbital_pointers.F:43
orbital_pointers::ncoset
integer, dimension(:), allocatable, public ncoset
Definition orbital_pointers.F:42
orbital_pointers::indco
integer, dimension(:, :), allocatable, public indco
Definition orbital_pointers.F:43
orbital_pointers::nso
integer, dimension(:), allocatable, public nso
Definition orbital_pointers.F:42
orbital_transformation_matrices
Calculation of the spherical harmonics and the corresponding orbital transformation matrices.
Definition orbital_transformation_matrices.F:17
orbital_transformation_matrices::orbtramat
type(orbtramat_type), dimension(:), pointer, public orbtramat
Definition orbital_transformation_matrices.F:47
paw_proj_set_types
Definition paw_proj_set_types.F:13
paw_proj_set_types::get_paw_proj_set
subroutine, public get_paw_proj_set(paw_proj_set, csprj, chprj, first_prj, first_prjs, last_prj, local_oce_sphi_h, local_oce_sphi_s, maxl, ncgauprj, nsgauprj, nsatbas, nsotot, nprj, o2nindex, n2oindex, rcprj, rzetprj, zisomin, zetprj)
Get informations about a paw projectors set.
Definition paw_proj_set_types.F:744
paw_proj_set_types::allocate_paw_proj_set
subroutine, public allocate_paw_proj_set(paw_proj_set)
Allocate projector type for GAPW.
Definition paw_proj_set_types.F:95
paw_proj_set_types::projectors
subroutine, public projectors(paw_proj, basis_1c, orb_basis, rc, qs_control, max_rad_local_type, force_env_section)
Initialize the projector-type set data set.
Definition paw_proj_set_types.F:180
paw_proj_set_types::deallocate_paw_proj_set
subroutine, public deallocate_paw_proj_set(paw_proj_set)
Deallocate a projector-type set data set.
Definition paw_proj_set_types.F:135
paw_proj_set_types::set_paw_proj_set
subroutine, public set_paw_proj_set(paw_proj_set, rzetprj, rcprj)
Set informations about a paw projectors set.
Definition paw_proj_set_types.F:790
basis_set_types::gto_basis_set_type
Definition basis_set_types.F:72
cp_control_types::qs_control_type
Definition cp_control_types.F:358
cp_log_handling::cp_logger_type
type of a logger, at the moment it contains just a print level starting at which level it should be l...
Definition cp_log_handling.F:140
input_section_types::section_vals_type
stores the values of a section
Definition input_section_types.F:127
paw_proj_set_types::paw_proj_set_type
Definition paw_proj_set_types.F:52