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 : !> \brief Some basic routines for atomic calculations
10 : !> \author jgh
11 : !> \date 01.04.2008
12 : !> \version 1.0
13 : !>
14 : ! **************************************************************************************************
15 : MODULE atom_utils
16 : USE ai_onecenter, ONLY: sg_overlap,&
17 : sto_overlap
18 : USE ai_overlap, ONLY: overlap_ab_s,&
19 : overlap_ab_sp
20 : USE ao_util, ONLY: exp_radius
21 : USE atom_types, ONLY: &
22 : CGTO_BASIS, GTO_BASIS, NUM_BASIS, STO_BASIS, atom_basis_type, atom_gthpot_type, &
23 : atom_hfx_type, atom_potential_type, atom_state, atom_type, ecp_pseudo, eri, gth_pseudo, &
24 : lmat, no_pseudo, sgp_pseudo, upf_pseudo
25 : USE basis_set_types, ONLY: srules
26 : USE cp_files, ONLY: close_file,&
27 : get_unit_number,&
28 : open_file
29 : USE input_constants, ONLY: do_rhf_atom,&
30 : do_rks_atom,&
31 : do_rohf_atom,&
32 : do_uhf_atom,&
33 : do_uks_atom
34 : USE kahan_sum, ONLY: accurate_dot_product
35 : USE kinds, ONLY: default_string_length,&
36 : dp
37 : USE lapack, ONLY: lapack_ssyev
38 : USE mathconstants, ONLY: dfac,&
39 : fac,&
40 : fourpi,&
41 : maxfac,&
42 : pi,&
43 : rootpi
44 : USE mathlib, ONLY: binomial_gen,&
45 : invmat_symm
46 : USE orbital_pointers, ONLY: deallocate_orbital_pointers,&
47 : init_orbital_pointers
48 : USE orbital_transformation_matrices, ONLY: deallocate_spherical_harmonics,&
49 : init_spherical_harmonics
50 : USE periodic_table, ONLY: nelem,&
51 : ptable
52 : USE physcon, ONLY: bohr
53 : USE qs_grid_atom, ONLY: grid_atom_type
54 : USE splines, ONLY: spline3ders
55 : USE string_utilities, ONLY: uppercase
56 : #include "./base/base_uses.f90"
57 :
58 : IMPLICIT NONE
59 :
60 : PRIVATE
61 :
62 : CHARACTER(len=*), PARAMETER, PRIVATE :: moduleN = 'atom_utils'
63 :
64 : PUBLIC :: atom_condnumber, atom_completeness, atom_basis_condnum
65 : PUBLIC :: atom_set_occupation, get_maxl_occ, get_maxn_occ
66 : PUBLIC :: atom_denmat, atom_density, atom_core_density
67 : PUBLIC :: integrate_grid, atom_trace, atom_solve
68 : PUBLIC :: coulomb_potential_numeric, coulomb_potential_analytic
69 : PUBLIC :: exchange_numeric, exchange_semi_analytic
70 : PUBLIC :: numpot_matrix, ceri_contract, eeri_contract, err_matrix
71 : PUBLIC :: slater_density, wigner_slater_functional, atom_local_potential
72 : PUBLIC :: atom_orbital_charge, atom_orbital_nodes, atom_consistent_method
73 : PUBLIC :: atom_orbital_max, atom_wfnr0, get_rho0
74 : PUBLIC :: contract2, contract4, contract2add
75 : ! ZMP added public subroutines
76 : PUBLIC :: atom_read_external_density
77 : PUBLIC :: atom_read_external_vxc
78 : PUBLIC :: atom_read_zmp_restart
79 : PUBLIC :: atom_write_zmp_restart
80 :
81 : !-----------------------------------------------------------------------------!
82 :
83 : INTERFACE integrate_grid
84 : MODULE PROCEDURE integrate_grid_function1, &
85 : integrate_grid_function2, &
86 : integrate_grid_function3
87 : END INTERFACE
88 :
89 : ! **************************************************************************************************
90 :
91 : CONTAINS
92 :
93 : ! **************************************************************************************************
94 : !> \brief Set occupation of atomic orbitals.
95 : !> \param ostring list of electronic states
96 : !> \param occupation ...
97 : !> \param wfnocc ...
98 : !> \param multiplicity ...
99 : !> \par History
100 : !> * 11.2009 unrestricted KS and HF methods [Juerg Hutter]
101 : !> * 08.2008 created [Juerg Hutter]
102 : ! **************************************************************************************************
103 470 : SUBROUTINE atom_set_occupation(ostring, occupation, wfnocc, multiplicity)
104 : CHARACTER(LEN=default_string_length), &
105 : DIMENSION(:), POINTER :: ostring
106 : REAL(Kind=dp), DIMENSION(0:lmat, 10) :: occupation, wfnocc
107 : INTEGER, INTENT(OUT), OPTIONAL :: multiplicity
108 :
109 : CHARACTER(len=2) :: elem
110 : CHARACTER(LEN=default_string_length) :: pstring
111 : INTEGER :: i, i1, i2, ielem, is, jd, jf, jp, js, k, &
112 : l, mult, n, no
113 : REAL(Kind=dp) :: e0, el, oo
114 :
115 470 : occupation = 0._dp
116 :
117 470 : CPASSERT(ASSOCIATED(ostring))
118 470 : CPASSERT(SIZE(ostring) > 0)
119 :
120 470 : no = SIZE(ostring)
121 :
122 470 : is = 1
123 : ! look for multiplicity
124 470 : mult = -1 !not specified
125 470 : IF (is <= no) THEN
126 470 : IF (INDEX(ostring(is), "(") /= 0) THEN
127 38 : i1 = INDEX(ostring(is), "(")
128 38 : i2 = INDEX(ostring(is), ")")
129 38 : CPASSERT((i2 - i1 - 1 > 0) .AND. (i2 - i1 - 1 < 3))
130 38 : elem = ostring(is) (i1 + 1:i2 - 1)
131 38 : IF (INDEX(elem, "HS") /= 0) THEN
132 24 : mult = -2 !High spin
133 14 : ELSE IF (INDEX(elem, "LS") /= 0) THEN
134 2 : mult = -3 !Low spin
135 : ELSE
136 12 : READ (elem, *) mult
137 : END IF
138 : is = is + 1
139 : END IF
140 : END IF
141 :
142 470 : IF (is <= no) THEN
143 470 : IF (INDEX(ostring(is), "CORE") /= 0) is = is + 1 !Pseudopotential detected
144 : END IF
145 470 : IF (is <= no) THEN
146 470 : IF (INDEX(ostring(is), "none") /= 0) is = is + 1 !no electrons, used with CORE
147 : END IF
148 470 : IF (is <= no) THEN
149 444 : IF (INDEX(ostring(is), "[") /= 0) THEN
150 : ! core occupation from element [XX]
151 308 : i1 = INDEX(ostring(is), "[")
152 308 : i2 = INDEX(ostring(is), "]")
153 308 : CPASSERT((i2 - i1 - 1 > 0) .AND. (i2 - i1 - 1 < 3))
154 308 : elem = ostring(is) (i1 + 1:i2 - 1)
155 308 : ielem = 0
156 9572 : DO k = 1, nelem
157 9572 : IF (elem == ptable(k)%symbol) THEN
158 : ielem = k
159 : EXIT
160 : END IF
161 : END DO
162 308 : CPASSERT(ielem /= 0)
163 1540 : DO l = 0, MIN(lmat, UBOUND(ptable(ielem)%e_conv, 1))
164 1232 : el = 2._dp*(2._dp*REAL(l, dp) + 1._dp)
165 1232 : e0 = ptable(ielem)%e_conv(l)
166 2846 : DO k = 1, 10
167 2538 : occupation(l, k) = MIN(el, e0)
168 2538 : e0 = e0 - el
169 2538 : IF (e0 <= 0._dp) EXIT
170 : END DO
171 : END DO
172 308 : is = is + 1
173 : END IF
174 :
175 : END IF
176 :
177 1200 : DO i = is, no
178 730 : pstring = ostring(i)
179 730 : CALL uppercase(pstring)
180 730 : js = INDEX(pstring, "S")
181 730 : jp = INDEX(pstring, "P")
182 730 : jd = INDEX(pstring, "D")
183 730 : jf = INDEX(pstring, "F")
184 730 : CPASSERT(js + jp + jd + jf > 0)
185 730 : IF (js > 0) THEN
186 386 : CPASSERT(jp + jd + jf == 0)
187 386 : READ (pstring(1:js - 1), *) n
188 386 : READ (pstring(js + 1:), *) oo
189 386 : CPASSERT(n > 0)
190 386 : CPASSERT(oo >= 0._dp)
191 386 : CPASSERT(occupation(0, n) == 0)
192 386 : occupation(0, n) = oo
193 : END IF
194 730 : IF (jp > 0) THEN
195 134 : CPASSERT(js + jd + jf == 0)
196 134 : READ (pstring(1:jp - 1), *) n
197 134 : READ (pstring(jp + 1:), *) oo
198 134 : CPASSERT(n > 1)
199 134 : CPASSERT(oo >= 0._dp)
200 134 : CPASSERT(occupation(1, n - 1) == 0)
201 134 : occupation(1, n - 1) = oo
202 : END IF
203 730 : IF (jd > 0) THEN
204 122 : CPASSERT(js + jp + jf == 0)
205 122 : READ (pstring(1:jd - 1), *) n
206 122 : READ (pstring(jd + 1:), *) oo
207 122 : CPASSERT(n > 2)
208 122 : CPASSERT(oo >= 0._dp)
209 122 : CPASSERT(occupation(2, n - 2) == 0)
210 122 : occupation(2, n - 2) = oo
211 : END IF
212 1200 : IF (jf > 0) THEN
213 88 : CPASSERT(js + jp + jd == 0)
214 88 : READ (pstring(1:jf - 1), *) n
215 88 : READ (pstring(jf + 1:), *) oo
216 88 : CPASSERT(n > 3)
217 88 : CPASSERT(oo >= 0._dp)
218 88 : CPASSERT(occupation(3, n - 3) == 0)
219 88 : occupation(3, n - 3) = oo
220 : END IF
221 :
222 : END DO
223 :
224 470 : wfnocc = 0._dp
225 3290 : DO l = 0, lmat
226 : k = 0
227 31490 : DO i = 1, 10
228 31020 : IF (occupation(l, i) /= 0._dp) THEN
229 2744 : k = k + 1
230 2744 : wfnocc(l, k) = occupation(l, i)
231 : END IF
232 : END DO
233 : END DO
234 :
235 : !Check for consistency with multiplicity
236 470 : IF (mult /= -1) THEN
237 : ! count open shells
238 : js = 0
239 266 : DO l = 0, lmat
240 228 : k = 2*(2*l + 1)
241 2546 : DO i = 1, 10
242 2508 : IF (wfnocc(l, i) /= 0._dp .AND. wfnocc(l, i) /= REAL(k, dp)) THEN
243 56 : js = js + 1
244 56 : i1 = l
245 56 : i2 = i
246 : END IF
247 : END DO
248 : END DO
249 :
250 38 : IF (js == 0 .AND. mult == -2) mult = 1
251 38 : IF (js == 0 .AND. mult == -3) mult = 1
252 38 : IF (js == 0) THEN
253 2 : CPASSERT(mult == 1)
254 : END IF
255 38 : IF (js == 1) THEN
256 16 : l = i1
257 16 : i = i2
258 16 : k = NINT(wfnocc(l, i))
259 16 : IF (k > (2*l + 1)) k = 2*(2*l + 1) - k
260 16 : IF (mult == -2) mult = k + 1
261 16 : IF (mult == -3) mult = MOD(k, 2) + 1
262 16 : CPASSERT(MOD(k + 1 - mult, 2) == 0)
263 : END IF
264 38 : IF (js > 1 .AND. mult /= -2) THEN
265 0 : CPASSERT(mult == -2)
266 : END IF
267 :
268 : END IF
269 :
270 470 : IF (PRESENT(multiplicity)) multiplicity = mult
271 :
272 470 : END SUBROUTINE atom_set_occupation
273 :
274 : ! **************************************************************************************************
275 : !> \brief Return the maximum orbital quantum number of occupied orbitals.
276 : !> \param occupation ...
277 : !> \return ...
278 : !> \par History
279 : !> * 08.2008 created [Juerg Hutter]
280 : ! **************************************************************************************************
281 13096 : PURE FUNCTION get_maxl_occ(occupation) RESULT(maxl)
282 : REAL(Kind=dp), DIMENSION(0:lmat, 10), INTENT(IN) :: occupation
283 : INTEGER :: maxl
284 :
285 : INTEGER :: l
286 :
287 13096 : maxl = 0
288 91672 : DO l = 0, lmat
289 877432 : IF (SUM(occupation(l, :)) /= 0._dp) maxl = l
290 : END DO
291 :
292 13096 : END FUNCTION get_maxl_occ
293 :
294 : ! **************************************************************************************************
295 : !> \brief Return the maximum principal quantum number of occupied orbitals.
296 : !> \param occupation ...
297 : !> \return ...
298 : !> \par History
299 : !> * 08.2008 created [Juerg Hutter]
300 : ! **************************************************************************************************
301 21809 : PURE FUNCTION get_maxn_occ(occupation) RESULT(maxn)
302 : REAL(Kind=dp), DIMENSION(0:lmat, 10), INTENT(IN) :: occupation
303 : INTEGER, DIMENSION(0:lmat) :: maxn
304 :
305 : INTEGER :: k, l
306 :
307 152663 : maxn = 0
308 152663 : DO l = 0, lmat
309 1461203 : DO k = 1, 10
310 1439394 : IF (occupation(l, k) /= 0._dp) maxn(l) = maxn(l) + 1
311 : END DO
312 : END DO
313 :
314 21809 : END FUNCTION get_maxn_occ
315 :
316 : ! **************************************************************************************************
317 : !> \brief Calculate a density matrix using atomic orbitals.
318 : !> \param pmat electron density matrix
319 : !> \param wfn atomic wavefunctions
320 : !> \param nbas number of basis functions
321 : !> \param occ occupation numbers
322 : !> \param maxl maximum angular momentum to consider
323 : !> \param maxn maximum principal quantum number for each angular momentum
324 : !> \par History
325 : !> * 08.2008 created [Juerg Hutter]
326 : ! **************************************************************************************************
327 78925 : PURE SUBROUTINE atom_denmat(pmat, wfn, nbas, occ, maxl, maxn)
328 : REAL(KIND=dp), DIMENSION(:, :, 0:), INTENT(INOUT) :: pmat
329 : REAL(KIND=dp), DIMENSION(:, :, 0:), INTENT(IN) :: wfn
330 : INTEGER, DIMENSION(0:lmat), INTENT(IN) :: nbas
331 : REAL(KIND=dp), DIMENSION(0:, :), INTENT(IN) :: occ
332 : INTEGER, INTENT(IN) :: maxl
333 : INTEGER, DIMENSION(0:lmat), INTENT(IN) :: maxn
334 :
335 : INTEGER :: i, j, k, l, n
336 :
337 37989991 : pmat = 0._dp
338 78925 : n = SIZE(wfn, 2)
339 226006 : DO l = 0, maxl
340 380833 : DO i = 1, MIN(n, maxn(l))
341 1635414 : DO k = 1, nbas(l)
342 33166297 : DO j = 1, nbas(l)
343 33011470 : pmat(j, k, l) = pmat(j, k, l) + occ(l, i)*wfn(j, i, l)*wfn(k, i, l)
344 : END DO
345 : END DO
346 : END DO
347 : END DO
348 :
349 78925 : END SUBROUTINE atom_denmat
350 :
351 : ! **************************************************************************************************
352 : !> \brief Map the electron density on an atomic radial grid.
353 : !> \param density computed electron density
354 : !> \param pmat electron density matrix
355 : !> \param basis atomic basis set
356 : !> \param maxl maximum angular momentum to consider
357 : !> \param typ type of the matrix to map:
358 : !> RHO -- density matrix;
359 : !> DER -- first derivatives of the electron density;
360 : !> KIN -- kinetic energy density;
361 : !> LAP -- second derivatives of the electron density.
362 : !> \param rr abscissa on the radial grid (required for typ == 'KIN')
363 : !> \par History
364 : !> * 08.2008 created [Juerg Hutter]
365 : ! **************************************************************************************************
366 159163 : SUBROUTINE atom_density(density, pmat, basis, maxl, typ, rr)
367 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: density
368 : REAL(KIND=dp), DIMENSION(:, :, 0:), INTENT(IN) :: pmat
369 : TYPE(atom_basis_type), INTENT(IN) :: basis
370 : INTEGER, INTENT(IN) :: maxl
371 : CHARACTER(LEN=*), OPTIONAL :: typ
372 : REAL(KIND=dp), DIMENSION(:), INTENT(IN), OPTIONAL :: rr
373 :
374 : CHARACTER(LEN=3) :: my_typ
375 : INTEGER :: i, j, l, n
376 : REAL(KIND=dp) :: ff
377 :
378 159163 : my_typ = "RHO"
379 159163 : IF (PRESENT(typ)) my_typ = typ(1:3)
380 159163 : IF (my_typ == "KIN") THEN
381 112 : CPASSERT(PRESENT(rr))
382 : END IF
383 :
384 64970837 : density = 0._dp
385 491474 : DO l = 0, maxl
386 332311 : n = basis%nbas(l)
387 2593900 : DO i = 1, n
388 20754596 : DO j = i, n
389 18319859 : ff = pmat(i, j, l)
390 18319859 : IF (i /= j) ff = 2._dp*pmat(i, j, l)
391 20422285 : IF (my_typ == "RHO") THEN
392 5257856398 : density(:) = density(:) + ff*basis%bf(:, i, l)*basis%bf(:, j, l)
393 5174076 : ELSE IF (my_typ == "DER") THEN
394 : density(:) = density(:) + ff*basis%dbf(:, i, l)*basis%bf(:, j, l) &
395 1966600830 : + ff*basis%bf(:, i, l)*basis%dbf(:, j, l)
396 194846 : ELSE IF (my_typ == "KIN") THEN
397 : density(:) = density(:) + 0.5_dp*ff*( &
398 : basis%dbf(:, i, l)*basis%dbf(:, j, l) + &
399 78133246 : REAL(l*(l + 1), dp)*basis%bf(:, i, l)*basis%bf(:, j, l)/rr(:))
400 0 : ELSE IF (my_typ == "LAP") THEN
401 : density(:) = density(:) + ff*basis%ddbf(:, i, l)*basis%bf(:, j, l) &
402 : + ff*basis%bf(:, i, l)*basis%ddbf(:, j, l) &
403 : + 2._dp*ff*basis%dbf(:, i, l)*basis%bf(:, j, l)/rr(:) &
404 0 : + 2._dp*ff*basis%bf(:, i, l)*basis%dbf(:, j, l)/rr(:)
405 : ELSE
406 0 : CPABORT("")
407 : END IF
408 : END DO
409 : END DO
410 : END DO
411 : ! this factor from the product of two spherical harmonics
412 64970837 : density = density/fourpi
413 :
414 159163 : END SUBROUTINE atom_density
415 :
416 : ! **************************************************************************************************
417 : !> \brief ZMP subroutine to write external restart file.
418 : !> \param atom information about the atomic kind
419 : !> \date 07.10.2013
420 : !> \author D. Varsano [daniele.varsano@nano.cnr.it]
421 : !> \version 1.0
422 : ! **************************************************************************************************
423 0 : SUBROUTINE atom_write_zmp_restart(atom)
424 : TYPE(atom_type), INTENT(IN) :: atom
425 :
426 : INTEGER :: extunit, i, k, l, n
427 :
428 0 : extunit = get_unit_number()
429 : CALL open_file(file_name=atom%zmp_restart_file, file_status="UNKNOWN", &
430 : file_form="FORMATTED", file_action="WRITE", &
431 0 : unit_number=extunit)
432 :
433 0 : n = SIZE(atom%orbitals%wfn, 2)
434 0 : WRITE (extunit, *) atom%basis%nbas
435 0 : DO l = 0, atom%state%maxl_occ
436 0 : DO i = 1, MIN(n, atom%state%maxn_occ(l))
437 0 : DO k = 1, atom%basis%nbas(l)
438 0 : WRITE (extunit, *) atom%orbitals%wfn(k, i, l)
439 : END DO
440 : END DO
441 : END DO
442 :
443 0 : CALL close_file(unit_number=extunit)
444 :
445 0 : END SUBROUTINE atom_write_zmp_restart
446 :
447 : ! **************************************************************************************************
448 : !> \brief ZMP subroutine to read external restart file.
449 : !> \param atom information about the atomic kind
450 : !> \param doguess flag that indicates that the restart file has not been read,
451 : !> so the initial guess is required
452 : !> \param iw output file unit
453 : !> \date 07.10.2013
454 : !> \author D. Varsano [daniele.varsano@nano.cnr.it]
455 : !> \version 1.0
456 : ! **************************************************************************************************
457 0 : SUBROUTINE atom_read_zmp_restart(atom, doguess, iw)
458 : TYPE(atom_type), INTENT(INOUT) :: atom
459 : LOGICAL, INTENT(INOUT) :: doguess
460 : INTEGER, INTENT(IN) :: iw
461 :
462 : INTEGER :: er, extunit, i, k, l, maxl, n
463 : INTEGER, DIMENSION(0:lmat) :: maxn, nbas
464 :
465 0 : INQUIRE (file=atom%zmp_restart_file, exist=atom%doread)
466 :
467 0 : IF (atom%doread) THEN
468 0 : WRITE (iw, fmt="(' ZMP | Restart file found')")
469 0 : extunit = get_unit_number()
470 : CALL open_file(file_name=atom%zmp_restart_file, file_status="OLD", &
471 : file_form="FORMATTED", file_action="READ", &
472 0 : unit_number=extunit)
473 :
474 0 : READ (extunit, *, IOSTAT=er) nbas
475 :
476 0 : IF (er .NE. 0) THEN
477 0 : WRITE (iw, fmt="(' ZMP | ERROR! Restart file unreadable')")
478 0 : WRITE (iw, fmt="(' ZMP | ERROR! Starting ZMP calculation form initial atomic guess')")
479 0 : doguess = .TRUE.
480 0 : atom%doread = .FALSE.
481 0 : ELSE IF (nbas(1) .NE. atom%basis%nbas(1)) THEN
482 0 : WRITE (iw, fmt="(' ZMP | ERROR! Restart file contains a different basis set')")
483 0 : WRITE (iw, fmt="(' ZMP | ERROR! Starting ZMP calculation form initial atomic guess')")
484 0 : doguess = .TRUE.
485 0 : atom%doread = .FALSE.
486 : ELSE
487 0 : nbas = atom%basis%nbas
488 0 : maxl = atom%state%maxl_occ
489 0 : maxn = atom%state%maxn_occ
490 0 : n = SIZE(atom%orbitals%wfn, 2)
491 0 : DO l = 0, atom%state%maxl_occ
492 0 : DO i = 1, MIN(n, atom%state%maxn_occ(l))
493 0 : DO k = 1, atom%basis%nbas(l)
494 0 : READ (extunit, *) atom%orbitals%wfn(k, i, l)
495 : END DO
496 : END DO
497 : END DO
498 0 : doguess = .FALSE.
499 : END IF
500 0 : CALL close_file(unit_number=extunit)
501 : ELSE
502 0 : WRITE (iw, fmt="(' ZMP | WARNING! Restart file not found')")
503 0 : WRITE (iw, fmt="(' ZMP | WARNING! Starting ZMP calculation form initial atomic guess')")
504 : END IF
505 0 : END SUBROUTINE atom_read_zmp_restart
506 :
507 : ! **************************************************************************************************
508 : !> \brief ZMP subroutine to read external density from linear grid of density matrix.
509 : !> \param density external density
510 : !> \param atom information about the atomic kind
511 : !> \param iw output file unit
512 : !> \date 07.10.2013
513 : !> \author D. Varsano [daniele.varsano@nano.cnr.it]
514 : !> \version 1.0
515 : ! **************************************************************************************************
516 0 : SUBROUTINE atom_read_external_density(density, atom, iw)
517 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: density
518 : TYPE(atom_type), INTENT(INOUT) :: atom
519 : INTEGER, INTENT(IN) :: iw
520 :
521 : CHARACTER(LEN=default_string_length) :: filename
522 : INTEGER :: extunit, ir, j, k, l, maxl_occ, maxnbas, &
523 : nbas, nr
524 : LOGICAL :: ldm
525 : REAL(KIND=dp) :: rr
526 0 : REAL(KIND=dp), ALLOCATABLE :: pmatread(:, :, :)
527 :
528 0 : filename = atom%ext_file
529 0 : ldm = atom%dm
530 0 : extunit = get_unit_number()
531 :
532 : CALL open_file(file_name=filename, file_status="OLD", &
533 : file_form="FORMATTED", file_action="READ", &
534 0 : unit_number=extunit)
535 :
536 0 : IF (.NOT. ldm) THEN
537 0 : READ (extunit, *) nr
538 :
539 0 : IF (nr .NE. atom%basis%grid%nr) THEN
540 0 : IF (iw > 0) WRITE (iw, fmt="(' ZMP | ERROR! External grid dimension ',I4,' differs from internal grid ',I4)") &
541 0 : nr, atom%basis%grid%nr
542 0 : IF (iw > 0) WRITE (iw, fmt="(' ZMP | ERROR! Stopping ZMP calculation')")
543 0 : CPABORT("")
544 : END IF
545 :
546 0 : DO ir = 1, nr
547 0 : READ (extunit, *) rr, density(ir)
548 0 : IF (ABS(rr - atom%basis%grid%rad(ir)) .GT. atom%zmpgrid_tol) THEN
549 0 : IF (iw > 0) WRITE (iw, fmt="(' ZMP | ERROR! Grid points do not coincide: ')")
550 0 : IF (iw > 0) WRITE (iw, fmt='(" ZMP |",T20,"R_out[bohr]",T36,"R_in[bohr]",T61,"R_diff[bohr]")')
551 0 : IF (iw > 0) WRITE (iw, fmt='(" ZMP |",T14,E24.15,T39,E24.15,T64,E24.15)') &
552 0 : rr, atom%basis%grid%rad(ir), ABS(rr - atom%basis%grid%rad(ir))
553 0 : CPABORT("")
554 : END IF
555 : END DO
556 0 : CALL close_file(unit_number=extunit)
557 : ELSE
558 0 : READ (extunit, *) maxl_occ
559 0 : maxnbas = MAXVAL(atom%basis%nbas)
560 0 : ALLOCATE (pmatread(maxnbas, maxnbas, 0:maxl_occ))
561 0 : pmatread = 0.0
562 0 : DO l = 0, maxl_occ
563 0 : nbas = atom%basis%nbas(l)
564 0 : READ (extunit, *) ! Read empty line
565 0 : DO k = 1, nbas
566 0 : READ (extunit, *) (pmatread(j, k, l), j=1, k)
567 0 : DO j = 1, k
568 0 : pmatread(k, j, l) = pmatread(j, k, l)
569 : END DO
570 : END DO
571 : END DO
572 :
573 0 : CALL close_file(unit_number=extunit)
574 :
575 0 : CALL atom_density(density, pmatread, atom%basis, maxl_occ, typ="RHO")
576 :
577 0 : extunit = get_unit_number()
578 : CALL open_file(file_name="rho_target.dat", file_status="UNKNOWN", &
579 0 : file_form="FORMATTED", file_action="WRITE", unit_number=extunit)
580 :
581 0 : IF (iw > 0) WRITE (iw, fmt="(' ZMP | Writing target density from density matrix')")
582 :
583 0 : WRITE (extunit, fmt='("# Target density built from density matrix : ",A20)') filename
584 0 : WRITE (extunit, fmt='("#",T10,"R[bohr]",T36,"Rho[au]")')
585 :
586 0 : nr = atom%basis%grid%nr
587 :
588 0 : DO ir = 1, nr
589 : WRITE (extunit, fmt='(T1,E24.15,T26,E24.15)') &
590 0 : atom%basis%grid%rad(ir), density(ir)
591 : END DO
592 0 : DEALLOCATE (pmatread)
593 :
594 0 : CALL close_file(unit_number=extunit)
595 :
596 : END IF
597 :
598 0 : END SUBROUTINE atom_read_external_density
599 :
600 : ! **************************************************************************************************
601 : !> \brief ZMP subroutine to read external v_xc in the atomic code.
602 : !> \param vxc external exchange-correlation potential
603 : !> \param atom information about the atomic kind
604 : !> \param iw output file unit
605 : !> \author D. Varsano [daniele.varsano@nano.cnr.it]
606 : ! **************************************************************************************************
607 0 : SUBROUTINE atom_read_external_vxc(vxc, atom, iw)
608 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: vxc
609 : TYPE(atom_type), INTENT(INOUT) :: atom
610 : INTEGER, INTENT(IN) :: iw
611 :
612 : CHARACTER(LEN=default_string_length) :: adum, filename
613 : INTEGER :: extunit, ir, nr
614 : REAL(KIND=dp) :: rr
615 :
616 0 : filename = atom%ext_vxc_file
617 0 : extunit = get_unit_number()
618 :
619 : CALL open_file(file_name=filename, file_status="OLD", &
620 : file_form="FORMATTED", file_action="READ", &
621 0 : unit_number=extunit)
622 :
623 0 : READ (extunit, *)
624 0 : READ (extunit, *) adum, nr
625 :
626 0 : IF (nr .NE. atom%basis%grid%nr) THEN
627 0 : IF (iw > 0) WRITE (iw, fmt="(' ZMP | ERROR! External grid dimension ',I4,' differs from internal grid ',I4)") &
628 0 : nr, atom%basis%grid%nr
629 0 : IF (iw > 0) WRITE (iw, fmt="(' ZMP | ERROR! Stopping ZMP calculation')")
630 0 : CPABORT("")
631 : END IF
632 0 : DO ir = 1, nr
633 0 : READ (extunit, *) rr, vxc(ir)
634 0 : IF (ABS(rr - atom%basis%grid%rad(ir)) .GT. atom%zmpvxcgrid_tol) THEN
635 0 : IF (iw > 0) WRITE (iw, fmt="(' ZMP | ERROR! Grid points do not coincide: ')")
636 0 : IF (iw > 0) WRITE (iw, fmt='(" ZMP |",T20,"R_out[bohr]",T36,"R_in[bohr]",T61,"R_diff[bohr]")')
637 0 : IF (iw > 0) WRITE (iw, fmt='(" ZMP |",T14,E24.15,T39,E24.15,T64,E24.15)') &
638 0 : rr, atom%basis%grid%rad(ir), ABS(rr - atom%basis%grid%rad(ir))
639 0 : CPABORT("")
640 : END IF
641 : END DO
642 :
643 0 : END SUBROUTINE atom_read_external_vxc
644 :
645 : ! **************************************************************************************************
646 : !> \brief ...
647 : !> \param charge ...
648 : !> \param wfn ...
649 : !> \param rcov ...
650 : !> \param l ...
651 : !> \param basis ...
652 : ! **************************************************************************************************
653 1374 : PURE SUBROUTINE atom_orbital_charge(charge, wfn, rcov, l, basis)
654 : REAL(KIND=dp), INTENT(OUT) :: charge
655 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: wfn
656 : REAL(KIND=dp), INTENT(IN) :: rcov
657 : INTEGER, INTENT(IN) :: l
658 : TYPE(atom_basis_type), INTENT(IN) :: basis
659 :
660 : INTEGER :: i, j, m, n
661 : REAL(KIND=dp) :: ff
662 1374 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: den
663 :
664 1374 : charge = 0._dp
665 1374 : m = SIZE(basis%bf, 1)
666 4122 : ALLOCATE (den(m))
667 1374 : n = basis%nbas(l)
668 520274 : den = 0._dp
669 30694 : DO i = 1, n
670 884334 : DO j = 1, n
671 853640 : ff = wfn(i)*wfn(j)
672 321059560 : den(1:m) = den(1:m) + ff*basis%bf(1:m, i, l)*basis%bf(1:m, j, l)
673 : END DO
674 : END DO
675 520274 : DO i = 1, m
676 520274 : IF (basis%grid%rad(i) > rcov) den(i) = 0._dp
677 : END DO
678 520274 : charge = SUM(den(1:m)*basis%grid%wr(1:m))
679 1374 : DEALLOCATE (den)
680 :
681 1374 : END SUBROUTINE atom_orbital_charge
682 :
683 : ! **************************************************************************************************
684 : !> \brief ...
685 : !> \param corden ...
686 : !> \param potential ...
687 : !> \param typ ...
688 : !> \param rr ...
689 : !> \par History
690 : !> * 01.2017 rewritten [Juerg Hutter]
691 : !> * 03.2010 extension of GTH pseudo-potential definition [Juerg Hutter]
692 : ! **************************************************************************************************
693 1692 : SUBROUTINE atom_core_density(corden, potential, typ, rr)
694 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: corden
695 : TYPE(atom_potential_type), INTENT(IN) :: potential
696 : CHARACTER(LEN=*), OPTIONAL :: typ
697 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rr
698 :
699 : CHARACTER(LEN=3) :: my_typ
700 : INTEGER :: i, j, m, n
701 : LOGICAL :: reverse
702 : REAL(KIND=dp) :: a, a2, cval, fb
703 1692 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: fe, rc, rhoc, rval
704 :
705 1692 : my_typ = "RHO"
706 1692 : IF (PRESENT(typ)) my_typ = typ(1:3)
707 :
708 3384 : SELECT CASE (potential%ppot_type)
709 : CASE (no_pseudo, ecp_pseudo)
710 : ! we do nothing
711 : CASE (gth_pseudo)
712 1692 : IF (potential%gth_pot%nlcc) THEN
713 1692 : m = SIZE(corden)
714 8460 : ALLOCATE (fe(m), rc(m))
715 1692 : n = potential%gth_pot%nexp_nlcc
716 3384 : DO i = 1, n
717 1692 : a = potential%gth_pot%alpha_nlcc(i)
718 1692 : a2 = a*a
719 : ! note that all terms are computed with rc, not rr
720 678492 : rc(:) = rr(:)/a
721 678492 : fe(:) = EXP(-0.5_dp*rc(:)*rc(:))
722 5112 : DO j = 1, potential%gth_pot%nct_nlcc(i)
723 1728 : cval = potential%gth_pot%cval_nlcc(j, i)
724 3420 : IF (my_typ == "RHO") THEN
725 368519 : corden(:) = corden(:) + fe(:)*rc**(2*j - 2)*cval
726 809 : ELSE IF (my_typ == "DER") THEN
727 324409 : corden(:) = corden(:) - fe(:)*rc**(2*j - 1)*cval/a
728 809 : IF (j > 1) THEN
729 0 : corden(:) = corden(:) + REAL(2*j - 2, dp)*fe(:)*rc**(2*j - 3)*cval/a
730 : END IF
731 0 : ELSE IF (my_typ == "LAP") THEN
732 0 : fb = 2._dp*cval/a
733 0 : corden(:) = corden(:) - fb*fe(:)/rr(:)*rc**(2*j - 1)
734 0 : corden(:) = corden(:) + fe(:)*rc**(2*j)*cval/a2
735 0 : IF (j > 1) THEN
736 0 : corden(:) = corden(:) + fb*REAL(2*j - 2, dp)*fe(:)/rr(:)*rc**(2*j - 3)
737 0 : corden(:) = corden(:) + REAL((2*j - 2)*(2*j - 3), dp)*fe(:)*rc**(2*j - 4)*cval/a2
738 0 : corden(:) = corden(:) - REAL(2*j - 2, dp)*fe(:)*rc**(2*j - 2)*cval/a2
739 : END IF
740 : ELSE
741 0 : CPABORT("")
742 : END IF
743 : END DO
744 : END DO
745 1692 : DEALLOCATE (fe, rc)
746 : END IF
747 : CASE (upf_pseudo)
748 0 : IF (potential%upf_pot%core_correction) THEN
749 0 : m = SIZE(corden)
750 0 : n = potential%upf_pot%mesh_size
751 0 : reverse = .FALSE.
752 0 : IF (rr(1) > rr(m)) reverse = .TRUE.
753 0 : ALLOCATE (rhoc(m), rval(m))
754 0 : IF (reverse) THEN
755 0 : DO i = 1, m
756 0 : rval(i) = rr(m - i + 1)
757 : END DO
758 : ELSE
759 0 : rval(1:m) = rr(1:m)
760 : END IF
761 0 : IF (my_typ == "RHO") THEN
762 : CALL spline3ders(potential%upf_pot%r(1:n), potential%upf_pot%rho_nlcc(1:n), rval(1:m), &
763 0 : ynew=rhoc(1:m))
764 0 : ELSE IF (my_typ == "DER") THEN
765 : CALL spline3ders(potential%upf_pot%r(1:n), potential%upf_pot%rho_nlcc(1:n), rval(1:m), &
766 0 : dynew=rhoc(1:m))
767 0 : ELSE IF (my_typ == "LAP") THEN
768 : CALL spline3ders(potential%upf_pot%r(1:n), potential%upf_pot%rho_nlcc(1:n), rval(1:m), &
769 0 : d2ynew=rhoc(1:m))
770 : ELSE
771 0 : CPABORT("")
772 : END IF
773 0 : IF (reverse) THEN
774 0 : DO i = 1, m
775 0 : rval(i) = rr(m - i + 1)
776 0 : corden(i) = corden(i) + rhoc(m - i + 1)
777 : END DO
778 : ELSE
779 0 : corden(1:m) = corden(1:m) + rhoc(1:m)
780 : END IF
781 0 : DEALLOCATE (rhoc, rval)
782 : END IF
783 : CASE (sgp_pseudo)
784 0 : IF (potential%sgp_pot%has_nlcc) THEN
785 0 : CPABORT("not implemented")
786 : END IF
787 : CASE DEFAULT
788 1692 : CPABORT("Unknown PP type")
789 : END SELECT
790 :
791 1692 : END SUBROUTINE atom_core_density
792 :
793 : ! **************************************************************************************************
794 : !> \brief ...
795 : !> \param locpot ...
796 : !> \param gthpot ...
797 : !> \param rr ...
798 : ! **************************************************************************************************
799 4 : PURE SUBROUTINE atom_local_potential(locpot, gthpot, rr)
800 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: locpot
801 : TYPE(atom_gthpot_type), INTENT(IN) :: gthpot
802 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rr
803 :
804 : INTEGER :: i, j, m, n
805 : REAL(KIND=dp) :: a
806 4 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: fe, rc
807 :
808 4 : m = SIZE(locpot)
809 20 : ALLOCATE (fe(m), rc(m))
810 2662 : rc(:) = rr(:)/gthpot%rc
811 2662 : DO i = 1, m
812 2662 : locpot(i) = -gthpot%zion*erf(rc(i)/SQRT(2._dp))/rr(i)
813 : END DO
814 4 : n = gthpot%ncl
815 2662 : fe(:) = EXP(-0.5_dp*rc(:)*rc(:))
816 12 : DO i = 1, n
817 5328 : locpot(:) = locpot(:) + fe(:)*rc**(2*i - 2)*gthpot%cl(i)
818 : END DO
819 4 : IF (gthpot%lpotextended) THEN
820 0 : DO j = 1, gthpot%nexp_lpot
821 0 : a = gthpot%alpha_lpot(j)
822 0 : rc(:) = rr(:)/a
823 0 : fe(:) = EXP(-0.5_dp*rc(:)*rc(:))
824 0 : n = gthpot%nct_lpot(j)
825 0 : DO i = 1, n
826 0 : locpot(:) = locpot(:) + fe(:)*rc**(2*i - 2)*gthpot%cval_lpot(i, j)
827 : END DO
828 : END DO
829 : END IF
830 4 : DEALLOCATE (fe, rc)
831 :
832 4 : END SUBROUTINE atom_local_potential
833 :
834 : ! **************************************************************************************************
835 : !> \brief ...
836 : !> \param rmax ...
837 : !> \param wfn ...
838 : !> \param rcov ...
839 : !> \param l ...
840 : !> \param basis ...
841 : ! **************************************************************************************************
842 80 : PURE SUBROUTINE atom_orbital_max(rmax, wfn, rcov, l, basis)
843 : REAL(KIND=dp), INTENT(OUT) :: rmax
844 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: wfn
845 : REAL(KIND=dp), INTENT(IN) :: rcov
846 : INTEGER, INTENT(IN) :: l
847 : TYPE(atom_basis_type), INTENT(IN) :: basis
848 :
849 : INTEGER :: i, m, n
850 : REAL(KIND=dp) :: ff
851 80 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dorb
852 :
853 80 : m = SIZE(basis%bf, 1)
854 240 : ALLOCATE (dorb(m))
855 80 : n = basis%nbas(l)
856 18180 : dorb = 0._dp
857 2714 : DO i = 1, n
858 2634 : ff = wfn(i)
859 600714 : dorb(1:m) = dorb(1:m) + ff*basis%dbf(1:m, i, l)
860 : END DO
861 80 : rmax = -1._dp
862 18100 : DO i = 1, m - 1
863 18100 : IF (basis%grid%rad(i) < 2*rcov) THEN
864 11962 : IF (dorb(i)*dorb(i + 1) < 0._dp) THEN
865 81 : rmax = MAX(rmax, basis%grid%rad(i))
866 : END IF
867 : END IF
868 : END DO
869 80 : DEALLOCATE (dorb)
870 :
871 80 : END SUBROUTINE atom_orbital_max
872 :
873 : ! **************************************************************************************************
874 : !> \brief ...
875 : !> \param node ...
876 : !> \param wfn ...
877 : !> \param rcov ...
878 : !> \param l ...
879 : !> \param basis ...
880 : ! **************************************************************************************************
881 590 : PURE SUBROUTINE atom_orbital_nodes(node, wfn, rcov, l, basis)
882 : INTEGER, INTENT(OUT) :: node
883 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: wfn
884 : REAL(KIND=dp), INTENT(IN) :: rcov
885 : INTEGER, INTENT(IN) :: l
886 : TYPE(atom_basis_type), INTENT(IN) :: basis
887 :
888 : INTEGER :: i, m, n
889 : REAL(KIND=dp) :: ff
890 590 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: orb
891 :
892 590 : node = 0
893 590 : m = SIZE(basis%bf, 1)
894 1770 : ALLOCATE (orb(m))
895 590 : n = basis%nbas(l)
896 230990 : orb = 0._dp
897 9670 : DO i = 1, n
898 9080 : ff = wfn(i)
899 3535270 : orb(1:m) = orb(1:m) + ff*basis%bf(1:m, i, l)
900 : END DO
901 230400 : DO i = 1, m - 1
902 230400 : IF (basis%grid%rad(i) < rcov) THEN
903 152990 : IF (orb(i)*orb(i + 1) < 0._dp) node = node + 1
904 : END IF
905 : END DO
906 590 : DEALLOCATE (orb)
907 :
908 590 : END SUBROUTINE atom_orbital_nodes
909 :
910 : ! **************************************************************************************************
911 : !> \brief ...
912 : !> \param value ...
913 : !> \param wfn ...
914 : !> \param basis ...
915 : ! **************************************************************************************************
916 296 : PURE SUBROUTINE atom_wfnr0(value, wfn, basis)
917 : REAL(KIND=dp), INTENT(OUT) :: value
918 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: wfn
919 : TYPE(atom_basis_type), INTENT(IN) :: basis
920 :
921 : INTEGER :: i, m, n
922 :
923 296 : value = 0._dp
924 115592 : m = MAXVAL(MINLOC(basis%grid%rad))
925 296 : n = basis%nbas(0)
926 5785 : DO i = 1, n
927 5785 : value = value + wfn(i)*basis%bf(m, i, 0)
928 : END DO
929 296 : END SUBROUTINE atom_wfnr0
930 :
931 : ! **************************************************************************************************
932 : !> \brief Solve the generalised eigenproblem for every angular momentum.
933 : !> \param hmat Hamiltonian matrix
934 : !> \param umat transformation matrix which reduces the overlap matrix to its unitary form
935 : !> \param orb atomic wavefunctions
936 : !> \param ener atomic orbital energies
937 : !> \param nb number of contracted basis functions
938 : !> \param nv number of linear-independent contracted basis functions
939 : !> \param maxl maximum angular momentum to consider
940 : !> \par History
941 : !> * 08.2008 created [Juerg Hutter]
942 : ! **************************************************************************************************
943 78529 : SUBROUTINE atom_solve(hmat, umat, orb, ener, nb, nv, maxl)
944 : REAL(KIND=dp), DIMENSION(:, :, 0:), INTENT(IN) :: hmat, umat
945 : REAL(KIND=dp), DIMENSION(:, :, 0:), INTENT(INOUT) :: orb
946 : REAL(KIND=dp), DIMENSION(:, 0:), INTENT(INOUT) :: ener
947 : INTEGER, DIMENSION(0:), INTENT(IN) :: nb, nv
948 : INTEGER, INTENT(IN) :: maxl
949 :
950 : INTEGER :: info, l, lwork, m, n
951 78529 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: w, work
952 78529 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: a
953 :
954 549703 : CPASSERT(ALL(nb >= nv))
955 :
956 6233677 : orb = 0._dp
957 232734 : DO l = 0, maxl
958 154205 : n = nb(l)
959 154205 : m = nv(l)
960 232734 : IF (n > 0 .AND. m > 0) THEN
961 150251 : lwork = 10*m
962 1202008 : ALLOCATE (a(n, n), w(n), work(lwork))
963 466792823 : a(1:m, 1:m) = MATMUL(TRANSPOSE(umat(1:n, 1:m, l)), MATMUL(hmat(1:n, 1:n, l), umat(1:n, 1:m, l)))
964 10416195 : CALL lapack_ssyev("V", "U", m, a(1:m, 1:m), m, w(1:m), work, lwork, info)
965 319802107 : a(1:n, 1:m) = MATMUL(umat(1:n, 1:m, l), a(1:m, 1:m))
966 :
967 150251 : m = MIN(m, SIZE(orb, 2))
968 2714493 : orb(1:n, 1:m, l) = a(1:n, 1:m)
969 360484 : ener(1:m, l) = w(1:m)
970 :
971 150251 : DEALLOCATE (a, w, work)
972 : END IF
973 : END DO
974 :
975 78529 : END SUBROUTINE atom_solve
976 :
977 : ! **************************************************************************************************
978 : !> \brief THIS FUNCTION IS NO LONGER IN USE.
979 : !> \param fun ...
980 : !> \param deps ...
981 : !> \return ...
982 : ! **************************************************************************************************
983 0 : FUNCTION prune_grid(fun, deps) RESULT(nc)
984 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: fun
985 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: deps
986 : INTEGER :: nc
987 :
988 : INTEGER :: i, nr
989 : REAL(KIND=dp) :: meps
990 :
991 0 : meps = 1.e-14_dp
992 0 : IF (PRESENT(deps)) meps = deps
993 :
994 0 : nr = SIZE(fun)
995 0 : nc = 0
996 0 : DO i = nr, 1, -1
997 0 : IF (ABS(fun(i)) > meps) THEN
998 : nc = i
999 : EXIT
1000 : END IF
1001 : END DO
1002 :
1003 0 : END FUNCTION prune_grid
1004 :
1005 : ! **************************************************************************************************
1006 : !> \brief Integrate a function on an atomic radial grid.
1007 : !> \param fun function to integrate
1008 : !> \param grid atomic radial grid
1009 : !> \return value of the integral
1010 : !> \par History
1011 : !> * 08.2008 created [Juerg Hutter]
1012 : ! **************************************************************************************************
1013 149936 : PURE FUNCTION integrate_grid_function1(fun, grid) RESULT(integral)
1014 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: fun
1015 : TYPE(grid_atom_type), INTENT(IN) :: grid
1016 : REAL(KIND=dp) :: integral
1017 :
1018 : INTEGER :: nc
1019 :
1020 149936 : nc = SIZE(fun)
1021 60898204 : integral = SUM(fun(1:nc)*grid%wr(1:nc))
1022 :
1023 149936 : END FUNCTION integrate_grid_function1
1024 :
1025 : ! **************************************************************************************************
1026 : !> \brief Integrate the product of two functions on an atomic radial grid.
1027 : !> \param fun1 first function
1028 : !> \param fun2 second function
1029 : !> \param grid atomic radial grid
1030 : !> \return value of the integral
1031 : !> \par History
1032 : !> * 08.2008 created [Juerg Hutter]
1033 : ! **************************************************************************************************
1034 693402 : PURE FUNCTION integrate_grid_function2(fun1, fun2, grid) RESULT(integral)
1035 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: fun1, fun2
1036 : TYPE(grid_atom_type), INTENT(IN) :: grid
1037 : REAL(KIND=dp) :: integral
1038 :
1039 : INTEGER :: nc
1040 :
1041 693402 : nc = MIN(SIZE(fun1), SIZE(fun2))
1042 278297426 : integral = SUM(fun1(1:nc)*fun2(1:nc)*grid%wr(1:nc))
1043 :
1044 693402 : END FUNCTION integrate_grid_function2
1045 :
1046 : ! **************************************************************************************************
1047 : !> \brief Integrate the product of three functions on an atomic radial grid.
1048 : !> \param fun1 first function
1049 : !> \param fun2 second function
1050 : !> \param fun3 third function
1051 : !> \param grid atomic radial grid
1052 : !> \return value of the integral
1053 : !> \par History
1054 : !> * 08.2008 created [Juerg Hutter]
1055 : ! **************************************************************************************************
1056 67547319 : PURE FUNCTION integrate_grid_function3(fun1, fun2, fun3, grid) RESULT(integral)
1057 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: fun1, fun2, fun3
1058 : TYPE(grid_atom_type), INTENT(IN) :: grid
1059 : REAL(KIND=dp) :: integral
1060 :
1061 : INTEGER :: nc
1062 :
1063 67547319 : nc = MIN(SIZE(fun1), SIZE(fun2), SIZE(fun3))
1064 26649269889 : integral = SUM(fun1(1:nc)*fun2(1:nc)*fun3(1:nc)*grid%wr(1:nc))
1065 :
1066 67547319 : END FUNCTION integrate_grid_function3
1067 :
1068 : ! **************************************************************************************************
1069 : !> \brief Numerically compute the Coulomb potential on an atomic radial grid.
1070 : !> \param cpot Coulomb potential on the radial grid
1071 : !> \param density electron density on the radial grid
1072 : !> \param grid atomic radial grid
1073 : !> \par History
1074 : !> * 08.2008 created [Juerg Hutter]
1075 : ! **************************************************************************************************
1076 86820 : SUBROUTINE coulomb_potential_numeric(cpot, density, grid)
1077 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: cpot
1078 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: density
1079 : TYPE(grid_atom_type), INTENT(IN) :: grid
1080 :
1081 : INTEGER :: i, nc
1082 : REAL(dp) :: int1, int2
1083 86820 : REAL(dp), DIMENSION(:), POINTER :: r, wr
1084 :
1085 86820 : nc = MIN(SIZE(cpot), SIZE(density))
1086 86820 : r => grid%rad
1087 86820 : wr => grid%wr
1088 :
1089 86820 : int1 = fourpi*integrate_grid(density, grid)
1090 86820 : int2 = 0._dp
1091 173640 : cpot(nc:) = int1/r(nc:)
1092 : ! IF (log_unit>0) WRITE(log_unit,"(A,2F10.8)") "r", int1, cpot(nc:)
1093 :
1094 : ! test that grid is decreasing
1095 86820 : CPASSERT(r(1) > r(nc))
1096 35216308 : DO i = 1, nc
1097 35129488 : cpot(i) = int1/r(i) + int2
1098 35129488 : int1 = int1 - fourpi*density(i)*wr(i)
1099 35216308 : int2 = int2 + fourpi*density(i)*wr(i)/r(i)
1100 : END DO
1101 :
1102 86820 : END SUBROUTINE coulomb_potential_numeric
1103 :
1104 : ! **************************************************************************************************
1105 : !> \brief Analytically compute the Coulomb potential on an atomic radial grid.
1106 : !> \param cpot Coulomb potential on the radial grid
1107 : !> \param pmat density matrix
1108 : !> \param basis atomic basis set
1109 : !> \param grid atomic radial grid
1110 : !> \param maxl maximum angular momentum to consider
1111 : !> \par History
1112 : !> * 08.2008 created [Juerg Hutter]
1113 : ! **************************************************************************************************
1114 38 : SUBROUTINE coulomb_potential_analytic(cpot, pmat, basis, grid, maxl)
1115 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: cpot
1116 : REAL(KIND=dp), DIMENSION(:, :, 0:), INTENT(IN) :: pmat
1117 : TYPE(atom_basis_type), INTENT(IN) :: basis
1118 : TYPE(grid_atom_type) :: grid
1119 : INTEGER, INTENT(IN) :: maxl
1120 :
1121 : INTEGER :: i, j, k, l, m, n
1122 : REAL(KIND=dp) :: a, b, ff, oab, sab
1123 38 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: erfa, expa, z
1124 38 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: unp
1125 :
1126 38 : m = SIZE(cpot)
1127 266 : ALLOCATE (erfa(1:m), expa(1:m), z(1:m))
1128 :
1129 18838 : cpot = 0._dp
1130 :
1131 172 : DO l = 0, maxl
1132 65958 : IF (MAXVAL(ABS(pmat(:, :, l))) < 1.e-14_dp) CYCLE
1133 38 : SELECT CASE (basis%basis_type)
1134 : CASE DEFAULT
1135 0 : CPABORT("")
1136 : CASE (GTO_BASIS)
1137 2266 : DO i = 1, basis%nbas(l)
1138 24774 : DO j = i, basis%nbas(l)
1139 22508 : IF (ABS(pmat(i, j, l)) < 1.e-14_dp) CYCLE
1140 22490 : ff = pmat(i, j, l)
1141 22490 : IF (i /= j) ff = 2._dp*ff
1142 22490 : a = basis%am(i, l)
1143 22490 : b = basis%am(j, l)
1144 22490 : sab = SQRT(a + b)
1145 22490 : oab = rootpi/(a + b)**(l + 1.5_dp)*ff
1146 12640090 : z(:) = sab*grid%rad(:)
1147 12640090 : DO k = 1, SIZE(erfa)
1148 12640090 : erfa(k) = oab*erf(z(k))/grid%rad(k)
1149 : END DO
1150 12640090 : expa(:) = EXP(-z(:)**2)*ff/(a + b)**(l + 1)
1151 2132 : SELECT CASE (l)
1152 : CASE DEFAULT
1153 0 : CPABORT("")
1154 : CASE (0)
1155 6153004 : cpot(:) = cpot(:) + 0.25_dp*erfa(:)
1156 : CASE (1)
1157 3956950 : cpot(:) = cpot(:) + 0.375_dp*erfa(:) - 0.25_dp*expa(:)
1158 : CASE (2)
1159 2096254 : cpot(:) = cpot(:) + 0.9375_dp*erfa(:) - expa(:)*(0.875_dp + 0.25_dp*z(:)**2)
1160 : CASE (3)
1161 455290 : cpot(:) = cpot(:) + 3.28125_dp*erfa(:) - expa(:)*(3.5625_dp + 1.375_dp*z(:)**2 + 0.25*z(:)**4)
1162 : END SELECT
1163 : END DO
1164 : END DO
1165 : CASE (CGTO_BASIS)
1166 0 : n = basis%nprim(l)
1167 0 : m = basis%nbas(l)
1168 0 : ALLOCATE (unp(n, n))
1169 :
1170 0 : unp(1:n, 1:n) = MATMUL(MATMUL(basis%cm(1:n, 1:m, l), pmat(1:m, 1:m, l)), &
1171 0 : TRANSPOSE(basis%cm(1:n, 1:m, l)))
1172 0 : DO i = 1, basis%nprim(l)
1173 0 : DO j = i, basis%nprim(l)
1174 0 : IF (ABS(unp(i, j)) < 1.e-14_dp) CYCLE
1175 0 : ff = unp(i, j)
1176 0 : IF (i /= j) ff = 2._dp*ff
1177 0 : a = basis%am(i, l)
1178 0 : b = basis%am(j, l)
1179 0 : sab = SQRT(a + b)
1180 0 : oab = rootpi/(a + b)**(l + 1.5_dp)*ff
1181 0 : z(:) = sab*grid%rad(:)
1182 0 : DO k = 1, SIZE(erfa)
1183 0 : erfa(k) = oab*erf(z(k))/grid%rad(k)
1184 : END DO
1185 0 : expa(:) = EXP(-z(:)**2)*ff/(a + b)**(l + 1)
1186 0 : SELECT CASE (l)
1187 : CASE DEFAULT
1188 0 : CPABORT("")
1189 : CASE (0)
1190 0 : cpot(:) = cpot(:) + 0.25_dp*erfa(:)
1191 : CASE (1)
1192 0 : cpot(:) = cpot(:) + 0.375_dp*erfa(:) - 0.25_dp*expa(:)
1193 : CASE (2)
1194 0 : cpot(:) = cpot(:) + 0.9375_dp*erfa(:) - expa(:)*(0.875_dp + 0.25_dp*z(:)**2)
1195 : CASE (3)
1196 0 : cpot(:) = cpot(:) + 3.28125_dp*erfa(:) - expa(:)*(3.5625_dp + 1.375_dp*z(:)**2 + 0.25*z(:)**4)
1197 : END SELECT
1198 : END DO
1199 : END DO
1200 :
1201 134 : DEALLOCATE (unp)
1202 : END SELECT
1203 : END DO
1204 38 : DEALLOCATE (erfa, expa, z)
1205 :
1206 38 : END SUBROUTINE coulomb_potential_analytic
1207 :
1208 : ! **************************************************************************************************
1209 : !> \brief Calculate the exchange potential numerically.
1210 : !> \param kmat Kohn-Sham matrix
1211 : !> \param state electronic state
1212 : !> \param occ occupation numbers
1213 : !> \param wfn wavefunctions
1214 : !> \param basis atomic basis set
1215 : !> \param hfx_pot potential parameters from Hartree-Fock section
1216 : !> \par History
1217 : !> * 01.2009 created [Juerg Hutter]
1218 : ! **************************************************************************************************
1219 20518 : SUBROUTINE exchange_numeric(kmat, state, occ, wfn, basis, hfx_pot)
1220 : REAL(KIND=dp), DIMENSION(:, :, 0:), INTENT(INOUT) :: kmat
1221 : TYPE(atom_state), INTENT(IN) :: state
1222 : REAL(KIND=dp), DIMENSION(0:, :), INTENT(IN) :: occ
1223 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: wfn
1224 : TYPE(atom_basis_type), INTENT(IN) :: basis
1225 : TYPE(atom_hfx_type), INTENT(IN) :: hfx_pot
1226 :
1227 : INTEGER :: i, ia, ib, k, lad, lbc, lh, ll, nbas, &
1228 : norb, nr, nu
1229 : REAL(KIND=dp) :: almn
1230 20518 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cpot, nai, nbi, pot
1231 20518 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: orb
1232 : REAL(KIND=dp), DIMENSION(0:maxfac) :: arho
1233 :
1234 20518 : arho = 0._dp
1235 20518 : DO ll = 0, maxfac, 2
1236 328288 : lh = ll/2
1237 328288 : arho(ll) = fac(ll)/fac(lh)**2
1238 : END DO
1239 :
1240 3946378 : kmat = 0._dp
1241 :
1242 20518 : nr = basis%grid%nr
1243 184662 : ALLOCATE (nai(nr), nbi(nr), cpot(nr), pot(nr))
1244 :
1245 43652 : DO lad = 0, state%maxl_calc
1246 73398 : DO lbc = 0, state%maxl_occ
1247 29746 : norb = state%maxn_occ(lbc)
1248 29746 : nbas = basis%nbas(lbc)
1249 : ! calculate orbitals for angmom lbc
1250 118948 : ALLOCATE (orb(nr, norb))
1251 20344532 : orb = 0._dp
1252 80532 : DO i = 1, norb
1253 595990 : DO k = 1, nbas
1254 205052644 : orb(:, i) = orb(:, i) + wfn(k, i, lbc)*basis%bf(:, k, lbc)
1255 : END DO
1256 : END DO
1257 65960 : DO nu = ABS(lad - lbc), lad + lbc, 2
1258 : almn = arho(-lad + lbc + nu)*arho(lad - lbc + nu)*arho(lad + lbc - nu)/(REAL(lad + lbc + nu + 1, dp) &
1259 36214 : *arho(lad + lbc + nu))
1260 36214 : almn = -0.5_dp*almn
1261 :
1262 350964 : DO ia = 1, basis%nbas(lad)
1263 1019040 : DO i = 1, norb
1264 277939422 : nai(:) = orb(:, i)*basis%bf(:, ia, lad)
1265 277939422 : cpot = 0.0_dp
1266 697822 : IF (hfx_pot%scale_coulomb /= 0.0_dp) THEN
1267 662850 : CALL potential_coulomb_numeric(pot, nai, nu, basis%grid)
1268 265802850 : cpot(:) = cpot(:) + pot(:)*hfx_pot%scale_coulomb
1269 : END IF
1270 697822 : IF (hfx_pot%scale_longrange /= 0.0_dp) THEN
1271 34972 : IF (hfx_pot%do_gh) THEN
1272 : CALL potential_longrange_numeric_gh(pot, nai, nu, basis%grid, hfx_pot%omega, &
1273 17486 : hfx_pot%kernel(:, :, nu))
1274 : ELSE
1275 : CALL potential_longrange_numeric(pot, nai, nu, basis%grid, hfx_pot%omega, &
1276 17486 : hfx_pot%kernel(:, :, nu))
1277 : END IF
1278 12136572 : cpot(:) = cpot(:) + pot(:)*hfx_pot%scale_longrange
1279 : END IF
1280 13514268 : DO ib = 1, basis%nbas(lad)
1281 : kmat(ia, ib, lad) = kmat(ia, ib, lad) + almn*occ(lbc, i)* &
1282 13229264 : integrate_grid(cpot, orb(:, i), basis%bf(:, ib, lad), basis%grid)
1283 : END DO
1284 : END DO
1285 : END DO
1286 :
1287 : END DO
1288 52880 : DEALLOCATE (orb)
1289 : END DO
1290 : END DO
1291 :
1292 20518 : DEALLOCATE (nai, nbi, cpot)
1293 :
1294 20518 : END SUBROUTINE exchange_numeric
1295 :
1296 : ! **************************************************************************************************
1297 : !> \brief Calculate the exchange potential semi-analytically.
1298 : !> \param kmat Kohn-Sham matrix
1299 : !> \param state electronic state
1300 : !> \param occ occupation numbers
1301 : !> \param wfn wavefunctions
1302 : !> \param basis atomic basis set
1303 : !> \param hfx_pot properties of the Hartree-Fock potential
1304 : !> \par History
1305 : !> * 01.2009 created [Juerg Hutter]
1306 : ! **************************************************************************************************
1307 191 : SUBROUTINE exchange_semi_analytic(kmat, state, occ, wfn, basis, hfx_pot)
1308 : REAL(KIND=dp), DIMENSION(:, :, 0:), INTENT(INOUT) :: kmat
1309 : TYPE(atom_state), INTENT(IN) :: state
1310 : REAL(KIND=dp), DIMENSION(0:, :), INTENT(IN) :: occ
1311 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: wfn
1312 : TYPE(atom_basis_type), INTENT(IN) :: basis
1313 : TYPE(atom_hfx_type), INTENT(IN) :: hfx_pot
1314 :
1315 : INTEGER :: i, ia, ib, k, lad, lbc, lh, ll, nbas, &
1316 : norb, nr, nu
1317 : REAL(KIND=dp) :: almn
1318 191 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cpot, nai, nbi
1319 191 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: orb
1320 191 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: pot
1321 : REAL(KIND=dp), DIMENSION(0:maxfac) :: arho
1322 :
1323 191 : arho = 0._dp
1324 191 : DO ll = 0, maxfac, 2
1325 3056 : lh = ll/2
1326 3056 : arho(ll) = fac(ll)/fac(lh)**2
1327 : END DO
1328 :
1329 574097 : kmat = 0._dp
1330 :
1331 191 : nr = basis%grid%nr
1332 1337 : nbas = MAXVAL(basis%nbas)
1333 955 : ALLOCATE (pot(nr, nbas, nbas))
1334 1337 : ALLOCATE (nai(nr), nbi(nr), cpot(nr))
1335 :
1336 764 : DO lad = 0, state%maxl_calc
1337 1910 : DO lbc = 0, state%maxl_occ
1338 1146 : norb = state%maxn_occ(lbc)
1339 1146 : nbas = basis%nbas(lbc)
1340 : ! calculate orbitals for angmom lbc
1341 4584 : ALLOCATE (orb(nr, norb))
1342 445170 : orb = 0._dp
1343 2370 : DO i = 1, norb
1344 29058 : DO k = 1, nbas
1345 9127512 : orb(:, i) = orb(:, i) + wfn(k, i, lbc)*basis%bf(:, k, lbc)
1346 : END DO
1347 : END DO
1348 2674 : DO nu = ABS(lad - lbc), lad + lbc, 2
1349 : almn = arho(-lad + lbc + nu)*arho(lad - lbc + nu) &
1350 1528 : *arho(lad + lbc - nu)/(REAL(lad + lbc + nu + 1, dp)*arho(lad + lbc + nu))
1351 1528 : almn = -0.5_dp*almn
1352 : ! calculate potential for basis function pair (lad,lbc)
1353 242333208 : pot = 0._dp
1354 1528 : CALL potential_analytic(pot, lad, lbc, nu, basis, hfx_pot)
1355 34098 : DO ia = 1, basis%nbas(lad)
1356 66794 : DO i = 1, norb
1357 11818242 : cpot = 0._dp
1358 801328 : DO k = 1, nbas
1359 249602528 : cpot(:) = cpot(:) + pot(:, ia, k)*wfn(k, i, lbc)
1360 : END DO
1361 813096 : DO ib = 1, basis%nbas(lad)
1362 : kmat(ia, ib, lad) = kmat(ia, ib, lad) + almn*occ(lbc, i)* &
1363 781672 : integrate_grid(cpot, orb(:, i), basis%bf(:, ib, lad), basis%grid)
1364 : END DO
1365 : END DO
1366 : END DO
1367 : END DO
1368 1719 : DEALLOCATE (orb)
1369 : END DO
1370 : END DO
1371 :
1372 191 : DEALLOCATE (pot)
1373 191 : DEALLOCATE (nai, nbi, cpot)
1374 :
1375 191 : END SUBROUTINE exchange_semi_analytic
1376 :
1377 : ! **************************************************************************************************
1378 : !> \brief Calculate the electrostatic potential using numerical differentiation.
1379 : !> \param cpot computed potential
1380 : !> \param density electron density on the atomic radial grid
1381 : !> \param nu integer parameter [ABS(la-lb) .. la+lb];
1382 : !> see potential_analytic() and exchange_numeric()
1383 : !> \param grid atomic radial grid
1384 : !> \par History
1385 : !> * 01.2009 created [Juerg Hutter]
1386 : ! **************************************************************************************************
1387 662850 : SUBROUTINE potential_coulomb_numeric(cpot, density, nu, grid)
1388 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: cpot
1389 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: density
1390 : INTEGER, INTENT(IN) :: nu
1391 : TYPE(grid_atom_type), INTENT(IN) :: grid
1392 :
1393 : INTEGER :: i, nc
1394 : REAL(dp) :: int1, int2
1395 662850 : REAL(dp), DIMENSION(:), POINTER :: r, wr
1396 :
1397 662850 : nc = MIN(SIZE(cpot), SIZE(density))
1398 662850 : r => grid%rad
1399 662850 : wr => grid%wr
1400 :
1401 265802850 : int1 = integrate_grid(density, r**nu, grid)
1402 662850 : int2 = 0._dp
1403 1325700 : cpot(nc:) = int1/r(nc:)**(nu + 1)
1404 :
1405 : ! test that grid is decreasing
1406 662850 : CPASSERT(r(1) > r(nc))
1407 265802850 : DO i = 1, nc
1408 265140000 : cpot(i) = int1/r(i)**(nu + 1) + int2*r(i)**nu
1409 265140000 : int1 = int1 - r(i)**(nu)*density(i)*wr(i)
1410 265802850 : int2 = int2 + density(i)*wr(i)/r(i)**(nu + 1)
1411 : END DO
1412 :
1413 662850 : END SUBROUTINE potential_coulomb_numeric
1414 :
1415 : ! **************************************************************************************************
1416 : !> \brief ...
1417 : !> \param cpot ...
1418 : !> \param density ...
1419 : !> \param nu ...
1420 : !> \param grid ...
1421 : !> \param omega ...
1422 : !> \param kernel ...
1423 : ! **************************************************************************************************
1424 17486 : SUBROUTINE potential_longrange_numeric(cpot, density, nu, grid, omega, kernel)
1425 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: cpot
1426 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: density
1427 : INTEGER, INTENT(IN) :: nu
1428 : TYPE(grid_atom_type), INTENT(IN) :: grid
1429 : REAL(KIND=dp), INTENT(IN) :: omega
1430 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: kernel
1431 :
1432 : INTEGER :: nc
1433 17486 : REAL(dp), DIMENSION(:), POINTER :: r, wr
1434 17486 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: work_inp, work_out
1435 :
1436 17486 : nc = MIN(SIZE(cpot), SIZE(density))
1437 17486 : r => grid%rad
1438 17486 : wr => grid%wr
1439 :
1440 87430 : ALLOCATE (work_inp(nc), work_out(nc))
1441 :
1442 6068286 : cpot = 0.0_dp
1443 :
1444 : ! First Bessel transform
1445 6068286 : work_inp(:nc) = density(:nc)*wr(:nc)
1446 17486 : CALL DSYMV('U', nc, 1.0_dp, kernel, nc, work_inp, 1, 0.0_dp, work_out, 1)
1447 :
1448 : ! Second Bessel transform
1449 6068286 : work_inp(:nc) = work_out(:nc)*EXP(-r(:nc)**2)/r(:nc)**2*wr(:nc)
1450 17486 : CALL DSYMV('U', nc, 1.0_dp, kernel, nc, work_inp, 1, 0.0_dp, work_out, 1)
1451 :
1452 6068286 : cpot(:nc) = work_out(:nc)*(2.0_dp*REAL(nu, dp) + 1.0_dp)*4.0_dp/pi*omega
1453 :
1454 34972 : END SUBROUTINE potential_longrange_numeric
1455 :
1456 : ! **************************************************************************************************
1457 : !> \brief ...
1458 : !> \param cpot ...
1459 : !> \param density ...
1460 : !> \param nu ...
1461 : !> \param grid ...
1462 : !> \param omega ...
1463 : !> \param kernel ...
1464 : ! **************************************************************************************************
1465 17486 : SUBROUTINE potential_longrange_numeric_gh(cpot, density, nu, grid, omega, kernel)
1466 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: cpot
1467 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: density
1468 : INTEGER, INTENT(IN) :: nu
1469 : TYPE(grid_atom_type), INTENT(IN) :: grid
1470 : REAL(KIND=dp), INTENT(IN) :: omega
1471 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: kernel
1472 :
1473 : INTEGER :: n_max, nc, nc_kernel, nr_kernel
1474 17486 : REAL(dp), DIMENSION(:), POINTER :: wr
1475 17486 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: work_inp, work_out
1476 :
1477 17486 : nc = MIN(SIZE(cpot), SIZE(density))
1478 17486 : wr => grid%wr
1479 :
1480 17486 : nc_kernel = SIZE(kernel, 1)
1481 17486 : nr_kernel = SIZE(kernel, 2)
1482 17486 : n_max = MAX(nc, nc_kernel, nr_kernel)
1483 :
1484 87430 : ALLOCATE (work_inp(n_max), work_out(n_max))
1485 :
1486 6068286 : cpot = 0.0_dp
1487 :
1488 : ! First Bessel transform
1489 6068286 : work_inp(:nc) = density(:nc)*wr(:nc)
1490 17486 : CALL DGEMV('T', nc_kernel, nr_kernel, 1.0_dp, kernel, nc_kernel, work_inp, 1, 0.0_dp, work_out, 1)
1491 :
1492 : ! Second Bessel transform
1493 1766086 : work_inp(:nr_kernel) = work_out(:nr_kernel)
1494 17486 : CALL DGEMV('N', nc_kernel, nr_kernel, 1.0_dp, kernel, nc_kernel, work_inp, 1, 0.0_dp, work_out, 1)
1495 :
1496 6068286 : cpot(:nc) = work_out(:nc)*(2.0_dp*REAL(nu, dp) + 1.0_dp)*4.0_dp/pi*omega
1497 :
1498 34972 : END SUBROUTINE potential_longrange_numeric_gh
1499 :
1500 : ! **************************************************************************************************
1501 : !> \brief Calculate the electrostatic potential using analytical expressions.
1502 : !> The interaction potential has the form
1503 : !> V(r)=scale_coulomb*1/r+scale_lr*erf(omega*r)/r
1504 : !> \param cpot computed potential
1505 : !> \param la angular momentum of the calculated state
1506 : !> \param lb angular momentum of the occupied state
1507 : !> \param nu integer parameter [ABS(la-lb) .. la+lb] with the parity of 'la+lb'
1508 : !> \param basis atomic basis set
1509 : !> \param hfx_pot properties of the Hartree-Fock potential
1510 : !> \par History
1511 : !> * 01.2009 created [Juerg Hutter]
1512 : ! **************************************************************************************************
1513 1528 : SUBROUTINE potential_analytic(cpot, la, lb, nu, basis, hfx_pot)
1514 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: cpot
1515 : INTEGER, INTENT(IN) :: la, lb, nu
1516 : TYPE(atom_basis_type), INTENT(IN) :: basis
1517 : TYPE(atom_hfx_type), INTENT(IN) :: hfx_pot
1518 :
1519 : INTEGER :: i, j, k, l, ll, m
1520 : REAL(KIND=dp) :: a, b
1521 1528 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: erfa, pot
1522 :
1523 1528 : m = SIZE(cpot, 1)
1524 4584 : ALLOCATE (erfa(1:m))
1525 :
1526 1528 : ll = la + lb
1527 :
1528 242333208 : cpot = 0._dp
1529 :
1530 1528 : SELECT CASE (basis%basis_type)
1531 : CASE DEFAULT
1532 0 : CPABORT("")
1533 : CASE (GTO_BASIS)
1534 32952 : DO i = 1, basis%nbas(la)
1535 715808 : DO j = 1, basis%nbas(lb)
1536 682856 : a = basis%am(i, la)
1537 682856 : b = basis%am(j, lb)
1538 :
1539 682856 : IF (hfx_pot%scale_coulomb /= 0.0_dp) THEN
1540 329160 : CALL potential_coulomb_analytic(erfa, a, b, ll, nu, basis%grid%rad)
1541 :
1542 112946760 : cpot(:, i, j) = cpot(:, i, j) + erfa(:)*hfx_pot%scale_coulomb
1543 : END IF
1544 :
1545 714280 : IF (hfx_pot%scale_longrange /= 0.0_dp) THEN
1546 353696 : CALL potential_longrange_analytic(erfa, a, b, ll, nu, basis%grid%rad, hfx_pot%omega)
1547 :
1548 119611296 : cpot(:, i, j) = cpot(:, i, j) + erfa(:)*hfx_pot%scale_longrange
1549 : END IF
1550 : END DO
1551 : END DO
1552 : CASE (CGTO_BASIS)
1553 0 : ALLOCATE (pot(1:m))
1554 :
1555 1528 : DO i = 1, basis%nprim(la)
1556 0 : DO j = 1, basis%nprim(lb)
1557 0 : a = basis%am(i, la)
1558 0 : b = basis%am(j, lb)
1559 :
1560 0 : pot = 0.0_dp
1561 :
1562 0 : IF (hfx_pot%scale_coulomb /= 0.0_dp) THEN
1563 0 : CALL potential_coulomb_analytic(erfa, a, b, ll, nu, basis%grid%rad)
1564 :
1565 0 : pot(:) = pot(:) + erfa(:)*hfx_pot%scale_coulomb
1566 : END IF
1567 :
1568 0 : IF (hfx_pot%scale_longrange /= 0.0_dp) THEN
1569 0 : CALL potential_longrange_analytic(erfa, a, b, ll, nu, basis%grid%rad, hfx_pot%omega)
1570 :
1571 0 : pot(:) = pot(:) + erfa(:)*hfx_pot%scale_longrange
1572 : END IF
1573 :
1574 0 : DO k = 1, basis%nbas(la)
1575 0 : DO l = 1, basis%nbas(lb)
1576 0 : cpot(:, k, l) = cpot(:, k, l) + pot(:)*basis%cm(i, k, la)*basis%cm(j, l, lb)
1577 : END DO
1578 : END DO
1579 : END DO
1580 : END DO
1581 :
1582 : END SELECT
1583 :
1584 1528 : DEALLOCATE (erfa)
1585 :
1586 1528 : END SUBROUTINE potential_analytic
1587 :
1588 : ! **************************************************************************************************
1589 : !> \brief ...
1590 : !> \param erfa Array will contain the potential
1591 : !> \param a Exponent of first Gaussian charge distribution
1592 : !> \param b Exponent of second Gaussian charge distribution
1593 : !> \param ll Sum of angular momentum quantum numbers building one charge distribution
1594 : !> \param nu Angular momentum of interaction, (ll-nu) should be even.
1595 : !> \param rad Radial grid
1596 : ! **************************************************************************************************
1597 329160 : SUBROUTINE potential_coulomb_analytic(erfa, a, b, ll, nu, rad)
1598 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: erfa
1599 : REAL(KIND=dp), INTENT(IN) :: a, b
1600 : INTEGER, INTENT(IN) :: ll, nu
1601 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rad
1602 :
1603 : INTEGER :: nr
1604 : REAL(KIND=dp) :: oab, sab
1605 329160 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: expa, z
1606 :
1607 329160 : nr = SIZE(rad)
1608 1645800 : ALLOCATE (expa(nr), z(nr))
1609 :
1610 329160 : sab = SQRT(a + b)
1611 329160 : oab = dfac(ll + nu + 1)*rootpi/sab**(ll + 2)/2._dp**((ll + nu)/2 + 2)
1612 112946760 : z(:) = sab*rad(:)
1613 112946760 : erfa(:) = oab*erf(z(:))/z(:)**(nu + 1)
1614 112946760 : expa(:) = EXP(-z(:)**2)/sab**(ll + 2)/2._dp**((ll + nu)/2 + 2)
1615 0 : SELECT CASE (ll)
1616 : CASE DEFAULT
1617 0 : CPABORT("")
1618 : CASE (0)
1619 43941 : CPASSERT(nu == 0)
1620 : CASE (1)
1621 84522 : CPASSERT(nu == 1)
1622 28685322 : erfa(:) = erfa(:) - 6._dp*expa(:)/z(:)
1623 : CASE (2)
1624 78570 : SELECT CASE (nu)
1625 : CASE DEFAULT
1626 0 : CPABORT("")
1627 : CASE (0)
1628 14043573 : erfa(:) = erfa(:) - 2._dp*expa(:)
1629 : CASE (2)
1630 28089327 : erfa(:) = erfa(:) - expa(:)*(20._dp + 30._dp/z(:)**2)
1631 : END SELECT
1632 : CASE (3)
1633 0 : SELECT CASE (nu)
1634 : CASE DEFAULT
1635 0 : CPABORT("")
1636 : CASE (1)
1637 13744485 : erfa(:) = erfa(:) - expa(:)*(12._dp*z(:) + 30._dp/z(:))
1638 : CASE (3)
1639 13783770 : erfa(:) = erfa(:) - expa(:)*(56._dp*z(:) + 140._dp/z(:) + 210._dp/z(:)**3)
1640 : END SELECT
1641 : CASE (4)
1642 0 : SELECT CASE (nu)
1643 : CASE DEFAULT
1644 0 : CPABORT("")
1645 : CASE (0)
1646 0 : erfa(:) = erfa(:) - expa(:)*(4._dp*z(:)**2 + 14._dp)
1647 : CASE (2)
1648 0 : erfa(:) = erfa(:) - expa(:)*(40._dp*z(:)**2 + 140._dp + 210._dp/z(:)**2)
1649 : CASE (4)
1650 0 : erfa(:) = erfa(:) - expa(:)*(144._dp*z(:)**2 + 504._dp + 1260._dp/z(:)**2 + 1890._dp/z(:)**4)
1651 : END SELECT
1652 : CASE (5)
1653 0 : SELECT CASE (nu)
1654 : CASE DEFAULT
1655 0 : CPABORT("")
1656 : CASE (1)
1657 0 : erfa(:) = erfa(:) - expa(:)*(24._dp*z(:)**3 + 108._dp*z(:) + 210._dp/z(:))
1658 : CASE (3)
1659 0 : erfa(:) = erfa(:) - expa(:)*(112._dp*z(:)**3 + 504._dp*z(:) + 1260._dp/z(:) + 1890._dp/z(:)**3)
1660 : CASE (5)
1661 : erfa(:) = erfa(:) - expa(:)*(352._dp*z(:)**3 + 1584._dp*z(:) + 5544._dp/z(:) + &
1662 0 : 13860._dp/z(:)**3 + 20790._dp/z(:)**5)
1663 : END SELECT
1664 : CASE (6)
1665 329160 : SELECT CASE (nu)
1666 : CASE DEFAULT
1667 0 : CPABORT("")
1668 : CASE (0)
1669 0 : erfa(:) = erfa(:) - expa(:)*(8._dp*z(:)**4 + 44._dp*z(:)**2 + 114._dp)
1670 : CASE (2)
1671 0 : erfa(:) = erfa(:) - expa(:)*(80._dp*z(:)**4 + 440._dp*z(:)**2 + 1260._dp + 1896._dp/z(:)**2)
1672 : CASE (4)
1673 : erfa(:) = erfa(:) - expa(:)*(288._dp*z(:)**4 + 1584._dp*z(:)**2 + 5544._dp + &
1674 0 : 13860._dp/z(:)**2 + 20790._dp/z(:)**4)
1675 : CASE (6)
1676 : erfa(:) = erfa(:) - expa(:)*(832._dp*z(:)**4 + 4576._dp*z(:)**2 + 20592._dp + &
1677 0 : 72072._dp/z(:)**2 + 180180._dp/z(:)**4 + 270270._dp/z(:)**6)
1678 : END SELECT
1679 : END SELECT
1680 :
1681 329160 : DEALLOCATE (expa, z)
1682 :
1683 329160 : END SUBROUTINE potential_coulomb_analytic
1684 :
1685 : ! **************************************************************************************************
1686 : !> \brief This routine calculates the longrange Coulomb potential of a product of two Gaussian with given angular momentum
1687 : !> \param erfa Array will contain the potential
1688 : !> \param a Exponent of first Gaussian charge distribution
1689 : !> \param b Exponent of second Gaussian charge distribution
1690 : !> \param ll Sum of angular momentum quantum numbers building one charge distribution
1691 : !> \param nu Angular momentum of interaction, (ll-nu) should be even.
1692 : !> \param rad Radial grid
1693 : !> \param omega Range-separation parameter
1694 : ! **************************************************************************************************
1695 353696 : PURE SUBROUTINE potential_longrange_analytic(erfa, a, b, ll, nu, rad, omega)
1696 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: erfa
1697 : REAL(KIND=dp), INTENT(IN) :: a, b
1698 : INTEGER, INTENT(IN) :: ll, nu
1699 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rad
1700 : REAL(KIND=dp), INTENT(IN) :: omega
1701 :
1702 : INTEGER :: i, lambda, nr
1703 : REAL(KIND=dp) :: ab, oab, pab, prel, sab
1704 353696 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: expa, z
1705 :
1706 353696 : IF (MOD(ll - nu, 2) == 0 .AND. ll >= nu .AND. nu >= 0) THEN
1707 353696 : nr = SIZE(rad)
1708 1768480 : ALLOCATE (expa(nr), z(nr))
1709 :
1710 353696 : lambda = (ll - nu)/2
1711 353696 : ab = a + b
1712 353696 : sab = SQRT(ab)
1713 353696 : pab = omega**2*ab/(omega**2 + ab)
1714 353696 : prel = pab/ab
1715 119611296 : z(:) = SQRT(pab)*rad(:)
1716 353696 : oab = fac(lambda)/SQRT(ab)**(ll + 2)/4.0_dp*SQRT(prel)**(nu + 1)*(2.0_dp*REAL(nu, KIND=dp) + 1.0_dp)
1717 119611296 : expa(:) = EXP(-z(:)**2)
1718 : lambda = (ll - nu)/2
1719 :
1720 353696 : IF (lambda > 0) THEN
1721 29379420 : erfa = 0.0_dp
1722 171640 : DO i = 1, lambda
1723 : erfa = erfa + (-prel)**i/REAL(i, KIND=dp)*binomial_gen(lambda + nu + 0.5_dp, lambda - i)* &
1724 29465240 : assoc_laguerre(z, REAL(nu, KIND=dp) + 0.5_dp, i - 1)
1725 : END DO
1726 29379420 : erfa = erfa*expa*z**nu
1727 :
1728 29379420 : erfa = erfa + 2.0_dp*binomial_gen(lambda + nu + 0.5_dp, lambda)*znFn(z, expa, nu)
1729 : ELSE
1730 90231876 : erfa = 2.0_dp*znFn(z, expa, nu)
1731 : END IF
1732 :
1733 119611296 : erfa = erfa*oab
1734 :
1735 353696 : DEALLOCATE (expa, z)
1736 : ELSE
1737 : ! Contribution to potential is zero (properties of spherical harmonics)
1738 0 : erfa = 0.0_dp
1739 : END IF
1740 :
1741 353696 : END SUBROUTINE potential_longrange_analytic
1742 :
1743 : ! **************************************************************************************************
1744 : !> \brief Boys' function times z**n
1745 : !> \param z ...
1746 : !> \param expa ...
1747 : !> \param n ...
1748 : !> \return ...
1749 : ! **************************************************************************************************
1750 119257600 : ELEMENTAL FUNCTION znFn(z, expa, n) RESULT(res)
1751 :
1752 : REAL(KIND=dp), INTENT(IN) :: z, expa
1753 : INTEGER, INTENT(IN) :: n
1754 : REAL(KIND=dp) :: res
1755 :
1756 : INTEGER :: i
1757 : REAL(KIND=dp) :: z_exp
1758 :
1759 119257600 : IF (n >= 0) THEN
1760 119257600 : IF (ABS(z) < 1.0E-20) THEN
1761 0 : res = z**n/(2.0_dp*REAL(n, KIND=dp) + 1.0_dp)
1762 119257600 : ELSE IF (n == 0) THEN
1763 30380000 : res = rootpi/2.0_dp*ERF(z)/z
1764 : ELSE
1765 88877600 : res = rootpi/4.0_dp*ERF(z)/z**2 - expa/z/2.0_dp
1766 88877600 : z_exp = -expa*0.5_dp
1767 :
1768 147420000 : DO i = 2, n
1769 58542400 : res = (REAL(i, KIND=dp) - 0.5_dp)*res/z + z_exp
1770 147420000 : z_exp = z_exp*z
1771 : END DO
1772 : END IF
1773 : ELSE ! Set it to zero (no Boys' function, to keep the ELEMENTAL keyword)
1774 : res = 0.0_dp
1775 : END IF
1776 :
1777 119257600 : END FUNCTION znFn
1778 :
1779 : ! **************************************************************************************************
1780 : !> \brief ...
1781 : !> \param z ...
1782 : !> \param a ...
1783 : !> \param n ...
1784 : !> \return ...
1785 : ! **************************************************************************************************
1786 29293600 : ELEMENTAL FUNCTION assoc_laguerre(z, a, n) RESULT(res)
1787 :
1788 : REAL(KIND=dp), INTENT(IN) :: z, a
1789 : INTEGER, INTENT(IN) :: n
1790 : REAL(KIND=dp) :: res
1791 :
1792 : INTEGER :: i
1793 : REAL(KIND=dp) :: f0, f1
1794 :
1795 29293600 : IF (n == 0) THEN
1796 : res = 1.0_dp
1797 0 : ELSE IF (n == 1) THEN
1798 0 : res = a + 1.0_dp - z
1799 0 : ELSE IF (n > 0) THEN
1800 0 : f0 = 1.0_dp
1801 0 : f1 = a + 1.0_dp - z
1802 :
1803 0 : DO i = 3, n
1804 0 : res = (2.0_dp + (a - 1.0_dp - z)/REAL(i, dp))*f1 - (1.0_dp + (a - 1.0_dp)/REAL(i, dp))*f0
1805 0 : f0 = f1
1806 0 : f1 = res
1807 : END DO
1808 : ELSE ! n is negative, set it zero (no polynomials, to keep the ELEMENTAL keyword)
1809 : res = 0.0_dp
1810 : END IF
1811 :
1812 29293600 : END FUNCTION assoc_laguerre
1813 :
1814 : ! **************************************************************************************************
1815 : !> \brief Compute Trace[opmat * pmat] == Trace[opmat * pmat^T] .
1816 : !> \param opmat operator matrix (e.g. Kohn-Sham matrix or overlap matrix)
1817 : !> \param pmat density matrix
1818 : !> \return value of trace
1819 : !> \par History
1820 : !> * 08.2008 created [Juerg Hutter]
1821 : ! **************************************************************************************************
1822 406298 : PURE FUNCTION atom_trace(opmat, pmat) RESULT(trace)
1823 : REAL(KIND=dp), DIMENSION(:, :, 0:), INTENT(IN) :: opmat, pmat
1824 : REAL(KIND=dp) :: trace
1825 :
1826 406298 : trace = accurate_dot_product(opmat, pmat)
1827 :
1828 406298 : END FUNCTION atom_trace
1829 :
1830 : ! **************************************************************************************************
1831 : !> \brief Calculate a potential matrix by integrating the potential on an atomic radial grid.
1832 : !> \param imat potential matrix
1833 : !> \param cpot potential on the atomic radial grid
1834 : !> \param basis atomic basis set
1835 : !> \param derivatives order of derivatives
1836 : !> \par History
1837 : !> * 08.2008 created [Juerg Hutter]
1838 : ! **************************************************************************************************
1839 189210 : SUBROUTINE numpot_matrix(imat, cpot, basis, derivatives)
1840 : REAL(KIND=dp), DIMENSION(:, :, 0:), INTENT(INOUT) :: imat
1841 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: cpot
1842 : TYPE(atom_basis_type), INTENT(INOUT) :: basis
1843 : INTEGER, INTENT(IN) :: derivatives
1844 :
1845 : INTEGER :: i, j, l, n
1846 :
1847 189210 : SELECT CASE (derivatives)
1848 : CASE (0)
1849 1153698 : DO l = 0, lmat
1850 988884 : n = basis%nbas(l)
1851 3628630 : DO i = 1, n
1852 29769297 : DO j = i, n
1853 : imat(i, j, l) = imat(i, j, l) + &
1854 26305481 : integrate_grid(cpot, basis%bf(:, i, l), basis%bf(:, j, l), basis%grid)
1855 28780413 : imat(j, i, l) = imat(i, j, l)
1856 : END DO
1857 : END DO
1858 : END DO
1859 : CASE (1)
1860 169694 : DO l = 0, lmat
1861 145452 : n = basis%nbas(l)
1862 1207081 : DO i = 1, n
1863 14942611 : DO j = i, n
1864 : imat(i, j, l) = imat(i, j, l) + &
1865 13759772 : integrate_grid(cpot, basis%dbf(:, i, l), basis%bf(:, j, l), basis%grid)
1866 : imat(i, j, l) = imat(i, j, l) + &
1867 13759772 : integrate_grid(cpot, basis%bf(:, i, l), basis%dbf(:, j, l), basis%grid)
1868 14797159 : imat(j, i, l) = imat(i, j, l)
1869 : END DO
1870 : END DO
1871 : END DO
1872 : CASE (2)
1873 1078 : DO l = 0, lmat
1874 924 : n = basis%nbas(l)
1875 24234 : DO i = 1, n
1876 459894 : DO j = i, n
1877 : imat(i, j, l) = imat(i, j, l) + &
1878 435814 : integrate_grid(cpot, basis%dbf(:, i, l), basis%dbf(:, j, l), basis%grid)
1879 458970 : imat(j, i, l) = imat(i, j, l)
1880 : END DO
1881 : END DO
1882 : END DO
1883 : CASE DEFAULT
1884 189210 : CPABORT("")
1885 : END SELECT
1886 :
1887 189210 : END SUBROUTINE numpot_matrix
1888 :
1889 : ! **************************************************************************************************
1890 : !> \brief Contract Coulomb Electron Repulsion Integrals.
1891 : !> \param jmat ...
1892 : !> \param erint ...
1893 : !> \param pmat ...
1894 : !> \param nsize ...
1895 : !> \param all_nu ...
1896 : !> \par History
1897 : !> * 08.2008 created [Juerg Hutter]
1898 : ! **************************************************************************************************
1899 1515 : PURE SUBROUTINE ceri_contract(jmat, erint, pmat, nsize, all_nu)
1900 : REAL(KIND=dp), DIMENSION(:, :, 0:), INTENT(INOUT) :: jmat
1901 : TYPE(eri), DIMENSION(:), INTENT(IN) :: erint
1902 : REAL(KIND=dp), DIMENSION(:, :, 0:), INTENT(IN) :: pmat
1903 : INTEGER, DIMENSION(0:), INTENT(IN) :: nsize
1904 : LOGICAL, INTENT(IN), OPTIONAL :: all_nu
1905 :
1906 : INTEGER :: i1, i2, ij1, ij2, j1, j2, l1, l2, ll, &
1907 : n1, n2
1908 : LOGICAL :: have_all_nu
1909 : REAL(KIND=dp) :: eint, f1, f2
1910 :
1911 1515 : IF (PRESENT(all_nu)) THEN
1912 0 : have_all_nu = all_nu
1913 : ELSE
1914 : have_all_nu = .FALSE.
1915 : END IF
1916 :
1917 4237113 : jmat(:, :, :) = 0._dp
1918 : ll = 0
1919 10605 : DO l1 = 0, lmat
1920 9090 : n1 = nsize(l1)
1921 42420 : DO l2 = 0, l1
1922 31815 : n2 = nsize(l2)
1923 31815 : ll = ll + 1
1924 31815 : ij1 = 0
1925 558425 : DO i1 = 1, n1
1926 6043326 : DO j1 = i1, n1
1927 5484901 : ij1 = ij1 + 1
1928 5484901 : f1 = 1._dp
1929 5484901 : IF (i1 /= j1) f1 = 2._dp
1930 5484901 : ij2 = 0
1931 119694873 : DO i2 = 1, n2
1932 1403690542 : DO j2 = i2, n2
1933 1284522279 : ij2 = ij2 + 1
1934 1284522279 : f2 = 1._dp
1935 1284522279 : IF (i2 /= j2) f2 = 2._dp
1936 1284522279 : eint = erint(ll)%int(ij1, ij2)
1937 1398205641 : IF (l1 == l2) THEN
1938 418306136 : jmat(i1, j1, l1) = jmat(i1, j1, l1) + f2*pmat(i2, j2, l2)*eint
1939 : ELSE
1940 866216143 : jmat(i1, j1, l1) = jmat(i1, j1, l1) + f2*pmat(i2, j2, l2)*eint
1941 866216143 : jmat(i2, j2, l2) = jmat(i2, j2, l2) + f1*pmat(i1, j1, l1)*eint
1942 : END IF
1943 : END DO
1944 : END DO
1945 : END DO
1946 : END DO
1947 40905 : IF (have_all_nu) THEN
1948 : ! skip integral blocks with nu/=0
1949 0 : ll = ll + l2
1950 : END IF
1951 : END DO
1952 : END DO
1953 10605 : DO l1 = 0, lmat
1954 9090 : n1 = nsize(l1)
1955 170209 : DO i1 = 1, n1
1956 1873302 : DO j1 = i1, n1
1957 1864212 : jmat(j1, i1, l1) = jmat(i1, j1, l1)
1958 : END DO
1959 : END DO
1960 : END DO
1961 :
1962 1515 : END SUBROUTINE ceri_contract
1963 :
1964 : ! **************************************************************************************************
1965 : !> \brief Contract exchange Electron Repulsion Integrals.
1966 : !> \param kmat ...
1967 : !> \param erint ...
1968 : !> \param pmat ...
1969 : !> \param nsize ...
1970 : !> \par History
1971 : !> * 08.2008 created [Juerg Hutter]
1972 : ! **************************************************************************************************
1973 547 : PURE SUBROUTINE eeri_contract(kmat, erint, pmat, nsize)
1974 : REAL(KIND=dp), DIMENSION(:, :, 0:), INTENT(INOUT) :: kmat
1975 : TYPE(eri), DIMENSION(:), INTENT(IN) :: erint
1976 : REAL(KIND=dp), DIMENSION(:, :, 0:), INTENT(IN) :: pmat
1977 : INTEGER, DIMENSION(0:), INTENT(IN) :: nsize
1978 :
1979 : INTEGER :: i1, i2, ij1, ij2, j1, j2, l1, l2, lh, &
1980 : ll, n1, n2, nu
1981 : REAL(KIND=dp) :: almn, eint, f1, f2
1982 : REAL(KIND=dp), DIMENSION(0:maxfac) :: arho
1983 :
1984 547 : arho = 0._dp
1985 547 : DO ll = 0, maxfac, 2
1986 8752 : lh = ll/2
1987 8752 : arho(ll) = fac(ll)/fac(lh)**2
1988 : END DO
1989 :
1990 1470517 : kmat(:, :, :) = 0._dp
1991 : ll = 0
1992 3829 : DO l1 = 0, lmat
1993 3282 : n1 = nsize(l1)
1994 15316 : DO l2 = 0, l1
1995 11487 : n2 = nsize(l2)
1996 45401 : DO nu = ABS(l1 - l2), l1 + l2, 2
1997 30632 : almn = arho(-l1 + l2 + nu)*arho(l1 - l2 + nu)*arho(l1 + l2 - nu)/(REAL(l1 + l2 + nu + 1, dp)*arho(l1 + l2 + nu))
1998 30632 : almn = -0.5_dp*almn
1999 30632 : ll = ll + 1
2000 30632 : ij1 = 0
2001 504257 : DO i1 = 1, n1
2002 5332602 : DO j1 = i1, n1
2003 4839832 : ij1 = ij1 + 1
2004 4839832 : f1 = 1._dp
2005 4839832 : IF (i1 /= j1) f1 = 2._dp
2006 4839832 : ij2 = 0
2007 104653578 : DO i2 = 1, n2
2008 1197184274 : DO j2 = i2, n2
2009 1092992834 : ij2 = ij2 + 1
2010 1092992834 : f2 = 1._dp
2011 1092992834 : IF (i2 /= j2) f2 = 2._dp
2012 1092992834 : eint = erint(ll)%int(ij1, ij2)
2013 1192344442 : IF (l1 == l2) THEN
2014 430803648 : kmat(i1, j1, l1) = kmat(i1, j1, l1) + f2*almn*pmat(i2, j2, l2)*eint
2015 : ELSE
2016 662189186 : kmat(i1, j1, l1) = kmat(i1, j1, l1) + f2*almn*pmat(i2, j2, l2)*eint
2017 662189186 : kmat(i2, j2, l2) = kmat(i2, j2, l2) + f1*almn*pmat(i1, j1, l1)*eint
2018 : END IF
2019 : END DO
2020 : END DO
2021 : END DO
2022 : END DO
2023 : END DO
2024 : END DO
2025 : END DO
2026 3829 : DO l1 = 0, lmat
2027 3282 : n1 = nsize(l1)
2028 58734 : DO i1 = 1, n1
2029 647742 : DO j1 = i1, n1
2030 644460 : kmat(j1, i1, l1) = kmat(i1, j1, l1)
2031 : END DO
2032 : END DO
2033 : END DO
2034 :
2035 547 : END SUBROUTINE eeri_contract
2036 :
2037 : ! **************************************************************************************************
2038 : !> \brief Calculate the error matrix for each angular momentum.
2039 : !> \param emat error matrix
2040 : !> \param demax the largest absolute value of error matrix elements
2041 : !> \param kmat Kohn-Sham matrix
2042 : !> \param pmat electron density matrix
2043 : !> \param umat transformation matrix which reduce overlap matrix to its unitary form
2044 : !> \param upmat transformation matrix which reduce density matrix to its unitary form
2045 : !> \param nval number of linear-independent contracted basis functions
2046 : !> \param nbs number of contracted basis functions
2047 : !> \par History
2048 : !> * 08.2008 created [Juerg Hutter]
2049 : ! **************************************************************************************************
2050 78874 : PURE SUBROUTINE err_matrix(emat, demax, kmat, pmat, umat, upmat, nval, nbs)
2051 : REAL(KIND=dp), DIMENSION(:, :, 0:), INTENT(OUT) :: emat
2052 : REAL(KIND=dp), INTENT(OUT) :: demax
2053 : REAL(KIND=dp), DIMENSION(:, :, 0:), INTENT(IN) :: kmat, pmat, umat, upmat
2054 : INTEGER, DIMENSION(0:), INTENT(IN) :: nval, nbs
2055 :
2056 : INTEGER :: l, m, n
2057 78874 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: tkmat, tpmat
2058 :
2059 37831654 : emat = 0._dp
2060 552118 : DO l = 0, lmat
2061 473244 : n = nval(l)
2062 473244 : m = nbs(l)
2063 552118 : IF (m > 0) THEN
2064 1355550 : ALLOCATE (tkmat(1:m, 1:m), tpmat(1:m, 1:m))
2065 13554160 : tkmat = 0._dp
2066 13554160 : tpmat = 0._dp
2067 936374169 : tkmat(1:m, 1:m) = MATMUL(TRANSPOSE(umat(1:n, 1:m, l)), MATMUL(kmat(1:n, 1:n, l), umat(1:n, 1:m, l)))
2068 936374169 : tpmat(1:m, 1:m) = MATMUL(TRANSPOSE(umat(1:n, 1:m, l)), MATMUL(pmat(1:n, 1:n, l), umat(1:n, 1:m, l)))
2069 689670701 : tpmat(1:m, 1:m) = MATMUL(upmat(1:m, 1:m, l), MATMUL(tpmat(1:m, 1:m), upmat(1:m, 1:m, l)))
2070 :
2071 382416872 : emat(1:m, 1:m, l) = MATMUL(tkmat(1:m, 1:m), tpmat(1:m, 1:m)) - MATMUL(tpmat(1:m, 1:m), tkmat(1:m, 1:m))
2072 :
2073 193650 : DEALLOCATE (tkmat, tpmat)
2074 : END IF
2075 : END DO
2076 37831654 : demax = MAXVAL(ABS(emat))
2077 :
2078 78874 : END SUBROUTINE err_matrix
2079 :
2080 : ! **************************************************************************************************
2081 : !> \brief Calculate Slater density on a radial grid.
2082 : !> \param density1 alpha-spin electron density
2083 : !> \param density2 beta-spin electron density
2084 : !> \param zcore nuclear charge
2085 : !> \param state electronic state
2086 : !> \param grid atomic radial grid
2087 : !> \par History
2088 : !> * 06.2018 bugfix [Rustam Khaliullin]
2089 : !> * 02.2010 unrestricted KS and HF methods [Juerg Hutter]
2090 : !> * 12.2008 created [Juerg Hutter]
2091 : !> \note An initial electron density to guess atomic orbitals.
2092 : ! **************************************************************************************************
2093 10850 : SUBROUTINE slater_density(density1, density2, zcore, state, grid)
2094 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: density1, density2
2095 : INTEGER, INTENT(IN) :: zcore
2096 : TYPE(atom_state), INTENT(IN) :: state
2097 : TYPE(grid_atom_type), INTENT(IN) :: grid
2098 :
2099 : INTEGER :: counter, i, l, mc, mm(0:lmat), mo, n, ns
2100 : INTEGER, DIMENSION(lmat+1, 20) :: ne
2101 : REAL(KIND=dp) :: a, pf
2102 :
2103 : ! fill out a table with the total number of electrons
2104 : ! core + valence. format ne(l,n)
2105 10850 : ns = SIZE(ne, 2)
2106 10850 : ne = 0
2107 10850 : mm = get_maxn_occ(state%core)
2108 75950 : DO l = 0, lmat
2109 65100 : mo = state%maxn_occ(l)
2110 726950 : IF (SUM(state%core(l, :)) == 0) THEN
2111 57763 : CPASSERT(ns >= l + mo)
2112 71865 : DO counter = 1, mo
2113 71865 : ne(l + 1, l + counter) = NINT(state%occ(l, counter))
2114 : END DO
2115 : ELSE
2116 7337 : mc = mm(l) ! number of levels in the core
2117 19800 : CPASSERT(SUM(state%occ(l, 1:mc)) == 0)
2118 7337 : CPASSERT(ns >= l + mc)
2119 19800 : DO counter = 1, mc
2120 19800 : ne(l + 1, l + counter) = NINT(state%core(l, counter))
2121 : END DO
2122 7337 : CPASSERT(ns >= l + mc + mo)
2123 13762 : DO counter = mc + 1, mc + mo
2124 13762 : ne(l + 1, l + counter) = NINT(state%occ(l, counter))
2125 : END DO
2126 : END IF
2127 : END DO
2128 :
2129 4338958 : density1 = 0._dp
2130 4338958 : density2 = 0._dp
2131 28993 : DO l = 0, state%maxl_occ
2132 210423 : DO i = 1, SIZE(state%occ, 2)
2133 199573 : IF (state%occ(l, i) > 0._dp) THEN
2134 20887 : n = i + l
2135 20887 : a = srules(zcore, ne, n, l)
2136 20887 : pf = 1._dp/SQRT(fac(2*n))*(2._dp*a)**(n + 0.5_dp)
2137 20887 : IF (state%multiplicity == -1) THEN
2138 8260935 : density1(:) = density1(:) + state%occ(l, i)/fourpi*(grid%rad(:)**(n - 1)*EXP(-a*grid%rad(:))*pf)**2
2139 : ELSE
2140 75786 : density1(:) = density1(:) + state%occa(l, i)/fourpi*(grid%rad(:)**(n - 1)*EXP(-a*grid%rad(:))*pf)**2
2141 75786 : density2(:) = density2(:) + state%occb(l, i)/fourpi*(grid%rad(:)**(n - 1)*EXP(-a*grid%rad(:))*pf)**2
2142 : END IF
2143 : END IF
2144 : END DO
2145 : END DO
2146 :
2147 10850 : END SUBROUTINE slater_density
2148 :
2149 : ! **************************************************************************************************
2150 : !> \brief Calculate the functional derivative of the Wigner (correlation) - Slater (exchange)
2151 : !> density functional.
2152 : !> \param rho electron density on a radial grid
2153 : !> \param vxc (output) exchange-correlation potential
2154 : !> \par History
2155 : !> * 12.2008 created [Juerg Hutter]
2156 : !> \note A model XC-potential to guess atomic orbitals.
2157 : ! **************************************************************************************************
2158 10946 : PURE SUBROUTINE wigner_slater_functional(rho, vxc)
2159 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rho
2160 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: vxc
2161 :
2162 : INTEGER :: i
2163 : REAL(KIND=dp) :: ec, ex, rs, vc, vx
2164 :
2165 4378254 : vxc = 0._dp
2166 4378254 : DO i = 1, SIZE(rho)
2167 4378254 : IF (rho(i) > 1.e-20_dp) THEN
2168 : ! 3/4 * (3/pi)^{1/3} == 0.7385588
2169 4092531 : ex = -0.7385588_dp*rho(i)**0.333333333_dp
2170 4092531 : vx = 1.333333333_dp*ex
2171 4092531 : rs = (3._dp/fourpi/rho(i))**0.333333333_dp
2172 4092531 : ec = -0.88_dp/(rs + 7.8_dp)
2173 4092531 : vc = ec*(1._dp + rs/(3._dp*(rs + 7.8_dp)))
2174 4092531 : vxc(i) = vx + vc
2175 : END IF
2176 : END DO
2177 :
2178 10946 : END SUBROUTINE wigner_slater_functional
2179 :
2180 : ! **************************************************************************************************
2181 : !> \brief Check that the atomic multiplicity is consistent with the electronic structure method.
2182 : !> \param method electronic structure method
2183 : !> \param multiplicity multiplicity
2184 : !> \return consistency status
2185 : !> \par History
2186 : !> * 11.2009 unrestricted KS and HF methods [Juerg Hutter]
2187 : ! **************************************************************************************************
2188 3106 : PURE FUNCTION atom_consistent_method(method, multiplicity) RESULT(consistent)
2189 : INTEGER, INTENT(IN) :: method, multiplicity
2190 : LOGICAL :: consistent
2191 :
2192 : ! multiplicity == -1 means it has not been specified explicitly;
2193 : ! see the source code of the subroutine atom_set_occupation() for further details.
2194 3106 : SELECT CASE (method)
2195 : CASE DEFAULT
2196 1737 : consistent = .FALSE.
2197 : CASE (do_rks_atom)
2198 1737 : consistent = (multiplicity == -1)
2199 : CASE (do_uks_atom)
2200 93 : consistent = (multiplicity /= -1)
2201 : CASE (do_rhf_atom)
2202 1268 : consistent = (multiplicity == -1)
2203 : CASE (do_uhf_atom)
2204 8 : consistent = (multiplicity /= -1)
2205 : CASE (do_rohf_atom)
2206 3106 : consistent = .FALSE.
2207 : END SELECT
2208 :
2209 3106 : END FUNCTION atom_consistent_method
2210 :
2211 : ! **************************************************************************************************
2212 : !> \brief Calculate the total electron density at R=0.
2213 : !> \param atom information about the atomic kind
2214 : !> \param rho0 (output) computed density
2215 : !> \par History
2216 : !> * 05.2016 created [Juerg Hutter]
2217 : ! **************************************************************************************************
2218 2176 : SUBROUTINE get_rho0(atom, rho0)
2219 : TYPE(atom_type), INTENT(IN) :: atom
2220 : REAL(KIND=dp), INTENT(OUT) :: rho0
2221 :
2222 : INTEGER :: m0, m1, m2, method, nr
2223 : LOGICAL :: nlcc, spinpol
2224 : REAL(KIND=dp) :: d0, d1, d2, r0, r1, r2, w0, w1
2225 2176 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: xfun
2226 2176 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: rho
2227 :
2228 2176 : method = atom%method_type
2229 : SELECT CASE (method)
2230 : CASE (do_rks_atom, do_rhf_atom)
2231 45 : spinpol = .FALSE.
2232 : CASE (do_uks_atom, do_uhf_atom)
2233 45 : spinpol = .TRUE.
2234 : CASE (do_rohf_atom)
2235 0 : CPABORT("")
2236 : CASE DEFAULT
2237 2176 : CPABORT("")
2238 : END SELECT
2239 :
2240 2176 : nr = atom%basis%grid%nr
2241 2176 : nlcc = .FALSE.
2242 2176 : IF (atom%potential%ppot_type == gth_pseudo) THEN
2243 1771 : nlcc = atom%potential%gth_pot%nlcc
2244 405 : ELSE IF (atom%potential%ppot_type == upf_pseudo) THEN
2245 2 : nlcc = atom%potential%upf_pot%core_correction
2246 403 : ELSE IF (atom%potential%ppot_type == sgp_pseudo) THEN
2247 0 : nlcc = atom%potential%sgp_pot%has_nlcc
2248 : END IF
2249 1773 : IF (nlcc) THEN
2250 24 : ALLOCATE (xfun(nr))
2251 : END IF
2252 :
2253 876208 : m0 = MAXVAL(MINLOC(atom%basis%grid%rad))
2254 2176 : IF (m0 == nr) THEN
2255 2176 : m1 = m0 - 1
2256 2176 : m2 = m0 - 2
2257 0 : ELSE IF (m0 == 1) THEN
2258 : m1 = 2
2259 : m2 = 3
2260 : ELSE
2261 0 : CPABORT("GRID Definition incompatible")
2262 : END IF
2263 2176 : r0 = atom%basis%grid%rad(m0)
2264 2176 : r1 = atom%basis%grid%rad(m1)
2265 2176 : r2 = atom%basis%grid%rad(m2)
2266 2176 : w0 = r1/(r1 - r0)
2267 2176 : w1 = 1 - w0
2268 :
2269 2176 : IF (spinpol) THEN
2270 135 : ALLOCATE (rho(nr, 2))
2271 45 : CALL atom_density(rho(:, 1), atom%orbitals%pmata, atom%basis, atom%state%maxl_occ, typ="RHO")
2272 45 : CALL atom_density(rho(:, 2), atom%orbitals%pmatb, atom%basis, atom%state%maxl_occ, typ="RHO")
2273 45 : IF (nlcc) THEN
2274 401 : xfun = 0.0_dp
2275 1 : CALL atom_core_density(xfun(:), atom%potential, typ="RHO", rr=atom%basis%grid%rad)
2276 401 : rho(:, 1) = rho(:, 1) + 0.5_dp*xfun(:)
2277 401 : rho(:, 2) = rho(:, 2) + 0.5_dp*xfun(:)
2278 : END IF
2279 18445 : rho(:, 1) = rho(:, 1) + rho(:, 2)
2280 : ELSE
2281 6393 : ALLOCATE (rho(nr, 1))
2282 2131 : CALL atom_density(rho(:, 1), atom%orbitals%pmat, atom%basis, atom%state%maxl_occ, typ="RHO")
2283 2131 : IF (nlcc) THEN
2284 7 : CALL atom_core_density(rho(:, 1), atom%potential, typ="RHO", rr=atom%basis%grid%rad)
2285 : END IF
2286 : END IF
2287 2176 : d0 = rho(m0, 1)
2288 2176 : d1 = rho(m1, 1)
2289 2176 : d2 = rho(m2, 1)
2290 :
2291 2176 : rho0 = w0*d0 + w1*d1
2292 2176 : rho0 = MAX(rho0, 0.0_dp)
2293 :
2294 2176 : DEALLOCATE (rho)
2295 2176 : IF (nlcc) THEN
2296 8 : DEALLOCATE (xfun)
2297 : END IF
2298 :
2299 2176 : END SUBROUTINE get_rho0
2300 :
2301 : ! **************************************************************************************************
2302 : !> \brief Print condition numbers of the given atomic basis set.
2303 : !> \param basis atomic basis set
2304 : !> \param crad atomic covalent radius
2305 : !> \param iw output file unit
2306 : !> \par History
2307 : !> * 05.2016 created [Juerg Hutter]
2308 : ! **************************************************************************************************
2309 5 : SUBROUTINE atom_condnumber(basis, crad, iw)
2310 : TYPE(atom_basis_type), POINTER :: basis
2311 : REAL(KIND=dp) :: crad
2312 : INTEGER, INTENT(IN) :: iw
2313 :
2314 : INTEGER :: i
2315 : REAL(KIND=dp) :: ci
2316 : REAL(KIND=dp), DIMENSION(10) :: cnum, rad
2317 :
2318 5 : WRITE (iw, '(/,A,F8.4)') " Basis Set Condition Numbers: 2*covalent radius=", 2*crad
2319 5 : CALL init_orbital_pointers(lmat)
2320 5 : CALL init_spherical_harmonics(lmat, 0)
2321 5 : cnum = 0.0_dp
2322 50 : DO i = 1, 9
2323 45 : ci = 2.0_dp*(0.85_dp + i*0.05_dp)
2324 45 : rad(i) = crad*ci
2325 45 : CALL atom_basis_condnum(basis, rad(i), cnum(i))
2326 45 : WRITE (iw, '(A,F15.3,T50,A,F14.4)') " Lattice constant:", &
2327 95 : rad(i), "Condition number:", cnum(i)
2328 : END DO
2329 5 : rad(10) = 0.01_dp
2330 5 : CALL atom_basis_condnum(basis, rad(10), cnum(10))
2331 5 : WRITE (iw, '(A,A,T50,A,F14.4)') " Lattice constant:", &
2332 10 : " Inf", "Condition number:", cnum(i)
2333 5 : CALL deallocate_orbital_pointers
2334 5 : CALL deallocate_spherical_harmonics
2335 :
2336 5 : END SUBROUTINE atom_condnumber
2337 :
2338 : ! **************************************************************************************************
2339 : !> \brief Estimate completeness of the given atomic basis set.
2340 : !> \param basis atomic basis set
2341 : !> \param zv atomic number
2342 : !> \param iw output file unit
2343 : ! **************************************************************************************************
2344 5 : SUBROUTINE atom_completeness(basis, zv, iw)
2345 : TYPE(atom_basis_type), POINTER :: basis
2346 : INTEGER, INTENT(IN) :: zv, iw
2347 :
2348 : INTEGER :: i, j, l, ll, m, n, nbas, nl, nr
2349 : INTEGER, DIMENSION(0:lmat) :: nelem, nlmax, nlmin
2350 : INTEGER, DIMENSION(lmat+1, 7) :: ne
2351 : REAL(KIND=dp) :: al, c1, c2, pf
2352 5 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: sfun
2353 5 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: bmat
2354 5 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: omat
2355 5 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: sint
2356 : REAL(KIND=dp), DIMENSION(0:lmat, 10) :: snl
2357 : REAL(KIND=dp), DIMENSION(2) :: sse
2358 : REAL(KIND=dp), DIMENSION(lmat+1, 7) :: sexp
2359 :
2360 5 : ne = 0
2361 5 : nelem = 0
2362 25 : nelem(0:3) = ptable(zv)%e_conv(0:3)
2363 35 : DO l = 0, lmat
2364 30 : ll = 2*(2*l + 1)
2365 64 : DO i = 1, 7
2366 89 : IF (nelem(l) >= ll) THEN
2367 22 : ne(l + 1, i) = ll
2368 22 : nelem(l) = nelem(l) - ll
2369 37 : ELSE IF (nelem(l) > 0) THEN
2370 7 : ne(l + 1, i) = nelem(l)
2371 7 : nelem(l) = 0
2372 : ELSE
2373 : EXIT
2374 : END IF
2375 : END DO
2376 : END DO
2377 :
2378 35 : nlmin = 1
2379 35 : nlmax = 1
2380 35 : DO l = 0, lmat
2381 30 : nlmin(l) = l + 1
2382 240 : DO i = 1, 7
2383 240 : IF (ne(l + 1, i) > 0) THEN
2384 29 : nlmax(l) = i + l
2385 : END IF
2386 : END DO
2387 35 : nlmax(l) = MAX(nlmax(l), nlmin(l) + 1)
2388 : END DO
2389 :
2390 : ! Slater exponents
2391 : sexp = 0.0_dp
2392 35 : DO l = 0, lmat
2393 30 : sse(1) = 0.05_dp
2394 30 : sse(2) = 10.0_dp
2395 165 : DO i = l + 1, 7
2396 135 : sexp(l + 1, i) = srules(zv, ne, i, l)
2397 165 : IF (ne(l + 1, i - l) > 0) THEN
2398 29 : sse(1) = MAX(sse(1), sexp(l + 1, i))
2399 29 : sse(2) = MIN(sse(2), sexp(l + 1, i))
2400 : END IF
2401 : END DO
2402 335 : DO i = 1, 10
2403 330 : snl(l, i) = ABS(2._dp*sse(1) - 0.5_dp*sse(2))/9._dp*REAL(i - 1, KIND=dp) + 0.5_dp*MIN(sse(1), sse(2))
2404 : END DO
2405 : END DO
2406 :
2407 35 : nbas = MAXVAL(basis%nbas)
2408 25 : ALLOCATE (omat(nbas, nbas, 0:lmat))
2409 5 : nr = SIZE(basis%bf, 1)
2410 30 : ALLOCATE (sfun(nr), sint(10, 2, nbas, 0:lmat))
2411 8453 : sint = 0._dp
2412 :
2413 : ! calculate overlaps between test functions and basis
2414 35 : DO l = 0, lmat
2415 335 : DO i = 1, 10
2416 300 : al = snl(l, i)
2417 300 : nl = nlmin(l)
2418 300 : pf = (2._dp*al)**nl*SQRT(2._dp*al/fac(2*nl))
2419 120300 : sfun(1:nr) = pf*basis%grid%rad(1:nr)**(nl - 1)*EXP(-al*basis%grid%rad(1:nr))
2420 1500 : DO j = 1, basis%nbas(l)
2421 481500 : sint(i, 1, j, l) = SUM(sfun(1:nr)*basis%bf(1:nr, j, l)*basis%grid%wr(1:nr))
2422 : END DO
2423 300 : nl = nlmax(l)
2424 300 : pf = (2._dp*al)**nl*SQRT(2._dp*al/fac(2*nl))
2425 120300 : sfun(1:nr) = pf*basis%grid%rad(1:nr)**(nl - 1)*EXP(-al*basis%grid%rad(1:nr))
2426 1530 : DO j = 1, basis%nbas(l)
2427 481500 : sint(i, 2, j, l) = SUM(sfun(1:nr)*basis%bf(1:nr, j, l)*basis%grid%wr(1:nr))
2428 : END DO
2429 : END DO
2430 : END DO
2431 :
2432 35 : DO l = 0, lmat
2433 30 : n = basis%nbas(l)
2434 30 : IF (n <= 0) CYCLE
2435 13 : m = basis%nprim(l)
2436 13 : SELECT CASE (basis%basis_type)
2437 : CASE DEFAULT
2438 0 : CPABORT("")
2439 : CASE (GTO_BASIS)
2440 10 : CALL sg_overlap(omat(1:n, 1:n, l), l, basis%am(1:n, l), basis%am(1:n, l))
2441 : CASE (CGTO_BASIS)
2442 12 : ALLOCATE (bmat(m, m))
2443 3 : CALL sg_overlap(bmat(1:m, 1:m), l, basis%am(1:m, l), basis%am(1:m, l))
2444 3 : CALL contract2(omat(1:n, 1:n, l), bmat(1:m, 1:m), basis%cm(1:m, 1:n, l))
2445 3 : DEALLOCATE (bmat)
2446 : CASE (STO_BASIS)
2447 : CALL sto_overlap(omat(1:n, 1:n, l), basis%ns(1:n, l), basis%as(1:n, l), &
2448 0 : basis%ns(1:n, l), basis%as(1:n, l))
2449 : CASE (NUM_BASIS)
2450 13 : CPABORT("")
2451 : END SELECT
2452 35 : CALL invmat_symm(omat(1:n, 1:n, l))
2453 : END DO
2454 :
2455 5 : WRITE (iw, '(/,A)') " Basis Set Completeness Estimates"
2456 35 : DO l = 0, lmat
2457 30 : n = basis%nbas(l)
2458 30 : IF (n <= 0) CYCLE
2459 13 : WRITE (iw, '(A,I3)') " L-quantum number: ", l
2460 13 : WRITE (iw, '(A,T31,A,I2,T61,A,I2)') " Slater Exponent", "Completeness n-qm=", nlmin(l), &
2461 26 : "Completeness n-qm=", nlmax(l)
2462 148 : DO i = 10, 1, -1
2463 21990 : c1 = DOT_PRODUCT(sint(i, 1, 1:n, l), MATMUL(omat(1:n, 1:n, l), sint(i, 1, 1:n, l)))
2464 21990 : c2 = DOT_PRODUCT(sint(i, 2, 1:n, l), MATMUL(omat(1:n, 1:n, l), sint(i, 2, 1:n, l)))
2465 160 : WRITE (iw, "(T6,F14.6,T41,F10.6,T71,F10.6)") snl(l, i), c1, c2
2466 : END DO
2467 : END DO
2468 :
2469 5 : DEALLOCATE (omat, sfun, sint)
2470 :
2471 5 : END SUBROUTINE atom_completeness
2472 :
2473 : ! **************************************************************************************************
2474 : !> \brief Calculate the condition number of the given atomic basis set.
2475 : !> \param basis atomic basis set
2476 : !> \param rad lattice constant (e.g. doubled atomic covalent radius)
2477 : !> \param cnum (output) condition number
2478 : ! **************************************************************************************************
2479 104 : SUBROUTINE atom_basis_condnum(basis, rad, cnum)
2480 : TYPE(atom_basis_type) :: basis
2481 : REAL(KIND=dp), INTENT(IN) :: rad
2482 : REAL(KIND=dp), INTENT(OUT) :: cnum
2483 :
2484 : INTEGER :: ia, ib, imax, info, ix, iy, iz, ja, jb, &
2485 : ka, kb, l, la, lb, lwork, na, nb, &
2486 : nbas, nna, nnb
2487 104 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: ibptr
2488 : REAL(KIND=dp) :: r1, r2, reps, rmax
2489 104 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: weig, work
2490 104 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: smat
2491 : REAL(KIND=dp), DIMENSION(2*lmat+1, 2*lmat+1) :: sab
2492 : REAL(KIND=dp), DIMENSION(3) :: rab
2493 104 : REAL(KIND=dp), DIMENSION(:), POINTER :: zeta, zetb
2494 :
2495 : ! Number of spherical Gaussian orbitals with angular momentum lmat: nso(lmat) = 2*lmat+1
2496 :
2497 : ! total number of basis functions
2498 104 : nbas = 0
2499 728 : DO l = 0, lmat
2500 728 : nbas = nbas + basis%nbas(l)*(2*l + 1)
2501 : END DO
2502 :
2503 624 : ALLOCATE (smat(nbas, nbas), ibptr(nbas, 0:lmat))
2504 175244 : smat = 0.0_dp
2505 20672 : ibptr = 0
2506 : na = 0
2507 728 : DO l = 0, lmat
2508 2276 : DO ia = 1, basis%nbas(l)
2509 1548 : ibptr(ia, l) = na + 1
2510 2172 : na = na + (2*l + 1)
2511 : END DO
2512 : END DO
2513 :
2514 104 : reps = 1.e-14_dp
2515 104 : IF (basis%basis_type == GTO_BASIS .OR. &
2516 : basis%basis_type == CGTO_BASIS) THEN
2517 728 : DO la = 0, lmat
2518 624 : na = basis%nprim(la)
2519 624 : nna = 2*la + 1
2520 624 : IF (na == 0) CYCLE
2521 238 : zeta => basis%am(:, la)
2522 1770 : DO lb = 0, lmat
2523 1428 : nb = basis%nprim(lb)
2524 1428 : nnb = 2*lb + 1
2525 1428 : IF (nb == 0) CYCLE
2526 566 : zetb => basis%am(:, lb)
2527 5768 : DO ia = 1, na
2528 56004 : DO ib = 1, nb
2529 49998 : IF (rad < 0.1_dp) THEN
2530 : imax = 0
2531 : ELSE
2532 46083 : r1 = exp_radius(la, zeta(ia), reps, 1.0_dp)
2533 46083 : r2 = exp_radius(lb, zetb(ib), reps, 1.0_dp)
2534 46083 : rmax = MAX(2._dp*r1, 2._dp*r2)
2535 46083 : imax = INT(rmax/rad) + 1
2536 : END IF
2537 50661 : IF (imax > 1) THEN
2538 36030 : CALL overlap_ab_sp(la, zeta(ia), lb, zetb(ib), rad, sab)
2539 36030 : IF (basis%basis_type == GTO_BASIS) THEN
2540 24238 : ja = ibptr(ia, la)
2541 24238 : jb = ibptr(ib, lb)
2542 194898 : smat(ja:ja + nna - 1, jb:jb + nnb - 1) = smat(ja:ja + nna - 1, jb:jb + nnb - 1) + sab(1:nna, 1:nnb)
2543 11792 : ELSEIF (basis%basis_type == CGTO_BASIS) THEN
2544 48574 : DO ka = 1, basis%nbas(la)
2545 181930 : DO kb = 1, basis%nbas(lb)
2546 133356 : ja = ibptr(ka, la)
2547 133356 : jb = ibptr(kb, lb)
2548 : smat(ja:ja + nna - 1, jb:jb + nnb - 1) = smat(ja:ja + nna - 1, jb:jb + nnb - 1) + &
2549 899042 : sab(1:nna, 1:nnb)*basis%cm(ia, ka, la)*basis%cm(ib, kb, lb)
2550 : END DO
2551 : END DO
2552 : END IF
2553 : ELSE
2554 48042 : DO ix = -imax, imax
2555 34074 : rab(1) = rad*ix
2556 142434 : DO iy = -imax, imax
2557 94392 : rab(2) = rad*iy
2558 403812 : DO iz = -imax, imax
2559 275346 : rab(3) = rad*iz
2560 275346 : CALL overlap_ab_s(la, zeta(ia), lb, zetb(ib), rab, sab)
2561 369738 : IF (basis%basis_type == GTO_BASIS) THEN
2562 269262 : ja = ibptr(ia, la)
2563 269262 : jb = ibptr(ib, lb)
2564 : smat(ja:ja + nna - 1, jb:jb + nnb - 1) = smat(ja:ja + nna - 1, jb:jb + nnb - 1) &
2565 1510034 : + sab(1:nna, 1:nnb)
2566 6084 : ELSEIF (basis%basis_type == CGTO_BASIS) THEN
2567 24030 : DO ka = 1, basis%nbas(la)
2568 79902 : DO kb = 1, basis%nbas(lb)
2569 55872 : ja = ibptr(ka, la)
2570 55872 : jb = ibptr(kb, lb)
2571 : smat(ja:ja + nna - 1, jb:jb + nnb - 1) = &
2572 : smat(ja:ja + nna - 1, jb:jb + nnb - 1) + &
2573 252778 : sab(1:nna, 1:nnb)*basis%cm(ia, ka, la)*basis%cm(ib, kb, lb)
2574 : END DO
2575 : END DO
2576 : END IF
2577 : END DO
2578 : END DO
2579 : END DO
2580 : END IF
2581 : END DO
2582 : END DO
2583 : END DO
2584 : END DO
2585 : ELSE
2586 0 : CPABORT("Condition number not available for this basis type")
2587 : END IF
2588 :
2589 104 : info = 0
2590 104 : lwork = MAX(nbas*nbas, nbas + 100)
2591 520 : ALLOCATE (weig(nbas), work(lwork))
2592 :
2593 104 : CALL lapack_ssyev("N", "U", nbas, smat, nbas, weig, work, lwork, info)
2594 104 : CPASSERT(info == 0)
2595 104 : IF (weig(1) < 0.0_dp) THEN
2596 14 : cnum = 100._dp
2597 : ELSE
2598 90 : cnum = ABS(weig(nbas)/weig(1))
2599 90 : cnum = LOG10(cnum)
2600 : END IF
2601 :
2602 104 : DEALLOCATE (smat, ibptr, weig, work)
2603 :
2604 104 : END SUBROUTINE atom_basis_condnum
2605 :
2606 : ! **************************************************************************************************
2607 : !> \brief Transform a matrix expressed in terms of a uncontracted basis set to a contracted one.
2608 : !> \param int (output) contracted matrix
2609 : !> \param omat uncontracted matrix
2610 : !> \param cm matrix of contraction coefficients
2611 : ! **************************************************************************************************
2612 64176 : SUBROUTINE contract2(int, omat, cm)
2613 : REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: int
2614 : REAL(dp), DIMENSION(:, :), INTENT(IN) :: omat, cm
2615 :
2616 : CHARACTER(len=*), PARAMETER :: routineN = 'contract2'
2617 :
2618 : INTEGER :: handle, m, n
2619 :
2620 64176 : CALL timeset(routineN, handle)
2621 :
2622 64176 : n = SIZE(int, 1)
2623 64176 : m = SIZE(omat, 1)
2624 :
2625 8662880 : INT(1:n, 1:n) = MATMUL(TRANSPOSE(cm(1:m, 1:n)), MATMUL(omat(1:m, 1:m), cm(1:m, 1:n)))
2626 :
2627 64176 : CALL timestop(handle)
2628 :
2629 64176 : END SUBROUTINE contract2
2630 :
2631 : ! **************************************************************************************************
2632 : !> \brief Same as contract2(), but add the new contracted matrix to the old one
2633 : !> instead of overwriting it.
2634 : !> \param int (input/output) contracted matrix to be updated
2635 : !> \param omat uncontracted matrix
2636 : !> \param cm matrix of contraction coefficients
2637 : ! **************************************************************************************************
2638 32140 : SUBROUTINE contract2add(int, omat, cm)
2639 : REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: int
2640 : REAL(dp), DIMENSION(:, :), INTENT(IN) :: omat, cm
2641 :
2642 : CHARACTER(len=*), PARAMETER :: routineN = 'contract2add'
2643 :
2644 : INTEGER :: handle, m, n
2645 :
2646 32140 : CALL timeset(routineN, handle)
2647 :
2648 32140 : n = SIZE(int, 1)
2649 32140 : m = SIZE(omat, 1)
2650 :
2651 3400530 : INT(1:n, 1:n) = INT(1:n, 1:n) + MATMUL(TRANSPOSE(cm(1:m, 1:n)), MATMUL(omat(1:m, 1:m), cm(1:m, 1:n)))
2652 :
2653 32140 : CALL timestop(handle)
2654 :
2655 32140 : END SUBROUTINE contract2add
2656 :
2657 : ! **************************************************************************************************
2658 : !> \brief Contract a matrix of Electron Repulsion Integrals (ERI-s).
2659 : !> \param eri (output) contracted ERI-s
2660 : !> \param omat uncontracted ERI-s
2661 : !> \param cm1 first matrix of contraction coefficients
2662 : !> \param cm2 second matrix of contraction coefficients
2663 : ! **************************************************************************************************
2664 1232 : SUBROUTINE contract4(eri, omat, cm1, cm2)
2665 : REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: eri
2666 : REAL(dp), DIMENSION(:, :), INTENT(IN) :: omat, cm1, cm2
2667 :
2668 : CHARACTER(len=*), PARAMETER :: routineN = 'contract4'
2669 :
2670 : INTEGER :: handle, i1, i2, m1, m2, mm1, mm2, n1, &
2671 : n2, nn1, nn2
2672 1232 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: amat, atran, bmat, btran, hint
2673 :
2674 1232 : CALL timeset(routineN, handle)
2675 :
2676 1232 : m1 = SIZE(cm1, 1)
2677 1232 : n1 = SIZE(cm1, 2)
2678 1232 : m2 = SIZE(cm2, 1)
2679 1232 : n2 = SIZE(cm2, 2)
2680 1232 : nn1 = SIZE(eri, 1)
2681 1232 : nn2 = SIZE(eri, 2)
2682 1232 : mm1 = SIZE(omat, 1)
2683 1232 : mm2 = SIZE(omat, 2)
2684 :
2685 7680 : ALLOCATE (amat(m1, m1), atran(n1, n1), bmat(m2, m2), btran(n2, n2))
2686 2844 : ALLOCATE (hint(mm1, nn2))
2687 :
2688 1894 : DO i1 = 1, mm1
2689 662 : CALL iunpack(bmat(1:m2, 1:m2), omat(i1, 1:mm2), m2)
2690 662 : CALL contract2(btran(1:n2, 1:n2), bmat(1:m2, 1:m2), cm2(1:m2, 1:n2))
2691 1894 : CALL ipack(btran(1:n2, 1:n2), hint(i1, 1:nn2), n2)
2692 : END DO
2693 :
2694 2330 : DO i2 = 1, nn2
2695 1098 : CALL iunpack(amat(1:m1, 1:m1), hint(1:mm1, i2), m1)
2696 1098 : CALL contract2(atran(1:n1, 1:n1), amat(1:m1, 1:m1), cm1(1:m1, 1:n1))
2697 2330 : CALL ipack(atran(1:n1, 1:n1), eri(1:nn1, i2), n1)
2698 : END DO
2699 :
2700 1232 : DEALLOCATE (amat, atran, bmat, btran)
2701 1232 : DEALLOCATE (hint)
2702 :
2703 1232 : CALL timestop(handle)
2704 :
2705 1232 : END SUBROUTINE contract4
2706 :
2707 : ! **************************************************************************************************
2708 : !> \brief Pack an n x n square real matrix into a 1-D vector.
2709 : !> \param mat matrix to pack
2710 : !> \param vec resulting vector
2711 : !> \param n size of the matrix
2712 : ! **************************************************************************************************
2713 1760 : PURE SUBROUTINE ipack(mat, vec, n)
2714 : REAL(dp), DIMENSION(:, :), INTENT(IN) :: mat
2715 : REAL(dp), DIMENSION(:), INTENT(INOUT) :: vec
2716 : INTEGER, INTENT(IN) :: n
2717 :
2718 : INTEGER :: i, ij, j
2719 :
2720 1760 : ij = 0
2721 4500 : DO i = 1, n
2722 10110 : DO j = i, n
2723 5610 : ij = ij + 1
2724 8350 : vec(ij) = mat(i, j)
2725 : END DO
2726 : END DO
2727 :
2728 1760 : END SUBROUTINE ipack
2729 :
2730 : ! **************************************************************************************************
2731 : !> \brief Unpack a 1-D real vector as a n x n square matrix.
2732 : !> \param mat resulting matrix
2733 : !> \param vec vector to unpack
2734 : !> \param n size of the matrix
2735 : ! **************************************************************************************************
2736 1760 : PURE SUBROUTINE iunpack(mat, vec, n)
2737 : REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: mat
2738 : REAL(dp), DIMENSION(:), INTENT(IN) :: vec
2739 : INTEGER, INTENT(IN) :: n
2740 :
2741 : INTEGER :: i, ij, j
2742 :
2743 1760 : ij = 0
2744 6541 : DO i = 1, n
2745 23252 : DO j = i, n
2746 16711 : ij = ij + 1
2747 16711 : mat(i, j) = vec(ij)
2748 21492 : mat(j, i) = vec(ij)
2749 : END DO
2750 : END DO
2751 :
2752 1760 : END SUBROUTINE iunpack
2753 :
2754 1021427 : END MODULE atom_utils
|
