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