Octopus
xc_ks_inversion.F90
Go to the documentation of this file.
1!! Copyright (C) 2010 H. Appel, N. Helbig
2!!
3!! This program is free software; you can redistribute it and/or modify
4!! it under the terms of the GNU General Public License as published by
5!! the Free Software Foundation; either version 2, or (at your option)
6!! any later version.
7!!
8!! This program is distributed in the hope that it will be useful,
9!! but WITHOUT ANY WARRANTY; without even the implied warranty of
10!! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
11!! GNU General Public License for more details.
12!!
13!! You should have received a copy of the GNU General Public License
14!! along with this program; if not, write to the Free Software
15!! Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
16!! 02110-1301, USA.
17!!
18
19#include "global.h"
20
22 use debug_oct_m
26 use global_oct_m
27 use grid_oct_m
29 use io_oct_m
32 use ions_oct_m
33 use, intrinsic :: iso_fortran_env
36 use lasers_oct_m
38 use mesh_oct_m
42 use parser_oct_m
44 use space_oct_m
49 use unit_oct_m
52 use xc_f03_lib_m
53 use xc_oct_m
54
55 implicit none
56
57 private
58 public :: &
67
69 integer, public, parameter :: &
70 XC_INV_METHOD_TWO_PARTICLE = 1, &
74
76 integer, public, parameter :: &
77 XC_KS_INVERSION_NONE = 1, &
80
82 integer, public, parameter :: &
83 XC_ASYMPTOTICS_NONE = 1, &
85
86 integer, parameter :: &
87 XC_FLAGS_NONE = 0
88
90 private
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
98 ! List with all the external partners
99 ! This will become a list of interactions in the future
100 type(partner_list_t) :: ext_partners
101
102 integer, public :: max_iter
103 end type xc_ks_inversion_t
104
105
106contains
107
108 ! ---------------------------------------------------------
109 subroutine xc_ks_inversion_init(ks_inv, namespace, gr, ions, st, xc, mc, space, kpoints)
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
119
120 class(lasers_t), pointer :: lasers
121
122 push_sub(xc_ks_inversion_init)
123
124 if(gr%parallel_in_domains) then
125 call messages_not_implemented("Kohn-Sham inversion in parallel", namespace=namespace)
126 end if
127
128 call messages_experimental("Kohn-Sham inversion")
129
130 !%Variable InvertKSmethod
131 !%Type integer
132 !%Default iter_godby
133 !%Section Calculation Modes::Invert KS
134 !%Description
135 !% Selects whether the exact two-particle method or the iterative scheme
136 !% is used to invert the density to get the KS potential.
137 !%Option two_particle 1
138 !% Exact two-particle scheme.
139 !%Option iterative 2
140 !% Iterative scheme for <math>v_s</math>.
141 !%Option iter_stella 3
142 !% Iterative scheme for <math>v_s</math> using Stella and Verstraete method.
143 !%Option iter_godby 4
144 !% Iterative scheme for <math>v_s</math> using power method from Rex Godby.
145 !%End
146 call parse_variable(namespace, 'InvertKSmethod', xc_inv_method_iter_godby, ks_inv%method)
147
148 if (ks_inv%method < xc_inv_method_two_particle &
149 .or. ks_inv%method > xc_inv_method_iter_godby) then
150 call messages_input_error(namespace, 'InvertKSmethod')
151 call messages_fatal(1, namespace=namespace)
152 end if
153
154 !%Variable KSInversionLevel
155 !%Type integer
156 !%Default ks_inversion_adiabatic
157 !%Section Calculation Modes::Invert KS
158 !%Description
159 !% At what level <tt>Octopus</tt> shall handle the KS inversion.
160 !%Option ks_inversion_none 1
161 !% Do not compute KS inversion.
162 !%Option ks_inversion_adiabatic 2
163 !% Compute exact adiabatic <math>v_{xc}</math>.
164 !%End
165 call messages_obsolete_variable(namespace, 'KS_Inversion_Level', 'KSInversionLevel')
166 call parse_variable(namespace, 'KSInversionLevel', xc_ks_inversion_adiabatic, ks_inv%level)
167 if (.not. varinfo_valid_option('KSInversionLevel', ks_inv%level)) call messages_input_error(namespace, 'KSInversionLevel')
168
169 !%Variable KSInversionAsymptotics
170 !%Type integer
171 !%Default xc_asymptotics_none
172 !%Section Calculation Modes::Invert KS
173 !%Description
174 !% Asymptotic correction applied to <math>v_{xc}</math>.
175 !%Option xc_asymptotics_none 1
176 !% Do not apply any correction in the asymptotic region.
177 !%Option xc_asymptotics_sc 2
178 !% Applies the soft-Coulomb decay of <math>-1/\sqrt{r^2+1}</math> to <math>v_{xc}</math> in the asymptotic region.
179 !%End
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
182 call messages_not_implemented("KSInversionAsymptotics /= xc_asymptotics_sc for dimensions higher than 1.")
183 end if
185 ! In order to get properly converged states, given a vxc, we might need to run more than
186 ! one eigensolver run.
187 ! This variable is documented in scf.F90
188 call parse_variable(namespace, 'MaximumIter', 50, ks_inv%max_iter)
190 if (ks_inv%level /= xc_ks_inversion_none) then
191 call states_elec_copy(ks_inv%aux_st, st, exclude_wfns = .true.)
192
193 ! initialize auxiliary random wavefunctions
194 call states_elec_allocate_wfns(ks_inv%aux_st, gr)
195 call states_elec_generate_random(ks_inv%aux_st, gr, kpoints)
196
197 ! initialize densities, hamiltonian and eigensolver
198 call states_elec_densities_init(ks_inv%aux_st, gr)
199
200 lasers => lasers_t(namespace)
202 call lasers_generate_potentials(lasers, gr, space, ions%latt)
203
204 if(lasers%no_lasers > 0) then
205 call ks_inv%ext_partners%add(lasers)
206 call lasers_check_symmetries(lasers, kpoints)
207 else
208 deallocate(lasers)
209 end if
210
211 call hamiltonian_elec_init(ks_inv%aux_hm, namespace, space, gr, ions, ks_inv%ext_partners, ks_inv%aux_st, &
212 kohn_sham_dft, xc, mc, kpoints)
213 call eigensolver_init(ks_inv%eigensolver, namespace, gr, ks_inv%aux_st, mc, space)
214 end if
215
216 pop_sub(xc_ks_inversion_init)
217 end subroutine xc_ks_inversion_init
218
219
220 ! ---------------------------------------------------------
221 subroutine xc_ks_inversion_end(ks_inv)
222 type(xc_ks_inversion_t), intent(inout) :: ks_inv
223
224 type(partner_iterator_t) :: iter
225 class(interaction_partner_t), pointer :: partner
226
227 push_sub(xc_ks_inversion_end)
228
229 if (ks_inv%level /= xc_ks_inversion_none) then
230 ! cleanup
231 call eigensolver_end(ks_inv%eigensolver)
232 call hamiltonian_elec_end(ks_inv%aux_hm)
233 call states_elec_end(ks_inv%aux_st)
234
235 call iter%start(ks_inv%ext_partners)
236 do while (iter%has_next())
237 partner => iter%get_next()
238 safe_deallocate_p(partner)
239 end do
240 call ks_inv%ext_partners%empty()
241
242 end if
243
244 pop_sub(xc_ks_inversion_end)
245 end subroutine xc_ks_inversion_end
246
247
248 ! ---------------------------------------------------------
249 subroutine xc_ks_inversion_write_info(ks_inversion, iunit, namespace)
250 type(xc_ks_inversion_t), intent(in) :: ks_inversion
251 integer, optional, intent(in) :: iunit
252 type(namespace_t), optional, intent(in) :: namespace
253
254 if (ks_inversion%level == xc_ks_inversion_none) return
255
257
258 call messages_print_var_option('KSInversionLevel', ks_inversion%level, iunit=iunit, namespace=namespace)
259
261 end subroutine xc_ks_inversion_write_info
262
263
264 ! ---------------------------------------------------------
270 subroutine invertks_2part(ks_inv, target_rho, nspin, aux_hm, gr, st, eigensolver, &
271 namespace, space, ext_partners)
272 type(xc_ks_inversion_t), intent(in) :: ks_inv
273 real(real64), intent(in) :: target_rho(:,:)
274 integer, intent(in) :: nspin
275 type(hamiltonian_elec_t), intent(inout) :: aux_hm
276 type(grid_t), intent(in) :: gr
277 type(states_elec_t), intent(inout) :: st
278 type(eigensolver_t), intent(inout) :: eigensolver
279 type(namespace_t), intent(in) :: namespace
280 class(space_t), intent(in) :: space
281 type(partner_list_t), intent(in) :: ext_partners
282
283 integer :: asym1, asym2, ip, ispin
284 integer :: np
285 real(real64) :: rr, shift
286 real(real64), parameter :: smalldensity = 5e-12_real64
287 real(real64), allocatable :: sqrtrho(:,:), laplace(:,:), vks(:,:), x_shifted(:)
288 type(spline_t) :: spl, lapl_spl
289
290 push_sub(invertks_2part)
291
292 assert(.not. gr%parallel_in_domains)
293
294 ! In periodic systems, the aymptotics needs to done differently.
295 ! In the KLI code, I used a formula from the litterature. The same logic could used here too - NTD
296 assert(space%periodic_dim == 0)
297
298 ! The treatment of the asymptotics is only valid for the 1D case
299 assert(space%dim == 1)
300
301 np = gr%np
302
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))
306
307 if (any(target_rho(:,:) < -m_epsilon)) then
308 write(message(1),*) "Target density has negative points. min value = ", minval(target_rho(:,:))
309 call messages_warning(1, namespace=namespace)
310 end if
311
312 ! Computing the square root of the density
313 !$omp parallel private(ip)
314 do ispin = 1, nspin
315 !$omp do simd
316 do ip = 1, gr%np
317 sqrtrho(ip, ispin) = sqrt(target_rho(ip, ispin))
318 end do
319 !$omp end do simd
320 end do
321 !$omp end parallel
322
323 ! Computing the Laplacian of \sqrt{\rho}
324 if (space%dim == 1) then
325 ! In the 1D case, we fit by a spline and compute the derivative with it
326 do ispin = 1, nspin
327 call spline_init(spl)
328 call spline_init(lapl_spl)
329 call spline_fit(gr%np, gr%x(:, 1), sqrtrho(:, ispin), spl)
330 ! Importantly, here we use a staggered grid to force interpolation
331 ! Without this, the low density region shows large oscilations, as for the case of finite differences.
332 ! To avoid this, we evaluate the derivative of the cubic spline on a different grid
333 ! than the one used in the calculation, fit this result and reevalute (interpolate) it on the actual grid.
334 ! This is needed as a cubic spline is by construction passing by the points on which it is fitted.
335 ! This leads to a much smoother result for Laplacian of the density in the low-density regions.
336
337 call spline_generate_shifted_grid(spl, x_shifted)
338 call spline_der2(spl, lapl_spl, m_zero, grid=x_shifted)
339 safe_deallocate_a(x_shifted)
340
341
342 do ip = 1, gr%np
343 laplace(ip, ispin) = spline_eval(lapl_spl, gr%x(ip, 1))
344 end do
345 call spline_end(spl)
346 call spline_end(lapl_spl)
347 end do
348
349 else
350 do ispin = 1, nspin
351 call dderivatives_lapl(gr%der, sqrtrho(:, ispin), laplace(:, ispin))
352 end do
353 end if
354
355 ! Used to check the asymptotics
356 ! Note that this does not work for open-shell systems where the up and down density are different.
357 ! This could cause problems for spin-polarized H atom for instance
358 asym1 = 0
359 asym2 = 0
360
361 ! Note that the code below does not works for grid parallelization
362 ! It also assumes a special ordering of the points, which which is only guarantied in 1D
363 ! A better approach would be to fix the constant with the eigenvalue (user input option or restarting)
364 do ispin = 1, nspin
365 !avoid division by zero and set parameters for asymptotics
366 !only for 1D potentials at the moment
367 !need to find a way to find all points from where asymptotics should start in 2D and 3D
368 do ip = 1, int(np/2)
369 if (target_rho(ip, ispin) < smalldensity) then
370 vks(ip, ispin) = aux_hm%ep%vpsl(ip) + aux_hm%ks_pot%vhartree(ip)
371 asym1 = ip
372 end if
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)
375 asym2 = np - ip + 1
376 end if
377 end do
378
379 ! The actual expression is \epsilon + 1/2 [\nabla^2 \sqrt{\rho}]/\sqrt{\rho} - v_{ext} - v_H
380 ! We first compute the first part
381 !$omp parallel do private(ip)
382 do ip = asym1+1, asym2-1
383 vks(ip, ispin) = m_half * laplace(ip, ispin) / (sqrtrho(ip, ispin) + m_tiny)
384 end do
385
386 ! And then compute the remaining part
387 !$omp parallel do private(ip)
388 do ip = 1, gr%np
389 aux_hm%ks_pot%vxc(ip, ispin) = vks(ip, ispin) - aux_hm%ep%vpsl(ip) - aux_hm%ks_pot%vhartree(ip)
390 end do
391 end do
392
393 !ensure correct asymptotic behavior, only for 1D potentials at the moment
394 !need to find a way to find all points from where asymptotics should start in 2D and 3D
395 if (ks_inv%asymptotics == xc_asymptotics_sc) then
396 do ispin = 1, nspin
397 !$omp parallel do private(rr, ip)
398 do ip = 1, asym1
399 call mesh_r(gr, ip, rr)
400 aux_hm%ks_pot%vxc(ip, ispin) = -st%qtot/sqrt(rr**2 + m_one)
401 end do
402 !$omp end parallel do
403
404 ! calculate constant shift for correct asymptotics and shift accordingly
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)
407
408 !$omp parallel do simd private(ip)
409 do ip = asym1+1, asym2-1
410 aux_hm%ks_pot%vxc(ip, ispin) = aux_hm%ks_pot%vxc(ip, ispin) - shift
411 end do
412 !$omp end parallel do simd
413
414
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)
417
418 !$omp parallel private(rr, ip)
419 !$omp do simd
420 do ip = 1, asym2-1
421 aux_hm%ks_pot%vxc(ip, ispin) = aux_hm%ks_pot%vxc(ip, ispin) - shift
422 end do
423 !$omp end do simd
424
425 !$omp do
426 do ip = asym2, np
427 call mesh_r(gr, ip, rr)
428 aux_hm%ks_pot%vxc(ip, ispin) = -st%qtot/sqrt(rr**2 + m_one)
429 end do
430 !$omp end do
431 !$omp end parallel
432 end do
433 end if !apply asymptotic correction
434
435 ! Remove a shift, imposing vxc to be continuous
436 if (ks_inv%asymptotics == xc_asymptotics_none) then
437 assert(asym1 > 0)
438 assert(asym2 > 0)
439 do ispin = 1, nspin
440 ! calculate constant shift to make potential continuous
441
442 ! Note that this will fail if the density does not vanish below small_density
443 shift = aux_hm%ks_pot%vxc(asym1+1, ispin)
444 !$omp parallel do simd
445 do ip = asym1+1, asym2-1
446 aux_hm%ks_pot%vxc(ip, ispin) = aux_hm%ks_pot%vxc(ip, ispin) - shift
447 end do
448 !$omp end parallel do simd
449
450 shift = aux_hm%ks_pot%vxc(asym2-1, ispin)
451 !$omp parallel do simd
452 do ip = 1, asym2-1
453 aux_hm%ks_pot%vxc(ip, ispin) = aux_hm%ks_pot%vxc(ip, ispin) - shift
454 end do
455 !$omp end parallel do simd
456 end do
457 end if
458
459 ! Store the final result
460 do ispin = 1, nspin
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)
462 end do
463
464
465 call invertks_update_hamiltonian(namespace, gr, space, ext_partners, eigensolver, st, aux_hm, ks_inv%max_iter)
466
467 safe_deallocate_a(sqrtrho)
468 safe_deallocate_a(laplace)
469 safe_deallocate_a(vks)
470
471 pop_sub(invertks_2part)
472 end subroutine invertks_2part
473
475 subroutine invertks_update_hamiltonian(namespace, gr, space, ext_partners, eigensolver, st, hm, max_iter)
476 type(namespace_t), intent(in) :: namespace
477 type(grid_t), intent(in) :: gr
478 class(space_t), intent(in) :: space
479 type(partner_list_t), intent(in) :: ext_partners
480 type(eigensolver_t), intent(inout) :: eigensolver
481 type(states_elec_t), intent(inout) :: st
482 type(hamiltonian_elec_t), intent(inout) :: hm
483 integer, intent(in) :: max_iter
484
485 logical :: converged
486 integer :: iter, nst_conv
487
489
490 call hm%update(gr, namespace, space, ext_partners)
491
492 !Poor-man estimate of the number of occupied states. Should be improved
493 nst_conv = floor(st%qtot/st%smear%el_per_state)
494
495 do iter = 1, max_iter
496 call eigensolver%run(namespace, gr, st, hm, 1, converged, nst_conv)
497
498 !TODO: Should we allow for updating the occupations? This might be necessary for open-shell or metallic systems
499 ! call states_elec_fermi(st, namespace, gr)
500
501 call write_iter()
502
503 if (converged) exit
504 end do
505
506 if (converged) then
507 message(1) = 'All states converged.'
508 call messages_info(1, namespace=namespace)
509 else
510 message(1) = 'Some of the states are not fully converged!'
511 call messages_info(1, namespace=namespace)
512 end if
513
514 write(message(1),'(a, e17.6)') 'Criterion = ', eigensolver%tolerance
515 call messages_info(1, namespace=namespace)
516 call messages_new_line()
517
518 ! Get the new guess density
519 call density_calc(st, gr, st%rho)
520
522
523 contains
524 ! ---------------------------------------------------------
525 subroutine write_iter()
526 character(len=50) :: str
527
529
530 write(str, '(a,i5)') 'Kohn-Sham inversion eigensolver iteration #', iter
531 call messages_print_with_emphasis(msg=trim(str), namespace=namespace)
532
533 write(message(1),'(a,i6,a,i6)') 'Converged states: ', minval(eigensolver%converged(1:st%nik))
534 call messages_info(1, namespace=namespace)
535
536 call states_elec_write_eigenvalues(st%nst, st, space, hm%kpoints, &
537 eigensolver%diff, compact = .true., namespace=namespace)
538
539 call messages_print_with_emphasis(namespace=namespace)
540
542 end subroutine write_iter
543 end subroutine invertks_update_hamiltonian
544
545 ! ---------------------------------------------------------
550 subroutine invertks_iter(ks_inv, target_rho, namespace, space, ext_partners, nspin, aux_hm, gr, &
551 st, eigensolver)
552 type(xc_ks_inversion_t), intent(in) :: ks_inv
553 real(real64), intent(in) :: target_rho(:,:)
554 type(grid_t), intent(in) :: gr
555 type(namespace_t), intent(in) :: namespace
556 class(space_t), intent(in) :: space
557 type(partner_list_t), intent(in) :: ext_partners
558 type(states_elec_t), intent(inout) :: st
559 type(hamiltonian_elec_t), intent(inout) :: aux_hm
560 type(eigensolver_t), intent(inout) :: eigensolver
561 integer, intent(in) :: nspin
562
563 integer :: ip, ispin, ierr, asym1, asym2
564 integer :: iunit, verbosity, counter, np
565 integer :: max_iter
566 integer :: imax
567 real(real64) :: rr, shift
568 real(real64) :: alpha, beta
569 real(real64) :: mu, npower, npower_in ! these constants are from Rex Godbys scheme
570 real(real64) :: convdensity, diffdensity
571 real(real64), allocatable :: vhxc(:,:)
572 real(real64), parameter :: small_density = 5e-12_real64
573
574 character(len=256) :: fname
575
576 push_sub(invertks_iter)
577
578 ! The treatment of the asymptotics is only valid for the 1D case
579 assert(space%dim == 1)
580
581 np = gr%np
582
583 !%Variable InvertKSConvAbsDens
584 !%Type float
585 !%Default 1e-5
586 !%Section Calculation Modes::Invert KS
587 !%Description
588 !% Absolute difference between the calculated and the target density in the KS
589 !% inversion. Has to be larger than the convergence of the density in the SCF run.
590 !%End
591 call parse_variable(namespace, 'InvertKSConvAbsDens', 1e-5_real64, convdensity)
592
593 !%Variable InvertKSStellaBeta
594 !%Type float
595 !%Default 1.0
596 !%Section Calculation Modes::Invert KS
597 !%Description
598 !% residual term in Stella iterative scheme to avoid 0 denominators
599 !%End
600 call parse_variable(namespace, 'InvertKSStellaBeta', 1e-6_real64, beta)
601
602 !%Variable InvertKSStellaAlpha
603 !%Type float
604 !%Default 0.05
605 !%Section Calculation Modes::Invert KS
606 !%Description
607 !% prefactor term in iterative scheme from L Stella
608 !%End
609 call parse_variable(namespace, 'InvertKSStellaAlpha', 0.25_real64, alpha)
610
611 !%Variable InvertKSGodbyMu
612 !%Type float
613 !%Default 1.0
614 !%Section Calculation Modes::Invert KS
615 !%Description
616 !% prefactor for iterative KS inversion convergence scheme from Godby based on van Leeuwen scheme
617 !%End
618 call parse_variable(namespace, 'InvertKSGodbyMu', m_one, mu)
619
620 !%Variable InvertKSGodbyPower
621 !%Type float
622 !%Default 0.05
623 !%Section Calculation Modes::Invert KS
624 !%Description
625 !% power to which density is elevated for iterative KS inversion convergence
626 !% scheme from Godby based on van Leeuwen scheme
627 !%End
628 call parse_variable(namespace, 'InvertKSGodbyPower', 0.05_real64, npower_in)
629 npower = npower_in
630
631 !%Variable InvertKSVerbosity
632 !%Type integer
633 !%Default 0
634 !%Section Calculation Modes::Invert KS
635 !%Description
636 !% Selects what is output during the calculation of the KS potential.
637 !%Option 0
638 !% Only outputs the converged density and KS potential.
639 !%Option 1
640 !% Same as 0 but outputs the maximum difference to the target density in each
641 !% iteration in addition.
642 !%Option 2
643 !% Same as 1 but outputs the density and the KS potential in each iteration in
644 !% addition.
645 !%End
646 call parse_variable(namespace, 'InvertKSVerbosity', 0, verbosity)
647 if (verbosity < 0 .or. verbosity > 2) then
648 call messages_input_error(namespace, 'InvertKSVerbosity')
649 call messages_fatal(1, namespace=namespace)
650 end if
651
652 !%Variable InvertKSMaxIter
653 !%Type integer
654 !%Default 200
655 !%Section Calculation Modes::Invert KS
656 !%Description
657 !% Selects how many iterations of inversion will be done in the iterative scheme
658 !%End
659 call parse_variable(namespace, 'InvertKSMaxIter', 200, max_iter)
660
661 safe_allocate(vhxc(1:np, 1:nspin))
662 call lalg_copy(np, nspin, aux_hm%ks_pot%vhxc, vhxc)
663
664 if (verbosity == 1 .or. verbosity == 2) then
665 iunit = io_open('InvertKSconvergence', namespace, action = 'write')
666 end if
667
668 counter = 0
669 imax = 0
670
671 diffdensity = m_zero
672 do ispin = 1, nspin
673 do ip = 1, np
674 if (abs(st%rho(ip, ispin)-target_rho(ip, ispin)) > diffdensity) then
675 diffdensity = abs(st%rho(ip, ispin)-target_rho(ip, ispin))
676 imax = ispin
677 end if
678 end do
679 end do
680
681
682 do while(diffdensity > convdensity .and. counter < max_iter)
683
684 counter = counter + 1
685
686 if (verbosity == 2) then
687 write(fname,'(i6.6)') counter
688 call dio_function_output(io_function_fill_how("AxisX"), ".", "vhxc"//fname, namespace, space, &
689 gr, aux_hm%ks_pot%vhxc(:,1), units_out%energy, ierr)
690 call dio_function_output(io_function_fill_how("AxisX"), ".", "rho"//fname, namespace, space, &
691 gr, st%rho(:,1), units_out%length**(-space%dim), ierr)
692 end if
693
694 ! Update the Hamiltonian
695 call invertks_update_hamiltonian(namespace, gr, space, ext_partners, eigensolver, st, aux_hm, ks_inv%max_iter)
696
697 ! Iterative inversion with fixed parameters in Stella Verstraete method
698 select case(ks_inv%method)
700 !$omp parallel private(ip)
701 do ispin = 1, nspin
702 !$omp do simd
703 do ip = 1, np
704 vhxc(ip, ispin) = vhxc(ip, ispin) &
705 + ((st%rho(ip, ispin) - target_rho(ip, ispin))/(target_rho(ip, ispin) + beta))*alpha
706 end do
707 !$omp end do simd
708 end do
709 !$omp end parallel
710
711 ! adaptative iterative method, with update of alpha and beta coefficients
712 ! based on residual in density
714 beta = diffdensity*1e-3_real64 !parameter to avoid numerical problems due to small denominator
715
716 ! proposition to increase convergence speed progressively
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
721 call messages_info(1, namespace=namespace)
722
723 !$omp parallel private(ip)
724 do ispin = 1, nspin
725 !$omp do simd
726 do ip = 1, np
727 vhxc(ip, ispin) = vhxc(ip, ispin) &
728 + ((st%rho(ip, ispin) - target_rho(ip, ispin))/(target_rho(ip, ispin) + beta))*alpha
729 end do
730 !$omp end do simd
731 end do
732 !$omp end parallel
733
735! ! below 1.e-3 start reducing power down to 0.01
736! if (diffdensity < 0.001_real64) then
737! npower = min(npower_in, diffdensity*50.0_real64)
738! end if
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
742 call messages_info(1, namespace=namespace)
743
744 !$omp parallel private(ip)
745 do ispin = 1, nspin
746 !$omp do simd
747 do ip = 1, np
748 vhxc(ip, ispin) = vhxc(ip, ispin) &
749 + (st%rho(ip, ispin)**npower - target_rho(ip, ispin)**npower)*mu
750 end do
751 !$omp end do simd
752 end do
753 !$omp end parallel
754
755 case default
756 assert(.false.)
757 end select
758
759 diffdensity = m_zero
760 !diffdensity = maxval(abs(st%rho(1:np,1:nspin)-target_rho(1:np,1:nspin)))
761
762 ! Note: The code below requires a reduction over the grid.
763 ! We now have routines for this, see mesh_minmaxloc - NTD
764 assert(.not. gr%parallel_in_domains)
765
766 do ispin = 1, nspin
767 do ip = 1, np
768 if (abs(st%rho(ip, ispin)-target_rho(ip, ispin)) > diffdensity) then
769 diffdensity = abs(st%rho(ip, ispin)-target_rho(ip, ispin))
770 imax = ispin
771 end if
772 end do
773 end do
774
775 if (verbosity == 1 .or. verbosity == 2) then
776 write(iunit,'(i6.6)', advance = 'no') counter
777 write(iunit,'(es18.10)') diffdensity
778#ifdef HAVE_FLUSH
779 call flush(iunit)
780#endif
781 end if
782
783 ! Store the final result
784 call lalg_copy(np, nspin, vhxc, aux_hm%ks_pot%vhxc)
785 do ispin = 1, nspin
786 aux_hm%ks_pot%vxc(1:np, ispin) = vhxc(1:np, ispin) - aux_hm%ks_pot%vhartree(1:np)
787 end do
788
789 end do ! end while statement on convergence
790
791 !ensure correct asymptotic behavior, only for 1D potentials at the moment
792 !need to find a way to find all points from where asymptotics should start in 2 and 3D
793 if (ks_inv%asymptotics == xc_asymptotics_sc) then
794 !Asympotics is only for vxc
795 call lalg_copy(np, nspin, aux_hm%ks_pot%vxc, vhxc)
796
797 do ispin = 1, nspin
798 ! In case the density is never smaller than 'small_density', we take the edge points
799 asym1 = 1
800 asym2 = np
801
802 do ip = 1, int(np/2)
803 if (target_rho(ip, ispin) < small_density) then
804 call mesh_r(gr, ip, rr)
805 vhxc(ip, ispin) = -st%qtot/sqrt(rr**2 + m_one)
806 asym1 = ip
807 end if
808 if (target_rho(np-ip+1, ispin) < small_density) then
809 asym2 = np - ip + 1
810 end if
811 end do
812
813 ! calculate constant shift for correct asymptotics and shift accordingly
814 call mesh_r(gr, asym1+1, rr)
815 shift = vhxc(asym1+1, ispin) + st%qtot/sqrt(rr**2 + m_one)
816
817 !$omp parallel do simd
818 do ip = asym1+1, asym2-1
819 vhxc(ip, ispin) = vhxc(ip, ispin) - shift
820 end do
821 !$omp end parallel do simd
822
823 call mesh_r(gr, asym2-1, rr)
824 shift = vhxc(asym2-1, ispin) + st%qtot/sqrt(rr**2 + m_one)
825
826 !$omp parallel private(rr)
827 !$omp do simd
828 do ip = 1, asym2-1
829 vhxc(ip, ispin) = vhxc(ip, ispin) - shift
830 end do
831 !$omp end do simd
832
833 !$omp do
834 do ip = asym2, np
835 call mesh_r(gr, ip, rr)
836 vhxc(ip, ispin) = -st%qtot/sqrt(rr**2 + m_one)
837 end do
838 !$omp end do
839 !$omp end parallel
840 end do
841
842 ! See the above comment, vhxc contains vxc here
843 call lalg_copy(np, nspin, vhxc, aux_hm%ks_pot%vxc)
844 do ispin = 1, nspin
845 aux_hm%ks_pot%vhxc(1:np, ispin) = vhxc(1:np, ispin) + aux_hm%ks_pot%vhartree(1:np)
846 end do
847 end if
848
849 !TODO: check that all arrays needed by hamiltonian update are sync`d in MPI fashion
850
851 !calculate final density
852 call invertks_update_hamiltonian(namespace, gr, space, ext_partners, eigensolver, st, aux_hm, ks_inv%max_iter)
853
854 write(message(1),'(a,I8)') "Iterative KS inversion, iterations needed:", counter
855 call messages_info(1, namespace=namespace)
856
857 call io_close(iunit)
858
859 safe_deallocate_a(vhxc)
860
861 pop_sub(invertks_iter)
862 end subroutine invertks_iter
863
864 ! ---------------------------------------------------------
865 subroutine xc_ks_inversion_calc(ks_inversion, namespace, space, gr, hm, ext_partners, st, vxc, time)
866 type(xc_ks_inversion_t), intent(inout) :: ks_inversion
867 type(namespace_t), intent(in) :: namespace
868 class(space_t), intent(in) :: space
869 type(grid_t), intent(in) :: gr
870 type(hamiltonian_elec_t), intent(in) :: hm
871 type(partner_list_t), intent(in) :: ext_partners
872 type(states_elec_t), intent(inout) :: st
873 real(real64), intent(inout) :: vxc(:,:)
874 real(real64), optional, intent(in) :: time
875
876 integer :: ispin
877
878 if (ks_inversion%level == xc_ks_inversion_none) return
879
880 push_sub(x(xc_ks_inversion_calc))
881
882 if (present(time) .and. time > m_zero) then
883 write(message(1),'(A,F18.12)') 'xc_ks_inversion_calc - time:', time
884 call messages_info(1, namespace=namespace)
885 end if
886
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
891
892 ks_inversion%aux_hm%ks_pot%vhartree = hm%ks_pot%vhartree
893
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)
896
897 else
898! no restart data available, start with vhxc = vh, which we know from exact input rho
899 do ispin = 1, st%d%nspin
900 ks_inversion%aux_hm%ks_pot%vxc(1:gr%np, ispin) = m_zero !hm%ep%vpsl(:)
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)
903 end do
904! TODO: restart data found. Use first KS orbital to invert equation and get starting vhxc
905! call invertks_2part(ks_inversion%aux_st%rho, st%d%nspin, ks_inversion%aux_hm, gr, &
906! ks_inversion%aux_st, ks_inversion%eigensolver, namespace, ks_inversion%asymp)
907
908 end if
909 !ks_inversion%aux_hm%ep%vpsl(:) = M_ZERO ! hm%ep%vpsl(:)
910 call lalg_copy(gr%np, hm%ep%vpsl, ks_inversion%aux_hm%ep%vpsl)
911
912 ! compute ks inversion, vhxc contains total KS potential
913 ! Update the Hamiltonian
914 call invertks_update_hamiltonian(namespace, gr, space, ext_partners, &
915 ks_inversion%eigensolver, ks_inversion%aux_st, ks_inversion%aux_hm, ks_inversion%max_iter)
916
917
918 ! these 2 routines need to be cleaned - they are not consistent in updating
919 ! the hamiltonian, states, etc...
920 select case (ks_inversion%method)
921 ! adiabatic ks inversion
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)
925
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)
929 end select
930
931 ! subtract Hartree potential
932 ! ATTENTION: subtracts true external potential not adiabatic one
933
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))
936 end do
937
938 vxc = ks_inversion%aux_hm%ks_pot%vxc
939
940 ! save vhxc for next step if we are running td
941 if (present(time)) then
942 ks_inversion%vhxc_previous_step = ks_inversion%aux_hm%ks_pot%vhxc
943 end if
944
945 pop_sub(x(xc_ks_inversion_calc))
946
947 end subroutine xc_ks_inversion_calc
948
949
950end module xc_ks_inversion_oct_m
951
952!! Local Variables:
953!! mode: f90
954!! coding: utf-8
955!! End:
constant times a vector plus a vector
Definition: lalg_basic.F90:171
Copies a vector x, to a vector y.
Definition: lalg_basic.F90:186
Some operations may be done for one spline-function, or for an array of them.
Definition: splines.F90:179
double floor(double __x) __attribute__((__nothrow__
This module implements a calculator for the density and defines related functions.
Definition: density.F90:120
subroutine, public density_calc(st, gr, density, istin)
Computes the density from the orbitals in st.
Definition: density.F90:609
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
Definition: global.F90:188
real(real64), parameter, public m_tiny
Definition: global.F90:205
real(real64), parameter, public m_epsilon
Definition: global.F90:204
real(real64), parameter, public m_half
Definition: global.F90:194
real(real64), parameter, public m_one
Definition: global.F90:189
This module implements the underlying real-space grid.
Definition: grid.F90:117
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...
Definition: io.F90:114
subroutine, public io_close(iunit, grp)
Definition: io.F90:418
integer function, public io_open(file, namespace, action, status, form, position, die, recl, grp)
Definition: io.F90:352
A module to handle KS potential, without the external potential.
integer, parameter, public kohn_sham_dft
subroutine, public lasers_check_symmetries(this, kpoints)
Definition: lasers.F90:555
subroutine, public lasers_parse_external_fields(this)
Definition: lasers.F90:244
subroutine, public lasers_generate_potentials(this, mesh, space, latt)
Definition: lasers.F90:444
This module defines the meshes, which are used in Octopus.
Definition: mesh.F90:118
pure subroutine, public mesh_r(mesh, ip, rr, origin, coords)
return the distance to the origin for a given grid point
Definition: mesh.F90:336
subroutine, public messages_print_with_emphasis(msg, iunit, namespace)
Definition: messages.F90:920
subroutine, public messages_not_implemented(feature, namespace)
Definition: messages.F90:1113
character(len=512), private msg
Definition: messages.F90:165
subroutine, public messages_warning(no_lines, all_nodes, namespace)
Definition: messages.F90:537
subroutine, public messages_obsolete_variable(namespace, name, rep)
Definition: messages.F90:1045
subroutine, public messages_new_line()
Definition: messages.F90:1134
character(len=256), dimension(max_lines), public message
to be output by fatal, warning
Definition: messages.F90:160
subroutine, public messages_fatal(no_lines, only_root_writes, namespace)
Definition: messages.F90:414
subroutine, public messages_input_error(namespace, var, details, row, column)
Definition: messages.F90:713
subroutine, public messages_experimental(name, namespace)
Definition: messages.F90:1085
subroutine, public messages_info(no_lines, iunit, debug_only, stress, all_nodes, namespace)
Definition: messages.F90:616
This module handles the communicators for the various parallelization strategies.
Definition: multicomm.F90:145
subroutine, public spline_fit(nrc, rofi, ffit, spl, threshold)
Definition: splines.F90:413
subroutine, public spline_der2(spl, dspl, threshold, grid)
Returns a spline that contains the second derivative of the original spline.
Definition: splines.F90:959
real(real64) function, public spline_eval(spl, x)
Definition: splines.F90:441
subroutine, public spline_generate_shifted_grid(this, x_shifted)
Definition: splines.F90:1171
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.
Definition: unit.F90:132
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.
Definition: grid.F90:168
abstract class for general interaction partners
the basic spline datatype
Definition: splines.F90:156
The states_elec_t class contains all electronic wave functions.
int true(void)
subroutine write_iter()