33 use,
intrinsic :: iso_fortran_env
69 integer,
public,
parameter :: &
70 XC_INV_METHOD_TWO_PARTICLE = 1, &
76 integer,
public,
parameter :: &
77 XC_KS_INVERSION_NONE = 1, &
82 integer,
public,
parameter :: &
83 XC_ASYMPTOTICS_NONE = 1, &
86 integer,
parameter :: &
91 integer,
public :: method
92 integer :: level = xc_ks_inversion_none
93 integer :: asymptotics
94 real(real64),
allocatable :: vhxc_previous_step(:,:)
95 type(states_elec_t),
public :: aux_st
96 type(hamiltonian_elec_t) :: aux_hm
97 type(eigensolver_t),
public :: eigensolver
100 type(partner_list_t) :: ext_partners
102 integer,
public :: max_iter
110 type(xc_ks_inversion_t),
intent(inout) :: ks_inv
111 type(namespace_t),
intent(in) :: namespace
112 type(grid_t),
intent(inout) :: gr
113 type(ions_t),
intent(inout) :: ions
114 type(states_elec_t),
intent(in) :: st
115 type(xc_t),
intent(in) :: xc
116 type(multicomm_t),
intent(in) :: mc
117 class(space_t),
intent(in) :: space
118 type(kpoints_t),
intent(in) :: kpoints
120 class(lasers_t),
pointer :: lasers
124 if(gr%parallel_in_domains)
then
148 if (ks_inv%method < xc_inv_method_two_particle &
180 call parse_variable(namespace,
'KSInversionAsymptotics', xc_asymptotics_none, ks_inv%asymptotics)
181 if (ks_inv%asymptotics /= xc_asymptotics_none .and. space%dim > 1)
then
190 if (ks_inv%level /= xc_ks_inversion_none)
then
204 if(lasers%no_lasers > 0)
then
205 call ks_inv%ext_partners%add(lasers)
211 call hamiltonian_elec_init(ks_inv%aux_hm, namespace, space, gr, ions, ks_inv%ext_partners, ks_inv%aux_st, &
213 call eigensolver_init(ks_inv%eigensolver, namespace, gr, ks_inv%aux_st, mc, space)
229 if (ks_inv%level /= xc_ks_inversion_none)
then
235 call iter%start(ks_inv%ext_partners)
236 do while (iter%has_next())
237 partner => iter%get_next()
238 safe_deallocate_p(partner)
240 call ks_inv%ext_partners%empty()
251 integer,
optional,
intent(in) :: iunit
252 type(
namespace_t),
optional,
intent(in) :: namespace
254 if (ks_inversion%level == xc_ks_inversion_none)
return
270 subroutine invertks_2part(ks_inv, target_rho, nspin, aux_hm, gr, st, eigensolver, &
271 namespace, space, ext_partners)
273 real(real64),
intent(in) :: target_rho(:,:)
274 integer,
intent(in) :: nspin
276 type(
grid_t),
intent(in) :: gr
280 class(
space_t),
intent(in) :: space
283 integer :: asym1, asym2, ip, ispin
285 real(real64) :: rr, shift
286 real(real64),
parameter :: smalldensity = 5e-12_real64
287 real(real64),
allocatable :: sqrtrho(:,:), laplace(:,:), vks(:,:), x_shifted(:)
292 assert(.not. gr%parallel_in_domains)
296 assert(space%periodic_dim == 0)
299 assert(space%dim == 1)
303 safe_allocate(sqrtrho(1:gr%np_part, 1:nspin))
304 safe_allocate(vks(1:np, 1:nspin))
305 safe_allocate(laplace(1:gr%np, 1:nspin))
307 if (any(target_rho(:,:) < -
m_epsilon))
then
308 write(
message(1),*)
"Target density has negative points. min value = ", minval(target_rho(:,:))
317 sqrtrho(ip, ispin) =
sqrt(target_rho(ip, ispin))
324 if (space%dim == 1)
then
329 call spline_fit(gr%np, gr%x(:, 1), sqrtrho(:, ispin), spl)
339 safe_deallocate_a(x_shifted)
343 laplace(ip, ispin) =
spline_eval(lapl_spl, gr%x(ip, 1))
369 if (target_rho(ip, ispin) < smalldensity)
then
370 vks(ip, ispin) = aux_hm%ep%vpsl(ip) + aux_hm%ks_pot%vhartree(ip)
373 if (target_rho(np-ip+1, ispin) < smalldensity)
then
374 vks(np-ip+1, ispin) = aux_hm%ep%vpsl(np-ip+1) + aux_hm%ks_pot%vhartree(np-ip+1)
382 do ip = asym1+1, asym2-1
383 vks(ip, ispin) =
m_half * laplace(ip, ispin) / (sqrtrho(ip, ispin) +
m_tiny)
389 aux_hm%ks_pot%vxc(ip, ispin) = vks(ip, ispin) - aux_hm%ep%vpsl(ip) - aux_hm%ks_pot%vhartree(ip)
400 aux_hm%ks_pot%vxc(ip, ispin) = -st%qtot/
sqrt(rr**2 +
m_one)
405 call mesh_r(gr, asym1+1, rr)
406 shift = aux_hm%ks_pot%vxc(asym1+1, ispin) + st%qtot/
sqrt(rr**2 +
m_one)
409 do ip = asym1+1, asym2-1
410 aux_hm%ks_pot%vxc(ip, ispin) = aux_hm%ks_pot%vxc(ip, ispin) - shift
415 call mesh_r(gr, asym2-1, rr)
416 shift = aux_hm%ks_pot%vxc(asym2-1, ispin) + st%qtot/
sqrt(rr**2 +
m_one)
421 aux_hm%ks_pot%vxc(ip, ispin) = aux_hm%ks_pot%vxc(ip, ispin) - shift
428 aux_hm%ks_pot%vxc(ip, ispin) = -st%qtot/
sqrt(rr**2 +
m_one)
436 if (ks_inv%asymptotics == xc_asymptotics_none)
then
443 shift = aux_hm%ks_pot%vxc(asym1+1, ispin)
445 do ip = asym1+1, asym2-1
446 aux_hm%ks_pot%vxc(ip, ispin) = aux_hm%ks_pot%vxc(ip, ispin) - shift
450 shift = aux_hm%ks_pot%vxc(asym2-1, ispin)
453 aux_hm%ks_pot%vxc(ip, ispin) = aux_hm%ks_pot%vxc(ip, ispin) - shift
461 aux_hm%ks_pot%vhxc(1:gr%np, ispin) = aux_hm%ks_pot%vxc(1:gr%np, ispin) + aux_hm%ks_pot%vhartree(1:gr%np)
467 safe_deallocate_a(sqrtrho)
468 safe_deallocate_a(laplace)
469 safe_deallocate_a(vks)
477 type(
grid_t),
intent(in) :: gr
478 class(
space_t),
intent(in) :: space
483 integer,
intent(in) :: max_iter
486 integer :: iter, nst_conv
490 call hm%update(gr, namespace, space, ext_partners)
493 nst_conv =
floor(st%qtot/st%smear%el_per_state)
495 do iter = 1, max_iter
496 call eigensolver%run(namespace, gr, st, hm, 1, converged, nst_conv)
507 message(1) =
'All states converged.'
510 message(1) =
'Some of the states are not fully converged!'
514 write(
message(1),
'(a, e17.6)')
'Criterion = ', eigensolver%tolerance
526 character(len=50) :: str
530 write(str,
'(a,i5)')
'Kohn-Sham inversion eigensolver iteration #', iter
533 write(
message(1),
'(a,i6,a,i6)')
'Converged states: ', minval(eigensolver%converged(1:st%nik))
537 eigensolver%diff, compact = .
true., namespace=namespace)
550 subroutine invertks_iter(ks_inv, target_rho, namespace, space, ext_partners, nspin, aux_hm, gr, &
553 real(real64),
intent(in) :: target_rho(:,:)
554 type(
grid_t),
intent(in) :: gr
556 class(
space_t),
intent(in) :: space
561 integer,
intent(in) :: nspin
563 integer :: ip, ispin, ierr, asym1, asym2
564 integer :: iunit, verbosity, counter, np
567 real(real64) :: rr, shift
568 real(real64) :: alpha, beta
569 real(real64) :: mu, npower, npower_in
570 real(real64) :: convdensity, diffdensity
571 real(real64),
allocatable :: vhxc(:,:)
572 real(real64),
parameter :: small_density = 5e-12_real64
574 character(len=256) :: fname
579 assert(space%dim == 1)
591 call parse_variable(namespace,
'InvertKSConvAbsDens', 1e-5_real64, convdensity)
600 call parse_variable(namespace,
'InvertKSStellaBeta', 1e-6_real64, beta)
609 call parse_variable(namespace,
'InvertKSStellaAlpha', 0.25_real64, alpha)
628 call parse_variable(namespace,
'InvertKSGodbyPower', 0.05_real64, npower_in)
647 if (verbosity < 0 .or. verbosity > 2)
then
661 safe_allocate(vhxc(1:np, 1:nspin))
662 call lalg_copy(np, nspin, aux_hm%ks_pot%vhxc, vhxc)
664 if (verbosity == 1 .or. verbosity == 2)
then
665 iunit =
io_open(
'InvertKSconvergence', namespace, action =
'write')
674 if (abs(st%rho(ip, ispin)-target_rho(ip, ispin)) > diffdensity)
then
675 diffdensity = abs(st%rho(ip, ispin)-target_rho(ip, ispin))
682 do while(diffdensity > convdensity .and. counter < max_iter)
684 counter = counter + 1
686 if (verbosity == 2)
then
687 write(fname,
'(i6.6)') counter
689 gr, aux_hm%ks_pot%vhxc(:,1),
units_out%energy, ierr)
691 gr, st%rho(:,1),
units_out%length**(-space%dim), ierr)
698 select case(ks_inv%method)
704 vhxc(ip, ispin) = vhxc(ip, ispin) &
705 + ((st%rho(ip, ispin) - target_rho(ip, ispin))/(target_rho(ip, ispin) + beta))*alpha
714 beta = diffdensity*1e-3_real64
717 alpha = max(0.05_real64,
m_half - diffdensity*100.0_real64*0.45_real64)
718 write(
message(1),
'(a,2E15.4,3I8, 2E15.4)') &
719 ' KSinversion: diffdensity, convdensity, imax, counter, max_iter, alpha, beta ', &
720 diffdensity, convdensity, imax, counter, max_iter, alpha, beta
727 vhxc(ip, ispin) = vhxc(ip, ispin) &
728 + ((st%rho(ip, ispin) - target_rho(ip, ispin))/(target_rho(ip, ispin) + beta))*alpha
739 write(
message(1),
'(a,2E15.4,3I8, 2E15.4)') &
740 ' KSinversion: diffdensity, convdensity, imax, counter, max_iter, power, mu ', &
741 diffdensity, convdensity, imax, counter, max_iter, npower, mu
748 vhxc(ip, ispin) = vhxc(ip, ispin) &
749 + (st%rho(ip, ispin)**npower - target_rho(ip, ispin)**npower)*mu
764 assert(.not. gr%parallel_in_domains)
768 if (abs(st%rho(ip, ispin)-target_rho(ip, ispin)) > diffdensity)
then
769 diffdensity = abs(st%rho(ip, ispin)-target_rho(ip, ispin))
775 if (verbosity == 1 .or. verbosity == 2)
then
776 write(iunit,
'(i6.6)', advance =
'no') counter
777 write(iunit,
'(es18.10)') diffdensity
784 call lalg_copy(np, nspin, vhxc, aux_hm%ks_pot%vhxc)
786 aux_hm%ks_pot%vxc(1:np, ispin) = vhxc(1:np, ispin) - aux_hm%ks_pot%vhartree(1:np)
795 call lalg_copy(np, nspin, aux_hm%ks_pot%vxc, vhxc)
803 if (target_rho(ip, ispin) < small_density)
then
805 vhxc(ip, ispin) = -st%qtot/
sqrt(rr**2 +
m_one)
808 if (target_rho(np-ip+1, ispin) < small_density)
then
814 call mesh_r(gr, asym1+1, rr)
815 shift = vhxc(asym1+1, ispin) + st%qtot/
sqrt(rr**2 +
m_one)
818 do ip = asym1+1, asym2-1
819 vhxc(ip, ispin) = vhxc(ip, ispin) - shift
823 call mesh_r(gr, asym2-1, rr)
824 shift = vhxc(asym2-1, ispin) + st%qtot/
sqrt(rr**2 +
m_one)
829 vhxc(ip, ispin) = vhxc(ip, ispin) - shift
836 vhxc(ip, ispin) = -st%qtot/
sqrt(rr**2 +
m_one)
843 call lalg_copy(np, nspin, vhxc, aux_hm%ks_pot%vxc)
845 aux_hm%ks_pot%vhxc(1:np, ispin) = vhxc(1:np, ispin) + aux_hm%ks_pot%vhartree(1:np)
854 write(
message(1),
'(a,I8)')
"Iterative KS inversion, iterations needed:", counter
859 safe_deallocate_a(vhxc)
865 subroutine xc_ks_inversion_calc(ks_inversion, namespace, space, gr, hm, ext_partners, st, vxc, time)
868 class(
space_t),
intent(in) :: space
869 type(
grid_t),
intent(in) :: gr
873 real(real64),
intent(inout) :: vxc(:,:)
874 real(real64),
optional,
intent(in) :: time
878 if (ks_inversion%level == xc_ks_inversion_none)
return
882 if (
present(time) .and. time >
m_zero)
then
883 write(
message(1),
'(A,F18.12)')
'xc_ks_inversion_calc - time:', time
887 ks_inversion%aux_hm%energy%intnvxc =
m_zero
888 ks_inversion%aux_hm%energy%hartree =
m_zero
889 ks_inversion%aux_hm%energy%exchange =
m_zero
890 ks_inversion%aux_hm%energy%correlation =
m_zero
892 ks_inversion%aux_hm%ks_pot%vhartree = hm%ks_pot%vhartree
894 if (
present(time) .and. time >
m_zero)
then
895 call lalg_copy(gr%np, st%d%nspin, ks_inversion%vhxc_previous_step, ks_inversion%aux_hm%ks_pot%vhxc)
899 do ispin = 1, st%d%nspin
900 ks_inversion%aux_hm%ks_pot%vxc(1:gr%np, ispin) =
m_zero
901 ks_inversion%aux_hm%ks_pot%vhxc(1:gr%np, ispin) = hm%ks_pot%vhartree(1:gr%np) &
902 + ks_inversion%aux_hm%ks_pot%vxc(1:gr%np, ispin)
910 call lalg_copy(gr%np, hm%ep%vpsl, ks_inversion%aux_hm%ep%vpsl)
915 ks_inversion%eigensolver, ks_inversion%aux_st, ks_inversion%aux_hm, ks_inversion%max_iter)
920 select case (ks_inversion%method)
922 case (xc_inv_method_two_particle)
923 call invertks_2part(ks_inversion, ks_inversion%aux_st%rho, st%d%nspin, ks_inversion%aux_hm, gr, &
924 ks_inversion%aux_st, ks_inversion%eigensolver, namespace, space, ext_partners)
927 call invertks_iter(ks_inversion, st%rho, namespace, space, ext_partners, st%d%nspin, ks_inversion%aux_hm, gr, &
928 ks_inversion%aux_st, ks_inversion%eigensolver)
934 do ispin = 1, st%d%nspin
935 call lalg_axpy(gr%np, -
m_one, hm%ks_pot%vhartree, ks_inversion%aux_hm%ks_pot%vhxc(:,ispin))
938 vxc = ks_inversion%aux_hm%ks_pot%vxc
941 if (
present(time))
then
942 ks_inversion%vhxc_previous_step = ks_inversion%aux_hm%ks_pot%vhxc
constant times a vector plus a vector
Copies a vector x, to a vector y.
Some operations may be done for one spline-function, or for an array of them.
double floor(double __x) __attribute__((__nothrow__
This module implements a calculator for the density and defines related functions.
subroutine, public density_calc(st, gr, density, istin)
Computes the density from the orbitals in st.
This module calculates the derivatives (gradients, Laplacians, etc.) of a function.
subroutine, public dderivatives_lapl(der, ff, op_ff, ghost_update, set_bc, factor)
apply the Laplacian to a mesh function
subroutine, public eigensolver_init(eigens, namespace, gr, st, mc, space)
subroutine, public eigensolver_end(eigens)
real(real64), parameter, public m_zero
real(real64), parameter, public m_tiny
real(real64), parameter, public m_epsilon
real(real64), parameter, public m_half
real(real64), parameter, public m_one
This module implements the underlying real-space grid.
subroutine, public hamiltonian_elec_end(hm)
subroutine, public hamiltonian_elec_init(hm, namespace, space, gr, ions, ext_partners, st, theory_level, xc, mc, kpoints, need_exchange, xc_photons)
This module defines classes and functions for interaction partners.
subroutine, public dio_function_output(how, dir, fname, namespace, space, mesh, ff, unit, ierr, pos, atoms, grp, root)
Top-level IO routine for functions defined on the mesh.
integer(int64) function, public io_function_fill_how(where)
Use this function to quickly plot functions for debugging purposes: call dio_function_output(io_funct...
subroutine, public io_close(iunit, grp)
integer function, public io_open(file, namespace, action, status, form, position, die, recl, grp)
A module to handle KS potential, without the external potential.
integer, parameter, public kohn_sham_dft
subroutine, public lasers_check_symmetries(this, kpoints)
subroutine, public lasers_parse_external_fields(this)
subroutine, public lasers_generate_potentials(this, mesh, space, latt)
This module defines the meshes, which are used in Octopus.
pure subroutine, public mesh_r(mesh, ip, rr, origin, coords)
return the distance to the origin for a given grid point
subroutine, public messages_print_with_emphasis(msg, iunit, namespace)
subroutine, public messages_not_implemented(feature, namespace)
character(len=512), private msg
subroutine, public messages_warning(no_lines, all_nodes, namespace)
subroutine, public messages_obsolete_variable(namespace, name, rep)
subroutine, public messages_new_line()
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
subroutine, public messages_input_error(namespace, var, details, row, column)
subroutine, public messages_experimental(name, namespace)
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
This module handles the communicators for the various parallelization strategies.
subroutine, public spline_fit(nrc, rofi, ffit, spl, threshold)
subroutine, public spline_der2(spl, dspl, threshold, grid)
Returns a spline that contains the second derivative of the original spline.
real(real64) function, public spline_eval(spl, x)
subroutine, public spline_generate_shifted_grid(this, x_shifted)
This module handles spin dimensions of the states and the k-point distribution.
This module defines routines to write information about states.
subroutine, public states_elec_write_eigenvalues(nst, st, space, kpoints, error, st_start, compact, iunit, namespace)
write the eigenvalues for some states to a file.
subroutine, public states_elec_densities_init(st, gr)
subroutine, public states_elec_end(st)
finalize the states_elec_t object
subroutine, public states_elec_allocate_wfns(st, mesh, wfs_type, skip, packed)
Allocates the KS wavefunctions defined within a states_elec_t structure.
subroutine, public states_elec_copy(stout, stin, exclude_wfns, exclude_eigenval, special)
make a (selective) copy of a states_elec_t object
subroutine, public states_elec_generate_random(st, mesh, kpoints, ist_start_, ist_end_, ikpt_start_, ikpt_end_, normalized)
randomize states
brief This module defines the class unit_t which is used by the unit_systems_oct_m module.
This module defines the unit system, used for input and output.
type(unit_system_t), public units_out
subroutine, public xc_ks_inversion_end(ks_inv)
integer, parameter, public xc_inv_method_iter_godby
subroutine, public invertks_2part(ks_inv, target_rho, nspin, aux_hm, gr, st, eigensolver, namespace, space, ext_partners)
Given a density, it performs the Kohn-Sham inversion, assuming a two-particle, one orbital case.
integer, parameter, public xc_ks_inversion_td_exact
integer, parameter, public xc_ks_inversion_adiabatic
integer, parameter, public xc_inv_method_vs_iter
integer, parameter, public xc_asymptotics_sc
subroutine, public xc_ks_inversion_write_info(ks_inversion, iunit, namespace)
subroutine, public invertks_iter(ks_inv, target_rho, namespace, space, ext_partners, nspin, aux_hm, gr, st, eigensolver)
Iterative inversion of KS potential from the density.
subroutine, public xc_ks_inversion_init(ks_inv, namespace, gr, ions, st, xc, mc, space, kpoints)
integer, parameter, public xc_inv_method_iter_stella
subroutine, public xc_ks_inversion_calc(ks_inversion, namespace, space, gr, hm, ext_partners, st, vxc, time)
subroutine, public invertks_update_hamiltonian(namespace, gr, space, ext_partners, eigensolver, st, hm, max_iter)
A small auxiliary function to perform the update of the Hamiltonian, run the eigensolver,...
Description of the grid, containing information on derivatives, stencil, and symmetries.
abstract class for general interaction partners
iterator for the list of partners
the basic spline datatype
The states_elec_t class contains all electronic wave functions.