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 1614 : SUBROUTINE allocate_paw_proj_set(paw_proj_set)
95 :
96 : TYPE(paw_proj_set_type), POINTER :: paw_proj_set
97 :
98 1614 : IF (ASSOCIATED(paw_proj_set)) CALL deallocate_paw_proj_set(paw_proj_set)
99 :
100 1614 : 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 1614 : 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 1614 : SUBROUTINE deallocate_paw_proj_set(paw_proj_set)
135 : TYPE(paw_proj_set_type), POINTER :: paw_proj_set
136 :
137 1614 : IF (ASSOCIATED(paw_proj_set)) THEN
138 :
139 1614 : IF (ASSOCIATED(paw_proj_set%zisomin)) DEALLOCATE (paw_proj_set%zisomin)
140 1614 : IF (ASSOCIATED(paw_proj_set%nprj)) DEALLOCATE (paw_proj_set%nprj)
141 1614 : IF (ASSOCIATED(paw_proj_set%lx)) DEALLOCATE (paw_proj_set%lx)
142 1614 : IF (ASSOCIATED(paw_proj_set%ly)) DEALLOCATE (paw_proj_set%ly)
143 1614 : IF (ASSOCIATED(paw_proj_set%lz)) DEALLOCATE (paw_proj_set%lz)
144 1614 : IF (ASSOCIATED(paw_proj_set%ll)) DEALLOCATE (paw_proj_set%ll)
145 1614 : IF (ASSOCIATED(paw_proj_set%m)) DEALLOCATE (paw_proj_set%m)
146 1614 : IF (ASSOCIATED(paw_proj_set%first_prj)) DEALLOCATE (paw_proj_set%first_prj)
147 1614 : IF (ASSOCIATED(paw_proj_set%last_prj)) DEALLOCATE (paw_proj_set%last_prj)
148 1614 : IF (ASSOCIATED(paw_proj_set%first_prjs)) DEALLOCATE (paw_proj_set%first_prjs)
149 1614 : IF (ASSOCIATED(paw_proj_set%zetprj)) DEALLOCATE (paw_proj_set%zetprj)
150 1614 : IF (ASSOCIATED(paw_proj_set%csprj)) DEALLOCATE (paw_proj_set%csprj)
151 1614 : IF (ASSOCIATED(paw_proj_set%chprj)) DEALLOCATE (paw_proj_set%chprj)
152 1614 : IF (ASSOCIATED(paw_proj_set%local_oce_sphi_h)) DEALLOCATE (paw_proj_set%local_oce_sphi_h)
153 1614 : IF (ASSOCIATED(paw_proj_set%local_oce_sphi_s)) DEALLOCATE (paw_proj_set%local_oce_sphi_s)
154 1614 : IF (ASSOCIATED(paw_proj_set%sphi_h)) DEALLOCATE (paw_proj_set%sphi_h)
155 1614 : IF (ASSOCIATED(paw_proj_set%sphi_s)) DEALLOCATE (paw_proj_set%sphi_s)
156 1614 : IF (ASSOCIATED(paw_proj_set%isoprj)) DEALLOCATE (paw_proj_set%isoprj)
157 1614 : IF (ASSOCIATED(paw_proj_set%rzetprj)) DEALLOCATE (paw_proj_set%rzetprj)
158 1614 : IF (ASSOCIATED(paw_proj_set%o2nindex)) DEALLOCATE (paw_proj_set%o2nindex)
159 1614 : IF (ASSOCIATED(paw_proj_set%n2oindex)) DEALLOCATE (paw_proj_set%n2oindex)
160 :
161 1614 : DEALLOCATE (paw_proj_set)
162 :
163 : END IF
164 :
165 1614 : 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 1614 : 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 1614 : eps_fit = qs_control%gapw_control%eps_fit
192 1614 : eps_iso = qs_control%gapw_control%eps_iso
193 1614 : eps_svd = qs_control%gapw_control%eps_svd
194 1614 : max_rad_local = qs_control%gapw_control%max_rad_local
195 1614 : IF (max_rad_local_type .LT. max_rad_local) THEN
196 1614 : max_rad_local = max_rad_local_type
197 : END IF
198 1614 : 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 1614 : rc, eps_orb, max_rad_local, force_env_section)
202 :
203 1614 : 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 1614 : 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 1614 : INTEGER, ALLOCATABLE, DIMENSION(:) :: IWORK
234 1614 : INTEGER, DIMENSION(:), POINTER :: lmax, lmin, npgf, nshell
235 1614 : INTEGER, DIMENSION(:, :), POINTER :: first_cgf, first_sgf, l, last_cgf, &
236 1614 : last_sgf
237 1614 : LOGICAL, ALLOCATABLE, DIMENSION(:) :: isoprj
238 : REAL(KIND=dp) :: expzet, my_error, prefac, radius, x, &
239 : zetmin, zetval
240 1614 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: S, Work_dgesdd
241 1614 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: U, VT
242 1614 : REAL(KIND=dp), DIMENSION(:), POINTER :: set_radius, zet, zetp
243 1614 : REAL(KIND=dp), DIMENSION(:, :), POINTER :: cprj_h, cprj_s, gcc, gcch, smat, sphi, &
244 1614 : work, zetb
245 1614 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: gcca, set_radius2
246 1614 : REAL(KIND=dp), DIMENSION(:, :, :, :), POINTER :: gcchprj, gccprj
247 : TYPE(cp_logger_type), POINTER :: logger
248 :
249 1614 : NULLIFY (logger)
250 3228 : logger => cp_get_default_logger()
251 :
252 1614 : NULLIFY (first_cgf, first_sgf, last_cgf, last_sgf, gcc, l, set_radius, set_radius2)
253 1614 : NULLIFY (sphi, lmax, lmin, npgf, nshell, zet, zetb, zetp, smat, work, gcca)
254 :
255 1614 : CPASSERT(ASSOCIATED(paw_proj))
256 1614 : CPASSERT(ASSOCIATED(orb_basis))
257 1614 : 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 1614 : zet=zetb, gcc=gcca)
267 :
268 1614 : paw_proj%maxl = maxl
269 1614 : CPASSERT(.NOT. ASSOCIATED(paw_proj%zisomin))
270 1614 : CPASSERT(.NOT. ASSOCIATED(paw_proj%nprj))
271 :
272 4842 : ALLOCATE (paw_proj%zisomin(0:maxl))
273 5466 : paw_proj%zisomin(0:maxl) = 0.0_dp
274 4842 : ALLOCATE (paw_proj%nprj(0:maxl))
275 5466 : paw_proj%nprj(0:maxl) = 0
276 :
277 : output_unit = cp_print_key_unit_nr(logger, force_env_section, &
278 1614 : "DFT%PRINT%GAPW%PROJECTORS", extension=".Log")
279 :
280 1614 : 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 4842 : ALLOCATE (set_radius(nset))
286 6760 : set_radius = 0.0_dp
287 6760 : DO iset = 1, nset
288 15048 : 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 13434 : rlow=set_radius(iset)))
293 : END DO ! is
294 : END DO ! iset
295 :
296 8070 : ALLOCATE (set_radius2(maxpgf, 0:maxl, nset))
297 89548 : set_radius2 = 0.0_dp
298 6760 : DO iset = 1, nset
299 13264 : DO lshell = lmin(iset), lmax(iset)
300 28116 : DO ip = 1, npgf(iset)
301 : set_radius2(ip, lshell, iset) = &
302 22970 : exp_radius(lshell, zetb(ip, iset), eps_orb, 1.0_dp)
303 : END DO
304 : END DO ! is
305 : END DO ! iset
306 :
307 1614 : maxnprj = 0
308 5466 : DO lshell = 0, maxl ! lshell
309 3852 : np = 0
310 17204 : DO iset = 1, nset
311 17204 : IF (lshell >= lmin(iset) .AND. lshell <= lmax(iset)) THEN
312 22970 : DO ip = 1, npgf(iset)
313 22970 : IF (set_radius2(ip, lshell, iset) < max_rad_local) THEN
314 16332 : np = np + 1
315 : END IF
316 : END DO
317 : END IF
318 : END DO
319 3852 : maxnprj = MAX(maxnprj, np)
320 3852 : paw_proj%nprj(lshell) = np
321 5466 : 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 6456 : ALLOCATE (paw_proj%zetprj(maxnprj, 0:maxl))
328 31018 : paw_proj%zetprj(1:maxnprj, 0:maxl) = 0.0_dp
329 6456 : ALLOCATE (paw_proj%rzetprj(maxnprj, 0:maxl))
330 31018 : paw_proj%rzetprj(1:maxnprj, 0:maxl) = 0.0_dp
331 6456 : ALLOCATE (paw_proj%isoprj(maxnprj, 0:maxl))
332 31018 : paw_proj%isoprj = .FALSE.
333 9684 : ALLOCATE (gccprj(maxnprj, maxpgf, 0:maxl, nset))
334 719338 : gccprj = 0.0_dp
335 9684 : ALLOCATE (gcchprj(maxnprj, maxpgf, 0:maxl, nset))
336 719338 : gcchprj = 0.0_dp
337 :
338 1614 : NULLIFY (zet, zetp, gcc, smat, work)
339 : ! Generate the projetor basis for each ang. mom. q.n.
340 5466 : DO lshell = 0, maxl ! lshell
341 :
342 3852 : np = paw_proj%nprj(lshell)
343 :
344 11556 : ALLOCATE (isoprj(np))
345 20184 : isoprj = .FALSE.
346 :
347 65484 : ALLOCATE (zet(np), zetp(np), gcc(np, np), gcch(np, np), smat(np, np), work(np, np))
348 :
349 20184 : zet(:) = 0.0_dp
350 20184 : zetp(:) = 0.0_dp
351 123784 : gcc(:, :) = 0.0_dp
352 123784 : gcch(:, :) = 0.0_dp
353 123784 : smat(:, :) = 0.0_dp
354 123784 : work(:, :) = 0.0_dp
355 :
356 3852 : npgfg = 0
357 : ! Collect all the exponent which contribute to lshell
358 17204 : DO iset = 1, nset ! iset
359 17204 : IF (lshell >= lmin(iset) .AND. lshell <= lmax(iset)) THEN
360 22970 : DO ip = 1, npgf(iset)
361 22970 : IF (set_radius2(ip, lshell, iset) < max_rad_local) THEN
362 16332 : npgfg = npgfg + 1
363 16332 : 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 3852 : 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 3852 : nisop = 0
376 20184 : DO ip = 1, np
377 : ! Check for equal exponents
378 59966 : DO ipp = 1, ip - 1
379 59966 : 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 20184 : IF (zet(ip) >= paw_proj%zisomin(lshell)) THEN
388 3438 : isoprj(ip) = .TRUE.
389 3438 : nisop = nisop + 1
390 : ELSE
391 12894 : isoprj(ip) = .FALSE.
392 : END IF
393 : END DO
394 :
395 : ! Smallest exp. due to eps_fit: where to start geometric progression
396 3852 : zetmin = gauss_exponent(lshell, rc, eps_fit, 1.0_dp)
397 :
398 : ! Generate the projectors by the geometric progression
399 3852 : IF (np - nisop - 1 > 2) THEN
400 2200 : x = (80.0_dp/zetmin)**(1.0_dp/REAL(np - nisop - 1, dp))
401 : ELSE
402 : x = 2.0_dp
403 : END IF
404 3852 : IF (x > 2.0_dp) x = 2.0_dp
405 :
406 3852 : zetval = zetmin
407 20184 : DO ip = np, 1, -1
408 20184 : IF (.NOT. isoprj(ip)) THEN
409 12894 : zetp(ip) = zetval
410 12894 : zetval = x*zetval
411 : END IF
412 : END DO
413 :
414 3852 : nisop = 0
415 20184 : DO ip = np, 1, -1
416 20184 : IF (isoprj(ip)) THEN
417 3438 : zetp(ip) = zetval
418 3438 : zetval = x*zetval
419 3438 : nisop = nisop + 1
420 : END IF
421 : END DO
422 :
423 : ! *** Build the overlap matrix: <projector|primitive> ***
424 3852 : prefac = 0.5_dp**(lshell + 2)*rootpi*dfac(2*lshell + 1)
425 3852 : expzet = REAL(lshell, dp) + 1.5_dp
426 :
427 20184 : DO ip = 1, np
428 20184 : IF (isoprj(ip)) THEN
429 36946 : DO jp = 1, np
430 36946 : IF (isoprj(jp)) THEN
431 17806 : smat(ip, jp) = prefac/(zetp(ip) + zet(jp))**expzet
432 : END IF
433 : END DO
434 : ELSE
435 82986 : DO jp = 1, np
436 82986 : IF (.NOT. isoprj(jp)) THEN
437 54390 : 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 3852 : IF (eps_svd .EQ. 0.0_dp) THEN
445 222 : CALL invert_matrix(smat, gcc, my_error, "T")
446 : ELSE
447 230358 : work = TRANSPOSE(smat)
448 : ! Workspace query
449 39930 : ALLOCATE (iwork(8*np), S(np), U(np, np), VT(np, np), work_dgesdd(1))
450 3630 : lwork = -1
451 3630 : CALL DGESDD('S', np, np, work, np, S, U, np, vt, np, work_dgesdd, lwork, iwork, info)
452 3630 : lwork = INT(work_dgesdd(1))
453 10890 : DEALLOCATE (work_dgesdd); ALLOCATE (work_dgesdd(lwork))
454 3630 : CALL DGESDD('S', np, np, work, np, S, U, np, vt, np, work_dgesdd, lwork, iwork, info)
455 : ! Construct the inverse
456 19062 : DO k = 1, np
457 : ! invert SV
458 15432 : IF (S(k) < eps_svd) THEN
459 1072 : S(k) = 0.0_dp
460 : ELSE
461 14360 : S(k) = 1.0_dp/S(k)
462 : END IF
463 116994 : VT(k, :) = VT(k, :)*S(k)
464 : END DO
465 3630 : CALL DGEMM('T', 'T', np, np, np, 1.0_dp, VT, np, U, np, 0.0_dp, gcc, np)
466 3630 : DEALLOCATE (iwork, S, U, VT, work_dgesdd)
467 : END IF
468 :
469 : ! Set the coefficient of the isolated projectors to 0
470 243716 : gcch(:, :) = gcc(:, :)
471 20184 : DO ip = 1, np
472 20184 : IF (isoprj(ip)) THEN
473 36946 : gcc(:, ip) = 0.0_dp
474 36946 : gcc(ip, :) = 0.0_dp
475 : END IF
476 : END DO
477 :
478 : ! Transfer data from local to global variables
479 :
480 36516 : paw_proj%zetprj(1:np, lshell) = zetp(1:np)
481 20184 : paw_proj%isoprj(1:np, lshell) = isoprj(1:np)
482 :
483 3852 : npgfg = 0
484 17204 : DO iset = 1, nset ! iset
485 17204 : IF (lshell >= lmin(iset) .AND. lshell <= lmax(iset)) THEN
486 22970 : DO ip = 1, npgf(iset)
487 22970 : IF (set_radius2(ip, lshell, iset) < max_rad_local) THEN
488 16332 : npgfg = npgfg + 1
489 223532 : gccprj(1:np, ip, lshell, iset) = gcc(1:np, npgfg)
490 223532 : 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 3852 : 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 3852 : DEALLOCATE (isoprj)
521 5466 : 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 1614 : "DFT%PRINT%GAPW%PROJECTORS")
526 :
527 : ! Release the working storage for the current value lshell
528 1614 : DEALLOCATE (set_radius)
529 1614 : DEALLOCATE (set_radius2)
530 :
531 : ! Count primitives basis functions for the projectors
532 1614 : paw_proj%ncgauprj = 0
533 1614 : paw_proj%nsgauprj = 0
534 5466 : DO lshell = 0, maxl
535 3852 : paw_proj%ncgauprj = paw_proj%ncgauprj + nco(lshell)*paw_proj%nprj(lshell)
536 5466 : paw_proj%nsgauprj = paw_proj%nsgauprj + nso(lshell)*paw_proj%nprj(lshell)
537 : END DO
538 :
539 1614 : ncgauprj = paw_proj%ncgauprj
540 1614 : nsgauprj = paw_proj%nsgauprj
541 1614 : CALL reallocate(paw_proj%lx, 1, ncgauprj)
542 1614 : CALL reallocate(paw_proj%ly, 1, ncgauprj)
543 1614 : CALL reallocate(paw_proj%lz, 1, ncgauprj)
544 1614 : CALL reallocate(paw_proj%first_prj, 0, maxl)
545 1614 : CALL reallocate(paw_proj%last_prj, 0, maxl)
546 1614 : CALL reallocate(paw_proj%ll, 1, nsgauprj)
547 1614 : CALL reallocate(paw_proj%m, 1, nsgauprj)
548 1614 : CALL reallocate(paw_proj%first_prjs, 0, maxl)
549 :
550 6456 : ALLOCATE (cprj_s(1:nsgauprj, 1:maxso*nset))
551 6456 : ALLOCATE (cprj_h(1:nsgauprj, 1:maxso*nset))
552 3003930 : cprj_s = 0.0_dp
553 3003930 : cprj_h = 0.0_dp
554 :
555 1614 : ncgauprj = 0
556 1614 : nsgauprj = 0
557 5466 : DO lshell = 0, maxl
558 3852 : np = paw_proj%nprj(lshell)
559 3852 : paw_proj%first_prj(lshell) = ncgauprj + 1
560 3852 : paw_proj%first_prjs(lshell) = nsgauprj + 1
561 3852 : paw_proj%last_prj(lshell) = ncgauprj + nco(lshell)*np
562 21798 : DO ip = 1, np
563 49694 : DO ico = ncoset(lshell - 1) + 1, ncoset(lshell)
564 33362 : ncgauprj = ncgauprj + 1
565 33362 : paw_proj%lx(ncgauprj) = indco(1, ico)
566 33362 : paw_proj%ly(ncgauprj) = indco(2, ico)
567 49694 : paw_proj%lz(ncgauprj) = indco(3, ico)
568 : END DO ! ico
569 52044 : DO iso = nsoset(lshell - 1) + 1, nsoset(lshell)
570 31860 : nsgauprj = nsgauprj + 1
571 31860 : paw_proj%ll(nsgauprj) = indso(1, iso)
572 48192 : paw_proj%m(nsgauprj) = indso(2, iso)
573 : END DO
574 : END DO ! ip
575 : END DO ! lshell
576 :
577 1614 : ms = 0
578 6760 : DO iset = 1, nset
579 5146 : ns = nsoset(lmax(iset))
580 11650 : DO lshell = lmin(iset), lmax(iset)
581 6504 : iprjfirst = paw_proj%first_prjs(lshell)
582 6504 : np = paw_proj%nprj(lshell)
583 28116 : DO ipgf = 1, npgf(iset)
584 127416 : DO ip = 1, np
585 293046 : DO il = 1, nso(lshell)
586 172134 : iprjs = iprjfirst - 1 + il + (ip - 1)*nso(lshell)
587 172134 : iso = nsoset(lshell - 1) + 1 + (lshell + paw_proj%m(iprjs))
588 :
589 172134 : iso = iso + (ipgf - 1)*ns + ms
590 172134 : cprj_s(iprjs, iso) = gccprj(ip, ipgf, lshell, iset)
591 276580 : 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 6760 : 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 1614 : nset=nseta, ncgf=ncgfo, nsgf=nsgfo)
603 :
604 6456 : ALLOCATE (paw_proj%local_oce_sphi_h(maxco, nsgfo))
605 444752 : paw_proj%local_oce_sphi_h = 0.0_dp
606 6456 : ALLOCATE (paw_proj%sphi_h(maxco, nsgfo))
607 444752 : paw_proj%sphi_h = 0.0_dp
608 :
609 6456 : ALLOCATE (paw_proj%local_oce_sphi_s(maxco, nsgfo))
610 444752 : paw_proj%local_oce_sphi_s = 0.0_dp
611 6456 : ALLOCATE (paw_proj%sphi_s(maxco, nsgfo))
612 444752 : paw_proj%sphi_s = 0.0_dp
613 :
614 : ! only use first nset of orb basis local projection!
615 6684 : DO iset = 1, nseta
616 5070 : n = ncoset(lmax(iset))
617 19046 : DO ipgf = 1, npgf(iset)
618 43384 : DO ishell = 1, nshell(iset)
619 25952 : lshell = l(ishell, iset)
620 25952 : icomin = ncoset(lshell - 1) + 1 + n*(ipgf - 1)
621 25952 : icomax = ncoset(lshell) + n*(ipgf - 1)
622 25952 : isgfmin = first_sgf(ishell, iset)
623 25952 : isgfmax = last_sgf(ishell, iset)
624 : radius = exp_radius(lshell, basis_1c%zet(ipgf, iset), &
625 25952 : eps_fit, 1.0_dp)
626 90294 : DO isgf = isgfmin, isgfmax
627 : paw_proj%sphi_h(icomin:icomax, isgf) = &
628 396996 : sphi(icomin:icomax, isgf)
629 77932 : IF (radius < rc) THEN
630 78668 : paw_proj%sphi_s(icomin:icomax, isgf) = 0.0_dp
631 : ELSE
632 : paw_proj%sphi_s(icomin:icomax, isgf) = &
633 258330 : 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 6684 : DO iset = 1, nseta
642 5070 : n = ncoset(lmax(iset))
643 5070 : ns = nsoset(lmax(iset))
644 19046 : DO ipgf = 1, npgf(iset)
645 43384 : DO ishell = 1, nshell(iset)
646 25952 : lshell = l(ishell, iset)
647 25952 : icomin = ncoset(lshell - 1) + 1 + n*(ipgf - 1)
648 25952 : icomax = ncoset(lshell) + n*(ipgf - 1)
649 25952 : isgfmin = first_sgf(ishell, iset)
650 25952 : isgfmax = last_sgf(ishell, iset)
651 25952 : isomin = nsoset(lshell - 1) + 1 + ns*(ipgf - 1)
652 90294 : DO is = 1, nso(lshell)
653 51980 : iso = isomin + is - 1
654 250440 : DO ic = 1, nco(lshell)
655 172508 : ico = icomin + ic - 1
656 172508 : lx = indco(1, ic + ncoset(lshell - 1))
657 172508 : ly = indco(2, ic + ncoset(lshell - 1))
658 172508 : lz = indco(3, ic + ncoset(lshell - 1))
659 940248 : DO isgf = isgfmin, isgfmax
660 : paw_proj%local_oce_sphi_h(iso, isgf) = &
661 : paw_proj%local_oce_sphi_h(iso, isgf) + &
662 715760 : 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 888268 : 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 4842 : ALLOCATE (paw_proj%o2nindex(maxso*nset))
675 4842 : ALLOCATE (paw_proj%n2oindex(maxso*nset))
676 67590 : paw_proj%o2nindex = 0
677 67590 : paw_proj%n2oindex = 0
678 1614 : ico = 1
679 6760 : DO iset = 1, nset
680 5146 : iso_set = (iset - 1)*maxso + 1
681 5146 : nsox = nsoset(lmax(iset))
682 19396 : DO ipgf = 1, npgf(iset)
683 12636 : iso_pgf = iso_set + (ipgf - 1)*nsox
684 12636 : iso = iso_pgf + nsoset(lmin(iset) - 1)
685 34248 : DO lx = lmin(iset), lmax(iset)
686 61252 : DO k = 1, nso(lx)
687 32150 : paw_proj%n2oindex(ico) = iso
688 32150 : paw_proj%o2nindex(iso) = ico
689 32150 : iso = iso + 1
690 48616 : ico = ico + 1
691 : END DO
692 : END DO
693 : END DO
694 : END DO
695 1614 : mp = ico - 1
696 1614 : paw_proj%nsatbas = mp
697 1614 : paw_proj%nsotot = nset*maxso
698 6456 : ALLOCATE (paw_proj%csprj(nsgauprj, mp))
699 1133032 : paw_proj%csprj = 0.0_dp
700 33764 : DO k = 1, mp
701 32150 : ico = paw_proj%n2oindex(k)
702 2232300 : paw_proj%csprj(:, k) = cprj_s(:, ico)
703 : END DO
704 6456 : ALLOCATE (paw_proj%chprj(nsgauprj, mp))
705 1133032 : paw_proj%chprj = 0.0_dp
706 33764 : DO k = 1, mp
707 32150 : ico = paw_proj%n2oindex(k)
708 2232300 : paw_proj%chprj(:, k) = cprj_h(:, ico)
709 : END DO
710 1614 : DEALLOCATE (cprj_s, cprj_h, gcchprj, gccprj)
711 :
712 6456 : 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 242392 : 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 242392 : IF (ASSOCIATED(paw_proj_set)) THEN
757 242392 : IF (PRESENT(csprj)) csprj => paw_proj_set%csprj
758 242392 : IF (PRESENT(chprj)) chprj => paw_proj_set%chprj
759 242392 : IF (PRESENT(local_oce_sphi_h)) local_oce_sphi_h => paw_proj_set%local_oce_sphi_h
760 242392 : IF (PRESENT(local_oce_sphi_s)) local_oce_sphi_s => paw_proj_set%local_oce_sphi_s
761 242392 : IF (PRESENT(first_prj)) first_prj => paw_proj_set%first_prj
762 242392 : IF (PRESENT(last_prj)) last_prj => paw_proj_set%last_prj
763 242392 : IF (PRESENT(first_prjs)) first_prjs => paw_proj_set%first_prjs
764 242392 : IF (PRESENT(maxl)) maxl = paw_proj_set%maxl
765 242392 : IF (PRESENT(ncgauprj)) ncgauprj = paw_proj_set%ncgauprj
766 242392 : IF (PRESENT(nsgauprj)) nsgauprj = paw_proj_set%nsgauprj
767 242392 : IF (PRESENT(nsatbas)) nsatbas = paw_proj_set%nsatbas
768 242392 : IF (PRESENT(nsotot)) nsotot = paw_proj_set%nsotot
769 242392 : IF (PRESENT(nprj)) nprj => paw_proj_set%nprj
770 242392 : IF (PRESENT(rcprj)) rcprj = paw_proj_set%rcprj
771 242392 : IF (PRESENT(rzetprj)) rzetprj => paw_proj_set%rzetprj
772 242392 : IF (PRESENT(zisomin)) zisomin => paw_proj_set%zisomin
773 242392 : IF (PRESENT(zetprj)) zetprj => paw_proj_set%zetprj
774 242392 : IF (PRESENT(o2nindex)) o2nindex => paw_proj_set%o2nindex
775 242392 : 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 242392 : 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 1642 : 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 1642 : IF (ASSOCIATED(paw_proj_set)) THEN
796 63052 : IF (PRESENT(rzetprj)) paw_proj_set%rzetprj(:, 0:) = rzetprj(:, 0:)
797 1642 : 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 1642 : END SUBROUTINE set_paw_proj_set
803 :
804 0 : END MODULE paw_proj_set_types
|
