Real Time Magnus

[pwborthwick]

Hi Josh, Sorry me again! In your /examples/real-time.py you have

      plt.plot(rt.time,np.asarray(rt.shape)*rt.field,label='Applied field')

however, at some point you've removed the

      import numpy as np

I have some questions, for example the pulse is set to None so the shape is zero always. This means that the dipole is never added to the field and the energy does not change which you can see by plotting rt.Energy. I'm wondering if the, rich in detail, dipole curves are an artifact of the calculation? The dipole is ~1e-15 and given that the matrix exponent is a Pade approximation...An estimate of the local error can be calculated as the norn of the difference between third and fourth order Magnus Numerical integrators based on the Magnus expansion for nonlinear dynamical systems, Hajiketabi & Casas . I notice that if you apply a pulse (your Gaussian one) then you get a smoothly oscillating dipole where both Magnus2 and Magnus4 coincide, I'm still trying to puzzle out why? I think this code is the only one on github using the Magnus expansion in this way and so is important and worth expanding (time permitting). Best Peter.

[pwborthwick]

I think the problem is that hydrogen molecule doesn't have a dipole so you're plotting numerical error. I've tried HF which does have a dipole and I get better results. Although you have to remove the plot of the applied field as that collapses the scale so detail in dipole is lost. Maybe change the example? With HF without a field - should the variation in dipole match the vibrational frequency we might calculate with BOMD? It seems it ought to so we could check the quantitative results possibly?

[jjgoings]

Hi Peter,

You are right, we need the import numpy as np!

I agree that is a poor example for RT-TDHF. It's a trivial example, but the important point is that energy is conserved, which is why it's a test in tests/test011.py. I should have had energy plotted, not dipole. 10^-15 is in the noise, and I'm surprised it had any structure at all! I think machine precision is around 10^-14.

Anyway, back to energy conservation: not all TDSE integrators are energy conserving, but those based off Magnus are, and are known as symplectic integrators. In other words, the error in the energy is bounded, unlike Runge-Kutta-4 which although technically more accurate, will show long-time energy drift. One of the first checks when verifying the implementation of an integrator for RT-TDSCF is energy conservation in the absence of external fields.

I think I mangled this example after trying to do too many experimental things with RT-TDHF a few years ago. I'll update soon with a better RT example and make sure any other bugs are worked out. Thanks!

[pwborthwick]

Hi Josh, I'm sure a new example when you have time would be good as many people use your program as a source of ideas and as a code reference. I think I mentioned symplectic integrators when we were discussing the velocity-Verlet earlier, although the Magnus is new to me having recently discovered it through your blog. I wonder how non-symplectic ones like the 4th-order Runge Kutta would compare over the range of this example? Of course the Magnus also preserves the unitarity unlike the Dyson series. I've run your example for HF (bond length 0.92A) and the energy plots as constant. The dipole oscillates gradually shifting its range... all in absence of external field. I've also run a calculation and reversed the time mid-way and that was fine. But I don't see why the dipole shifts as it does in the example above?

I noticed that you had intended, in addition to dipole, to plot 'angular'. I assume this was going to use the angular integral? I'm guessing this is the angular momentum? I ran a comparison between mmd and psi4 (using water) and got different matrices. I used mints.ao_angular_momentum() with sto-3g basis, scf_type = direct and puream = False. Energy and dipole agree with mmd, am I comparing like with like? I had a very quick look at code and I can't see anything obviously wrong. It's just a combination of a dipole and a nabla integral. If this is a real issue let me know and I'll elaborate with some numbers from the calculations. Thanks for all your efforts @jjgoings

[pwborthwick]

Hi @jjgoings , I have revisited this issue to add some observations and ask for your comments. Firstly, the seeming downward movement of the hydrogen fluoride dipole (in absence of external fields) is due to a continuous drop in the molecular energy during the calculation. The energy is constant upto the 11th decimal place but that small energy loss causes the dipole range to fall accordingly. So it's not something that happens in the real world! As for the dipole oscillation itself. HF has a real dipole and a BOMD calculation using mmd shows a frequency of ~4000 cm^-1 which is in agreement with the experimental value of 4461 cm^-1. However, your TDHF (I can't see any conversion factors so I assume working is in atomic time units) gives a frequency a factor 10 greater. The oscillation's origin is in the density matrix (you can see this by printing out Tr(ρ²) which shows the oscillation (again in 11th decimal place)). So a) is the dipole behaviour real, in which case what does the period of the oscillation represent? Or b) is it just noise, we've seen with H₂ noise can produce a seeming oscillation. I can't find any reference calculation for Magnus propogating a molecular dipole, but https://nwchemgit.github.io/RT-TDDFT.html#Hints_and_Tricks has some results that may be reproducible if I implement a 'kick' - I'll try to do it and see what happens. But I still don't know if without an applied field the graph is 'real'? Best Peter

[jjgoings]

Okay, I (think) I have a good minimal example for RT on H2/3-21G with commit 37514399a98e43

With a narrow Gaussian envelope (essentially a delta-kick), we get the non-trivial z-dipole propagation like so:

And working up the spectra with the utilities in mmd/utils/spectrum.py; e.g. Pade transform variant of FFT:

Note that Gaussian16 with LR-TDHF gives peaks with non-zero oscillator strength at

 Excitation energies and oscillator strengths:

 Excited State   1:      Singlet-SGU   15.4997 eV   79.99 nm  f=0.6685  <S**2>=0.000
       1 ->  2         0.70960
 This state for optimization and/or second-order correction.
 Total Energy, E(TD-HF/TD-DFT) = -0.553338323725
 Copying the excited state density for this state as the 1-particle RhoCI density.

 Excited State   2:      Singlet-SGG   31.9771 eV   38.77 nm  f=0.0000  <S**2>=0.000
       1 ->  3         0.70794

 Excited State   3:      Singlet-SGU   46.4065 eV   26.72 nm  f=0.0381  <S**2>=0.000
       1 ->  4         0.70804

So we have the right number of peaks at the right energies and at the right oscillator strengths. If you integrate over each peak you should get the oscillator strength f.

To get the right sign I had to add the static electric field in the addField routine, which I was always under the impression was supposed to be subtracted. So I may have missed a sign somewhere and want to have this choice documented in case I need to fix it later! In practice, this is not such a big deal for absorption cross sections. Usually the absolute value of the y-axis (intensity) is re-scaled. This may cause issues with sign-sensitive spectroscopies like electronic circular dichroism, so I need to think about this more carefully.

[jjgoings]

@pwborthwick :

In regards to your observation, the origin of the oscillating dipole moment between BOMD and RT-TDHF is a bit subtle. In BOMD, the electrons are assumed to move adiabatically in the sense that we converge a new SCF at each time step. So the dipole change is due to nuclear vibrations on the ground state. Getting the ~4000 cm-1 for HF is a good sign that we are getting the nuclear gradients correct.

In contrast, with RT-TDHF the nuclei are fixed and it is the electronic density that oscillates around the "fixed points" of the nuclear potential. You correctly observe this. These oscillations should be much more rapid than the nuclear vibrations and give rise to the electronic absorption spectrum (at higher energies).

Put another way, a time-independent optimization + frequency calculation is analogous to a time-dependent BOMD calculation, whereas a "time-independent" linear response LR-TDHF is analogous to an explicitly time-dependent RT-TDHF calculation. The purposes might vary, but in general the spectroscopic predictions should agree between methods at a given level of theory.

I think that the dipole oscillations around the 11th decimal place are effectively in the noise due to machine precision and convergence criteria, but the oscillations may have some physical connection. The connection is due to the fact that the time derivative of the density matrix is given by the commutator [F,P] but the convergence criteria for the SCF is also [F,P]. So if the SCF is not so tightly converged that [F,P] is strictly zero, there may be some residual oscillations in the time evolution, which it appears you can plot and view. This also make sense in light that you see some energy loss -- the SCF is not sufficiently converged if you want to look at this high of precision.

(Side note that is also with mentioning: It is also possible to converge the SCF equations with the RT-TDHF equations by propagating in imaginary time. The change from real to imaginary time converts the problem from a dynamics to a statics problem. The connections between time evolution of the SCF and the convergence of the SCF are quite intricate!)

[pwborthwick]

Hi @jjgoings , Wow, thank you for your reply which as usual contains a wealth of information. Sorry it took me so long to get back to TDHF but you will remember we got into angular momemtum and gauge centers which lead to coding nabla, quadrupole and electric field and potential integrals which in turn lead to a foray into restrained ESP. I also decided that although I understood the Magnus algorithm there were subtilties in the implementation that I was missing and needed to think about. TDHF I see, now you've explained, that the overweight hadrons will not move as fast as the 'gym-going' svelte leptons so the difference between BOMD and TDHF oscillations are explained. The TDHF dipole (no field) oscillation is ghostly as it lies tantalisingly near the noise - sometimes it seems real, sometimes not. I suspect it may be real but with a period of tenths of femtoseconds you'd need a sharp eye and a good watch to measure it. Also if you apply a field the resulting oscillation looks to be suspiciously of the same period as the 'ghost', so maybe it is something fundemental. I'd like to think about this some more and look at other cases. The example you now have is great, I love the spectrum. I find TDHF interesting and will go back to my coding it to explore further. I love the way you add little remarks like 'imaginary time' and 'negative energy' (your relativistic HF article), they get the mind working overtime. Incidently a negative energy would make Magnum propogate into the past rather than the future.

Gauge I decided on nuclear charge center as default in my code and added flexibility to select other options. You will want to evaluate a ESP at an arbitrary point anyway. This won't be the same as the center of charge although the electron distribution is centered on the nuclei. But if there are eg lone-pairs then the center of nuclear charge may be away from the center of charge. But finding the center of charge is not trivial. With the Born-Oppenheumer approximation we have a fixed set of nuclei but only an averaged electron distribution. There is a paper https://www.sciencedirect.com/science/article/abs/pii/0009261471800817 and the abstract says 'The positions of the effective centers of electronic charge have been calculated for the atoms in a group of homonuclear diatomic molecules. For all of the molecules studied except Li2, these centers of charge were found to be outside the internuclear region.' Electric Field I thought the positive direction of a static electric field was away from positive charge and towards negative charge. The dipole vector points from negative to positive! Since we are dealing with dipoles does this explain the sign? Integrals As I said I've been coding a few integrals and I notice that with the derivatives there are savings to be made using translational invariance. Eg with water the form of the sto-3g gradient overlap matrix is	Os	Os	Op	Op	Op	H1	H2
0	0	0	0	0	x	x
0	0	0	0	0	x	x
0	0	0	0	0	x	x
0	0	0	0	0	x	x
0	0	0	0	0	x	x
x	x	x	x	x	0	x
x	x	x	x	x	x	0

This structure occurs because if the orbitals (μ,ν) are on the same center then S_μν will be zero by translational invariance. In fact, in general, S^ax_μν = -S^bx_μν for centers a and b and similarly for y,z. Should be true for Kinetic integrals too but not Coulomb. Might save a bit of time in calculating forces?

Thanks again, I hope this isn't taking up too much of your time? I'm retired and in lock-down because of Corvid so not much else to do really but you have a job so I appreciate your time spent. Best Peter

[pwborthwick]

This is interesting... HF with same pulse as in your test example sto-3g... and now nothing changed except basis set -> 3-21g... Clearly doesn't happen with H₂, so because HF has an intrinsic dipole ? I think avoid sto-3g for these computations I get strange results for spectra with water with sto-3g as well. Peter

[pwborthwick]

Hi @jjgoings In the magnus2 and 4 routines you are getting the dipole from scf.computeDipole(). Doesn't this return the dipole in Debyes - the last statement is

        self.mu *= 2.541765

These 'Debye' values go into dipole list which is then used directly as mt.dipole in plot with axis label

        ax1.set_ylabel('z-dipole / au',color='tab:blue')

Shouldn't that be 'z-dipole / Debye' ? It might explain why I'm getting values 2.5 times larger than yours in a version of HDHF I'm trying to write...but it could be I'm missing something. Best Peter

[jjgoings]

@pwborthwick Yeah, good catch. That's an easy fix -- already done.

[pwborthwick]

@jjgoings Reproduced the nwchem example I referenced earlier for water in 6-31g basis with mmd. Has kick of 0.0001 at t=0 and result looks good.

If you use a Gaussian centered away from t=0 (t=2) eg exp(-10(t-2)(t-2)) then plotting the resultant dipole you get

Peter

[jjgoings]

@pwborthwick Not sure if this was a question or not, but in the second image it looks like you are plotting the Fourier transform of the dipole, not the dipole itself.

An idealized (decaying) oscillating dipole after a narrow Gaussian electric field "kick" at time = 0 is

where theta is the Heaviside step function. When the oscillating dipole is Fourier transformed, and you just look at the imaginary component, you get a Lorentzian (Cauchy) function centered at the period,

where "a" is related to the FWHM of the Lorentizian. If you perturb the dipole at a later time, but still Fourier Transform over t = [0,+inf) you will get phase issues and the real and imaginary components "mix". That is, you take the Fourier Transform of

(I'm not sure of the analytic result here, but I know it will be a phase mixing of the real and imag components of the t = 0 case.) So in this respect you plot looks OK, though it looks unlike how the experimental spectra should appear. Plotting the abs of the Fourier Transform can eliminate the phase, although the line shape and intensity are for a power spectrum, not an absorption cross section. The energies will be at the correct locations, however.

[pwborthwick]

Hi,   Sorry to clarify the second plot uses same parameters as first plot, but instead of a kick at t=0 has a narrow gaussian at t=2 as shown. It's not what I expected, I think similar things happen for your H2 example with a pulse at t>0 - I'll check tomorrow. My feeling was something wasn't right but I need to spend more time experimenting. The data plotted were the same variables in both diagrams. I'll go back to H2 example tomorrow and try to sort out what's happening.      Peter 

[pwborthwick]

These are the plots you get by using respectively in genSpectra,

      numerator = np.imag(fw)  : np.real(fw)  : np.abs(fw)

with a Gaussian pulse at t=2

compare with kick (narrow Gaussian) at t=0 as you say the t=2 imaginary plot is a sort of mix of t=0 real and imaginary (below is a simple Im + Re (t=0) plot)...

I guess all the plots are correct and you just have to chose which is right for your problem. Maybe genSpectra could have an option to use Re, Im or Abs? I think I might have been guilty of pushing the code beyond what you intended as an example of TDHF?

I tried to verify that the area under the peaks is the oscillator strength. By eye the peaks are between 10-22eV and 44-48eV. These are the spec array locations n*20000/2*27.211, so using np.trapz(spec[start:end]) I get a ratio of ~7:1 which isn't anywhere near 0.6685:0.0381. I must have misunderstood how trapz works ! I don't want to bother you with this any longer so if you're happy then we'll leave this @jjgoings. Peter