Version 1.4.0
@@ -26,9 +31,30 @@
The interaction between laser and matter has attracted widespread interest since the invention of laser technology decades ago.
To study the interaction between an ultrafast and intense laser pulse and atoms/molecules, where the electrons are ionized from the targets through multi-photon or tunneling/over-barrier processes, a time-dependent Schrödinger equation (TDSE) simulation is usually required to be carried out.
-However, its high demand in computational resources and limited application scope (atoms and simple molecules) prevent it from its extensive application.
+However, its high demand in computational resources and limited application scope (atoms and simple molecules) prevents it from its extensive application.
+
+To overcome the shortcomings of TDSE, Corkum *et al.* [^Corkum_1989] proposed a scheme, where the electron is first ionized from the target through the tunneling mechanism, and then acts as a classical electron in the laser field.
+This scheme was further developed by Hu *et al.* [^Hu_1997], in which the initial conditions of the classical electrons and the Coulomb potential of the parent ion are more appropriatedly taken account.
+This scheme is named after the *Classical Trajectory Monte-Carlo (CTMC)* method, which has been widely adopted for research in interaction between high-intensity ultra-fast laser pulses and atoms/molecules.
+Compared with TDSE, trajectory simulation schemes including CTMC and its variants, are less demanding in computational resources, which, in addition, provides a clear physical picture of strong-field ionization.
+
+The essence of the trajectory simulation scheme lies in two aspects:
+(1) The initial conditions of the classical electron samples at the beginning of the classical trajectories, which consists of initial position $\vec{r}_0$ (i.e., the tunneling exit position), initial momenta $\vec{p}_0$, and the corresponding ionization probability $W$ carried by the electron sample.
+(2) The quantum phase property of classical trajectories, while the full classical trajectory (i.e., the CTMC) is widely adopted, there are schemes (e.g., QTMC and SCTS, which would be discussed further in the documentation) which introduce quantum phases in the electron trajectories and develop a semiclassical method for trajectory simulations.
+
+After decades of accumulation of research and development, the trajectory simulation has grown to a complete solution of research on strong-field ionization of atoms and molecules. Developing a library with implementation of existing methods, efficiency of calculation, extensibility for future development and ease of maintenance would provide great convenience for theoretical research on strong-field ionization. With such aim, here we present *SemiclassicalSFI.jl*, a program package written in julia language, which provides a general, efficient and out-of-box solution of performing trajectory simulations.
+
+[^Corkum_1989]: P. B. Corkum *et al.*, Above-Threshold Ionization in the Long-Wavelength Limit. *Phys. Rev. Lett.* **62**(11), 1259–1262 (1989). DOI: [10.1103/PhysRevLett.62.1259](https://dx.doi.org/10.1103/PhysRevLett.62.1259)
+
+[^Hu_1997]: B. Hu *et al.*, Plateau in Above-Threshold-Ionization Spectra and Chaotic Behavior in Rescattering Processes. *Phys. Lett. A* **236**, 533–542 (1997). DOI: [10.1016/S0375-9601(97)00811-6](https://dx.doi.org/10.1016/S0375-9601(97)00811-6)
+
+---------------------------
-As an alternative, semiclassical/classical electron trajectory simulation is widely used in numerical simulation in studies of strong-field ionization because it is less demanding in computational resources, which, in addition, provides a clear physical picture of strong-field ionization. This library written in julia aims to provide a general, efficient and out-of-box solution to perform trajectory simulations.
+## Features
+
+- *Versatile* : *SemiclassicalSFI.jl* supports a wide range of functions. As for initial conditions (rate method), the library supports (for atoms) *ADK*, *SFA* and *SFA-AE*, (for molecules) *MOADK* and *WFAT*. As for the trajectory phase method, the library supports *CTMC*, *QTMC* and *SCTS*. Non-dipole effects can also be included during the trajectory simulation.
+- *Out-of-box* : The usage of *SemiclassicalSFI.jl* is simple and straightforward.
+- *Extensible* : *SemiclassicalSFI.jl* has a well-defined structure, which makes it easy to include new features.
---------------------------
@@ -36,36 +62,52 @@ As an alternative, semiclassical/classical electron trajectory simulation is wid
### Prerequisites
-Some prerequisites are listed below:
-- [ **Hardware** ] A supported graphic card (optional, if you need to use GPU acceleration)
-- [ **Platform** ] *Linux* & Windows (Linux is suggested because `PySCFMolecularCalculator` which is used to calculate molecular structure factors is incompatible with Windows)
-- [ **Environment** ] Julia 1.7 & Python 3 (Python 3 is optional, but is required if you need to calculate the molecular structure factors)
+- *Minimum prerequisites* : Julia ≥1.7
+
+- *GPU acceleration of traj. simulation* : a supported graphic card (NVIDIA)
+
+- *MOADK and WFAT features* : Linux or macOS platform, Python 3 with the [PySCF](https://github.com/pyscf/pyscf) python package installed and the [PyCall.jl](https://github.com/JuliaPy/PyCall.jl) package successfully built.
+
+### Installing the package
-### Installing dependency PySCF (optional)
+This package is currently not in julia's general registry, but can be added through the repository URL:
-Currently the calculation of molecular structure factors relies on the [PySCF](https://github.com/pyscf/pyscf) python package, which only supports the Linux platform. This library calls the PySCF using the julia [PyCall.jl](https://github.com/JuliaPy/PyCall.jl) package, and there are two ways to set up the Python environment used by PyCall, here we suggest using your local Python environment for convenience:
+```julia
+using Pkg
+Pkg.add(url="https://github.com/TheStarAlight/SemiclassicalSFI.jl.git")
+# In pkg mode of REPL:
+# (@v1.8) pkg> add https://github.com/TheStarAlight/SemiclassicalSFI.jl.git
```
-$ julia
-julia> ENV["PYTHON"] = ""
-(@v1.7) pkg> add PyCall # if you haven't installed it yet
-(@v1.7) pkg> build PyCall
+
+It is suggested to test the package to check if the functions check if some special features (e.g., GPU acceleration and molecular calculation) work on your platform:
+
+```julia
+Pkg.test("SemiclassicalSFI")
+# In pkg mode of REPL:
+# (@v1.8) pkg> test SemiclassicalSFI
```
-Remember to install PySCF in your python !!
-The alternative choice is to install a private Python environment for julia via [Conda.jl](https://github.com/Luthaf/Conda.jl), for more information, please refer to [Conda.jl](https://github.com/JuliaPy/Conda.jl) and [PyCall.jl](https://github.com/JuliaPy/PyCall.jl).
+### Configuring Python and PySCF
-### Installing the package
+Currently the MO-ADK and WFAT features related to molecules rely on the [PySCF](https://github.com/pyscf/pyscf) python package. *SemiclassicalSFI.jl* calls the PySCF using the [PyCall.jl](https://github.com/JuliaPy/PyCall.jl) package. There are two ways to set up the Python environment used by PyCall, here we suggest using your local Python environment for convenience.
+
+To correctly set up the configuration of PyCall, first, set the `PYTHON` environment variable to your Python executable, and build the PyCall package:
+
+```julia
+ENV["PYTHON"] = "path/to/python_exec"
+using Pkg
+Pkg.build("PyCall")
+```
+
+And don't forget to install PySCF in your Python via pip:
-The package is still under development, to install the package, it's better to use the dev mode in the julia pkg manager (\ refers to a location where you want to store the source code, \ refers to the branch or tag version of the library you want to install):
```
-$ cd
-$ git clone https://github.com/TheStarAlight/SemiclassicalSFI.jl.git
-$ cd ./SemiclassicalSFI.jl
-$ git switch
-$ julia
-(@v1.7) pkg> dev .
+$ pip3 install pyscf
```
-When you want to update or use another version, switch to the version you want using `git`, and "`dev`" it in julia's pkg manager.
+
+
+**Note**: Since the PySCF does not support the Windows, the molecular calculation must be performed on a Linux or macOS platform.
+However, for Windows users, they may install the [WSL](https://learn.microsoft.com/en-us/windows/wsl/install) (Windows Subsystem for Linux), which supports the PySCF.
---------------------------
@@ -74,48 +116,50 @@ When you want to update or use another version, switch to the version you want u
The usage of the program is simple and straightforward.
The main function of the program is wrapped in a method `performSFI`, with all simulation parameters provided, the program would automatically generate initial electron samples and perform trajectory simulations, collecting the electrons and mapping them to the momentum spectrum, and finally save the spectrum to the specified location.
-The following are parameters of the method `performSFI` (the more detailed documentation is coming soon):
-
-### Required params. for all methods
-- `ionRateMethod = <:ADK|:SFA|:SFA_AE|:WFAT>` : Method of determining ionization rate. Currently supports `:ADK`, `:SFA`, `:SFA_AE` for atoms and `:WFAT` for molecules.
-- `laser::Laser` : Parameters of the laser field.
-- `target::Target` : Parameters of the target.
-- `sample_tSpan = (start,stop)` : Time span in which electrons are sampled.
-- `sample_tSampleNum` : Number of time samples.
-- `simu_tFinal` : Time when every trajectory simulation ends.
-- `finalMomentum_pMax = (pxMax,pyMax,pzMax)` : Boundaries of final momentum spectrum collecting in three dimensions.
-- `finalMomentum_pNum = (pxNum,pyNum,pzNum)` : Numbers of final momentum spectrum collecting in three dimensions.
-
-### Required params. for step-sampling methods
-- `ss_kdMax` : Boundary of kd (momentum's component along transverse direction (in xy plane)) samples.
-- `ss_kdNum` : Number of kd (momentum's component along transverse direction (in xy plane)) samples.
-- `ss_kzMax` : Boundary of kz (momentum's component along propagation direction (z ax.)) samples.
-- `ss_kzNum` : Number of kz (momentum's component along propagation direction (z ax.)) samples.
-
-### Required params. for Monte-Carlo-sampling methods
-- `mc_tBatchSize` : Number of electron samples in a single time sample.
-- `mc_ktMax` : Maximum value of momentum's transversal component (perpendicular to field direction).
-
-### Optional params. for all methods
-- `save_fileName` : Output HDF5 file name.
-- `save_3D_momentumSpec = false` : Determines whether 3D momentum spectrum is saved.
-- `simu_phaseMethod = <:CTMC|:QTMC|:SCTS>` : Method of classical trajectories' phase.
-- `simu_relTol = 1e-6` : Relative error tolerance when solving classical trajectories.
-- `simu_nondipole = false` : Determines whether non-dipole effect is taken account in the simulation.
-- `simu_GPU = false` : Determines whether GPU acceleration in trajectory simulation is used.
-- `rate_monteCarlo = false` : Determines whether Monte-Carlo sampling is used when generating electron samples.
-- `rate_ionRatePrefix = <:ExpRate|:ExpPre|:ExpJac|:Full>` : Prefix of the exponential term in the ionization rate. For MOADK & WFAT, `:ExpRate` & `:ExpPre` are the same.
-- `rydberg_collect = false` : Determines whether rydberg final states are collected.
-- `rydberg_prinQNMax` : Maximum principle quantum number n to be collected.
-
-### Optional params. for target `Molecule`
-- `mol_ionOrbitRelHOMO = 0` : Index of the ionizing orbit relative to the HOMO (e.g., 0 indicates HOMO, and -1 indicates HOMO-1) (default 0).
-
-### Optional params. for MOADK method:
-- `moadk_ionOrbit_m = 0` : Magnetic quantum number m of the ionizing orbital along the z axis. m = 0,1,2 indicate σ, π and δ respectively.
-
-### Optional params. for ADK method
-- `adk_ADKTunExit = <:IpF|:FDM|:Para>` : Tunneling exit method for ADK methods (when `ionRateMethod==:ADK`).
+The following are parameters of the method `performSFI` (for more detailed documentation, see [here](https://TheStarAlight.github.io/SemiclassicalSFI.jl/stable/manual3_main_method/)):
+
+### Required params. for all methods:
+- `init_cond_method = <:ADK|:SFA|:SFAAE|:WFAT|:MOADK>` : Method of electrons' initial conditions. Currently supports `:ADK`, `:SFA`, `:SFAAE` for atoms and `:WFAT`, `:MOADK` for molecules.
+- `laser::Laser` : A `Lasers.Laser` object containing information of the laser field.
+- `target::Target` : A `Targets.Target` object containing information of the atom/molecule target.
+- `sample_t_intv = (start,stop)` : Time interval in which the initial electrons are sampled.
+- `sample_t_num` : Number of time samples.
+- `traj_t_final` : Time when every trajectory simulation ends.
+- `final_p_max = (pxMax,pyMax,pzMax)` : Boundaries of final momentum spectrum collected in three dimensions.
+- `final_p_num = (pxNum,pyNum,pzNum)` : Numbers of final momentum spectrum collected in three dimensions.
+
+### Required params. for step-sampling methods:
+- `ss_kd_max` : Boundary of kd (momentum's component along transverse direction (in xy plane)) samples.
+- `ss_kd_num` : Number of kd (momentum's component along transverse direction (in xy plane)) samples.
+- `ss_kz_max` : Boundary of kz (momentum's component along propagation direction (z ax.)) samples.
+- `ss_kz_num` : Number of kz (momentum's component along propagation direction (z ax.)) samples.
+
+### Required params. for Monte-Carlo-sampling methods:
+- `mc_kp_num` : Number of kp (initial momentum which is perpendicular to field direction, two dimensional) samples in a single time sample.
+- `mc_kp_max` : Maximum value of momentum's transversal component (perpendicular to field direction).
+
+### Optional params. for all methods:
+- `save_path` : Output HDF5 file path.
+- `save_3D_spec = false` : Determines whether to save the 3D momentum spectrum (otherwise 2D) (default `false`).
+- `traj_phase_method = <:CTMC|:QTMC|:SCTS>` : Method of classical trajectories' phase (default `:CTMC`). Currently `:QTMC` and `:SCTS` only support atomic cases.
+- `traj_rtol = 1e-6` : Relative error tolerance when solving classical trajectories using adaptive methods (default `1e-6`).
+- `traj_nondipole = false` : Determines whether the non-dipole effect is taken account in the simulation (default `false`).
+- `traj_GPU = false` : [Experimental] Determines whether to enable GPU acceleration in trajectory simulation (default `false`).
+- `sample_monte_carlo = false` : Determines whether Monte-Carlo sampling is used when generating electron samples (default `false`). Currently only supports ADK.
+- `final_ryd_collect = false` : Determines whether the rydberg final states are collected (default `false`).
+- `final_ryd_n_max` : Determines the maximum principle quantum number n for rydberg final states to be collected.
+
+### Optional params. for atomic SFA, SFA-AE and ADK methods:
+- `rate_prefix = <:ExpRate|:ExpPre|:ExpJac|:Full>` : Prefix of the exponential term in the ionization rate (default `:ExpRate`).
+
+### Optional params. for target `Molecule`:
+- `mol_orbit_idx = 0` : Index of the ionizing orbit relative to the HOMO (e.g., `0` indicates HOMO, and `-1` indicates HOMO-1) (default `0`).
+
+### Optional params. for MO-ADK method:
+- `moadk_orbit_m = 0` : Magnetic quantum number m of the ionizing orbital along the z axis. m = 0,1,2 indicate σ, π and δ respectively (default `0`).
+
+### Optional params. for ADK method:
+- `adk_tun_exit = <:IpF|:FDM|:Para>` : Tunneling exit method for ADK methods (when `init_cond_method==:ADK`) (default `:IpF`).
---------------------------
@@ -132,35 +176,35 @@ filename = "SCSFI_HydLike_Ip_0.5662.h5" # output file name
# The target is a Hydrogen-like atom with ionization potential of 0.5662 a.u. and single nuclear charge.
t = SemiclassicalSFI.Targets.HydrogenLikeAtom(Ip=0.5662, Z=1)
-# The laser is a 800nm (NIR) 6-cycle circularly polarized pulse with a cos^4-shaped envelope (propagating in z direction).
-l = SemiclassicalSFI.Lasers.Cos4Laser(peakInt=3e14, waveLen=800, cycNum=6, ellip=1.0)
+# The laser is a 800nm (NIR) 2-cycle circularly polarized pulse with a cos^4-shaped envelope (propagating in z direction).
+l = SemiclassicalSFI.Lasers.Cos4Laser(peak_int=4e14, wave_len=800.0, cyc_num=2, ellip=1.0)
# Invokes the main method.
SemiclassicalSFI.performSFI(
# Using ADK as the rate method to give initial conditions of electron samples.
- ionRateMethod = :ADK,
+ init_cond_method = :ADK,
laser = l,
target = t,
- # Electrons ejected from the target between time range -300 and 300 a.u. would be sampled, and 20000 time samples would be picked at equal distances in the above time range.
- sample_tSpan = (-300,300),
- sample_tSampleNum = 20000,
- # Trajectory simulation of the ejected electrons would end at t=400 a.u. (a short span after the laser ends)
- simu_tFinal = 400,
- # Parameters of the final momentum spectrum. Electrons with momentum within (±3,±3,±3) a.u. would be collected, and mapped to a 800×800×1 cartesian grid (grid size in z direction is 1 and save_3D is false, which indicates that the simulation only saves 2D momentum spectrum).
- finalMomentum_pMax = (3,3,3),
- finalMomentum_pNum = (800,800,1),
- save_3D_momentumSpec = false,
+ # Electrons ejected from the target between time range -80 and 80 a.u. would be sampled, and 5000 time samples would be picked at equal distances in the above time range.
+ sample_t_intv = (-80,80),
+ sample_t_num = 5000,
+ # Trajectory simulation of the ejected electrons would end at t=120 a.u. (a short span after the laser ends)
+ traj_t_final = 120,
+ # Parameters of the final momentum spectrum. Electrons with momentum within (±2,±2,±2) a.u. would be collected, and mapped to a 800×800×1 cartesian grid (grid size in z direction is 1 and save_3D is false, which indicates that the simulation only saves 2D momentum spectrum).
+ final_p_max = (2,2,2),
+ final_p_num = (500,500,1),
+ save_3D_spec = false,
# Step-sampling parameters of electrons' initial momentum. `kd` refers to the momentum's component along transverse direction (in xy plane), while `kz` refers to the momentum's component along propagation direction (z ax.).
- ss_kdMax = 2.,
- ss_kdNum = 500,
- ss_kzMax = 2.,
- ss_kzNum = 150,
+ ss_kd_max = 2.,
+ ss_kd_num = 500,
+ ss_kz_max = 2.,
+ ss_kz_num = 150,
# Classical Trajectory Monte-Carlo (CTMC), which indicates no account of phase.
- simu_phaseMethod = :CTMC,
+ traj_phase_method = :CTMC,
# GPU acceleration in trajectory simulation is enabled.
- simu_GPU = true,
+ traj_GPU = true,
# output file name.
- save_fileName = filename
+ save_path = filename
)
```
@@ -168,28 +212,26 @@ SemiclassicalSFI.performSFI(
## Troubleshooting
-### Precompilation failure
-Sometimes the precompilation of the package and its dependencies fails, which usually happens on SciML's packages, while no action in the pkg manager works.
+### Possible solution to precompilation failure
-Under such circumstances, try to delete the compiled julia code (usually stored in ~/.julia/compiled/\) and precompile again.
-
-If the problem still exists after precompiling from scratch, you may try switching the SciML dependencies' versions in the julia.
-As for the author, *OrdinaryDiffEq@6.51* and *DiffEqGPU@1.26* runs well on Windows 10, while for OrdinaryDiffEq@6.20/37/41, the precompilation never succeeded.
-
-### Runtime warning: electrons with anomalous momentum
-Sometimes (or usually, for some unlucky users) you may encounter such warning during the simulation if you use GPU acceleration:
-```
-[Ensemble Simulation] Found electron (#<...> in the batch) with anomalous momentum <...>.
+Sometimes the precompilation of the package and its dependencies fails, which usually happens on SciML's packages.
+Under such circumstances, try to delete the compiled julia code (usually stored in `~/.julia/compiled/`) and precompile again.
+If the problem still exists after precompiling from scratch, you may try switching the dependencies' versions in the julia, which is done by specifying the version when adding the packages:
+```julia
+using Pkg
+Pkg.add(name="package_name", version="1.0")
+# In pkg mode of REPL:
+# (@v1.8) pkg> add package_name@1.0
```
-This is possibly resulted from the DiffEqGPU package for some unknown reasons.
-To solve this problem, try switching to another DiffEqGPU version (v1.26 is suggested).
+
+It is shown that *OrdinaryDiffEq@6.53.3* and *DiffEqGPU@2.4.1* runs well on Windows 10 (10.0.19044) and WSL Ubuntu (22.04.1 LTS).
---------------------------
## Contributors
- [Mingyu Zhu](https://github.com/TheStarAlight) @ ECNU
-- Hongcheng Ni @ ECNU
+- [Hongcheng Ni](https://faculty.ecnu.edu.cn/_s29/nhc_en/main.psp) @ ECNU
---------------------------
diff --git a/docs/make.jl b/docs/make.jl
new file mode 100644
index 0000000..640eae2
--- /dev/null
+++ b/docs/make.jl
@@ -0,0 +1,26 @@
+using Documenter, SemiclassicalSFI
+
+makedocs(
+ sitename = "SemiclassicalSFI.jl",
+ pages = [
+ "Home" => "index.md",
+ "Theory" => [
+ "Initial Conditions" => "theory1_initial_conditions.md",
+ "Trajectory Simulation and Phase Methods" => "theory2_trajectory_simulation_phase_methods.md"
+ ],
+ "Manual" => [
+ "Targets" => "manual1_targets.md",
+ "Lasers" => "manual2_lasers.md",
+ "Main Method" => "manual3_main_method.md"
+ ],
+ "Examples" => [
+ "Attoclock and Initial Condition Methods" => "example1_attoclock.md",
+ "Phase Methods" => "example2_phase_methods.md"
+ ]
+ ],
+ format = Documenter.HTML(
+ ansicolor=true,
+ edit_link=nothing,
+ footer="· *SemiclassicalSFI.jl* Documentation · by *Mingyu Zhu* and other contributors"
+ )
+ )
\ No newline at end of file
diff --git a/docs/run_doc_server.jl b/docs/run_doc_server.jl
new file mode 100644
index 0000000..d9e8839
--- /dev/null
+++ b/docs/run_doc_server.jl
@@ -0,0 +1,4 @@
+# run this script to start a local web server that hosts the documentation.
+
+using LiveServer
+serve(dir="./docs/build")
\ No newline at end of file
diff --git a/docs/src/assets/example_attoclock.png b/docs/src/assets/example_attoclock.png
new file mode 100644
index 0000000..579dcf2
Binary files /dev/null and b/docs/src/assets/example_attoclock.png differ
diff --git a/docs/src/assets/example_phase_methods.png b/docs/src/assets/example_phase_methods.png
new file mode 100644
index 0000000..0713b9d
Binary files /dev/null and b/docs/src/assets/example_phase_methods.png differ
diff --git a/docs/src/assets/manual2_lasers_azimuth.svg b/docs/src/assets/manual2_lasers_azimuth.svg
new file mode 100644
index 0000000..38dd14c
--- /dev/null
+++ b/docs/src/assets/manual2_lasers_azimuth.svg
@@ -0,0 +1,3674 @@
+
+
+
diff --git a/docs/src/assets/manual2_lasers_cep.svg b/docs/src/assets/manual2_lasers_cep.svg
new file mode 100644
index 0000000..ea5598b
--- /dev/null
+++ b/docs/src/assets/manual2_lasers_cep.svg
@@ -0,0 +1,5072 @@
+
+
+
diff --git a/docs/src/assets/manual2_lasers_ellip.svg b/docs/src/assets/manual2_lasers_ellip.svg
new file mode 100644
index 0000000..8f29a13
--- /dev/null
+++ b/docs/src/assets/manual2_lasers_ellip.svg
@@ -0,0 +1,5050 @@
+
+
+
diff --git a/docs/src/assets/manual2_lasers_env_cos2.svg b/docs/src/assets/manual2_lasers_env_cos2.svg
new file mode 100644
index 0000000..2dd71c6
--- /dev/null
+++ b/docs/src/assets/manual2_lasers_env_cos2.svg
@@ -0,0 +1,6018 @@
+
+
+
diff --git a/docs/src/assets/manual2_lasers_env_cos4.svg b/docs/src/assets/manual2_lasers_env_cos4.svg
new file mode 100644
index 0000000..de528b2
--- /dev/null
+++ b/docs/src/assets/manual2_lasers_env_cos4.svg
@@ -0,0 +1,5971 @@
+
+
+
diff --git a/docs/src/assets/manual2_lasers_env_gaussian.svg b/docs/src/assets/manual2_lasers_env_gaussian.svg
new file mode 100644
index 0000000..86f6ac2
--- /dev/null
+++ b/docs/src/assets/manual2_lasers_env_gaussian.svg
@@ -0,0 +1,9237 @@
+
+
+
diff --git a/docs/src/assets/manual2_lasers_env_trapezoidal.svg b/docs/src/assets/manual2_lasers_env_trapezoidal.svg
new file mode 100644
index 0000000..035162c
--- /dev/null
+++ b/docs/src/assets/manual2_lasers_env_trapezoidal.svg
@@ -0,0 +1,5922 @@
+
+
+
diff --git a/docs/src/example1_attoclock.md b/docs/src/example1_attoclock.md
new file mode 100644
index 0000000..609a756
--- /dev/null
+++ b/docs/src/example1_attoclock.md
@@ -0,0 +1,46 @@
+# Example: Attoclock and Initial Condition Methods
+
+With the aim of studying the influence of non-adiabatic effects on the attoclock signal, we employ ADK, SFA-AE and SFA to provide the initial conditions of the electrons, and perform trajectory simulations.
+
+The simulation parameters are taken from Ma *et al.* [*J. Phys. B* **54**, 144001 (2021)] [^Ma_2021]
+
+[^Ma_2021]: Y. Ma *et al.*, Influence of Nonadiabatic, Nondipole and Quantum Effects on the Attoclock Signal. *J. Phys. B: At. Mol. Opt. Phys.* **54**, 144001 (2021). DOI: [10.1088/1361-6455/ac0d3e](https://dx.doi.org/10.1088/1361-6455/ac0d3e)
+
+```julia
+@info "SemiclassicalSFI Example - Attoclock"
+
+@info "Loading Packages..."
+using SemiclassicalSFI
+using SemiclassicalSFI.Targets
+using SemiclassicalSFI.Lasers
+
+t = HeAtom()
+l = Cos4Laser(peak_int=4e14, wave_len=800, cyc_num=2, ellip=1.0)
+
+init_cond_list = [:ADK, :SFAAE, :SFA]
+
+for init_cond in init_cond_list
+ @info "Running $(init_cond)..."
+ filename = "SCSFI_He_4e14_800nm_2cyc_CP_$(init_cond)_CTMC.h5"
+ performSFI(
+ target = t,
+ laser = l,
+ init_cond_method = init_cond,
+ sample_t_intv = (-80,80),
+ sample_t_num = 20000,
+ ss_kd_max = 1.5,
+ ss_kd_num = 1000,
+ ss_kz_max = 1.5,
+ ss_kz_num = 50,
+ traj_phase_method = :CTMC,
+ traj_t_final = 120,
+ final_p_max = (2,2,2),
+ final_p_num = (500,500,1),
+ save_path = filename
+ )
+end
+```
+
+Below shows the momentum spectrum squashed on the polarization plane.
+
+
diff --git a/docs/src/example2_phase_methods.md b/docs/src/example2_phase_methods.md
new file mode 100644
index 0000000..8eb5128
--- /dev/null
+++ b/docs/src/example2_phase_methods.md
@@ -0,0 +1,93 @@
+# Example: Phase Methods
+
+To investigate the difference between the phase methods QTMC and SCTS, we use an 8-cycle linearly polarized laser pulse, employing the QTMC and SCTS as the phase methods during the trajectory simulation.
+
+The simulation parameters are brought from [*Phys. Rev. A* **94**, 013415 (2016)] [^ShvetsovShilovski_2016]
+
+[^ShvetsovShilovski_2016]: N. I. Shvetsov-Shilovski *et al.* Semiclassical Two-Step Model for Strong-Field Ionization. *Phys. Rev. A* **94**, 013415 (2016). DOI: [10.1103/PhysRevA.94.013415](https://dx.doi.org/10.1103/PhysRevA.94.013415)
+
+## QTMC Code
+
+```julia
+@info "SemiclassicalSFI Example - Phase Methods - QTMC"
+
+@info "Loading Packages..."
+using SemiclassicalSFI
+using SemiclassicalSFI.Targets
+using SemiclassicalSFI.Lasers
+
+t = HAtom()
+l = Cos4Laser(peak_int=0.9e14, wave_len=800, cyc_num=8, ellip=0.0)
+
+filename = "SCSFI_H_0.9e14_800nm_8cyc_LP_ADKPara_ExpRate_QTMC.h5"
+
+performSFI(
+ target = t,
+ laser = l,
+ init_cond_method = :ADK,
+ adk_tun_exit = :Para,
+ sample_t_intv = (-300,300),
+ sample_t_num = 50000,
+ sample_cutoff_limit = 1e-16,
+ ss_kd_max = 1.0,
+ ss_kd_num = 400,
+ ss_kz_max = 1.0,
+ ss_kz_num = 400,
+ rate_prefix = :ExpRate,
+ traj_phase_method = :QTMC,
+ traj_t_final = 450,
+ traj_rtol = 1e-6,
+ save_3D_spec = true,
+ final_p_max = (1,1,1),
+ final_p_num = (400,400,400),
+ save_path = filename
+ )
+```
+
+## SCTS Code
+
+```julia
+@info "SemiclassicalSFI Example - Phase Methods - SCTS"
+
+@info "Loading Packages..."
+using SemiclassicalSFI
+using SemiclassicalSFI.Targets
+using SemiclassicalSFI.Lasers
+
+t = HAtom()
+l = Cos4Laser(peak_int=0.9e14, wave_len=800, cyc_num=8, ellip=0.0)
+
+filename = "SCSFI_H_0.9e14_800nm_8cyc_LP_ADKPara_ExpRate_SCTS.h5"
+
+performSFI(
+ target = t,
+ laser = l,
+ init_cond_method = :ADK,
+ adk_tun_exit = :Para,
+ sample_t_intv = (-300,300),
+ sample_t_num = 50000,
+ sample_cutoff_limit = 1e-16,
+ ss_kd_max = 1.0,
+ ss_kd_num = 400,
+ ss_kz_max = 1.0,
+ ss_kz_num = 400,
+ rate_prefix = :ExpRate,
+ traj_phase_method = :SCTS,
+ traj_t_final = 450,
+ traj_rtol = 1e-6,
+ save_3D_spec = true,
+ final_p_max = (1,1,1),
+ final_p_num = (400,400,400),
+ save_path = filename
+ )
+```
+
+Each piece of code takes 8.5 hours using 4 threads on an AMD Ryzen 9 7950X CPU at Manjaro Linux (Uranos 23.0.0).
+
+!!! note "Note"
+ - This example takes a very long time (~10h each task for 4 threads) to finish.
+ - In this example, the 3D momentum spectrum is collected, which is very memory-consuming, consider using fewer threads.
+ - A 16-thread task is much less efficient than a 4-thread task and consumes much more memory. Therefore, it's suggested to use fewer threads (e.g., 4) per task.
+
+The final momentum spectrum (the slice on the ``p_z-p_x`` plane near ``p_y=0``) is displayed below.
+
\ No newline at end of file
diff --git a/docs/src/index.md b/docs/src/index.md
new file mode 100644
index 0000000..e530d0f
--- /dev/null
+++ b/docs/src/index.md
@@ -0,0 +1,109 @@
+# 🎆SemiclassicalSFI.jl
+
+*Implementation of classical/semiclassical methods in strong-field ionization of atoms and molecules.*
+
+```@contents
+Pages = ["index.md"]
+```
+
+## Background
+
+The interaction between laser and matter has attracted widespread interest since the invention of laser technology decades ago.
+To study the interaction between an ultrafast and intense laser pulse and atoms/molecules, where the electrons are ionized from the targets through multi-photon or tunneling/over-barrier processes, a time-dependent Schrödinger equation (TDSE) simulation is usually required to be carried out.
+However, its high demand in computational resources and limited application scope (atoms and simple molecules) prevents it from its extensive application.
+
+To overcome the shortcomings of TDSE, Corkum *et al.* [^Corkum_1989] proposed a scheme, where the electron is first ionized from the target through the tunneling mechanism, and then acts as a classical electron in the laser field.
+This scheme was further developed by Hu *et al.* [^Hu_1997], in which the initial conditions of the classical electrons and the Coulomb potential of the parent ion are more appropriatedly taken account.
+This scheme is named after the *Classical Trajectory Monte-Carlo (CTMC)* method, which has been widely adopted for research in interaction between high-intensity ultra-fast laser pulses and atoms/molecules.
+Compared with TDSE, trajectory simulation schemes including CTMC and its variants, are less demanding in computational resources, which, in addition, provides a clear physical picture of strong-field ionization.
+
+The essence of the trajectory simulation scheme lies in two aspects:
+(1) The initial conditions of the classical electron samples at the beginning of the classical trajectories, which consists of initial position $\bm{r}_0$ (i.e., the tunneling exit position), initial momenta $\bm{p}_0$, and the corresponding ionization probability $W$ carried by the electron sample.
+(2) The quantum phase property of classical trajectories, while the full classical trajectory (i.e., the CTMC) is widely adopted, there are schemes (e.g., QTMC and SCTS, which would be discussed further in the documentation) which introduce quantum phases in the electron trajectories and develop a semiclassical method for trajectory simulations.
+
+After decades of accumulation of research and development, the trajectory simulation has grown to a complete solution of research on strong-field ionization of atoms and molecules. Developing a library with implementation of existing methods, efficiency of calculation, extensibility for future development and ease of maintenance would provide great convenience for theoretical research on strong-field ionization. With such aim, here we present *SemiclassicalSFI.jl*, a program package written in julia language, which provides a general, efficient and out-of-box solution of performing trajectory simulations.
+
+[^Corkum_1989]: P. B. Corkum *et al.*, Above-Threshold Ionization in the Long-Wavelength Limit. *Phys. Rev. Lett.* **62**(11), 1259–1262 (1989). DOI: [10.1103/PhysRevLett.62.1259](https://dx.doi.org/10.1103/PhysRevLett.62.1259)
+
+[^Hu_1997]: B. Hu *et al.*, Plateau in Above-Threshold-Ionization Spectra and Chaotic Behavior in Rescattering Processes. *Phys. Lett. A* **236**, 533–542 (1997). DOI: [10.1016/S0375-9601(97)00811-6](https://dx.doi.org/10.1016/S0375-9601(97)00811-6)
+
+## Features
+
+- *Versatile* : *SemiclassicalSFI.jl* supports a wide range of functions. As for initial conditions (rate method), the library supports (for atoms) *ADK*, *SFA* and *SFA-AE*, (for molecules) *MOADK* and *WFAT*. As for the trajectory phase method, the library supports *CTMC*, *QTMC* and *SCTS*. Non-dipole effects can also be included during the trajectory simulation.
+- *Out-of-box* : The usage of *SemiclassicalSFI.jl* is simple and straightforward.
+- *Extensible* : *SemiclassicalSFI.jl* has a well-defined structure, which makes it easy to include new features.
+
+## Installation
+
+### Prerequisites
+
+- *Minimum prerequisites* : Julia ≥1.7
+
+- *GPU acceleration of traj. simulation* : a supported graphic card (NVIDIA)
+
+- *MOADK and WFAT features* : Linux or macOS platform, Python 3 with the [PySCF](https://github.com/pyscf/pyscf) python package installed and the [PyCall.jl](https://github.com/JuliaPy/PyCall.jl) package successfully built.
+
+### Installing the package
+
+This package is currently not in julia's general registry, but can be added through the repository URL:
+
+```julia
+using Pkg
+Pkg.add(url="https://github.com/TheStarAlight/SemiclassicalSFI.jl.git")
+# In pkg mode of REPL:
+# (@v1.8) pkg> add https://github.com/TheStarAlight/SemiclassicalSFI.jl.git
+```
+
+It is suggested to test the package to check if the functions check if some special features (e.g., GPU acceleration and molecular calculation) work on your platform:
+
+```julia
+Pkg.test("SemiclassicalSFI")
+# In pkg mode of REPL:
+# (@v1.8) pkg> test SemiclassicalSFI
+```
+
+!!! note "Possible solution to precompilation failure"
+
+ Sometimes the precompilation of the package and its dependencies fails, which usually happens on SciML's packages.
+ Under such circumstances, try to delete the compiled julia code (usually stored in `~/.julia/compiled/`) and precompile again.
+ If the problem still exists after precompiling from scratch, you may try switching the dependencies' versions in the julia, which is done by specifying the version when adding the packages:
+ ```julia
+ using Pkg
+ Pkg.add(name="package_name", version="1.0")
+ # In pkg mode of REPL:
+ # (@v1.8) pkg> add package_name@1.0
+ ```
+
+ It is shown that *OrdinaryDiffEq@6.53.3* and *DiffEqGPU@2.4.1* runs well on Windows 10 (10.0.19044), WSL Ubuntu (22.04.1 LTS) and Manjaro Linux (Uranos 23.0.0).
+
+### Configuring Python and PySCF
+
+Currently the [MO-ADK](@ref MOADK) and [WFAT](@ref WFAT) features related to molecules rely on the [PySCF](https://github.com/pyscf/pyscf) python package. *SemiclassicalSFI.jl* calls the PySCF using the [PyCall.jl](https://github.com/JuliaPy/PyCall.jl) package. There are two ways to set up the Python environment used by PyCall, here we suggest using your local Python environment for convenience.
+
+To correctly set up the configuration of PyCall, first, set the `PYTHON` environment variable to your Python executable, and build the PyCall package:
+
+```julia
+ENV["PYTHON"] = "path/to/python_exec"
+using Pkg
+Pkg.build("PyCall")
+```
+
+And don't forget to install PySCF in your Python via pip:
+
+```
+$ pip3 install pyscf
+```
+
+!!! tip "Note"
+ Since the PySCF does not support the Windows, the molecular calculation must be performed on a Linux or macOS platform.
+ However, for Windows users, they may install the [WSL](https://learn.microsoft.com/en-us/windows/wsl/install) (Windows Subsystem for Linux), which supports the PySCF.
+
+
+## Contributors
+
+- [Mingyu Zhu](https://github.com/TheStarAlight) @ ECNU
+- [Hongcheng Ni](https://faculty.ecnu.edu.cn/_s29/nhc_en/main.psp) @ ECNU
+
+## License
+
+This package is licensed under the Apache 2.0 license, and is copyrighted by Mingyu Zhu, Hongcheng Ni and the other contributors.
diff --git a/docs/src/manual1_targets.md b/docs/src/manual1_targets.md
new file mode 100644
index 0000000..4fe3e07
--- /dev/null
+++ b/docs/src/manual1_targets.md
@@ -0,0 +1,201 @@
+# [Targets](@id targets_doc)
+
+*This section provides information of available targets (atoms/molecules) in the library.*
+
+A target interacts with the laser field and release the electron through multi-photon or tunneling processes.
+Here we list available targets implemented in the [`Targets`](@ref) module of the library.
+
+```@docs
+Targets
+```
+
+```@contents
+Pages = ["manual1_targets.md"]
+Depth = 3
+```
+
+```@meta
+CurrentModule = SemiclassicalSFI.Targets
+```
+
+```@setup manual_targets
+using SemiclassicalSFI
+using SemiclassicalSFI.Targets
+```
+
+## List of Available Properties and Methods
+
+The available properties and methods of the targets are listed below.
+To obtain a property of the target, invoke the property as a method and pass the target object as an argument. The following shows an example.
+
+```@repl manual_targets
+t = HeAtom()
+IonPotential(t)
+V = TargetPotential(t)
+V(1.0,0.0,0.0)
+```
+
+### Atoms
+
+| |[`HydrogenLikeAtom`](@ref) | [`SAEAtom`](@ref) |
+|:------------------|:-:|:-:|
+|`IonPotential` | ✔ | ✔ |
+|`AsympNuclCharge` | ✔ | ✔ |
+|`SoftCore` | ✔ | |
+|`TargetName` | ✔ | ✔ |
+|`TargetPotential` | ✔ | ✔ |
+|`TargetForce` | ✔ | ✔ |
+
+### Molecule
+
+#### Properties
+
+| property name | params | example |
+|:------------------------------|:------------------------------|:----------------------------------|
+|`MolAtoms` | -- |`["H","H"]` |
+|`MolAtomCoords` | -- |`[0.0 0.0 -0.375; 0.0 0.0 0.375]` |
+|`MolCharge` | -- |`0` |
+|`MolEnergyDataAvailable` | -- |`true` |
+|`MolEnergyLevels` | -- |`[-0.590975, 0.174188, ...]` |
+|`MolHOMOIndex` | -- |`1` |
+|`MolHOMOEnergy` | -- |`-0.590975` |
+|`MolWFATAvailableIndices` | -- |`Set([0])` |
+|`MolWFATData` |`orbitIdx_relHOMO=0` | -- |
+|`MolWFATStructureFactor_G` |`orbitIdx_relHOMO=0,nξ,m,β,γ` |`1.91369 - 1.46476e-33im` |
+|`MolAsympCoeffAvailableIndices`| -- |`Set([0])` |
+|`MolAsympCoeff` |`orbitIdx_relHOMO=0` | -- |
+|`MolMOADKStructureFactor_B` |`orbitIdx_relHOMO=0,m_,β,γ` |`2.21496 + 0.0im` |
+|[`MolRotation`](@ref) | -- |`(0.0,0.0,0.0)` |
+|`IonPotential` |`orbitIdx_relHOMO=0` |`0.590975` |
+|`AsympNuclCharge` | -- |`1` |
+|`TargetName` | -- |`"Hydrogen"` |
+|`TargetPotential` | -- | -- |
+|`TargetForce` | -- | -- |
+
+#### Methods
+
+| method name | params |
+|:------------------------------|:--------------------------------------|
+|[`SetMolRotation`](@ref) |`α,β,γ` or `(α,β,γ)` |
+|`MolCalcEnergyData!` |`orbitIdx_relHOMO=0, MCType, kwargs...`|
+|[`MolCalcWFATData!`](@ref) |`orbitIdx_relHOMO=0, MCType, kwargs...`|
+|[`MolCalcAsympCoeff!`](@ref) |`orbitIdx_relHOMO=0, MCType, kwargs...`|
+
+## Hydrogen-Like Atom
+
+A hydrogen-like atom has a potential of the form
+```math
+V(r) = -\frac{Z}{\sqrt{r^2+a}},
+```
+where ``Z`` is the nuclear charge number;
+``a`` denotes the soft-core parameter, which is applied to avoid singularity of the potential and can be adjusted to fit the actual ionization potential of the atom (obtained by TDSE).
+
+The hydrogen-like atom is implemented in the library as [`HydrogenLikeAtom`](@ref).
+```@docs
+Targets.HydrogenLikeAtom
+```
+
+
+## Single-Active-Electron (SAE) Atom
+
+The single-active-electron (SAE) atom is an implementation of the empirical atomic SAE model potential proposed by Tong *et al.* [^Tong_2005]
+The model potential of the SAE atom has the form
+```math
+V(r) = -\frac{Z + a_1 \mathrm{e}^{-b_1 r} + a_2 r \mathrm{e}^{-b_2 r} + a_3 \mathrm{e}^{-b_3 r}}{r},
+```
+where the ``a_i`` and ``b_i`` are tunable model potential parameters [^note].
+
+The SAE atom is implemented in the library as [`SAEAtom`](@ref).
+```@docs
+Targets.SAEAtom
+```
+
+[^Tong_2005]: X. M. Tong *et al.*, Empirical Formula for Static Field Ionization Rates of Atoms and Molecules by Lasers in the Barrier-Suppression Regime. *J. Phys. B: At. Mol. Opt. Phys.* **38**, 2593–2600. DOI: [10.1088/0953-4075/38/15/001](https://dx.doi.org/10.1088/0953-4075/38/15/001)
+[^note]: The symbols of the parameters are different from that in the original article. ``a_1, b_1, a_2, b_2, a_3, b_3`` correspond to ``a_1, a_2, a_3, a_4, a_5, a_6`` in the original article respectively.
+
+
+## Preset Atom Library
+
+The library provides some preset commonly-used atoms or atomic ions for convenience.
+
+- [`HydrogenLikeAtom`](@ref) : H, He⁺, Li²⁺
+
+- [`SAEAtom`](@ref) : He, Ne, Ne⁺, Ne²⁺, Ar, Ar⁺, Ar²⁺, V, Ni, Kr, Kr⁺, Rb, Nb, Pd, Xe, Xe⁺, Ta
+
+These atoms/ions can be obtained by invoking `Targets.**Atom()` (for neutral atoms) or `Targets.**#pAtom()` (for positive atomic ions), where `**` denotes the symbol of the nucleus and `#` denotes the positive charge the ion carries.
+
+Example:
+
+```@setup manual_targets
+using SemiclassicalSFI
+```
+```@repl manual_targets
+t1 = Targets.HAtom()
+t2 = Targets.Xe1pAtom()
+```
+
+
+## Molecule
+
+The `Molecule` object represents a generic molecule, which is implemented in the library as [`Molecule`](@ref).
+The structure of `Molecule` is much more complex than that of atoms because the [Molecular ADK (MO-ADK) theory](@ref MOADK) and [Weak-Field Asymptotic Theory (WFAT)](@ref WFAT) features for molecular strong-field ionization require a number of coefficients, which are saved to files for convenience.
+
+### Initialization, saving and loading
+
+The `Molecule` object can be initialized either by providing necessary information of the molecule (mainly atoms, coordinates of the atoms and the charge of the molecule) or from external data (stored in the HDF5 format), cf. the documentation of [`Molecule`](@ref):
+
+```@docs
+Targets.Molecule
+```
+
+The `Molecule` object, after modification, can be manually saved to a HDF5 file via [`MolSaveDataAs`](@ref).
+
+```@docs
+Targets.MolSaveDataAs
+```
+
+### Molecular-SFI Data Preparation
+
+To use the molecular strong-field ionization theories such as the MO-SFA, MO-ADK and WFAT to provide the intitial conditions of the electrons, the structure coefficients of the `Molecule` have to be calculated beforehand and stored in the object.
+Cf. the documentation of [`MolCalcAsympCoeff!`](@ref) and [`MolCalcWFATData!`](@ref).
+
+Evaluation of the structure coefficients depends on the external quantum chemistry packages.
+The [`Targets.MolecularCalculators`](@ref) module undertakes the task of communication with the external quantum packages.
+Currently only the [`PySCFMolecularCalculator`](@ref) is implemented.
+
+!!! tip "Customized calculation parameters"
+ When invoking `MolCalcAsympCoeff!` and `MolCalcWFATData!` to perform calculation of structure coefficients, customized calculation parameters can be passed to the `kwargs` of these methods.
+ These parameters would be passed to the constructor method of the `MolecularCalculator` (e.g., the `basis` parameter of the [`PySCFMolecularCalculator`](@ref)), as well as the [`MolecularCalculators.calc_asymp_coeff`](@ref), [`MolecularCalculators.calc_WFAT_data`](@ref) methods.
+ Refer to their documentation below for more information.
+
+```@docs
+Targets.MolCalcAsympCoeff!
+Targets.MolCalcWFATData!
+```
+
+```@docs
+Targets.MolecularCalculators
+Targets.MolecularCalculators.PySCFMolecularCalculator
+```
+
+```@docs
+Targets.MolecularCalculators.calc_asymp_coeff
+Targets.MolecularCalculators.calc_WFAT_data
+```
+
+### Molecule's Orientation
+
+The molecule's orientation is described by a set of Euler angles (``z-y'-z''`` convention), which defines a rotational transformation from the molecular frame (MF) to the lab frame (LF).
+This property of `Molecule` is NOT included in the saved file and thus needs to be specified each time upon initialization of the `Molecule` object from external files.
+
+!!! note "Note"
+ Here the three Euler angles `(α,β,γ)` that describe the `Molecule`'s orientation are completely different from that of the Euler angles `(β',γ')` in the [WFAT](@ref WFAT) and [MO-ADK](@ref MOADK) theory.
+ These theories' "lab frame" is chosen for convenience of theoretical formulation, where the electric field is assumed to be static, pointing towards the ``+z`` direction,
+ and has no relation with the lab frame mentioned above.
+
+The orientation of the molecule can be obtained and set via the [`MolRotation`](@ref) and [`SetMolRotation`](@ref) methods.
+
+```@docs
+Targets.MolRotation
+Targets.SetMolRotation
+```
diff --git a/docs/src/manual2_lasers.md b/docs/src/manual2_lasers.md
new file mode 100644
index 0000000..b4dfaa2
--- /dev/null
+++ b/docs/src/manual2_lasers.md
@@ -0,0 +1,208 @@
+# [Lasers](@id lasers_doc)
+
+*This section provides information of available lasers in the library.*
+
+In this section we list available lasers implemented in the [`Lasers`](@ref) module of the library.
+
+```@docs
+Lasers
+```
+
+```@contents
+Pages = ["manual2_lasers.md"]
+Depth = 3
+```
+
+```@meta
+CurrentModule = SemiclassicalSFI.Lasers
+```
+
+```@setup manual_lasers
+using SemiclassicalSFI
+using SemiclassicalSFI.Lasers
+```
+
+
+## Basic Properties of Monochromatic Lasers
+
+A monochromatic laser is composed of the carrier wave ``\cos{(\omega t+\phi)}`` and the envelope ``f_{\mathrm{env}}(t)``.
+Given the amplitude of the vector potential ``A_0``, the time-dependent vector potential of the laser,
+which we assume to propagate in ``z`` direction and have ``x`` axis as the principle axis of polarization, reads
+```math
+\bm{A}(t) = A_0 f_{\mathrm{env}}(t) \left[ \cos{(\omega t+\phi)} \bm{e}_x + \varepsilon \sin{(\omega t+\phi)} \bm{e}_y \right],
+```
+where ``\varepsilon`` is the ellipticity.
+
+### List of Available Properties
+
+Currently the monochromatic lasers implemented in the library include
+[`Cos4Laser`](@ref), [`Cos2Laser`](@ref), [`GaussianLaser`](@ref) and [`TrapezoidalLaser`](@ref).
+
+The available properties of the laser fields are listed below.
+To obtain a property of the laser field, invoke the property as a method and pass the laser object as an argument. The following shows an example.
+
+```@repl manual_lasers
+l = Cos4Laser(peak_int=1e14, wave_len=800.0, cyc_num=10, ellip=0)
+LaserA0(l)
+Ax = LaserAx(l)
+Ax(0.0)
+```
+
+| |[`Cos4Laser`](@ref) | [`Cos2Laser`](@ref) | [`GaussianLaser`](@ref) | [`TrapezoidalLaser`](@ref) |
+|:--------------|:-:|:-:|:-:|:-:|
+|`PeakInt` | ✔ | ✔ | ✔ | ✔ |
+|`WaveLen` | ✔ | ✔ | ✔ | ✔ |
+|`CycNum` | ✔ | ✔ | | |
+|`SpreadCycNum` | | | ✔ | |
+|`SpreadDuration`| | | ✔ | |
+|`FWHM_Duration`| | | ✔ | |
+|`CycNumTotal` | | | | ✔ |
+|`CycNumTurnOn` | | | | ✔ |
+|`CycNumTurnOff`| | | | ✔ |
+|`CycNumConst` | | | | ✔ |
+|`Ellipticity` | ✔ | ✔ | ✔ | ✔ |
+|`Azimuth` | ✔ | ✔ | ✔ | ✔ |
+|`AngFreq` | ✔ | ✔ | ✔ | ✔ |
+|`Period` | ✔ | ✔ | ✔ | ✔ |
+|`CEP` | ✔ | ✔ | ✔ | ✔ |
+|`TimeShift` | ✔ | ✔ | ✔ | ✔ |
+|`LaserF0` | ✔ | ✔ | ✔ | ✔ |
+|`LaserA0` | ✔ | ✔ | ✔ | ✔ |
+|`LaserFx` | ✔ | ✔ | ✔ | ✔ |
+|`LaserFy` | ✔ | ✔ | ✔ | ✔ |
+|`LaserAx` | ✔ | ✔ | ✔ | ✔ |
+|`LaserAy` | ✔ | ✔ | ✔ | ✔ |
+
+### Ellipticity
+
+The ellipticity ``\varepsilon`` defines the polarization type of the laser field.
+For special cases, `0` indicates linear polarization and `±1` indicates circular polarization. The electric field rotates clockwise for positive ellipticities and counter-clockwise for negative ones.
+
+
+
+
+### Azimuth of Principle Axis
+
+The azimuth angle ``\varphi`` of the principle axis defines a clockwise rotation of the laser field in the polarization plane.
+
+
+
+
+### Carrier-Envelope-Phase (CEP)
+
+The carrier-envelope-phase (CEP) ``\phi`` is the difference between the optical phase of the carrier wave and the envelope position. For few-cycle laser pulses, the influence of the CEP to the laser-matter interaction becomes significant.
+
+
+
+
+## Cos⁴-envelope Laser
+
+A [`Cos4Laser`](@ref) has a cos⁴-shaped-envelope:
+```math
+f_{\mathrm{env}}(t) =
+ \begin{cases}
+ \cos^4{\left[ \omega (t-t_0)/2N \right]}, & -NT/2 \leq t-t_0 \leq NT/2, \\
+ 0, & \mathrm{otherwise,} \\
+ \end{cases}
+```
+where ``\omega`` is the angular frequency of the laser field, ``N`` is the cycle number, ``T`` the period, and ``t_0`` denotes the peak time (corresponding to the `time_shift` in the constructor method).
+
+For the detailed usage, cf. the documentation of [`Cos4Laser`](@ref).
+
+```@docs
+Lasers.Cos4Laser
+```
+
+The following shows an example of [`Cos4Laser`](@ref) and its envelope shape.
+
+```@repl manual_lasers
+Cos4Laser(peak_int=1e14, wave_len=800.0, cyc_num=10, ellip=0)
+```
+
+
+
+
+## Cos²-envelope Laser
+
+A [`Cos2Laser`](@ref) has a cos²-shaped-envelope, similiar with the `Cos4Laser`:
+```math
+f_{\mathrm{env}}(t) =
+ \begin{cases}
+ \cos^2{\left[ \omega (t-t_0)/2N \right]}, & -NT/2 \leq t-t_0 \leq NT/2, \\
+ 0, & \mathrm{otherwise.} \\
+ \end{cases}
+```
+
+For the detailed usage, cf. the documentation of [`Cos2Laser`](@ref).
+
+```@docs
+Lasers.Cos2Laser
+```
+
+The following shows an example of [`Cos2Laser`](@ref) and its envelope shape.
+
+```@repl manual_lasers
+Cos2Laser(peak_int=1e14, wave_len=800.0, cyc_num=10, ellip=0)
+```
+
+
+
+
+## Gaussian-envelope Laser
+
+A [`GaussianLaser`](@ref) has a Gaussian-shaped-envelope:
+```math
+f_{\mathrm{env}}(t) = \exp{\left[ -(t-t_0)^2/\sigma^2 \right]} = \exp{\left[ -8\ln{2}(t-t_0)^2/\tau_{\mathrm{FWHM}}^2 \right]},
+```
+where ``\sigma`` is the temporal width of the laser (relating to the `spread_duration` in the constructor method) and ``\tau_{\mathrm{FWHM}}=2\sqrt{2\ln{2}}\sigma`` denotes the laser's temporal FWHM (Full-Width at Half Maximum).
+
+For the detailed usage, cf. the documentation of [`GaussianLaser`](@ref).
+
+```@docs
+Lasers.GaussianLaser
+```
+
+An example of [`GaussianLaser`](@ref) and its envelope shape are shown as follows.
+
+```@repl manual_lasers
+GaussianLaser(peak_int=1e14, wave_len=800.0, FWHM_duration=1103.2, ellip=0)
+```
+
+
+
+
+## Trapezoidal-envelope Laser
+
+A [`TrapezoidalLaser`](@ref) has a trapezoidal-shaped-envelope:
+```math
+f_{\mathrm{env}}(t) =
+ \begin{cases}
+ (t-t_0) / N_{\mathrm{on}}T,
+ & 0 < t-t_0 \leq N_{\mathrm{on}}T, \\
+ 1,
+ & N_{\mathrm{on}}T < t-t_0 \leq (N_{\mathrm{on}}+N_{\mathrm{const}})T, \\
+ 1-[t-t_0-(N_{\mathrm{on}}+N_{\mathrm{const}})T] / N_{\mathrm{off}}T,
+ & (N_{\mathrm{on}}+N_{\mathrm{const}})T < t-t_0 \leq (N_{\mathrm{on}}+N_{\mathrm{const}}+N_{\mathrm{off}})T, \\
+ 0,
+ & \mathrm{otherwise.}
+ \end{cases}
+```
+where ``N_{\mathrm{on}}, N_{\mathrm{const}}, N_{\mathrm{off}}`` are cycle numbers during the turn-on, constant, and turn-off stages. Note that for the `TrapezoidalLaser`, the ``t_0`` denotes the *time of rise* instead of time of peak, in contrast to the previous lasers.
+
+For the detailed usage, cf. the documentation of [`TrapezoidalLaser`](@ref).
+
+```@docs
+Lasers.TrapezoidalLaser
+```
+
+In the following is an example of [`TrapezoidalLaser`](@ref) and its envelope shape.
+
+```@repl manual_lasers
+l = TrapezoidalLaser(
+ peak_int=1e14, wave_len=800.0,
+ cyc_num_turn_on=3, cyc_num_turn_off=3, cyc_num_const=4,
+ ellip=0, t_shift=-551.6)
+```
+
+
+
diff --git a/docs/src/manual3_main_method.md b/docs/src/manual3_main_method.md
new file mode 100644
index 0000000..8742370
--- /dev/null
+++ b/docs/src/manual3_main_method.md
@@ -0,0 +1,335 @@
+# Main Method `performSFI`
+
+*This section introduces the main method* [`performSFI`](@ref)*, which performs the trajectory simulation and saves the final electron momentum spectrum.*
+
+```@contents
+Pages = ["manual3_main_method.md"]
+Depth = 3
+```
+
+```@meta
+CurrentModule = SemiclassicalSFI
+```
+
+
+## Brief Documentation
+
+```@docs
+performSFI
+```
+
+
+
+## Lasers & Targets
+
+### Lasers
+
+- `laser::Laser`
+
+A `Lasers.Laser` object containing information of the laser field.
+Cf. the documentation for [Lasers](@ref lasers_doc).
+
+### Targets
+
+- `target::Target`
+
+A `Targets.Target` object containing information of the atom/molecule target.
+Cf. the documentation for [Targets](@ref targets_doc).
+
+### Workflow for Preparation of the `Molecule` Target
+
+To use the [`Molecule`](@ref) as the target and the relating initial condition methods [MO-ADK](@ref MOADK) or [WFAT](@ref WFAT), some coefficients need to be calculated in advance.
+Here we present a workflow for preparation of the `Molecule` target before invoking the `performSFI` method, taking the carbon monoxide molecule as an example.
+
+!!! note "Note"
+ Currently only the PySCF is implemented as the library's molecular computation interface, which only supports the Linux platform.
+ What's more, in current version, only *close-shell* molecules are supported!
+
+
+#### Initialization
+
+First of all, initialize a `Molecule` object, and provide the necessary information of the molecule. Cf. the documentation of [`Molecule`](@ref).
+
+```julia
+using SemiclassicalSFI.Targets
+mol = Molecule(atoms=["C","O"], atom_coords=[0 0 -0.180; 0 0 0.950],
+ charge=0, name="Carbon Monoxide",
+ data_path="./Molecule_CarbonMonoxide.h5")
+```
+
+!!! warning "Unit of the atom coordinates"
+ The input of the atom coordinates is in **Angstrom (Å)**, which is not atomic unit (Bohr). 1 Bohr = 0.53 Å.
+
+!!! tip "Data saving of the Molecule object"
+ If the user specifies `data_path` in the constructor method of `Molecule`, the data would be automatically saved each time the user invokes the [`MolCalcAsympCoeff!`](@ref) and [`MolCalcWFATData!`](@ref).
+ However, if doesn't specify (in case the user does not wish to save the data), the data would not be saved, and the user has to manually invoke [`MolSaveDataAs`](@ref) to save the data afterwards.
+
+
+#### Calculate Asymptotic Coefficients
+
+To calculate the asymptotic coefficients of the `Molecule` object which is used in MO-SFA and [MO-ADK](@ref MOADK), invoke the [`MolCalcAsympCoeff!`](@ref) method:
+```julia
+MolCalcAsympCoeff!(mol)
+```
+
+For typical small molecules, using default parameters usually gives satisfactory results.
+However, for special demands, the user may refer to the documentation of [`MolCalcAsympCoeff!`](@ref) and [`Targets.MolecularCalculators.calc_asymp_coeff`](@ref) for more calculation parameters.
+
+#### Calculate WFAT Data
+
+To calculate the data necessary for the [WFAT](@ref WFAT) of the `Molecule` object, invoke the [`MolCalcWFATData!`](@ref) method:
+```julia
+MolCalcWFATData!(mol, orbitIdx_relHOMO = 0)
+```
+
+To obtain the data of other orbitals besides HOMO, the user may alter the `orbitIdx_relHOMO` parameter.
+For custom calculation parameters, refer to the documentation of [`MolCalcWFATData!`](@ref) and [`Targets.MolecularCalculators.calc_WFAT_data`](@ref).
+
+#### Setting the Molecule's Orientation
+
+The orientation of the molecule can be specified by invoking the [`SetMolRotation`](@ref) method:
+```julia
+SetMolRotation(mol, 0.0,π/2,π/3)
+```
+and can be obtained through the [`MolRotation`](@ref) method:
+```julia
+MolRotation(mol)
+```
+
+!!! tip "Testing molecular calculation setup"
+ The test sets of this library include molecular calculations, the user may run the tests of this library to check if the environment is correctly set up for molecular calculation:
+ ```julia
+ using Pkg
+ Pkg.test("SemiclassicalSFI")
+ # In pkg mode of REPL:
+ # (@v1.8) pkg> test SemiclassicalSFI
+ ```
+
+
+## Initial Condition Methods
+
+- `init_cond_method = <:ADK|:SFA|:SFAAE|:WFAT|:MOADK>`
+
+Method of electrons' initial conditions.
+Currently supports `:ADK`, `:SFA`, `:SFAAE` for atoms, and `:WFAT`, `:MOADK` for molecules.
+
+For more information about the theories, cf. [Theory - Initial Conditions](@ref theory_init_cond).
+
+!!! note "Note"
+ For initial condition methods [ADK](@ref ADK), [SFA](@ref SFA) or [SFA-AE](@ref SFAAE), you must specify an atom target of types [`HydrogenLikeAtom`](@ref) or [`SAEAtom`](@ref);
+ For initial condition methods [WFAT](@ref WFAT) or [MO-ADK](@ref MOADK), you must specify a molecule target of type [`Molecule`](@ref).
+
+ | | [ADK](@ref ADK) | [SFA](@ref SFA) | [SFA-AE](@ref SFAAE) | [WFAT](@ref WFAT) | [MO-ADK](@ref MOADK) |
+ | :------------------------- |:-:|:-:|:-:|:-:|:-:|
+ | [`HydrogenLikeAtom`](@ref) | ✔ | ✔ | ✔ | | |
+ | [`SAEAtom`](@ref) | ✔ | ✔ | ✔ | | |
+ | [`Molecule`](@ref) | | | | ✔ | ✔ |
+
+### Atomic Rate Prefix
+
+- `rate_prefix = <:ExpRate|:ExpPre|:ExpJac|:Full>`
+
+Prefix of the exponential term in the ionization rate (default `:ExpRate`).
+
+For atomic [ADK](@ref ADK), [SFA](@ref SFA) and [SFA-AE](@ref SFAAE), we obtained the ionization probability in the following form:
+```math
+\mathrm{d}W/\mathrm{d}\bm{k}_\perp \mathrm{d}t_{\mathrm{r}} = \bm{J}(k_d,t_{\mathrm{r}}) \lvert P_{\bm{p}}(t_{\mathrm{s}}) \rvert^2 \exp(2\ \mathrm{Im}\ \Phi_{\mathrm{tun}}),
+```
+where the ``\lvert P_{\bm{p}}(t_{\mathrm{s}}) \rvert^2`` denotes the prefactor, which has different expressions for different theories:
+
+| [SFA](@ref SFA) | [SFA-AE](@ref SFAAE) | [ADK](@ref ADK) |
+| :-: | :-: | :-: |
+| ``\{ [\bm{p}+\bm{A}(t_{\mathrm{s}})] \cdot \bm{F}(t_{\mathrm{s}}) \}^{-\alpha}`` | ``\left[ (k_\perp^2+2I_{\mathrm{p}})(F^2-\bm{k}_\perp \cdot \bm{F}') \right]^{-\alpha/2}`` | ``\left[ (k_\perp^2+2I_{\mathrm{p}})F^2\right]^{-\alpha/2}`` |
+
+The ``\bm{J}(k_d,t_{\mathrm{r}})`` denotes the Jacobian which arises from the coordinate transformation.
+
+For initial condition methods [ADK](@ref ADK), [SFA](@ref SFA) and [SFA-AE](@ref SFAAE),
+specifying `ExpRate` would not add any prefix besides the exponential term;
+`ExpPre` would include the prefactor ``\lvert P_{\bm{p}}(t_{\mathrm{s}}) \rvert^2``;
+`ExpJac` would include the Jacobian ``\bm{J}(k_d,t_{\mathrm{r}})``;
+`Full` indicates inclusion of both the prefactor and Jacobian.
+
+### Tunneling Exit Methods For Atomic ADK
+
+- `adk_tun_exit = <:IpF|:FDM|:Para>`
+
+Tunneling exit method for atomic ADK methods (when `init_cond_method==:ADK`) (default `:IpF`).
+Cf. [Tunneling Exit Methods For Atomic ADK](@ref tun_exit_atomic_adk).
+
+
+
+## Sampling Methods and Parameters
+
+In *SemiclassicalSFI.jl* the initial electrons are sampled in the ``(t_{\mathrm{r}},k_d,k_z)`` coordinate, where ``k_d`` denotes the initial momentum's component in the ``xy`` plane (which is perpendicular to the electric field).
+There are two ways of sampling in these coordinates, namely *step sampling* and *Monte-Carlo sampling*.
+
+- `sample_monte_carlo = false` : Determines whether Monte-Carlo sampling is used when generating electron samples (default `false`).
+
+- `sample_t_intv = (start,stop)` : Time interval in which the initial electrons are sampled.
+
+- `sample_t_num` : Number of time samples.
+
+- `sample_cutoff_limit = 1e-16` : The cut-off limit of the probability of the sampled electron, electrons with probabilities lower than the limit would be discarded (default `1e-16`).
+
+### Step Sampling
+
+In the step sampling scheme, the `sample_t_num` time samples are uniformly distributed in the interval `sample_t_intv`.
+In each time sample, a batch of electrons of different initial conditions are launched and collected, whose initial momenta ``\bm{k}_\perp`` are distributed on a Cartesian grid ``(k_d,k_z)``.
+The Cartesian grid is defined by the following parameters:
+
+- `ss_kd_max`, `ss_kd_num`, `ss_kz_max`, `ss_kz_num`
+
+In the ``k_d`` dimension, `ss_kd_num` samples distribute uniformly in the interval (-`ss_kd_max`,`ss_kd_max`);
+and in the ``k_z`` dimension, there are `ss_kz_num` equidistant samples in the interval (-`ss_kz_max`,`ss_kz_max`).
+
+The step sampling method is supported for all initial condition methods.
+
+!!! warning "ss_kz_num should be an even number"
+ Under the dipole case, electrons with zero initial ``k_z`` would remain zero ``k_z``, which may induce some problems in the step-length choice of the ODE solver and result in "final electron with anomalously large momentum" warning. Thus the sampling on the ``k_z`` direction should avoid zero. Setting `ss_kz_num` to an even number is a good choice.
+
+### Monte-Carlo Sampling
+
+In the Monte-Carlo sampling scheme, the `sample_t_num` time samples are randomly chosen in the `sample_t_intv`;
+A batch containing `mc_kp_num` electrons would be sampled in a single time sample, the electrons' initial momenta ``\bm{k}_\perp`` are also randomly sampled inside a circle ``k_d^2+k_z^2 \leq k_{\perp\mathrm{max}}^2``, where the ``k_{\perp\mathrm{max}}`` is defined in the parameter as `mc_kp_max`.
+
+- `mc_kp_num` : Number of kp (initial momentum which is perpendicular to field direction, two dimensional) samples in a single time sample.
+- `mc_kp_max` : Maximum value of momentum's transversal component (perpendicular to field direction).
+
+Currently the Monte-Carlo sampling method is only supported for the ADK initial condition method.
+
+
+
+## Trajectory Simulation
+
+After preparation of the initial electrons, the electrons evolve classically in the combined potential of the nucleus and laser field, and the trajectory simulation terminates at `traj_t_final`.
+
+- `traj_t_final` : Time when every trajectory simulation ends.
+
+### Phase Methods
+
+- `traj_phase_method = <:CTMC|:QTMC|:SCTS>`
+
+Method of classical trajectories' phase (default `CTMC`).
+Currently `:QTMC` and `:SCTS` only support atomic cases.
+
+| | [ADK](@ref ADK) | [SFA](@ref SFA) | [SFA-AE](@ref SFAAE) | [WFAT](@ref WFAT) | [MO-ADK](@ref MOADK) |
+| :---------------- |:-:|:-:|:-:|:-:|:-:|
+| [CTMC](@ref CTMC) | ✔ | ✔ | ✔ | ✔ | ✔ |
+| [QTMC](@ref QTMC) | ✔ | ✔ | ✔ | | |
+| [SCTS](@ref SCTS) | ✔ | ✔ | ✔ | | |
+
+For more information about phase methods, cf. [Theory - Trajectory Simulation and Phase Methods](@ref theory_traj_phase).
+
+### Non-dipole Effect
+
+- `traj_nondipole = false`
+
+Determines whether the non-dipole effect is taken account in the simulation (default `false`).
+
+For more information about the non-dipole effects, cf. [Theory - Non-dipole Effects on the Trajectory Motion](@ref traj_nondipole).
+
+Currently, all targets support the inclusion of non-dipole effects.
+
+### Accuracy Control
+
+- `traj_rtol = 1e-6`
+
+Relative error tolerance when solving classical trajectories using adaptive methods (default `1e-6`).
+
+The classical trajectories of the electrons are obtained by solving ordinary differential equations using the [OrdinaryDiffEq.jl](https://github.com/SciML/OrdinaryDiffEq.jl) package, where the solver would adjust the time-step length according to the relative error.
+Stricter error tolerance is favorable if the electron moves close to the nucleus, which usually happens in a linearly polarized laser field.
+
+### GPU Acceleration (Experimental)
+
+- `traj_GPU = false`
+
+Determines whether to enable GPU acceleration in trajectory simulation (default `false`).
+
+The ordinary differential equations related to the electrons' classical trajectories can also be solved by GPU via the [DiffEqGPU.jl](https://github.com/SciML/DiffEqGPU.jl) package.
+The GPU acceleration now only supports the NVIDIA graphic cards, which also requires the CUDA driver to be correctly installed.
+
+!!! tip "Testing GPU capability"
+ The test sets of this library include GPU tests, the user may run the tests of this library to check if the graphic card is ready for GPU acceleration:
+ ```julia
+ using Pkg
+ Pkg.test("SemiclassicalSFI")
+ # In pkg mode of REPL:
+ # (@v1.8) pkg> test SemiclassicalSFI
+ ```
+
+!!! compat "Note: Experimental feature"
+ GPU acceleration is an experimental feature of the package and the API may change in the near future.
+
+!!! note "Note: GPU acceleration efficiency"
+ The GPU's advantage of speed over CPU may not be so obvious in the trajectory simulation.
+ Sometimes the program is even slower when using GPU.
+ Thus it's better to benchmark for specific hardware and choose the most efficient option.
+
+## Final Electron Collecting & Saving
+
+After the trajectory simulation ends, the electrons would be analyzed and collected.
+Those with positive energies finally become free electrons and would reach the detectors;
+while those with negative energies finally fall on Rydberg states.
+The collected momentum spectra and Rydberg spectra would be saved in an HDF5 file together with the simulation abstract.
+This library provides some parameters to customize the collecting and saving procedure.
+
+### 2D/3D Momentum Spectrum Collecting
+
+- `final_p_max = (pxMax,pyMax,pzMax)` : Boundaries of final momentum spectrum collected in three dimensions.
+- `final_p_num = (pxNum,pyNum,pzNum)` : Numbers of final momentum spectrum collected in three dimensions.
+- `save_3D_spec = false` : Determines whether the 3D momentum spectrum is saved (if not, will only save 2D by flattening on the xy plane) (default `false`).
+
+Electrons with positive final energies would be collected and placed on the 3D momentum grid determined by `final_p_max` and `final_p_num`.
+When setting `save_3D_spec = false`, the three-dimensional final momentum spectrum would be squashed into two-dimensional ones (by summing over the z axis).
+
+The ``p_x``, ``p_y`` and ``p_z`` grids would be saved in entries that are named after `px`, `py` and `pz` respectively in the output file, and the final momentum spectrum would be saved in the `momentum_spec_2D` or `momentum_spec_3D` entry (compressed to reduce the file size).
+
+### Rydberg Final State Collecting
+
+- `final_ryd_collect = false` : Determines whether the rydberg final states are collected (default `false`).
+- `final_ryd_n_max` : Determines the maximum principle quantum number n for rydberg final states to be collected.
+
+Rydberg final states would be collected if `final_ryd_collect` is set to `true`.
+Only Rydberg states with principle quantum number that below `final_ryd_n_max` would be collected.
+
+The Rydberg spectrum is saved as a three-dimensional array in the entry named after `ryd_spec` in the output file.
+To get the probability of Rydberg state at ``(n,l,m)``, index with the indices `(n,l+1,m+final_ryd_n_max)`.
+
+### Output File
+
+- `save_path` : Output HDF5 file path.
+
+The output HDF5 file would be saved in the `save_path`.
+If the parameter is left unspecified or an error occurs when trying to write to the specified path, the output path would be set to `./SCSFI-yyyymmdd-hhmmss.h5`.
+
+Apart from the information related to the momentum and Rydberg spectrum, an abstract encoded in YAML, which contains necessary input parameters, is also saved in the `abstract` entry of the output file.
+
+The output file can be opened and accessed using the [HDF5.jl](https://github.com/JuliaIO/HDF5.jl) or [JLD2.jl](https://github.com/JuliaIO/JLD2.jl) packages.
+Below shows the structure of typical output files.
+```
+# saving only 2D momentum spectrum
+🗂️ HDF5.File
+├─ 🔢 abstract
+├─ 🔢 info
+├─ 🔢 ion_prob
+├─ 🔢 ion_prob_uncollected
+├─ 🔢 momentum_spec_2D
+├─ 🔢 num_effective_traj
+├─ 🔢 px
+└─ 🔢 py
+```
+```
+# saving 3D momentum spectrum
+🗂️ HDF5.File
+├─ 🔢 abstract
+├─ 🔢 info
+├─ 🔢 ion_prob
+├─ 🔢 ion_prob_uncollected
+├─ 🔢 momentum_spec_2D
+├─ 🔢 momentum_spec_3D
+├─ 🔢 num_effective_traj
+├─ 🔢 px
+├─ 🔢 py
+└─ 🔢 pz
+```
+
diff --git a/docs/src/plot_manual2.jl b/docs/src/plot_manual2.jl
new file mode 100644
index 0000000..bddfa73
--- /dev/null
+++ b/docs/src/plot_manual2.jl
@@ -0,0 +1,203 @@
+# code to generate plots used in manual2_lasers.md
+
+using SemiclassicalSFI
+using SemiclassicalSFI.Lasers
+using Plots
+pyplot()
+
+# =======================
+
+l = Cos4Laser(peak_int=1e14, wave_len=800.0, cyc_num=6, ellip=0)
+F0 = LaserF0(l)
+Fx = LaserFx(l)
+Fy = LaserFy(l)
+t = -400:1:400
+Fxt = Fx.(t)/F0
+Fyt = Fy.(t)/F0
+color = RGB(0/255,154/255,250/255)
+color_trans = RGBA(0/255,154/255,250/255,0.3)
+f1=plot(t,Fxt,Fyt, label="ε=0.0", color=color)
+plot!(repeat([-400],length(t)),Fxt,Fyt, label=:none, color=color_trans)
+plot!(t,Fxt,repeat([-1.0],length(t)), label=:none, color=color_trans)
+plot!(t,repeat([1.0],length(t)),Fyt, label=:none, color=color_trans)
+plot!(xlim=(-400,400),ylim=(-1,1),zlim=(-1,1))
+plot!(xlabel=raw"$t$/a.u.", ylabel=raw"$F_x/F_0$", zlabel=raw"$F_y/F_0$")
+
+l = Cos4Laser(peak_int=1e14, wave_len=800.0, cyc_num=6, ellip=0.5)
+F0 = LaserF0(l)
+Fx = LaserFx(l)
+Fy = LaserFy(l)
+t = -400:1:400
+Fxt = Fx.(t)/F0
+Fyt = Fy.(t)/F0
+color = RGB(255/255,127/255,24/255)
+color_trans = RGBA(255/255,127/255,24/255,0.3)
+f2=plot(t,Fxt,Fyt, label="ε=0.5", color=color)
+plot!(repeat([-400],length(t)),Fxt,Fyt, label=:none, color=color_trans)
+plot!(t,Fxt,repeat([-1.0],length(t)), label=:none, color=color_trans)
+plot!(t,repeat([1.0],length(t)),Fyt, label=:none, color=color_trans)
+plot!(xlim=(-400,400),ylim=(-1,1),zlim=(-1,1))
+plot!(xlabel=raw"$t$/a.u.", ylabel=raw"$F_x/F_0$", zlabel=raw"$F_y/F_0$")
+
+l = Cos4Laser(peak_int=1e14, wave_len=800.0, cyc_num=6, ellip=1.0)
+F0 = LaserF0(l)
+Fx = LaserFx(l)
+Fy = LaserFy(l)
+t = -400:1:400
+Fxt = Fx.(t)/F0
+Fyt = Fy.(t)/F0
+color = RGB(44/255,160/255,44/255)
+color_trans = RGBA(44/255,160/255,44/255,0.3)
+f3=plot(t,Fxt,Fyt, label="ε=1.0", color=color)
+plot!(repeat([-400],length(t)),Fxt,Fyt, label=:none, color=color_trans)
+plot!(t,Fxt,repeat([-1.0],length(t)), label=:none, color=color_trans)
+plot!(t,repeat([1.0],length(t)),Fyt, label=:none, color=color_trans)
+plot!(xlim=(-400,400),ylim=(-1,1),zlim=(-1,1))
+plot!(xlabel=raw"$t$/a.u.", ylabel=raw"$F_x/F_0$", zlabel=raw"$F_y/F_0$")
+
+l = Cos4Laser(peak_int=1e14, wave_len=800.0, cyc_num=6, ellip=-1.0)
+F0 = LaserF0(l)
+Fx = LaserFx(l)
+Fy = LaserFy(l)
+t = -400:1:400
+Fxt = Fx.(t)/F0
+Fyt = Fy.(t)/F0
+color = RGB(214/255,39/255,40/255)
+color_trans = RGBA(214/255,39/255,40/255,0.3)
+f4=plot(t,Fxt,Fyt, label="ε=-1.0", color=color)
+plot!(repeat([-400],length(t)),Fxt,Fyt, label=:none, color=color_trans)
+plot!(t,Fxt,repeat([-1.0],length(t)), label=:none, color=color_trans)
+plot!(t,repeat([1.0],length(t)),Fyt, label=:none, color=color_trans)
+plot!(xlim=(-400,400),ylim=(-1,1),zlim=(-1,1))
+plot!(xlabel=raw"$t$/a.u.", ylabel=raw"$F_x/F_0$", zlabel=raw"$F_y/F_0$")
+
+f=plot(f1,f2,f3,f4, layout=(2,2), size=(1000,800), fmt=:svg)
+savefig("./docs/src/assets/manual2_lasers_ellip.svg")
+
+# =======================
+
+l = Cos4Laser(peak_int=1e14, wave_len=800.0, cyc_num=6, ellip=0, azi=0π)
+F0 = LaserF0(l)
+Fx = LaserFx(l)
+Fy = LaserFy(l)
+t = -400:1:400
+Fxt = Fx.(t)/F0
+Fyt = Fy.(t)/F0
+color = RGB(0/255,154/255,250/255)
+color_trans = RGBA(0/255,154/255,250/255,0.3)
+f1=plot(t,Fxt,Fyt, label="azi=0", color=color)
+plot!(repeat([-400],length(t)),Fxt,Fyt, label=:none, color=color_trans)
+plot!(t,Fxt,repeat([-1.0],length(t)), label=:none, color=color_trans)
+plot!(t,repeat([1.0],length(t)),Fyt, label=:none, color=color_trans)
+plot!(xlim=(-400,400),ylim=(-1,1),zlim=(-1,1))
+plot!(xlabel=raw"$t$/a.u.", ylabel=raw"$F_x/F_0$", zlabel=raw"$F_y/F_0$")
+
+l = Cos4Laser(peak_int=1e14, wave_len=800.0, cyc_num=6, ellip=0, azi=π/6)
+F0 = LaserF0(l)
+Fx = LaserFx(l)
+Fy = LaserFy(l)
+t = -400:1:400
+Fxt = Fx.(t)/F0
+Fyt = Fy.(t)/F0
+color = RGB(255/255,127/255,24/255)
+color_trans = RGBA(255/255,127/255,24/255,0.3)
+f2=plot(t,Fxt,Fyt, label="azi=π/6 (30°)", color=color)
+plot!(repeat([-400],length(t)),Fxt,Fyt, label=:none, color=color_trans)
+plot!(t,Fxt,repeat([-1.0],length(t)), label=:none, color=color_trans)
+plot!(t,repeat([1.0],length(t)),Fyt, label=:none, color=color_trans)
+plot!(xlim=(-400,400),ylim=(-1,1),zlim=(-1,1))
+plot!(xlabel=raw"$t$/a.u.", ylabel=raw"$F_x/F_0$", zlabel=raw"$F_y/F_0$")
+
+l = Cos4Laser(peak_int=1e14, wave_len=800.0, cyc_num=6, ellip=0, azi=π/3)
+F0 = LaserF0(l)
+Fx = LaserFx(l)
+Fy = LaserFy(l)
+t = -400:1:400
+Fxt = Fx.(t)/F0
+Fyt = Fy.(t)/F0
+color = RGB(44/255,160/255,44/255)
+color_trans = RGBA(44/255,160/255,44/255,0.3)
+f3=plot(t,Fxt,Fyt, label="azi=π/3 (60°)", color=color)
+plot!(repeat([-400],length(t)),Fxt,Fyt, label=:none, color=color_trans)
+plot!(t,Fxt,repeat([-1.0],length(t)), label=:none, color=color_trans)
+plot!(t,repeat([1.0],length(t)),Fyt, label=:none, color=color_trans)
+plot!(xlim=(-400,400),ylim=(-1,1),zlim=(-1,1))
+plot!(xlabel=raw"$t$/a.u.", ylabel=raw"$F_x/F_0$", zlabel=raw"$F_y/F_0$")
+
+l = Cos4Laser(peak_int=1e14, wave_len=800.0, cyc_num=6, ellip=0, azi=π/2)
+F0 = LaserF0(l)
+Fx = LaserFx(l)
+Fy = LaserFy(l)
+t = -400:1:400
+Fxt = Fx.(t)/F0
+Fyt = Fy.(t)/F0
+color = RGB(214/255,39/255,40/255)
+color_trans = RGBA(214/255,39/255,40/255,0.3)
+f4=plot(t,Fxt,Fyt, label="azi=π/2 (90°)", color=color)
+plot!(repeat([-400],length(t)),Fxt,Fyt, label=:none, color=color_trans)
+plot!(t,Fxt,repeat([-1.0],length(t)), label=:none, color=color_trans)
+plot!(t,repeat([1.0],length(t)),Fyt, label=:none, color=color_trans)
+plot!(xlim=(-400,400),ylim=(-1,1),zlim=(-1,1))
+plot!(xlabel=raw"$t$/a.u.", ylabel=raw"$F_x/F_0$", zlabel=raw"$F_y/F_0$")
+
+f=plot(f1,f2,f3,f4, layout=(2,2), size=(1000,800), fmt=:svg)
+savefig("./docs/src/assets/manual2_lasers_azimuth.svg")
+
+# =======================
+
+t = -400:1:400
+l = Cos4Laser(peak_int=1e14, wave_len=800.0, cyc_num=6, ellip=0, cep=0)
+plot(t, UnitEnvelope(l).(t), color=:gray, fill=0, α=0.1, label=:none)
+plot!(t, -1 .* UnitEnvelope(l).(t), color=:gray, fill=0, α=0.1, label=:none)
+
+l = Cos4Laser(peak_int=1e14, wave_len=800.0, cyc_num=6, ellip=0, cep=0)
+plot!(t, LaserAx(l).(t) ./ LaserA0(l), label="CEP=0", color=RGB(64/255,183/255,173/255))
+l = Cos4Laser(peak_int=1e14, wave_len=800.0, cyc_num=6, ellip=0, cep=π/3)
+plot!(t, LaserAx(l).(t) ./ LaserA0(l), label="CEP=π/3", color=RGB(52/255,143/255,167/255))
+l = Cos4Laser(peak_int=1e14, wave_len=800.0, cyc_num=6, ellip=0, cep=2π/3)
+plot!(t, LaserAx(l).(t) ./ LaserA0(l), label="CEP=2π/3", color=RGB(55/255,101/255,158/255))
+l = Cos4Laser(peak_int=1e14, wave_len=800.0, cyc_num=6, ellip=0, cep=π)
+plot!(t, LaserAx(l).(t) ./ LaserA0(l), label="CEP=π", color=RGB(65/255,61/255,163/255))
+
+f=plot!(xlabel=raw"$t$/a.u.", ylabel=raw"$f_\mathrm{env}$", xlim=(-300,300), ylim=(-1.05,1.05), size=(500,400), framestyle=:box, fmt=:svg)
+savefig("./docs/src/assets/manual2_lasers_cep.svg")
+
+# =======================
+
+t = -600:1:600
+l = Cos4Laser(peak_int=1e14, wave_len=800.0, cyc_num=10, ellip=0)
+plot(t, UnitEnvelope(l).(t), color=:gray, fill=0, α=0.2)
+plot!(t, -1 .* UnitEnvelope(l).(t), color=:gray, fill=0, α=0.2)
+plot!(t, LaserAx(l).(t) ./ LaserA0(l), color=RGB(0/255,154/255,250/255))
+plot!(xlabel=raw"$t$/a.u.", ylabel=raw"$f_\mathrm{env}$", xlim=(-600,600), ylim=(-1.05,1.05), legend=:none, size=(500,400), framestyle=:box, fmt=:svg)
+savefig("./docs/src/assets/manual2_lasers_env_cos4.svg")
+
+# =======================
+
+t = -600:1:600
+l = Cos2Laser(peak_int=1e14, wave_len=800.0, cyc_num=10, ellip=0)
+plot(t, UnitEnvelope(l).(t), color=:gray, fill=0, α=0.2)
+plot!(t, -1 .* UnitEnvelope(l).(t), color=:gray, fill=0, α=0.2)
+plot!(t, LaserAx(l).(t) ./ LaserA0(l), color=RGB(0/255,154/255,250/255))
+plot!(xlabel=raw"$t$/a.u.", ylabel=raw"$f_\mathrm{env}$", xlim=(-600,600), ylim=(-1.05,1.05), legend=:none, size=(500,400), framestyle=:box, fmt=:svg)
+savefig("./docs/src/assets/manual2_lasers_env_cos2.svg")
+
+# =======================
+
+t = -1000:1:1000
+l = GaussianLaser(peak_int=1e14, wave_len=800.0, FWHM_duration=1103.2, ellip=0)
+plot(t, UnitEnvelope(l).(t), color=:gray, fill=0, α=0.2)
+plot!(t, -1 .* UnitEnvelope(l).(t), color=:gray, fill=0, α=0.2)
+plot!(t, LaserAx(l).(t) ./ LaserA0(l), color=RGB(0/255,154/255,250/255))
+plot!(xlabel=raw"$t$/a.u.", ylabel=raw"$f_\mathrm{env}$", xlim=(-1000,1000), ylim=(-1.05,1.05), legend=:none, size=(500,400), framestyle=:box, fmt=:svg)
+savefig("./docs/src/assets/manual2_lasers_env_gaussian.svg")
+
+# =======================
+
+t = -600:1:600
+l = TrapezoidalLaser(peak_int=1e14, wave_len=800.0, cyc_num_turn_on=3, cyc_num_turn_off=3, cyc_num_const=4, ellip=0, t_shift=-551.6)
+plot(t, UnitEnvelope(l).(t), color=:gray, fill=0, α=0.2)
+plot!(t, -1 .* UnitEnvelope(l).(t), color=:gray, fill=0, α=0.2)
+plot!(t, LaserAx(l).(t) ./ LaserA0(l), color=RGB(0/255,154/255,250/255))
+plot!(xlabel=raw"$t$/a.u.", ylabel=raw"$f_\mathrm{env}$", xlim=(-600,600), ylim=(-1.05,1.05), legend=:none, size=(500,400), framestyle=:box, fmt=:svg)
+savefig("./docs/src/assets/manual2_lasers_env_trapezoidal.svg")
diff --git a/docs/src/theory1_initial_conditions.md b/docs/src/theory1_initial_conditions.md
new file mode 100644
index 0000000..3b1724e
--- /dev/null
+++ b/docs/src/theory1_initial_conditions.md
@@ -0,0 +1,344 @@
+# [Theory: Initial Conditions](@id theory_init_cond)
+
+*This section reviews commonly-used theories used to provide initial conditions in the trajectory simulations.*
+
+A number of theories can be adapted to provide initial conditions of the classical electrons in the trajectory simulation scheme.
+The initial condition usually consists of three properties:
+
+- Initial position ``\bm{r}_0`` (i.e., the tunneling exit position);
+- Initial momentum ``\bm{p}_0``, we note that in the trajectory simulation schemes, initial momenta are usually denoted using ``\bm{k}_0``;
+- The corresponding ionization probability ``W`` carried by each electron sample, depending on the time-dependent laser field and the properties of the target atoms/molecules.
+
+In the following we will give a brief review on the available theories we implemented in *SemiclassicalSFI.jl*.
+Atomic units (a.u.) are used throughout unless stated otherwise.
+
+```@contents
+Pages = ["theory1_initial_conditions.md"]
+```
+
+## [Strong-Field Approximation (SFA)](@id SFA)
+
+The Strong-Field Approximation (SFA) [^Popruzhenko_2014] is originated from the Keldysh theory of strong-field ionization.
+Compared with the pertubative methods and adiabatic tunneling theories, the SFA is able to predict both the multi-photon and the tunneling process during the laser-atom interaction, as well as high-order non-pertubative phenomenona such as the above-threshold ionization (ATI) because it fully includes the non-adiabatic effect of the laser-atom interaction.
+The broad scope of SFA has contributed to its widespread application in theoretical investigations of strong-field ionization.
+
+Considering the electron evolving under a combined field of the Coulomb field ``V(\bm{r})`` and the laser field ``\bm{F}(t)=-\partial_t \bm{A}(t)``, under the length gauge (LG), its Hamiltonian reads
+```math
+H = \frac12 \bm{p}^2 + V(\bm{r}) + \bm{F}(t)\cdot\bm{r}.
+```
+Denoting ``\ket{\Psi_0} = \ket{\psi_0} \mathrm{e}^{\mathrm{i}I_{\mathrm{p}}t}`` as the unperturbed initial state with ionization potential of ``I_{\mathrm{p}}``, ``\ket{\Psi_{\bm{p}}}`` as the continuum state of momentum ``\bm{p}``, and
+```math
+U(t_{\mathrm{f}},t_0) = \exp \left[ -\mathrm{i} \int_{t_0}^{t_{\mathrm{f}}} H(\tau) \mathrm{d}\tau \right]
+```
+the time-evolution operator, the transition amplitude between the initial state (at ``t_0``) and the final state of momentum ``\bm{p}`` (at ``t_{\mathrm{f}}``) is written as
+```math
+M_{\bm{p}} = \braket{ \Psi_{\bm{p}} | U(t_{\mathrm{f}},t_0) | \Psi_0 }.
+```
+
+Here lies the key idea of SFA: when the influence of the Coulomb field to the ionized electrons is weak compared with that of the external laser field, we may neglect the influence of the Coulomb field in the expression of ``M_{\bm{p}}`` by replacing the time-evolution operator with a Coulomb-free one ``U_{\mathrm{f}}``, and meanwhile replacing the continuum state with the Volkov state ``\ket{\Psi^{\mathrm{V}}_{\bm{p}}}`` which represents a free electron evolving under the same laser field:
+```math
+M_{\bm{p}} = \braket{ \Psi^{\mathrm{V}}_{\bm{p}} | U_{\mathrm{f}}(t_{\mathrm{f}},t_{\mathrm{0}}) | \Psi_0 },
+```
+where the Volkov state under the LG is the product of a plane wave and a phase factor:
+```math
+\ket{ \Psi^{\mathrm{V}}_{\bm{p}} } = \ket{ \bm{p}+\bm{A}(t) } \mathrm{e}^{\mathrm{i}S_{\bm{p}}(t)},
+```
+and the phase has the expression:
+```math
+S_{\bm{p}}(t) = - \int^{t} \frac12 [\bm{p}+\bm{A}(\tau)]^2 \mathrm{d}\tau.
+```
+In this way the ``M_{\bm{p}}`` is expressed as
+```math
+M_{\bm{p}} = -\mathrm{i} \int_{t_{\mathrm{0}}}^{t_{\mathrm{f}}} \braket{ \bm{p}+\bm{A}(\tau) | \bm{F}(\tau)\cdot\bm{r} | \psi_0 } \mathrm{e}^{-\mathrm{i}\tilde{S}_{\bm{p}}(\tau)} \mathrm{d}\tau,
+```
+and we note that here we have extracted the phase factor of ``\ket{\Psi_0}`` and combined it with the former ``\mathrm{e}^{-\mathrm{i}S_{\bm{p}}(t)}``, giving
+```math
+\tilde{S}_{\bm{p}}(t) = - \int^{t} \left[ \frac12 [\bm{p}+\bm{A}(\tau)]^2 + I_{\mathrm{p}} \right] \mathrm{d}\tau.
+```
+
+An additional saddle-point approximation facilitates preparation of initial conditions of the electron trajectories.
+The variation of ``\tilde{S}_{\bm{p}}(t)`` is much more sensitive than that of ``\braket{ \bm{p}+\bm{A}(t) | \bm{F}(t)\cdot\bm{r} | \psi_0 }`` as ``t`` varies, which leads to a fact that the whole integrand in our latest expression of ``M_{\bm{p}}`` oscillates in its complex phase and its values cancel out each other in most cases, except when the variation of the phase ``\tilde{S}_{\bm{p}}(t)`` becomes stable, i.e., at the saddle points of ``\tilde{S}_{\bm{p}}(t)``. The saddle points ``t_{\mathrm{s}}=t_{\mathrm{r}}+\mathrm{i}t_{\mathrm{i}}`` are the zeroes of the derivative of the complex function ``\tilde{S}_{\bm{p}}(t)``, which satisfy
+```math
+-\partial_t \tilde{S}_{\bm{p}}(t) |_{t=t_{\mathrm{s}}} = \frac12 [\bm{p}+\bm{A}(t_{\mathrm{s}})]^2 + I_{\mathrm{p}} = 0.
+```
+The integral can be approximated by a summation over the saddle points:
+```math
+M_{\bm{p}} \approx \sum_{t_{\mathrm{s}}} P_{\bm{p}}(t_{\mathrm{s}}) \mathrm{e}^{-\mathrm{i}\tilde{S}_{\bm{p}}(t_{\mathrm{s}})},
+```
+where ``P_{\bm{p}}(t_{\mathrm{s}})`` denotes the prefactor.
+Here we use a modified version of SFA which takes account of the Coulomb tail [^Kjeldsen_2006] [^Milosevic_2006], which gives the prefactor
+```math
+P_{\bm{p}}(t_{\mathrm{s}}) = \{ [\bm{p}+\bm{A}(t_{\mathrm{s}})] \cdot \bm{F}(t_{\mathrm{s}}) \}^{-\alpha/2},
+```
+where ``\alpha = 1+Z/\sqrt{2I_{\mathrm{p}}}``, and ``Z`` is the asymptotic charge of the nucleus.
+The phase ``\tilde{S}_{\bm{p}}(t_{\mathrm{s}})`` is obtained by solving the integral
+```math
+\begin{aligned}
+ \tilde{S}_{\bm{p}}(t_{\mathrm{s}})
+ &= \int_{t_{\mathrm{s}}}^{\infty} \left[ \frac12 [\bm{p}+\bm{A}(\tau)]^2 + I_{\mathrm{p}} \right] \mathrm{d}\tau \\
+ &= \left( \int_{t_{\mathrm{s}}}^{t_{\mathrm{r}}} + \int_{t_{\mathrm{r}}}^{\infty} \right) \left[ \frac12 [\bm{p}+\bm{A}(\tau)]^2 + I_{\mathrm{p}} \right] \mathrm{d}\tau \\
+ &= \Phi_{\mathrm{tun}} + \Phi_{\mathrm{traj}},
+\end{aligned}
+```
+where ``\Phi_{\mathrm{tun}}`` represents the complex phase accumulation during the tunneling process, whose real part denotes the quantum phase, while its imaginary part, is related to the ionization probability; the ``\Phi_{\mathrm{traj}}``, is the phase accumulation during the electron trajectory motion in the continuum.
+
+The SFA provides the final momentum distribution, while the trajectory simulation requires initial conditions of the eletrons.
+To utilize the SFA to give initial conditions, we suppose that the classical electron is ejected at time ``t_{\mathrm{r}}`` at tunneling exit ``\bm{r}_0`` with momentum ``\bm{k}_0``.
+The initial momentum ``\bm{k}_0``, neglecting the Coulomb potential, is related to the final momentum ``\bm{p}`` through
+```math
+\bm{p} = \bm{k}_0 - \int_{t_{\mathrm{r}}}^{\infty} \bm{F}(\tau) \mathrm{d}\tau = \bm{k}_0 - \bm{A}(t_{\mathrm{r}}).
+```
+The initial position ``\bm{r}_0``, i.e., the tunneling exit, is found by constructing a quantum tunneling trajectory.
+The beginning of the trajectory, i.e., the tunneling entrance, has a real part of zero; the electron tunnels through the barrier during the time interval ``t_{\mathrm{s}}`` to ``t_{\mathrm{r}}`` and emerges as a classical electron at the tunneling exit with real position and momentum.
+In this way we obtain the expression of the initial position:
+```math
+\bm{r}_0 = \mathrm{Re} \int_{t_{\mathrm{s}}}^{t_{\mathrm{r}}} \bm{A}(\tau) \mathrm{d}\tau = \mathrm{Im} \int_{0}^{t_{\mathrm{i}}} \bm{A}(t_{\mathrm{r}}+\mathrm{i}\tau) \mathrm{d}\tau.
+```
+The probablity density (in the final momentum space) carried by the electron sample is
+```math
+\mathrm{d}W/\mathrm{d}\bm{p} = \lvert P_{\bm{p}}(t_{\mathrm{s}}) \rvert^2 \exp(2\ \mathrm{Im}\ \Phi_{\mathrm{tun}}).
+```
+
+We note that the ionization probability is expressed in the coordinate of final momentum ``(p_x,p_y,p_z)``.
+However, in the trajectory simulation, the initial electrons are sampled in the coordinate ``(t_{\mathrm{r}},k_d,k_z)``, where ``k_d`` denotes the initial momentum's component in the ``xy`` plane (which is perpendicular to the electric field).
+Thus, adding a Jacobian in the prefix of the ionization probability is required if we sample the initial electrons within such coordinate, the transformed expression reads
+```math
+\mathrm{d}W/\mathrm{d}\bm{k}_\perp \mathrm{d}t_{\mathrm{r}} = \bm{J}(k_d,t_{\mathrm{r}}) \lvert P_{\bm{p}}(t_{\mathrm{s}}) \rvert^2 \exp(2\ \mathrm{Im}\ \Phi_{\mathrm{tun}}),
+```
+where the Jacobian is
+```math
+\bm{J}(k_d,t_{\mathrm{r}}) = \frac{\partial(p_x,p_y)}{\partial(k_d,t_{\mathrm{r}})} =
+\begin{vmatrix}
+ \partial p_x/\partial k_d & \partial p_x/\partial t_{\mathrm{r}} \\
+ \partial p_y/\partial k_d & \partial p_y/\partial t_{\mathrm{r}}
+\end{vmatrix}
+.
+```
+
+[^Popruzhenko_2014]: S. V. Popruzhenko, Keldysh Theory of Strong Field Ionization: History, Applications, Difficulties and Perspectives. *J. Phys. B: At. Mol. Opt. Phys.* **47**, 204001 (2014). DOI:[10.1088/0953-4075/47/20/204001](https://dx.doi.org/10.1088/0953-4075/47/20/204001)
+[^Kjeldsen_2006]: T. K. Kjeldsen *et al.*, Strong-Field Ionization of Atoms and Molecules: The Two-Term Saddle-Point Method. *Phys. Rev. A* **74**, 023407 (2006). DOI:[10.1103/PhysRevA.74.023407](https://dx.doi.org/10.1103/PhysRevA.74.023407)
+[^Milosevic_2006]: D. B. Milošević *et al.*, Above-Threshold Ionization by Few-Cycle Pulses. *J. Phys. B: At. Mol. Opt. Phys.* **39**, R203–R262 (2006). DOI: [10.1088/0953-4075/39/14/R01](https://dx.doi.org/10.1088/0953-4075/39/14/R01)
+
+
+## [SFA with Adiabatic Expansion (SFA-AE)](@id SFAAE)
+
+For small Keldysh parameter ``\gamma``, the non-adiabatic effect is not significant, thus an adiabatic expansion scheme can be carried out to develop a modified theory based on the SFA, which is named after the SFA with adiabatic expansion (SFA-AE) [^Ni_2018]. It includes the non-adiabatic effect to a large extent and is capable of giving similar results compared with that given by the SFA under small Keldysh parameters.
+
+The SFA-AE is applicable when the Keldysh parameter is small or the non-adiabatic effect is insignificant, and we recall that in the SFA there is a corresponding quantity ``t_{\mathrm{i}}`` which quantifies the non-adiabacity of tunneling.
+For small ``t_{\mathrm{i}}``, we expand the vector potential ``\bm{A}(t_{\mathrm{r}} + \mathrm{i}t_{\mathrm{i}})`` at ``t_{\mathrm{r}}``, up to the second order of ``t_{\mathrm{i}}``:
+```math
+\bm{A}(t_{\mathrm{r}} + \mathrm{i}t_{\mathrm{i}}) = \bm{A}(t_{\mathrm{r}}) - \mathrm{i}t_{\mathrm{i}}\bm{F}(t_{\mathrm{r}}) + \frac12 t_{\mathrm{i}}^2 \bm{F}'(t_{\mathrm{r}}) + o(t_{\mathrm{i}}^2).
+```
+Inserting the above expression into the saddle-point equation in the SFA gives
+```math
+\bm{k}(t_{\mathrm{r}}) \cdot \bm{F}(t_{\mathrm{r}}) = 0,
+```
+and
+```math
+t_{\mathrm{i}} = \sqrt{\frac{k^2(t_{\mathrm{r}})+2I_{\mathrm{p}}}{F^2(t_{\mathrm{r}})-\bm{k}(t_{\mathrm{r}}) \cdot \bm{F}'(t_{\mathrm{r}})}}.
+```
+
+The ``\mathrm{Im}\ \Phi_{\mathrm{tun}}`` related to the ionization rate, in the SFA-AE, is
+```math
+\mathrm{Im}\ \Phi_{\mathrm{tun}} \approx -\frac13 \frac{[k^2(t_{\mathrm{r}})+2I_{\mathrm{p}}]^{3/2}}{\sqrt{F^2(t_{\mathrm{r}})-\bm{k}(t_{\mathrm{r}}) \cdot \bm{F}'(t_{\mathrm{r}})}},
+```
+and we obtain
+```math
+\begin{aligned}
+ \mathrm{d}W/\mathrm{d}\bm{p}
+ &= \lvert P_{\bm{p}}(t_{\mathrm{s}}) \rvert^2 \exp(2\ \mathrm{Im}\ \Phi_{\mathrm{tun}}) \\
+ &= \lvert P_{\bm{p}}(t_{\mathrm{s}}) \rvert^2 \exp \left[ -\frac23 \frac{(k_\perp^2+2I_{\mathrm{p}})^{3/2}}{\sqrt{F^2-\bm{k}_\perp \cdot \bm{F}'}} \right],
+\end{aligned}
+```
+where the ``\bm{k}_{\perp}`` denotes the transverse momentum at the tunneling exit, which is actually equivalent to ``\bm{k}(t_{\mathrm{r}})`` in the SFA-AE because the above saddle-point equation requires ``\bm{k}(t_{\mathrm{r}}) \cdot \bm{F}(t_{\mathrm{r}}) = 0``. We note that the initial momentum, ``\bm{k}_0``, is exactly ``\bm{k}_{\perp}``.
+The prefactor ``P_{\bm{p}}(t_{\mathrm{s}})`` which included the Coulomb tail correction, in the SFA-AE, has the expression [^Frolov_2017]:
+```math
+P_{\bm{p}}(t_{\mathrm{s}}) = \left[ (k_\perp^2+2I_{\mathrm{p}})(F^2-\bm{k}_\perp \cdot \bm{F}') \right]^{-\alpha/4}.
+```
+
+The initial position has the expression
+```math
+\bm{r}_0 = \mathrm{Im} \int_{0}^{t_{\mathrm{i}}} \bm{A}(t_{\mathrm{r}}+\mathrm{i}\tau) \mathrm{d}\tau = -\frac{\bm{F}}{2} \frac{k_\perp^2+2I_{\mathrm{p}}}{F^2-\bm{k}_\perp \cdot \bm{F}'}.
+```
+
+[^Ni_2018]: H. Ni *et al.*, Tunneling Criteria and a Nonadiabatic Term for Strong-Field Ionization. *Phys. Rev. A* **98**, 013411 (2018). DOI:[10.1103/PhysRevA.98.013411](https://dx.doi.org/10.1103/PhysRevA.98.013411)
+[^Frolov_2017]: Frolov *et al.*, Adiabatic-Limit Coulomb Factors for Photoelectron and High-Order-Harmonic Spectra. *Phys. Rev. A* **96**, 023406 (2017). DOI:[10.1103/PhysRevA.96.023406](https://dx.doi.org/10.1103/PhysRevA.96.023406)
+
+
+
+## [Ammosov-Delone-Krainov (ADK)](@id ADK)
+
+The Ammosov-Delone-Krainov (ADK) theory [^Ammosov_1986] [^Delone_1998] is used to study the adiabatic tunneling in the strong-field ionization, and is, in a sense, the adiabatic limit of the SFA.
+
+In the adiabatic limit, the laser field can be treated as static, thus we have ``\bm{F}'(t)=\bm{0}`` (higher order derivatives of ``\bm{F}(t)`` remains zero as well).
+Substuting it into the expressions of SFA-AE yields the ADK rate
+```math
+\mathrm{d}W/\mathrm{d}\bm{p} = \lvert P_{\bm{p}}(t_{\mathrm{s}}) \rvert^2 \exp \left[ -\frac23 \frac{[k_\perp^2+2I_{\mathrm{p}}]^{3/2}}{F} \right],
+```
+where the prefactor reads
+```math
+\lvert P_{\bm{p}}(t_{\mathrm{s}}) \rvert^2 = \left[ (k_\perp^2+2I_{\mathrm{p}})F^2\right]^{-\alpha/2}.
+```
+The tunneling exit position can be obtained in the same approach:
+```math
+\bm{r}_0 = \mathrm{Im} \int_{0}^{t_{\mathrm{i}}} \bm{A}(t_{\mathrm{r}}+\mathrm{i}\tau) \mathrm{d}\tau = -\frac{\bm{F}}{2} \frac{k_\perp^2+2I_{\mathrm{p}}}{F^2}.
+```
+
+[^Ammosov_1986]: M. V. Ammosov *et al.*, Tunnel Ionization of Complex Atoms and of Atomic Ions in an Alternating Electromagnetic Field. *Sov. Phys. JETP* **64**, 1191 (1986).
+[^Delone_1998]: N. B. Delone *et al.*, Tunneling and Barrier-Suppression Ionization of Atoms and Ions in a Laser Radiation Field. *Phys.-Usp.* **41**, 469–485. DOI: [10.1070/PU1998v041n05ABEH000393](https://dx.doi.org/10.1070/PU1998v041n05ABEH000393)
+
+### [Tunneling Exit Methods for Atomic ADK](@id tun_exit_atomic_adk)
+
+The tunneling exit, i.e., the initial position ``\bm{r}_0`` of the classical trajectories, remains a controversial problem.
+Below we briefly review three methods for determination of the tunneling exit, which apply for the ADK method.
+
+(1) The ``I_{\mathrm{p}}/F`` model.
+
+The ``I_{\mathrm{p}}/F`` model is initially derived from the SFA theory and is the adiabatic limit of it.
+In this model the distance ``r_0`` between the tunneling exit and the nucleus is simply expressed as the quotient of the ionization potential ``I_{\mathrm{p}}`` and the electric field strength ``F``.
+Taking account of the initial kinetic energy ``k_\perp^2/2``, ``r_0`` reads
+```math
+r_0 = \frac{I_{\mathrm{p}}+k_\perp^2/2}{F}.
+```
+However, the shortcoming of the ``I_{\mathrm{p}}/F`` model is obvious: derived from the SFA which assumes a short-range potential (i.e., a delta potential), this model neglects the complexity of the Coulomb characteristics of the potential and thus the conclusion lacks persuasiveness.
+
+(2) The field-direction model (FDM).
+
+The field-direction model (FDM) determines the tunneling exit by treating the problem as a one-dimensional tunneling problem: the tunneling exit is found in the one-dimensional cut of the combined potential of the laser and nucleus along the field direction.
+We assume that the field points towards the ``-z`` direction, and the tunneling exit ``r_0`` is found from
+```math
+V(r_0 \bm{e}_z) - F r_0 = - I_{\mathrm{p}}.
+```
+For sufficiently large ``r_0``, the potential ``V(r)`` behaves as ``-Z/r``, and the tunneling exit is expressed as
+```math
+r_0 = \frac{I_{\mathrm{p}}+\sqrt{I_{\mathrm{p}}^2-4FZ}}{2F}.
+```
+The FDM approach takes the Coulomb field into account, however, not in such a scientific way because this problem is actually not independent of the transverse dimensions.
+
+(3) The parabolic-coordinate model.
+
+For a static electric field towards the ``-z`` direction, the parabolic-coordinate model introduces the parabolic coordinate ``\xi = r+z``, ``\eta = r-z``, ``\phi = \tan^{-1}{y/x}``, in which the Schrödinger equation becomes separable in the three dimensions.
+The parabolic-coordinate approach gives expression of the tunneling exit:
+```math
+r_0 = \frac{I_{\mathrm{p}}+\sqrt{I_{\mathrm{p}}^2-4F \left[ Z-(1+|m|)\sqrt{I_{\mathrm{p}}/2} \right] }}{2F},
+```
+where ``m`` is the magnetic quantum number along the ``z`` axis (due to the isotropic angular distribution of the atomic wavefunction, we set ``m=0``).
+
+
+## [Molecular ADK (MO-ADK)](@id MOADK)
+
+The molecular ADK (MO-ADK) theory generalizes the original ADK theory by extending the application scope from atomic to simple linear molecules [^Tong_2002].
+
+In the MO-ADK theory, the wavefunction of a linear molecule's ionizing orbital behaves asymptotically as
+```math
+\psi_0(\bm{r}) \sim \sum_{lm} F_{lm}(r) Y_{lm}(\theta,\phi)
+```
+in the molecular frame (MF) when ``r\rightarrow\infty``.
+Assigning ``\kappa=\sqrt{2I_{\mathrm{p}}}``, the ``F_{lm}(r)`` has the following asymptotic behavior when ``r\rightarrow\infty``:
+```math
+F_{lm}(r) \sim C_{lm} r^{Z/\kappa-1} \mathrm{e}^{-\kappa r}.
+```
+In numerical implementation we obtain the parameters ``C_{lm}`` by fitting the above expression [^Zhang_2015], and the ``F_{lm}(r)`` is found by the spherical-harmonic expansion of the wavefunction:
+```math
+F_{lm}(r) = \int \mathrm{d}\bm{\Omega} Y_{lm}^{*}(\bm{\Omega}) \psi_0(\bm{r}).
+```
+
+[^Zhang_2015]: Zhang, B. *et al.*, SLIMP: Strong Laser Interaction Model Package for Atoms and Molecules. *Comp. Phys. Comm.* **192**, 330–341 (2015). DOI: [10.1016/j.cpc.2015.02.031](https://dx.doi.org/10.1016/j.cpc.2015.02.031)
+
+We assume the electric field is pointing towards the ``z`` axis in the laboratory frame (LF).
+The angle-dependent tunneling ionization rate in the MO-ADK theory reads
+```math
+\Gamma(\beta,\gamma) = \mathrm{d}W/\mathrm{d}t = \sum_{m'} \frac{|B_{m'}(\beta,\gamma)|^2}{2^{|m'|}|m'|!} \kappa^{-|m'|} \left(\frac{2\kappa^2}{F}\right)^{2Z/\kappa-|m'|-1} \mathrm{e}^{-2\kappa^3/3F},
+```
+where the molecule's orientation is described using a set of Euler angles ``\hat{\bm{R}} = (\alpha,\beta,\gamma)`` (``z-y'-z''`` convention), which represents the rotational transformation from the MF to the LF; ``B_{m'}(\beta,\gamma)`` are the structural parameters which depend on the molecule's orbital wavefunction (here we omitted the ``\alpha`` dependence because the structural parameters are independent of ``\alpha``).
+The structural parameters ``B_{m'}(\beta,\gamma)`` have the following expression:
+```math
+B_{m'}(\beta,\gamma) = \sum_{lm} C_{lm} d_{m' m}^{l}(\beta) \mathrm{e}^{-\mathrm{i}m\gamma} Q_{l m'},
+```
+with ``d_{m' m}^{l}(\beta)`` being the Wigner-``d`` rotation matrix, and
+```math
+Q_{l m'} = (-1)^{(m'+|m'|)/2} \sqrt{\frac{(2l+1)(l+|m'|)!}{2(l-|m'|)!}}.
+```
+
+To ultilize the MO-ADK theory to provide the initial conditions in the trajectory simulation, we simply adopt the result of the atomic ADK theory:
+```math
+\bm{r}_0 = -\frac{\bm{F}}{2} \frac{k_\perp^2+2I_{\mathrm{p}}}{F^2}.
+```
+As for the ionization probability, we include the influence of the initial kinetic energy ``k_\perp^2/2`` by replacing the ``\kappa=\sqrt{2I_{\mathrm{p}}}`` with ``\kappa'(k_\perp)=\sqrt{2I_{\mathrm{p}}+k_\perp^2}`` in the exponential term of the ionization probability in the MO-ADK theory, giving
+```math
+\mathrm{d}W/\mathrm{d}\bm{k}_\perp \mathrm{d}t = \sum_{m'} \frac{|B_{m'}(\beta,\gamma)|^2}{2^{|m'|}|m'|!} \kappa^{-|m'|} \left(\frac{2\kappa^2}{F}\right)^{2Z/\kappa-|m'|-1} \mathrm{e}^{-2\kappa'^3(k_\perp)/3F}.
+```
+
+[^Tong_2002]: X. M. Tong *et al.*, Theory of Molecular Tunneling Ionization. *Phys. Rev. A* **66**, 033402 (2002). DOI: [10.1103/PhysRevA.66.033402](https://dx.doi.org/10.1103/PhysRevA.66.033402)
+
+
+
+## [Weak-Field Asymptotic Theory (WFAT)](@id WFAT)
+
+The weak-field asymptotic theory (WFAT) generalizes the tunneling ionization from isotropic atomic potentials to arbitrary molecular potentials [^Tolstikhin_2011]. Compared with the MO-ADK theory, the WFAT accounts for the influence of the molecules' permanent dipole moment, and is applicable for complex molecules other than simple linear molecules.
+
+The formulation of the WFAT is based on the expansion in the parabolic coordinates.
+The total ionization rate ``\Gamma(\beta,\gamma) = \mathrm{d}W/\mathrm{d}t``, is split into different parabolic channels:
+```math
+\Gamma(\beta,\gamma) = \sum_\nu \Gamma_\nu(\beta,\gamma),
+```
+where ``\Gamma_\nu(\beta,\gamma)`` are partial rates of parabolic quantum number indices ``\nu=(n_\xi,m)``, and ``n_\xi=0,1,2,\cdots``, ``m=0,\pm 1,\pm 2,\cdots``.
+In the leading-order approximation of the WFAT, the partial rates can be separated into two factors, namely the structural part ``|G_\nu(\beta,\gamma)|^2`` and the field part ``W_\nu(F)``:
+```math
+\Gamma_\nu(\beta,\gamma) = |G_\nu(\beta,\gamma)|^2 W_\nu(F).
+```
+The field factor is expressed as
+```math
+W_\nu(F) = \frac{\kappa}{2} \left(\frac{4\kappa^2}{F}\right)^{2Z/\kappa-2n_\xi-|m|-1} \mathrm{e}^{-2\kappa^3/3F}.
+```
+The structure factor, in the integral representation of the WFAT [^Dnestryan_2018], is given as an integral:
+```math
+G_\nu (\beta,\gamma) = \mathrm{e}^{-\kappa\mu_z} \int \Omega_\nu^* \left(\hat{\bm{R}}^{-1} \bm{r}\right) \hat{V}_{\mathrm{c}}(\bm{r}) \psi_0(\bm{r}) \mathrm{d} \bm{r},
+```
+where ``\psi_0`` is the wavefunction of the ionizing orbital,
+```math
+\bm{\mu} = \int \psi_0^*(\bm{r}) \bm{r} \psi_0(\bm{r}) \mathrm{d} \bm{r}
+```
+denotes the orbital dipole moment in the LF, with ``\mu_z`` being its component along the field direction;
+```math
+\Omega_\nu(\bm{r}) = \sum_{l=|m|}^{\infty} \Omega^\nu_{lm}(\bm{r}) = \sum_{l=|m|}^{\infty} R_l^\nu(r) Y_{lm}(\theta, \phi)
+```
+is a reference function which can be expanded into spherical harmonics, its radial part is expressed as
+```math
+R_l^\nu(r)=\omega_l^\nu \ (\kappa r)^l \ \mathrm{e}^{-\kappa r} \ \mathrm{M}(l+1-Z/\kappa, 2l+2, 2 \kappa r),
+```
+with ``\mathrm{M}(a,b,x)`` being the confluent hyper-geometric function and
+```math
+\begin{aligned}
+ \omega_l^\nu = & \ (-1)^{l+(|m|-m)/2+1}\ 2^{l+3/2}\ \kappa^{Z/\kappa-(|m|+1)/2-n_\xi} \\
+ & \times \sqrt{(2l+1)(l+m)!(l-m)!(|m|+n_\xi)!n_\xi!}\ \frac{l!}{(2l+1)!} \\
+ & \times \!\!\!\!\!\! \sum_{k=0}^{\min{(n_\xi,l-|m|)}} \!\!\!\!\!\!\!\!\!\! \frac{\Gamma(l+1-Z/\kappa+n_\xi-k)}{k!(l-k)!(|m|+k)!(l-|m|-k)!(n_\xi-k)!}
+\end{aligned}
+```
+the normalization coefficient;
+``\hat{V}_{\mathrm{c}}(\bm{r})=\hat{V}(\bm{r})+Z/r`` is the core potential with the Coulomb tail removed, where ``Z`` is the asymptotic charge of the residual ion.
+
+The effective potential ``\hat{V}(\bm{r})`` describes the interaction between the ionizing electron and the residual parent ion.
+Under the framework of the Hartree-Fock method, the effective potential consists of three parts, namely the nuclear Coulomb potential (``V_{\mathrm{nuc}}``), the direct (``V_{\mathrm{d}}``) and exchange (``V_{\mathrm{ex}}``) parts of inter-electron interactions:
+```math
+\hat{V}(\bm{r}) = V_{\mathrm{nuc}}(\bm{r}) + V_{\mathrm{d}}(\bm{r}) + \hat{V}_{\mathrm{ex}}(\bm{r}),
+```
+and
+```math
+\begin{aligned}
+ V_{\mathrm{nuc}}(\bm{r}) &= -\sum_{A=1}^{N_{\mathrm{atm}}} \frac{Z_A}{\left|\bm{r}-\bm{R}_A\right|}, \\
+ V_{\mathrm{d}}(\bm{r}) &= \sum_{i=1}^N \int \frac{\psi_i^*(\bm{r}') \psi_i(\bm{r}')}{|\bm{r}-\bm{r}'|} \mathrm{d} \bm{r}', \\
+ \hat{V}_{\mathrm{ex}}(\bm{r}) \psi_0(\bm{r}) &= -\sum_{i=1}^N \psi_i(\bm{r}) \int \frac{\psi_i^*(\bm{r}') \psi_0(\bm{r}')}{|\bm{r}-\bm{r}'|} \braket{\sigma_i | \sigma_0} \mathrm{d} \bm{r}',
+\end{aligned}
+```
+where ``N`` is the number of electrons, ``N_{\mathrm{atm}}`` is the number of nuclei, ``\psi_i(\bm{r})`` and ``\sigma_i`` denote the molecular orbital and the spin state of the electron of index ``i`` (``\braket{\sigma_i|\sigma_j}=1`` if electrons of index ``i`` and ``j`` have the same spin, and ``\braket{\sigma_i|\sigma_j}=0`` otherwise), ``Z_A`` and ``\bm{R}_A`` are the nuclear charge and position of atom of index ``A``.
+
+As the WFAT provides only the ionization rate ``\Gamma = \mathrm{d}W/\mathrm{d}t`` as the MO-ADK does, we adopt the same procedure as we did in the MO-ADK theory to provide the initial conditions for the trajectory simulation.
+The initial position ``\bm{r}_0`` is the same as that in the MO-ADK theory, and the ionization rate reads
+```math
+\mathrm{d}W/\mathrm{d}\bm{k}_\perp \mathrm{d}t = \sum_\nu |G_\nu(\beta,\gamma)|^2 \cdot \frac{\kappa}{2} \left(\frac{4\kappa^2}{F}\right)^{2Z/\kappa-2n_\xi-|m|-1} \mathrm{e}^{-2\kappa'^3(k_\perp)/3F}.
+```
+
+[^Tolstikhin_2011]: O. I. Tolstikhin *et al.*, Theory of Tunneling Ionization of Molecules: Weak-Field Asymptotics Including Dipole Effects. *Phys. Rev. A* **84**, 053423 (2011). DOI: [10.1103/PhysRevA.84.053423](https://dx.doi.org/10.1103/PhysRevA.84.053423)
+[^Dnestryan_2018]: A. I. Dnestryan *et al.*, Structure Factors for Tunneling Ionization Rates of Molecules: General Grid-Based Methodology and Convergence Studies. *J. Chem. Phys.* **149**, 164107. DOI: [10.1063/1.5046902](https://dx.doi.org/10.1063/1.5046902)
+
diff --git a/docs/src/theory2_trajectory_simulation_phase_methods.md b/docs/src/theory2_trajectory_simulation_phase_methods.md
new file mode 100644
index 0000000..fed411e
--- /dev/null
+++ b/docs/src/theory2_trajectory_simulation_phase_methods.md
@@ -0,0 +1,116 @@
+# [Theory: Trajectory Simulation and Phase Methods](@id theory_traj_phase)
+
+*This section reviews the trajectory simulation procedure and the phase methods within.*
+
+Given the initial conditions, the tunneled electrons evolve classically in the combined field of Coulomb and laser, following a classical trajectory, the scheme is named after the *Classical Trajectory Monte Carlo (CTMC)*.
+Apart from the position and momentum, phase methods like the *Quantum Trajectory Monte Carlo (QTMC)* and *Semiclassical Two-Step (SCTS) Model* give an additional quantum phase property to the classical trajectories, which are capable of reproducing more details in the final momentum spectrum than the full-classical CTMC.
+
+In the following we review the scheme of trajectory simulation and introduce the quantum phase methods available at present.
+Note that atomic units (a.u.) is used throughout.
+
+```@contents
+Pages = ["theory2_trajectory_simulation_phase_methods.md"]
+Depth = 3
+```
+
+
+## [Classical Trajectory Monte Carlo (CTMC)](@id CTMC)
+
+In the CTMC, each sample electron carries a probability ``W``, following a classical trajectory, and finally ends up with a final momentum ``\bm{p}_\infty = \bm{p}|_{t=\infty}``, which is our interested physical quantity.
+
+The tunneled electrons, each having different tunneling time, initial positions and momenta, evolve under the Hamiltonian equation:
+```math
+\dot{\bm{r}} = \bm{\nabla}_{\bm{p}}H, \quad \dot{\bm{p}} = -\bm{\nabla}_{\bm{r}}H.
+```
+The Hamiltonian reads
+```math
+H = \frac12 \left[ \bm{p}+\bm{A}(t) \right]^2 + V(\bm{r}),
+```
+where ``V(\bm{r})`` denotes the potential of the parent ion.
+
+After the laser ends, the electron interacts only with the residual parent ion.
+At a distance from the parent ion, the electron interacts with the potential's Coulomb tail, and its Runge-Lenz vector ``\bm{a} = \bm{p}\times\bm{L} - Z\bm{r}/r`` can be viewed as approximately conserved. Taking advantage of this conserved quantity, combining with the conservation of angular momentum and energy, we obtain the expression of the final momentum [^ShvetsovShilovski_2012]:
+```math
+\begin{aligned}
+ \bm{p}_\infty &= p_\infty \frac{p_\infty(\bm{L}\times\bm{a})-\bm{a}}{1+p_\infty^2 L^2}, \\
+ p_\infty^2/2 &= p^2/2 - Z/r, \\
+ \bm{L} &= \bm{r}\times\bm{p}, \\
+ \bm{a} &= \bm{p}\times\bm{L} - Z\bm{r}/r, \\
+\end{aligned}
+```
+where ``\bm{r},\bm{p}`` are quantities of the electron at any time after the laser ends.
+This scheme applies for electrons with positive energy, which are able to finally escape the parent ion and reach the detector.
+For electrons with negative energy, we assume that they finally become Rydberg states.
+
+Finally, electrons with similar final momenta (i.e., in the same small box of the final momentum grid) would be collected by summing up the probabilities they carry: ``W_{\bm{p}} = \sum_i{W_i}``, and the final momentum spectrum is given by ``W_{\bm{p}}``.
+
+[^ShvetsovShilovski_2012]: N. I. Shvetsov-Shilovski *et al.*, Ionization in elliptically polarized pulses: Multielectron polarization effects and asymmetry of photoelectron momentum distributions, *Phys. Rev. A* **85**, 023428 (2012). DOI: [10.1103/PhysRevA.85.023428](https://dx.doi.org/10.1103/PhysRevA.85.023428)
+
+
+### [Non-dipole Effects on the Trajectory Motion](@id traj_nondipole)
+
+Dipole approximation is usually applied in the study of laser-matter interaction by neglecting the spatial dependence of the laser field, i.e., we let ``\bm{A}(\bm{r},t) = \bm{A}(t)``.
+However, for lasers with high intensity or frequency, the spatial dependence of the laser becomes noticeable, the dipole approximation breaks down, and we have to take non-dipole effects into account.
+
+To include the first-order non-dipole effects in trajectory simulation, we first refer to the spatial dependent vector potential ``\bm{A}(\bm{r},t)``, giving its first-order expansion in space coordinates:
+```math
+\bm{A}(\bm{r},t) = \bm{A}(t) \mathrm{e}^{\mathrm{i}\bm{k}\cdot\bm{r}} \approx \bm{A}(t) + \frac{z}{c} \bm{F}(t),
+```
+where we assume that the laser propagates in ``z`` direction.
+In this way we obtained the Hamiltonian which includes the first-order non-dipole effects:
+```math
+H = \frac12 \left[ \bm{p}+\bm{A}(t)+\frac{z}{c}\bm{F}(t) \right]^2 + V(\bm{r}).
+```
+
+
+## [Quantum Trajectory Monte Carlo (QTMC)](@id QTMC)
+
+Compared with the CTMC, the QTMC scheme endows each electron trajectory with a quantum phase ``\Phi`` based on the Feynman path-integral approach [^Li_2014].
+The phase gets acculmulated during the electron's excursion and is expressed as
+```math
+\Phi = - \int_{t_0}^\infty \left[ \frac{k^2}{2} + V(\bm{r}) + I_{\mathrm{p}} \right] \mathrm{d}t.
+```
+where ``t_0`` is the time when the electron tunneled, and ``\bm{k}=\dot{\bm{r}}`` denotes the momentum. It's worth noting that in the velocity gauge, we distinguish the actual momentum/velocity ``\bm{k}`` from ``\bm{p}=\bm{k}-\bm{A}`` which denotes the canonical momentum when the laser is present, and use ``\bm{p}`` for the momentum after the laser ends (the vector potential ``\bm{A}`` vanishes and we have ``\bm{p}=\bm{k}``).
+Finally the momentum spectrum is given by coherently summing up the probability amplitude, and taking the square modulus of the summation result:
+```math
+W_{\bm{p}} = \left| \sum_i \sqrt{W_i}\ \mathrm{e}^{\mathrm{i}\Phi_i} \right|^2.
+```
+
+It's also worthwhile noting that in practical implementation, the upper limit of the integral of the quantum phase ``\Phi`` doesn't have to be infinity.
+Since electrons which arrived at the same final momentum share the same energy after the laser ends (at ``t_{\mathrm{f}}``), the integral
+```math
+\int_{t_{\mathrm{f}}}^\infty \left[ \frac{k^2}{2} + V(\bm{r}) + I_{\mathrm{p}} \right] \mathrm{d}t
+```
+is same for electrons with the same final momentum.
+Therefore, in numerical implementation, the upper limit of the phase integral can be simply set as the end of the laser, i.e., the ``t_{\mathrm{f}}``.
+
+[^Li_2014]: M. Li *et al.*, Classical-quantum correspondence for above-threshold ionization, *Phys. Rev. Lett.* **112**, 113002 (2014). DOI: [10.1103/PhysRevLett.112.113002](https://dx.doi.org/10.1103/PhysRevLett.112.113002)
+
+
+## [Semiclassical Two-Step (SCTS) Model](@id SCTS)
+
+The SCTS model [^ShvetsovShilovski_2016] improves the quantum phase in the QTMC scheme, giving
+```math
+\Phi = - \bm{k}_0\cdot\bm{r}_0 - \int_{t_0}^\infty \left[ \frac{k^2}{2} + V(\bm{r}) - \bm{r}\cdot\bm{\nabla}V(\bm{r}) + I_{\mathrm{p}} \right] \mathrm{d}t.
+```
+The difference between the SCTS phase and the QTMC lies in two aspects:
+The first is the initial phase ``\bm{k}_0\cdot\bm{r}_0``, which is non-zero for non-zero initial longitudinal momentum ``k_\parallel`` w.r.t. the non-adiabatic tunneling process.
+The second is the ``\bm{r}\cdot\bm{\nabla}V(\bm{r})`` term in the integrand which is omitted in the QTMC scheme.
+
+For the SCTS model, the phase integral in the interval ``[t_{\mathrm{f}},\infty)`` cannot be simply neglected due to the presence of the ``\bm{r}\cdot\bm{\nabla}V(\bm{r})`` term in the integrand.
+However, the integral of this term can be reduced to an analytical expression in case of Coulomb potential (``V(r)=Z/r``):
+```math
+\begin{aligned}
+ \Phi_{\mathrm{f}}^{\mathrm{C}}(t_{\mathrm{f}})
+ &= \int_{t_{\mathrm{f}}}^\infty \bm{r}\cdot\bm{\nabla}V(\bm{r}) \mathrm{d}t \\
+ &= Z \int_{t_{\mathrm{f}}}^\infty \frac{\mathrm{d}t}{r} \\
+ &= - \frac{Z}{\kappa} \left[ \ln{g} + \sinh^{-1} \left( \frac{\kappa}{g}\bm{r}_{\mathrm{f}}\cdot\bm{p}_{\mathrm{f}} \right) \right],
+\end{aligned}
+```
+where ``\bm{r}_{\mathrm{f}}=\bm{r}(t_{\mathrm{f}})``, ``\bm{p}_{\mathrm{f}}=\bm{p}(t_{\mathrm{f}})`` and ``g = \sqrt{1+2\kappa^2 L^2} = \sqrt{1+2\kappa^2 (\bm{r}_{\mathrm{f}}\times\bm{p}_{\mathrm{f}})^2}``.
+In this way we obtain the expression of the SCTS phase that is suitable for numerical implementation:
+```math
+\Phi = - \bm{k}_0\cdot\bm{r}_0 + I_{\mathrm{p}}t_0 - \int_{t_0}^{t_{\mathrm{f}}} \left[ \frac{k^2}{2} + V(\bm{r}) - \bm{r}\cdot\bm{\nabla}V(\bm{r}) \right] \mathrm{d}t + \Phi_{\mathrm{f}}^{\mathrm{C}}(t_{\mathrm{f}}).
+```
+
+[^ShvetsovShilovski_2016]: N. I. Shvetsov-Shilovski *et al.*, Semiclassical two-step model for strong-field ionization, *Phys. Rev. A* **94**, 013415 (2016). DOI: [10.1103/PhysRevA.94.013415](https://dx.doi.org/10.1103/PhysRevA.94.013415)
diff --git a/examples/example_attoclock.jl b/examples/example_attoclock.jl
new file mode 100644
index 0000000..60ad13b
--- /dev/null
+++ b/examples/example_attoclock.jl
@@ -0,0 +1,44 @@
+# === example_attoclock.jl ===
+#
+# With the aim of studying the influence of non-adiabatic effects on the attoclock signal, we employ ADK, SFA-AE and SFA to provide the initial conditions of the electrons, and perform trajectory simulations.
+# The parameters are taken from [J. Phys. B 54, 144001 (2020)].
+#
+# To run the example, a julia intepreter (version ≥ 1.7) with the *SemiclassicalSFI* package (and its dependencies) installed is required.
+#
+# Usage: place example_attoclock.jl in the working directory, and execute:
+# julia -t example_attoclock.jl
+# where the `n_threads` should be replaced with the desired number of threads allocated for julia.
+
+
+@info "SemiclassicalSFI Example - Attoclock"
+
+@info "Loading Packages..."
+using SemiclassicalSFI
+using SemiclassicalSFI.Targets
+using SemiclassicalSFI.Lasers
+
+t = HeAtom()
+l = Cos4Laser(peak_int=4e14, wave_len=800, cyc_num=2, ellip=1.0)
+
+init_cond_list = [:ADK, :SFAAE, :SFA]
+
+for init_cond in init_cond_list
+ @info "Running $(init_cond)..."
+ filename = "SCSFI_He_4e14_800nm_2cyc_CP_$(init_cond)_CTMC.h5"
+ performSFI(
+ target = t,
+ laser = l,
+ init_cond_method = init_cond,
+ sample_t_intv = (-80,80),
+ sample_t_num = 20000,
+ ss_kd_max = 1.5,
+ ss_kd_num = 1000,
+ ss_kz_max = 1.5,
+ ss_kz_num = 50,
+ traj_phase_method = :CTMC,
+ traj_t_final = 120,
+ final_p_max = (2,2,2),
+ final_p_num = (500,500,1),
+ save_path = filename
+ )
+end
diff --git a/examples/example_attoclock_plot.jl b/examples/example_attoclock_plot.jl
new file mode 100644
index 0000000..ffe7719
--- /dev/null
+++ b/examples/example_attoclock_plot.jl
@@ -0,0 +1,32 @@
+# === example_attoclock_plot.jl ===
+#
+# To correctly generate the plots using the code in the example, the HDF5, Colors, ColorSchemes, Plots and PyPlot (which is successfully precompiled) packages are also required to be installed via the package manager of julia.
+
+
+@info "Plotting..."
+using HDF5
+using Plots
+pyplot() # PyPlot backend, need to install the PyPlot package manually. Switching to other backends is okay but the output might be slightly different.
+include("ijet.jl") # modified jet colormap which has white color for the bottom value.
+clim = (-5,0)
+
+function plot_heatmap(h5file; kwargs...)
+ ion_prob = read(h5file, "ion_prob")
+ px = read(h5file, "px")
+ py = read(h5file, "py")
+ spec = read(h5file, "momentum_spec_2D")
+ spec ./= ion_prob
+ @. spec = max(1e-100, spec) # remove zero value because the figure is displayed in logarithmic scale.
+ return heatmap(px,py,log10.(spec'), xlim=(-2,2), ylim=(-2,2), c=:ijet, clim=clim,
+ aspectratio=:equal, tick_direction=:out, bordercolor=:black, framestyle=:box; kwargs...)
+end
+
+figure_list = Vector{Any}(nothing,3)
+figure_list[1] = plot_heatmap(h5open("SCSFI_He_4e14_800nm_2cyc_CP_ADK_CTMC.h5"); title="ADK", xlabel=raw"$p_x$ (a.u.)", ylabel=raw"$p_y$ (a.u.)", legend=:none)
+figure_list[2] = plot_heatmap(h5open("SCSFI_He_4e14_800nm_2cyc_CP_SFAAE_CTMC.h5"); title="SFAAE", xlabel=raw"$p_x$ (a.u.)", legend=:none)
+figure_list[3] = plot_heatmap(h5open("SCSFI_He_4e14_800nm_2cyc_CP_SFA_CTMC.h5"); title="SFA", xlabel=raw"$p_x$ (a.u.)")
+
+l = @layout [a{0.32w} b{0.32w} c{0.36w}]
+plot(figure_list..., size=(900,300), layout=l, fmt=:svg)
+savefig("example_attoclock.pdf")
+@info "Plot saved as \"example_attoclock.pdf\"."
\ No newline at end of file
diff --git a/examples/example_phase_methods_QTMC.jl b/examples/example_phase_methods_QTMC.jl
new file mode 100644
index 0000000..e30cb42
--- /dev/null
+++ b/examples/example_phase_methods_QTMC.jl
@@ -0,0 +1,48 @@
+# === example_phase_methods_QTMC.jl ===
+#
+# To investigate the difference between the phase methods QTMC and SCTS, we use an 8-cycle linearly polarized laser pulse, employing the QTMC and SCTS as the phase methods during the trajectory simulation.
+# The parameters are taken from [Phys. Rev. A 94, 013415 (2016)].
+# This is the QTMC task.
+#
+# To run the example, a julia intepreter (version ≥ 1.7) with the *SemiclassicalSFI* package (and its dependencies) installed is required.
+#
+# Usage: place example_phase_methods_QTMC.jl in the working directory, and execute:
+# julia -t example_phase_methods_QTMC.jl
+# where the `n_threads` should be replaced with the desired number of threads allocated for julia.
+#
+# Note:
+# 1. This example takes a very long time (~10h each task for 4 threads) to finish.
+# 2. In this example, the 3D momentum spectrum is collected, which is very memory-consuming, consider using fewer threads.
+# 3. A 16-thread task is much less efficient than a 4-thread task and consumes much more memory. Therefore, it's suggested to use fewer threads (e.g., 4) per task.
+
+@info "SemiclassicalSFI Example - Phase Methods - QTMC"
+
+@info "Loading Packages..."
+using SemiclassicalSFI
+using SemiclassicalSFI.Targets
+using SemiclassicalSFI.Lasers
+
+t = HAtom()
+l = Cos4Laser(peak_int=0.9e14, wave_len=800, cyc_num=8, ellip=0.0)
+filename = "SCSFI_H_0.9e14_800nm_8cyc_LP_ADKPara_ExpRate_QTMC.h5"
+performSFI(
+ target = t,
+ laser = l,
+ init_cond_method = :ADK,
+ adk_tun_exit = :Para,
+ sample_t_intv = (-300,300),
+ sample_t_num = 50000,
+ sample_cutoff_limit = 1e-16,
+ ss_kd_max = 1.0,
+ ss_kd_num = 400,
+ ss_kz_max = 1.0,
+ ss_kz_num = 400,
+ rate_prefix = :ExpRate,
+ traj_phase_method = :QTMC,
+ traj_t_final = 450,
+ traj_rtol = 1e-6,
+ save_3D_spec = true,
+ final_p_max = (1,1,1),
+ final_p_num = (400,400,400),
+ save_path = filename
+ )
\ No newline at end of file
diff --git a/examples/example_phase_methods_SCTS.jl b/examples/example_phase_methods_SCTS.jl
new file mode 100644
index 0000000..b046803
--- /dev/null
+++ b/examples/example_phase_methods_SCTS.jl
@@ -0,0 +1,48 @@
+# === example_phase_methods_SCTS.jl ===
+#
+# To investigate the difference between the phase methods QTMC and SCTS, we emplopyed the QTMC and SCTS as the phase methods during the trajectory simulation.
+# The parameters are taken from [Phys. Rev. A 94, 013415 (2016)].
+# This is the SCTS task.
+#
+# To run the example, a julia intepreter (version ≥ 1.7) with the *SemiclassicalSFI* package (and its dependencies) installed is required.
+#
+# Usage: place example_phase_methods_SCTS.jl in the working directory, and execute:
+# julia -t example_phase_methods_SCTS.jl
+# where the `n_threads` should be replaced with the desired number of threads allocated for julia.
+#
+# Note:
+# 1. This example takes a very long time (~10h each task for 4 threads) to finish.
+# 2. In this example, the 3D momentum spectrum is collected, which is very memory-consuming, consider using fewer threads.
+# 3. A 16-thread task is much less efficient than a 4-thread task and consumes much more memory. Therefore, it's suggested to use fewer threads (e.g., 4) per task.
+
+@info "SemiclassicalSFI Example - Phase Methods - SCTS"
+
+@info "Loading Packages..."
+using SemiclassicalSFI
+using SemiclassicalSFI.Targets
+using SemiclassicalSFI.Lasers
+
+t = HAtom()
+l = Cos4Laser(peak_int=0.9e14, wave_len=800, cyc_num=8, ellip=0.0)
+filename = "SCSFI_H_0.9e14_800nm_8cyc_LP_ADKPara_ExpRate_SCTS.h5"
+performSFI(
+ target = t,
+ laser = l,
+ init_cond_method = :ADK,
+ adk_tun_exit = :Para,
+ sample_t_intv = (-300,300),
+ sample_t_num = 50000,
+ sample_cutoff_limit = 1e-16,
+ ss_kd_max = 1.0,
+ ss_kd_num = 400,
+ ss_kz_max = 1.0,
+ ss_kz_num = 400,
+ rate_prefix = :ExpRate,
+ traj_phase_method = :SCTS,
+ traj_t_final = 450,
+ traj_rtol = 1e-6,
+ save_3D_spec = true,
+ final_p_max = (1,1,1),
+ final_p_num = (400,400,400),
+ save_path = filename
+ )
\ No newline at end of file
diff --git a/examples/example_phase_methods_plot.jl b/examples/example_phase_methods_plot.jl
new file mode 100644
index 0000000..325b4f5
--- /dev/null
+++ b/examples/example_phase_methods_plot.jl
@@ -0,0 +1,33 @@
+# === example_phase_methods_plot.jl ===
+#
+# To correctly generate the plots using the code in the example, the HDF5, Colors, ColorSchemes, Plots and PyPlot (which is successfully precompiled) packages are also required to be installed via the package manager of julia.
+
+
+@info "Plotting..."
+using HDF5
+using Plots
+pyplot() # PyPlot backend, need to install the PyPlot package manually. Switching to other backends is okay but the output might be slightly different.
+include("ijet.jl") # modified jet colormap which has white color for the bottom value.
+clim = (2,6)
+
+function plot_heatmap(h5file; kwargs...)
+ ion_prob = read(h5file, "ion_prob")
+ px = read(h5file, "px")
+ py = read(h5file, "py")
+ pz = read(h5file, "pz")
+ spec = read(h5file, "momentum_spec_3D")
+ spec ./= ion_prob # normalize the total ionization yield
+ @. spec = max(1e-100, spec) # remove zero value because the figure is displayed in logarithmic scale.
+ # plots the slice on the pz-px plane near py=0.
+ return heatmap(pz,px,log10.(mapreduce(i->spec[:,i,:],+,length(py)÷2-5:length(py)÷2+5))', xlim=(-1,1), ylim=(-1,1), c=:ijet, clim=clim,
+ aspectratio=:equal, tick_direction=:out, bordercolor=:black, framestyle=:box; kwargs...)
+end
+
+figure_list = Vector{Any}(nothing,2)
+figure_list[1] = plot_heatmap(h5open("SCSFI_H_0.9e14_800nm_8cyc_LP_ADKPara_ExpRate_QTMC.h5"); title="QTMC", xlabel=raw"$p_z$ (a.u.)", ylabel=raw"$p_x$ (a.u.)", legend=:none)
+figure_list[2] = plot_heatmap(h5open("SCSFI_H_0.9e14_800nm_8cyc_LP_ADKPara_ExpRate_SCTS.h5"); title="SCTS", xlabel=raw"$p_z$ (a.u.)")
+
+l = @layout [a{0.45w} b{0.55w}]
+plot(figure_list..., size=(600,300), layout=l, fmt=:svg)
+savefig("example_phase_methods.pdf")
+@info "Plot saved as \"example_phase_methods.pdf\"."
\ No newline at end of file
diff --git a/examples/ijet.jl b/examples/ijet.jl
new file mode 100644
index 0000000..d8673ec
--- /dev/null
+++ b/examples/ijet.jl
@@ -0,0 +1,10 @@
+using Colors
+using ColorSchemes
+
+r = cat(range(1.0,0.0,32),zeros(64),range(0.0,1.0,64),ones(64),range(1.0,0.5,32); dims=1)
+g = cat(range(1.0,0.0,32),range(0.0,1.0,64),ones(64),range(1.0,0.0,64),zeros(32); dims=1)
+b = cat(ones(96),range(1.0,0.0,64),zeros(96); dims=1)
+
+RGB_list = map((r,g,b)->RGB(r,g,b), r,g,b)
+
+loadcolorscheme(:ijet, RGB_list)
\ No newline at end of file
diff --git a/src/Lasers/Cos2Laser.jl b/src/Lasers/Cos2Laser.jl
index b374db0..3b2a933 100644
--- a/src/Lasers/Cos2Laser.jl
+++ b/src/Lasers/Cos2Laser.jl
@@ -1,5 +1,27 @@
-"Represents a monochromatic elliptically polarized laser field with Cos2-shape envelope propagating in z direction."
+"""
+```
+struct Cos2Laser <: MonochromaticLaser
+```
+Represents a monochromatic elliptically polarized laser field with Cos2-shape envelope propagating in z direction.
+
+An instance of `Cos2Laser` can be initialized via the constructor method:
+```julia
+Cos2Laser(peak_int, wave_len|ang_freq, cyc_num|duration, ellip, azi=0.0, cep=0.0, t_shift=0.0)
+```
+
+# Parameters
+- `peak_int` : Peak intensity of the laser field (in W/cm²).
+- `wave_len` : Wave length of the laser field (in nm). Must specify either `wave_len` or `ang_freq`.
+- `ang_freq` : Angular frequency of the laser field (in a.u.). Must specify either `wave_len` or `ang_freq`.
+- `cyc_num` : Number of cycles of the laser field. Must specify either `cyc_num` or `duration`.
+- `duration` : Duration of the laser field (in a.u.). Must specify either `cyc_num` or `duration`.
+- `ellip` : Ellipticity of the laser field [-1≤ε≤1, 0 indicates linear polarization and ±1 indicates circular polarization].
+- `azi` : Azimuth angle of the laser's polarization's principle axis relative to x axis (in radians) (optional, default 0).
+- `cep` : Carrier-Envelope-Phase of the laser field (optional, default 0).
+- `t_shift` : Time shift of the laser (in a.u.) relative to the peak (optional, default 0).
+
+"""
struct Cos2Laser <: MonochromaticLaser
"Peak intensity of the laser field (in W/cm^2)."
peak_int;
diff --git a/src/Lasers/Cos4Laser.jl b/src/Lasers/Cos4Laser.jl
index d538b10..a4d2ebd 100644
--- a/src/Lasers/Cos4Laser.jl
+++ b/src/Lasers/Cos4Laser.jl
@@ -1,5 +1,27 @@
-"Represents a monochromatic elliptically polarized laser field with Cos4-shape envelope propagating in z direction."
+"""
+```
+struct Cos4Laser <: MonochromaticLaser
+```
+Represents a monochromatic elliptically polarized laser field with Cos4-shape envelope propagating in z direction.
+
+An instance of `Cos4Laser` can be initialized via the constructor method:
+```julia
+Cos4Laser(peak_int, wave_len|ang_freq, cyc_num|duration, ellip, azi=0.0, cep=0.0, t_shift=0.0)
+```
+
+# Parameters
+- `peak_int` : Peak intensity of the laser field (in W/cm²).
+- `wave_len` : Wave length of the laser field (in nm). Must specify either `wave_len` or `ang_freq`.
+- `ang_freq` : Angular frequency of the laser field (in a.u.). Must specify either `wave_len` or `ang_freq`.
+- `cyc_num` : Number of cycles of the laser field. Must specify either `cyc_num` or `duration`.
+- `duration` : Duration of the laser field (in a.u.). Must specify either `cyc_num` or `duration`.
+- `ellip` : Ellipticity of the laser field [-1≤ε≤1, 0 indicates linear polarization and ±1 indicates circular polarization].
+- `azi` : Azimuth angle of the laser's polarization's principle axis relative to x axis (in radians) (optional, default 0).
+- `cep` : Carrier-Envelope-Phase of the laser field (optional, default 0).
+- `t_shift` : Time shift of the laser (in a.u.) relative to the peak (optional, default 0).
+
+"""
struct Cos4Laser <: MonochromaticLaser
"Peak intensity of the laser field (in W/cm^2)."
peak_int;
diff --git a/src/Lasers/GaussianLaser.jl b/src/Lasers/GaussianLaser.jl
index 1974e7f..dd0812f 100644
--- a/src/Lasers/GaussianLaser.jl
+++ b/src/Lasers/GaussianLaser.jl
@@ -1,11 +1,34 @@
-"Represents a monochromatic elliptically polarized laser field with Gaussian-shape envelope propagating in z direction."
+"""
+```
+struct GaussianLaser <: MonochromaticLaser
+```
+Represents a monochromatic elliptically polarized laser field with Gaussian-shape envelope propagating in z direction.
+
+An instance of `GaussianLaser` can be initialized via the constructor method:
+```julia
+GaussianLaser(peak_int, wave_len|ang_freq, spread_cyc_num|spread_duration|FWHM_duration, ellip, azi=0., cep=0., t_shift=0.)
+```
+
+# Parameters
+- `peak_int` : Peak intensity of the laser field (in W/cm²).
+- `wave_len` : Wave length of the laser field (in nm). Must specify either `wave_len` or `ang_freq`.
+- `ang_freq` : Angular frequency of the laser field (in a.u.). Must specify either `wave_len` or `ang_freq`.
+- `spread_cyc_num` : Temporal width (converting to cycle numbers) of the laser field (in a.u.), namely σ. Must specify one in `spread_cyc_num`, `spread_duration` and `FWHM_duration`.
+- `spread_duration` : Temporal width of the laser field (in a.u.). Must specify one in `spread_cyc_num`, `spread_duration` and `FWHM_duration`.
+- `FWHM_duration` : Temporal FWHM (Full Width at Half Maximum) of the laser field (in a.u.). Must specify one in `spread_cyc_num`, `spread_duration` and `FWHM_duration`.
+- `ellip` : Ellipticity of the laser field [-1≤ε≤1, 0 indicates linear polarization and ±1 indicates circular polarization].
+- `azi` : Azimuth angle of the laser's polarization's principle axis relative to x axis (in radians) (optional, default 0).
+- `cep` : Carrier-Envelope-Phase of the laser field (optional, default 0).
+- `t_shift` : Time shift of the laser (in a.u.) relative to the peak (optional, default 0).
+
+"""
struct GaussianLaser <: MonochromaticLaser
"Peak intensity of the laser field (in W/cm^2)."
peak_int;
"Wavelength of the laser field (in NANOMETER)."
wave_len;
- "Temporal width (converting to cycle numbers) of the laser field, namely σ."
+ "Temporal width (converting to cycle numbers) of the laser field (in a.u.), namely σ."
spread_cyc_num;
"Ellipticity of the laser field."
ellip;
@@ -21,7 +44,7 @@ struct GaussianLaser <: MonochromaticLaser
# Parameters
- `peak_int` : Peak intensity of the laser field (in W/cm²).
- `wave_len` : Wavelength of the laser field (in nm).
- - `spread_cyc_num` : Temporal width (converting to cycle numbers) of the laser field, namely σ.
+ - `spread_cyc_num` : Temporal width (converting to cycle numbers) of the laser field (in a.u.), namely σ.
- `ellip` : Ellipticity of the laser field [-1≤e≤1, 0 indicates linear polarization and ±1 indicates circular polarization].
- `azi` : Azimuth angle of the laser's polarization's principle axis relative to x axis (in radians) (optional, default 0).
- `cep` : Carrier-Envelope-Phase of the laser field (optional, default 0).
@@ -41,9 +64,9 @@ struct GaussianLaser <: MonochromaticLaser
- `peak_int` : Peak intensity of the laser field (in W/cm²).
- `wave_len` : Wave length of the laser field (in nm). Must specify either `wave_len` or `ang_freq`.
- `ang_freq` : Angular frequency of the laser field (in a.u.). Must specify either `wave_len` or `ang_freq`.
- - `spread_cyc_num` : Temporal width (converting to cycle numbers) of the laser field, namely σ. Must specify one in `spread_cyc_num`, `spread_duration` and `FWHM_duration`.
+ - `spread_cyc_num` : Temporal width (converting to cycle numbers) of the laser field (in a.u.), namely σ. Must specify one in `spread_cyc_num`, `spread_duration` and `FWHM_duration`.
- `spread_duration` : Temporal width of the laser field (in a.u.). Must specify one in `spread_cyc_num`, `spread_duration` and `FWHM_duration`.
- - `FWHM_duration` : Temporal FWHM (Full Width at Half Maxima) of the laser field (in a.u.). Must specify one in `spread_cyc_num`, `spread_duration` and `FWHM_duration`.
+ - `FWHM_duration` : Temporal FWHM (Full Width at Half Maximum) of the laser field (in a.u.). Must specify one in `spread_cyc_num`, `spread_duration` and `FWHM_duration`.
- `ellip` : Ellipticity of the laser field [-1≤ε≤1, 0 indicates linear polarization and ±1 indicates circular polarization].
- `azi` : Azimuth angle of the laser's polarization's principle axis relative to x axis (in radians) (optional, default 0).
- `cep` : Carrier-Envelope-Phase of the laser field (optional, default 0).
diff --git a/src/Lasers/Lasers.jl b/src/Lasers/Lasers.jl
index 69a7a1f..9fb6edd 100644
--- a/src/Lasers/Lasers.jl
+++ b/src/Lasers/Lasers.jl
@@ -1,6 +1,6 @@
"""
-The Laser module provides information about laser fields.
+The `Lasers` module provides information about laser fields.
"""
module Lasers
diff --git a/src/Lasers/TrapezoidalLaser.jl b/src/Lasers/TrapezoidalLaser.jl
index b01edc7..9bcc93d 100644
--- a/src/Lasers/TrapezoidalLaser.jl
+++ b/src/Lasers/TrapezoidalLaser.jl
@@ -1,5 +1,27 @@
-"Represents a monochromatic elliptically polarized laser field with Trapezoidal-shape envelope propagating in z direction."
+"""
+```
+struct TrapezoidalLaser <: MonochromaticLaser
+```
+Represents a monochromatic elliptically polarized laser field with Trapezoidal-shape envelope propagating in z direction.
+
+An instance of `TrapezoidalLaser` can be initialized via the constructor method:
+```julia
+TrapezoidalLaser(peak_int, wave_len|ang_freq, cyc_num_turn_on, cyc_num_turn_off, cyc_num_const, ellip, azi=0.0, cep=0.0, t_shift=0.0)
+```
+
+# Parameters
+- `peak_int` : Peak intensity of the laser field (in W/cm²).
+- `wave_len` : Wavelength of the laser field (in nm). Must specify either `wave_len` or `ang_freq`.
+- `ang_freq` : Angular frequency of the laser field (in a.u.). Must specify either `wave_len` or `ang_freq`.
+- `cyc_num_turn_on` : Number of cycles of the laser field in the turn-on stage.
+- `cyc_num_turn_off`: Number of cycles of the laser field in the turn-off stage.
+- `cyc_num_const` : Number of cycles of the laser field in the constant-intensity stage.
+- `ellip` : Ellipticity of the laser field [-1≤ε≤1, 0 indicates linear polarization and ±1 indicates circular polarization].
+- `azi` : Azimuth angle of the laser's polarization's principle axis relative to x axis (in radians) (optional, default 0).
+- `t_shift` : Time shift of the laser (in a.u.) relative to the beginning of TURN-ON (optional, default 0).
+
+"""
struct TrapezoidalLaser <: MonochromaticLaser
"Peak intensity of the laser field (in W/cm^2)."
peak_int;
@@ -24,9 +46,9 @@ struct TrapezoidalLaser <: MonochromaticLaser
# Parameters
- `peak_int` : Peak intensity of the laser field (in W/cm²).
- `wave_len` : Wavelength of the laser field (in nm).
- - `cyc_num_turn_on` : Number of cycles of the laser field in the turn-on.
- - `cyc_num_turn_off`: Number of cycles of the laser field in the turn-off.
- - `cyc_num_const` : Number of cycles of the laser field in the constant-intensity.
+ - `cyc_num_turn_on` : Number of cycles of the laser field in the turn-on stage.
+ - `cyc_num_turn_off`: Number of cycles of the laser field in the turn-off stage.
+ - `cyc_num_const` : Number of cycles of the laser field in the constant-intensity stage.
- `ellip` : Ellipticity of the laser field [-1≤ε≤1, 0 indicates linear polarization and ±1 indicates circular polarization].
- `azi` : Azimuth angle of the laser's polarization's principle axis relative to x axis (in radians) (optional, default 0).
- `t_shift` : Time shift of the laser (in a.u.) relative to the beginning of TURN-ON (optional, default 0).
@@ -44,9 +66,9 @@ struct TrapezoidalLaser <: MonochromaticLaser
- `peak_int` : Peak intensity of the laser field (in W/cm²).
- `wave_len` : Wavelength of the laser field (in nm). Must specify either `wave_len` or `ang_freq`.
- `ang_freq` : Angular frequency of the laser field (in a.u.). Must specify either `wave_len` or `ang_freq`.
- - `cyc_num_turn_on` : Number of cycles of the laser field in the turn-on.
- - `cyc_num_turn_off`: Number of cycles of the laser field in the turn-off.
- - `cyc_num_const` : Number of cycles of the laser field in the constant-intensity.
+ - `cyc_num_turn_on` : Number of cycles of the laser field in the turn-on stage.
+ - `cyc_num_turn_off`: Number of cycles of the laser field in the turn-off stage.
+ - `cyc_num_const` : Number of cycles of the laser field in the constant-intensity stage.
- `ellip` : Ellipticity of the laser field [-1≤ε≤1, 0 indicates linear polarization and ±1 indicates circular polarization].
- `azi` : Azimuth angle of the laser's polarization's principle axis relative to x axis (in radians) (optional, default 0).
- `t_shift` : Time shift of the laser (in a.u.) relative to the beginning of TURN-ON (optional, default 0).
diff --git a/src/SemiclassicalSFI.jl b/src/SemiclassicalSFI.jl
index d234e87..4829c38 100644
--- a/src/SemiclassicalSFI.jl
+++ b/src/SemiclassicalSFI.jl
@@ -32,8 +32,8 @@ Performs a semiclassical simulation with given parameters.
## Required params. for all:
- `init_cond_method = <:ADK|:SFA|:SFAAE|:WFAT|:MOADK>` : Method of electrons' initial conditions. Currently supports `:ADK`, `:SFA`, `:SFAAE` for atoms and `:WFAT`, `:MOADK` for molecules.
-- `laser::Laser` : Parameters of the laser field.
-- `target::Target` : Parameters of the target.
+- `laser::Laser` : A `Lasers.Laser` object containing information of the laser field.
+- `target::Target` : A `Targets.Target` object containing information of the atom/molecule target.
- `sample_t_intv = (start,stop)` : Time interval in which the initial electrons are sampled.
- `sample_t_num` : Number of time samples.
- `traj_t_final` : Time when every trajectory simulation ends.
@@ -53,7 +53,7 @@ Performs a semiclassical simulation with given parameters.
## Optional params. for all:
- `save_path` : Output HDF5 file path.
- `save_3D_spec = false` : Determines whether the 3D momentum spectrum is saved (if not, will only save 2D by flattening on the xy plane) (default `false`).
-- `traj_phase_method = <:CTMC|:QTMC|:SCTS>` : Method of classical trajectories' phase (default `CTMC`). Currently `:QTMC` and `:SCTS` only supports atom targets.
+- `traj_phase_method = <:CTMC|:QTMC|:SCTS>` : Method of classical trajectories' phase (default `:CTMC`). Currently `:QTMC` and `:SCTS` only supports atom targets.
- `traj_rtol = 1e-6` : Relative error tolerance when solving classical trajectories using adaptive methods (default `1e-6`).
- `traj_nondipole = false` : Determines whether the non-dipole effect is taken account in the simulation (default `false`).
- `traj_GPU = false` : [Experimental] Determines whether to enable GPU acceleration in trajectory simulation (default `false`).
diff --git a/src/Targets/HydrogenLikeAtom.jl b/src/Targets/HydrogenLikeAtom.jl
index 19cbfcd..03a9ca5 100644
--- a/src/Targets/HydrogenLikeAtom.jl
+++ b/src/Targets/HydrogenLikeAtom.jl
@@ -1,6 +1,23 @@
using StaticArrays
-"Represents a Hydrogen-like atom."
+"""
+```
+struct HydrogenLikeAtom <: SAEAtomBase
+```
+
+Represents a hydrogen-like atom.
+
+An instance of `HydrogenLikeAtom` can be initialized via the constructor method:
+```julia
+HydrogenLikeAtom(Ip, Z, soft_core=1.0, name="[NA]")
+```
+
+# Example:
+```jldoctest
+julia> t = Targets.HydrogenLikeAtom(Ip=0.5, Z=1.0, soft_core=1.0, name="H")
+[HydrogenLikeAtom] Atom H, Ip=0.5, Z=1.0, SoftCore=1.0
+```
+"""
struct HydrogenLikeAtom <: SAEAtomBase
"Ionization potential of the atom."
Ip;
diff --git a/src/Targets/MolecularCalculators/PySCFMolecularCalculator.jl b/src/Targets/MolecularCalculators/PySCFMolecularCalculator.jl
index feaed45..bb1adf7 100644
--- a/src/Targets/MolecularCalculators/PySCFMolecularCalculator.jl
+++ b/src/Targets/MolecularCalculators/PySCFMolecularCalculator.jl
@@ -10,7 +10,22 @@ using HDF5
using Dates
using ProgressMeter
-"An interface of molecular calculation using PySCF."
+"""
+```
+struct PySCFMolecularCalculator <: MolecularCalculatorBase
+```
+An interface of molecular calculation using PySCF.
+
+An instance of `PySCFMolecularCalculator` can be initialized via the following constructor method:
+
+```
+PySCFMolecularCalculator(; mol, basis::String="pc-1", kwargs...)
+```
+
+# Parameters
+- `mol::Molecule` : The `Molecule` to be calculated.
+- `basis::String` : Basis set used for calculation (default `pc-1`).
+"""
mutable struct PySCFMolecularCalculator <: MolecularCalculatorBase
"The molecule to be calculated."
# mol::Molecule; # linter reported error when specifying type Molecule.
@@ -99,7 +114,21 @@ function DipoleMomentum(mc::PySCFMolecularCalculator)
end
"""
-Calculates the data used in WFAT structure factor calculation of the given molecule.
+```
+calc_WFAT_data(;
+ mc::PySCFMolecularCalculator,
+ orbitIdx_relHOMO::Int = 0,
+ grid_rNum::Int = 200,
+ grid_rMax::Real = 10.,
+ grid_θNum::Int = 60,
+ grid_ϕNum::Int = 60,
+ sf_nξMax ::Int = 5,
+ sf_mMax ::Int = 5,
+ sf_lMax ::Int = 10,
+ kwargs...)
+```
+
+Calculates the DATA used in structure factor calculation in WFAT of the given molecule.
# Returns
`(μ, int_data)` : Orbital dipole momentum and the array which stores the integrals.
@@ -109,11 +138,11 @@ Calculates the data used in WFAT structure factor calculation of the given molec
- `orbitIdx_relHOMO`: Index of selected orbit relative to the HOMO (e.g., 0 indicates HOMO and -1 indicates HOMO-1) (default 0).
- `grid_rNum` : The number of radial grid (default 200).
- `grid_rMax` : The maximum radius of the radial grid (default 10.0).
-- `grid_θNum` : The number of angular grid in the θ direction (default 60).
-- `grid_ϕNum` : The number of angular grid in the ϕ direction (default 60).
-- `sf_nξMax` : The maximum number of nξ used in calculation (default 3).
-- `sf_mMax` : The maximum number of |m| used in calculation (default 3).
-- `sf_lMax` : The maximum angular quantum number l used in calculation (default 6).
+- `grid_θNum` : The number of angular grid in the ``θ`` direction (default 60).
+- `grid_ϕNum` : The number of angular grid in the ``ϕ`` direction (default 60).
+- `sf_nξMax` : The maximum number of ``n_ξ`` used in calculation (default 5).
+- `sf_mMax` : The maximum number of ``|m|`` used in calculation (default 5).
+- `sf_lMax` : The maximum angular quantum number ``l`` used in calculation (default 10).
"""
function calc_WFAT_data(;
mc::PySCFMolecularCalculator,
@@ -122,9 +151,9 @@ function calc_WFAT_data(;
grid_rMax::Real = 10.,
grid_θNum::Int = 60,
grid_ϕNum::Int = 60,
- sf_nξMax ::Int = 5,
- sf_mMax ::Int = 5,
- sf_lMax ::Int = 10,
+ sf_nξMax ::Int = 5,
+ sf_mMax ::Int = 5,
+ sf_lMax ::Int = 10,
kwargs...)
# == PROCEDURE ==
# 0. Obtain the coefficients (finished in the initialization).
@@ -459,6 +488,18 @@ function calc_WFAT_data(;
end
"""
+```
+calc_asymp_coeff(;
+ mc::PySCFMolecularCalculator,
+ orbitIdx_relHOMO::Int = 0,
+ grid_rNum::Int = 200,
+ grid_rReg::Tuple{<:Real,<:Real} = (3,8),
+ grid_θNum::Int = 60,
+ grid_ϕNum::Int = 60,
+ l_max::Int = 6,
+ kwargs...)
+```
+
Calculates the asymptotic coefficients (used in MO-ADK and MO-SFA) of the given molecule.
# Parameters
diff --git a/src/Targets/Molecule.jl b/src/Targets/Molecule.jl
index 39567ec..d2570a6 100644
--- a/src/Targets/Molecule.jl
+++ b/src/Targets/Molecule.jl
@@ -5,7 +5,50 @@ using Rotations
using Dates
using HDF5
-"Represents a molecule."
+"""
+```
+struct Molecule <: Target
+```
+
+Represents a generic molecule.
+
+# Fresh initialization
+
+A new instance of `Molecule` can be initialized via the constructor method:
+```julia
+Molecule(atoms, atom_coords, charge::Integer=0, name::String="[NA]", data_path::String="", calc_energy::Bool=false, rot_α=0.0, rot_β=0.0, rot_γ=0.0)
+```
+
+## Parameters:
+- `atoms` : Atoms in the molecule, stored as a vector of String.
+- `atom_coords` : Atoms' coordinates in the molecule (in **Angstrom**), stored as a N×3 matrix.
+- `charge` : Total charge of the molecule (ion) (optional, default 0).
+- `name` : Name of the molecule (optional).
+- `data_path` : Path to the molecule's data (default empty). Specifying an empty string indicates no saving (but can still be saved later by calling method `MolSaveDataAs`).
+- `calc_energy` : Indicates whether to calculate the energy data of the molecule upon initialization (default false).
+- `rot_α`,`rot_β`,`rot_γ` : Euler angles (ZYZ convention) specifying the molecule's orientation (optional, default 0).
+
+## Example:
+The following example creates a new instance of Hydrogen `Molecule`, and saves the data related to the `Molecule` to the specified path.
+```jldoctest
+julia> mol = Targets.Molecule(atoms=["H","H"], atom_coords=[0 0 -0.375; 0 0 0.375], charge=0, name="Hydrogen", data_path="./Molecule_Hydrogen.h5")
+[ Info: [Molecule] Data saved for molecule Hydrogen at "./Molecule_Hydrogen.h5".
+Molecule [Hydrogen]
+```
+
+# Initialization from existing data
+
+To initialize an instance of `Molecule` from existing external data, invoke
+```julia
+Molecule(data_path::String, rot_α=0.0, rot_β=0.0, rot_γ=0.0)
+```
+
+## Example:
+```jldoctest
+julia> mol = Targets.Molecule("./Molecule_Hydrogen.h5")
+Molecule [Hydrogen]
+```
+"""
mutable struct Molecule <: Target
"Path to the data file that stores information and data about the molecule."
@@ -63,7 +106,7 @@ mutable struct Molecule <: Target
Initializes a new instance of `Molecule` with given parameters.
# Parameters
- `atoms` : Atoms in the molecule, stored as a vector of String.
- - `atom_coords` : Atoms' coordinates in the molecule, stored as a N×3 matrix.
+ - `atom_coords` : Atoms' coordinates in the molecule (in **Angstrom**), stored as a N×3 matrix.
- `charge` : Total charge of the molecule (ion) (optional, default 0).
- `name` : Name of the molecule (optional).
- `data_path` : Path to the molecule's data (default empty). Specifying an empty string indicates no saving (but can still be saved later by calling method `MolSaveDataAs`).
@@ -357,12 +400,28 @@ function MolAsympCoeff_mMax(mol::Molecule, orbitIdx_relHOMO::Integer)
return size(MolAsympCoeff(mol, orbitIdx_relHOMO), 2) - 1
end
-"Gets the Euler angles (ZYZ convention) specifying the molecule's orientation in format (α,β,γ)."
+"""
+```
+MolRotation(mol::Molecule)
+```
+Gets the Euler angles (ZYZ convention) specifying the molecule's orientation in format (α,β,γ).
+"""
MolRotation(mol::Molecule) = (mol.rot_α,mol.rot_β,mol.rot_γ)
-"Sets the Euler angles (ZYZ convention) specifying the molecule's orientation in format (α,β,γ)."
+"""
+```
+SetMolRotation(mol::Molecule, α,β,γ)
+```
+Sets the Euler angles (ZYZ convention) specifying the molecule's orientation in format (α,β,γ).
+"""
function SetMolRotation(mol::Molecule, α,β,γ)
mol.rot_α = α; mol.rot_β = β; mol.rot_γ = γ;
end
+"""
+```
+SetMolRotation(mol::Molecule, (α,β,γ))
+```
+Sets the Euler angles (ZYZ convention) specifying the molecule's orientation in format (α,β,γ).
+"""
function SetMolRotation(mol::Molecule, (α,β,γ))
SetMolRotation(mol, α,β,γ)
end
@@ -423,10 +482,31 @@ function _MolSaveEnergyData(mol::Molecule)
end
"""
-Calculates the WFAT data of the molecule and saves the data.
-- `MCType` : Type of `MolecularCalculator` if it is not initialized. `PySCFMolecularCalculator` if `MC` is not specified.
+```
+MolCalcWFATData!(mol::Molecule,
+ orbitIdx_relHOMO::Integer = 0,
+ MCType::Type = PySCFMolecularCalculator;
+ kwargs...)
+```
+Calculates the WFAT data of the `Molecule` and saves the data.
+- `MCType` : Type of `MolecularCalculator` if the one for this `Molecule` is not initialized before. Default is `PySCFMolecularCalculator` if the `MCType` is not specified.
- `orbitIdx_relHOMO` : Index of selected orbit relative to the HOMO (e.g., 0 indicates HOMO, and -1 indicates HOMO-1) (default 0).
- `kwargs...` : Keyword arguments to pass to the `MolecularCalculator` and the `calcStructFactorData` method, e.g. `basis`, `grid_rNum`, `grid_rMax`, `sf_lMax`, ⋯
+
+# Example:
+```jldoctest
+julia> mol = Targets.Molecule(["H","H"], [0 0 -0.375; 0 0 0.375], 0, "Hydrogen")
+Molecule [Hydrogen]
+
+julia> Targets.MolCalcWFATData!(mol, orbitIdx_relHOMO=0, basis="6-31g")
+[ Info: [PySCFMolecularCalculator] Running molecular calculation...
+[ Info: Finished initialization [taking 0.0466409 second(s)].
+[ Info: [PySCFMolecularCalculator] Running calculation of structure factor data... (ionizing orbital 0 relative to HOMO)
+✓ Calculating the effective potential... (720000 pts) Time: 0:00:22
+Progress: 100%[●●●●●●●●●●●●●●●●●●●●●●●●●] Time: 0:00:22 (31.08 μs/it)
+✓ Calculating the integrals... (7986 integrals) Time: 0:02:48
+Progress: 100%[●●●●●●●●●●●●●●●●●●●●●●●●●] Time: 0:02:48 (21.08 ms/it)
+```
"""
function MolCalcWFATData!(mol::Molecule, orbitIdx_relHOMO::Integer = 0, MCType::Type = PySCFMolecularCalculator; kwargs...)
if isnothing(mol.mol_calc)
@@ -482,10 +562,25 @@ function _MolSaveWFATData(mol::Molecule, orbitIdx_relHOMO::Integer)
@info "[Molecule] WFAT data saved for molecule $(mol.name) at \"$(mol.data_path)\"."
end
"""
-Calculates the asymptotic coefficients of the molecule and saves the data.
-- `MCType` : Type of `MolecularCalculator` if it is not initialized. `PySCFMolecularCalculator` if `MC` is not specified.
+```
+MolCalcAsympCoeff!(mol::Molecule,
+ orbitIdx_relHOMO::Integer = 0,
+ MCType::Type = PySCFMolecularCalculator;
+ kwargs...)
+```
+Calculates the asymptotic coefficients of the `Molecule` and saves the data.
+- `MCType` : Type of `MolecularCalculator` if the one for this `Molecule` is not initialized before. Default is `PySCFMolecularCalculator` if `MCType` is not specified.
- `orbitIdx_relHOMO` : Index of selected orbit relative to the HOMO (e.g., 0 indicates HOMO, and -1 indicates HOMO-1) (default 0).
- `kwargs...` : Keyword arguments to pass to the `MolecularCalculator`, e.g. `grid_rNum`, `l_max`.
+
+## Example:
+```jldoctest
+julia> mol = Molecule(["H","H"], [0 0 -0.375; 0 0 0.375], 0, "Hydrogen")
+Molecule [Hydrogen]
+
+julia> MolCalcAsympCoeff!(mol, basis="6-31g")
+[ Info: [PySCFMolecularCalculator] Running calculation of asymptotic coefficients... (ionizing orbital 0 relative to HOMO)
+```
"""
function MolCalcAsympCoeff!(mol::Molecule, orbitIdx_relHOMO::Integer = 0, MCType::Type = PySCFMolecularCalculator; kwargs...)
if isnothing(mol.mol_calc)
@@ -538,7 +633,12 @@ function _MolSaveAsympCoeff(mol::Molecule, orbitIdx_relHOMO::Integer)
@info "[Molecule] Asymptotic coefficients of orbital index $(orbitIdx_relHOMO) saved for molecule $(mol.name) at \"$(mol.data_path)\"."
end
-"Saves the data of the `Molecule` to the `data_path` (will change the `Molecule`'s inner field `data_path`)."
+"""
+```
+MolSaveDataAs(mol::Molecule, data_path::String)
+```
+Saves the data of the `Molecule` to the `data_path` (will change the `Molecule`'s inner field `data_path`).
+"""
function MolSaveDataAs(mol::Molecule, data_path::String)
function defaultFileName()
Y,M,D = yearmonthday(now())
diff --git a/src/Targets/SAEAtom.jl b/src/Targets/SAEAtom.jl
index c055c70..7ef11a9 100644
--- a/src/Targets/SAEAtom.jl
+++ b/src/Targets/SAEAtom.jl
@@ -1,6 +1,23 @@
using StaticArrays
-"Represents an atom under single-active-electron (SAE) approximation."
+"""
+```
+struct SAEAtom <: SAEAtomBase
+```
+
+Represents an atom under single-active-electron (SAE) approximation.
+
+An instance of `SAEAtom` can be initialized via the constructor method:
+```julia
+SAEAtom(Ip, Z, a1=0.0, b1=0.0, a2=0.0, b2=0.0, a3=0.0, b3=0.0, name="[NA]")
+```
+
+# Example:
+```jldoctest
+julia> t = Targets.SAEAtom(Ip=0.9035698802, Z=1.0, a1=1.230723, b1=0.6620055, a2=-1.325040, b2=1.236224, a3=-0.2307230, b3=0.4804286, name="He")
+[SAEAtom] Atom He, Ip=0.9035698802, Z=1.0
+```
+"""
struct SAEAtom <: SAEAtomBase
# should implement TargetPotential, TargetForce, TrajectoryFunction, ADKRateExp.
"Ionization potential of the atom."
diff --git a/src/Targets/Targets.jl b/src/Targets/Targets.jl
index d856ca2..2d0c27f 100644
--- a/src/Targets/Targets.jl
+++ b/src/Targets/Targets.jl
@@ -1,7 +1,7 @@
"""
-The Targets module provides information about the targeting atoms or molecules.
+The `Targets` module provides information about the targeting atoms or molecules.
"""
module Targets