(git:b279b6b)
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-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  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 
804 END 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:195
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)
...
Definition: basis_set_types.F:615
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:803
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:1013
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