Skip to content

[vasp] first draft of VASP dft driver incl example#2

Open
the-hampel wants to merge 15 commits into
TRIQS:unstablefrom
the-hampel:vasp_driver
Open

[vasp] first draft of VASP dft driver incl example#2
the-hampel wants to merge 15 commits into
TRIQS:unstablefrom
the-hampel:vasp_driver

Conversation

@the-hampel

Copy link
Copy Markdown
Member

still in draft mode

@Wentzell

Copy link
Copy Markdown
Member

Let us know once this is ready for review

@the-hampel

the-hampel commented May 26, 2026

Copy link
Copy Markdown
Member Author

Hi @Wentzell ,

in principle it would be but I am waiting for an email response of @harrisonlabollita @GingrasO . The code seems to hit a bug in modest and I wanted them to confirm.

Let me paste the issues here so we can work on it:

I had to change in

  1. I had to remove line 90 in solver_interface/ctseg.py:_prepare_solver() (https://github.com/TRIQS/modest/blob/solver_interfaces/python/triqs_modest/solver_interfaces/ctseg.py)

    S.h_int = h_int

Which is a bug I think? That does not belong there.

  1. Then my example crashes with:
  File "/Users/ahampel/git/triqs/modest/vasp_csc_example/vasp_modest_csc.py", line 87, in <module>
    solver_results = impurity_solve(Delta, Eimp, h_int, **solver_params)
  File "/Users/ahampel/pyvenv/devpy/lib/python3.14/site-packages/triqs_modest/solver_interfaces/ctseg.py", line 155, in solve
    S, solver_params, post_proc_params = _prepare_solver(
                                         ~~~~~~~~~~~~~~~^
        Delta_iw, h_loc0, solver_interface_params
        ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
    )
    ^
  File "/Users/ahampel/pyvenv/devpy/lib/python3.14/site-packages/triqs_modest/solver_interfaces/ctseg.py", line 90, in _prepare_solver
    S.h_loc0_mat = block_matrix_from_op(h_loc0, gf_struct)
                   ~~~~~~~~~~~~~~~~~~~~^^^^^^^^^^^^^^^^^^^
TypeError: Error: no suitable C++ overload found in implementation of function extractors.block_matrix_from_op

block_matrix_from_op(many_body_operator H, gf_struct_t gf_struct, bool ignore_irrelevant) -> block_matrix_t 
 failed with the error : 
  Error message was not convertible to raw string.

This is already quite far into the script just before the solve call. So VASP is invoked correctly, the converter runs correctly (see output log) and there I stopped working.

  1. another question: When running modest with multiple mpi ranks I see that the chemical potential finder is not running mpi parallelized is that correct?

I put my example in the PR. You just have to replace the POTCAR dummy in the PR with the correct one.

Best,
Alex
log.txt

@cnyeh

cnyeh commented May 27, 2026

Copy link
Copy Markdown

Hi @the-hampel,

Point 1 is definitely a bug. I am not yet sure about the cause of point 2, I will take a closer look at it while fixing the bug int point 1.

@the-hampel

Copy link
Copy Markdown
Member Author

Made progress: TRIQS/modest#9 and I pushed some updates to this branch as well. Now I am facing some problem with the complex part of the hybridization I think. There are some changes that I do not quite fully understand.

@the-hampel

the-hampel commented Jun 9, 2026

Copy link
Copy Markdown
Member Author

While running the VASP CSC example (doc/examples/vasp_csc_svo/vasp_modest_csc.py) end-to-end I hit two separate bugs that block a full charge self-consistent loop. Here is a summary from my little helper claude:

1. dftkit — band-count mismatch in band_energy_and_write_charge_update

python/triqs_dftkit/vasp/driver.py:299 dots the full N_k[ik][spin] (shape n_max_bands × n_max_bands) against Hk[ik, idx, :nb, :nb] (sliced to the per-k band count nb = n_bands_per_k[ik, idx]). At any k-point where nb < n_max_bands this raises:

ValueError: shapes (5,5) and (4,4) not aligned: 5 (dim 1) != 4 (dim 0)

N_k is allocated at n_max_bands, so it must also be sliced to :nb (line 298 already does the right thing via np.diag_indices(nb)):

-                band_energy_correction += np.dot(N_k[ik][spin],          Hk[ik, idx, :nb, :nb]).trace().real * bz_weights[ik]
+                band_energy_correction += np.dot(N_k[ik][spin][:nb, :nb], Hk[ik, idx, :nb, :nb]).trace().real * bz_weights[ik]

With this one-line fix the band energy correction computes fine (Band Energy Correction = 0.475 eV in my run).

2. ctseg — real_or_complex throws on a real-but-complex-typed h_loc0

After the dftkit fix, solve_generic(Delta, Eimp, h_int, ...) fails inside S.solve:

Triqs runtime error at triqs/utility/real_or_complex.hpp:49
Logic error : the number is not real, it is complex

Root cause: real_or_complex::_is_real is set by which constructor was used, not by whether the imaginary part is zero — building from std::complex<double> flags it complex even when imag == 0. In triqs_ctseg/work_data.cpp the interaction path is tolerant (uses real(...)):

auto U_full = dict_to_matrix(extract_U_dict2(p.h_int), p.gf_struct);
U           = nda::matrix<double>{real(U_full)};                  // OK

but the h_loc0 → mu path casts implicitly:

mu          = nda::zeros<double>(n_color);
auto h_loc0 = dict_to_matrix(extract_h_dict(p.h_loc0), p.gf_struct);
for (auto const &col : range(n_color)) { mu(col) = -h_loc0(col, col); }   // real_or_complex -> double, throws

Our h_loc0 (Eimp) comes from op_from_block_matrix of complex128 matrices, so every coefficient is complex-flagged even though max|Im| ~ 1e-18. The interaction path silently drops the imaginary part; the h_loc0 path now throws.

Two ways to address it:

  • ctseg (proper, fixes all callers): mirror the U path — mu(col) = -real(h_loc0(col, col)); in work_data.cpp.
  • caller workaround (what I used to keep going): pass a real h_loc0 to solve_generic, e.g. Eimp_real = [h.real.copy() for h in Eimp]. Since impurity levels are Hermitian/real this just drops ~1e-18 noise.

@the-hampel

the-hampel commented Jun 23, 2026

Copy link
Copy Markdown
Member Author

I think this is ready now. See also TRIQS/modest#11 . @harrisonlabollita @GingrasO please have a look. Up for discussion: Should there be an example in doc/examples or only the modest tutorial? I vote for some examples here but this is of course a bit harder to maintain.

@harrisonlabollita

Copy link
Copy Markdown
Collaborator

Hi @the-hampel, I agree this looks good to me! Thank you for this PR!

@Wentzell, can we go ahead and merge this?

band_energy_and_write_charge_update only stored delta_N in the seedname
h5 (dft_update group), which VASP's ADD_GAMMA_FROM_FILE does not read, so
the CSC loop aborted with 'no GAMMA or vaspgamma.h5 file found'. Also write
the correction to vaspgamma.h5 in the layout VASP expects (top-level
band_window + deltaN/{up,down}, one nb x nb block per IBZ k-point), mirroring
triqs_dft_tools SumkDFT.calc_density_correction(dm_type='vasp').

Handles the non-magnetic / spin-averaged case (up == down); SP!=0 and SOC
are not written by this h5 path (SO needs the legacy GAMMA text file).
band_energy_and_write_charge_update now branches on the SO flag in the
seedname h5: for spin-orbit it writes the single spinor-band channel to
deltaN/ud (no spin average); otherwise the spin-averaged collinear blocks
to deltaN/up and deltaN/down. This matches what VASP's ADD_GAMMA_FROM_FILE
reads and mirrors triqs_dft_tools calc_density_correction(dm_type='vasp').

Add test/python/vasp/charge_update with collinear, spin-orbit and IBZ-slicing
cases: the expected correction is built independently in numpy, written to a
reference vaspgamma.h5, and compared to the writer output with triqs h5diff;
the band-energy correction value is checked separately.
Report the global DFT+DMFT iteration, each inner DMFT iteration, and the
point just before the VASP charge update / DFT driver is triggered.
vasp.Driver landed on this branch; link it in the interface table
instead of marking it as in development.
The upstream ctseg fix (real(...) in work_data.cpp mu extraction) means
solve_generic now accepts the complex-flagged Hermitian impurity levels
directly; no need to strip the imaginary noise with h.real.copy().
Skip the final VASP charge update so the run always terminates on a DMFT
step instead of a now-unused DFT charge update, and add a header comment
plus inline comments explaining the setup and the outer/inner loops.
Call mpi.barrier with poll_msec=100 in the vasp, qe and abinit drivers
(and the plovasp sc_dmft loop) so idle ranks sleep instead of busy-
polling every 1 ms while the external DFT code occupies the CPUs.
Remove a leftover vasp.lock (in addition to STOPCAR) at the start of
run_initial_stage: a stale lock would make the driver believe VASP is
already running. Also clarify in run_update_stage that running several
VASP steps per DMFT iteration is left to the caller.
@the-hampel

Copy link
Copy Markdown
Member Author

@Wentzell @harrisonlabollita the more I think about it I come to believe that it is best to remove the doc/example files - they would really be duplicate between here and modest. I think this should live on the modest site no? Even though it might be hard to find for the user. We should put strong cross links between the docs to point to the tutorial I would say. Should I maybe draft this in here, i.e. doc/tutorials contains a section pointing to modest and dft_tools with tutorials and guides?

@harrisonlabollita

Copy link
Copy Markdown
Collaborator

@the-hampel, when I was drafting the modest docs I had a similar thought because the initial step is running the converter. So it felt odd to me that we would point people to modest/dft_tools, then immediately point them to dftkit to really get started.

So are you suggesting that we embed the documentation of dftkit within modest and dft_tools? Meaning the dftkit doc is duplicated within the modest and dft tools repositories?

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment

Labels

None yet

Projects

None yet

Development

Successfully merging this pull request may close these issues.

4 participants