Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

OptaDOS photoemission module #38

Open
wants to merge 3 commits into
base: develop
Choose a base branch
from

Conversation

victorchanglee
Copy link

Dear all,

This includes a new module for simulating photoemission in OptaDOS.

Details of the implementation are included in documents/photoemission

A .cell .param and .odi example files are also included in examples/Photoemission

Please, let me know about any issue or question about how to use the module.

Kind regards,
Victor

This pull request includes a new module for simulating photoemission in OptaDOS.

Details of the implementation are included in documents/photoemission

A .cell .param and .odi example files are also included in examples/Photoemission
@ajm143
Copy link
Contributor

ajm143 commented May 20, 2021

@victorchanglee I'm confused by a ddome_bin file? What does it contain, and where does it come from? An ome_bin is the optical matrix elements. A dome_bin is the diagonal optical matrix elements, so what does ddome mean? And is that a consistent naming scheme?

Copy link
Contributor

@ajm143 ajm143 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Looks really nice and useful.

These are some initial thoughts on the pull request. I will do more later, but I think it's good for @victorchanglee to have a look. And @Rebecca-Nicholls also needs to see the comments.

integer :: max_layer
real(kind=dp) :: e_fermi
real(kind=dp) :: q_weight
real(kind=dp), parameter :: epsilon_0 = 8.8541878176E-12_dp !A^2 s^4 kg^-1 m^-3
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Constants need to live in constants.f90 and be consistent throughout the code, please.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Thanks, will move that to the appropriate file. It was initially coded so that the original optados code is modified as little as possible

optados/src/photo.f90 Show resolved Hide resolved
@@ -0,0 +1,50 @@
@article{DFT-HK,
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This shouldn't be checked in.

@@ -0,0 +1,46 @@
# CELL file written 10:16:14 (GMT+2.0) 18th July 2020 from run Si2
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This shouldn't be checked in.

@@ -0,0 +1,20 @@
TASK : photoemission
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Might be worth commenting this file more, to explain what the new keywords do here.

@@ -0,0 +1,970 @@
\documentclass[a4paper,11pt,twoside]{book}
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not convinced that we should have a separate manual. I'm quite happy for an extra chapter in the manual.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Our longer term plan is to have an html / online manual. But for the moment, I think this should just be a chapter of the manual.

SPIN_POLARIZED : FALSE
WRITE_FORMATTED_POTENTIAL : TRUE
WRITE_FORMATTED_DENSITY : TRUE
MIXING_SCHEME : Pulay
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Are all of these parameters necessary? It would be helpful not to have red-herrings in here.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

You are right, they are not a requirement. They has to do more with other results in my work.

+============================================================================+
+ Projected Dispersion Curve Calculation +
+============================================================================+

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

By all means have a section here that says, "Photoemission by V.C. N.H. et al, please cite..."

@@ -0,0 +1,17 @@
module od_build
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Please don't check this in.

return

101 call io_error('Error: Problem opening cst_vel file in read_band_curvature')
102 call io_error('Error: Problem opening dome_bin file in read_band_curvature')
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Surely not a dome_bin ?

@Rebecca-Nicholls
Copy link
Contributor

@victorchanglee Hi Victor - thank you for all of this! I'm in the middle of teaching at the moment but should be able to look at it in ~three weeks time. Andrew has sent you some comments in the meantime. Best wishes, Rebecca

@victorchanglee
Copy link
Author

@victorchanglee I'm confused by a ddome_bin file? What does it contain, and where does it come from? An ome_bin is the optical matrix elements. A dome_bin is the diagonal optical matrix elements, so what does ddome mean? And is that a consistent naming scheme?

@ajm143 It contains the second order kp perturbation, which is the band curvature/effective mass. This is currently under "testing"

It is related to how the parallel momentum in ARPES is interpreted. In ARPES, the angle is determined by the parallel momentum of the electrons inside of a crystal. It is usually assign to the crystal momentum or wavevector k. However, the nature of the Bloch theorem makes k only preserved for a defined unit cell or in other words, if we define a supercell, the bands will be folded and thus changing the interpretation of k.

Related to this issue, several band unfolding schemes has being developed (for example, CASTEP has a band unfolding tool), however, in my thesis I took a different look at this momentum interpretation of the band structure.

I used two different approaches, the first one is applying the momentum operator to get the momentum expectation value, which is in fact, the kp first order perturbation already computed in CASTEP in .dome

The second approach I used is a semi-classical momentum, for which I used the band gradient to obtain the electron velocity and the band curvature to get the effective mass. The band curvature is computed as the second order perturbation in the kp method, for which Phil Hasnip and Peter Byrnes provided me with a piece a code to modify Castep.

There is a key word "momentum" that switch between one momentum interpretation or the other. It is under testing since both three approaches gives different answers.

@ajm143
Copy link
Contributor

ajm143 commented May 25, 2021

@pjh503 This is technically a CASTEP issue. But seems to fit better with this PR.

We seem now to have a ddome filename. I think that name could be confusing.

  • ome are the optical matrix elements j,i
  • dome is the diagonal optical matrix elements I=j. Hence d stands for diagonal
  • ddome second order k.p perturbation, I=j. Hence the first d stands for derivative, and the second d for diagonal.

I'm worried that this might get messy if we use d to mean two different things. What are your thoughts?

@pjh503
Copy link
Collaborator

pjh503 commented May 28, 2021

I agree that ddome_bin would be confusing, the question is what it should be called. I prefer "d" to mean "derivative", but unfortunately we've already used it for "diagonal".

I would also prefer the name to reflect the fact that it is a derivative with respect to "dk", as we can imagine cases whether we look at the derivative with respect to other things (e.g. strain, phonon mode).

@ajm143
Copy link
Contributor

ajm143 commented May 28, 2021

.dome_dk ? I agree we're using d to mean two different things here... maybe we get away from omes and call it de_dk or similar?
.bands_dk ?

Copy link
Contributor

@ajm143 ajm143 left a comment

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

The looks good, and is almost totally in the OptaDOS coding style (despite us never publishing such).

There seems to me to be a lot of redundant code, in photo.f90 that could just call other functions instead. You shouldn't feel bad about generalising functions in other modules to make them do what we both want. It will be much easier for us all to maintain.

How do we test your code? What version of CASTEP generates the necessary OptaDOS input files?

I'm very keen to get this code into the stable branch.


epsilon_zero = 55.26349406_dp
ge = 10.24633_dp
kB=8.617333262E-5_dp
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Put constants in the constant module please.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I moved all the constants to the constant module.

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

How do we test your code? What version of CASTEP generates the necessary OptaDOS input files?

I will prepare some input/output examples for the code. There are two input files required to run some of the options in the code that are generated by a modified version of CASTEP 18.1

One is the free electron matrix, which I generated by modifying the spectral.f90 in CASTEP, in which the final state is substituted by a free electron planewave.

the second file is the band curvature which I also modified the spectral.f90 in CASTEP for this purpose.

I think ideally both files should be written as a tool file in CASTEP but at the time, it was faster and easier for me to just modify the spectral.f90

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think the CASTEP code freeze is coming up. When will your modifications to CASTEP be in a release ?

width = (1.0_dp/11604.45_dp)*temp
qe_factor=1.0_dp/(4*pi**2*s_area)
norm_vac = gaussian(0.0_dp,width,0.0_dp)
kB=8.617333262E-5_dp
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

put contestants in constants module please.

optados/src/photo.f90 Outdated Show resolved Hide resolved

real(kind=dp) :: z,epsilon_zero

epsilon_zero = 55.26349406_dp
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

please put the constants in constants.f90

@@ -0,0 +1,970 @@
\documentclass[a4paper,11pt,twoside]{book}
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Our longer term plan is to have an html / online manual. But for the moment, I think this should just be a chapter of the manual.


!***************************************************************
subroutine calc_absorp_photo
!***************************************************************
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Copy ?

!***************************************************************
subroutine calc_reflect_photo
!***************************************************************
! This subroutine calculates the reflection coefficient
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Copy?


!===============================================================================
subroutine setup_energy_scale
!===============================================================================
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This is a copy? Why not call the original ?

!===============================================================================
subroutine jdos_utils_calculate_delta(delta_temp)
!===============================================================================
! Copy from the jdos_utils.f90
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Why not just call jdos_utlis ?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This function is slightly trickier to call.

In the original jdos_utils, the delta function is computed on the fly (it is not allocated) in a variable called dos_temp to generate the weighted_jdos which is integrated over the BZ.

In photoemission, I want to generate the photoemission intensities as function of k, so I need the Dirac delta function to be allocated.

I think the solution is to allocate the dos_temp in calculate_jdos.

end subroutine transverse_energy_spread

subroutine binding_energy_spread
! This subroutine applies a Gaussian broadenning to the binding energy
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This looks like a copy of the function above? Why not make a function that does both ?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I made changes so now it only generates a plot of the photoemission intensities (QE) as a function of the binding energy, which is what experimentalist usually obtain from UPS/XPS/ARPES.

Several changes:

- Constants moved to constants.f90

-Merged the optical properties in the photoemission properties with the now public functions in optics.f90.

There are still a few functions that are copies with small changes.

I included all the input requires to run the photoemission module for a graphene sheet, including the free electron matrix.

I will try to update the code with further small changes.
@@ -405,7 +405,7 @@ subroutine elec_read_band_curvature

if(on_root) close (unit=curvature_unit)

! Convert all band curvatures to eV Ang
! Convert all band curvatures to eV Ang^2
band_curvature=band_curvature*bohr2ang*bohr2ang*H2eV
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Be very careful here. OptaDOS does a double conversion here bohr2ang twice because of an error in the units of the optical matrix elements coming out of CASTEP. Are you sure that the band curvature has the same error?

Copy link
Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it has same error. I compared the calculated band curvature from a near-free electron band calculated as a chain of Li atoms in CASTEP with a polynomial fitting of the same band.

After unit conversion, the difference in curvature with both methods is around 0.02, which I think it is more likely to be an issue of the polynomial fitting.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

Successfully merging this pull request may close these issues.

None yet

5 participants