Fast recombination time using QE

[Shudong Wang]

Dear developers and users,

I reproduce the results of Long's paper discribing the phosphorene recombinationas mentioned in JPCL,7,2016,653-659, where the recombination is 0.108 ns. But muy repeat result

is about 3 ps. I used single unit cell of phosphorene, QE for ab initio MD and relaxation. 

My input file is :

# Simulation type
params["runtype"] = "namd"                   # Type of calculation to perform. Possible values:
                                             # "namd" - to do NA-MD calculations, "no-namd"(or any other) - to
                                             # perform only pre-processing steps - this will create the files with
                                             # the energies of basis states and will output some useful information,
                                             # it may be particularly helpful for preparing your input
params["decoherence"] = 1                    # Do you want to include decoherence via DISH? Possible values:
                                             # 0 - no, 1 - yes
params["is_field"] = 0                       # Do you want to include laser excitation via explicit light-matter
                                             # interaction Hamiltonian? Possible values: 0 - no, 1 - yes

# Integrator parameters
params["elec_dt"] = 0.1                      # Electronic integration time step, fs
params["nucl_dt"] = 1.0                      # Nuclear integration time step, fs (this parameter comes from
                                             # you x.md.in file)
params["integrator"] = 0                     # Integrator to solve TD-SE. Possible values: 0, 10,11, 2

# NA-MD trajectory and SH control
params["namdtime"] = 6000                    # Trajectory time, fs
params["num_sh_traj"] = 1500                 # Number of stochastic realizations for each initial condition
params["boltz_flag"] = 1                     # Boltzmann flag (set to 1 anyways)
params["Temp"] = 300.0                       # Temperature of the system
params["alp_bet"] = 0                        # How to treat spin. Possible values: 0 - alpha and beta spins are not
                                             # coupled to each other, 1 - don't care about spins, only orbitals matter

params["debug_flag"] = 0                     # If you want extra output. Possible values: 0, 1, 2, ...
                                             # as the number increases the amount of the output increases too
                                             # Be carefull - it may result in a huge output!

# Parameters of the field (if it is included)
params["field_dir"] = "xyz"                 # Direction of the field. Possible values: "x","y","z","xy","xz","yz","xyz"
params["field_protocol"] = 1                # Envelope function. Possible values: 1 - step function, 2 - saw-tooth
params["field_Tm"] = 25.0                   # Middle of the time interval during which the field is active
params["field_T"] = 25.0                    # The period (duration) of the field pulse
params["field_freq"] = 3.0                  # The frequency of the field radiation = energy of the photons
params["field_freq_units"] = "eV"           # Units of the above quantity. Possible values: "eV", "nm","1/fs","rad/fs"
params["field_fluence"] = 1.0               # Defines the light radiation intensity (fluence), mJ/cm^2

###############################################################

# Set active space and the basis states
params["active_space"] =  range(1,30) #[175,176,...,185,186]

params["states"] = []
params["states"].append(["GS",[10,-10],0.00])  # ground state
params["states"].append(["S0",[11,-10],0.69])  # excited state 

################################################################

############ Initial conditions choice:222222
nmicrost = len(params["states"])
ic = []
i = 0
while i<100:
    j = 0
    while j<nmicrost:
        ic.append([20*i,j])
        j = j + 1
    i = i + 1

params["iconds"] = ic

and here is the input for NAC

# on the NAC calculation step
python -c "from PYXAID import *
params = { }
params[\"NP\"]=16
params[\"EXE\"]=\"$exe_qespresso\"
params[\"EXE_EXPORT\"]=\"$exe_export\"
params[\"start_indx\"]=$param1
params[\"stop_indx\"]=$param2
params[\"wd\"]=\"wd_test\"
params[\"rd\"]=\"$res\"
params[\"minband\"]=1
params[\"nocc\"]=10
params[\"maxband\"]=64
params[\"nac_method\"]=0
params[\"wfc_preprocess\"]=\"complete\"
params[\"do_complete\"]=1
params[\"prefix0\"]=\"x0.scf\"
params[\"prefix1\"]=\"x1.scf\"
params[\"compute_Hprime\"]=1
print params
runMD.runMD(params)
"

Why the two results differ so much?

Thanks!

Shudong

[Nijamudheen A]

Hi Shudong,

I think, using a simple unit cell (4 P atoms) may not be sufficient to calculate the relaxation timescales. In PYXAID we use gamma point approximation in the NAC calculations. So, a small unit cell would introduce large errors than using a converged large super cell. With this small unit cell, if you have used gamma point calculations for the ground state MD, this would introduce errors in the vibrational properties. During NAMD calculations, there would be significantly large concentration of charge, another source of error. In Long's paper, 4x4x1 super cell of phosphorene is used. Could you please repeat the calculations with a similar supercell? Please, also verify that the pseudopotentials used, various cutoffs are all good.

Regards,

Nijam

[Shudong Wang]

Thank you Nijam,

I used ultrasoft-pseudopotential and the cutoff is 40 Ry. I will take the  same supercell as Long's and let you konw the rsults soon...

By ther way, how we should  take the size of the unitcell? I mean how big the supercell we should take in the calculations.

Thank you agian!

[brendanqhd]

Hi, Shudong,

Based on Nijam's response, the super cell should be 4x4x1 in size. You can easily do this by taking a given unit cell, and extending it to your desired size (4x4x1 in this case), in VESTA, and then using those coordinates for your calculation. I believe the minimum energy cutoff for ultra soft pseudo potentials (PP) is ~ 20 (Ry) , so with 40 (Ry), you may be alright. However for different PP's , such as PAW, etc, the minimum suggested cutoff changes. As always, a convergence study will suggest to you which cutoffs are the most optimal with regards to accuracy / computational cost.

Best,
Brendan

[Nijamudheen A]

   Ideally, size of the supercell should be converged with respect to the gamma point electronic structure calculations, since we use only a gamma point for NAC calculations in PYXAID. I think, PYXAID2 can be used for running NAC calculations with multiple k-points. By using a converged k-point mesh for adiabatic MD simulations, one can fix the error in nuclear dynamics. Finally, during NA-MD simulations, a large supercell may be necessary to mimic the small concentration of charge carriers in experiments. 

[Alexey Akimov]

Dear Shudong,

I'd recommend you to make sure you use exactly the same configuration as Run did. 

As a general comment, I keep seeing that there is a disturbing issue with reproducibility of the dynamics computed with Pyxaid. I don't think this is just a Pyxaid-related thing, more an issue of the precision with which the computational procedures are described. Many things may depend on the choice of the unit cell size, DFT functional, nuclear dynamics protocol (e.g. thermostat or not), perhaps, there is another layer of complexity with the VASP vs QE details of electronic structure implementations. Another factor could be the choice of the active space - e.g. just HOMO-LUMO or also other nearby orbitals. 

All,

As many of you do the Pyxaid-related calculations, I ask you all to provide as much details of your calculations (in the SI or online - e.g. via GitHub) in the next paper you publish. Although the routine-type calculations such as all sort of convergences are not regarded as highly as "physical insights" these days, they are important to ensure consistency of all reported data. If everyone adopts this way, it'll help the entire community.

Best,

Alexey

[Shudong Wang]

Dear Alexey and all,

Firstly, thank you all of your reply for my stupid questions.

I am reproducing the published data just so as to correct my own work, or solve the issue I am concerned. I put the input file about NAMD and other related

files just searching for any suggestions about it, since I am not the expert of NAMD.

I did the calculations the same configuration as the exmple did, as possibly as I can obtain from the code itself, the forum here...

Anyway, thanks again.

Shudong 

[Tong Yang]

Dear all,

I am also trying to reproduce the electron-hole recombination rate in phosphorus monolayer in Long's JPCL work. I have obtained the 6000-step MD trajectory using QE and already obtained the wavefunctions of the first 1000 steps. I averaged the Hamiltonian matrix and the transition dipole moments over the first 1000 steps using py-scr6.py. However, the average of the real part of the Hamiltonian matrix, the average transition dipole moments and the absorption intensity are all zero except the average of the imaginary part of the Hamiltonian matrix. It seems that there must be something wrong that I didn't realize. Because the PBE band gap of monolayer phosphorus is around 0.9 eV, I should observe peaks at.0.9 eV or higher energy in the absorption spectrum. What could be the cause of such unphysical picture?

calculation details:

1. I use a 4x4 supercell

2. Only the wavefunction at the Gamma point is calculated for every step

3. The calculation is non-spin-polarized.

4. There are 320 electrons in the system, namely 160 occupied bands. I include 20 empty bands by setting nbnd=180

5. I use pyxaid to calculate the Hamiltonian matrix, where I include 10 valence bands and 10 empty bands (minband=150, maxband=170)

6. I checked the Hamiltonian matrix of every step. For 0_Ham_xxx_re, only the diagonal elements are nonzero, while for 0_Ham_xxx_im, the diagonal entries vanish. 

7. I use py-scr6.py to average the Hamiltonian matrix, the transition dipole moments and estimate the absorption intensity.

YANG Tong

[Alexey Akimov]

Hi Yang Tong,

This is just a convention used in the Pyxaid - use the imaginary part. Please also see other posts in this forum - I think, it has been discussed several times. 

Best,

Alexey