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