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