23 gto_basis,
sto_basis, atom_basis_type, atom_gthpot_type, atom_integrals, atom_p_type, &
50 #include "./base/base_uses.f90"
56 CHARACTER(len=*),
PARAMETER,
PRIVATE :: moduleN =
'atom_fit'
61 REAL(KIND=
dp),
DIMENSION(:, :, :),
POINTER :: wfn => null()
77 TYPE(atom_type),
POINTER ::
atom
78 INTEGER,
INTENT(IN) :: num_gto, norder, iunit
79 TYPE(section_vals_type),
OPTIONAL,
POINTER :: powell_section
80 REAL(kind=
dp),
DIMENSION(:),
OPTIONAL :: results
83 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: co
84 REAL(kind=
dp),
DIMENSION(2) :: x
85 TYPE(opgrid_type),
POINTER :: density
88 ALLOCATE (co(num_gto))
93 atom%state%maxl_occ,
atom%state%maxn_occ)
96 density%op =
fourpi*density%op
98 density%op = density%op*
atom%basis%grid%rad**norder
105 IF (
PRESENT(powell_section))
THEN
110 ostate%rhoend = 1.e-8_dp
111 ostate%rhobeg = 5.e-2_dp
119 WRITE (iunit,
'(/," POWELL| Start optimization procedure")')
121 n10 = ostate%maxfun/10
125 IF (ostate%state == 2)
THEN
126 CALL density_fit(density,
atom, num_gto, x(1), x(2), co, ostate%f)
129 IF (ostate%state == -1)
EXIT
133 IF (mod(ostate%nf, n10) == 0 .AND. iunit > 0)
THEN
134 WRITE (iunit,
'(" POWELL| Reached",i4,"% of maximal function calls")') &
135 int(real(ostate%nf,
dp)/real(ostate%maxfun,
dp)*100._dp)
142 CALL density_fit(density,
atom, num_gto, x(1), x(2), co, ostate%f)
147 WRITE (iunit,
'(" POWELL| Number of function evaluations",T71,I10)') ostate%nf
148 WRITE (iunit,
'(" POWELL| Final value of function",T61,G20.10)') ostate%fopt
149 WRITE (iunit,
'(" Optimized representation of density using a Geometrical Gaussian basis")')
150 WRITE (iunit,
'(A,I3,/,T10,A,F10.6,T48,A,F10.6)')
" Number of Gaussians:", num_gto, &
151 "Starting exponent:", x(1),
"Proportionality factor:", x(2)
152 WRITE (iunit,
'(A)')
" Expansion coefficients"
153 WRITE (iunit,
'(4F20.10)') co(1:num_gto)
156 IF (
PRESENT(results))
THEN
157 cpassert(
SIZE(results) >= num_gto + 2)
160 results(3:2 + num_gto) = co(1:num_gto)
175 TYPE(atom_p_type),
DIMENSION(:, :),
POINTER :: atom_info
176 TYPE(atom_basis_type),
POINTER :: basis
177 LOGICAL,
INTENT(IN) :: pptype
178 INTEGER,
INTENT(IN) :: iunit
179 TYPE(section_vals_type),
POINTER :: powell_section
181 INTEGER :: i, j, k, l, ll, m, n, n10, nl, nr, zval
182 INTEGER,
ALLOCATABLE,
DIMENSION(:, :) :: xtob
183 LOGICAL :: explicit, mult, penalty
184 REAL(kind=
dp) :: al, crad, ear, fopt, pf, rk
185 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: x
186 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: wem
187 REAL(kind=
dp),
DIMENSION(:),
POINTER :: w
191 SELECT CASE (basis%basis_type)
195 IF (basis%geometrical)
THEN
198 x(1) = sqrt(basis%aval)
199 x(2) = sqrt(basis%cval)
201 ll = maxval(basis%nprim(:))
202 ALLOCATE (xtob(ll, 0:
lmat))
204 ll = sum(basis%nprim(:))
209 DO i = 1, basis%nprim(l)
212 IF (abs(basis%am(i, l) - x(k)) < 1.e-6_dp)
THEN
219 x(ll) = basis%am(i, l)
225 DO i = 1, ostate%nvar
226 x(i) = sqrt(log(1._dp + x(i)))
231 ll = maxval(basis%nbas(:))
232 ALLOCATE (xtob(ll, 0:
lmat))
234 ll = sum(basis%nbas(:))
239 DO i = 1, basis%nbas(l)
241 x(ll) = basis%as(i, l)
246 DO i = 1, ostate%nvar
247 x(i) = sqrt(log(1._dp + x(i)))
255 n =
SIZE(atom_info, 1)
256 m =
SIZE(atom_info, 2)
262 DO i = 1, min(
SIZE(w), n)
263 wem(i, :) = w(i)*wem(i, :)
269 DO i = 1, min(
SIZE(w), m)
270 wem(:, i) = w(i)*wem(:, i)
274 DO i = 1,
SIZE(atom_info, 1)
275 DO j = 1,
SIZE(atom_info, 2)
276 atom_info(i, j)%atom%weight = wem(i, j)
287 WRITE (iunit,
'(/," POWELL| Start optimization procedure")')
288 WRITE (iunit,
'(/," POWELL| Total number of parameters in optimization",T71,I10)') ostate%nvar
290 n10 = max(ostate%maxfun/100, 1)
296 IF (ostate%state == 2)
THEN
297 SELECT CASE (basis%basis_type)
301 IF (basis%geometrical)
THEN
304 DO i = 1, basis%nbas(l)
305 ll = i - 1 + basis%start(l)
306 basis%am(i, l) = x(1)*x(1)*(x(2)*x(2))**(ll)
309 basis%aval = x(1)*x(1)
310 basis%cval = x(2)*x(2)
313 DO i = 1, basis%nprim(l)
314 al = x(xtob(i, l))**2
315 basis%am(i, l) = exp(al) - 1._dp
324 DO i = 1, basis%nbas(l)
327 rk = basis%grid%rad(k)
328 ear = exp(-al*basis%grid%rad(k)**2)
329 basis%bf(k, i, l) = rk**l*ear
330 basis%dbf(k, i, l) = (real(l,
dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear
331 basis%ddbf(k, i, l) = (real(l*(l - 1),
dp)*rk**(l - 2) - &
332 2._dp*al*real(2*l + 1,
dp)*rk**(l) + 4._dp*al*rk**(l + 2))*ear
338 DO i = 1, basis%nbas(l)
339 al = x(xtob(i, l))**2
340 basis%as(i, l) = exp(al) - 1._dp
348 DO i = 1, basis%nbas(l)
351 pf = (2._dp*al)**nl*sqrt(2._dp*al/
fac(2*nl))
353 rk = basis%grid%rad(k)
354 ear = rk**(nl - 1)*exp(-al*rk)
355 basis%bf(k, i, l) = pf*ear
356 basis%dbf(k, i, l) = pf*(real(nl - 1,
dp)/rk - al)*ear
357 basis%ddbf(k, i, l) = pf*(real((nl - 2)*(nl - 1),
dp)/rk/rk &
358 - al*real(2*(nl - 1),
dp)/rk + al*al)*ear
363 CALL basis_fit(atom_info, basis, pptype, ostate%f, 0, penalty)
364 fopt = min(fopt, ostate%f)
367 IF (ostate%state == -1)
EXIT
371 IF (ostate%nf == 2 .AND. iunit > 0)
THEN
372 WRITE (iunit,
'(" POWELL| Initial value of function",T61,F20.10)') ostate%f
374 IF (mod(ostate%nf, n10) == 0 .AND. iunit > 0)
THEN
375 WRITE (iunit,
'(" POWELL| Reached",i4,"% of maximal function calls",T61,F20.10)') &
376 int(real(ostate%nf,
dp)/real(ostate%maxfun,
dp)*100._dp), fopt
385 WRITE (iunit,
'(" POWELL| Number of function evaluations",T71,I10)') ostate%nf
386 WRITE (iunit,
'(" POWELL| Final value of function",T61,F20.10)') ostate%fopt
388 SELECT CASE (basis%basis_type)
392 IF (basis%geometrical)
THEN
395 DO i = 1, basis%nbas(l)
396 ll = i - 1 + basis%start(l)
397 basis%am(i, l) = x(1)*x(1)*(x(2)*x(2))**(ll)
400 basis%aval = x(1)*x(1)
401 basis%cval = x(2)*x(2)
404 DO i = 1, basis%nprim(l)
405 al = x(xtob(i, l))**2
406 basis%am(i, l) = exp(al) - 1._dp
415 DO i = 1, basis%nbas(l)
418 rk = basis%grid%rad(k)
419 ear = exp(-al*basis%grid%rad(k)**2)
420 basis%bf(k, i, l) = rk**l*ear
421 basis%dbf(k, i, l) = (real(l,
dp)*rk**(l - 1) - 2._dp*al*rk**(l + 1))*ear
422 basis%ddbf(k, i, l) = (real(l*(l - 1),
dp)*rk**(l - 2) - &
423 2._dp*al*real(2*l + 1,
dp)*rk**(l) + 4._dp*al*rk**(l + 2))*ear
429 DO i = 1, basis%nprim(l)
430 al = x(xtob(i, l))**2
431 basis%as(i, l) = exp(al) - 1._dp
439 DO i = 1, basis%nbas(l)
442 pf = (2._dp*al)**nl*sqrt(2._dp*al/
fac(2*nl))
444 rk = basis%grid%rad(k)
445 ear = rk**(nl - 1)*exp(-al*rk)
446 basis%bf(k, i, l) = pf*ear
447 basis%dbf(k, i, l) = pf*(real(nl - 1,
dp)/rk - al)*ear
448 basis%ddbf(k, i, l) = pf*(real((nl - 2)*(nl - 1),
dp)/rk/rk &
449 - al*real(2*(nl - 1),
dp)/rk + al*al)*ear
460 IF (
ALLOCATED(xtob))
THEN
464 SELECT CASE (basis%basis_type)
468 zval = atom_info(1, 1)%atom%z
488 TYPE(atom_p_type),
DIMENSION(:, :),
POINTER :: atom_info, atom_refs
489 TYPE(atom_potential_type),
POINTER :: ppot
490 INTEGER,
INTENT(IN) :: iunit
491 TYPE(section_vals_type),
POINTER :: powell_section
493 LOGICAL,
PARAMETER :: debug = .false.
495 CHARACTER(len=2) :: pc1, pc2, pct
496 INTEGER :: i, i1, i2, iw, j, j1, j2, k, l, m, n, &
497 n10, np, nre, nreinit, ntarget
498 INTEGER,
ALLOCATABLE,
DIMENSION(:, :) :: xtob
499 INTEGER,
DIMENSION(0:lmat) :: ncore
500 LOGICAL :: explicit, noopt_nlcc, preopt_nlcc
501 REAL(kind=
dp) :: afun, charge, de, deig, drho, dx, eig, fopt, oc, pchg, peig, pv, rcm, rcov, &
502 rmax, semicore_level, step_size_scaling, t_ener, t_psir0, t_semi, t_valence, t_virt, &
503 w_ener, w_node, w_psir0, w_semi, w_valence, w_virt, wtot
504 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: x, xi
505 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: wem
506 REAL(kind=
dp),
ALLOCATABLE, &
507 DIMENSION(:, :, :, :, :) :: dener, pval
508 REAL(kind=
dp),
DIMENSION(:),
POINTER :: w
509 TYPE(atom_type),
POINTER ::
atom
511 TYPE(wfn_init),
DIMENSION(:, :),
POINTER :: wfn_guess
540 n =
SIZE(atom_info, 1)
541 m =
SIZE(atom_info, 2)
544 ALLOCATE (pval(8, 10, 0:
lmat, m, n))
545 ALLOCATE (dener(2, m, m, n, n))
548 ALLOCATE (wfn_guess(n, m))
551 atom => atom_info(i, j)%atom
552 NULLIFY (wfn_guess(i, j)%wfn)
554 i1 =
SIZE(
atom%orbitals%wfn, 1)
555 i2 =
SIZE(
atom%orbitals%wfn, 2)
556 ALLOCATE (wfn_guess(i, j)%wfn(i1, i2, 0:
lmat))
564 DO i = 1, min(
SIZE(w), n)
565 wem(i, :) = w(i)*wem(i, :)
571 DO i = 1, min(
SIZE(w), m)
572 wem(:, i) = w(i)*wem(:, i)
577 CALL open_file(file_name=
"POWELL_RESULT", file_status=
"UNKNOWN", file_action=
"WRITE", unit_number=iw)
580 IF (ppot%gth_pot%nlcc)
THEN
581 CALL opt_nlcc_param(atom_info, atom_refs, ppot%gth_pot, iunit, preopt_nlcc)
584 preopt_nlcc = .false.
589 CALL get_pseudo_param(xi, ostate%nvar, ppot%gth_pot, noopt_nlcc)
590 ALLOCATE (x(ostate%nvar))
591 x(1:ostate%nvar) = xi(1:ostate%nvar)
600 WRITE (iunit,
'(/," POWELL| Start optimization procedure")')
601 WRITE (iunit,
'(/," POWELL| Total number of parameters in optimization",T71,I10)') ostate%nvar
603 n10 = max(ostate%maxfun/100, 1)
605 rcov =
ptable(atom_refs(1, 1)%atom%z)%covalent_radius*
bohr*rcm
609 DO i = 1,
SIZE(atom_info, 1)
610 DO j = 1,
SIZE(atom_info, 2)
611 atom => atom_info(i, j)%atom
613 dener(2, j, j, i, i) = atom_refs(i, j)%atom%energy%etot
614 IF (
atom%state%multiplicity == -1)
THEN
616 atom%weight = wem(i, j)
619 np =
atom%state%maxn_calc(l)
621 CALL atom_orbital_max(rmax, atom_refs(i, j)%atom%orbitals%wfn(:, ncore(l) + k, l), &
622 rcov, l, atom_refs(i, j)%atom%basis)
623 atom%orbitals%rcmax(k, l, 1) = max(rcov, rmax)
625 atom%orbitals%rcmax(k, l, 1), l, atom_refs(i, j)%atom%basis)
626 atom%orbitals%refene(k, l, 1) = atom_refs(i, j)%atom%orbitals%ener(ncore(l) + k, l)
627 atom%orbitals%refchg(k, l, 1) = charge
628 IF (k >
atom%state%maxn_occ(l))
THEN
629 IF (k <=
atom%state%maxn_occ(l) + 1)
THEN
630 atom%orbitals%wrefene(k, l, 1) = w_virt
631 atom%orbitals%wrefchg(k, l, 1) = w_virt/100._dp
632 atom%orbitals%crefene(k, l, 1) = t_virt
633 atom%orbitals%reftype(k, l, 1) =
"U1"
634 ntarget = ntarget + 2
635 wtot = wtot +
atom%weight*(w_virt + w_virt/100._dp)
637 atom%orbitals%wrefene(k, l, 1) = w_virt/100._dp
638 atom%orbitals%wrefchg(k, l, 1) = 0._dp
639 atom%orbitals%crefene(k, l, 1) = t_virt*10._dp
640 atom%orbitals%reftype(k, l, 1) =
"U2"
641 ntarget = ntarget + 1
642 wtot = wtot +
atom%weight*w_virt/100._dp
644 ELSEIF (k <
atom%state%maxn_occ(l))
THEN
645 atom%orbitals%wrefene(k, l, 1) = w_semi
646 atom%orbitals%wrefchg(k, l, 1) = w_semi/100._dp
647 atom%orbitals%crefene(k, l, 1) = t_semi
648 atom%orbitals%reftype(k, l, 1) =
"SC"
649 ntarget = ntarget + 2
650 wtot = wtot +
atom%weight*(w_semi + w_semi/100._dp)
652 IF (abs(
atom%state%occupation(l, k) - real(4*l + 2, kind=
dp)) < 0.01_dp .AND. &
653 abs(
atom%orbitals%refene(k, l, 1)) > semicore_level)
THEN
655 atom%orbitals%wrefene(k, l, 1) = w_semi
656 atom%orbitals%wrefchg(k, l, 1) = w_semi/100._dp
657 atom%orbitals%crefene(k, l, 1) = t_semi
658 atom%orbitals%reftype(k, l, 1) =
"SC"
659 wtot = wtot +
atom%weight*(w_semi + w_semi/100._dp)
661 atom%orbitals%wrefene(k, l, 1) = w_valence
662 atom%orbitals%wrefchg(k, l, 1) = w_valence/100._dp
663 atom%orbitals%crefene(k, l, 1) = t_valence
664 atom%orbitals%reftype(k, l, 1) =
"VA"
665 wtot = wtot +
atom%weight*(w_valence + w_valence/100._dp)
668 atom%orbitals%tpsir0(k, 1) = t_psir0
669 atom%orbitals%wpsir0(k, 1) = w_psir0
670 wtot = wtot +
atom%weight*w_psir0
672 ntarget = ntarget + 2
676 atom%orbitals%refnod(k, l, 1) = real(k - 1, kind=
dp)
678 IF (k == 1 .AND.
atom%state%occupation(l, k) /= 0._dp)
THEN
679 atom%orbitals%wrefnod(k, l, 1) = w_node
680 wtot = wtot +
atom%weight*w_node
686 atom%weight = wem(i, j)
689 np =
atom%state%maxn_calc(l)
691 CALL atom_orbital_max(rmax, atom_refs(i, j)%atom%orbitals%wfna(:, ncore(l) + k, l), &
692 rcov, l, atom_refs(i, j)%atom%basis)
693 atom%orbitals%rcmax(k, l, 1) = max(rcov, rmax)
694 CALL atom_orbital_max(rmax, atom_refs(i, j)%atom%orbitals%wfnb(:, ncore(l) + k, l), &
695 rcov, l, atom_refs(i, j)%atom%basis)
696 atom%orbitals%rcmax(k, l, 2) = max(rcov, rmax)
698 atom%orbitals%rcmax(k, l, 1), l, atom_refs(i, j)%atom%basis)
699 atom%orbitals%refene(k, l, 1) = atom_refs(i, j)%atom%orbitals%enera(ncore(l) + k, l)
700 atom%orbitals%refchg(k, l, 1) = charge
702 atom%orbitals%rcmax(k, l, 1), l, atom_refs(i, j)%atom%basis)
703 atom%orbitals%refene(k, l, 2) = atom_refs(i, j)%atom%orbitals%enerb(ncore(l) + k, l)
704 atom%orbitals%refchg(k, l, 2) = charge
706 IF (k >
atom%state%maxn_occ(l))
THEN
707 IF (k <=
atom%state%maxn_occ(l) + 1)
THEN
708 atom%orbitals%wrefene(k, l, 1:2) = w_virt
709 atom%orbitals%wrefchg(k, l, 1:2) = w_virt/100._dp
710 atom%orbitals%crefene(k, l, 1:2) = t_virt
711 atom%orbitals%reftype(k, l, 1:2) =
"U1"
712 ntarget = ntarget + 4
713 wtot = wtot +
atom%weight*2._dp*(w_virt + w_virt/100._dp)
715 atom%orbitals%wrefene(k, l, 1:2) = w_virt/100._dp
716 atom%orbitals%wrefchg(k, l, 1:2) = 0._dp
717 atom%orbitals%crefene(k, l, 1:2) = t_virt*10.0_dp
718 atom%orbitals%reftype(k, l, 1:2) =
"U2"
719 wtot = wtot +
atom%weight*2._dp*w_virt/100._dp
720 ntarget = ntarget + 2
722 ELSEIF (k <
atom%state%maxn_occ(l))
THEN
723 atom%orbitals%wrefene(k, l, 1:2) = w_semi
724 atom%orbitals%wrefchg(k, l, 1:2) = w_semi/100._dp
725 atom%orbitals%crefene(k, l, 1:2) = t_semi
726 atom%orbitals%reftype(k, l, 1:2) =
"SC"
727 ntarget = ntarget + 4
728 wtot = wtot +
atom%weight*2._dp*(w_semi + w_semi/100._dp)
730 IF (abs(
atom%state%occupation(l, k) - real(2*l + 1, kind=
dp)) < 0.01_dp .AND. &
731 abs(
atom%orbitals%refene(k, l, 1)) > semicore_level)
THEN
732 atom%orbitals%wrefene(k, l, 1:2) = w_semi
733 atom%orbitals%wrefchg(k, l, 1:2) = w_semi/100._dp
734 atom%orbitals%crefene(k, l, 1:2) = t_semi
735 atom%orbitals%reftype(k, l, 1:2) =
"SC"
736 wtot = wtot +
atom%weight*2._dp*(w_semi + w_semi/100._dp)
738 atom%orbitals%wrefene(k, l, 1:2) = w_valence
739 atom%orbitals%wrefchg(k, l, 1:2) = w_valence/100._dp
740 atom%orbitals%crefene(k, l, 1:2) = t_valence
741 atom%orbitals%reftype(k, l, 1:2) =
"VA"
742 wtot = wtot +
atom%weight*2._dp*(w_valence + w_valence/100._dp)
744 ntarget = ntarget + 4
746 atom%orbitals%tpsir0(k, 1:2) = t_psir0
747 atom%orbitals%wpsir0(k, 1:2) = w_psir0
748 wtot = wtot +
atom%weight*2._dp*w_psir0
753 atom%orbitals%refnod(k, l, 1:2) = real(k - 1, kind=
dp)
755 IF (k == 1 .AND.
atom%state%occa(l, k) /= 0._dp)
THEN
756 atom%orbitals%wrefnod(k, l, 1) = w_node
757 wtot = wtot +
atom%weight*w_node
759 IF (k == 1 .AND.
atom%state%occb(l, k) /= 0._dp)
THEN
760 atom%orbitals%wrefnod(k, l, 2) = w_node
761 wtot = wtot +
atom%weight*w_node
771 DO j1 = 1,
SIZE(atom_info, 2)
772 DO j2 = 1,
SIZE(atom_info, 2)
773 DO i1 = 1,
SIZE(atom_info, 1)
774 DO i2 = 1,
SIZE(atom_info, 1)
775 IF ((j1 > j2) .OR. (j1 == j2 .AND. i1 >= i2)) cycle
776 dener(2, j1, j2, i1, i2) = dener(2, j1, j1, i1, i1) - dener(2, j2, j2, i2, i2)
785 ALLOCATE (xi(ostate%nvar))
788 CALL put_pseudo_param(x, ppot%gth_pot, noopt_nlcc)
789 CALL pseudo_fit(atom_info, wfn_guess, ppot, ostate%f, wtot, pval, dener, w_ener, t_ener, .true.)
791 WRITE (iunit,
'(/," POWELL| Initial errors of target values")')
793 DO i = 1,
SIZE(atom_info, 1)
794 DO j = 1,
SIZE(atom_info, 2)
795 atom => atom_info(i, j)%atom
798 wfn_guess(i, j)%wfn =
atom%orbitals%wfn
800 WRITE (iunit,
'(/," Reference configuration ",T31,i5,T50," Method number ",T76,i5)') i, j
801 IF (
atom%state%multiplicity == -1)
THEN
803 WRITE (iunit,
'(" L N Occupation Eigenvalue [eV] dE [eV] dCharge ")')
805 np =
atom%state%maxn_calc(l)
808 oc =
atom%state%occupation(l, k)
809 eig =
atom%orbitals%ener(k, l)
810 deig = eig -
atom%orbitals%refene(k, l, 1)
811 peig = pval(1, k, l, j, i)/afun*100._dp
812 IF (pval(5, k, l, j, i) > 0.5_dp)
THEN
815 WRITE (pc1,
"(I2)") nint(peig)
818 atom%orbitals%rcmax(k, l, 1), l,
atom%basis)
819 drho = charge -
atom%orbitals%refchg(k, l, 1)
820 pchg = pval(2, k, l, j, i)/afun*100._dp
821 IF (pval(6, k, l, j, i) > 0.5_dp)
THEN
824 WRITE (pc2,
"(I2)") min(nint(pchg), 99)
826 pct =
atom%orbitals%reftype(k, l, 1)
827 WRITE (iunit,
'(I5,I5,F14.2,F21.10,A3,F11.6,"[",A2,"]",F13.6,"[",A2,"]")') &
828 l, k, oc, eig*
evolt, pct, deig*
evolt, pc1, drho, pc2
834 WRITE (iunit,
'(" L N Spin Occupation Eigenvalue [eV] dE [eV] dCharge ")')
836 np =
atom%state%maxn_calc(l)
839 oc =
atom%state%occa(l, k)
840 eig =
atom%orbitals%enera(k, l)
841 deig = eig -
atom%orbitals%refene(k, l, 1)
842 peig = pval(1, k, l, j, i)/afun*100._dp
843 IF (pval(5, k, l, j, i) > 0.5_dp)
THEN
846 WRITE (pc1,
"(I2)") nint(peig)
849 charge,
atom%orbitals%wfna(:, k, l),
atom%orbitals%rcmax(k, l, 1), l,
atom%basis)
850 drho = charge -
atom%orbitals%refchg(k, l, 1)
851 pchg = pval(2, k, l, j, i)/afun*100._dp
852 IF (pval(6, k, l, j, i) > 0.5_dp)
THEN
855 WRITE (pc2,
"(I2)") min(nint(pchg), 99)
857 pct =
atom%orbitals%reftype(k, l, 1)
858 WRITE (iunit,
'(I5,I5,2X,A5,F11.2,F19.10,A3,F10.6,"[",A2,"]",F12.6,"[",A2,"]")') &
859 l, k,
"alpha", oc, eig*
evolt, pct, deig*
evolt, pc1, drho, pc2
860 oc =
atom%state%occb(l, k)
861 eig =
atom%orbitals%enerb(k, l)
862 deig = eig -
atom%orbitals%refene(k, l, 2)
863 peig = pval(3, k, l, j, i)/afun*100._dp
864 IF (pval(7, k, l, j, i) > 0.5_dp)
THEN
867 WRITE (pc1,
"(I2)") nint(peig)
870 charge,
atom%orbitals%wfnb(:, k, l),
atom%orbitals%rcmax(k, l, 2), l,
atom%basis)
871 drho = charge -
atom%orbitals%refchg(k, l, 2)
872 pchg = pval(4, k, l, j, i)/afun*100._dp
873 IF (pval(8, k, l, j, i) > 0.5_dp)
THEN
876 WRITE (pc2,
"(I2)") min(nint(pchg), 99)
878 pct =
atom%orbitals%reftype(k, l, 2)
879 WRITE (iunit,
'(I5,I5,2X,A5,F11.2,F19.10,A3,F10.6,"[",A2,"]",F12.6,"[",A2,"]")') &
880 l, k,
" beta", oc, eig*
evolt, pct, deig*
evolt, pc1, drho, pc2
891 WRITE (iunit,
'(" Method/Econf 1 "," Method/Econf 2"," Delta Energy "," Error Energy "," Target")')
892 DO j1 = 1,
SIZE(atom_info, 2)
893 DO j2 = 1,
SIZE(atom_info, 2)
894 DO i1 = 1,
SIZE(atom_info, 1)
895 DO i2 = 1,
SIZE(atom_info, 1)
896 IF ((j1 > j2) .OR. (j1 == j2 .AND. i1 >= i2)) cycle
897 IF (abs(dener(1, j1, j1, i1, i1)) < 0.000001_dp) cycle
898 IF (abs(dener(2, j1, j1, i1, i1)) < 0.000001_dp) cycle
899 IF (abs(dener(1, j2, j2, i2, i2)) < 0.000001_dp) cycle
900 IF (abs(dener(2, j2, j2, i2, i2)) < 0.000001_dp) cycle
901 de = dener(2, j1, j2, i1, i2) - dener(1, j1, j2, i1, i2)
902 WRITE (iunit,
'(i6,i6,i10,i6,5X,F16.6,F19.6,F12.6)') &
903 j1, i1, j2, i2, dener(2, j1, j2, i1, i2), de, t_ener
911 WRITE (iunit,
'(" Number of target values reached: ",T69,i5," of ",i3)') &
912 int(sum(pval(5:8, :, :, :, :))), ntarget
917 WRITE (iunit,
'(" POWELL| Start optimization",I4," of total",I4,T60,"rhobeg = ",F12.8)') &
918 nre, nreinit, ostate%rhobeg
924 IF (ostate%state == 2)
THEN
925 CALL put_pseudo_param(x, ppot%gth_pot, noopt_nlcc)
926 CALL pseudo_fit(atom_info, wfn_guess, ppot, ostate%f, wtot, pval, dener, w_ener, t_ener, .false.)
927 fopt = min(fopt, ostate%f)
930 IF (ostate%state == -1)
EXIT
934 IF (ostate%nf == 2 .AND. iunit > 0)
THEN
935 WRITE (iunit,
'(" POWELL| Initial value of function",T61,F20.10)') ostate%f
937 IF (mod(ostate%nf, n10) == 0 .AND. iunit > 0 .AND. ostate%nf > 2)
THEN
938 WRITE (iunit,
'(" POWELL| Reached",i4,"% of maximal function calls",T61,F20.10)') &
939 int(real(ostate%nf,
dp)/real(ostate%maxfun,
dp)*100._dp), fopt
940 CALL put_pseudo_param(ostate%xopt, ppot%gth_pot, noopt_nlcc)
944 IF (fopt < 1.e-12_dp)
EXIT
947 WRITE (iw, *) ostate%nf, ostate%f, x(1:ostate%nvar)
952 dx = sqrt(sum((ostate%xopt(:) - xi(:))**2)/real(ostate%nvar, kind=
dp))
954 WRITE (iunit,
'(" POWELL| RMS average of variables",T69,F12.10)') dx
958 CALL put_pseudo_param(x, ppot%gth_pot, noopt_nlcc)
961 IF ((dx*dx) < ostate%rhoend)
EXIT
962 ostate%rhobeg = step_size_scaling*ostate%rhobeg
967 WRITE (iunit,
'(" POWELL| Number of function evaluations",T71,I10)') ostate%nf
968 WRITE (iunit,
'(" POWELL| Final value of function",T61,F20.10)') ostate%fopt
970 CALL put_pseudo_param(x, ppot%gth_pot, noopt_nlcc)
971 CALL pseudo_fit(atom_info, wfn_guess, ppot, ostate%f, wtot, pval, dener, w_ener, t_ener, .false.)
974 WRITE (iunit,
'(/," POWELL| Final errors of target values")')
975 DO i = 1,
SIZE(atom_info, 1)
976 DO j = 1,
SIZE(atom_info, 2)
977 atom => atom_info(i, j)%atom
979 WRITE (iunit,
'(/," Reference configuration ",T31,i5,T50," Method number ",T76,i5)') i, j
980 IF (
atom%state%multiplicity == -1)
THEN
982 WRITE (iunit,
'(" L N Occupation Eigenvalue [eV] dE [eV] dCharge ")')
984 np =
atom%state%maxn_calc(l)
987 oc =
atom%state%occupation(l, k)
988 eig =
atom%orbitals%ener(k, l)
989 deig = eig -
atom%orbitals%refene(k, l, 1)
990 peig = pval(1, k, l, j, i)/afun*100._dp
991 IF (pval(5, k, l, j, i) > 0.5_dp)
THEN
994 WRITE (pc1,
"(I2)") nint(peig)
997 charge,
atom%orbitals%wfn(:, k, l),
atom%orbitals%rcmax(k, l, 1), l,
atom%basis)
998 drho = charge -
atom%orbitals%refchg(k, l, 1)
999 pchg = pval(2, k, l, j, i)/afun*100._dp
1000 IF (pval(6, k, l, j, i) > 0.5_dp)
THEN
1003 WRITE (pc2,
"(I2)") min(nint(pchg), 99)
1005 pct =
atom%orbitals%reftype(k, l, 1)
1006 WRITE (iunit,
'(I5,I5,F14.2,F21.10,A3,F11.6,"[",A2,"]",F13.6,"[",A2,"]")') &
1007 l, k, oc, eig*
evolt, pct, deig*
evolt, pc1, drho, pc2
1011 np =
atom%state%maxn_calc(0)
1014 IF (abs(
atom%orbitals%tpsir0(k, 1)) > 1.e-14_dp)
THEN
1015 IF (abs(pv) > abs(
atom%orbitals%tpsir0(k, 1)))
THEN
1018 pv = 10._dp*(abs(pv) - abs(
atom%orbitals%tpsir0(k, 1)))
1020 pchg =
atom%weight*
atom%orbitals%wpsir0(k, 1)*pv*pv/afun
1022 pchg =
atom%weight*
atom%orbitals%wpsir0(k, 1)*pv*pv/afun*100._dp
1024 WRITE (iunit,
'(" s-states"," N=",I5,T40,"Wavefunction at r=0:",T64,F13.6,"[",I2,"]")') &
1029 WRITE (iunit,
'(" L N Spin Occupation Eigenvalue [eV] dE [eV] dCharge ")')
1031 np =
atom%state%maxn_calc(l)
1034 oc =
atom%state%occa(l, k)
1035 eig =
atom%orbitals%enera(k, l)
1036 deig = eig -
atom%orbitals%refene(k, l, 1)
1037 peig = pval(1, k, l, j, i)/afun*100._dp
1038 IF (pval(5, k, l, j, i) > 0.5_dp)
THEN
1041 WRITE (pc1,
"(I2)") nint(peig)
1044 charge,
atom%orbitals%wfna(:, k, l),
atom%orbitals%rcmax(k, l, 1), l,
atom%basis)
1045 drho = charge -
atom%orbitals%refchg(k, l, 1)
1046 pchg = pval(2, k, l, j, i)/afun*100._dp
1047 IF (pval(6, k, l, j, i) > 0.5_dp)
THEN
1050 WRITE (pc2,
"(I2)") min(nint(pchg), 99)
1052 pct =
atom%orbitals%reftype(k, l, 1)
1053 WRITE (iunit,
'(I5,I5,A,F11.2,F20.10,A3,F10.6,"[",A2,"]",F11.6,"[",A2,"]")') &
1054 l, k,
" alpha", oc, eig*
evolt, pct, deig*
evolt, pc1, drho, pc2
1055 oc =
atom%state%occb(l, k)
1056 eig =
atom%orbitals%enerb(k, l)
1057 deig = eig -
atom%orbitals%refene(k, l, 2)
1058 peig = pval(3, k, l, j, i)/afun*100._dp
1059 IF (pval(7, k, l, j, i) > 0.5_dp)
THEN
1062 WRITE (pc1,
"(I2)") nint(peig)
1065 charge,
atom%orbitals%wfnb(:, k, l),
atom%orbitals%rcmax(k, l, 2), l,
atom%basis)
1066 drho = charge -
atom%orbitals%refchg(k, l, 2)
1067 pchg = pval(4, k, l, j, i)/afun*100._dp
1068 IF (pval(8, k, l, j, i) > 0.5_dp)
THEN
1071 WRITE (pc2,
"(I2)") min(nint(pchg), 99)
1073 pct =
atom%orbitals%reftype(k, l, 2)
1074 WRITE (iunit,
'(I5,I5,A,F11.2,F20.10,A3,F10.6,"[",A2,"]",F11.6,"[",A2,"]")') &
1075 l, k,
" beta", oc, eig*
evolt, pct, deig*
evolt, pc1, drho, pc2
1079 np =
atom%state%maxn_calc(0)
1082 IF (abs(
atom%orbitals%tpsir0(k, 1)) > 1.e-14_dp)
THEN
1083 IF (abs(pv) > abs(
atom%orbitals%tpsir0(k, 1)))
THEN
1086 pv = 10._dp*(abs(pv) - abs(
atom%orbitals%tpsir0(k, 1)))
1088 pchg =
atom%weight*
atom%orbitals%wpsir0(k, 1)*pv*pv/afun
1090 pchg =
atom%weight*
atom%orbitals%wpsir0(k, 1)*pv*pv/afun*100._dp
1092 WRITE (iunit,
'(" s-states"," N=",I5,T35,"Alpha Wavefunction at r=0:",T64,F13.6,"[",I2,"]")') &
1096 IF (abs(
atom%orbitals%tpsir0(k, 2)) > 1.e-14_dp)
THEN
1097 IF (abs(pv) > abs(
atom%orbitals%tpsir0(k, 2)))
THEN
1100 pv = 10._dp*(abs(pv) - abs(
atom%orbitals%tpsir0(k, 2)))
1102 pchg =
atom%weight*
atom%orbitals%wpsir0(k, 2)*pv*pv/afun
1104 pchg =
atom%weight*
atom%orbitals%wpsir0(k, 2)*pv*pv/afun*100._dp
1106 WRITE (iunit,
'(" s-states"," N=",I5,T36,"Beta Wavefunction at r=0:",T64,F13.6,"[",I2,"]")') &
1116 WRITE (iunit,
'(" Method/Econf 1 "," Method/Econf 2"," Delta Energy "," Error Energy "," Target")')
1117 DO j1 = 1,
SIZE(atom_info, 2)
1118 DO j2 = 1,
SIZE(atom_info, 2)
1119 DO i1 = 1,
SIZE(atom_info, 1)
1120 DO i2 = 1,
SIZE(atom_info, 1)
1121 IF ((j1 > j2) .OR. (j1 == j2 .AND. i1 >= i2)) cycle
1122 IF (abs(dener(1, j1, j1, i1, i1)) < 0.000001_dp) cycle
1123 IF (abs(dener(2, j1, j1, i1, i1)) < 0.000001_dp) cycle
1124 IF (abs(dener(1, j2, j2, i2, i2)) < 0.000001_dp) cycle
1125 IF (abs(dener(2, j2, j2, i2, i2)) < 0.000001_dp) cycle
1126 de = dener(2, j1, j2, i1, i2) - dener(1, j1, j2, i1, i2)
1127 WRITE (iunit,
'(i6,i6,i10,i6,5X,F16.6,F19.6,F12.6)') j1, i1, j2, i2, dener(2, j1, j2, i1, i2), de, t_ener
1135 WRITE (iunit,
'(/," Number of target values reached: ",T69,i5," of ",i3)') int(sum(pval(5:8, :, :, :, :))), ntarget
1139 DEALLOCATE (x, pval, dener)
1141 DO i = 1,
SIZE(wfn_guess, 1)
1142 DO j = 1,
SIZE(wfn_guess, 2)
1143 IF (
ASSOCIATED(wfn_guess(i, j)%wfn))
THEN
1144 DEALLOCATE (wfn_guess(i, j)%wfn)
1148 DEALLOCATE (wfn_guess)
1150 IF (
ALLOCATED(xtob))
THEN
1168 SUBROUTINE opt_nlcc_param(atom_info, atom_refs, gthpot, iunit, preopt_nlcc)
1169 TYPE(atom_p_type),
DIMENSION(:, :),
POINTER :: atom_info, atom_refs
1170 TYPE(atom_gthpot_type),
INTENT(inout) :: gthpot
1171 INTEGER,
INTENT(in) :: iunit
1172 LOGICAL,
INTENT(in) :: preopt_nlcc
1174 INTEGER :: i, im, j, k, method
1175 REAL(kind=
dp) :: rcov, zcore, zrc, zrch
1176 TYPE(atom_type),
POINTER :: aref,
atom
1177 TYPE(opgrid_type),
POINTER :: corden, den, den1, den2, density
1178 TYPE(opmat_type),
POINTER :: denmat, dma, dmb
1180 cpassert(gthpot%nlcc)
1183 WRITE (iunit,
'(/," Core density information")')
1184 WRITE (iunit,
'(A,T37,A,T55,A,T75,A)')
" State Density:",
"Full",
"Rcov/2",
"Rcov/4"
1187 rcov =
ptable(atom_refs(1, 1)%atom%z)%covalent_radius*
bohr
1188 atom => atom_info(1, 1)%atom
1193 DO i = 1,
SIZE(atom_info, 1)
1194 DO j = 1,
SIZE(atom_info, 2)
1195 atom => atom_info(i, j)%atom
1196 aref => atom_refs(i, j)%atom
1198 NULLIFY (den, denmat)
1202 method =
atom%method_type
1203 SELECT CASE (method)
1206 atom%basis%nbas,
atom%state%core, &
1207 aref%state%maxl_occ, aref%state%maxn_occ)
1213 atom%basis%nbas,
atom%state%core, &
1214 aref%state%maxl_occ, aref%state%maxn_occ)
1216 atom%basis%nbas,
atom%state%core, &
1217 aref%state%maxl_occ, aref%state%maxn_occ)
1218 denmat%op = 0.5_dp*(dma%op + dmb%op)
1228 CALL atom_density(den%op, denmat%op,
atom%basis, aref%state%maxl_occ, typ=
"RHO")
1229 density%op = density%op + den%op
1230 zcore = integrate_grid(den%op,
atom%basis%grid)
1232 NULLIFY (den1, den2)
1237 DO k = 1,
atom%basis%grid%nr
1238 IF (
atom%basis%grid%rad(k) < 0.50_dp*rcov)
THEN
1239 den1%op(k) = den%op(k)
1241 IF (
atom%basis%grid%rad(k) < 0.25_dp*rcov)
THEN
1242 den2%op(k) = den%op(k)
1245 zrc = integrate_grid(den1%op,
atom%basis%grid)
1247 zrch = integrate_grid(den2%op,
atom%basis%grid)
1254 WRITE (iunit,
'(2I5,T31,F10.4,T51,F10.4,T71,F10.4)') i, j, zcore, zrc, zrch
1259 density%op = density%op/real(im, kind=
dp)
1261 IF (preopt_nlcc)
THEN
1262 gthpot%nexp_nlcc = 1
1264 gthpot%cval_nlcc = 0._dp
1265 gthpot%alpha_nlcc = 0._dp
1266 gthpot%nct_nlcc(1) = 1
1267 gthpot%cval_nlcc(1, 1) = 1.0_dp
1268 gthpot%alpha_nlcc(1) = gthpot%rc
1272 DO k = 1,
atom%basis%grid%nr
1273 IF (
atom%basis%grid%rad(k) > 0.25_dp*rcov)
THEN
1274 corden%op(k) = 0.0_dp
1277 zrc = integrate_grid(corden%op,
atom%basis%grid)
1279 gthpot%cval_nlcc(1, 1) = zrch/zrc
1284 END SUBROUTINE opt_nlcc_param
1297 SUBROUTINE density_fit(density, atom, n, aval, cval, co, aerr)
1298 TYPE(opgrid_type),
POINTER :: density
1299 TYPE(atom_type),
POINTER ::
atom
1300 INTEGER,
INTENT(IN) :: n
1301 REAL(
dp),
INTENT(IN) :: aval, cval
1302 REAL(
dp),
DIMENSION(:),
INTENT(INOUT) :: co
1303 REAL(
dp),
INTENT(OUT) :: aerr
1305 INTEGER :: i, info, ip, iq, k, nr
1306 INTEGER,
ALLOCATABLE,
DIMENSION(:) :: ipiv
1307 REAL(
dp) :: a, rk, zval
1308 REAL(
dp),
ALLOCATABLE,
DIMENSION(:) :: den, pe, uval
1309 REAL(
dp),
ALLOCATABLE,
DIMENSION(:, :) :: bf, smat, tval
1312 nr =
atom%basis%grid%nr
1313 ALLOCATE (bf(nr, n), den(nr))
1317 pe(i) = aval*cval**i
1320 rk =
atom%basis%grid%rad(k)
1321 bf(k, i) = exp(-a*rk**2)
1326 zval = integrate_grid(density%op,
atom%basis%grid)
1329 ALLOCATE (tval(n + 1, 1), uval(n), smat(n + 1, n + 1))
1331 uval(i) = (
pi/pe(i))**1.5_dp
1332 tval(i, 1) = integrate_grid(density%op, bf(:, i),
atom%basis%grid)
1334 tval(n + 1, 1) = zval
1338 smat(ip, iq) = (
pi/(pe(ip) + pe(iq)))**1.5_dp
1341 smat(1:n, n + 1) = uval(1:n)
1342 smat(n + 1, 1:n) = uval(1:n)
1343 smat(n + 1, n + 1) = 0._dp
1345 ALLOCATE (ipiv(n + 1))
1346 CALL lapack_sgesv(n + 1, 1, smat, n + 1, ipiv, tval, n + 1, info)
1349 co(1:n) = tval(1:n, 1)
1354 den(:) = den(:) + co(i)*bf(:, i)
1357 den(:) = (den(:) - density%op(:))**2
1358 aerr = sqrt(integrate_grid(den,
atom%basis%grid))
1360 DEALLOCATE (pe, bf, den)
1362 DEALLOCATE (tval, uval, smat)
1364 END SUBROUTINE density_fit
1374 SUBROUTINE basis_fit(atom_info, basis, pptype, afun, iw, penalty)
1375 TYPE(atom_p_type),
DIMENSION(:, :),
POINTER :: atom_info
1376 TYPE(atom_basis_type),
POINTER :: basis
1377 LOGICAL,
INTENT(IN) :: pptype
1378 REAL(
dp),
INTENT(OUT) :: afun
1379 INTEGER,
INTENT(IN) :: iw
1380 LOGICAL,
INTENT(IN) :: penalty
1382 INTEGER :: do_eric, do_erie, i, im, in, l, nm, nn, &
1384 LOGICAL :: eri_c, eri_e
1385 REAL(kind=
dp) :: amin
1386 TYPE(atom_integrals),
POINTER :: atint
1387 TYPE(atom_potential_type),
POINTER :: pot
1388 TYPE(atom_type),
POINTER ::
atom
1392 nn =
SIZE(atom_info, 1)
1393 nm =
SIZE(atom_info, 2)
1402 atom => atom_info(in, im)%atom
1403 IF (pptype .EQV.
atom%pp_calc)
THEN
1404 pot =>
atom%potential
1405 zval = atom_info(in, im)%atom%z
1406 reltyp = atom_info(in, im)%atom%relativistic
1407 do_eric = atom_info(in, im)%atom%coulomb_integral_type
1408 do_erie = atom_info(in, im)%atom%exchange_integral_type
1414 IF (
ASSOCIATED(pot))
EXIT
1417 CALL atom_int_setup(atint, basis, potential=pot, eri_coulomb=eri_c, eri_exchange=eri_e)
1421 NULLIFY (atint%tzora, atint%hdkh)
1431 atom => atom_info(in, im)%atom
1433 IF (pptype .EQV.
atom%pp_calc)
THEN
1437 afun = afun +
atom%energy%etot*
atom%weight
1446 DO i = 1, basis%nbas(l) - 1
1447 amin = minval(abs(basis%am(i, l) - basis%am(i + 1:basis%nbas(l), l)))
1448 amin = amin/basis%am(i, l)
1449 afun = afun + 10._dp*exp(-(20._dp*amin)**4)
1460 END SUBROUTINE basis_fit
1474 SUBROUTINE pseudo_fit(atom_info, wfn_guess, ppot, afun, wtot, pval, dener, wen, ten, init)
1475 TYPE(atom_p_type),
DIMENSION(:, :),
POINTER :: atom_info
1476 TYPE(wfn_init),
DIMENSION(:, :),
POINTER :: wfn_guess
1477 TYPE(atom_potential_type),
POINTER :: ppot
1478 REAL(
dp),
INTENT(OUT) :: afun
1479 REAL(
dp),
INTENT(IN) :: wtot
1480 REAL(
dp),
DIMENSION(:, :, 0:, :, :),
INTENT(INOUT) :: pval
1481 REAL(
dp),
DIMENSION(:, :, :, :, :),
INTENT(INOUT) :: dener
1482 REAL(
dp),
INTENT(IN) :: wen, ten
1483 LOGICAL,
INTENT(IN) :: init
1485 INTEGER :: i, i1, i2, j, j1, j2, k, l, n, node
1486 LOGICAL :: converged, noguess, shift
1487 REAL(kind=
dp) :: charge, de, fde, pv, rcov, rcov1, rcov2, &
1489 TYPE(atom_integrals),
POINTER :: pp_int
1490 TYPE(atom_type),
POINTER ::
atom
1494 dener(1, :, :, :, :) = 0._dp
1496 noguess = .NOT. init
1498 pp_int => atom_info(1, 1)%atom%integrals
1502 DO i = 1,
SIZE(atom_info, 1)
1503 DO j = 1,
SIZE(atom_info, 2)
1504 atom => atom_info(i, j)%atom
1507 atom%orbitals%wfn = wfn_guess(i, j)%wfn
1511 IF (.NOT. converged)
THEN
1513 IF (.NOT. shift)
THEN
1514 atom%orbitals%ener(:, :) = 1.5_dp*
atom%orbitals%ener(:, :) + 0.5_dp
1517 dener(1, j, j, i, i) =
atom%energy%etot
1518 DO l = 0,
atom%state%maxl_calc
1519 n =
atom%state%maxn_calc(l)
1521 IF (
atom%state%multiplicity == -1)
THEN
1523 rcov =
atom%orbitals%rcmax(k, l, 1)
1524 tv =
atom%orbitals%crefene(k, l, 1)
1525 de = abs(
atom%orbitals%ener(k, l) -
atom%orbitals%refene(k, l, 1))
1526 fde = get_error_value(de, tv)
1527 IF (fde < 1.e-8) pval(5, k, l, j, i) = 1._dp
1528 pv =
atom%weight*
atom%orbitals%wrefene(k, l, 1)*fde
1530 pval(1, k, l, j, i) = pv
1531 IF (
atom%orbitals%wrefchg(k, l, 1) > 0._dp)
THEN
1533 de = abs(charge -
atom%orbitals%refchg(k, l, 1))
1534 fde = get_error_value(de, 25._dp*tv)
1535 IF (fde < 1.e-8) pval(6, k, l, j, i) = 1._dp
1536 pv =
atom%weight*
atom%orbitals%wrefchg(k, l, 1)*fde
1538 pval(2, k, l, j, i) = pv
1540 IF (
atom%orbitals%wrefnod(k, l, 1) > 0._dp)
THEN
1542 afun = afun +
atom%weight*
atom%orbitals%wrefnod(k, l, 1)* &
1543 abs(real(node,
dp) -
atom%orbitals%refnod(k, l, 1))
1546 IF (
atom%orbitals%wpsir0(k, 1) > 0._dp)
THEN
1548 IF (abs(
atom%orbitals%tpsir0(k, 1)) > 1.e-14_dp)
THEN
1549 IF (abs(pv) > abs(
atom%orbitals%tpsir0(k, 1)))
THEN
1552 pv = 10._dp*(abs(pv) - abs(
atom%orbitals%tpsir0(k, 1)))
1554 pv =
atom%weight*
atom%orbitals%wpsir0(k, 1)*pv*pv
1556 pv =
atom%weight*
atom%orbitals%wpsir0(k, 1)*pv*pv*100._dp
1563 rcov1 =
atom%orbitals%rcmax(k, l, 1)
1564 rcov2 =
atom%orbitals%rcmax(k, l, 2)
1565 tv =
atom%orbitals%crefene(k, l, 1)
1566 de = abs(
atom%orbitals%enera(k, l) -
atom%orbitals%refene(k, l, 1))
1567 fde = get_error_value(de, tv)
1568 IF (fde < 1.e-8) pval(5, k, l, j, i) = 1._dp
1569 pv =
atom%weight*
atom%orbitals%wrefene(k, l, 1)*fde
1571 pval(1, k, l, j, i) = pv
1572 tv =
atom%orbitals%crefene(k, l, 2)
1573 de = abs(
atom%orbitals%enerb(k, l) -
atom%orbitals%refene(k, l, 2))
1574 fde = get_error_value(de, tv)
1575 IF (fde < 1.e-8) pval(7, k, l, j, i) = 1._dp
1576 pv =
atom%weight*
atom%orbitals%wrefene(k, l, 2)*fde
1578 pval(3, k, l, j, i) = pv
1579 IF (
atom%orbitals%wrefchg(k, l, 1) > 0._dp)
THEN
1581 de = abs(charge -
atom%orbitals%refchg(k, l, 1))
1582 fde = get_error_value(de, 25._dp*tv)
1583 IF (fde < 1.e-8) pval(6, k, l, j, i) = 1._dp
1584 pv =
atom%weight*
atom%orbitals%wrefchg(k, l, 1)*fde
1586 pval(2, k, l, j, i) = pv
1588 IF (
atom%orbitals%wrefchg(k, l, 2) > 0._dp)
THEN
1590 de = abs(charge -
atom%orbitals%refchg(k, l, 2))
1591 fde = get_error_value(de, 25._dp*tv)
1592 IF (fde < 1.e-8) pval(8, k, l, j, i) = 1._dp
1593 pv =
atom%weight*
atom%orbitals%wrefchg(k, l, 2)*fde
1595 pval(4, k, l, j, i) = pv
1597 IF (
atom%orbitals%wrefnod(k, l, 1) > 0._dp)
THEN
1599 afun = afun +
atom%weight*
atom%orbitals%wrefnod(k, l, 1)* &
1600 abs(real(node,
dp) -
atom%orbitals%refnod(k, l, 1))
1602 IF (
atom%orbitals%wrefnod(k, l, 2) > 0._dp)
THEN
1604 afun = afun +
atom%weight*
atom%orbitals%wrefnod(k, l, 2)* &
1605 abs(real(node,
dp) -
atom%orbitals%refnod(k, l, 2))
1608 IF (
atom%orbitals%wpsir0(k, 1) > 0._dp)
THEN
1610 IF (abs(
atom%orbitals%tpsir0(k, 1)) > 1.e-14_dp)
THEN
1611 IF (abs(pv) > abs(
atom%orbitals%tpsir0(k, 1)))
THEN
1614 pv = 10._dp*(abs(pv) - abs(
atom%orbitals%tpsir0(k, 1)))
1616 pv =
atom%weight*
atom%orbitals%wpsir0(k, 1)*pv*pv
1618 pv =
atom%weight*
atom%orbitals%wpsir0(k, 1)*pv*pv*100._dp
1622 IF (
atom%orbitals%wpsir0(k, 2) > 0._dp)
THEN
1624 IF (abs(
atom%orbitals%tpsir0(k, 2)) > 1.e-14_dp)
THEN
1625 IF (abs(pv) > abs(
atom%orbitals%tpsir0(k, 2)))
THEN
1628 pv = 10._dp*(abs(pv) - abs(
atom%orbitals%tpsir0(k, 2)))
1630 pv =
atom%weight*
atom%orbitals%wpsir0(k, 2)*pv*pv
1632 pv =
atom%weight*
atom%orbitals%wpsir0(k, 2)*pv*pv*100._dp
1645 DO j1 = 1,
SIZE(atom_info, 2)
1646 DO j2 = 1,
SIZE(atom_info, 2)
1647 DO i1 = 1,
SIZE(atom_info, 1)
1648 DO i2 = i1 + 1,
SIZE(atom_info, 1)
1649 IF ((j1 > j2) .OR. (j1 == j2 .AND. i1 >= i2)) cycle
1650 dener(1, j1, j2, i1, i2) = dener(1, j1, j1, i1, i1) - dener(1, j2, j2, i2, i2)
1651 de = abs(dener(2, j1, j2, i1, i2) - dener(1, j1, j2, i1, i2))
1652 fde = get_error_value(de, ten)
1653 afun = afun + wen*fde
1662 END SUBROUTINE pseudo_fit
1669 FUNCTION get_error_value(fval, ftarget)
RESULT(errval)
1670 REAL(kind=
dp),
INTENT(in) :: fval, ftarget
1671 REAL(kind=
dp) :: errval
1673 cpassert(fval >= 0.0_dp)
1675 IF (fval <= ftarget)
THEN
1678 errval = (fval - ftarget)/max(ftarget, 1.e-10_dp)
1679 errval = errval*errval
1682 END FUNCTION get_error_value
1690 SUBROUTINE get_pseudo_param(pvec, nval, gthpot, noopt_nlcc)
1691 REAL(kind=
dp),
DIMENSION(:),
INTENT(out) :: pvec
1692 INTEGER,
INTENT(out) :: nval
1693 TYPE(atom_gthpot_type),
INTENT(in) :: gthpot
1694 LOGICAL,
INTENT(in) :: noopt_nlcc
1696 INTEGER :: i, ival, j, l, n
1698 IF (gthpot%lsdpot)
THEN
1701 DO j = 1, gthpot%nexp_lsd
1703 pvec(ival) = rcpro(-1, gthpot%alpha_lsd(j))
1704 DO i = 1, gthpot%nct_lsd(j)
1706 pvec(ival) = gthpot%cval_lsd(i, j)
1712 pvec(ival) = rcpro(-1, gthpot%rc)
1713 DO i = 1, gthpot%ncl
1715 pvec(ival) = gthpot%cl(i)
1717 IF (gthpot%lpotextended)
THEN
1718 DO j = 1, gthpot%nexp_lpot
1720 pvec(ival) = rcpro(-1, gthpot%alpha_lpot(j))
1721 DO i = 1, gthpot%nct_lpot(j)
1723 pvec(ival) = gthpot%cval_lpot(i, j)
1727 IF (gthpot%nlcc .AND. (.NOT. noopt_nlcc))
THEN
1728 DO j = 1, gthpot%nexp_nlcc
1730 pvec(ival) = rcpro(-1, gthpot%alpha_nlcc(j))
1731 DO i = 1, gthpot%nct_nlcc(j)
1733 pvec(ival) = gthpot%cval_nlcc(i, j)
1738 IF (gthpot%nl(l) > 0)
THEN
1740 pvec(ival) = rcpro(-1, gthpot%rcnl(l))
1744 IF (gthpot%nl(l) > 0)
THEN
1749 pvec(ival) = gthpot%hnl(i, j, l)
1757 END SUBROUTINE get_pseudo_param
1765 SUBROUTINE put_pseudo_param(pvec, gthpot, noopt_nlcc)
1766 REAL(kind=
dp),
DIMENSION(:),
INTENT(in) :: pvec
1767 TYPE(atom_gthpot_type),
INTENT(inout) :: gthpot
1768 LOGICAL,
INTENT(in) :: noopt_nlcc
1770 INTEGER :: i, ival, j, l, n
1772 IF (gthpot%lsdpot)
THEN
1774 DO j = 1, gthpot%nexp_lsd
1776 gthpot%alpha_lsd(j) = rcpro(1, pvec(ival))
1777 DO i = 1, gthpot%nct_lsd(j)
1779 gthpot%cval_lsd(i, j) = pvec(ival)
1784 gthpot%rc = rcpro(1, pvec(ival))
1785 DO i = 1, gthpot%ncl
1787 gthpot%cl(i) = pvec(ival)
1789 IF (gthpot%lpotextended)
THEN
1790 DO j = 1, gthpot%nexp_lpot
1792 gthpot%alpha_lpot(j) = rcpro(1, pvec(ival))
1793 DO i = 1, gthpot%nct_lpot(j)
1795 gthpot%cval_lpot(i, j) = pvec(ival)
1799 IF (gthpot%nlcc .AND. (.NOT. noopt_nlcc))
THEN
1800 DO j = 1, gthpot%nexp_nlcc
1802 gthpot%alpha_nlcc(j) = rcpro(1, pvec(ival))
1803 DO i = 1, gthpot%nct_nlcc(j)
1805 gthpot%cval_nlcc(i, j) = pvec(ival)
1810 IF (gthpot%nl(l) > 0)
THEN
1812 gthpot%rcnl(l) = rcpro(1, pvec(ival))
1816 IF (gthpot%nl(l) > 0)
THEN
1821 gthpot%hnl(i, j, l) = pvec(ival)
1828 END SUBROUTINE put_pseudo_param
1838 FUNCTION rcpro(id, xval)
RESULT(yval)
1839 INTEGER,
INTENT(IN) :: id
1840 REAL(
dp),
INTENT(IN) :: xval
1846 yval = 2.0_dp*tanh(0.1_dp*xval)**2
1847 ELSE IF (id == -1)
THEN
1848 x1 = sqrt(xval/2.0_dp)
1849 cpassert(x1 <= 1._dp)
1850 x2 = 0.5_dp*log((1._dp + x1)/(1._dp - x1))
1867 TYPE(atom_type),
POINTER ::
atom
1868 INTEGER,
INTENT(IN) :: num_gau, num_pol, iunit
1869 TYPE(section_vals_type),
OPTIONAL,
POINTER :: powell_section
1870 REAL(kind=
dp),
DIMENSION(:),
OPTIONAL :: results
1872 REAL(kind=
dp),
PARAMETER :: t23 = 2._dp/3._dp
1874 INTEGER :: i, ig, ip, iw, j, n10
1875 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: x
1876 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:, :) :: co
1877 TYPE(opgrid_type),
POINTER :: density
1882 cpassert(num_pol*num_gau > 0)
1884 ALLOCATE (co(num_pol + 1, num_gau), x(num_pol*num_gau + num_gau))
1891 atom%state%maxl_occ,
atom%state%maxn_occ)
1894 density%op = t23*(0.3_dp*(3.0_dp*
pi*
pi)**t23)*density%op**t23
1898 ostate%nvar = num_pol*num_gau + num_gau
1900 co(1, i) = 0.5_dp + real(i - 1, kind=
dp)
1902 DO j = 3, num_pol + 1
1906 CALL putvar(x, co, num_pol, num_gau)
1908 IF (
PRESENT(powell_section))
THEN
1913 ostate%rhoend = 1.e-8_dp
1914 ostate%rhobeg = 5.e-2_dp
1915 ostate%maxfun = 1000
1922 WRITE (iunit,
'(/," ",13("*")," Approximated Nonadditive Kinetic Energy Functional ",14("*"))')
1923 WRITE (iunit,
'(" POWELL| Start optimization procedure")')
1929 IF (ostate%state == 2)
THEN
1930 CALL getvar(x, co, num_pol, num_gau)
1931 CALL kgpot_fit(density, num_gau, num_pol, co, ostate%f)
1934 IF (ostate%state == -1)
EXIT
1938 IF (mod(ostate%nf, n10) == 0 .AND. iunit > 0)
THEN
1939 WRITE (iunit,
'(" POWELL| Reached",i4,"% of maximal function calls",T66,G15.7)') &
1940 int(real(ostate%nf,
dp)/real(ostate%maxfun,
dp)*100._dp), ostate%fopt
1947 CALL getvar(x, co, num_pol, num_gau)
1952 WRITE (iunit,
'(" POWELL| Number of function evaluations",T71,I10)') ostate%nf
1953 WRITE (iunit,
'(" POWELL| Final value of function",T61,G20.10)') ostate%fopt
1954 WRITE (iunit,
'(" Optimized local potential of approximated nonadditive kinetic energy functional")')
1956 WRITE (iunit,
'(I2,T15,"Gaussian polynomial expansion",T66,"Rc=",F12.4)') ig, co(1, ig)
1957 WRITE (iunit,
'(T15,"Coefficients",T33,4F12.4)') (co(1 + ip, ig), ip=1, num_pol)
1961 CALL open_file(file_name=
"FIT_RESULT", file_status=
"UNKNOWN", file_action=
"WRITE", unit_number=iw)
1963 WRITE (iw, *) num_gau, num_pol
1965 WRITE (iw,
'(T10,F12.4,6X,4F12.4)') (co(ip, ig), ip=1, num_pol + 1)
1969 IF (
PRESENT(results))
THEN
1970 cpassert(
SIZE(results) >=
SIZE(x))
1986 SUBROUTINE kgpot_fit(kgpot, ng, np, cval, aerr)
1987 TYPE(opgrid_type),
POINTER :: kgpot
1988 INTEGER,
INTENT(IN) :: ng, np
1989 REAL(
dp),
DIMENSION(:, :),
INTENT(INOUT) :: cval
1990 REAL(
dp),
INTENT(OUT) :: aerr
1992 INTEGER :: i, ig, ip, n
1993 REAL(kind=
dp) :: pc, rc
1994 REAL(kind=
dp),
ALLOCATABLE,
DIMENSION(:) :: pval
2001 rc = kgpot%grid%rad(i)/cval(1, ig)
2004 pc = pc + cval(ip + 1, ig)*rc**(2*ip - 2)
2006 pval(i) = pval(i) + pc*exp(-0.5_dp*rc*rc)
2009 pval(1:n) = (pval(1:n) - kgpot%op(1:n))**2
2010 aerr =
fourpi*sum(pval(1:n)*kgpot%grid%wr(1:n))
2014 END SUBROUTINE kgpot_fit
2023 PURE SUBROUTINE getvar(xvar, cvar, np, ng)
2024 REAL(kind=
dp),
DIMENSION(:),
INTENT(in) :: xvar
2025 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(inout) :: cvar
2026 INTEGER,
INTENT(IN) :: np, ng
2028 INTEGER :: ig, ii, ip
2033 cvar(1, ig) = xvar(ii)
2036 cvar(ip + 1, ig) = xvar(ii)**2
2040 END SUBROUTINE getvar
2049 PURE SUBROUTINE putvar(xvar, cvar, np, ng)
2050 REAL(kind=
dp),
DIMENSION(:),
INTENT(inout) :: xvar
2051 REAL(kind=
dp),
DIMENSION(:, :),
INTENT(in) :: cvar
2052 INTEGER,
INTENT(IN) :: np, ng
2054 INTEGER :: ig, ii, ip
2059 xvar(ii) = cvar(1, ig)
2062 xvar(ii) = sqrt(abs(cvar(ip + 1, ig)))
2066 END SUBROUTINE putvar
subroutine, public calculate_atom(atom, iw, noguess, converged)
General routine to perform electronic structure atomic calculations.
routines that fit parameters for /from atomic calculations
subroutine, public atom_fit_pseudo(atom_info, atom_refs, ppot, iunit, powell_section)
...
subroutine, public atom_fit_density(atom, num_gto, norder, iunit, powell_section, results)
Fit the atomic electron density using a geometrical Gaussian basis set.
subroutine, public atom_fit_kgpot(atom, num_gau, num_pol, iunit, powell_section, results)
...
subroutine, public atom_fit_basis(atom_info, basis, pptype, iunit, powell_section)
...
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).
Routines that print various information about an atomic kind.
subroutine, public atom_print_basis(atom_basis, iw, title)
Print atomic basis set.
subroutine, public atom_print_basis_file(atom_basis, wfn)
Print the optimized atomic basis set into a file.
subroutine, public atom_write_pseudo_param(gthpot, iunit)
Print GTH pseudo-potential parameters.
Define the atom type and its sub types.
integer, parameter, public gto_basis
integer, parameter, public sto_basis
subroutine, public release_opmat(opmat)
...
subroutine, public release_opgrid(opgrid)
...
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 create_opmat(opmat, n, lmax)
...
subroutine, public create_opgrid(opgrid, grid)
...
Some basic routines for atomic calculations.
subroutine, public atom_condnumber(basis, crad, iw)
Print condition numbers of the given atomic basis set.
pure subroutine, public atom_denmat(pmat, wfn, nbas, occ, maxl, maxn)
Calculate a density matrix using atomic orbitals.
pure subroutine, public atom_orbital_charge(charge, wfn, rcov, l, basis)
...
pure logical function, public atom_consistent_method(method, multiplicity)
Check that the atomic multiplicity is consistent with the electronic structure method.
subroutine, public atom_core_density(corden, potential, typ, rr)
...
pure integer function, dimension(0:lmat), public get_maxn_occ(occupation)
Return the maximum principal quantum number of occupied orbitals.
subroutine, public atom_density(density, pmat, basis, maxl, typ, rr)
Map the electron density on an atomic radial grid.
pure subroutine, public atom_orbital_max(rmax, wfn, rcov, l, basis)
...
subroutine, public atom_completeness(basis, zv, iw)
Estimate completeness of the given atomic basis set.
pure subroutine, public atom_orbital_nodes(node, wfn, rcov, l, basis)
...
pure subroutine, public atom_wfnr0(value, wfn, basis)
...
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
Interface to the LAPACK F77 library.
Definition of mathematical constants and functions.
real(kind=dp), parameter, public pi
real(kind=dp), parameter, public fourpi
real(kind=dp), dimension(0:maxfac), parameter, public fac
Periodic Table related data definitions.
type(atom), dimension(0:nelem), public ptable
Definition of physical constants:
real(kind=dp), parameter, public evolt
real(kind=dp), parameter, public bohr
subroutine, public powell_optimize(n, x, optstate)
...