[Octopus-users] ionization probabilities by TDKLI

[Nam A. Nguyen]

Well, after many sleepless nights, I finally found an aswer to my own 
question. I think that the octopus developers also detected the error a 
couple of days ago and probably they have corrected it. I'm posting here 
nevertherless in case some of the users want to do OEP calculation, then 
they should be aware of a possible numerical error.

Problem:
Propagate a system with ZERO external field with TD-KLI and the density 
and energy will change. They should not do so.

Answer:
the local xc potential, oep%lxc should be complex in the time dependet case

Solution
1- modify the real oep%lxc to  complex. This would involve declaring 
oep%zlxc and recompiling octopus
2- download the latest version

Happy TDDFT number crunching everyone !
Nam

>Dear octopus users,
>
>I've been experimenting with TD-OEP/KLI. (see the attached input file). 
>Basically I propagate  the H2 system in ZERO fields (i.e. intensity=0). 
>Then I calculate the ionization probabilities by integrating the 
>electron density in a certain inner region of the space.
>
>Here's the problem:
>
>at time=0, the single-electron ionization probability is about 10^{-12} 
>which makes perfect sense.
>However, at time=T (about some femto seconds), the probability becomes: 
>  0.1, which does NOT make sense because in a ZERO field there should 
>not be any ionization at all.
>
>I also did the same calculation with LDA, not using KLI scheme, and got 
>a correct result.
>
>So my conclusion is that the finite ionization probability in TD-KLI is 
>a numerical error originated from the propagation. Any idea how I can 
>reduce this error ? smaller timestep ?.
>
>Comments, thoughts ?
>
>Best regards
>Nam
>
>-----------------------------
>verbose = 1130
>DebugLevel = 3
>SystemName = "H2 OEP calculation"
>CalculationMode = 5
>Dimensions=1
>Units = "a.u"
>UnitsInput  = a.u
>UnitsOutput = a.u
>
>%Species
>  "H"  | 1.0 | 1 | "-1/sqrt(x^2+0.7)"
>%
>
>%Coordinates
>  "H" | 0.85 | 0.0 | 0.0 | no
>  "H" | -0.85 | 0.0 | 0.0 | no
>%
>
>#BoxShape=4
>#%Lsize
>#100|10|10
>#%
>#%Spacing
>#0.4|0.4|0.4
>#%
>
>#PoissonSolver=3
>
>Radius=1000
>Spacing=0.4
>
>XFunctional='SIC'
>CFunctional='SIC'
>
>#SpinComponents=2
>#ExtraStates=10
>LCAOStart=no
>#AbsorbingBoundaries=2
>
>%MakeStates
>1|1|1|1|0.9|0.7|0.0
>1|2|1|1|0.9|-0.7|0.0
>%
>
>
>%TDLasers
>1|0|0|0.0000|0.05841455|3|214|321|214
>%
>
>#FFTOptimize=no
>
>TDEvolutionMethod=magnus
>TDExponentialMethod=split
>
>TDTimeStep=0.03
>TDMaximumIter=25000
>
>#TDOutputOccAnalysis=1
>#UnOccNumberStates=30
>UnoccMaximumIter = 25000
>UnoccConv = 0.000001
>GuessDensitySpin=2
>
>OutputDensity = yes
>OutputWfs=yes
>OutputKSPotential=yes
>OutputWfsNumber="1-5"
>#OutputDX=yes
>#OutputPlaneZ=yes
>OutputAxisX=yes
>OutputEvery=5000
>
>SpecMaxEnergy=300
>HSPolarization='x'
>