53 #include "./base/base_uses.f90"
58 TYPE(atom_basis_type),
POINTER :: basis => null()
64 CHARACTER(len=*),
PARAMETER,
PRIVATE :: modulen =
'atom_grb'
79 TYPE(atom_p_type),
DIMENSION(:, :),
POINTER :: atom_info
80 TYPE(section_vals_type),
POINTER :: atom_section
81 INTEGER,
INTENT(IN) :: iw
83 CHARACTER(len=default_string_length) :: abas, basname
84 CHARACTER(len=default_string_length),
DIMENSION(1) :: basline
85 CHARACTER(len=default_string_length),
DIMENSION(3) :: headline
86 INTEGER :: i, ider, is, iunit, j, k, l, lhomo, ll, &
87 lval, m, maxl, mb, method, mo, n, &
88 nder, ngp, nhomo, nr, num_gto, &
89 num_pol, quadtype, s1, s2
90 INTEGER,
DIMENSION(0:7) :: nbas
91 INTEGER,
DIMENSION(0:lmat) :: next_bas, next_prim
92 INTEGER,
DIMENSION(:),
POINTER :: num_bas
93 REAL(kind=
dp) :: al, amin, aval, cnum, crad, cradx, cval, delta, dene, ear, emax, &
94 energy_ex(0:
lmat), energy_ref, energy_vb(0:
lmat), expzet, fhomo, o, prefac, rconf, rk, &
95 rmax, scon, zeta, zval
96 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: ale, alp, rho
97 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: amat
98 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :) :: ebasis, pbasis, qbasis, rbasis
99 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :, :, :) :: wfn
100 REAL(kind=
dp),
DIMENSION(:, :, :),
POINTER :: ovlp
101 TYPE(atom_basis_type),
POINTER :: basis, basis_grb, basis_ref, basis_vrb
102 TYPE(atom_integrals),
POINTER :: atint
103 TYPE(atom_orbitals),
POINTER :: orbitals
104 TYPE(atom_state),
POINTER :: state
105 TYPE(atom_type),
POINTER ::
atom, atom_ref, atom_test
107 TYPE(section_vals_type),
POINTER :: grb_section, powell_section
109 IF (iw > 0)
WRITE (iw,
'(/," ",79("*"),/,T28,A,/," ",79("*"))')
"GEOMETRICAL RESPONSE BASIS"
112 NULLIFY (vbasis(i)%basis)
116 IF (iw > 0 .AND. is > 1)
THEN
117 WRITE (iw,
'(/,A,/)')
" WARNING: Only use first electronic structure/method for basis set generation"
119 atom_ref => atom_info(1, 1)%atom
122 method = atom_ref%method_type
127 cpabort(
"Unrestricted methods not allowed for GRB generation")
138 CALL copy_atom_basics(atom_ref,
atom, state=.true., potential=.true., optimization=.true.,
xc=.true.)
140 atom%potential%confinement = .true.
142 atom%potential%acon = 200._dp
143 atom%potential%rcon = 4._dp
145 atom%potential%scon = scon
147 basis_ref => atom_ref%basis
149 NULLIFY (basis%am, basis%cm, basis%as, basis%ns, basis%bf, basis%dbf, basis%ddbf)
157 cpabort(
"# point radial grid < 0")
161 maxl =
atom%state%maxl_occ
165 basis%nbas(0:maxl) = num_gto
166 basis%nprim = basis%nbas
169 m = maxval(basis%nbas)
170 ALLOCATE (basis%am(m, 0:
lmat))
173 DO i = 1, basis%nbas(l)
175 basis%am(i, l) = aval*cval**(ll)
179 basis%eps_eig = basis_ref%eps_eig
180 basis%geometrical = .true.
187 m = maxval(basis%nbas)
188 ALLOCATE (basis%bf(nr, m, 0:
lmat))
189 ALLOCATE (basis%dbf(nr, m, 0:
lmat))
190 ALLOCATE (basis%ddbf(nr, m, 0:
lmat))
195 DO i = 1, basis%nbas(l)
198 rk = basis%grid%rad(k)
199 ear = exp(-al*basis%grid%rad(k)**2)
200 basis%bf(k, i, l) = rk**l*ear
201 basis%dbf(k, i, l) = (real(l,
dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear
202 basis%ddbf(k, i, l) = (real(l*(l - 1),
dp)*rk**(l - 2) - &
203 2._dp*al*real(2*l + 1,
dp)*rk**(l) + 4._dp*al*rk**(l + 2))*ear
209 mo = maxval(
atom%state%maxn_calc)
210 mb = maxval(basis%nbas)
215 CALL atom_fit_grb(
atom, basis, iw, powell_section)
222 WRITE (iw,
'(/,A,T76,I5)')
" Generate Response Basis Sets with Order ", nder
230 DO l = 0, state%maxl_occ
231 DO i = 1, state%maxn_occ(l)
232 IF (
atom%orbitals%ener(i, l) > emax)
THEN
235 emax =
atom%orbitals%ener(i, l)
236 fhomo = state%occupation(l, i)
241 s1 =
SIZE(
atom%orbitals%wfn, 1)
242 s2 =
SIZE(
atom%orbitals%wfn, 2)
243 ALLOCATE (wfn(s1, s2, 0:
lmat, -nder:nder))
244 s2 = maxval(state%maxn_occ) + nder
245 ALLOCATE (rbasis(s1, s2, 0:
lmat), qbasis(s1, s2, 0:
lmat))
251 CALL atom_int_setup(atint, basis, potential=
atom%potential, eri_coulomb=.false., eri_exchange=.false.)
253 IF (
atom%pp_calc)
THEN
254 NULLIFY (atint%tzora, atint%hdkh)
262 DO ider = -nder, nder
263 dene = real(ider, kind=
dp)*delta
264 cpassert(fhomo > abs(dene))
265 state%occupation(lhomo, nhomo) = fhomo + dene
267 wfn(:, :, :, ider) =
atom%orbitals%wfn
268 state%occupation(lhomo, nhomo) = fhomo
271 WRITE (iw,
'(A,T76,I5)')
" Total number of electronic structure calculations ", 2*nder + 1
274 ovlp =>
atom%integrals%ovlp
276 DO l = 0, state%maxl_occ
278 WRITE (iw,
'(A,T76,I5)')
" Response derivatives for l quantum number ", l
281 DO i = 1, max(state%maxn_occ(l), 1)
282 rbasis(:, i, l) = wfn(:, i, l, 0)
286 i = max(state%maxn_occ(l), 1)
289 rbasis(:, i + 1, l) = 0.5_dp*(wfn(:, i, l, 1) - wfn(:, i, l, -1))/delta
291 rbasis(:, i + 2, l) = 0.25_dp*(wfn(:, i, l, 2) - 2._dp*wfn(:, i, l, 0) + wfn(:, i, l, -2))/delta**2
293 rbasis(:, i + 3, l) = 0.125_dp*(wfn(:, i, l, 3) - 3._dp*wfn(:, i, l, 1) &
294 + 3._dp*wfn(:, i, l, -1) - wfn(:, i, l, -3))/delta**3
301 n = state%maxn_occ(l) + nder
302 m =
atom%basis%nbas(l)
305 o = dot_product(rbasis(1:m, j, l), reshape(matmul(ovlp(1:m, 1:m, l), rbasis(1:m, i:i, l)), (/m/)))
306 rbasis(1:m, i, l) = rbasis(1:m, i, l) - o*rbasis(1:m, j, l)
308 o = dot_product(rbasis(1:m, i, l), reshape(matmul(ovlp(1:m, 1:m, l), rbasis(1:m, i:i, l)), (/m/)))
309 rbasis(1:m, i, l) = rbasis(1:m, i, l)/sqrt(o)
313 ALLOCATE (amat(n, n))
314 amat(1:n, 1:n) = matmul(transpose(rbasis(1:m, 1:n, l)), matmul(ovlp(1:m, 1:m, l), rbasis(1:m, 1:n, l)))
316 amat(i, i) = amat(i, i) - 1._dp
318 IF (maxval(abs(amat)) > 1.e-12)
THEN
319 IF (iw > 0)
WRITE (iw,
'(A,G20.10)')
" Orthogonality error ", maxval(abs(amat))
324 expzet = 0.25_dp*real(2*l + 3,
dp)
325 prefac = sqrt(
rootpi/2._dp**(l + 2)*
dfac(2*l + 1))
327 zeta = (2._dp*
atom%basis%am(i, l))**expzet
328 qbasis(i, 1:n, l) = rbasis(i, 1:n, l)*prefac/zeta
334 IF (iw > 0)
WRITE (iw,
'(/,A)')
" Condition Number of Valence Response Basis Sets"
340 NULLIFY (basis_vrb%am, basis_vrb%cm, basis_vrb%as, basis_vrb%ns, basis_vrb%bf, &
341 basis_vrb%dbf, basis_vrb%ddbf)
343 NULLIFY (basis_vrb%grid)
346 basis_vrb%grid%nr = ngp
348 basis_vrb%eps_eig = basis_ref%eps_eig
349 basis_vrb%geometrical = .false.
351 basis_vrb%nprim = basis%nprim
353 DO l = 0, state%maxl_occ
354 basis_vrb%nbas(l) = state%maxn_occ(l) + ider
356 m = maxval(basis_vrb%nprim)
357 n = maxval(basis_vrb%nbas)
358 ALLOCATE (basis_vrb%am(m, 0:
lmat))
359 basis_vrb%am = basis%am
361 ALLOCATE (basis_vrb%cm(m, n, 0:
lmat))
362 DO l = 0, state%maxl_occ
363 m = basis_vrb%nprim(l)
364 n = basis_vrb%nbas(l)
365 basis_vrb%cm(1:m, 1:n, l) = rbasis(1:m, 1:n, l)
369 nr = basis_vrb%grid%nr
370 m = maxval(basis_vrb%nbas)
371 ALLOCATE (basis_vrb%bf(nr, m, 0:
lmat))
372 ALLOCATE (basis_vrb%dbf(nr, m, 0:
lmat))
373 ALLOCATE (basis_vrb%ddbf(nr, m, 0:
lmat))
375 basis_vrb%dbf = 0._dp
376 basis_vrb%ddbf = 0._dp
378 DO i = 1, basis_vrb%nprim(l)
379 al = basis_vrb%am(i, l)
381 rk = basis_vrb%grid%rad(k)
382 ear = exp(-al*basis_vrb%grid%rad(k)**2)
383 DO j = 1, basis_vrb%nbas(l)
384 basis_vrb%bf(k, j, l) = basis_vrb%bf(k, j, l) + rk**l*ear*basis_vrb%cm(i, j, l)
385 basis_vrb%dbf(k, j, l) = basis_vrb%dbf(k, j, l) &
386 + (real(l,
dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear*basis_vrb%cm(i, j, l)
387 basis_vrb%ddbf(k, j, l) = basis_vrb%ddbf(k, j, l) + &
388 (real(l*(l - 1),
dp)*rk**(l - 2) - 2._dp*al*real(2*l + 1,
dp)*rk**(l) + &
389 4._dp*al*rk**(l + 2))*ear*basis_vrb%cm(i, j, l)
396 CALL basis_label(abas, basis_vrb%nprim, basis_vrb%nbas)
397 WRITE (iw,
'(A,A)')
" Basis set ", trim(abas)
402 IF (iw > 0)
WRITE (iw,
'(T5,A,F15.3,T50,A,F14.4)')
" Lattice constant:", cradx,
"Condition number:", cnum
405 IF (iw > 0)
WRITE (iw,
'(T5,A,F15.3,T50,A,F14.4)')
" Lattice constant:", cradx,
"Condition number:", cnum
408 IF (iw > 0)
WRITE (iw,
'(T5,A,F15.3,T50,A,F14.4)')
" Lattice constant:", cradx,
"Condition number:", cnum
409 vbasis(ider)%basis => basis_vrb
415 ALLOCATE (rho(basis%grid%nr))
418 n = sum(maxloc(rho(:)))
419 rmax = basis%grid%rad(n)
420 IF (rmax < 0.1_dp) rmax = 1.0_dp
424 maxl =
atom%state%maxl_occ
427 IF (num_gto > 0)
THEN
429 WRITE (iw,
'(/,A)')
" Polarization basis set "
431 ALLOCATE (pbasis(num_gto, num_gto, 0:7), alp(num_gto))
435 zval = sqrt(real(2*lval + 2,
dp))*real(lval + 1,
dp)/(2._dp*rmax)
436 aval =
atom%basis%am(1, 0)
438 rconf =
atom%potential%scon
439 CALL atom_fit_pol(zval, rconf, lval, aval, cval, num_gto, iw, powell_section)
442 alp(i) = aval*cval**(i - 1)
444 ALLOCATE (rho(num_gto))
445 DO l = maxl + 1, min(maxl + num_gto, 7)
446 zval = sqrt(real(2*l + 2,
dp))*real(l + 1,
dp)/(2._dp*rmax)
447 CALL hydrogenic(zval, rconf, l, alp, num_gto, rho, pbasis(:, :, l))
448 IF (iw > 0)
WRITE (iw,
'(T5,A,i5,T66,A,F10.4)') &
449 " Polarization basis set contraction for lval=", l,
"zval=", zval
455 maxl =
atom%state%maxl_occ
459 IF (num_bas(1) == -1)
THEN
461 next_bas(l) = maxl - l + 1
464 n = min(
SIZE(num_bas, 1), 4)
465 next_bas(0:n - 1) = num_bas(1:n)
469 IF (next_bas(l) > 0) next_prim(l) = num_gto
472 CALL basis_label(abas, next_prim, next_bas)
473 WRITE (iw,
'(/,A,A)')
" Extension basis set ", trim(abas)
475 n = maxval(next_prim)
477 ALLOCATE (ebasis(n, n, 0:
lmat), ale(n))
478 basis_vrb => vbasis(0)%basis
479 amin =
atom%basis%aval/
atom%basis%cval**1.5_dp
481 ale(i) = amin*
atom%basis%cval**(i - 1)
485 rconf = 2.0_dp*
atom%potential%scon
487 IF (next_bas(l) < 1) cycle
488 zval = sqrt(real(2*l + 2,
dp))*real(l + 1,
dp)/(2._dp*rmax)
489 CALL hydrogenic(zval, rconf, l, ale, n, rho, ebasis(:, :, l))
490 IF (iw > 0)
WRITE (iw,
'(T5,A,i5,T66,A,F10.4)') &
491 " Extension basis set contraction for lval=", l,
"zval=", zval
495 IF (iw > 0)
WRITE (iw,
'(/,A)')
" Condition Number of Extended Basis Sets"
501 NULLIFY (basis_vrb%am, basis_vrb%cm, basis_vrb%as, basis_vrb%ns, basis_vrb%bf, &
502 basis_vrb%dbf, basis_vrb%ddbf)
504 NULLIFY (basis_vrb%grid)
507 basis_vrb%grid%nr = ngp
509 basis_vrb%eps_eig = basis_ref%eps_eig
510 basis_vrb%geometrical = .false.
512 basis_vrb%nprim = basis%nprim + next_prim
514 DO l = 0, state%maxl_occ
515 basis_vrb%nbas(l) = state%maxn_occ(l) + ider + next_bas(l)
517 m = maxval(basis_vrb%nprim)
518 ALLOCATE (basis_vrb%am(m, 0:
lmat))
520 m =
SIZE(basis%am, 1)
521 basis_vrb%am(1:m, :) = basis%am(1:m, :)
523 DO l = 0, state%maxl_occ
524 basis_vrb%am(m + 1:m + n, l) = ale(1:n)
527 m = maxval(basis_vrb%nprim)
528 n = maxval(basis_vrb%nbas)
529 ALLOCATE (basis_vrb%cm(m, n, 0:
lmat))
530 basis_vrb%cm = 0.0_dp
531 DO l = 0, state%maxl_occ
533 n = state%maxn_occ(l) + ider
534 basis_vrb%cm(1:m, 1:n, l) = rbasis(1:m, 1:n, l)
535 basis_vrb%cm(m + 1:m + next_prim(l), n + 1:n + next_bas(l), l) = ebasis(1:next_prim(l), 1:next_bas(l), l)
539 nr = basis_vrb%grid%nr
540 m = maxval(basis_vrb%nbas)
541 ALLOCATE (basis_vrb%bf(nr, m, 0:
lmat))
542 ALLOCATE (basis_vrb%dbf(nr, m, 0:
lmat))
543 ALLOCATE (basis_vrb%ddbf(nr, m, 0:
lmat))
545 basis_vrb%dbf = 0._dp
546 basis_vrb%ddbf = 0._dp
548 DO i = 1, basis_vrb%nprim(l)
549 al = basis_vrb%am(i, l)
551 rk = basis_vrb%grid%rad(k)
552 ear = exp(-al*basis_vrb%grid%rad(k)**2)
553 DO j = 1, basis_vrb%nbas(l)
554 basis_vrb%bf(k, j, l) = basis_vrb%bf(k, j, l) + rk**l*ear*basis_vrb%cm(i, j, l)
555 basis_vrb%dbf(k, j, l) = basis_vrb%dbf(k, j, l) &
556 + (real(l,
dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear*basis_vrb%cm(i, j, l)
557 basis_vrb%ddbf(k, j, l) = basis_vrb%ddbf(k, j, l) + &
558 (real(l*(l - 1),
dp)*rk**(l - 2) - 2._dp*al*real(2*l + 1,
dp)*rk**(l) + &
559 4._dp*al*rk**(l + 2))*ear*basis_vrb%cm(i, j, l)
566 CALL basis_label(abas, basis_vrb%nprim, basis_vrb%nbas)
567 WRITE (iw,
'(A,A)')
" Basis set ", trim(abas)
572 IF (iw > 0)
WRITE (iw,
'(T5,A,F15.3,T50,A,F14.4)')
" Lattice constant:", cradx,
"Condition number:", cnum
575 IF (iw > 0)
WRITE (iw,
'(T5,A,F15.3,T50,A,F14.4)')
" Lattice constant:", cradx,
"Condition number:", cnum
578 IF (iw > 0)
WRITE (iw,
'(T5,A,F15.3,T50,A,F14.4)')
" Lattice constant:", cradx,
"Condition number:", cnum
579 vbasis(nder + 1 + ider)%basis => basis_vrb
585 energy_ref = atom_ref%energy%etot
586 IF (iw > 0)
WRITE (iw,
'(/,A,A)')
" Basis set tests "
587 IF (iw > 0)
WRITE (iw,
'(T10,A,T59,F22.9)')
" Reference Energy [a.u.] ", energy_ref
588 DO ider = 0, 2*nder + 1
592 CALL copy_atom_basics(atom_ref, atom_test, state=.true., potential=.true., &
593 optimization=.true.,
xc=.true.)
594 basis_grb => vbasis(ider)%basis
596 mo = maxval(atom_test%state%maxn_calc)
597 mb = maxval(basis_grb%nbas)
599 CALL set_atom(atom_test, orbitals=orbitals, basis=basis_grb)
602 CALL atom_int_setup(atint, basis_grb, potential=atom_test%potential, eri_coulomb=.false., eri_exchange=.false.)
604 IF (atom_test%pp_calc)
THEN
605 NULLIFY (atint%tzora, atint%hdkh)
608 CALL atom_relint_setup(atint, basis_grb, atom_test%relativistic, zcore=real(atom_test%z,
dp))
610 CALL set_atom(atom_test, integrals=atint)
613 IF (ider <= nder)
THEN
614 energy_vb(ider) = atom_test%energy%etot
615 IF (iw > 0)
WRITE (iw,
'(T10,A,i1,A,T40,F13.9,T59,F22.9)')
" GRB (VB)", ider,
" Energy [a.u.] ", &
616 energy_ref - energy_vb(ider), energy_vb(ider)
619 energy_ex(i) = atom_test%energy%etot
620 IF (iw > 0)
WRITE (iw,
'(T10,A,i1,A,T40,F13.9,T59,F22.9)')
" GRB (EX)", i,
" Energy [a.u.] ", &
621 energy_ref - energy_ex(i), energy_ex(i)
626 DEALLOCATE (atom_test%state, atom_test%potential, atint)
632 expzet = 0.25_dp*real(2*l + 3,
dp)
633 prefac = sqrt(
rootpi/2._dp**(l + 2)*
dfac(2*l + 1))
635 zeta = (2._dp*alp(i))**expzet
636 pbasis(i, 1:num_pol, l) = pbasis(i, 1:num_pol, l)*prefac/zeta
641 expzet = 0.25_dp*real(2*l + 3,
dp)
642 prefac = sqrt(
rootpi/2._dp**(l + 2)*
dfac(2*l + 1))
643 DO i = 1, next_prim(l)
644 zeta = (2._dp*ale(i))**expzet
645 ebasis(i, 1:next_bas(l), l) = ebasis(i, 1:next_bas(l), l)*prefac/zeta
651 CALL open_file(file_name=
"GRB_BASIS", file_status=
"UNKNOWN", file_action=
"WRITE", unit_number=iunit)
655 headline(2) =
"# Generated with CP2K Atom Code"
657 CALL grb_print_basis(
header=headline, iunit=iunit)
660 WRITE (basline(1),
"(T2,A)") adjustl(
ptable(atom_ref%z)%symbol)
663 WRITE (basline(1),
"(T2,A,T5,A,I1)") adjustl(
ptable(atom_ref%z)%symbol), trim(adjustl(basname))//
"-VAL-", ider
664 CALL grb_print_basis(
header=basline, nprim=vbasis(ider)%basis%nprim(0), nbas=vbasis(ider)%basis%nbas, &
665 al=vbasis(ider)%basis%am(:, 0), gcc=qbasis, iunit=iunit)
668 maxl = atom_ref%state%maxl_occ
669 DO l = maxl + 1, min(maxl + num_pol, 7)
676 WRITE (basline(1),
"(T2,A,T5,A,I1)") adjustl(
ptable(atom_ref%z)%symbol), trim(adjustl(basname))//
"-POL-", i
677 CALL grb_print_basis(
header=basline, nprim=num_pol, nbas=nbas, al=alp, gcc=pbasis, iunit=iunit)
680 IF (sum(next_bas) > 0)
THEN
682 WRITE (basline(1),
"(T2,A,T5,A)") adjustl(
ptable(atom_ref%z)%symbol), trim(adjustl(basname))//
"-EXT"
683 CALL grb_print_basis(
header=basline, nprim=next_prim(0), nbas=next_bas, al=ale, gcc=ebasis, iunit=iunit)
689 IF (
ALLOCATED(pbasis))
DEALLOCATE (pbasis)
690 IF (
ALLOCATED(alp))
DEALLOCATE (alp)
691 IF (
ALLOCATED(ebasis))
DEALLOCATE (ebasis)
692 DEALLOCATE (wfn, rbasis, qbasis, ale)
695 IF (
ASSOCIATED(vbasis(ider)%basis))
THEN
697 DEALLOCATE (vbasis(ider)%basis)
705 DEALLOCATE (
atom%potential,
atom%state,
atom%integrals, basis)
708 IF (iw > 0)
WRITE (iw,
'(" ",79("*"))')
724 SUBROUTINE grb_print_basis(header, nprim, nbas, al, gcc, iunit)
725 CHARACTER(len=*),
DIMENSION(:),
INTENT(IN), &
727 INTEGER,
INTENT(IN),
OPTIONAL :: nprim
728 INTEGER,
DIMENSION(0:),
INTENT(IN),
OPTIONAL :: nbas
729 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN),
OPTIONAL :: al
730 REAL(kind=
dp),
DIMENSION(:, :, 0:),
INTENT(IN), &
732 INTEGER,
INTENT(IN) :: iunit
734 INTEGER :: i, j, l, lmax, lmin, nval
738 IF (
header(i) .NE.
"")
THEN
739 WRITE (iunit,
"(A)") trim(
header(i))
744 IF (
PRESENT(nprim))
THEN
746 cpassert(
PRESENT(nbas))
747 cpassert(
PRESENT(al))
748 cpassert(
PRESENT(gcc))
750 DO i = lbound(nbas, 1), ubound(nbas, 1)
751 IF (nbas(i) > 0)
THEN
756 DO i = ubound(nbas, 1), lbound(nbas, 1), -1
757 IF (nbas(i) > 0)
THEN
764 WRITE (iunit, *)
" 1"
765 WRITE (iunit,
"(40I3)") nval, lmin, lmax, nprim, (nbas(l), l=lmin, lmax)
767 WRITE (iunit,
"(G20.12)", advance=
"no") al(i)
770 WRITE (iunit,
"(F16.10)", advance=
"no") gcc(i, j, l)
779 END SUBROUTINE grb_print_basis
790 SUBROUTINE basis_label(label, np, nb)
791 CHARACTER(len=*),
INTENT(out) :: label
792 INTEGER,
DIMENSION(0:),
INTENT(in) :: np, nb
794 INTEGER :: i, l, lmax
795 CHARACTER(len=1),
DIMENSION(0:7),
PARAMETER :: &
796 lq = (/
"s",
"p",
"d",
"f",
"g",
"h",
"i",
"k"/)
799 lmax = min(ubound(np, 1), ubound(nb, 1), 7)
806 WRITE (label(i:i + 1),
"(I2)") np(l)
809 WRITE (label(i:i),
"(I1)") np(l)
816 label(i:i + 6) =
") --> ["
822 WRITE (label(i:i + 1),
"(I2)") nb(l)
825 WRITE (label(i:i),
"(I1)") nb(l)
834 END SUBROUTINE basis_label
845 SUBROUTINE grb_fit(atom, basis, afun, iw)
846 TYPE(atom_type),
POINTER ::
atom
847 TYPE(atom_basis_type),
POINTER :: basis
848 REAL(
dp),
INTENT(OUT) :: afun
849 INTEGER,
INTENT(IN) :: iw
851 INTEGER :: do_eric, do_erie, reltyp, zval
852 LOGICAL :: eri_c, eri_e
853 TYPE(atom_integrals),
POINTER :: atint
854 TYPE(atom_potential_type),
POINTER :: pot
861 pot =>
atom%potential
863 reltyp =
atom%relativistic
864 do_eric =
atom%coulomb_integral_type
865 do_erie =
atom%exchange_integral_type
869 CALL atom_int_setup(atint, basis, potential=pot, eri_coulomb=eri_c, eri_exchange=eri_e)
872 IF (
atom%pp_calc)
THEN
873 NULLIFY (atint%tzora, atint%hdkh)
881 afun =
atom%energy%etot
886 END SUBROUTINE grb_fit
899 SUBROUTINE copy_atom_basics(atom_ref, atom, state, potential, optimization, xc)
900 TYPE(atom_type),
POINTER :: atom_ref,
atom
901 LOGICAL,
INTENT(IN),
OPTIONAL :: state, potential, optimization,
xc
904 atom%zcore = atom_ref%zcore
905 atom%pp_calc = atom_ref%pp_calc
906 atom%method_type = atom_ref%method_type
907 atom%relativistic = atom_ref%relativistic
908 atom%coulomb_integral_type = atom_ref%coulomb_integral_type
909 atom%exchange_integral_type = atom_ref%exchange_integral_type
914 IF (
PRESENT(state))
THEN
916 ALLOCATE (
atom%state)
917 atom%state = atom_ref%state
921 IF (
PRESENT(potential))
THEN
923 ALLOCATE (
atom%potential)
924 atom%potential = atom_ref%potential
928 IF (
PRESENT(optimization))
THEN
929 IF (optimization)
THEN
930 atom%optimization = atom_ref%optimization
934 IF (
PRESENT(
xc))
THEN
936 atom%xc_section => atom_ref%xc_section
940 END SUBROUTINE copy_atom_basics
951 SUBROUTINE atom_fit_grb(atom, basis, iunit, powell_section)
952 TYPE(atom_type),
POINTER ::
atom
953 TYPE(atom_basis_type),
POINTER :: basis
954 INTEGER,
INTENT(IN) :: iunit
955 TYPE(section_vals_type),
POINTER :: powell_section
957 INTEGER :: i, k, l, ll, n10, nr
958 REAL(kind=
dp) :: al, cnum, crad, cradx, ear, fopt, rk
959 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: x
962 cpassert(basis%geometrical)
970 x(1) = sqrt(basis%aval)
971 x(2) = sqrt(basis%cval)
979 WRITE (iunit,
'(/," POWELL| Start optimization procedure")')
980 WRITE (iunit,
'(" POWELL| Total number of parameters in optimization",T71,I10)') ostate%nvar
982 n10 = max(ostate%maxfun/100, 1)
988 IF (ostate%state == 2)
THEN
991 DO i = 1, basis%nbas(l)
992 ll = i - 1 + basis%start(l)
993 basis%am(i, l) = x(1)*x(1)*(x(2)*x(2))**(ll)
996 basis%aval = x(1)*x(1)
997 basis%cval = x(2)*x(2)
1003 DO i = 1, basis%nbas(l)
1006 rk = basis%grid%rad(k)
1007 ear = exp(-al*basis%grid%rad(k)**2)
1008 basis%bf(k, i, l) = rk**l*ear
1009 basis%dbf(k, i, l) = (real(l,
dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear
1010 basis%ddbf(k, i, l) = (real(l*(l - 1),
dp)*rk**(l - 2) - &
1011 2._dp*al*real(2*l + 1,
dp)*rk**(l) + 4._dp*al*rk**(l + 2))*ear
1015 CALL grb_fit(
atom, basis, ostate%f, 0)
1016 fopt = min(fopt, ostate%f)
1019 IF (ostate%state == -1)
EXIT
1023 IF (ostate%nf == 2 .AND. iunit > 0)
THEN
1024 WRITE (iunit,
'(" POWELL| Initial value of function",T61,F20.10)') ostate%f
1026 IF (mod(ostate%nf, n10) == 0 .AND. iunit > 0)
THEN
1027 WRITE (iunit,
'(" POWELL| Reached",i4,"% of maximal function calls",T61,F20.10)') &
1028 int(real(ostate%nf,
dp)/real(ostate%maxfun,
dp)*100._dp), fopt
1037 WRITE (iunit,
'(" POWELL| Number of function evaluations",T71,I10)') ostate%nf
1038 WRITE (iunit,
'(" POWELL| Final value of function",T61,F20.10)') ostate%fopt
1043 DO i = 1, basis%nbas(l)
1044 ll = i - 1 + basis%start(l)
1045 basis%am(i, l) = x(1)*x(1)*(x(2)*x(2))**(ll)
1048 basis%aval = x(1)*x(1)
1049 basis%cval = x(2)*x(2)
1055 DO i = 1, basis%nbas(l)
1058 rk = basis%grid%rad(k)
1059 ear = exp(-al*basis%grid%rad(k)**2)
1060 basis%bf(k, i, l) = rk**l*ear
1061 basis%dbf(k, i, l) = (real(l,
dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear
1062 basis%ddbf(k, i, l) = (real(l*(l - 1),
dp)*rk**(l - 2) - &
1063 2._dp*al*real(2*l + 1,
dp)*rk**(l) + 4._dp*al*rk**(l + 2))*ear
1072 WRITE (iunit,
'(/,A)')
" Optimized Geometrical GTO basis set"
1073 WRITE (iunit,
'(A,F15.8,T41,A,F15.8)')
" Initial exponent: ", basis%aval, &
1074 " Proportionality factor: ", basis%cval
1076 WRITE (iunit,
'(T41,A,I2,T76,I5)')
" Number of exponents for l=", l, basis%nbas(l)
1080 IF (iunit > 0)
WRITE (iunit,
'(/,A)')
" Condition number of uncontracted basis set"
1084 cradx = crad*1.00_dp
1086 IF (iunit > 0)
WRITE (iunit,
'(T5,A,F15.3,T50,A,F14.4)')
" Lattice constant:", cradx,
"Condition number:", cnum
1087 cradx = crad*1.10_dp
1089 IF (iunit > 0)
WRITE (iunit,
'(T5,A,F15.3,T50,A,F14.4)')
" Lattice constant:", cradx,
"Condition number:", cnum
1090 cradx = crad*1.20_dp
1092 IF (iunit > 0)
WRITE (iunit,
'(T5,A,F15.3,T50,A,F14.4)')
" Lattice constant:", cradx,
"Condition number:", cnum
1096 END SUBROUTINE atom_fit_grb
1111 SUBROUTINE atom_fit_pol(zval, rconf, lval, aval, cval, nbas, iunit, powell_section)
1112 REAL(kind=
dp),
INTENT(IN) :: zval, rconf
1113 INTEGER,
INTENT(IN) :: lval
1114 REAL(kind=
dp),
INTENT(INOUT) :: aval, cval
1115 INTEGER,
INTENT(IN) :: nbas, iunit
1116 TYPE(section_vals_type),
POINTER :: powell_section
1119 REAL(kind=
dp) :: fopt, x(2)
1120 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: am, ener
1121 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: orb
1124 ALLOCATE (am(nbas), ener(nbas), orb(nbas, nbas))
1140 WRITE (iunit,
'(/," POWELL| Start optimization procedure")')
1141 WRITE (iunit,
'(" POWELL| Total number of parameters in optimization",T71,I10)') ostate%nvar
1143 n10 = max(ostate%maxfun/100, 1)
1149 IF (ostate%state == 2)
THEN
1153 am(i) = aval*cval**(i - 1)
1155 CALL hydrogenic(zval, rconf, lval, am, nbas, ener, orb)
1157 fopt = min(fopt, ostate%f)
1160 IF (ostate%state == -1)
EXIT
1164 IF (ostate%nf == 2 .AND. iunit > 0)
THEN
1165 WRITE (iunit,
'(" POWELL| Initial value of function",T61,F20.10)') ostate%f
1167 IF (mod(ostate%nf, n10) == 0 .AND. iunit > 0)
THEN
1168 WRITE (iunit,
'(" POWELL| Reached",i4,"% of maximal function calls",T61,F20.10)') &
1169 int(real(ostate%nf,
dp)/real(ostate%maxfun,
dp)*100._dp), fopt
1178 WRITE (iunit,
'(" POWELL| Number of function evaluations",T71,I10)') ostate%nf
1179 WRITE (iunit,
'(" POWELL| Final value of function",T61,F20.10)') ostate%fopt
1187 WRITE (iunit,
'(/,A,T51,A,T76,I5)')
" Optimized Polarization basis set", &
1188 " Number of exponents:", nbas
1189 WRITE (iunit,
'(A,F15.8,T41,A,F15.8)')
" Initial exponent: ", aval, &
1190 " Proportionality factor: ", cval
1193 DEALLOCATE (am, ener, orb)
1195 END SUBROUTINE atom_fit_pol
1209 SUBROUTINE hydrogenic(zval, rconf, lval, am, nbas, ener, orb)
1210 REAL(kind=
dp),
INTENT(IN) :: zval, rconf
1211 INTEGER,
INTENT(IN) :: lval
1212 REAL(kind=
dp),
DIMENSION(:),
INTENT(IN) :: am
1213 INTEGER,
INTENT(IN) :: nbas
1214 REAL(kind=
dp),
DIMENSION(:),
INTENT(OUT) :: ener
1215 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(OUT) :: orb
1217 INTEGER :: info, k, lwork, n
1219 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: w, work
1220 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: confmat, hmat, potmat, smat, tmat
1223 ALLOCATE (smat(n, n), tmat(n, n), potmat(n, n), confmat(n, n), hmat(n, n))
1225 CALL sg_overlap(smat(1:n, 1:n), lval, am(1:n), am(1:n))
1227 CALL sg_kinetic(tmat(1:n, 1:n), lval, am(1:n), am(1:n))
1229 CALL sg_nuclear(potmat(1:n, 1:n), lval, am(1:n), am(1:n))
1233 CALL sg_conf(confmat, rconf, k, lval, am(1:n), am(1:n))
1235 hmat(1:n, 1:n) = tmat(1:n, 1:n) - zval*potmat(1:n, 1:n) + cf*confmat(1:n, 1:n)
1238 ALLOCATE (w(n), work(lwork))
1239 CALL lapack_ssygv(1,
"V",
"U", n, hmat, n, smat, n, w, work, lwork, info)
1241 orb(1:n, 1:n) = hmat(1:n, 1:n)
1243 DEALLOCATE (w, work)
1244 DEALLOCATE (smat, tmat, potmat, confmat, hmat)
1246 END SUBROUTINE hydrogenic
subroutine, public sg_nuclear(umat, l, pa, pb)
...
subroutine, public sg_kinetic(kmat, l, pa, pb)
...
subroutine, public sg_conf(gmat, rc, k, l, pa, pb)
...
subroutine, public sg_overlap(smat, l, pa, pb)
...
subroutine, public calculate_atom(atom, iw, noguess, converged)
General routine to perform electronic structure atomic calculations.
subroutine, public atom_grb_construction(atom_info, atom_section, iw)
Construct geometrical response basis set.
Calculate the atomic operator matrices.
subroutine, public atom_ppint_release(integrals)
Release memory allocated for atomic integrals (core electrons).
subroutine, public atom_int_setup(integrals, basis, potential, eri_coulomb, eri_exchange, all_nu)
Set up atomic integrals.
subroutine, public atom_relint_setup(integrals, basis, reltyp, zcore, alpha)
...
subroutine, public atom_relint_release(integrals)
Release memory allocated for atomic integrals (relativistic effects).
subroutine, public atom_ppint_setup(integrals, basis, potential)
...
subroutine, public atom_int_release(integrals)
Release memory allocated for atomic integrals (valence electrons).
Define the atom type and its sub types.
subroutine, public create_atom_type(atom)
...
integer, parameter, public cgto_basis
integer, parameter, public gto_basis
subroutine, public release_atom_type(atom)
...
integer, parameter, public lmat
subroutine, public set_atom(atom, basis, state, integrals, orbitals, potential, zcore, pp_calc, do_zmp, doread, read_vxc, method_type, relativistic, coulomb_integral_type, exchange_integral_type, fmat)
...
subroutine, public release_atom_basis(basis)
...
subroutine, public create_atom_orbs(orbs, mbas, mo)
...
Some basic routines for atomic calculations.
subroutine, public atom_density(density, pmat, basis, maxl, typ, rr)
Map the electron density on an atomic radial grid.
subroutine, public atom_basis_condnum(basis, rad, cnum)
Calculate the condition number of the given atomic basis set.
Utility routines to open and close files. Tracking of preconnections.
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.
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.
Defines the basic variable types.
integer, parameter, public dp
integer, parameter, public default_string_length
Interface to the LAPACK F77 library.
Definition of mathematical constants and functions.
real(kind=dp), dimension(-1:2 *maxfac+1), parameter, public dfac
real(kind=dp), parameter, public rootpi
Provides Cartesian and spherical orbital pointers and indices.
subroutine, public init_orbital_pointers(maxl)
Initialize or update the orbital pointers.
subroutine, public deallocate_orbital_pointers()
Deallocate the orbital pointers.
Periodic Table related data definitions.
type(atom), dimension(0:nelem), public ptable
Definition of physical constants:
real(kind=dp), parameter, public bohr
subroutine, public powell_optimize(n, x, optstate)
...
subroutine, public allocate_grid_atom(grid_atom)
Initialize components of the grid_atom_type structure.
subroutine, public create_grid_atom(grid_atom, nr, na, llmax, ll, quadrature)
...
Exchange and Correlation functional calculations.