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