MatSQ provides a tutorial video for users who are not familiar with the service.
How would you simulate materials used in your research with DFT? DFT can compute electronic structures for only a few to dozens of models because of the limitations of computing performance. Periodic boundary conditions were introduced as a way of overcoming these limitations and simulating a bulk condition similar to the macroscopic system. In periodic boundary conditions, a material is defined as follows:
In periodic boundary conditions, a structure is repeated infinitely, so a silicon unit cell with eight atoms, a silicon supercell with 32 atoms, and a silicon bulk are theoretically identical. Moreover, because atoms are repeated by the space-group rules within the unit cell, silicon unit cells can be modeled only with the following information:
The computer recognizes a calculation model with a view that the structure shown on Structure Builder is repeated infinitely. If the cell size increases because of increased vacuum, the vacuum is recognized as repetitive, and the calculation is performed in the vacuum space too, so the amount of computation also increases. Roughly speaking, computational cost is proportional to V N (2 < N < 3), as the number of basis functions (plane waves) proportional to the volume of periodic unit cell. If you add a vacuum to model an atom, molecule, or slab structure, it is appropriate to have 10 to 15 Å clearance with the next repetitive model. You can check the repetitive structure by ticking "Ghost" in the right-click menu in the Visualizer Canvas of Structure Builder.
You can obtain the information or structure file required for unit cell modeling at the following link:
The following section gives a brief introduction to Quantum Espresso (QE), one of the DFT codes available on Materials Squire to run materials simulations. To understand in detail how Quantum Espresso works and what it can do, we recommend reading the documentation provided at www.quantum-espresso.org.
Quantum Espresso is an integrated suite of open-source computer codes for electronic-structure calculations and materials modeling at the nanoscale. It is based on the density-functional theory (DFT), plane waves, and pseudopotentials.
- The Quantum Espresso package is an expandable distribution of related packages. Two core packages for DFT electronic structure calculations, PWscf (Plane Wave self-consistent field), and CP (Car-Parrinello Molecular Dynamics) are supplemented by various packages for specialized applications as well as plug-ins. For more information on all packages, refer to the Quantum Espresso official user guide .
- Quantum Espresso provides open-source software packages. This means that everyone can study, expand, and modify the source code, allowing for continuous improvement.
- The theory behind Quantum Espresso’s algorithm for electronic structure calculations and materials modeling is determined by the Kohn-Sham density functional theory (KS-DFT). In other words, the code implements the common iterative self-consistent method to solve the Kohn-Sham equations.
- To perform materials simulations by a computer, we must expand functions like the wave function or the electron density in a basis set. Quantum Espresso uses plane waves as a basis set under the Bloch’s theorem. For the Gamma point, it is a Fourier series expansion of functions. A finite number of expansion coefficients, which is required for computation, can be achieved by an energy cutoff (ecutwfc, ecutrho).
- Pseudopotentials are used to smoothen the Coulomb potential by atomic nucleis, resulting in fewer plane waves, without affecting the result too much. Quantum Espresso allows the use of various pseudopotentials (norm-conserving, ultrasoft, PAW). Choosing the pseudopotential for a given system is a science in itself. Refer to the Quantum Espresso documentation for further information.
Based on DFT, Quantum Espresso has numerous applications, ranging from ground-state energy calculations and structural optimization to molecular dynamics as well as the modeling of response and spectroscopic properties. DFT simulations can be done on any crystal structure or supercell, which features some form of periodicity. Furthermore, Quantum Espresso works for insulators, semiconductors, and metals, offering different options for k-point sampling as well as smearing of energy states. To speed up computations, Quantum Espresso can deal with various pseudopotentials and approximate exchange-correlation functionals. For more information, visit the Quantum Espresso website .
The Quantum Espresso input script contains all information about the system of interest and defines the calculation process. The information consists of the namelist and input_cards. The following figure shows the general syntax of the input script.
There are three mandatory namelists in the PWscf package:
- It lists input variables that control the calculation process or determine the number of I/O. Examples include the calculation type, information amount (verbosity), and directory.
- It includes input variables that determine the calculation system such as the number of atoms, Bravais lattice index, cutoff energies, and smearing methods.
- It controls the algorithms used to reach the self-consistent solution of Kohn-Sham equations. Examples include the convergence threshold for self-consistency and mixing beta.
- If the nuclei of the system are allowed to move in a calculation, this namelist contains necessary variables to control their motion. Variable atomic positions occur in molecular dynamics or structural relaxation computations.
- Similar to &IONS, this namelist must be included for calculations with variable cell dimensions.
For some input data, such as atomic coordinates, it is inconvenient to write it in the namelist syntax. To make life easier, Quantum Espresso, therefore, features input_cards that allow you to enter data in a more practical format. There are three mandatory input_cards to be entered:
- Lists the name, mass, and PP of the atomic species included in a calculation model
- Lists the name and coordinates of all atoms in a calculation model
- Refers to the k-point grid and shift information, which are used to determine the number of k-points to be sampled in each lattice direction
The following links might be helpful to learn more about Quantum Espresso:
| Parameter | Value | Description | |
|---|---|---|---|
| Calculation type | scf | It performs self-consistent field calculations without affecting the atoms’ position. Namelists &IONS and &CELL are ignored in calculations. An iterative solution process runs to calculate the total energy, forces, and stresses. | |
| relax | In a relax calculation, the atoms are allowed to move to find their minimum energy (structural optimization). Includes geometric optimization steps and iterative self-consistent field calculations. | ||
| vc-relax | It optimizes the structure for the atom position and the cell. The cell shape (angles, lattice constants) may change to find the optimized structure. It also includes geometric optimization steps and iterative self-consistent field calculations. | ||
| (vc-)md | It calculates molecular dynamics with DFT. (ab-initio MD, AIMD) | ||
| restart | nscf | It performs non-scf calculations. Using this scheme, you make a single step calculation with the superposition of atomic orbitals. In contrast with the scf calculation, unoccupied electron states are also considered. Therefore, an nscf calculation is an economical choice for calculations that require a large k-point sampling. | |
| bands | It calculates only the Kohn-Sham states for the given set of k-points. | ||
| Max scf steps | It determines the maximum number of scf algorithm runs until convergence is reached (scf is fixed to 1 unconditionally). | ||
| Information amount | low | Default | |
| high | It adds detailed information about k-points or the character table to a job.stdout file. | ||
| Force threshold | It is the convergence threshold for the force to an ionic minimization. Any force for all elements must be less than this value (3.8E-4 Ry/Bohr = 0.01 eV/Å ). | ||
| Time step | It sets the time step for an MD simulation (atomic unit, 1 a.u. = 4.8378*10^-17 s). | ||
| Parameter | Value | Description |
|---|---|---|
| occupations | smearing | It performs Gaussian smearing of occupation numbers on the assumption that an electron would occupy up to a point slightly above the balance band (suitable for metals). |
| fixed | It calculates without smearing on the assumption that the structure is an insulator. | |
| tetrahedra | The tetrahedron method (P.E. Bloechl, PRB 49, 16223 (1994)) is used, which is suitable for DOS calculations. It must use a Monkhorst’s pack (automatic) k-point. | |
| Ecut(wfc) | It is the kinetic energy cutoffs for wave functions. | |
| Ecut(rho) | It is the kinetic energy cutoffs for charge densities and potentials. The default value is "ecutrho = 4*ecutwfc". Norm-conserving potentials and ultrasoft pseudopotentials require a higher ecut (rho) (about 8 to 12 times ecutwfc). | |
| Gaussian broadeing | It is the value required for Gaussian spreading in the Brillouin zone (broadening similar to '0.002 Ry = 27E-3eV = approx.300 K'). | |
| Number of electron spin | 1 (all up) | It calculates without considering the spin. |
| 2 (up, down) | It calculates taking into account the spin polarization. The amount of calculation doubles the case of "nspin=1." | |
| 4 (Noncolinear) | It is used in noncollinear cases. | |
| Van der Walls correction | It is used to compensate data for models that are significantly affected by the Van der Waals force, such as a layered structure. | |
| grimme-d2 (DFT-D) | It corrects by the semi-empirical Grimme’s DFT-D2 method (S. Grimme, J. Comp. Chem. 27, 1787 (2006), V. Barone et al., J. Comp. Chem. 30, 934 (2009)). | |
| grimme-d3 (DFT-D3) | It corrects by the semi-empirical Grimme’s DFT-D3 method (S. Grimme et al., J. Chem. Phys 132, 154104 (2010)). | |
| tkatchenko-scheffler | It corrects the Tkatchenko-Scheffler dispersion with the C6 coefficient induced by the first principles (A. Tkatchenko and M. Scheffler, PRL 102, 073005 (2009)). | |
| XDM | It performs corrections with the exchange-hole dipole-moment model (A. D. Becke et al., J. Chem. Phys. 127, 154108 (2007), A. Otero de la Roza et al., J. Chem. Phys. 136, 174109 (2012)). | |
| Hubbard_U | the U parameter for the element of interest. | |
| Parameter | Value | Description |
|---|---|---|
| Max iteration step | Within an scf step, it sets the maximum step value for iteration, which continues until the convergence is completed. It is recommended to increase this value for structures with poor convergence. | |
| Mixing beta | The rate at which the final electron density is mixed with the initial electron density under the scf algorithm. It is recommended to decrease this value for structures with poor convergence. | |
| Convergence threshold | The value that sets the convergence threshold, which is the limit of the energy difference before and after an scf step. | |
| Mixing mode | plain | The Broyden charge density mixing method |
| TF | It adds a simple Thomas-Fermi screening (applicable to highly consistent systems). | |
| local-TF | It screens TF depending on the local density (applicable to highly consistent systems). | |
| Starting wavefunction | atomic | It starts wave function calculations from the superposition of atomic orbitals. Calculations are performed normally in most cases, but some fail occasionally. |
| atomic+random | In addition to the superposition of atomic orbitals, it considers random wave functions. | |
| random | It starts a calculation with random wave functions. It has a slower but safer start of scf. | |
| Parameter | Value | Description | |
|---|---|---|---|
| Ion dynamics | It specifies the algorithm in Structural Relaxation to consider atomic movement. | ||
| relax | bfgs | It uses the BFGS quasi-newton algorithm for structural relaxation; cell_dynamics must be "bfgs" too. (Default). | |
| damp | It uses damped (quick-min Verlet) dynamics for structural relaxation. | ||
| vc-relax | bfgs | It uses BFGS quasi-newton algorithm; cell_dynamics must be "bfgs" too. | |
| damp | It uses damped (Beeman) dynamics for structural relaxation | ||
| md | verlet | It uses the Verlet algorithm to integrate Newton’s equation (Default). | |
| langevin | The ion dynamics is the overdamped Langevin. | ||
| langevin-smc | The overdamped Langevin with Smart Monte Carlo (R.J. Rossky, JCP, 69, 4628 (1978)). | ||
| vc-md | beeman | It uses Beeman’s algorithm to integrate Newton’s equations (Default). | |
| upscale | It reduces conv_thr by conv_thr/upscale during structural optimization to increase accuracy when the relaxation approaches convergence (available in the bfgs option only). | ||
| Ion temperature | (vc-)md | not_controlled | It does not control ionic temperatures (Default). |
| rescaling | It controls ionic temperatures via velocity rescaling (first method). | ||
| md | rescale-v | It controls ionic temperatures via velocity rescaling (second method). | |
| rescale-T | It controls ionic temperatures via velocity rescaling (third method). | ||
| reduce-T | It reduces ionic temperatures for every nraise step by the ΔT value. | ||
| berendsen | It controls ionic temperatures with "soft" velocity rescaling. | ||
| andersen | It controls ionic temperatures with the Andersen thermostat method. | ||
| initial | It initializes ion temperatures to the starting temperature and leaves uncontrolled further on. | ||
| Starting temperature (K) | The starting temperature in MD simulations. | ||
| ΔT | Ion temperature = "rescale-T": At each step, the instantaneous temperature is multiplied by ΔT.
Ion temperature = "reduce-T": In every "nraise" step, the instantaneous temperature is reduced by ΔT. |
||
| nraise | It rescales the instantaneous temperature (https://www.quantum-espresso.org/Doc/INPUT_PW.html#idm938). | ||
| Parameter | Value | Description |
|---|---|---|
| Cell dynamics | It specifies the type of algorithms for variable cell relaxations to control the cell size. | |
| bfgs | The BFGS quasi-newton algorithm; ion_dynamics must be "bfgs" too. | |
| none | no dynamics | |
| sd | It is the steepest descent (not implemented). | |
| damp-pr | It is the damped (Beeman) dynamics of the Parrinello-Rahman extended lagrangian. | |
| damp-w | It is the damped (Beeman) dynamics of the new Wentzcovitch extended lagrangian. | |
| Cell factor | It is the maximum strain ratio of the cell size. | |
| Press threshold | It is the convergence threshold of the pressure applied to cells (Kbar). | |
| Parameter | Value | Description | |
|---|---|---|---|
| Sampling | automatic | It is an option to sample k-points in the Monkhorst-Pack method. It also distributes the k-point grid evenly to the supercell. | |
| GAMMA | It is similar to automatic 1x1x1 in that only one k-point is sampled, but there is a difference: the k-point is recognized as real rather than a complex number. It has the benefit of fast calculation. | ||
| crystal(_b) | It designates the k-point as the relative coordinates to the reciprocal lattice vector. If the last column of data is {crystal}, it represents the weight for each k-point. If the last column of data is {crystal_b}, it represents the k-point number to be sampled by the next crystal coordination. | ||
| # k-points | It is the number of k-points in the direction of the three lattice vectors, respectively. | ||
| shift | It shifts the k-point grid with respect to the origin. Depending on the symmetry of the supercell, shifting the k-points could lead to better results. | ||
| crystal(_b) | Path | You can use a high-symmetric point to set the path for k-point sampling. You have to enter the weight. | |
| weight | crystal: the weight for each k-point to have
crystal_b: the number of k-points to be sampled until the next k-point |
||
This is the entire list of functions available in Quantum Espresso. Put 'Short name' in the DFT Functional input field.
| Short name | Complete name | Short description | References |
|---|---|---|---|
| pz | sla+pz | Perdew-Zunger LDA | J.P.Perdew and A.Zunger, PRB 23, 5048 (1981) |
| bp | b88+p86 | Becke-Perdew grad.corr. | |
| pw91 | sla+pw+ggx+ggc | PW91 (aka GGA) | J.P.Perdew and Y. Wang, PRB 46, 6671 (1992) |
| blyp | sla+b88+lyp+blyp | BLYP | C.Lee, W.Yang, R.G.Parr, PRB 37, 785 (1988) |
| pbe | sla+pw+pbx+pbc | PBE | J.P.Perdew, K.Burke, M.Ernzerhof, PRL 77, 3865 (1996) |
| revpbe | sla+pw+rpb+pbc | revPBE (Zhang-Yang) | Zhang and Yang, PRL 80, 890 (1998) |
| pw86pbe | sla+pw+pw86+pbc | PW86 exchange + PBE correlation | |
| b86bpbe | sla+pw+b86b+pbc | B86b exchange + PBE correlation | |
| pbesol | sla+pw+psx+psc | PBEsol | J.P. Perdew et al., PRL 100, 136406 (2008) |
| q2d | sla+pw+q2dx+q2dc | PBEQ2D | L. Chiodo et al., PRL 108, 126402 (2012) |
| hcth | nox+noc+hcth+hcth | HCTH/120 | Handy et al, JCP 109, 6264 (1998) |
| olyp | nox+lyp+optx+blyp | OLYP | Handy et al, JCP 116, 5411 (2002) |
| wc | sla+pw+wcx+pbc | Wu-Cohen | Z. Wu and R. E. Cohen, PRB 73, 235116 (2006) |
| sogga | sla+pw+sox+pbec | SOGGA | Y. Zhao and D. G. Truhlar, JCP 128, 184109 (2008) |
| optbk88 | sla+pw+obk8+p86 | optB88 | |
| optb86b | sla+pw+ob86+p86 | optB86 | |
| ev93 | sla+pw+evx+nogc | Engel-Vosko | Engel-Vosko, Phys. Rev. B 47, 13164 (1993) |
| tpss | sla+pw+tpss+tpss | TPSS Meta-GGA | J.Tao, J.P.Perdew, V.N.Staroverov, G.E. Scuseria, PRL 91, 146401 (2003) |
| m06l | nox+noc+m6lx+m6lc | M06L Meta-GGA | Y. Zhao and D. G. Truhlar, JCP 125, 194101 (2006) |
| tb09 | sla+pw+tb09+tb09 | TB09 Meta-GGA | F Tran and P Blaha, Phys.Rev.Lett. 102, 226401 (2009) |
| pbe0 | pb0x+pw+pb0x+pbc | PBE0 | J.P.Perdew, M. Ernzerhof, K.Burke, JCP 105, 9982 (1996) |
| b86bx | pb0x+pw+b86x+pbc | B86bPBE hybrid | |
| bhahlyp | pb0x+pw+b88x+blyp | Becke half-and-half LYP | |
| hse | sla+pw+hse+pbc | Heyd-Scuseria-Ernzerhof (HSE 06, see note below) | Heyd, Scuseria, Ernzerhof, J. Chem. Phys. 118, 8207 (2003); Heyd, Scuseria, Ernzerhof, J. Chem. Phys. 124, 219906 (2006). |
| b3lyp | b3lp+b3lp+b3lp+b3lp | B3LYP | P.J. Stephens,F.J. Devlin,C.F. Chabalowski,M.J. Frisch, J.Phys.Chem 98, 11623 (1994) |
| b3lypv1r | b3lp+b3lpv1r+b3lp+b3lp | B3LYP-VWN1-RPA | |
| x3lyp | x3lp+x3lp+x3lp+x3lp | X3LYP | X. Xu, W.A Goddard III, PNAS 101, 2673 (2004) |
| vwn-rpa | sla+vwn-rpa | VWN LDA using vwn1-rpa parametriz | |
| gaupbe | sla+pw+gaup+pbc | Gau-PBE (also "gaup") | |
| vdw-df | sla+pw+rpb +vdw1 | vdW-DF1 | M. Dion et al., PRL 92, 246401 (2004); T. Thonhauser et al., PRL 115, 136402 (2015) |
| vdw-df2 | sla+pw+rw86+vdw2 | vdW-DF2 | Lee et al., Phys. Rev. B 82, 081101 (2010) |
| vdw-df-c09 | sla+pw+c09x+vdw1 | vdW-DF-C09 | |
| vdw-df2-c09 | sla+pw+c09x+vdw2 | vdW-DF2-C09 | |
| vdw-df-obk8 | sla+pw+obk8+vdw1 | vdW-DF-obk8 (optB88-vdW) | Klimes et al, J. Phys. Cond. Matter, 22, 022201 (2010) |
| vdw-df-ob86 | sla+pw+ob86+vdw1 | vdW-DF-ob86 (optB86b-vdW) | Klimes et al, Phys. Rev. B, 83, 195131 (2011) |
| vdw-df2-b86r | sla+pw+b86r+vdw2 | vdW-DF2-B86R (rev-vdw-df2) | |
| vdw-df-cx | sla+pw+cx13+vdW1 | vdW-DF-cx | K. Berland and P. Hyldgaard, PRB 89, 035412 (2014) |
| vdw-df-cx0 | sla+pw+cx13+vdW1+HF/4 | vdW-DF-cx-0 | K. Berland, Y. Jiao, J.-H. Lee, T. Rangel, J. B. Neaton and P. Hyldgaard, J. Chem. Phys. 146, 234106 (2017) |
| vdw-df2-0 | sla+pw+rw86+vdw2+HF/4 | vdW-DF2-0 | |
| vdw-df2-br0 | sla+pw+b86r+vdW2+HF/4 | vdW-DF2-b86r-0 | |
| vdw-df-c090 | sla+pw+c09x+vdw1+HF/4 | vdW-DF-C09-0 | |
| vdw-df-x | sla+pw+????+vdwx | vdW-DF-x, reserved Thonhauser, not implemented | |
| vdw-df-y | sla+pw+????+vdwy | vdW-DF-y, reserved Thonhauser, not implemented | |
| vdw-df-z | sla+pw+????+vdwz | vdW-DF-z, reserved Thonhauser, not implemented | |
| rvv10 | sla+pw+rw86+pbc+vv10 | rVV10 | R. Sabatini et al. Phys. Rev. B 87, 041108(R) (2013) |
This is a brief introduction to LAMMPS, the molecular dynamics code available in Material Square. For more information on how LAMMPS works and possible applications, refer the official manual at lammps.sandia.gov..
LAMMPS (Large-scale Atomic/Molecular Massively Parallel Simulator) is a classical molecular dynamics simulation code, developed at the Sandia National Laboratory, it designed to be calculated effectively with parallel computing. Base on Newton's equation of motion, LAMMPS can calculate a system consists of several hundred to several billion atoms using forcefield for fast speed. It can calculate not only the stability at the given temperature which is a basic property in the materials study, but also mechanical, thermal and chemical properties of material. Recently, the effectivity of parallel computing using the Graphics Processing Unit (GPU) was increased, and machine learning can be used for the expanding interatomic potential.
| Parameter | Value | Description |
|---|---|---|
| Reactive Forcefield | ID | The identification value of each reactive forcefield |
| Type | The type of reactive forcefield | |
| Elements | Elements for the reactive forcefield to consider | |
| Author | The author of the reactive forcefield | |
| DOI | The paper on the corresponding reactive forcefield | |
| Ensemble | NVT | An NVT ensemble is also known as the canonical ensemble. It assumes a calculation model as an isolated system, which has a fixed temperature, volume, and the number of atoms. |
| NVE | An NVE ensemble is also known as the microcanonical ensemble. This assumes a calculation model as an isolated system, which cannot exchange energy or particles with its environment because the energy of the system is fixed. | |
| After Relax | Before starting an MD simulation, a structural relaxation is performed at room temperature (200 K – 300 K) first for structural stabilization. It can reduce the probability that the calculation would fail because of structural instability. | |
| Dump | It saves the trajectory in a file for every step selected. For calculations with long simulation times, you can reduce the overall file size. | |
| Temperature | Begin (K) | The temperature at the start of a simulation |
| Final (K) | The target temperature during simulation | |
| Damping (Step) | Resets the temperature for each damping step | |
| Time (ps) | Specifies the total simulation time (1 ns = 1000 ps) | |
| Initial Velocity (Å/fs) | The initial speed of the atom group selected | |
| Force (Kcal/mole-Å) | The force per fs applied to the selected atom group | |
| Move (Linear, Å/fs) | The distance traveled per fs for the selected atom group | |
To cite the Materials Square, add citation by the URL.