Line data Source code
1 : !--------------------------------------------------------------------------------------------------!
2 : ! CP2K: A general program to perform molecular dynamics simulations !
3 : ! Copyright 2000-2026 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 480 : 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 480 : occupation = 0._dp
115 :
116 480 : CPASSERT(ASSOCIATED(ostring))
117 480 : CPASSERT(SIZE(ostring) > 0)
118 :
119 480 : no = SIZE(ostring)
120 :
121 480 : is = 1
122 : ! look for multiplicity
123 480 : mult = -1 !not specified
124 480 : IF (is <= no) THEN
125 480 : 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 480 : IF (is <= no) THEN
142 480 : IF (INDEX(ostring(is), "CORE") /= 0) is = is + 1 !Pseudopotential detected
143 : END IF
144 480 : IF (is <= no) THEN
145 480 : IF (INDEX(ostring(is), "none") /= 0) is = is + 1 !no electrons, used with CORE
146 : END IF
147 480 : IF (is <= no) THEN
148 454 : IF (INDEX(ostring(is), "[") /= 0) THEN
149 : ! core occupation from element [XX]
150 310 : i1 = INDEX(ostring(is), "[")
151 310 : i2 = INDEX(ostring(is), "]")
152 310 : CPASSERT((i2 - i1 - 1 > 0) .AND. (i2 - i1 - 1 < 3))
153 310 : elem = ostring(is) (i1 + 1:i2 - 1)
154 310 : ielem = 0
155 9644 : DO k = 1, nelem
156 9644 : IF (elem == ptable(k)%symbol) THEN
157 : ielem = k
158 : EXIT
159 : END IF
160 : END DO
161 310 : CPASSERT(ielem /= 0)
162 1550 : DO l = 0, MIN(lmat, UBOUND(ptable(ielem)%e_conv, 1))
163 1240 : el = 2._dp*(2._dp*REAL(l, dp) + 1._dp)
164 1240 : e0 = ptable(ielem)%e_conv(l)
165 2866 : DO k = 1, 10
166 2556 : occupation(l, k) = MIN(el, e0)
167 2556 : e0 = e0 - el
168 2556 : IF (e0 <= 0._dp) EXIT
169 : END DO
170 : END DO
171 310 : is = is + 1
172 : END IF
173 :
174 : END IF
175 :
176 1228 : DO i = is, no
177 748 : pstring = ostring(i)
178 748 : CALL uppercase(pstring)
179 748 : js = INDEX(pstring, "S")
180 748 : jp = INDEX(pstring, "P")
181 748 : jd = INDEX(pstring, "D")
182 748 : jf = INDEX(pstring, "F")
183 748 : CPASSERT(js + jp + jd + jf > 0)
184 748 : IF (js > 0) THEN
185 396 : CPASSERT(jp + jd + jf == 0)
186 396 : READ (pstring(1:js - 1), *) n
187 396 : READ (pstring(js + 1:), *) oo
188 396 : CPASSERT(n > 0)
189 396 : CPASSERT(oo >= 0._dp)
190 396 : CPASSERT(occupation(0, n) == 0)
191 396 : occupation(0, n) = oo
192 : END IF
193 748 : IF (jp > 0) THEN
194 136 : CPASSERT(js + jd + jf == 0)
195 136 : READ (pstring(1:jp - 1), *) n
196 136 : READ (pstring(jp + 1:), *) oo
197 136 : CPASSERT(n > 1)
198 136 : CPASSERT(oo >= 0._dp)
199 136 : CPASSERT(occupation(1, n - 1) == 0)
200 136 : occupation(1, n - 1) = oo
201 : END IF
202 748 : IF (jd > 0) THEN
203 126 : CPASSERT(js + jp + jf == 0)
204 126 : READ (pstring(1:jd - 1), *) n
205 126 : READ (pstring(jd + 1:), *) oo
206 126 : CPASSERT(n > 2)
207 126 : CPASSERT(oo >= 0._dp)
208 126 : CPASSERT(occupation(2, n - 2) == 0)
209 126 : occupation(2, n - 2) = oo
210 : END IF
211 1228 : IF (jf > 0) THEN
212 90 : CPASSERT(js + jp + jd == 0)
213 90 : READ (pstring(1:jf - 1), *) n
214 90 : READ (pstring(jf + 1:), *) oo
215 90 : CPASSERT(n > 3)
216 90 : CPASSERT(oo >= 0._dp)
217 90 : CPASSERT(occupation(3, n - 3) == 0)
218 90 : occupation(3, n - 3) = oo
219 : END IF
220 :
221 : END DO
222 :
223 480 : wfnocc = 0._dp
224 3360 : DO l = 0, lmat
225 : k = 0
226 32160 : DO i = 1, 10
227 31680 : IF (occupation(l, i) /= 0._dp) THEN
228 2774 : k = k + 1
229 2774 : wfnocc(l, k) = occupation(l, i)
230 : END IF
231 : END DO
232 : END DO
233 :
234 : !Check for consistency with multiplicity
235 480 : 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 : 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 480 : IF (PRESENT(multiplicity)) multiplicity = mult
270 :
271 480 : 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 15576 : 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 15576 : maxl = 0
287 109032 : DO l = 0, lmat
288 1043592 : IF (SUM(occupation(l, :)) /= 0._dp) maxl = l
289 : END DO
290 :
291 15576 : 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 25607 : 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 179249 : maxn = 0
307 179249 : DO l = 0, lmat
308 1715669 : DO k = 1, 10
309 1690062 : IF (occupation(l, k) /= 0._dp) maxn(l) = maxn(l) + 1
310 : END DO
311 : END DO
312 :
313 25607 : 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 84345 : 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 39159051 : pmat = 0._dp
337 84345 : n = SIZE(wfn, 2)
338 243818 : DO l = 0, maxl
339 410991 : DO i = 1, MIN(n, maxn(l))
340 1700482 : DO k = 1, nbas(l)
341 33517043 : DO j = 1, nbas(l)
342 33349870 : 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 84345 : 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 173859 : 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 173859 : my_typ = "RHO"
378 173859 : IF (PRESENT(typ)) my_typ = typ(1:3)
379 173859 : IF (my_typ == "KIN") THEN
380 112 : CPASSERT(PRESENT(rr))
381 : END IF
382 :
383 70726189 : density = 0._dp
384 536973 : DO l = 0, maxl
385 363114 : n = basis%nbas(l)
386 2754472 : DO i = 1, n
387 21372367 : DO j = i, n
388 18791754 : ff = pmat(i, j, l)
389 18791754 : IF (i /= j) ff = 2._dp*pmat(i, j, l)
390 21009253 : IF (my_typ == "RHO") THEN
391 5425705847 : density(:) = density(:) + ff*basis%bf(:, i, l)*basis%bf(:, j, l)
392 5268602 : ELSE IF (my_typ == "DER") THEN
393 : density(:) = density(:) + ff*basis%dbf(:, i, l)*basis%bf(:, j, l) &
394 2002103356 : + 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("Unknown matrix type specified. Check the code!")
406 : END IF
407 : END DO
408 : END DO
409 : END DO
410 : ! this factor from the product of two spherical harmonics
411 70726189 : density = density/fourpi
412 :
413 173859 : 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 /= 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) /= 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 /= 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("Unable to continue reading external density file")
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)) > 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("Unable to continue reading external density file")
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 /= 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("Unable to continue reading external v_xc file")
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)) > 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("Unable to continue reading external v_xc file")
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 1395 : 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 1395 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: den
662 :
663 1395 : charge = 0._dp
664 1395 : m = SIZE(basis%bf, 1)
665 4185 : ALLOCATE (den(m))
666 1395 : n = basis%nbas(l)
667 1395 : den = 0._dp
668 31555 : DO i = 1, n
669 918795 : DO j = 1, n
670 887240 : ff = wfn(i)*wfn(j)
671 334534000 : 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 528695 : DO i = 1, m
675 528695 : IF (basis%grid%rad(i) > rcov) den(i) = 0._dp
676 : END DO
677 528695 : charge = SUM(den(1:m)*basis%grid%wr(1:m))
678 1395 : DEALLOCATE (den)
679 :
680 1395 : 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 1745 : 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 1745 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: fe, rc, rhoc, rval
703 :
704 1745 : my_typ = "RHO"
705 1745 : IF (PRESENT(typ)) my_typ = typ(1:3)
706 :
707 3490 : SELECT CASE (potential%ppot_type)
708 : CASE (no_pseudo, ecp_pseudo)
709 : ! we do nothing
710 : CASE (gth_pseudo)
711 1745 : IF (potential%gth_pot%nlcc) THEN
712 1745 : m = SIZE(corden)
713 6980 : ALLOCATE (fe(m), rc(m))
714 1745 : n = potential%gth_pot%nexp_nlcc
715 3490 : DO i = 1, n
716 1745 : a = potential%gth_pot%alpha_nlcc(i)
717 1745 : a2 = a*a
718 : ! note that all terms are computed with rc, not rr
719 699745 : rc(:) = rr(:)/a
720 699745 : fe(:) = EXP(-0.5_dp*rc(:)*rc(:))
721 5271 : DO j = 1, potential%gth_pot%nct_nlcc(i)
722 1781 : cval = potential%gth_pot%cval_nlcc(j, i)
723 3526 : IF (my_typ == "RHO") THEN
724 389772 : corden(:) = corden(:) + fe(:)*rc**(2*j - 2)*cval
725 809 : ELSE IF (my_typ == "DER") THEN
726 324409 : corden(:) = corden(:) - fe(:)*rc**(2*j - 1)*cval/a
727 809 : 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 : CALL cp_abort(__LOCATION__, &
741 : "Only <RHO>, <DER>, <LAP> are supported as <my_typ> "// &
742 0 : "in atom_core_density, found <"//TRIM(my_typ)//">")
743 : END IF
744 : END DO
745 : END DO
746 1745 : DEALLOCATE (fe, rc)
747 : END IF
748 : CASE (upf_pseudo)
749 0 : IF (potential%upf_pot%core_correction) THEN
750 0 : m = SIZE(corden)
751 0 : n = potential%upf_pot%mesh_size
752 0 : reverse = .FALSE.
753 0 : IF (rr(1) > rr(m)) reverse = .TRUE.
754 0 : ALLOCATE (rhoc(m), rval(m))
755 0 : IF (reverse) THEN
756 0 : DO i = 1, m
757 0 : rval(i) = rr(m - i + 1)
758 : END DO
759 : ELSE
760 0 : rval(1:m) = rr(1:m)
761 : END IF
762 0 : IF (my_typ == "RHO") THEN
763 : CALL spline3ders(potential%upf_pot%r(1:n), potential%upf_pot%rho_nlcc(1:n), rval(1:m), &
764 0 : ynew=rhoc(1:m))
765 0 : ELSE IF (my_typ == "DER") THEN
766 : CALL spline3ders(potential%upf_pot%r(1:n), potential%upf_pot%rho_nlcc(1:n), rval(1:m), &
767 0 : dynew=rhoc(1:m))
768 0 : ELSE IF (my_typ == "LAP") THEN
769 : CALL spline3ders(potential%upf_pot%r(1:n), potential%upf_pot%rho_nlcc(1:n), rval(1:m), &
770 0 : d2ynew=rhoc(1:m))
771 : ELSE
772 : CALL cp_abort(__LOCATION__, &
773 : "Only <RHO>, <DER>, <LAP> are supported as <my_typ> "// &
774 0 : "in atom_core_density, found <"//TRIM(my_typ)//">")
775 : END IF
776 0 : IF (reverse) THEN
777 0 : DO i = 1, m
778 0 : rval(i) = rr(m - i + 1)
779 0 : corden(i) = corden(i) + rhoc(m - i + 1)
780 : END DO
781 : ELSE
782 0 : corden(1:m) = corden(1:m) + rhoc(1:m)
783 : END IF
784 0 : DEALLOCATE (rhoc, rval)
785 : END IF
786 : CASE (sgp_pseudo)
787 0 : IF (potential%sgp_pot%has_nlcc) THEN
788 0 : CPABORT("not implemented")
789 : END IF
790 : CASE DEFAULT
791 1745 : CPABORT("Unknown PP type")
792 : END SELECT
793 :
794 1745 : END SUBROUTINE atom_core_density
795 :
796 : ! **************************************************************************************************
797 : !> \brief ...
798 : !> \param locpot ...
799 : !> \param gthpot ...
800 : !> \param rr ...
801 : ! **************************************************************************************************
802 9 : PURE SUBROUTINE atom_local_potential(locpot, gthpot, rr)
803 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: locpot
804 : TYPE(atom_gthpot_type), INTENT(IN) :: gthpot
805 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rr
806 :
807 : INTEGER :: i, j, m, n
808 : REAL(KIND=dp) :: a
809 9 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: fe, rc
810 :
811 9 : m = SIZE(locpot)
812 36 : ALLOCATE (fe(m), rc(m))
813 6783 : rc(:) = rr(:)/gthpot%rc
814 6783 : DO i = 1, m
815 6783 : locpot(i) = -gthpot%zion*erf(rc(i)/SQRT(2._dp))/rr(i)
816 : END DO
817 9 : n = gthpot%ncl
818 6783 : fe(:) = EXP(-0.5_dp*rc(:)*rc(:))
819 25 : DO i = 1, n
820 11715 : locpot(:) = locpot(:) + fe(:)*rc**(2*i - 2)*gthpot%cl(i)
821 : END DO
822 9 : IF (gthpot%lpotextended) THEN
823 0 : DO j = 1, gthpot%nexp_lpot
824 0 : a = gthpot%alpha_lpot(j)
825 0 : rc(:) = rr(:)/a
826 0 : fe(:) = EXP(-0.5_dp*rc(:)*rc(:))
827 0 : n = gthpot%nct_lpot(j)
828 0 : DO i = 1, n
829 0 : locpot(:) = locpot(:) + fe(:)*rc**(2*i - 2)*gthpot%cval_lpot(i, j)
830 : END DO
831 : END DO
832 : END IF
833 9 : DEALLOCATE (fe, rc)
834 :
835 9 : END SUBROUTINE atom_local_potential
836 :
837 : ! **************************************************************************************************
838 : !> \brief ...
839 : !> \param rmax ...
840 : !> \param wfn ...
841 : !> \param rcov ...
842 : !> \param l ...
843 : !> \param basis ...
844 : ! **************************************************************************************************
845 84 : PURE SUBROUTINE atom_orbital_max(rmax, wfn, rcov, l, basis)
846 : REAL(KIND=dp), INTENT(OUT) :: rmax
847 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: wfn
848 : REAL(KIND=dp), INTENT(IN) :: rcov
849 : INTEGER, INTENT(IN) :: l
850 : TYPE(atom_basis_type), INTENT(IN) :: basis
851 :
852 : INTEGER :: i, m, n
853 : REAL(KIND=dp) :: ff
854 84 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: dorb
855 :
856 84 : m = SIZE(basis%bf, 1)
857 252 : ALLOCATE (dorb(m))
858 84 : n = basis%nbas(l)
859 84 : dorb = 0._dp
860 2878 : DO i = 1, n
861 2794 : ff = wfn(i)
862 664878 : dorb(1:m) = dorb(1:m) + ff*basis%dbf(1:m, i, l)
863 : END DO
864 84 : rmax = -1._dp
865 19700 : DO i = 1, m - 1
866 19700 : IF (basis%grid%rad(i) < 2*rcov) THEN
867 12798 : IF (dorb(i)*dorb(i + 1) < 0._dp) THEN
868 102 : rmax = MAX(rmax, basis%grid%rad(i))
869 : END IF
870 : END IF
871 : END DO
872 84 : DEALLOCATE (dorb)
873 :
874 84 : END SUBROUTINE atom_orbital_max
875 :
876 : ! **************************************************************************************************
877 : !> \brief ...
878 : !> \param node ...
879 : !> \param wfn ...
880 : !> \param rcov ...
881 : !> \param l ...
882 : !> \param basis ...
883 : ! **************************************************************************************************
884 593 : PURE SUBROUTINE atom_orbital_nodes(node, wfn, rcov, l, basis)
885 : INTEGER, INTENT(OUT) :: node
886 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: wfn
887 : REAL(KIND=dp), INTENT(IN) :: rcov
888 : INTEGER, INTENT(IN) :: l
889 : TYPE(atom_basis_type), INTENT(IN) :: basis
890 :
891 : INTEGER :: i, m, n
892 : REAL(KIND=dp) :: ff
893 593 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: orb
894 :
895 593 : node = 0
896 593 : m = SIZE(basis%bf, 1)
897 1779 : ALLOCATE (orb(m))
898 593 : n = basis%nbas(l)
899 593 : orb = 0._dp
900 9793 : DO i = 1, n
901 9200 : ff = wfn(i)
902 3583393 : orb(1:m) = orb(1:m) + ff*basis%bf(1:m, i, l)
903 : END DO
904 231600 : DO i = 1, m - 1
905 231600 : IF (basis%grid%rad(i) < rcov) THEN
906 153617 : IF (orb(i)*orb(i + 1) < 0._dp) node = node + 1
907 : END IF
908 : END DO
909 593 : DEALLOCATE (orb)
910 :
911 593 : END SUBROUTINE atom_orbital_nodes
912 :
913 : ! **************************************************************************************************
914 : !> \brief ...
915 : !> \param value ...
916 : !> \param wfn ...
917 : !> \param basis ...
918 : ! **************************************************************************************************
919 299 : PURE SUBROUTINE atom_wfnr0(value, wfn, basis)
920 : REAL(KIND=dp), INTENT(OUT) :: value
921 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: wfn
922 : TYPE(atom_basis_type), INTENT(IN) :: basis
923 :
924 : INTEGER :: i, m, n
925 :
926 299 : value = 0._dp
927 116798 : m = MAXVAL(MINLOC(basis%grid%rad))
928 299 : n = basis%nbas(0)
929 5908 : DO i = 1, n
930 5908 : value = value + wfn(i)*basis%bf(m, i, 0)
931 : END DO
932 299 : END SUBROUTINE atom_wfnr0
933 :
934 : ! **************************************************************************************************
935 : !> \brief Solve the generalised eigenproblem for every angular momentum.
936 : !> \param hmat Hamiltonian matrix
937 : !> \param umat transformation matrix which reduces the overlap matrix to its unitary form
938 : !> \param orb atomic wavefunctions
939 : !> \param ener atomic orbital energies
940 : !> \param nb number of contracted basis functions
941 : !> \param nv number of linear-independent contracted basis functions
942 : !> \param maxl maximum angular momentum to consider
943 : !> \par History
944 : !> * 08.2008 created [Juerg Hutter]
945 : ! **************************************************************************************************
946 83924 : SUBROUTINE atom_solve(hmat, umat, orb, ener, nb, nv, maxl)
947 : REAL(KIND=dp), DIMENSION(:, :, 0:), INTENT(IN) :: hmat, umat
948 : REAL(KIND=dp), DIMENSION(:, :, 0:), INTENT(INOUT) :: orb
949 : REAL(KIND=dp), DIMENSION(:, 0:), INTENT(INOUT) :: ener
950 : INTEGER, DIMENSION(0:), INTENT(IN) :: nb, nv
951 : INTEGER, INTENT(IN) :: maxl
952 :
953 : INTEGER :: info, l, lwork, m, n
954 83924 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: w, work
955 83924 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: a
956 :
957 587468 : CPASSERT(ALL(nb >= nv))
958 :
959 6447800 : orb = 0._dp
960 250518 : DO l = 0, maxl
961 166594 : n = nb(l)
962 166594 : m = nv(l)
963 250518 : IF (n > 0 .AND. m > 0) THEN
964 162640 : lwork = 10*m
965 1301120 : ALLOCATE (a(n, n), w(n), work(lwork))
966 475330458 : 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)))
967 10525496 : CALL dsyev("V", "U", m, a(1:m, 1:m), m, w(1:m), work, lwork, info)
968 179859303 : a(1:n, 1:m) = MATMUL(umat(1:n, 1:m, l), a(1:m, 1:m))
969 :
970 162640 : m = MIN(m, SIZE(orb, 2))
971 2796539 : orb(1:n, 1:m, l) = a(1:n, 1:m)
972 387558 : ener(1:m, l) = w(1:m)
973 :
974 162640 : DEALLOCATE (a, w, work)
975 : END IF
976 : END DO
977 :
978 83924 : END SUBROUTINE atom_solve
979 :
980 : ! **************************************************************************************************
981 : !> \brief THIS FUNCTION IS NO LONGER IN USE.
982 : !> \param fun ...
983 : !> \param deps ...
984 : !> \return ...
985 : ! **************************************************************************************************
986 0 : FUNCTION prune_grid(fun, deps) RESULT(nc)
987 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: fun
988 : REAL(KIND=dp), INTENT(IN), OPTIONAL :: deps
989 : INTEGER :: nc
990 :
991 : INTEGER :: i, nr
992 : REAL(KIND=dp) :: meps
993 :
994 0 : meps = 1.e-14_dp
995 0 : IF (PRESENT(deps)) meps = deps
996 :
997 0 : nr = SIZE(fun)
998 0 : nc = 0
999 0 : DO i = nr, 1, -1
1000 0 : IF (ABS(fun(i)) > meps) THEN
1001 : nc = i
1002 : EXIT
1003 : END IF
1004 : END DO
1005 :
1006 0 : END FUNCTION prune_grid
1007 :
1008 : ! **************************************************************************************************
1009 : !> \brief Integrate a function on an atomic radial grid.
1010 : !> \param fun function to integrate
1011 : !> \param grid atomic radial grid
1012 : !> \return value of the integral
1013 : !> \par History
1014 : !> * 08.2008 created [Juerg Hutter]
1015 : ! **************************************************************************************************
1016 164216 : PURE FUNCTION integrate_grid_function1(fun, grid) RESULT(integral)
1017 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: fun
1018 : TYPE(grid_atom_type), INTENT(IN) :: grid
1019 : REAL(KIND=dp) :: integral
1020 :
1021 : INTEGER :: nc
1022 :
1023 164216 : nc = SIZE(fun)
1024 66543226 : integral = SUM(fun(1:nc)*grid%wr(1:nc))
1025 :
1026 164216 : END FUNCTION integrate_grid_function1
1027 :
1028 : ! **************************************************************************************************
1029 : !> \brief Integrate the product of two functions on an atomic radial grid.
1030 : !> \param fun1 first function
1031 : !> \param fun2 second function
1032 : !> \param grid atomic radial grid
1033 : !> \return value of the integral
1034 : !> \par History
1035 : !> * 08.2008 created [Juerg Hutter]
1036 : ! **************************************************************************************************
1037 685958 : PURE FUNCTION integrate_grid_function2(fun1, fun2, grid) RESULT(integral)
1038 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: fun1, fun2
1039 : TYPE(grid_atom_type), INTENT(IN) :: grid
1040 : REAL(KIND=dp) :: integral
1041 :
1042 : INTEGER :: nc
1043 :
1044 685958 : nc = MIN(SIZE(fun1), SIZE(fun2))
1045 275312382 : integral = SUM(fun1(1:nc)*fun2(1:nc)*grid%wr(1:nc))
1046 :
1047 685958 : END FUNCTION integrate_grid_function2
1048 :
1049 : ! **************************************************************************************************
1050 : !> \brief Integrate the product of three functions on an atomic radial grid.
1051 : !> \param fun1 first function
1052 : !> \param fun2 second function
1053 : !> \param fun3 third function
1054 : !> \param grid atomic radial grid
1055 : !> \return value of the integral
1056 : !> \par History
1057 : !> * 08.2008 created [Juerg Hutter]
1058 : ! **************************************************************************************************
1059 68806843 : PURE FUNCTION integrate_grid_function3(fun1, fun2, fun3, grid) RESULT(integral)
1060 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: fun1, fun2, fun3
1061 : TYPE(grid_atom_type), INTENT(IN) :: grid
1062 : REAL(KIND=dp) :: integral
1063 :
1064 : INTEGER :: nc
1065 :
1066 68806843 : nc = MIN(SIZE(fun1), SIZE(fun2), SIZE(fun3))
1067 27192720293 : integral = SUM(fun1(1:nc)*fun2(1:nc)*fun3(1:nc)*grid%wr(1:nc))
1068 :
1069 68806843 : END FUNCTION integrate_grid_function3
1070 :
1071 : ! **************************************************************************************************
1072 : !> \brief Numerically compute the Coulomb potential on an atomic radial grid.
1073 : !> \param cpot Coulomb potential on the radial grid
1074 : !> \param density electron density on the radial grid
1075 : !> \param grid atomic radial grid
1076 : !> \par History
1077 : !> * 08.2008 created [Juerg Hutter]
1078 : ! **************************************************************************************************
1079 93778 : SUBROUTINE coulomb_potential_numeric(cpot, density, grid)
1080 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: cpot
1081 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: density
1082 : TYPE(grid_atom_type), INTENT(IN) :: grid
1083 :
1084 : INTEGER :: i, nc
1085 : REAL(dp) :: int1, int2
1086 93778 : REAL(dp), DIMENSION(:), POINTER :: r, wr
1087 :
1088 93778 : nc = MIN(SIZE(cpot), SIZE(density))
1089 93778 : r => grid%rad
1090 93778 : wr => grid%wr
1091 :
1092 93778 : int1 = fourpi*integrate_grid(density, grid)
1093 93778 : int2 = 0._dp
1094 187556 : cpot(nc:) = int1/r(nc:)
1095 : ! IF (log_unit>0) WRITE(log_unit,"(A,2F10.8)") "r", int1, cpot(nc:)
1096 :
1097 : ! test that grid is decreasing
1098 93778 : CPASSERT(r(1) > r(nc))
1099 37966538 : DO i = 1, nc
1100 37872760 : cpot(i) = int1/r(i) + int2
1101 37872760 : int1 = int1 - fourpi*density(i)*wr(i)
1102 37966538 : int2 = int2 + fourpi*density(i)*wr(i)/r(i)
1103 : END DO
1104 :
1105 93778 : END SUBROUTINE coulomb_potential_numeric
1106 :
1107 : ! **************************************************************************************************
1108 : !> \brief Analytically compute the Coulomb potential on an atomic radial grid.
1109 : !> \param cpot Coulomb potential on the radial grid
1110 : !> \param pmat density matrix
1111 : !> \param basis atomic basis set
1112 : !> \param grid atomic radial grid
1113 : !> \param maxl maximum angular momentum to consider
1114 : !> \par History
1115 : !> * 08.2008 created [Juerg Hutter]
1116 : ! **************************************************************************************************
1117 38 : SUBROUTINE coulomb_potential_analytic(cpot, pmat, basis, grid, maxl)
1118 : REAL(KIND=dp), DIMENSION(:), INTENT(INOUT) :: cpot
1119 : REAL(KIND=dp), DIMENSION(:, :, 0:), INTENT(IN) :: pmat
1120 : TYPE(atom_basis_type), INTENT(IN) :: basis
1121 : TYPE(grid_atom_type) :: grid
1122 : INTEGER, INTENT(IN) :: maxl
1123 :
1124 : INTEGER :: i, j, k, l, m, n
1125 : REAL(KIND=dp) :: a, b, ff, oab, sab
1126 38 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: erfa, expa, z
1127 38 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: unp
1128 :
1129 38 : m = SIZE(cpot)
1130 190 : ALLOCATE (erfa(1:m), expa(1:m), z(1:m))
1131 :
1132 18838 : cpot = 0._dp
1133 :
1134 172 : DO l = 0, maxl
1135 65958 : IF (MAXVAL(ABS(pmat(:, :, l))) < 1.e-14_dp) CYCLE
1136 38 : SELECT CASE (basis%basis_type)
1137 : CASE DEFAULT
1138 0 : CPABORT("Unknown basis type for coulomb_potential_analytic")
1139 : CASE (GTO_BASIS)
1140 2266 : DO i = 1, basis%nbas(l)
1141 24774 : DO j = i, basis%nbas(l)
1142 22508 : IF (ABS(pmat(i, j, l)) < 1.e-14_dp) CYCLE
1143 22490 : ff = pmat(i, j, l)
1144 22490 : IF (i /= j) ff = 2._dp*ff
1145 22490 : a = basis%am(i, l)
1146 22490 : b = basis%am(j, l)
1147 22490 : sab = SQRT(a + b)
1148 22490 : oab = rootpi/(a + b)**(l + 1.5_dp)*ff
1149 12640090 : z(:) = sab*grid%rad(:)
1150 12640090 : DO k = 1, SIZE(erfa)
1151 12640090 : erfa(k) = oab*erf(z(k))/grid%rad(k)
1152 : END DO
1153 12640090 : expa(:) = EXP(-z(:)**2)*ff/(a + b)**(l + 1)
1154 2132 : SELECT CASE (l)
1155 : CASE (0)
1156 6153004 : cpot(:) = cpot(:) + 0.25_dp*erfa(:)
1157 : CASE (1)
1158 3956950 : cpot(:) = cpot(:) + 0.375_dp*erfa(:) - 0.25_dp*expa(:)
1159 : CASE (2)
1160 2096254 : cpot(:) = cpot(:) + 0.9375_dp*erfa(:) - expa(:)*(0.875_dp + 0.25_dp*z(:)**2)
1161 : CASE (3)
1162 433882 : cpot(:) = cpot(:) + 3.28125_dp*erfa(:) - expa(:)*(3.5625_dp + 1.375_dp*z(:)**2 + 0.25*z(:)**4)
1163 : CASE DEFAULT
1164 22508 : CPABORT("Invalid l number for GTO specified. Check the code!")
1165 : END SELECT
1166 : END DO
1167 : END DO
1168 : CASE (CGTO_BASIS)
1169 0 : n = basis%nprim(l)
1170 0 : m = basis%nbas(l)
1171 0 : ALLOCATE (unp(n, n))
1172 :
1173 0 : unp(1:n, 1:n) = MATMUL(MATMUL(basis%cm(1:n, 1:m, l), pmat(1:m, 1:m, l)), &
1174 0 : TRANSPOSE(basis%cm(1:n, 1:m, l)))
1175 0 : DO i = 1, basis%nprim(l)
1176 0 : DO j = i, basis%nprim(l)
1177 0 : IF (ABS(unp(i, j)) < 1.e-14_dp) CYCLE
1178 0 : ff = unp(i, j)
1179 0 : IF (i /= j) ff = 2._dp*ff
1180 0 : a = basis%am(i, l)
1181 0 : b = basis%am(j, l)
1182 0 : sab = SQRT(a + b)
1183 0 : oab = rootpi/(a + b)**(l + 1.5_dp)*ff
1184 0 : z(:) = sab*grid%rad(:)
1185 0 : DO k = 1, SIZE(erfa)
1186 0 : erfa(k) = oab*erf(z(k))/grid%rad(k)
1187 : END DO
1188 0 : expa(:) = EXP(-z(:)**2)*ff/(a + b)**(l + 1)
1189 0 : SELECT CASE (l)
1190 : CASE (0)
1191 0 : cpot(:) = cpot(:) + 0.25_dp*erfa(:)
1192 : CASE (1)
1193 0 : cpot(:) = cpot(:) + 0.375_dp*erfa(:) - 0.25_dp*expa(:)
1194 : CASE (2)
1195 0 : cpot(:) = cpot(:) + 0.9375_dp*erfa(:) - expa(:)*(0.875_dp + 0.25_dp*z(:)**2)
1196 : CASE (3)
1197 0 : cpot(:) = cpot(:) + 3.28125_dp*erfa(:) - expa(:)*(3.5625_dp + 1.375_dp*z(:)**2 + 0.25*z(:)**4)
1198 : CASE DEFAULT
1199 0 : CPABORT("Invalid l number for CGTO specified. Check the code!")
1200 : END SELECT
1201 : END DO
1202 : END DO
1203 :
1204 134 : DEALLOCATE (unp)
1205 : END SELECT
1206 : END DO
1207 38 : DEALLOCATE (erfa, expa, z)
1208 :
1209 38 : END SUBROUTINE coulomb_potential_analytic
1210 :
1211 : ! **************************************************************************************************
1212 : !> \brief Calculate the exchange potential numerically.
1213 : !> \param kmat Kohn-Sham matrix
1214 : !> \param state electronic state
1215 : !> \param occ occupation numbers
1216 : !> \param wfn wavefunctions
1217 : !> \param basis atomic basis set
1218 : !> \param hfx_pot potential parameters from Hartree-Fock section
1219 : !> \par History
1220 : !> * 01.2009 created [Juerg Hutter]
1221 : ! **************************************************************************************************
1222 18038 : SUBROUTINE exchange_numeric(kmat, state, occ, wfn, basis, hfx_pot)
1223 : REAL(KIND=dp), DIMENSION(:, :, 0:), INTENT(INOUT) :: kmat
1224 : TYPE(atom_state), INTENT(IN) :: state
1225 : REAL(KIND=dp), DIMENSION(0:, :), INTENT(IN) :: occ
1226 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: wfn
1227 : TYPE(atom_basis_type), INTENT(IN) :: basis
1228 : TYPE(atom_hfx_type), INTENT(IN) :: hfx_pot
1229 :
1230 : INTEGER :: i, ia, ib, k, lad, lbc, lh, ll, nbas, &
1231 : norb, nr, nu
1232 : REAL(KIND=dp) :: almn
1233 18038 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cpot, nai, nbi, pot
1234 18038 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: orb
1235 : REAL(KIND=dp), DIMENSION(0:maxfac) :: arho
1236 :
1237 18038 : arho = 0._dp
1238 18038 : DO ll = 0, maxfac, 2
1239 288608 : lh = ll/2
1240 288608 : arho(ll) = fac(ll)/fac(lh)**2
1241 : END DO
1242 :
1243 4004426 : kmat = 0._dp
1244 :
1245 18038 : nr = basis%grid%nr
1246 108228 : ALLOCATE (nai(nr), nbi(nr), cpot(nr), pot(nr))
1247 :
1248 38714 : DO lad = 0, state%maxl_calc
1249 65994 : DO lbc = 0, state%maxl_occ
1250 27280 : norb = state%maxn_occ(lbc)
1251 27280 : nbas = basis%nbas(lbc)
1252 : ! calculate orbitals for angmom lbc
1253 109084 : ALLOCATE (orb(nr, norb))
1254 27280 : orb = 0._dp
1255 75576 : DO i = 1, norb
1256 585504 : DO k = 1, nbas
1257 202832624 : orb(:, i) = orb(:, i) + wfn(k, i, lbc)*basis%bf(:, k, lbc)
1258 : END DO
1259 : END DO
1260 61024 : DO nu = ABS(lad - lbc), lad + lbc, 2
1261 : almn = arho(-lad + lbc + nu)*arho(lad - lbc + nu)*arho(lad + lbc - nu)/(REAL(lad + lbc + nu + 1, dp) &
1262 33744 : *arho(lad + lbc + nu))
1263 33744 : almn = -0.5_dp*almn
1264 :
1265 340554 : DO ia = 1, basis%nbas(lad)
1266 1005554 : DO i = 1, norb
1267 275717080 : nai(:) = orb(:, i)*basis%bf(:, ia, lad)
1268 692280 : cpot = 0.0_dp
1269 692280 : IF (hfx_pot%scale_coulomb /= 0.0_dp) THEN
1270 657308 : CALL potential_coulomb_numeric(pot, nai, nu, basis%grid)
1271 263580508 : cpot(:) = cpot(:) + pot(:)*hfx_pot%scale_coulomb
1272 : END IF
1273 692280 : IF (hfx_pot%scale_longrange /= 0.0_dp) THEN
1274 34972 : IF (hfx_pot%do_gh) THEN
1275 : CALL potential_longrange_numeric_gh(pot, nai, nu, basis%grid, hfx_pot%omega, &
1276 17486 : hfx_pot%kernel(:, :, nu))
1277 : ELSE
1278 : CALL potential_longrange_numeric(pot, nai, nu, basis%grid, hfx_pot%omega, &
1279 17486 : hfx_pot%kernel(:, :, nu))
1280 : END IF
1281 12136572 : cpot(:) = cpot(:) + pot(:)*hfx_pot%scale_longrange
1282 : END IF
1283 13563642 : DO ib = 1, basis%nbas(lad)
1284 : kmat(ia, ib, lad) = kmat(ia, ib, lad) + almn*occ(lbc, i)* &
1285 13284112 : integrate_grid(cpot, orb(:, i), basis%bf(:, ib, lad), basis%grid)
1286 : END DO
1287 : END DO
1288 : END DO
1289 :
1290 : END DO
1291 47956 : DEALLOCATE (orb)
1292 : END DO
1293 : END DO
1294 :
1295 18038 : DEALLOCATE (nai, nbi, cpot)
1296 :
1297 18038 : END SUBROUTINE exchange_numeric
1298 :
1299 : ! **************************************************************************************************
1300 : !> \brief Calculate the exchange potential semi-analytically.
1301 : !> \param kmat Kohn-Sham matrix
1302 : !> \param state electronic state
1303 : !> \param occ occupation numbers
1304 : !> \param wfn wavefunctions
1305 : !> \param basis atomic basis set
1306 : !> \param hfx_pot properties of the Hartree-Fock potential
1307 : !> \par History
1308 : !> * 01.2009 created [Juerg Hutter]
1309 : ! **************************************************************************************************
1310 191 : SUBROUTINE exchange_semi_analytic(kmat, state, occ, wfn, basis, hfx_pot)
1311 : REAL(KIND=dp), DIMENSION(:, :, 0:), INTENT(INOUT) :: kmat
1312 : TYPE(atom_state), INTENT(IN) :: state
1313 : REAL(KIND=dp), DIMENSION(0:, :), INTENT(IN) :: occ
1314 : REAL(KIND=dp), DIMENSION(:, :, :), POINTER :: wfn
1315 : TYPE(atom_basis_type), INTENT(IN) :: basis
1316 : TYPE(atom_hfx_type), INTENT(IN) :: hfx_pot
1317 :
1318 : INTEGER :: i, ia, ib, k, lad, lbc, lh, ll, nbas, &
1319 : norb, nr, nu
1320 : REAL(KIND=dp) :: almn
1321 191 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: cpot, nai, nbi
1322 191 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: orb
1323 191 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: pot
1324 : REAL(KIND=dp), DIMENSION(0:maxfac) :: arho
1325 :
1326 191 : arho = 0._dp
1327 191 : DO ll = 0, maxfac, 2
1328 3056 : lh = ll/2
1329 3056 : arho(ll) = fac(ll)/fac(lh)**2
1330 : END DO
1331 :
1332 574097 : kmat = 0._dp
1333 :
1334 191 : nr = basis%grid%nr
1335 1337 : nbas = MAXVAL(basis%nbas)
1336 955 : ALLOCATE (pot(nr, nbas, nbas))
1337 955 : ALLOCATE (nai(nr), nbi(nr), cpot(nr))
1338 :
1339 764 : DO lad = 0, state%maxl_calc
1340 1910 : DO lbc = 0, state%maxl_occ
1341 1146 : norb = state%maxn_occ(lbc)
1342 1146 : nbas = basis%nbas(lbc)
1343 : ! calculate orbitals for angmom lbc
1344 4584 : ALLOCATE (orb(nr, norb))
1345 1146 : orb = 0._dp
1346 2370 : DO i = 1, norb
1347 29058 : DO k = 1, nbas
1348 9127512 : orb(:, i) = orb(:, i) + wfn(k, i, lbc)*basis%bf(:, k, lbc)
1349 : END DO
1350 : END DO
1351 2674 : DO nu = ABS(lad - lbc), lad + lbc, 2
1352 : almn = arho(-lad + lbc + nu)*arho(lad - lbc + nu) &
1353 1528 : *arho(lad + lbc - nu)/(REAL(lad + lbc + nu + 1, dp)*arho(lad + lbc + nu))
1354 1528 : almn = -0.5_dp*almn
1355 : ! calculate potential for basis function pair (lad,lbc)
1356 1528 : pot = 0._dp
1357 1528 : CALL potential_analytic(pot, lad, lbc, nu, basis, hfx_pot)
1358 34098 : DO ia = 1, basis%nbas(lad)
1359 66794 : DO i = 1, norb
1360 33842 : cpot = 0._dp
1361 801328 : DO k = 1, nbas
1362 249602528 : cpot(:) = cpot(:) + pot(:, ia, k)*wfn(k, i, lbc)
1363 : END DO
1364 813096 : DO ib = 1, basis%nbas(lad)
1365 : kmat(ia, ib, lad) = kmat(ia, ib, lad) + almn*occ(lbc, i)* &
1366 781672 : integrate_grid(cpot, orb(:, i), basis%bf(:, ib, lad), basis%grid)
1367 : END DO
1368 : END DO
1369 : END DO
1370 : END DO
1371 1719 : DEALLOCATE (orb)
1372 : END DO
1373 : END DO
1374 :
1375 191 : DEALLOCATE (pot)
1376 191 : DEALLOCATE (nai, nbi, cpot)
1377 :
1378 191 : END SUBROUTINE exchange_semi_analytic
1379 :
1380 : ! **************************************************************************************************
1381 : !> \brief Calculate the electrostatic potential using numerical differentiation.
1382 : !> \param cpot computed potential
1383 : !> \param density electron density on the atomic radial grid
1384 : !> \param nu integer parameter [ABS(la-lb) .. la+lb];
1385 : !> see potential_analytic() and exchange_numeric()
1386 : !> \param grid atomic radial grid
1387 : !> \par History
1388 : !> * 01.2009 created [Juerg Hutter]
1389 : ! **************************************************************************************************
1390 657308 : SUBROUTINE potential_coulomb_numeric(cpot, density, nu, grid)
1391 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: cpot
1392 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: density
1393 : INTEGER, INTENT(IN) :: nu
1394 : TYPE(grid_atom_type), INTENT(IN) :: grid
1395 :
1396 : INTEGER :: i, nc
1397 : REAL(dp) :: int1, int2
1398 657308 : REAL(dp), DIMENSION(:), POINTER :: r, wr
1399 :
1400 657308 : nc = MIN(SIZE(cpot), SIZE(density))
1401 657308 : r => grid%rad
1402 657308 : wr => grid%wr
1403 :
1404 263580508 : int1 = integrate_grid(density, r**nu, grid)
1405 657308 : int2 = 0._dp
1406 1314616 : cpot(nc:) = int1/r(nc:)**(nu + 1)
1407 :
1408 : ! test that grid is decreasing
1409 657308 : CPASSERT(r(1) > r(nc))
1410 263580508 : DO i = 1, nc
1411 262923200 : cpot(i) = int1/r(i)**(nu + 1) + int2*r(i)**nu
1412 262923200 : int1 = int1 - r(i)**(nu)*density(i)*wr(i)
1413 263580508 : int2 = int2 + density(i)*wr(i)/r(i)**(nu + 1)
1414 : END DO
1415 :
1416 657308 : END SUBROUTINE potential_coulomb_numeric
1417 :
1418 : ! **************************************************************************************************
1419 : !> \brief ...
1420 : !> \param cpot ...
1421 : !> \param density ...
1422 : !> \param nu ...
1423 : !> \param grid ...
1424 : !> \param omega ...
1425 : !> \param kernel ...
1426 : ! **************************************************************************************************
1427 17486 : SUBROUTINE potential_longrange_numeric(cpot, density, nu, grid, omega, kernel)
1428 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: cpot
1429 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: density
1430 : INTEGER, INTENT(IN) :: nu
1431 : TYPE(grid_atom_type), INTENT(IN) :: grid
1432 : REAL(KIND=dp), INTENT(IN) :: omega
1433 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: kernel
1434 :
1435 : INTEGER :: nc
1436 17486 : REAL(dp), DIMENSION(:), POINTER :: r, wr
1437 17486 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: work_inp, work_out
1438 :
1439 17486 : nc = MIN(SIZE(cpot), SIZE(density))
1440 17486 : r => grid%rad
1441 17486 : wr => grid%wr
1442 :
1443 69944 : ALLOCATE (work_inp(nc), work_out(nc))
1444 :
1445 6068286 : cpot = 0.0_dp
1446 :
1447 : ! First Bessel transform
1448 6068286 : work_inp(:nc) = density(:nc)*wr(:nc)
1449 17486 : CALL dsymv('U', nc, 1.0_dp, kernel, nc, work_inp, 1, 0.0_dp, work_out, 1)
1450 :
1451 : ! Second Bessel transform
1452 6068286 : work_inp(:nc) = work_out(:nc)*EXP(-r(:nc)**2)/r(:nc)**2*wr(:nc)
1453 17486 : CALL dsymv('U', nc, 1.0_dp, kernel, nc, work_inp, 1, 0.0_dp, work_out, 1)
1454 :
1455 6068286 : cpot(:nc) = work_out(:nc)*(2.0_dp*REAL(nu, dp) + 1.0_dp)*4.0_dp/pi*omega
1456 :
1457 34972 : END SUBROUTINE potential_longrange_numeric
1458 :
1459 : ! **************************************************************************************************
1460 : !> \brief ...
1461 : !> \param cpot ...
1462 : !> \param density ...
1463 : !> \param nu ...
1464 : !> \param grid ...
1465 : !> \param omega ...
1466 : !> \param kernel ...
1467 : ! **************************************************************************************************
1468 17486 : SUBROUTINE potential_longrange_numeric_gh(cpot, density, nu, grid, omega, kernel)
1469 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: cpot
1470 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: density
1471 : INTEGER, INTENT(IN) :: nu
1472 : TYPE(grid_atom_type), INTENT(IN) :: grid
1473 : REAL(KIND=dp), INTENT(IN) :: omega
1474 : REAL(KIND=dp), DIMENSION(:, :), INTENT(IN) :: kernel
1475 :
1476 : INTEGER :: n_max, nc, nc_kernel, nr_kernel
1477 17486 : REAL(dp), DIMENSION(:), POINTER :: wr
1478 17486 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: work_inp, work_out
1479 :
1480 17486 : nc = MIN(SIZE(cpot), SIZE(density))
1481 17486 : wr => grid%wr
1482 :
1483 17486 : nc_kernel = SIZE(kernel, 1)
1484 17486 : nr_kernel = SIZE(kernel, 2)
1485 17486 : n_max = MAX(nc, nc_kernel, nr_kernel)
1486 :
1487 69944 : ALLOCATE (work_inp(n_max), work_out(n_max))
1488 :
1489 6068286 : cpot = 0.0_dp
1490 :
1491 : ! First Bessel transform
1492 6068286 : work_inp(:nc) = density(:nc)*wr(:nc)
1493 17486 : CALL DGEMV('T', nc_kernel, nr_kernel, 1.0_dp, kernel, nc_kernel, work_inp, 1, 0.0_dp, work_out, 1)
1494 :
1495 : ! Second Bessel transform
1496 1766086 : work_inp(:nr_kernel) = work_out(:nr_kernel)
1497 17486 : CALL DGEMV('N', nc_kernel, nr_kernel, 1.0_dp, kernel, nc_kernel, work_inp, 1, 0.0_dp, work_out, 1)
1498 :
1499 6068286 : cpot(:nc) = work_out(:nc)*(2.0_dp*REAL(nu, dp) + 1.0_dp)*4.0_dp/pi*omega
1500 :
1501 34972 : END SUBROUTINE potential_longrange_numeric_gh
1502 :
1503 : ! **************************************************************************************************
1504 : !> \brief Calculate the electrostatic potential using analytical expressions.
1505 : !> The interaction potential has the form
1506 : !> V(r)=scale_coulomb*1/r+scale_lr*erf(omega*r)/r
1507 : !> \param cpot computed potential
1508 : !> \param la angular momentum of the calculated state
1509 : !> \param lb angular momentum of the occupied state
1510 : !> \param nu integer parameter [ABS(la-lb) .. la+lb] with the parity of 'la+lb'
1511 : !> \param basis atomic basis set
1512 : !> \param hfx_pot properties of the Hartree-Fock potential
1513 : !> \par History
1514 : !> * 01.2009 created [Juerg Hutter]
1515 : ! **************************************************************************************************
1516 1528 : SUBROUTINE potential_analytic(cpot, la, lb, nu, basis, hfx_pot)
1517 : REAL(KIND=dp), DIMENSION(:, :, :), INTENT(OUT) :: cpot
1518 : INTEGER, INTENT(IN) :: la, lb, nu
1519 : TYPE(atom_basis_type), INTENT(IN) :: basis
1520 : TYPE(atom_hfx_type), INTENT(IN) :: hfx_pot
1521 :
1522 : INTEGER :: i, j, k, l, ll, m
1523 : REAL(KIND=dp) :: a, b
1524 1528 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: erfa, pot
1525 :
1526 1528 : m = SIZE(cpot, 1)
1527 4584 : ALLOCATE (erfa(1:m))
1528 :
1529 1528 : ll = la + lb
1530 :
1531 242333208 : cpot = 0._dp
1532 :
1533 1528 : SELECT CASE (basis%basis_type)
1534 : CASE DEFAULT
1535 0 : CPABORT("Unknown basis type for potential_analytic")
1536 : CASE (GTO_BASIS)
1537 32952 : DO i = 1, basis%nbas(la)
1538 715808 : DO j = 1, basis%nbas(lb)
1539 682856 : a = basis%am(i, la)
1540 682856 : b = basis%am(j, lb)
1541 :
1542 682856 : IF (hfx_pot%scale_coulomb /= 0.0_dp) THEN
1543 329160 : CALL potential_coulomb_analytic(erfa, a, b, ll, nu, basis%grid%rad)
1544 :
1545 112946760 : cpot(:, i, j) = cpot(:, i, j) + erfa(:)*hfx_pot%scale_coulomb
1546 : END IF
1547 :
1548 714280 : IF (hfx_pot%scale_longrange /= 0.0_dp) THEN
1549 353696 : CALL potential_longrange_analytic(erfa, a, b, ll, nu, basis%grid%rad, hfx_pot%omega)
1550 :
1551 119611296 : cpot(:, i, j) = cpot(:, i, j) + erfa(:)*hfx_pot%scale_longrange
1552 : END IF
1553 : END DO
1554 : END DO
1555 : CASE (CGTO_BASIS)
1556 0 : ALLOCATE (pot(1:m))
1557 :
1558 1528 : DO i = 1, basis%nprim(la)
1559 0 : DO j = 1, basis%nprim(lb)
1560 0 : a = basis%am(i, la)
1561 0 : b = basis%am(j, lb)
1562 :
1563 0 : pot = 0.0_dp
1564 :
1565 0 : IF (hfx_pot%scale_coulomb /= 0.0_dp) THEN
1566 0 : CALL potential_coulomb_analytic(erfa, a, b, ll, nu, basis%grid%rad)
1567 :
1568 0 : pot(:) = pot(:) + erfa(:)*hfx_pot%scale_coulomb
1569 : END IF
1570 :
1571 0 : IF (hfx_pot%scale_longrange /= 0.0_dp) THEN
1572 0 : CALL potential_longrange_analytic(erfa, a, b, ll, nu, basis%grid%rad, hfx_pot%omega)
1573 :
1574 0 : pot(:) = pot(:) + erfa(:)*hfx_pot%scale_longrange
1575 : END IF
1576 :
1577 0 : DO k = 1, basis%nbas(la)
1578 0 : DO l = 1, basis%nbas(lb)
1579 0 : cpot(:, k, l) = cpot(:, k, l) + pot(:)*basis%cm(i, k, la)*basis%cm(j, l, lb)
1580 : END DO
1581 : END DO
1582 : END DO
1583 : END DO
1584 :
1585 : END SELECT
1586 :
1587 1528 : DEALLOCATE (erfa)
1588 :
1589 1528 : END SUBROUTINE potential_analytic
1590 :
1591 : ! **************************************************************************************************
1592 : !> \brief ...
1593 : !> \param erfa Array will contain the potential
1594 : !> \param a Exponent of first Gaussian charge distribution
1595 : !> \param b Exponent of second Gaussian charge distribution
1596 : !> \param ll Sum of angular momentum quantum numbers building one charge distribution
1597 : !> \param nu Angular momentum of interaction, (ll-nu) should be even.
1598 : !> \param rad Radial grid
1599 : ! **************************************************************************************************
1600 329160 : SUBROUTINE potential_coulomb_analytic(erfa, a, b, ll, nu, rad)
1601 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: erfa
1602 : REAL(KIND=dp), INTENT(IN) :: a, b
1603 : INTEGER, INTENT(IN) :: ll, nu
1604 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rad
1605 :
1606 : INTEGER :: nr
1607 : REAL(KIND=dp) :: oab, sab
1608 329160 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: expa, z
1609 :
1610 329160 : nr = SIZE(rad)
1611 1316640 : ALLOCATE (expa(nr), z(nr))
1612 :
1613 329160 : sab = SQRT(a + b)
1614 329160 : oab = dfac(ll + nu + 1)*rootpi/sab**(ll + 2)/2._dp**((ll + nu)/2 + 2)
1615 112946760 : z(:) = sab*rad(:)
1616 112946760 : erfa(:) = oab*erf(z(:))/z(:)**(nu + 1)
1617 112946760 : expa(:) = EXP(-z(:)**2)/sab**(ll + 2)/2._dp**((ll + nu)/2 + 2)
1618 0 : SELECT CASE (ll)
1619 : CASE DEFAULT
1620 : CALL cp_abort(__LOCATION__, &
1621 0 : "Only 0, 1, 2, 3, 4, 5, 6 are supported as the value of ll")
1622 : CASE (0)
1623 43941 : CPASSERT(nu == 0)
1624 : CASE (1)
1625 84522 : CPASSERT(nu == 1)
1626 28685322 : erfa(:) = erfa(:) - 6._dp*expa(:)/z(:)
1627 : CASE (2)
1628 78570 : SELECT CASE (nu)
1629 : CASE DEFAULT
1630 : CALL cp_abort(__LOCATION__, &
1631 0 : "Only 0, 2 are supported as the value of nu when ll = 2")
1632 : CASE (0)
1633 14043573 : erfa(:) = erfa(:) - 2._dp*expa(:)
1634 : CASE (2)
1635 28089327 : erfa(:) = erfa(:) - expa(:)*(20._dp + 30._dp/z(:)**2)
1636 : END SELECT
1637 : CASE (3)
1638 0 : SELECT CASE (nu)
1639 : CASE DEFAULT
1640 : CALL cp_abort(__LOCATION__, &
1641 0 : "Only 1, 3 are supported as the value of nu when ll = 3")
1642 : CASE (1)
1643 13744485 : erfa(:) = erfa(:) - expa(:)*(12._dp*z(:) + 30._dp/z(:))
1644 : CASE (3)
1645 13783770 : erfa(:) = erfa(:) - expa(:)*(56._dp*z(:) + 140._dp/z(:) + 210._dp/z(:)**3)
1646 : END SELECT
1647 : CASE (4)
1648 0 : SELECT CASE (nu)
1649 : CASE DEFAULT
1650 : CALL cp_abort(__LOCATION__, &
1651 0 : "Only 0, 2, 4 are supported as the value of nu when ll = 4")
1652 : CASE (0)
1653 0 : erfa(:) = erfa(:) - expa(:)*(4._dp*z(:)**2 + 14._dp)
1654 : CASE (2)
1655 0 : erfa(:) = erfa(:) - expa(:)*(40._dp*z(:)**2 + 140._dp + 210._dp/z(:)**2)
1656 : CASE (4)
1657 0 : erfa(:) = erfa(:) - expa(:)*(144._dp*z(:)**2 + 504._dp + 1260._dp/z(:)**2 + 1890._dp/z(:)**4)
1658 : END SELECT
1659 : CASE (5)
1660 0 : SELECT CASE (nu)
1661 : CASE DEFAULT
1662 : CALL cp_abort(__LOCATION__, &
1663 0 : "Only 1, 3, 5 are supported as the value of nu when ll = 5")
1664 : CASE (1)
1665 0 : erfa(:) = erfa(:) - expa(:)*(24._dp*z(:)**3 + 108._dp*z(:) + 210._dp/z(:))
1666 : CASE (3)
1667 0 : erfa(:) = erfa(:) - expa(:)*(112._dp*z(:)**3 + 504._dp*z(:) + 1260._dp/z(:) + 1890._dp/z(:)**3)
1668 : CASE (5)
1669 : erfa(:) = erfa(:) - expa(:)*(352._dp*z(:)**3 + 1584._dp*z(:) + 5544._dp/z(:) + &
1670 0 : 13860._dp/z(:)**3 + 20790._dp/z(:)**5)
1671 : END SELECT
1672 : CASE (6)
1673 329160 : SELECT CASE (nu)
1674 : CASE DEFAULT
1675 : CALL cp_abort(__LOCATION__, &
1676 0 : "Only 0, 2, 4, 6 are supported as the value of nu when ll = 6")
1677 : CASE (0)
1678 0 : erfa(:) = erfa(:) - expa(:)*(8._dp*z(:)**4 + 44._dp*z(:)**2 + 114._dp)
1679 : CASE (2)
1680 0 : erfa(:) = erfa(:) - expa(:)*(80._dp*z(:)**4 + 440._dp*z(:)**2 + 1260._dp + 1896._dp/z(:)**2)
1681 : CASE (4)
1682 : erfa(:) = erfa(:) - expa(:)*(288._dp*z(:)**4 + 1584._dp*z(:)**2 + 5544._dp + &
1683 0 : 13860._dp/z(:)**2 + 20790._dp/z(:)**4)
1684 : CASE (6)
1685 : erfa(:) = erfa(:) - expa(:)*(832._dp*z(:)**4 + 4576._dp*z(:)**2 + 20592._dp + &
1686 0 : 72072._dp/z(:)**2 + 180180._dp/z(:)**4 + 270270._dp/z(:)**6)
1687 : END SELECT
1688 : END SELECT
1689 :
1690 329160 : DEALLOCATE (expa, z)
1691 :
1692 329160 : END SUBROUTINE potential_coulomb_analytic
1693 :
1694 : ! **************************************************************************************************
1695 : !> \brief This routine calculates the longrange Coulomb potential of a product of two Gaussian with given angular momentum
1696 : !> \param erfa Array will contain the potential
1697 : !> \param a Exponent of first Gaussian charge distribution
1698 : !> \param b Exponent of second Gaussian charge distribution
1699 : !> \param ll Sum of angular momentum quantum numbers building one charge distribution
1700 : !> \param nu Angular momentum of interaction, (ll-nu) should be even.
1701 : !> \param rad Radial grid
1702 : !> \param omega Range-separation parameter
1703 : ! **************************************************************************************************
1704 353696 : PURE SUBROUTINE potential_longrange_analytic(erfa, a, b, ll, nu, rad, omega)
1705 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: erfa
1706 : REAL(KIND=dp), INTENT(IN) :: a, b
1707 : INTEGER, INTENT(IN) :: ll, nu
1708 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rad
1709 : REAL(KIND=dp), INTENT(IN) :: omega
1710 :
1711 : INTEGER :: i, lambda, nr
1712 : REAL(KIND=dp) :: ab, oab, pab, prel, sab
1713 353696 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: expa, z
1714 :
1715 353696 : IF (MOD(ll - nu, 2) == 0 .AND. ll >= nu .AND. nu >= 0) THEN
1716 353696 : nr = SIZE(rad)
1717 1414784 : ALLOCATE (expa(nr), z(nr))
1718 :
1719 353696 : lambda = (ll - nu)/2
1720 353696 : ab = a + b
1721 : sab = SQRT(ab)
1722 353696 : pab = omega**2*ab/(omega**2 + ab)
1723 353696 : prel = pab/ab
1724 119611296 : z(:) = SQRT(pab)*rad(:)
1725 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)
1726 119611296 : expa(:) = EXP(-z(:)**2)
1727 : lambda = (ll - nu)/2
1728 :
1729 353696 : IF (lambda > 0) THEN
1730 29379420 : erfa = 0.0_dp
1731 171640 : DO i = 1, lambda
1732 : erfa = erfa + (-prel)**i/REAL(i, KIND=dp)*binomial_gen(lambda + nu + 0.5_dp, lambda - i)* &
1733 29465240 : assoc_laguerre(z, REAL(nu, KIND=dp) + 0.5_dp, i - 1)
1734 : END DO
1735 29379420 : erfa = erfa*expa*z**nu
1736 :
1737 29379420 : erfa = erfa + 2.0_dp*binomial_gen(lambda + nu + 0.5_dp, lambda)*znFn(z, expa, nu)
1738 : ELSE
1739 90231876 : erfa = 2.0_dp*znFn(z, expa, nu)
1740 : END IF
1741 :
1742 119611296 : erfa = erfa*oab
1743 :
1744 353696 : DEALLOCATE (expa, z)
1745 : ELSE
1746 : ! Contribution to potential is zero (properties of spherical harmonics)
1747 0 : erfa = 0.0_dp
1748 : END IF
1749 :
1750 353696 : END SUBROUTINE potential_longrange_analytic
1751 :
1752 : ! **************************************************************************************************
1753 : !> \brief Boys' function times z**n
1754 : !> \param z ...
1755 : !> \param expa ...
1756 : !> \param n ...
1757 : !> \return ...
1758 : ! **************************************************************************************************
1759 119257600 : ELEMENTAL FUNCTION znFn(z, expa, n) RESULT(res)
1760 :
1761 : REAL(KIND=dp), INTENT(IN) :: z, expa
1762 : INTEGER, INTENT(IN) :: n
1763 : REAL(KIND=dp) :: res
1764 :
1765 : INTEGER :: i
1766 : REAL(KIND=dp) :: z_exp
1767 :
1768 119257600 : IF (n >= 0) THEN
1769 119257600 : IF (ABS(z) < 1.0E-20) THEN
1770 0 : res = z**n/(2.0_dp*REAL(n, KIND=dp) + 1.0_dp)
1771 119257600 : ELSE IF (n == 0) THEN
1772 30380000 : res = rootpi/2.0_dp*ERF(z)/z
1773 : ELSE
1774 88877600 : res = rootpi/4.0_dp*ERF(z)/z**2 - expa/z/2.0_dp
1775 88877600 : z_exp = -expa*0.5_dp
1776 :
1777 147420000 : DO i = 2, n
1778 58542400 : res = (REAL(i, KIND=dp) - 0.5_dp)*res/z + z_exp
1779 147420000 : z_exp = z_exp*z
1780 : END DO
1781 : END IF
1782 : ELSE ! Set it to zero (no Boys' function, to keep the ELEMENTAL keyword)
1783 : res = 0.0_dp
1784 : END IF
1785 :
1786 119257600 : END FUNCTION znFn
1787 :
1788 : ! **************************************************************************************************
1789 : !> \brief ...
1790 : !> \param z ...
1791 : !> \param a ...
1792 : !> \param n ...
1793 : !> \return ...
1794 : ! **************************************************************************************************
1795 29293600 : ELEMENTAL FUNCTION assoc_laguerre(z, a, n) RESULT(res)
1796 :
1797 : REAL(KIND=dp), INTENT(IN) :: z, a
1798 : INTEGER, INTENT(IN) :: n
1799 : REAL(KIND=dp) :: res
1800 :
1801 : INTEGER :: i
1802 : REAL(KIND=dp) :: f0, f1
1803 :
1804 29293600 : IF (n == 0) THEN
1805 : res = 1.0_dp
1806 0 : ELSE IF (n == 1) THEN
1807 0 : res = a + 1.0_dp - z
1808 0 : ELSE IF (n > 0) THEN
1809 0 : f0 = 1.0_dp
1810 0 : f1 = a + 1.0_dp - z
1811 :
1812 0 : DO i = 3, n
1813 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
1814 0 : f0 = f1
1815 0 : f1 = res
1816 : END DO
1817 : ELSE ! n is negative, set it zero (no polynomials, to keep the ELEMENTAL keyword)
1818 : res = 0.0_dp
1819 : END IF
1820 :
1821 29293600 : END FUNCTION assoc_laguerre
1822 :
1823 : ! **************************************************************************************************
1824 : !> \brief Compute Trace[opmat * pmat] == Trace[opmat * pmat^T] .
1825 : !> \param opmat operator matrix (e.g. Kohn-Sham matrix or overlap matrix)
1826 : !> \param pmat density matrix
1827 : !> \return value of trace
1828 : !> \par History
1829 : !> * 08.2008 created [Juerg Hutter]
1830 : ! **************************************************************************************************
1831 439604 : PURE FUNCTION atom_trace(opmat, pmat) RESULT(trace)
1832 : REAL(KIND=dp), DIMENSION(:, :, 0:), INTENT(IN) :: opmat, pmat
1833 : REAL(KIND=dp) :: trace
1834 :
1835 439604 : trace = accurate_dot_product(opmat, pmat)
1836 :
1837 439604 : END FUNCTION atom_trace
1838 :
1839 : ! **************************************************************************************************
1840 : !> \brief Calculate a potential matrix by integrating the potential on an atomic radial grid.
1841 : !> \param imat potential matrix
1842 : !> \param cpot potential on the atomic radial grid
1843 : !> \param basis atomic basis set
1844 : !> \param derivatives order of derivatives
1845 : !> \par History
1846 : !> * 08.2008 created [Juerg Hutter]
1847 : ! **************************************************************************************************
1848 208416 : SUBROUTINE numpot_matrix(imat, cpot, basis, derivatives)
1849 : REAL(KIND=dp), DIMENSION(:, :, 0:), INTENT(INOUT) :: imat
1850 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: cpot
1851 : TYPE(atom_basis_type), INTENT(INOUT) :: basis
1852 : INTEGER, INTENT(IN) :: derivatives
1853 :
1854 : INTEGER :: i, j, l, n
1855 :
1856 208416 : SELECT CASE (derivatives)
1857 : CASE (0)
1858 1280412 : DO l = 0, lmat
1859 1097496 : n = basis%nbas(l)
1860 3902082 : DO i = 1, n
1861 30808019 : DO j = i, n
1862 : imat(i, j, l) = imat(i, j, l) + &
1863 27088853 : integrate_grid(cpot, basis%bf(:, i, l), basis%bf(:, j, l), basis%grid)
1864 29710523 : imat(j, i, l) = imat(i, j, l)
1865 : END DO
1866 : END DO
1867 : END DO
1868 : CASE (1)
1869 177422 : DO l = 0, lmat
1870 152076 : n = basis%nbas(l)
1871 1241677 : DO i = 1, n
1872 15182941 : DO j = i, n
1873 : imat(i, j, l) = imat(i, j, l) + &
1874 13966610 : integrate_grid(cpot, basis%dbf(:, i, l), basis%bf(:, j, l), basis%grid)
1875 : imat(i, j, l) = imat(i, j, l) + &
1876 13966610 : integrate_grid(cpot, basis%bf(:, i, l), basis%dbf(:, j, l), basis%grid)
1877 15030865 : imat(j, i, l) = imat(i, j, l)
1878 : END DO
1879 : END DO
1880 : END DO
1881 : CASE (2)
1882 1078 : DO l = 0, lmat
1883 924 : n = basis%nbas(l)
1884 24234 : DO i = 1, n
1885 459894 : DO j = i, n
1886 : imat(i, j, l) = imat(i, j, l) + &
1887 435814 : integrate_grid(cpot, basis%dbf(:, i, l), basis%dbf(:, j, l), basis%grid)
1888 458970 : imat(j, i, l) = imat(i, j, l)
1889 : END DO
1890 : END DO
1891 : END DO
1892 : CASE DEFAULT
1893 208416 : CPABORT("Only 0, 1, 2 are supported as the value of derivatives")
1894 : END SELECT
1895 :
1896 208416 : END SUBROUTINE numpot_matrix
1897 :
1898 : ! **************************************************************************************************
1899 : !> \brief Contract Coulomb Electron Repulsion Integrals.
1900 : !> \param jmat ...
1901 : !> \param erint ...
1902 : !> \param pmat ...
1903 : !> \param nsize ...
1904 : !> \param all_nu ...
1905 : !> \par History
1906 : !> * 08.2008 created [Juerg Hutter]
1907 : ! **************************************************************************************************
1908 1515 : PURE SUBROUTINE ceri_contract(jmat, erint, pmat, nsize, all_nu)
1909 : REAL(KIND=dp), DIMENSION(:, :, 0:), INTENT(INOUT) :: jmat
1910 : TYPE(eri), DIMENSION(:), INTENT(IN) :: erint
1911 : REAL(KIND=dp), DIMENSION(:, :, 0:), INTENT(IN) :: pmat
1912 : INTEGER, DIMENSION(0:), INTENT(IN) :: nsize
1913 : LOGICAL, INTENT(IN), OPTIONAL :: all_nu
1914 :
1915 : INTEGER :: i1, i2, ij1, ij2, j1, j2, l1, l2, ll, &
1916 : n1, n2
1917 : LOGICAL :: have_all_nu
1918 : REAL(KIND=dp) :: eint, f1, f2
1919 :
1920 1515 : IF (PRESENT(all_nu)) THEN
1921 0 : have_all_nu = all_nu
1922 : ELSE
1923 : have_all_nu = .FALSE.
1924 : END IF
1925 :
1926 4237113 : jmat(:, :, :) = 0._dp
1927 : ll = 0
1928 10605 : DO l1 = 0, lmat
1929 9090 : n1 = nsize(l1)
1930 42420 : DO l2 = 0, l1
1931 31815 : n2 = nsize(l2)
1932 31815 : ll = ll + 1
1933 31815 : ij1 = 0
1934 558425 : DO i1 = 1, n1
1935 6043326 : DO j1 = i1, n1
1936 5484901 : ij1 = ij1 + 1
1937 5484901 : f1 = 1._dp
1938 5484901 : IF (i1 /= j1) f1 = 2._dp
1939 5484901 : ij2 = 0
1940 119694873 : DO i2 = 1, n2
1941 1403690542 : DO j2 = i2, n2
1942 1284522279 : ij2 = ij2 + 1
1943 1284522279 : f2 = 1._dp
1944 1284522279 : IF (i2 /= j2) f2 = 2._dp
1945 1284522279 : eint = erint(ll)%int(ij1, ij2)
1946 1398205641 : IF (l1 == l2) THEN
1947 418306136 : jmat(i1, j1, l1) = jmat(i1, j1, l1) + f2*pmat(i2, j2, l2)*eint
1948 : ELSE
1949 866216143 : jmat(i1, j1, l1) = jmat(i1, j1, l1) + f2*pmat(i2, j2, l2)*eint
1950 866216143 : jmat(i2, j2, l2) = jmat(i2, j2, l2) + f1*pmat(i1, j1, l1)*eint
1951 : END IF
1952 : END DO
1953 : END DO
1954 : END DO
1955 : END DO
1956 40905 : IF (have_all_nu) THEN
1957 : ! skip integral blocks with nu/=0
1958 0 : ll = ll + l2
1959 : END IF
1960 : END DO
1961 : END DO
1962 10605 : DO l1 = 0, lmat
1963 9090 : n1 = nsize(l1)
1964 170209 : DO i1 = 1, n1
1965 1873302 : DO j1 = i1, n1
1966 1864212 : jmat(j1, i1, l1) = jmat(i1, j1, l1)
1967 : END DO
1968 : END DO
1969 : END DO
1970 :
1971 1515 : END SUBROUTINE ceri_contract
1972 :
1973 : ! **************************************************************************************************
1974 : !> \brief Contract exchange Electron Repulsion Integrals.
1975 : !> \param kmat ...
1976 : !> \param erint ...
1977 : !> \param pmat ...
1978 : !> \param nsize ...
1979 : !> \par History
1980 : !> * 08.2008 created [Juerg Hutter]
1981 : ! **************************************************************************************************
1982 547 : PURE SUBROUTINE eeri_contract(kmat, erint, pmat, nsize)
1983 : REAL(KIND=dp), DIMENSION(:, :, 0:), INTENT(INOUT) :: kmat
1984 : TYPE(eri), DIMENSION(:), INTENT(IN) :: erint
1985 : REAL(KIND=dp), DIMENSION(:, :, 0:), INTENT(IN) :: pmat
1986 : INTEGER, DIMENSION(0:), INTENT(IN) :: nsize
1987 :
1988 : INTEGER :: i1, i2, ij1, ij2, j1, j2, l1, l2, lh, &
1989 : ll, n1, n2, nu
1990 : REAL(KIND=dp) :: almn, eint, f1, f2
1991 : REAL(KIND=dp), DIMENSION(0:maxfac) :: arho
1992 :
1993 547 : arho = 0._dp
1994 547 : DO ll = 0, maxfac, 2
1995 8752 : lh = ll/2
1996 8752 : arho(ll) = fac(ll)/fac(lh)**2
1997 : END DO
1998 :
1999 1470517 : kmat(:, :, :) = 0._dp
2000 : ll = 0
2001 3829 : DO l1 = 0, lmat
2002 3282 : n1 = nsize(l1)
2003 15316 : DO l2 = 0, l1
2004 11487 : n2 = nsize(l2)
2005 45401 : DO nu = ABS(l1 - l2), l1 + l2, 2
2006 30632 : almn = arho(-l1 + l2 + nu)*arho(l1 - l2 + nu)*arho(l1 + l2 - nu)/(REAL(l1 + l2 + nu + 1, dp)*arho(l1 + l2 + nu))
2007 30632 : almn = -0.5_dp*almn
2008 30632 : ll = ll + 1
2009 30632 : ij1 = 0
2010 504257 : DO i1 = 1, n1
2011 5332602 : DO j1 = i1, n1
2012 4839832 : ij1 = ij1 + 1
2013 4839832 : f1 = 1._dp
2014 4839832 : IF (i1 /= j1) f1 = 2._dp
2015 4839832 : ij2 = 0
2016 104653578 : DO i2 = 1, n2
2017 1197184274 : DO j2 = i2, n2
2018 1092992834 : ij2 = ij2 + 1
2019 1092992834 : f2 = 1._dp
2020 1092992834 : IF (i2 /= j2) f2 = 2._dp
2021 1092992834 : eint = erint(ll)%int(ij1, ij2)
2022 1192344442 : IF (l1 == l2) THEN
2023 430803648 : kmat(i1, j1, l1) = kmat(i1, j1, l1) + f2*almn*pmat(i2, j2, l2)*eint
2024 : ELSE
2025 662189186 : kmat(i1, j1, l1) = kmat(i1, j1, l1) + f2*almn*pmat(i2, j2, l2)*eint
2026 662189186 : kmat(i2, j2, l2) = kmat(i2, j2, l2) + f1*almn*pmat(i1, j1, l1)*eint
2027 : END IF
2028 : END DO
2029 : END DO
2030 : END DO
2031 : END DO
2032 : END DO
2033 : END DO
2034 : END DO
2035 3829 : DO l1 = 0, lmat
2036 3282 : n1 = nsize(l1)
2037 58734 : DO i1 = 1, n1
2038 647742 : DO j1 = i1, n1
2039 644460 : kmat(j1, i1, l1) = kmat(i1, j1, l1)
2040 : END DO
2041 : END DO
2042 : END DO
2043 :
2044 547 : END SUBROUTINE eeri_contract
2045 :
2046 : ! **************************************************************************************************
2047 : !> \brief Calculate the error matrix for each angular momentum.
2048 : !> \param emat error matrix
2049 : !> \param demax the largest absolute value of error matrix elements
2050 : !> \param kmat Kohn-Sham matrix
2051 : !> \param pmat electron density matrix
2052 : !> \param umat transformation matrix which reduce overlap matrix to its unitary form
2053 : !> \param upmat transformation matrix which reduce density matrix to its unitary form
2054 : !> \param nval number of linear-independent contracted basis functions
2055 : !> \param nbs number of contracted basis functions
2056 : !> \par History
2057 : !> * 08.2008 created [Juerg Hutter]
2058 : ! **************************************************************************************************
2059 84272 : PURE SUBROUTINE err_matrix(emat, demax, kmat, pmat, umat, upmat, nval, nbs)
2060 : REAL(KIND=dp), DIMENSION(:, :, 0:), INTENT(OUT) :: emat
2061 : REAL(KIND=dp), INTENT(OUT) :: demax
2062 : REAL(KIND=dp), DIMENSION(:, :, 0:), INTENT(IN) :: kmat, pmat, umat, upmat
2063 : INTEGER, DIMENSION(0:), INTENT(IN) :: nval, nbs
2064 :
2065 : INTEGER :: l, m, n
2066 84272 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: tkmat, tpmat
2067 :
2068 38934032 : emat = 0._dp
2069 589904 : DO l = 0, lmat
2070 505632 : n = nval(l)
2071 505632 : m = nbs(l)
2072 589904 : IF (m > 0) THEN
2073 1274040 : ALLOCATE (tkmat(1:m, 1:m), tpmat(1:m, 1:m))
2074 212340 : tkmat = 0._dp
2075 212340 : tpmat = 0._dp
2076 948153152 : 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)))
2077 948153152 : 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)))
2078 697200990 : tpmat(1:m, 1:m) = MATMUL(upmat(1:m, 1:m, l), MATMUL(tpmat(1:m, 1:m), upmat(1:m, 1:m, l)))
2079 :
2080 381391479 : 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))
2081 :
2082 212340 : DEALLOCATE (tkmat, tpmat)
2083 : END IF
2084 : END DO
2085 38934032 : demax = MAXVAL(ABS(emat))
2086 :
2087 84272 : END SUBROUTINE err_matrix
2088 :
2089 : ! **************************************************************************************************
2090 : !> \brief Calculate Slater density on a radial grid.
2091 : !> \param density1 alpha-spin electron density
2092 : !> \param density2 beta-spin electron density
2093 : !> \param zcore nuclear charge
2094 : !> \param state electronic state
2095 : !> \param grid atomic radial grid
2096 : !> \par History
2097 : !> * 06.2018 bugfix [Rustam Khaliullin]
2098 : !> * 02.2010 unrestricted KS and HF methods [Juerg Hutter]
2099 : !> * 12.2008 created [Juerg Hutter]
2100 : !> \note An initial electron density to guess atomic orbitals.
2101 : ! **************************************************************************************************
2102 12528 : SUBROUTINE slater_density(density1, density2, zcore, state, grid)
2103 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: density1, density2
2104 : INTEGER, INTENT(IN) :: zcore
2105 : TYPE(atom_state), INTENT(IN) :: state
2106 : TYPE(grid_atom_type), INTENT(IN) :: grid
2107 :
2108 : INTEGER :: counter, i, l, mc, mm(0:lmat), mo, n, ns
2109 : INTEGER, DIMENSION(lmat+1, 20) :: ne
2110 : REAL(KIND=dp) :: a, pf
2111 :
2112 : ! fill out a table with the total number of electrons
2113 : ! core + valence. format ne(l,n)
2114 12528 : ns = SIZE(ne, 2)
2115 12528 : ne = 0
2116 12528 : mm = get_maxn_occ(state%core)
2117 87696 : DO l = 0, lmat
2118 75168 : mo = state%maxn_occ(l)
2119 839376 : IF (SUM(state%core(l, :)) == 0) THEN
2120 66789 : CPASSERT(ns >= l + mo)
2121 82679 : DO counter = 1, mo
2122 82679 : ne(l + 1, l + counter) = NINT(state%occ(l, counter))
2123 : END DO
2124 : ELSE
2125 8379 : mc = mm(l) ! number of levels in the core
2126 22354 : CPASSERT(SUM(state%occ(l, 1:mc)) == 0)
2127 8379 : CPASSERT(ns >= l + mc)
2128 22354 : DO counter = 1, mc
2129 22354 : ne(l + 1, l + counter) = NINT(state%core(l, counter))
2130 : END DO
2131 8379 : CPASSERT(ns >= l + mc + mo)
2132 15778 : DO counter = mc + 1, mc + mo
2133 15778 : ne(l + 1, l + counter) = NINT(state%occ(l, counter))
2134 : END DO
2135 : END IF
2136 : END DO
2137 :
2138 5013238 : density1 = 0._dp
2139 5013238 : density2 = 0._dp
2140 33363 : DO l = 0, state%maxl_occ
2141 241713 : DO i = 1, SIZE(state%occ, 2)
2142 229185 : IF (state%occ(l, i) > 0._dp) THEN
2143 23659 : n = i + l
2144 23659 : a = srules(zcore, ne, n, l)
2145 23659 : pf = 1._dp/SQRT(fac(2*n))*(2._dp*a)**(n + 0.5_dp)
2146 23659 : IF (state%multiplicity == -1) THEN
2147 9369671 : density1(:) = density1(:) + state%occ(l, i)/fourpi*(grid%rad(:)**(n - 1)*EXP(-a*grid%rad(:))*pf)**2
2148 : ELSE
2149 81400 : density1(:) = density1(:) + state%occa(l, i)/fourpi*(grid%rad(:)**(n - 1)*EXP(-a*grid%rad(:))*pf)**2
2150 81400 : density2(:) = density2(:) + state%occb(l, i)/fourpi*(grid%rad(:)**(n - 1)*EXP(-a*grid%rad(:))*pf)**2
2151 : END IF
2152 : END IF
2153 : END DO
2154 : END DO
2155 :
2156 12528 : END SUBROUTINE slater_density
2157 :
2158 : ! **************************************************************************************************
2159 : !> \brief Calculate the functional derivative of the Wigner (correlation) - Slater (exchange)
2160 : !> density functional.
2161 : !> \param rho electron density on a radial grid
2162 : !> \param vxc (output) exchange-correlation potential
2163 : !> \par History
2164 : !> * 12.2008 created [Juerg Hutter]
2165 : !> \note A model XC-potential to guess atomic orbitals.
2166 : ! **************************************************************************************************
2167 12632 : PURE SUBROUTINE wigner_slater_functional(rho, vxc)
2168 : REAL(KIND=dp), DIMENSION(:), INTENT(IN) :: rho
2169 : REAL(KIND=dp), DIMENSION(:), INTENT(OUT) :: vxc
2170 :
2171 : INTEGER :: i
2172 : REAL(KIND=dp) :: ec, ex, rs, vc, vx
2173 :
2174 5055742 : vxc = 0._dp
2175 5055742 : DO i = 1, SIZE(rho)
2176 5055742 : IF (rho(i) > 1.e-20_dp) THEN
2177 : ! 3/4 * (3/pi)^{1/3} == 0.7385588
2178 4732413 : ex = -0.7385588_dp*rho(i)**0.333333333_dp
2179 4732413 : vx = 1.333333333_dp*ex
2180 4732413 : rs = (3._dp/fourpi/rho(i))**0.333333333_dp
2181 4732413 : ec = -0.88_dp/(rs + 7.8_dp)
2182 4732413 : vc = ec*(1._dp + rs/(3._dp*(rs + 7.8_dp)))
2183 4732413 : vxc(i) = vx + vc
2184 : END IF
2185 : END DO
2186 :
2187 12632 : END SUBROUTINE wigner_slater_functional
2188 :
2189 : ! **************************************************************************************************
2190 : !> \brief Check that the atomic multiplicity is consistent with the electronic structure method.
2191 : !> \param method electronic structure method
2192 : !> \param multiplicity multiplicity
2193 : !> \return consistency status
2194 : !> \par History
2195 : !> * 11.2009 unrestricted KS and HF methods [Juerg Hutter]
2196 : ! **************************************************************************************************
2197 3083 : PURE FUNCTION atom_consistent_method(method, multiplicity) RESULT(consistent)
2198 : INTEGER, INTENT(IN) :: method, multiplicity
2199 : LOGICAL :: consistent
2200 :
2201 : ! multiplicity == -1 means it has not been specified explicitly;
2202 : ! see the source code of the subroutine atom_set_occupation() for further details.
2203 3083 : SELECT CASE (method)
2204 : CASE DEFAULT
2205 1810 : consistent = .FALSE.
2206 : CASE (do_rks_atom)
2207 1810 : consistent = (multiplicity == -1)
2208 : CASE (do_uks_atom)
2209 93 : consistent = (multiplicity /= -1)
2210 : CASE (do_rhf_atom)
2211 1172 : consistent = (multiplicity == -1)
2212 : CASE (do_uhf_atom)
2213 8 : consistent = (multiplicity /= -1)
2214 : CASE (do_rohf_atom)
2215 3083 : consistent = .FALSE.
2216 : END SELECT
2217 :
2218 3083 : END FUNCTION atom_consistent_method
2219 :
2220 : ! **************************************************************************************************
2221 : !> \brief Calculate the total electron density at R=0.
2222 : !> \param atom information about the atomic kind
2223 : !> \param rho0 (output) computed density
2224 : !> \par History
2225 : !> * 05.2016 created [Juerg Hutter]
2226 : ! **************************************************************************************************
2227 2554 : SUBROUTINE get_rho0(atom, rho0)
2228 : TYPE(atom_type), INTENT(IN) :: atom
2229 : REAL(KIND=dp), INTENT(OUT) :: rho0
2230 :
2231 : INTEGER :: m0, m1, m2, method, nr
2232 : LOGICAL :: nlcc, spinpol
2233 : REAL(KIND=dp) :: d0, d1, d2, r0, r1, r2, w0, w1
2234 2554 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: xfun
2235 2554 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: rho
2236 :
2237 2554 : method = atom%method_type
2238 : SELECT CASE (method)
2239 : CASE (do_rks_atom, do_rhf_atom)
2240 48 : spinpol = .FALSE.
2241 : CASE (do_uks_atom, do_uhf_atom)
2242 48 : spinpol = .TRUE.
2243 : CASE (do_rohf_atom)
2244 0 : CPABORT("ROHF not yet implemented for get_rho0")
2245 : CASE DEFAULT
2246 2554 : CPABORT("Unknown method for get_rho0")
2247 : END SELECT
2248 :
2249 2554 : nr = atom%basis%grid%nr
2250 2554 : nlcc = .FALSE.
2251 2554 : IF (atom%potential%ppot_type == gth_pseudo) THEN
2252 2103 : nlcc = atom%potential%gth_pot%nlcc
2253 451 : ELSE IF (atom%potential%ppot_type == upf_pseudo) THEN
2254 2 : nlcc = atom%potential%upf_pot%core_correction
2255 449 : ELSE IF (atom%potential%ppot_type == sgp_pseudo) THEN
2256 0 : nlcc = atom%potential%sgp_pot%has_nlcc
2257 : END IF
2258 2105 : IF (nlcc) THEN
2259 33 : ALLOCATE (xfun(nr))
2260 : END IF
2261 :
2262 1028164 : m0 = MAXVAL(MINLOC(atom%basis%grid%rad))
2263 2554 : IF (m0 == nr) THEN
2264 2554 : m1 = m0 - 1
2265 2554 : m2 = m0 - 2
2266 0 : ELSE IF (m0 == 1) THEN
2267 : m1 = 2
2268 : m2 = 3
2269 : ELSE
2270 0 : CPABORT("GRID Definition incompatible")
2271 : END IF
2272 2554 : r0 = atom%basis%grid%rad(m0)
2273 2554 : r1 = atom%basis%grid%rad(m1)
2274 2554 : r2 = atom%basis%grid%rad(m2)
2275 2554 : w0 = r1/(r1 - r0)
2276 2554 : w1 = 1 - w0
2277 :
2278 2554 : IF (spinpol) THEN
2279 144 : ALLOCATE (rho(nr, 2))
2280 48 : CALL atom_density(rho(:, 1), atom%orbitals%pmata, atom%basis, atom%state%maxl_occ, typ="RHO")
2281 48 : CALL atom_density(rho(:, 2), atom%orbitals%pmatb, atom%basis, atom%state%maxl_occ, typ="RHO")
2282 48 : IF (nlcc) THEN
2283 1 : xfun = 0.0_dp
2284 1 : CALL atom_core_density(xfun(:), atom%potential, typ="RHO", rr=atom%basis%grid%rad)
2285 401 : rho(:, 1) = rho(:, 1) + 0.5_dp*xfun(:)
2286 401 : rho(:, 2) = rho(:, 2) + 0.5_dp*xfun(:)
2287 : END IF
2288 19648 : rho(:, 1) = rho(:, 1) + rho(:, 2)
2289 : ELSE
2290 7518 : ALLOCATE (rho(nr, 1))
2291 2506 : CALL atom_density(rho(:, 1), atom%orbitals%pmat, atom%basis, atom%state%maxl_occ, typ="RHO")
2292 2506 : IF (nlcc) THEN
2293 10 : CALL atom_core_density(rho(:, 1), atom%potential, typ="RHO", rr=atom%basis%grid%rad)
2294 : END IF
2295 : END IF
2296 2554 : d0 = rho(m0, 1)
2297 2554 : d1 = rho(m1, 1)
2298 2554 : d2 = rho(m2, 1)
2299 :
2300 2554 : rho0 = w0*d0 + w1*d1
2301 2554 : rho0 = MAX(rho0, 0.0_dp)
2302 :
2303 2554 : DEALLOCATE (rho)
2304 2554 : IF (nlcc) THEN
2305 11 : DEALLOCATE (xfun)
2306 : END IF
2307 :
2308 2554 : END SUBROUTINE get_rho0
2309 :
2310 : ! **************************************************************************************************
2311 : !> \brief Print condition numbers of the given atomic basis set.
2312 : !> \param basis atomic basis set
2313 : !> \param crad atomic covalent radius
2314 : !> \param iw output file unit
2315 : !> \par History
2316 : !> * 05.2016 created [Juerg Hutter]
2317 : ! **************************************************************************************************
2318 5 : SUBROUTINE atom_condnumber(basis, crad, iw)
2319 : TYPE(atom_basis_type), POINTER :: basis
2320 : REAL(KIND=dp) :: crad
2321 : INTEGER, INTENT(IN) :: iw
2322 :
2323 : INTEGER :: i
2324 : REAL(KIND=dp) :: ci
2325 : REAL(KIND=dp), DIMENSION(10) :: cnum, rad
2326 :
2327 5 : WRITE (iw, '(/,A,F8.4)') " Basis Set Condition Numbers: 2*covalent radius=", 2*crad
2328 5 : CALL init_orbital_pointers(lmat)
2329 5 : CALL init_spherical_harmonics(lmat, 0)
2330 5 : cnum = 0.0_dp
2331 50 : DO i = 1, 9
2332 45 : ci = 2.0_dp*(0.85_dp + i*0.05_dp)
2333 45 : rad(i) = crad*ci
2334 45 : CALL atom_basis_condnum(basis, rad(i), cnum(i))
2335 45 : WRITE (iw, '(A,F15.3,T50,A,F14.4)') " Lattice constant:", &
2336 95 : rad(i), "Condition number:", cnum(i)
2337 : END DO
2338 5 : rad(10) = 0.01_dp
2339 5 : CALL atom_basis_condnum(basis, rad(10), cnum(10))
2340 5 : WRITE (iw, '(A,A,T50,A,F14.4)') " Lattice constant:", &
2341 10 : " Inf", "Condition number:", cnum(i)
2342 5 : CALL deallocate_orbital_pointers
2343 5 : CALL deallocate_spherical_harmonics
2344 :
2345 5 : END SUBROUTINE atom_condnumber
2346 :
2347 : ! **************************************************************************************************
2348 : !> \brief Estimate completeness of the given atomic basis set.
2349 : !> \param basis atomic basis set
2350 : !> \param zv atomic number
2351 : !> \param iw output file unit
2352 : ! **************************************************************************************************
2353 5 : SUBROUTINE atom_completeness(basis, zv, iw)
2354 : TYPE(atom_basis_type), POINTER :: basis
2355 : INTEGER, INTENT(IN) :: zv, iw
2356 :
2357 : INTEGER :: i, j, l, ll, m, n, nbas, nl, nr
2358 : INTEGER, DIMENSION(0:lmat) :: nelem, nlmax, nlmin
2359 : INTEGER, DIMENSION(lmat+1, 7) :: ne
2360 : REAL(KIND=dp) :: al, c1, c2, pf
2361 5 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: sfun
2362 5 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: bmat
2363 5 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :) :: omat
2364 5 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :, :, :) :: sint
2365 : REAL(KIND=dp), DIMENSION(0:lmat, 10) :: snl
2366 : REAL(KIND=dp), DIMENSION(2) :: sse
2367 : REAL(KIND=dp), DIMENSION(lmat+1, 7) :: sexp
2368 :
2369 5 : ne = 0
2370 5 : nelem = 0
2371 25 : nelem(0:3) = ptable(zv)%e_conv(0:3)
2372 35 : DO l = 0, lmat
2373 30 : ll = 2*(2*l + 1)
2374 64 : DO i = 1, 7
2375 89 : IF (nelem(l) >= ll) THEN
2376 22 : ne(l + 1, i) = ll
2377 22 : nelem(l) = nelem(l) - ll
2378 37 : ELSE IF (nelem(l) > 0) THEN
2379 7 : ne(l + 1, i) = nelem(l)
2380 7 : nelem(l) = 0
2381 : ELSE
2382 : EXIT
2383 : END IF
2384 : END DO
2385 : END DO
2386 :
2387 35 : nlmin = 1
2388 35 : nlmax = 1
2389 35 : DO l = 0, lmat
2390 30 : nlmin(l) = l + 1
2391 240 : DO i = 1, 7
2392 240 : IF (ne(l + 1, i) > 0) THEN
2393 29 : nlmax(l) = i + l
2394 : END IF
2395 : END DO
2396 35 : nlmax(l) = MAX(nlmax(l), nlmin(l) + 1)
2397 : END DO
2398 :
2399 : ! Slater exponents
2400 : sexp = 0.0_dp
2401 35 : DO l = 0, lmat
2402 30 : sse(1) = 0.05_dp
2403 30 : sse(2) = 10.0_dp
2404 165 : DO i = l + 1, 7
2405 135 : sexp(l + 1, i) = srules(zv, ne, i, l)
2406 165 : IF (ne(l + 1, i - l) > 0) THEN
2407 29 : sse(1) = MAX(sse(1), sexp(l + 1, i))
2408 29 : sse(2) = MIN(sse(2), sexp(l + 1, i))
2409 : END IF
2410 : END DO
2411 335 : DO i = 1, 10
2412 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))
2413 : END DO
2414 : END DO
2415 :
2416 35 : nbas = MAXVAL(basis%nbas)
2417 25 : ALLOCATE (omat(nbas, nbas, 0:lmat))
2418 5 : nr = SIZE(basis%bf, 1)
2419 30 : ALLOCATE (sfun(nr), sint(10, 2, nbas, 0:lmat))
2420 5 : sint = 0._dp
2421 :
2422 : ! calculate overlaps between test functions and basis
2423 35 : DO l = 0, lmat
2424 335 : DO i = 1, 10
2425 300 : al = snl(l, i)
2426 300 : nl = nlmin(l)
2427 300 : pf = (2._dp*al)**nl*SQRT(2._dp*al/fac(2*nl))
2428 120300 : sfun(1:nr) = pf*basis%grid%rad(1:nr)**(nl - 1)*EXP(-al*basis%grid%rad(1:nr))
2429 1500 : DO j = 1, basis%nbas(l)
2430 481500 : sint(i, 1, j, l) = SUM(sfun(1:nr)*basis%bf(1:nr, j, l)*basis%grid%wr(1:nr))
2431 : END DO
2432 300 : nl = nlmax(l)
2433 300 : pf = (2._dp*al)**nl*SQRT(2._dp*al/fac(2*nl))
2434 120300 : sfun(1:nr) = pf*basis%grid%rad(1:nr)**(nl - 1)*EXP(-al*basis%grid%rad(1:nr))
2435 1530 : DO j = 1, basis%nbas(l)
2436 481500 : sint(i, 2, j, l) = SUM(sfun(1:nr)*basis%bf(1:nr, j, l)*basis%grid%wr(1:nr))
2437 : END DO
2438 : END DO
2439 : END DO
2440 :
2441 35 : DO l = 0, lmat
2442 30 : n = basis%nbas(l)
2443 30 : IF (n <= 0) CYCLE
2444 13 : m = basis%nprim(l)
2445 13 : SELECT CASE (basis%basis_type)
2446 : CASE DEFAULT
2447 0 : CPABORT("Unknown basis type for atom_completeness")
2448 : CASE (GTO_BASIS)
2449 10 : CALL sg_overlap(omat(1:n, 1:n, l), l, basis%am(1:n, l), basis%am(1:n, l))
2450 : CASE (CGTO_BASIS)
2451 12 : ALLOCATE (bmat(m, m))
2452 3 : CALL sg_overlap(bmat(1:m, 1:m), l, basis%am(1:m, l), basis%am(1:m, l))
2453 3 : CALL contract2(omat(1:n, 1:n, l), bmat(1:m, 1:m), basis%cm(1:m, 1:n, l))
2454 3 : DEALLOCATE (bmat)
2455 : CASE (STO_BASIS)
2456 : CALL sto_overlap(omat(1:n, 1:n, l), basis%ns(1:n, l), basis%as(1:n, l), &
2457 0 : basis%ns(1:n, l), basis%as(1:n, l))
2458 : CASE (NUM_BASIS)
2459 13 : CPABORT("Numerical basis not yet implemented for atom_completeness")
2460 : END SELECT
2461 35 : CALL invmat_symm(omat(1:n, 1:n, l))
2462 : END DO
2463 :
2464 5 : WRITE (iw, '(/,A)') " Basis Set Completeness Estimates"
2465 35 : DO l = 0, lmat
2466 30 : n = basis%nbas(l)
2467 30 : IF (n <= 0) CYCLE
2468 13 : WRITE (iw, '(A,I3)') " L-quantum number: ", l
2469 13 : WRITE (iw, '(A,T31,A,I2,T61,A,I2)') " Slater Exponent", "Completeness n-qm=", nlmin(l), &
2470 26 : "Completeness n-qm=", nlmax(l)
2471 148 : DO i = 10, 1, -1
2472 20790 : c1 = DOT_PRODUCT(sint(i, 1, 1:n, l), MATMUL(omat(1:n, 1:n, l), sint(i, 1, 1:n, l)))
2473 20790 : c2 = DOT_PRODUCT(sint(i, 2, 1:n, l), MATMUL(omat(1:n, 1:n, l), sint(i, 2, 1:n, l)))
2474 160 : WRITE (iw, "(T6,F14.6,T41,F10.6,T71,F10.6)") snl(l, i), c1, c2
2475 : END DO
2476 : END DO
2477 :
2478 5 : DEALLOCATE (omat, sfun, sint)
2479 :
2480 5 : END SUBROUTINE atom_completeness
2481 :
2482 : ! **************************************************************************************************
2483 : !> \brief Calculate the condition number of the given atomic basis set.
2484 : !> \param basis atomic basis set
2485 : !> \param rad lattice constant (e.g. doubled atomic covalent radius)
2486 : !> \param cnum (output) condition number
2487 : ! **************************************************************************************************
2488 104 : SUBROUTINE atom_basis_condnum(basis, rad, cnum)
2489 : TYPE(atom_basis_type) :: basis
2490 : REAL(KIND=dp), INTENT(IN) :: rad
2491 : REAL(KIND=dp), INTENT(OUT) :: cnum
2492 :
2493 : INTEGER :: ia, ib, imax, info, ix, iy, iz, ja, jb, &
2494 : ka, kb, l, la, lb, lwork, na, nb, &
2495 : nbas, nna, nnb
2496 104 : INTEGER, ALLOCATABLE, DIMENSION(:, :) :: ibptr
2497 : REAL(KIND=dp) :: r1, r2, reps, rmax
2498 104 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:) :: weig, work
2499 104 : REAL(KIND=dp), ALLOCATABLE, DIMENSION(:, :) :: smat
2500 : REAL(KIND=dp), DIMENSION(2*lmat+1, 2*lmat+1) :: sab
2501 : REAL(KIND=dp), DIMENSION(3) :: rab
2502 104 : REAL(KIND=dp), DIMENSION(:), POINTER :: zeta, zetb
2503 :
2504 : ! Number of spherical Gaussian orbitals with angular momentum lmat: nso(lmat) = 2*lmat+1
2505 :
2506 : ! total number of basis functions
2507 104 : nbas = 0
2508 728 : DO l = 0, lmat
2509 728 : nbas = nbas + basis%nbas(l)*(2*l + 1)
2510 : END DO
2511 :
2512 624 : ALLOCATE (smat(nbas, nbas), ibptr(nbas, 0:lmat))
2513 104 : smat = 0.0_dp
2514 104 : ibptr = 0
2515 104 : na = 0
2516 728 : DO l = 0, lmat
2517 2276 : DO ia = 1, basis%nbas(l)
2518 1548 : ibptr(ia, l) = na + 1
2519 2172 : na = na + (2*l + 1)
2520 : END DO
2521 : END DO
2522 :
2523 104 : reps = 1.e-14_dp
2524 104 : IF (basis%basis_type == GTO_BASIS .OR. &
2525 : basis%basis_type == CGTO_BASIS) THEN
2526 728 : DO la = 0, lmat
2527 624 : na = basis%nprim(la)
2528 624 : nna = 2*la + 1
2529 624 : IF (na == 0) CYCLE
2530 238 : zeta => basis%am(:, la)
2531 1770 : DO lb = 0, lmat
2532 1428 : nb = basis%nprim(lb)
2533 1428 : nnb = 2*lb + 1
2534 1428 : IF (nb == 0) CYCLE
2535 566 : zetb => basis%am(:, lb)
2536 5768 : DO ia = 1, na
2537 56004 : DO ib = 1, nb
2538 49998 : IF (rad < 0.1_dp) THEN
2539 : imax = 0
2540 : ELSE
2541 46083 : r1 = exp_radius(la, zeta(ia), reps, 1.0_dp)
2542 46083 : r2 = exp_radius(lb, zetb(ib), reps, 1.0_dp)
2543 46083 : rmax = MAX(2._dp*r1, 2._dp*r2)
2544 46083 : imax = INT(rmax/rad) + 1
2545 : END IF
2546 50661 : IF (imax > 1) THEN
2547 36030 : CALL overlap_ab_sp(la, zeta(ia), lb, zetb(ib), rad, sab)
2548 36030 : IF (basis%basis_type == GTO_BASIS) THEN
2549 24238 : ja = ibptr(ia, la)
2550 24238 : jb = ibptr(ib, lb)
2551 194898 : smat(ja:ja + nna - 1, jb:jb + nnb - 1) = smat(ja:ja + nna - 1, jb:jb + nnb - 1) + sab(1:nna, 1:nnb)
2552 11792 : ELSEIF (basis%basis_type == CGTO_BASIS) THEN
2553 48574 : DO ka = 1, basis%nbas(la)
2554 181930 : DO kb = 1, basis%nbas(lb)
2555 133356 : ja = ibptr(ka, la)
2556 133356 : jb = ibptr(kb, lb)
2557 : smat(ja:ja + nna - 1, jb:jb + nnb - 1) = smat(ja:ja + nna - 1, jb:jb + nnb - 1) + &
2558 899042 : sab(1:nna, 1:nnb)*basis%cm(ia, ka, la)*basis%cm(ib, kb, lb)
2559 : END DO
2560 : END DO
2561 : END IF
2562 : ELSE
2563 48042 : DO ix = -imax, imax
2564 34074 : rab(1) = rad*ix
2565 142434 : DO iy = -imax, imax
2566 94392 : rab(2) = rad*iy
2567 403812 : DO iz = -imax, imax
2568 275346 : rab(3) = rad*iz
2569 275346 : CALL overlap_ab_s(la, zeta(ia), lb, zetb(ib), rab, sab)
2570 369738 : IF (basis%basis_type == GTO_BASIS) THEN
2571 269262 : ja = ibptr(ia, la)
2572 269262 : jb = ibptr(ib, lb)
2573 : smat(ja:ja + nna - 1, jb:jb + nnb - 1) = smat(ja:ja + nna - 1, jb:jb + nnb - 1) &
2574 1510034 : + sab(1:nna, 1:nnb)
2575 6084 : ELSEIF (basis%basis_type == CGTO_BASIS) THEN
2576 24030 : DO ka = 1, basis%nbas(la)
2577 79902 : DO kb = 1, basis%nbas(lb)
2578 55872 : ja = ibptr(ka, la)
2579 55872 : jb = ibptr(kb, lb)
2580 : smat(ja:ja + nna - 1, jb:jb + nnb - 1) = &
2581 : smat(ja:ja + nna - 1, jb:jb + nnb - 1) + &
2582 252778 : sab(1:nna, 1:nnb)*basis%cm(ia, ka, la)*basis%cm(ib, kb, lb)
2583 : END DO
2584 : END DO
2585 : END IF
2586 : END DO
2587 : END DO
2588 : END DO
2589 : END IF
2590 : END DO
2591 : END DO
2592 : END DO
2593 : END DO
2594 : ELSE
2595 0 : CPABORT("Condition number not available for this basis type")
2596 : END IF
2597 :
2598 104 : info = 0
2599 104 : lwork = MAX(nbas*nbas, nbas + 100)
2600 520 : ALLOCATE (weig(nbas), work(lwork))
2601 :
2602 104 : CALL dsyev("N", "U", nbas, smat, nbas, weig, work, lwork, info)
2603 104 : CPASSERT(info == 0)
2604 104 : IF (weig(1) < 0.0_dp) THEN
2605 14 : cnum = 100._dp
2606 : ELSE
2607 90 : cnum = ABS(weig(nbas)/weig(1))
2608 90 : cnum = LOG10(cnum)
2609 : END IF
2610 :
2611 104 : DEALLOCATE (smat, ibptr, weig, work)
2612 :
2613 104 : END SUBROUTINE atom_basis_condnum
2614 :
2615 : ! **************************************************************************************************
2616 : !> \brief Transform a matrix expressed in terms of a uncontracted basis set to a contracted one.
2617 : !> \param int (output) contracted matrix
2618 : !> \param omat uncontracted matrix
2619 : !> \param cm matrix of contraction coefficients
2620 : ! **************************************************************************************************
2621 76732 : SUBROUTINE contract2(int, omat, cm)
2622 : REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: int
2623 : REAL(dp), DIMENSION(:, :), INTENT(IN) :: omat, cm
2624 :
2625 : CHARACTER(len=*), PARAMETER :: routineN = 'contract2'
2626 :
2627 : INTEGER :: handle, m, n
2628 :
2629 76732 : CALL timeset(routineN, handle)
2630 :
2631 76732 : n = SIZE(int, 1)
2632 76732 : m = SIZE(omat, 1)
2633 :
2634 8877022 : INT(1:n, 1:n) = MATMUL(TRANSPOSE(cm(1:m, 1:n)), MATMUL(omat(1:m, 1:m), cm(1:m, 1:n)))
2635 :
2636 76732 : CALL timestop(handle)
2637 :
2638 76732 : END SUBROUTINE contract2
2639 :
2640 : ! **************************************************************************************************
2641 : !> \brief Same as contract2(), but add the new contracted matrix to the old one
2642 : !> instead of overwriting it.
2643 : !> \param int (input/output) contracted matrix to be updated
2644 : !> \param omat uncontracted matrix
2645 : !> \param cm matrix of contraction coefficients
2646 : ! **************************************************************************************************
2647 38208 : SUBROUTINE contract2add(int, omat, cm)
2648 : REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: int
2649 : REAL(dp), DIMENSION(:, :), INTENT(IN) :: omat, cm
2650 :
2651 : CHARACTER(len=*), PARAMETER :: routineN = 'contract2add'
2652 :
2653 : INTEGER :: handle, m, n
2654 :
2655 38208 : CALL timeset(routineN, handle)
2656 :
2657 38208 : n = SIZE(int, 1)
2658 38208 : m = SIZE(omat, 1)
2659 :
2660 3265676 : 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)))
2661 :
2662 38208 : CALL timestop(handle)
2663 :
2664 38208 : END SUBROUTINE contract2add
2665 :
2666 : ! **************************************************************************************************
2667 : !> \brief Contract a matrix of Electron Repulsion Integrals (ERI-s).
2668 : !> \param eri (output) contracted ERI-s
2669 : !> \param omat uncontracted ERI-s
2670 : !> \param cm1 first matrix of contraction coefficients
2671 : !> \param cm2 second matrix of contraction coefficients
2672 : ! **************************************************************************************************
2673 1232 : SUBROUTINE contract4(eri, omat, cm1, cm2)
2674 : REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: eri
2675 : REAL(dp), DIMENSION(:, :), INTENT(IN) :: omat, cm1, cm2
2676 :
2677 : CHARACTER(len=*), PARAMETER :: routineN = 'contract4'
2678 :
2679 : INTEGER :: handle, i1, i2, m1, m2, mm1, mm2, n1, &
2680 : n2, nn1, nn2
2681 1232 : REAL(dp), ALLOCATABLE, DIMENSION(:, :) :: amat, atran, bmat, btran, hint
2682 :
2683 1232 : CALL timeset(routineN, handle)
2684 :
2685 1232 : m1 = SIZE(cm1, 1)
2686 1232 : n1 = SIZE(cm1, 2)
2687 1232 : m2 = SIZE(cm2, 1)
2688 1232 : n2 = SIZE(cm2, 2)
2689 1232 : nn1 = SIZE(eri, 1)
2690 1232 : nn2 = SIZE(eri, 2)
2691 1232 : mm1 = SIZE(omat, 1)
2692 1232 : mm2 = SIZE(omat, 2)
2693 :
2694 7680 : ALLOCATE (amat(m1, m1), atran(n1, n1), bmat(m2, m2), btran(n2, n2))
2695 2844 : ALLOCATE (hint(mm1, nn2))
2696 :
2697 1894 : DO i1 = 1, mm1
2698 662 : CALL iunpack(bmat(1:m2, 1:m2), omat(i1, 1:mm2), m2)
2699 662 : CALL contract2(btran(1:n2, 1:n2), bmat(1:m2, 1:m2), cm2(1:m2, 1:n2))
2700 1894 : CALL ipack(btran(1:n2, 1:n2), hint(i1, 1:nn2), n2)
2701 : END DO
2702 :
2703 2330 : DO i2 = 1, nn2
2704 1098 : CALL iunpack(amat(1:m1, 1:m1), hint(1:mm1, i2), m1)
2705 1098 : CALL contract2(atran(1:n1, 1:n1), amat(1:m1, 1:m1), cm1(1:m1, 1:n1))
2706 2330 : CALL ipack(atran(1:n1, 1:n1), eri(1:nn1, i2), n1)
2707 : END DO
2708 :
2709 1232 : DEALLOCATE (amat, atran, bmat, btran)
2710 1232 : DEALLOCATE (hint)
2711 :
2712 1232 : CALL timestop(handle)
2713 :
2714 1232 : END SUBROUTINE contract4
2715 :
2716 : ! **************************************************************************************************
2717 : !> \brief Pack an n x n square real matrix into a 1-D vector.
2718 : !> \param mat matrix to pack
2719 : !> \param vec resulting vector
2720 : !> \param n size of the matrix
2721 : ! **************************************************************************************************
2722 1760 : PURE SUBROUTINE ipack(mat, vec, n)
2723 : REAL(dp), DIMENSION(:, :), INTENT(IN) :: mat
2724 : REAL(dp), DIMENSION(:), INTENT(INOUT) :: vec
2725 : INTEGER, INTENT(IN) :: n
2726 :
2727 : INTEGER :: i, ij, j
2728 :
2729 1760 : ij = 0
2730 4500 : DO i = 1, n
2731 10110 : DO j = i, n
2732 5610 : ij = ij + 1
2733 8350 : vec(ij) = mat(i, j)
2734 : END DO
2735 : END DO
2736 :
2737 1760 : END SUBROUTINE ipack
2738 :
2739 : ! **************************************************************************************************
2740 : !> \brief Unpack a 1-D real vector as a n x n square matrix.
2741 : !> \param mat resulting matrix
2742 : !> \param vec vector to unpack
2743 : !> \param n size of the matrix
2744 : ! **************************************************************************************************
2745 1760 : PURE SUBROUTINE iunpack(mat, vec, n)
2746 : REAL(dp), DIMENSION(:, :), INTENT(INOUT) :: mat
2747 : REAL(dp), DIMENSION(:), INTENT(IN) :: vec
2748 : INTEGER, INTENT(IN) :: n
2749 :
2750 : INTEGER :: i, ij, j
2751 :
2752 1760 : ij = 0
2753 6541 : DO i = 1, n
2754 23252 : DO j = i, n
2755 16711 : ij = ij + 1
2756 16711 : mat(i, j) = vec(ij)
2757 21492 : mat(j, i) = vec(ij)
2758 : END DO
2759 : END DO
2760 :
2761 1760 : END SUBROUTINE iunpack
2762 :
2763 1127200 : END MODULE atom_utils
|
