79 TYPE(
atom_p_type),
DIMENSION(:, :),
POINTER :: atom_info
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
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("*"))')