[REVIEW]: SHGYield

[whedon]

Submitting author: roguephysicist (Sean M. Anderson)
Repository: https://github.com/roguephysicist/SHGYield
Version: v1.0.1
Editor: @kyleniemeyer
Reviewer: @nicoguaro
Archive: 10.5281/zenodo.803483

Status

Status badge code:

HTML: <a href="http://joss.theoj.org/papers/2988b11612a03c1d15eab588f817106d"><img src="http://joss.theoj.org/papers/2988b11612a03c1d15eab588f817106d/status.svg"></a>
Markdown: [![status](http://joss.theoj.org/papers/2988b11612a03c1d15eab588f817106d/status.svg)](http://joss.theoj.org/papers/2988b11612a03c1d15eab588f817106d)

Reviewers and authors:

Please avoid lengthy details of difficulties in the review thread. Instead, please create a new issue in the target repository and link to those issues (especially acceptance-blockers) in the review thread below. (For completists: if the target issue tracker is also on GitHub, linking the review thread in the issue or vice versa will create corresponding breadcrumb trails in the link target.)

Reviewer questions

Conflict of interest

• As the reviewer I confirm that there are no conflicts of interest for me to review this work (such as being a major contributor to the software).

General checks

• Repository: Is the source code for this software available at the repository url?
• License: Does the repository contain a plain-text LICENSE file with the contents of an OSI approved software license?
• Version: Does the release version given match the GitHub release (v1.0.1)?
• Authorship: Has the submitting author (roguephysicist) made major contributions to the software?

Functionality

• Installation: Does installation proceed as outlined in the documentation?
• Functionality: Have the functional claims of the software been confirmed?
• Performance: Have any performance claims of the software been confirmed?

Documentation

• A statement of need: Do the authors clearly state what problems the software is designed to solve and who the target audience is?
• Installation instructions: Is there a clearly-stated list of dependencies? Ideally these should be handled with an automated package management solution.
• Example usage: Do the authors include examples of how to use the software (ideally to solve real-world analysis problems).
• Functionality documentation: Is the core functionality of the software documented to a satisfactory level (e.g. API method documentation)?
• Automated tests: Are there automated tests or manual steps described so that the function of the software can be verified?
• Community guidelines: Are there clear guidelines for third parties wishing to 1) Contribute to the software 2) Report issues or problems with the software 3) Seek support

Software paper

• Authors: Does the paper.md file include a list of authors with their affiliations?
• A statement of need: Do the authors clearly state what problems the software is designed to solve and who the target audience is?
• References: Do all archival references that should have a DOI list one (e.g. papers, datasets, software)?

[whedon]

Hello human, I'm @whedon. I'm here to help you with some common editorial tasks for JOSS. @nicoguaro it looks like you're currently assigned as the reviewer for this paper 🎉.

⭐ Important ⭐

If you haven't already, you should seriously consider unsubscribing from GitHub notifications for this (https://github.com/openjournals/joss-reviews) repository. As as reviewer, you're probably currently watching this repository which means for GitHub's default behaviour you will receive notifications (emails) for all JOSS reviews 😿

To fix this do the following two things:

1. Set yourself as 'Not watching' https://github.com/openjournals/joss-reviews:

2. You may also like to change your default settings for this watching repositories in your GitHub profile here: https://github.com/settings/notifications

For a list of things I can do to help you, just type:

@whedon commands

[roguephysicist]

Thanks @nicoguaro and @kyleniemeyer. I am available for any questions you might have, though the theory behind the program is completely described in the references listed in the README (particularly Phys. Rev. B 94, 115314 (2016)). The most recent publication using this software has been accepted but is still undergoing final processing. You can find a provisional version here.

[nicoguaro]

As a first comment, I would say that the main file should have docstrings to understand better the functions that appear there. As a second one, a description of the input and output files is missing, what makes really difficult to use the software for other cases different to the provided example.

[roguephysicist]

I can definitely see that. I will improve the documentation in the main file and add documentation to the input file.

[roguephysicist]

I have added many comments to both the main and input files. Hopefully that will make things a lot easier to understand.

[nicoguaro]

I still don't get the sample-data input files thing. For example

## chi1 & chi2 components. Can be path to files, or number (like 0).

How come this can be a single number like 0, or a whole table like in the file sample-data/SiBulk-chi1-xx_yy_zz.

I understand that you mention this

You must first calculate the different components of the nonlinear susceptibility tensor. The theory surrounding this problem is still being developed, and there are many ways to go about it. We leave it to the reader to find the best method for their particular problem.

But let's suppose that I can "magically" know this tensor (that should have 27(?) components for every energy value). I still don't know how to format that information into inputs for the whole process.

[roguephysicist]

Some basic theory on the susceptibility tensor

chi2 (chi1) is the nonlinear (linear) susceptibility tensor.

In particular, chi2 determines the nonlinear polarizability of a material and is responsible for second-harmonic generation. This relationship is expressed as

where a, b, and c are crystallographic directions that depend on how you orient your crystalline structure. We can see that a material can produce a polarization response in direction a from two incident fields (E) in directions b and c, by means of chi2^{abc}.

chi2 is a second-order tensor, and thus has 27 possible components (unique combinations of a, b, and c; for instance, aaa, aab, and so on.). Second-harmonic generation implies that the incoming fields are identical (two photons of equal energy in, one photon of double-energy out) so it is also implied that

for this particular phenomenon. This reduces 9 of the possible combinations, reducing to 18 unique components. It is very convenient to express the crystallographic directions in terms of x, y, and z; therefore, we can express chi2 with all 18 components as

Symmetry relations are very important for determining chi2. A given crystal symmetry can greatly reduce the complexity of the problem by eliminating many of the components. For instance, for the (001) face of cubic crystals, we have that

which has only 3 independent components. There are many articles and books with tables and extensive discussion that describe the form that this tensor should have for a given symmetry.

As far as the program is concerned, you can indicate in the input file that you have all 18 components for your structure (maybe the symmetry of your material is unknown, so you calculate all of them), or that some of them are 0 (the program assumes that any of them NOT listed are 0).

The case for chi1 is considerable simpler. chi1 is directly related to the dielectric function of the material

which is directly related to the index of refraction as

The chi1 spectra should obviously have non-zero regions; otherwise, the problem is not very interesting.

[roguephysicist]

Concerning the calculation of the nonlinear susceptibility tensor

I think that the other part of your question concerns how we actually calculate chi2. That is definitely a non-trivial problem, and as I mentioned in the README, there are different ways to approach it. The method that we use is derived in full detail here, although that reference does not really touch upon the relevant computational aspects.

This response ended up a little long-winded, so the last section essentially summarizes what I think you are asking.

On calculating the electronic properties of semiconductors

In general, there are several ab initio formalisms that can be used to obtain the wavefunctions/electron densities/energies of crystalline semiconductor materials. The most common are:

• all-electron methods (older, declining in popularity)
• density functional theory with the local-density approximation (DFT-LDA) (solid workhorse with known problems, used in this work)
• time-dependent density functional theory (TDDFT) (good for optical responses, but hard to implement)
• many-body perturbation theory (MBPT) (state-of-the-art, but extreme computational expense)

Each of these has their pros and cons that mainly relate to accuracy vs. difficulty of the theory vs. computer resources. There are many free and open-source codes available for download, that are under active development by thousands of researchers and groups.

On calculating the susceptibility tensors and optical properties

Once the initial electronic properties of the material are determined, we can then proceed to calculate the optical properties that include the linear and nonlinear responses. In general, nonlinear optics is now pretty well understood within the DFT-LDA (see ref. above and refs. within) and the TDDFT frameworks. MBPT is currently (AFAIK) still at the linear optics level, but has the highest level of accuracy available. Of course, you can combine different methods to exploit the strengths that each has to offer.

Unfortunately, most groups do not publish the software codes used to calculate the susceptibility tensor. The code that we use is called TINIBA; however, it is VERY unpolished and requires a major overhaul before it could ever be "officially" released/published. That said, we are an active research group that uses this software every day to produce high-quality scientific work, and are constantly improving and adding features to it. Publishing TINIBA is definitely a long-term goal that we have.

There are, however, some codes available, such as RT-SIESTA which works with TDDFT, and even ABINIT has built-in utilities for calculating the susceptibility tensors.

On calculating the SHG yield

In summary, SHGYield.py produces the SHG radiation (in reflectance) produced from a semiconductor surface. As you mentioned, it requires the susceptibility tensors that are calculated with the above considerations. The theory is developed considering a reflectance model with three distinct regions that allows the user to readily simulate the SHG response of thin-films over bulk substrates, or of any crystalline surface with any symmetry considerations (see my previous reply to this thread). The program does NOT care how you produced the susceptibility tensors; you can use MBPT for producing the linear susceptibility tensor and then use TDDFT to produce the chi2 components, etc.

This is why I think that this program has value and is worth being published - it tackles a research problem (the SHG radiation from a surface), but allows the researcher to use their data that can be produced using a variety of frameworks and theoretical developments.

[nicoguaro]

@roguephysicist, these are good explanations, and I think that you should include them in the documentation of your software. I little detail, though, I suppose that you meant third-order tensor, since chi2^{abc} has three free indices. On the other hand, chi1^{ab} is a second-order tensor, since it has two free indices.

So, you answered part of my question. But some part is still missing, and is the formatting of these files. And that was my main concern, because right now is not clear on how to format the information of the tensors so that the program read the properly.

This is why I think that this program has value and is worth being published - it tackles a research problem (the SHG radiation from a surface), but allows the researcher to use their data that can be produced using a variety of frameworks and theoretical developments.

I also think that it is worth being published. But it is important that the documentation is good enough for people to use it.

[roguephysicist]

I suppose that you meant third-order tensor, since chi2^{abc} has three free indices. On the other hand, chi1^{ab} is a second-order tensor, since it has two free indices.

I thought one thing and wrote another! chi2 is a second-order nonlinear susceptibility, but is a third-rank tensor. Thanks for spotting that.

So, you answered part of my question. But some part is still missing, and is the formatting of these files. And that was my main concern, because right now is not clear on how to format the information of the tensors so that the program read the properly.

Sorry for missing that. The chi2 components are dealt with in the shgcomp() function on line 76. The file is organized in 5 columns:

1. 1w energy from 0-20 eV in 0.01 eV steps,
2. real part of the 1w resonances
3. imaginary part of the 1w resonances
4. real part of the 2w resonances
5. imaginary part of the 2w resonances

We separate the 1w and 2w spectra because they are computed from separate equations, but in the end chi2^{abc} = Re[chi2(w) + chi2(2w] + i*Im[chi2(w) + chi2(2w] (line 96). The function also takes care of the scale and units that are produced in TINIBA (line 97), and can easily accommodate any set of units the user might want to use.

chi1 is handled in the epsilon() function. We calculate all three directions of chi1 (xx, yy, and zz) in the same file. Those are organized as:

1. 1w energy from 0-20 eV in 0.01 eV steps,
2. real part xx
3. imaginary part xx
4. real part yy
5. imaginary part yy
6. real part of zz
7. imaginary part of zz

chi1 gets averaged and converted to epsilon per the equations mentioned in my previous reply. This all happens on lines 70 -73.

I also think that it is worth being published. But it is important that the documentation is good enough for people to use it.

You are exactly right, thanks. I will add some of this discussion to the README to really beef up the documentation.

[roguephysicist]

Hey @nicoguaro, I just updated the repo with the following changes.

I expanded the README with a lot of the information I posted above. I changed the input file to YAML which greatly improves the parsing abilities, and added a lot of comments to both the main file and the input file. I improved the way that the main file interprets the chi2 components, so it's smarter and more flexible now.

Let me know what you think.

[nicoguaro]

I think that it looks cleaner now, and the choice of YAML seems to be aligned in the same direction, as well. It could be useful to do the same with the output, a Jupyter Notebook reproducing the plots from above might be good for that purpose.

Regarding automatic testing, it seems that it still does not cover this. The same goes for the community guidelines for contribution.

[roguephysicist]

Thanks for these good suggestions.

I added a minimal test script in the test directory, and am working on a full fledged example with a Jupyter tutorial that includes a couple of different calculations. I'm still working on the Jupyter notebook (my first time), but most of the other stuff is already in place. I have no experience writing test programs, so I appreciate any feedback you might have. I also added a CONTRIBUTING.md file.

I have a few questions about the "JOSS paper" format. Should I enhance it with some of the stuff in the main README, or is it best to keep it simple? Also, does the LaTeX conversion happen automatically, or is that something I need to add in manually?

[kyleniemeyer]

@roguephysicist you are welcome to edit your paper, especially if @nicoguaro thinks it could use more content. I took a quick look, and agree it could use a bit more. Please look at some examples of published JOSS papers.

The LaTeX conversion does happen automatically, so you don't need to worry about that.

[nicoguaro]

@roguephysicist, you could check steps 5 - 7 from this link. Basically, it is good to automate the testing and use something like nosetest.

Regarding your tests this is the output I get

Loading reference data.
Calculating SHG yield for comparison.
Traceback (most recent call last):

  File "<ipython-input-1-bc5c09d2d0d0>", line 1, in <module>
    runfile('D:/workspace/SHGYield/test/test.py', wdir='D:/workspace/SHGYield/test')

  File "C:\Anaconda2\lib\site-packages\spyder\utils\site\sitecustomize.py", line 866, in runfile
    execfile(filename, namespace)

  File "C:\Anaconda2\lib\site-packages\spyder\utils\site\sitecustomize.py", line 87, in execfile
    exec(compile(scripttext, filename, 'exec'), glob, loc)

  File "D:/workspace/SHGYield/test/test.py", line 13, in <module>
    test = np.loadtxt('test-output.dat', unpack=True)

  File "C:\Anaconda2\lib\site-packages\numpy\lib\npyio.py", line 858, in loadtxt
    fh = iter(open(fname, 'U'))

IOError: [Errno 2] No such file or directory: 'test-output.dat'

I also think that you could add some details to the paper. Particularly, I would add the images.

[roguephysicist]

@nicoguaro, I will take a look and improve the tests. I will need to try the program on a windows machine to figure out the problem. So far, it has worked fine for me under Linux and macOS. I will hopefully be able to do this next week.

I will enhance the JOSS paper as well according to both your suggestions. Thanks!

[tacaswell]

@tacaswell pinging my self so I get updates. I will be a second reviewer on this.

I agree on automated tests. It looks like you already have a good reference result (the hard part!), using np.assert_almost_equal instead of plotting should be easy. I advise against nose, the authors have declared it effectively un-maintained, I suggest pytest instead.

[roguephysicist]

@nicoguaro @tacaswell, sorry for the absence -- work and long weekends have gotten in the way.

@tacaswell, Thanks for the suggestions -- I definitely agree that a numeric comparison is much simpler than plotting. I already implemented a small test using pytest and it was indeed very easy. I will produce a simple test suite encompassing the most important components of the program.

@nicoguaro, I will start working on a tutorial using a Jupyter notebook. It will walk the user through the different parts of the calculation and compare against the math featured in the references. It should be pretty informative for people looking to learn what this is all about.

I will hopefully be able to complete these tasks next week. Thanks again.

[kyleniemeyer]

@roguephysicist those plans sound great!

[roguephysicist]

Hi all,

I have finally finished the aforementioned additions to the program.

@nicoguaro, I added a Jupyter notebook tutorial in the examples/ directory. I think it contains most of the relevant theory and code, but it was hard to keep it from getting too long with all the math.

@tacaswell, I added a test script in tests/ that works with pytest under both MacOS and Linux (Windows pending). It runs 6 tests that try out different aspects of the program, ending with a numerical comparison against reference data. The run-time on my machine is around 0.50 seconds.

@kyleniemeyer, I added some content from the main README file to the JOSS paper. I looked at other examples, and kept it brief.

Anyway, looking forward to your comments and suggestions. Thanks again.

[nicoguaro]

@roguephysicist, I just ran your test under Windows, this is the output:

============================= test session starts =============================
platform win32 -- Python 3.6.0, pytest-3.0.5, py-1.4.32, pluggy-0.4.0
rootdir: C:\Users\nguarinz\workspace\SHGYield\tests, inifile:
collected 0 items / 1 errors

=================================== ERRORS ====================================
________________________ ERROR collecting self_test.py ________________________
self_test.py:112: in <module>
    EPS = test_epsilon() # Creates the dictionary with the epsilons
self_test.py:24: in test_epsilon
    epsl = read_eps('tests/data/SiH1x1-epsilon') # Epsilon from chi1, layered, normalized
self_test.py:11: in read_eps
    data = np.loadtxt(in_file, unpack=True, skiprows=1)
..\..\..\Anaconda3\lib\site-packages\numpy\lib\npyio.py:805: in loadtxt
    fh = iter(open(fname))
E   FileNotFoundError: [Errno 2] No such file or directory: 'tests/data/SiH1x1-epsilon'
!!!!!!!!!!!!!!!!!!! Interrupted: 1 errors during collection !!!!!!!!!!!!!!!!!!!
=========================== 1 error in 0.47 seconds ===========================

The rest seems to be fine!

[roguephysicist]

@nicoguaro, you have to run "pytest" in the root/ directory, not in the tests/ directory.

[nicoguaro]

@roguephysicist, you are right. It works nicely now.

You have my blessing.

[roguephysicist]

Thanks @nicoguaro. I don't know why I said `root/` but I'm glad it worked out anyway 😅.

[kyleniemeyer]

@roguephysicist Your changes look great, between the paper, example Jupyter notebook, and addition of tests. Thanks! Could you now give the DOI for an archived version of the (accepted) software, via (e.g.) Zenodo?

@nicoguaro for completeness, could you check off the review items at the top of the issue?

@arfon Do we need to update the version of the software in this issue? The software is now at v1.1, vs. v1.0.1 when submitted. (I can't remember if this is important or not...)

[roguephysicist]

Thanks @kyleniemeyer. I went ahead and added one more release with all these changes (v1.1.1). The Zenodo DOI for this release is 10.5281/zenodo.803483. It does not resolve yet, but should be working in a few hours.

Let me know if anything else is needed on my part. Thanks again!

[kyleniemeyer]

@whedon commands

[whedon]

Here are some things you can ask me to do:

# List all of Whedon's capabilities
@whedon commands

# Assign a GitHub user as the reviewer of this submission
@whedon assign @username as reviewer

# List the GitHub usernames of the JOSS editors
@whedon list editors

# List of JOSS reviewers together with programming language preferences and domain expertise
@whedon list reviewers

# Change editorial assignment
@whedon assign @username as editor

# Set the software archive DOI at the top of the issue e.g.
@whedon set 10.0000/zenodo.00000 as archive

# Open the review issue
@whedon start review

🚧 Important 🚧

This is all quite new. Please make sure you check the top of the issue after running a @whedon command (you might also need to refresh the page to see the issue update).

[kyleniemeyer]

@whedon set 10.5281/zenodo.803483 as archive

[whedon]

OK. 10.5281/zenodo.803483 is the archive.

[kyleniemeyer]

@arfon ok, this one is good to go, other than the potential version change (which may not be necessary)!

[arfon]

@roguephysicist - Could you move the references you currently have in the paper.md file into a paper.bib file and cite them directly please? (You can read how to do that here)

[roguephysicist]

@arfon, I made an attempt at this, but I don't have any way to test it. Let me know if I need to do something different.

[arfon]

@nicoguaro many thanks for your review here and to @kyleniemeyer for editing this one ✨

@roguephysicist - your paper is now accepted into JOSS and your DOI is http://dx.doi.org/10.21105/joss.00242 ⚡ 🚀 💥

[kyleniemeyer]

🎉

[roguephysicist]

Thanks guys! 🎉