[Octopus-users] Time propogation problems

[Robertson Burgess]

Dear Octopus Users,

Firstly I apologise in advance, I expect this email will be quite long.
I'm trying to calculate the optical response of large gold clusters. As a result of what I learnt at the School in Benasque last January, I've decided to go back to the tetrahedral 20 atom cluster and attempt to reproduce the results of Castro et. al. JCP 144110 (2008) before moving back to larger clusters. I also wanted to take advantage of the symmetry of the cluster and hence only use one calculation. The problem with this for the tetrahedron is that the equivalent axes are not orthogonal and so calculating the symmetry operations was a bit tricky. Unfortunately after I was finished my absorption coefficient is non-physical, fluctuating between negative and positive values. There's a couple of places where I could have gone wrong, but I'm not sure.

I used the APE code to generate pseudopotentials, attempting to follow the scheme set out by Engel et. al. Phys Rev B, 125121 (2001). My inp.ape file (including both variables for the all-electron run and the PP generation) looks like this:

Title = "Au"
CalculationMode = pp
Verbose = 30

WaveEquation = dirac

NuclearCharge=79

ConvAbsDens = 1e-4
Mixing=0.3
MaximumIter=300

PPOutputFileFormat=upf

PPScheme = tm

%PPComponents
 0 | 2.2
 1 | 3.0
 2 | 1.2
%

%Orbitals
 "Xe"
 4 | 3 |  7 | 7
 5 | 2 |  5 | 5
 6 | 0 |  1 | 0
 6 | 1 |  0 | 0
%


Where I tried calculating the PP using both PPScheme = tm and PPScheme = rtm to investigate the effects.
After this I ran a ground state run of octopus, using spin-polarised for the PPScheme = rtm PP and spin unpolarised for the PPScheme = tm PP. In terms of ground state energy levels the two runs barely differ. They are significantly different, however from my earlier ground state calculation using the HGH pseudopotential that came with octopus. This is why I think that maybe my problem is with my pseudopotential.

My octopus input file is as follows:

calculationmode=td
DevelVersion=True

units=2
%species
'Au'|196.9666|spec_ps_upf|79|3|0
%

extrastates=10

XYZCoordinates='Au.ANI'

ConvAbsDens=5e-6
EigensolverFinalTolerance=1.e-7
EigensolverInitTolerance=1.e-5
BoxShape=3

Radius=3.5

SpinComponents=2
Spacing=0.20
TypeOfMixing=2
What2Mix=1

FromScratch=True

periodicdimensions=0

LRConvAbsDens = 1e-3
LinearSolverTol = 1e-5
Mixing=0.1
ParallelizationStrategy=2
EigenSolverMaxIter=200

SmearingFunction = fermi_dirac
Smearing = 0.01

MaximumIter=100

TDDeltaStrength=0.01
%TDPolarization
  0.01868 |  0.99958 | -0.02206
 -0.63010 | -0.33635 | -0.69989
  0.91751 | -0.35555 | -0.17820
%
TDPolarizationDirection=1
TDPolarizationEquivAxis=3

TDExponentialMethod=lanczos
TDExpOrder=100

tmax=30
TDEvolutionMethod= aetrs
TDTimeStep =0.010
TDMaximumIter=tmax/TDTimeStep

And my structure file, Au.ANI, is as follows:

20
 ****************
      Au                   6.701665    7.644312  -12.291109
      Au                   9.262230    7.511184  -11.585020
      Au                  11.744520    7.477735  -10.747932
      Au                   7.661670    9.919073  -11.303238
      Au                  10.323387    9.844710  -10.544820
      Au                   8.701111   12.063421  -10.218144
      Au                  14.212748    7.553470   -9.758498
      Au                  11.457858   12.030503   -9.284862
      Au                  12.899834    9.855944   -9.535751
      Au                   9.847969   14.129839   -9.003509
      Au                   8.269814    7.786824   -4.520125
      Au                  12.315103    7.523030   -7.896841
      Au                  10.911745    9.891688   -7.614757
      Au                   9.278478   12.113773   -7.365947
      Au                  10.350900    7.598816   -6.165846
      Au                   8.759664   10.016064   -5.884352
      Au                   7.596525    7.639833   -7.089721
      Au                   9.723153    7.450242   -8.856109
      Au                   7.081880    7.591330   -9.658726
      Au                   8.082900    9.928859   -8.567959


At the time that I ran the time propogation, I had not yet calculated my TDPolarizationWprime and so hadn't entered it in yet. I calculated it later and entered it manually into the multipoles file before running oct-cross-section.

My three polarisation vectors were calculated as follows:
vector 1 is the vector running from atom 10 through atom 18, vector 2 runs from atom 1 through atom 13, and vector 3 runs from atom 7 to atom 20. Each of these axes runs though one of the corner atoms and out though the atom in the middle of the opposite face.
The symmetry operation A I defined to be a 180 degree rotation about the axis that bisects the angle made by vector 1 and vector 2. This axis therefore bisects the edge running between the corners that vector 1 and 2 come out of, and proceeds to bisect the edge opposite. This rotation axis I calculated to be
[-0.541861290107313542, 0.574110693038497377, -0.614262905137877180]

And so the Rotation Matrix A is calculated as
[-0.413087071435233044, -0.621843623431909354, 0.665334186091203206]
[-0.621843623473128049, -0.341146746808180135, -0.704932198794475884]
[0.665334186012490614, -0.704932198777862728, -0.245766181141780810]

Given that the matrix represents a rotation of 180 degrees one would expect that A^2 would be the identity matrix, and such is the case, within several many significant figures. As such the inverse of Matrix A is itself.

Applying Matrix A to vector 3 I obtained Wprime as

[-0.276478173280484906,-0.323634804281765331,0.904886917420508530]

I then put this into my multipoles input file and ran oct-cross-section. As I said earlier, the resulting (1/3)*Tr(sigma) fluctuated about 0, going to quite large negative values (<-100). Unfortunately I don't understand what I did wrong and why my results differ so dramatically from what I've been able to calculate before. As far as I can see, the problem must come down to one of two things;

My pseudopotentials are bad
My axes and/or Wprime are bad.

These are really the only two things that I have changed from before. I have adopted the Lanczos scheme this time and increased my propagation step, but I checked these my running a time propagation without a perturbation and the energy remained stable over 30 iterations.

I do know that this is a long email and I apologise. Anything at all to help me track down the source of these errors would be greatly appreciated. What can result in highly negative absorption coefficients? Previously for this structure I calculated all three components separately, but especially as I move to larger structures it will be very useful to be able to correctly use the symmetry to reduce my calculation times.

Sincerely,
Robertson Burgess
University of Newcastle, Australia