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