Question: Velocity gauge RT-TDDFT - Reproducing Figure
Details
Hello, thank you for sharing this interesting and important project.
I am currently trying to reproduce Figure 3a from the paper "Hybrid Gauge Approach for Accurate Real-Time TDDFT Simulations with Numerical Atomic Orbitals", specifically focusing on the Velocity Gauge (VG) only. I’ve slightly modified the simulation parameters as follows:
Electric field amplitude: Reduced from the paper’s value of 8×10⁻⁵ au to 2×10⁻⁵ au ~ (0.001eV/A)
Time step (dt): Changed from 1.93 as to 2 as
Delta function excitation: Used a Gaussian shape instead
I ran the simulation for 100 steps, but the result does not resemble the behavior shown in Figure 3a of the paper. I am wondering if these parameter changes may significantly impact the outcome, or if I might have missed something important in the setup.
I’m attaching my input files as Si_VG.tar for reference. Could you kindly check if there’s anything wrong in the configuration, or suggest what might be causing the difference?
Thanks in advance for your help!
Have you read FAQ on the online manual http://abacus.deepmodeling.com/en/latest/community/faq.html
- [x] Yes, I have read the FAQ part on online manual.
Task list for Issue attackers (only for developers)
- [ ] Understand the problem or question described by the user.
- [ ] Check if the issue is a known problem or has been addressed in the documentation.
- [ ] Test the issue or problem on a similar system or environment, if possible.
- [ ] Identify the root cause or provide clarification on the user's question.
- [ ] Provide a step-by-step guide, including any necessary resources, to resolve the issue or answer the question.
- [ ] If the issue is related to documentation, update the documentation to prevent future confusion (optional).
- [ ] If the issue is related to code, consider implementing a fix or improvement (optional).
- [ ] Review and incorporate any relevant feedback from users or developers.
- [ ] Ensure the user's issue is resolved or their question is answered and close the ticket.
Thanks for your question, however, we haven't uploaded the new codes in the paper yet, we will do this asap. Please see the release note for the upcoming feature.
One possible reason is that you used a system that is too small. The paper you mentioned used an 8-atom cubic conventional cell with a 16×16×16 k-point mesh. But your simulation only uses a 2-atom silicon primitive cell and a 4×4×4 k-point mesh, which can significantly affect the results.
Thanks for your suggestion. Following your advice, I adjusted the cell size and k-points, but the current is still extremely small. I also tried increasing the number of atoms to 16, but observed similar results.
I think I see the problem. It's a known issue with the workflow — the propagation doesn't really start until the second step. So if you apply the electric field in the first step, the system isn't properly excited. This issue has been fixed, but the related update hasn't been uploaded yet. The fix will be uploaded later once the hybrid gauge code is uploaded. For now, if you want correct results, try using td_gauss_t0 2 or td_gauss_t0 3.
Thank you for the follow-up.
I re-ran the calculation using td_gauss_t0 3.
Upon reviewing the setup, I noticed that the delta-function-like Gaussian pulse did not produce a vector potential on the same scale as in the paper (8×10⁻⁵ a.u.).
To match that scale more closely (approximately 5.36×10⁻⁵ a.u.), I increased the td_gauss_amp from 0.002 to 0.1.
This adjustment led to a larger current, but it is still significantly smaller than the figure shown in the paper (only around 10⁻² a.u., rather than ~1 a.u.).
From inspecting the source code, it seems that the values written to current_total.dat are in atomic units — could you confirm if this is correct?
The vector potential and current output are all in atomic units. But the current shown in the paper is not as large as $$1$$ a.u. — if you check Fig. 3A, their amplitudes are on the order of $$1e^{-6}$$ a.u. To obtain this result, you also need to subtract the current value from the first step. The fact that the current in the first step is not zero is actually another bug, which will also be fixed when the hybrid gauge code is uploaded. Although this bug causes the current to be offset by a constant, our tests show that the oscillations of current are correct, so you can still obtain the desired results.
Finally, if you want to get the same order of magnitude as in the paper, you also need to divide by the volume of the system, which is around $$1000 Bohr^{3}$$ for an 8-atom Si system. You can check the running log for the exact value. This volume normalization will be included in the output in future versions, but for now, you’ll need to do it manually.