[Octopus-users] Time propogation problems
Robertson.Burgess at newcastle.edu.au
Tue Feb 2 00:35:17 WET 2010
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
ConvAbsDens = 1e-4
PPScheme = tm
0 | 2.2
1 | 3.0
2 | 1.2
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:
LRConvAbsDens = 1e-3
LinearSolverTol = 1e-5
SmearingFunction = fermi_dirac
Smearing = 0.01
0.01868 | 0.99958 | -0.02206
-0.63010 | -0.33635 | -0.69989
0.91751 | -0.35555 | -0.17820
And my structure file, Au.ANI, is as follows:
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
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.
University of Newcastle, Australia
More information about the Octopus-users