Enhanced readme, changed input file to YAML, and improved functionali…

…ty when interpreting the chi2 components

README.md

@@ -3,21 +3,31 @@ SHG Yield for Semiconductor Surfaces
 
 [![DOI](https://zenodo.org/badge/11697217.svg)](https://zenodo.org/badge/latestdoi/11697217)
 
-SHGYield is a python program for calculating the second-harmonic generation (SHG) yield for semiconductor surfaces. This is useful for conducting theoretical second-harmonic spectroscopy of any crystalline semiconductor or nanostructure. The physical theory behind this phenomenon is presented in the references below. We have applied this software to several silicon test cases and have achieved very accurate results that compare well with experiment.
+SHGYield is a python program for calculating the surface second-harmonic generation (SSHG) yield (in reflectance) for semiconductor surfaces. 
 
-These calculations can predict effects in novel materials that are still under study, or save valuable time and resources for the experimentalist. For example, the figure below is an overview of the angular dependence of the reflected SHG Yield from the Si(111)(1x1)H surface. Experimentalists will find this very useful, as they can plan the experiment accordingly in order to optimize the output signal strength and polarization.
 
-![An overview of the angular dependence of the SHG Yield for the Si(111)(1x1)H surface](paper/3D-Si1x1.png)
+Introduction
+------------------------------------
+
+SSHG is an effective, nondestructive, and noninvasive probe for studying surface and interface properties, and even for characterizing buried interfaces and nanostructures. The high surface sensitivity of SSHG spectroscopy is due to the fact that within the dipole approximation, the bulk SHG response in centrosymmetric materials is identically zero. The SHG process can occur only at the surface where the inversion symmetry is broken. 
+
+This program has several potential applications and uses:
+* determining and analyzing the physical origin of SSHG spectra
+* predicting and characterizing the radiated SSHG for interesting new materials
+* characterizing thin films based on measured SH spectra
+* allowing the experimenter to calculate and analyze the SSHG yield to optimize experiments
 
-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. As an example, we use [ABINIT](http://www.abinit.org) to calculate the electron density and matrix elements and then [TINIBA](https://github.com/bemese/tiniba) to calculate the tensor components. 
+For example, the figure below is an overview of the angular dependence of the reflected SHG Yield from the Si(111)(1x1)H surface. Experimentalists will find this very useful, as they can plan the experiment accordingly in order to optimize the output signal strength and polarization.
+
+![An overview of the angular dependence of the SHG Yield for the Si(111)(1x1)H surface](paper/figures/3D-Si1x1.png)
 
 
 References
 ------------------------------------
 
-The complete theory is derived step-by-step in my [Ph.D. thesis](https://github.com/roguephysicist/thesis-phd). This software has been used in the following publications:
+The complete theory is derived step-by-step in [Phys. Rev. B 94, 115314 (2016)](https://doi.org/10.1103/PhysRevB.94.115314). This software has been developed and used in the following publications:
 
-* [Front. Mater. 4:12 (2017)](http://journal.frontiersin.org/article/10.3389/fmats.2017.00012/abstract)
+* [Front. Mater. 4:12 (2017)](https://doi.org/10.3389/fmats.2017.00012)
 * [Phys. Rev. B 94, 115314 (2016)](https://doi.org/10.1103/PhysRevB.94.115314)
 * [Phys. Rev. B 93, 235304 (2016)](https://doi.org/10.1103/PhysRevB.93.235304)
 * [Phys. Rev. B 91, 075302 (2015)](https://doi.org/10.1103/PhysRevB.91.075302)
@@ -28,15 +38,81 @@ The complete theory is derived step-by-step in my [Ph.D. thesis](https://github.
 Installation
 ------------------------------------
 
-SHGYield has been tested with Python 2 and 3, and Anaconda Python 4+ on Mac OS X and Linux. It should work on any system with the required Python packages installed. 
+SHGYield has been tested with Python 2 and 3, and Anaconda Python 4+ on both macOS and Linux. It should work on any system (including Windows) with the required Python packages installed.
 
 Python requirements:
-`sys`, `math`, `numpy`, `scipy`
+`sys`, `yaml`, `numpy`, `scipy`
 
 Usage:
-`python shgyield.py <sample.in>`
+`python shgyield.py <input.yml>`
+
+A sample input file `input.yml` and data set `sample-data/` are included for your convenience. Both the input file and the program itself are extensively documented.
+
+
+Theory
+------------------------------------
+
+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). 
+
+In order to calculate the SHG yield, you must first calculate the linear and nonlinear susceptibility tensor (<b><i>χ</i></b>(−2ω; ω, ω), <b><i>χ</i></b> for short) for the material of interest. 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. As an example, we use [ABINIT](http://www.abinit.org) to calculate the electron density/wavefunction/energies and then [TINIBA](https://github.com/bemese/tiniba) to calculate the relevant matrix elements and the <b><i>χ</i></b> tensor components. The exact method that we use is derived in full detail in [Phys. Rev. B 91, 075302 (2015)](https://doi.org/10.1103/PhysRevB.91.075302). The program does NOT care how you have produced the susceptibility tensors; you can use different frameworks (such as MBPT, DFT-LDA, TDDFT, etc.) for producing the linear and nonlinear susceptibility tensors. 
+
+
+### The nonlinear susceptibility tensor, <i>χ</i>(−2ω; ω, ω)
+
+
+<b><i>χ</i></b> determines the nonlinear polarizability of a material and is responsible for second-harmonic generation. This relationship is expressed as
+
+<img src="paper/figures/polar.png" height="25">
+
+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 (<i>E</i>) in directions b and c, by means of <i>χ</i><sup>abc</sup>.
+
+<b><i>χ</i></b> 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 
+
+<img src="paper/figures/chi2-shg.png" height="25">
+
+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 <b><i>χ</i></b> with all 18 components as
+
+<img src="paper/figures/chi2-full.png" height="90">
+
+Symmetry relations are very important for determining <b><i>χ</i></b>. 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
+
+<img src="paper/figures/chi2-001.png" height="90">
+
+which has only 3 independent components. There are many [articles](https://journals.aps.org/prb/abstract/10.1103/PhysRevB.35.1129) and [books](https://books.google.com/books?id=bdFju3af2FsC) with tables and extensive discussion that describe the form that this tensor should have for a given symmetry.
+
+
+### The linear susceptibility tensor
+
+The case for <i>χ</i><sup>ab</sup>(ω) is considerable simpler. <i>χ</i><sup>ab</sup>(ω) is directly related to the dielectric function of the material
+
+<img src="paper/figures/epsilon.png" height="25">
+
+which is directly related to the index of refraction as
+
+<img src="paper/figures/refrac.png" height="45">
+
+The <i>χ</i><sup>ab</sup>(ω) spectra should obviously have non-zero regions; otherwise, the problem is not very interesting.
+
+
+### Calculating the electronic properties of semiconductors
+
+In general, there are several *ab initio* formalisms that can be used to obtain the electron density/wavefunction/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)
+* 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](http://psi-k.net/software/) available for download that are under active development by thousands of researchers and groups.
+
+
+### 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 references within [Phys. Rev. B 91, 075302 (2015)](https://doi.org/10.1103/PhysRevB.91.075302)) and the [TDDFT](http://dx.doi.org/10.1063/1.2790014) 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. 
+
+The code that we use is called [TINIBA](https://github.com/roguephysicist/tiniba); however, we have not yet created an official release for production use. 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.
 
-A sample input file `sample.in` and data set `sample-data/` are included for your convenience.
+There are, however, some codes available, such as [RT-SIESTA](http://monalisa.phys.washington.edu/feffproject-rtsiesta.html) which works with TDDFT, and even [ABINIT](http://www.abinit.org/doc/helpfiles/for-v8.2/users/optic_help.html) has built-in utilities for calculating the susceptibility tensors.
 
 
 License

input.yml

@@ -0,0 +1,65 @@
+# Sample input file for calculating SHG yield
+# This is in YAML, so please preserve indentation below!
+
+# Establish reflectance model, see PRB 93, 235304 (2016). Available choices:
+# 3-layer, 2-layer-fresnel, 2-layer-bulk, 2-layer-vacuum, and 3-layer-hybrid.
+# Anything other than '3-layer' will exclude multiple reflections. When in 
+# doubt, use '3-layer'.
+mode:       3-layer
+
+# Required calculation parameters
+parameters:
+    theta: 65           # angle of incidence in degrees
+    phi:   30           # azimuthal angle in degrees
+    sigma: 0.10         # gaussian output broadening in eV
+
+# Multiple reflection framework
+multiref:
+    enable:    yes      # if multiple reflections are taken into account. Can
+                        # be 'yes' or 'no'; if 'no', next values are ignored.
+    thickness: 10       # thickness "d" of the thin layer, in nm
+    depth:     average  # depth "d2" of the nonlinear polarization, in nm. can
+                        # be a float or 'average'. See Eqs (21) and 22 of 
+                        # PRB 93, 235304 (2016).
+
+# chi1: linear susceptibility tensor.
+# Paths to files for chib (bulk) and chil (layered), and normalization of chil.
+# For information about normalizing, see PRB 92, 245308 (2015). File columns are
+# Energy(1w) Re[chi_xx] Im[chi_xx] Re[chi_yy] Im[chi_yy] Re[chi_zz] Im[chi_zz].
+chi1:
+    chib: sample-data/SiBulk-chi1-xx_yy_zz
+    chil: sample-data/SiH1x1-chi1-xx_yy_zz
+    norm: 1.2659296143
+
+# chi2: nonlinear susceptibility tensor 
+# Components should be in separate files, and listed separately. Symmetry 
+# relations can be represented here. There are 4 options for each component:
+#
+# 1. ZERO if it is not listed or commented out (i.e. '# xxx: ...')
+# 2. ZERO if you give it a value of 0, i.e     (i.e. 'xxx: 0')         
+# 3. Loaded from file if path is given         (i.e. 'xxx: ../SiH1x1-chi2-xxx')
+# 4. Has some relation to another component    (i.e. 'xyy: -xxx')
+#
+# File columns should be organized as Energy(1w) Re[1w] Im[1w] Re[2w] Im[2w].
+chi2:
+    xxx: sample-data/SiH1x1-chi2-xxx
+    xyy: sample-data/SiH1x1-chi2-xyy
+    xzz: sample-data/SiH1x1-chi2-xzz
+    xyz: sample-data/SiH1x1-chi2-xyz
+    xxz: sample-data/SiH1x1-chi2-xxz
+    xxy: sample-data/SiH1x1-chi2-xxy
+    yxx: sample-data/SiH1x1-chi2-yxx
+    yyy: sample-data/SiH1x1-chi2-yyy
+    yzz: sample-data/SiH1x1-chi2-yzz
+    yyz: sample-data/SiH1x1-chi2-yyz
+    yxz: sample-data/SiH1x1-chi2-yxz
+    yxy: sample-data/SiH1x1-chi2-yxy
+    zxx: sample-data/SiH1x1-chi2-zxx
+    zyy: sample-data/SiH1x1-chi2-zyy
+    zzz: sample-data/SiH1x1-chi2-zzz
+    zyz: sample-data/SiH1x1-chi2-zyz
+    zxz: sample-data/SiH1x1-chi2-zxz
+    zxy: sample-data/SiH1x1-chi2-zxy
+
+# Path to desired output file
+output:     sample.out

paper/paper.bib

@@ -6,6 +6,18 @@ @Misc{shgyieldrepo
   url    = {https://github.com/roguephysicist/SHGYield},
 }
 
+@article{andersonFMATS17,
+  author  = {Anderson, S. M. and Mendoza, B. S.},
+  title   = {Depth-Dependent Three-Layer Model for the Surface Second-Harmonic Generation Yield},
+  journal = {Frontiers in Materials},
+  volume  = {4},
+  pages   = {12},
+  year    = {2017},   
+  url     = {http://journal.frontiersin.org/article/10.3389/fmats.2017.00012},
+  issn    = {2296-8016},
+  doi     = {10.3389/fmats.2017.00012}
+}
+
 @article{andersonPRB16b,
   doi = {10.1103/PhysRevB.94.115314},
   year  = {2016},