Density Functional Theory calculations at the University of Texas al El Paso
User login
Recent Pictures
- 1 of 4
- next ›
Monthly archive
- January 2023 (1)
- November 2022 (1)
- April 2022 (3)
- February 2022 (1)
- November 2021 (1)
- October 2021 (1)
- September 2021 (4)
- July 2021 (2)
- May 2021 (3)
- March 2021 (2)
Pages
FAQ
FREQUENTLY ASKED QUESTIONS ABOUT NRLMOL
Tunna Baruah, Mark R. Pederson
- Which molecular properties are calculated in NRLMOL ?
The molecular properties one can get from NRLMOL are :
Electronic structure, density of states, bond lengths, optimized geometry, spin magnetic moment, magnetic anisotropy energy, dipole moment, harmonic vibrational frequencies, infrared and Raman spectra,
electronic density, orbital density, electron localization function, ionization potential, electron affinity, optical absorption spectra, joint density of states, atomic charge, screened and unscreened polarizabilities, vibrationa polarizability etc. - Which exchange-correlation functionals can be used in NRLMOL ?
The LDA exchange-correlation functionals are:
VWN
CA
Perdew-Zunger 81
RPA
Wigner Interpolation Functional
Kohn-Sham exchange only
Gunnarsson-LundquistThe GGA functionals are PW91 and PBE.
- How can I start a spin-polarized calculation ?
If you specify a spin magnetic moment in the CLUSTER file, the calculation will automatically be spin-polarized. For a system with spin S=0, to run a spin-polarized calculation you have to make changes in the SYMBOL file. In the SYMBOL file, the number of ELECTRON with down spin should be specified by a '-' sign. e.g. in methane it should be
ELECTRONS = 5.00000000 -5.00000000
- What is the basis set used in NRLMOL ?
The basis set is specified in the ISYMGEN file which is created when you run NRLMOL. This file contains the basis set for each identity atom. For each orbital of a given atom, same set of primitive Gaussians is used. For example, the default basis set for carbon is written as below :
6 6 ELECTRONIC AND NUCLEAR CHARGE
ALL ALL-ELECTRON ATOM TYPE
3 NUMBER OF ATOMS OF TYPE CAR
ALL-CAR001
ALL-CAR002
ALL-CAR003
EXTRABASIS CONTROLS USAGE OF SUPPLEMENTARY BASIS FUNCTIONS
12 NUMBER OF BARE GAUSSIANS
5 4 3 NUMBER OF S,P,D FUNCTIONS
0 0 1 SUPPLEMENTARY S,P,D FUNCTIONS
0.22213361D+05 0.33317370D+04 0.75790135D+03
0.21454372D+03 0.69924889D+02 0.25086135D+02
0.95910418D+01 0.38024557D+01 0.14891854D+01
0.57487653D+00 0.21494732D+00 0.77209650D-01
0.19792249D+00 0.36998977D+00 0.63644615D+00
0.10124931D+01 0.14480787D+01 0.17173689D+01
0.14931932D+01 0.68987161D+00 0.86072247D-01
-0.16566695D-02 0.37766033D-03 -0.47105343D-04
-0.45005260D-01 -0.84621052D-01 -0.14496564D+00
-0.23535601D+00 -0.34215368D+00 -0.44595124D+00
-0.45263971D+00 -0.32216414D+00 -0.12988420D-01
0.20135471D+00 0.12769913D+00 0.14135467D-01
0.00000000D+00 0.00000000D+00 0.00000000D+00
0.00000000D+00 0.00000000D+00 0.00000000D+00
0.00000000D+00 0.00000000D+00 0.00000000D+00
0.10000000D+01 0.00000000D+00 0.00000000D+00
0.00000000D+00 0.00000000D+00 0.00000000D+00
0.00000000D+00 0.00000000D+00 0.00000000D+00
0.00000000D+00 0.00000000D+00 0.00000000D+00
0.00000000D+00 0.10000000D+01 0.00000000D+00
0.00000000D+00 0.00000000D+00 0.00000000D+00
0.00000000D+00 0.00000000D+00 0.00000000D+00
0.00000000D+00 0.00000000D+00 0.00000000D+00
0.00000000D+00 0.00000000D+00 0.10000000D+01
0.23138630D-01 0.42649133D-01 0.74658851D-01
0.12024115D+00 0.18351176D+00 0.24706804D+00
0.30714219D+00 0.31372706D+00 0.26726340D+00
0.14756585D+00 0.47585576D-01 0.72796459D-02
0.00000000D+00 0.00000000D+00 0.00000000D+00
0.00000000D+00 0.00000000D+00 0.00000000D+00
0.00000000D+00 0.00000000D+00 0.00000000D+00
0.10000000D+01 0.00000000D+00 0.00000000D+00
......
......
......Here, the first line specifies the nuclear and electronic charges in the atom. The nuclear and electronic charges specifies the actual atom (e.g. Nuc charge = 6 for carbon) but the electronic charge depends on
whether all electron or pseudo potential calculations are used (e.g. for carbon it would be 6 for all electron and 4 for a pseudopotential calculation)The second line specifies the type of calculation :
ALL for all-electron BHS for BHS pseudopotential TAB for tabulated user-supplied pseudopotential
The third and fourth lines specify the number of total such atoms in the geometry and their symbols in the SYMBOL file. EXTRABASIS=1 in SYMBOL file will signal the program to use the supplementary basis functions. Then comes the number of primitive Gaussians followed by the number of the s, p and d -type contracted Gaussians. The number of supplementary functions of s, p and d type are written next. These informations are followed by some blocks of numbers. The first block lists the exponents of the primitive Gaussians. This is then followed by Ns blocks where Ns is the number of contracted s-type Gaussians. The first block is the coefficients multiplying the primitive Gaussians for the 1s contracted Gaussian, the second block is for 2s Gaussian and so on. After the Ns number of such blocks, comes the Np blocks corresponding to the p-type contracted Gaussians followed by similar Nd number of d-type Gaussians. In the example for carbon atom above, the 1s and 2s contracted Gaussian is a linear combination of the all the primitive Gaussians whereas the higher unoccupied s orbitals are taken as single long-range Gaussians. Similarly for p orbitals where only the 2p orbital is occupied. These are then followed by similar blocks corresponding to the supplementary functions. The supplementary functions are used only when
EXTRABASIS = 1
in the SYMBOL file. These are generally polarization functions and are usually only used for the calculation of IR and Raman intensities and possibly for dipole moments and polarizabilities.
The carbon atom basis set in the example might be called a "12-1211+++G(3d)" basis in Pople/Gaussian parlance if we consider the last s,p,d orbitals as diffuse functions. This means 12 primitive gaussians contracted to the core 1s orbital, then a triple-zeta valence set (three sets of s and three sets of p orbitals) represented by 12, 1 and 1 primitive gaussians, three sets of d polarization functions and three sets of diffuse functions (s,p,d) each represented by a single primitive gaussian. This is an example of a "generally contracted" basis in which the same set of primitives contributes to all the orbitals. If the diffuse functions are considered as valence, then this basis has a quadruple-zeta functions.
- NRLMOL writes out a bunch of files. Are all of them output files ?
NRLMOL writes out the output in the screen file as well as in particular output files. Here is a list of the output files (you will get almost all information from the screen file as well ):
EVALUES : The eigenvalues, their spin, symmetry representation, degeneracy and occupancy,
Fermi energy for the current iteration.
EVALNNN : The eigenvalues, their spin, symmetry representation, degeneracy and occupancy,
Fermi energy for each iteration number NNN.
SUMMARY : Total energy, electronic charge, kinetic energy and trace of hamiltonian for each
iteration.
MAJNNN : Majority spin density of states for the atom number NNN.
MINNNN : Minority spin density of states for the atom number NNN.
GEOCNVRG : Convergence criteria, Total energy , largest force, information about geometry
optimization. Is written after SCF.
FRCOUT : Total energy, forces on each atom, dipole moment, applied electric field .
XMOL.DAT : Goemetry in the xyz format. The values are in Angstrom.
DIPOLE : Dipole moment.
ATOMSPHNN: Charge and spin charge in each inequivalent atom integrated over a sphere of
specified radius.
SPNRES: Gamma matrices, anisotropy energy, energy eigenvalue with spin-orbit coupling
(more description later).
JNTOUT : Joint density of states and absoprtion spectra.
RHOTOT : Total density on a specified grid. Is written in Gaussian cubic format.
RHOSPN : Spin density on a specified grid. Is written in Gaussian cubic format.
WFHOMONN: Orbital density on a specified grid. Is written in Gaussian cubic format. NN is
00 for HOMO and is positive for states above HOMO. NN is negative for states below HOMO. - In the self-consistency cycle, the energy oscillates. How can I reach self-consistency ?
This happens when the energy levels near the Fermi energy are degenerate and half-filled. The cure is to increase the electronic temperature in TMPTRE from the default value of 0.0001 Hartree. Doing so will smear the occupancies over all the levels and will help in reaching a self-consistent solution.
- The geometry optimization seems to get stuck in the same geometry. How do I get out ?
This seems to happen mostly when the LBFGS scheme is used for geometry optimization. The way out is either delete the .LBF files or use conjugate-gradient method. The later can be started by replacing the LBFGS in the first line of SYMBOL file by CONJUGATE-GRADIENT. The conjugate gradient option seems to be most reliable but often slower than LBFGS by a factor of 2. The user should decide whether to use LBFGS or conjugate-gradient depending on the size of the system and also the smallness of the forces. Using conjugate-gradient towards the end of the optimization seems to be more useful.
- The geometry changes very little in the optimization and forces are small. But the code does not say that convergence is reached. What do I do ?
This is numerical problem and can be solved by refining the mesh used in the calculation. To refine a mesh see following question.
- Is it possible to change the mesh by the user ?
Yes, the parameters used in creating the mesh are listed in the MESHDAT file. The default MESHDAT file contains the following parameters :
F line 1
0.10000E-06 1.2000 line 2
6 line 3
0.20000 0.40000 0.60000 1.0000 1.6000 line 4
4 line 5
2.1000 10.100 18.100 line 6
2 1 3 5 5 7 9 11 19 21 line 7
4 1 3 5 5 7 9 11 19 21 line 8
4 1 3 6 5 7 9 11 19 21 line 9
6 1 3 6 7 7 9 11 19 21 line 10
0.20000E-06 1.2000 line 11
2.0000 8 line 12
2.0000 line 13Explanation: There are two types of meshes used in the calculation. One is the radial mesh use for calculations within atomic spheres and the other is the interstitial mesh.
Line 1:
Line 2: The numbers in this line pertains to the radial mesh. The radial mesh is constructed to give integrations involving various Gaussians from short to long range. The small number is the error allowed in integration by the mesh. The second number 1.2 separates the exponentials of the Gaussians which are tested i.e. the n+1st exponential is 1.2 times nth exponent. Decreasing the error tolerance will result in increase in number of mesh points and better integrals.
Line 3: The atomic sphere is divided into different concentric regions with different numbers of mesh points. The number of such radial zones is written in line 3.
Line 4: The outer radii of the radial zones.
Line 5: The mesh has to be different for different atoms. Hence the periodic table (upto Z=56) is divided in 4 types of meshes.
Line 6: This line shows the 4 types of meshes for atoms : first one for Z<2.1, second one for Z<10.1 and third on for Z<18.1 and the last one for Z>18.
Line 7,8,9,10: The four lines contain the parameters for each type of atomic mesh. The numbers are
NPATS, NPIST, NTHET, NPHI, LMAX FOR EACH RADIAL ZONE
NPATS :
NPIST, NTHET, NPHI : Before creating the mesh, the space is divided into boxes such that each atom is contained in one box. Then a sphere around the atom is assumed and the radial mesh is created. The space at the box corners between the atomic sphere and the rectangular box is divided in a different mesh. The parameters NPIST, NTHET and NPHI correspond to the mesh at the box corners.
LMAX : LMAX for the each radial zone. This mesh will integrate a function of the type
r lmax exp(-α r2) within the given error limit in this region. Increasing LMAX results in a more refined mesh.
Line 11 : This line contains the same parameters as line 2 but for the interstitial mesh. Decreasing the error limit results in more accurate integrals in the interstitial region.
Line 12: Cutfac, mx1d
CUTFAC: If a box transforms into itself due to symmetry, it will be split if it is larger than cutfac times the distance to the closest atom.
MX1D: max. number of points in a one-dimensional interstitial partition
Line 13: Splrat
SPLRAT: largest allowed ratio for: the size of an atomic box divided by the distance of another atom to the box boundary
- What are the symmetry groups recognized by NRLMOL ?
Only a handful of symmetry groups are hard-coded in NRLMOL. These are :
Ih, Th, Oh, Td, X, Y, Z, XY, YZ, XZ, C3v, ...
Caution: C3v,... assumes the highest symmetry rotational axis along the z direction.
It however allows use of any(most) symmetry group provided corresponding GRPMAT file containing symmetry opration matrices is provided. e.g. Point groups D20, D30 can also be used in NRLMOL. - How can I know the symmetry of the energy levels ?
The eigen values are not labeled by its symmetry in NRLMOL. However, it can be easily identified by the representation it belongs to and the character of the representation which is written near the beginning in the screen output.
- I am doing calculations on a singly or doubly charged anion. The energy eigenvalues near the top of the valence band are positive. Does it mean that these electrons are not bound ?
Not necessarily. The binding energy of an electron is given by the difference of the total energy E(n)-E(n+1) where n is the number of electron. The positive eigenvalues is an artifact of the inaccuracy in the exchange-correlation functional in the asymptotic region and can be cured by inclusion of self-interaction correction (SIC) in the calculation. However, the SIC version of NRLMOL is not yet fully operational.
- Can I optimize the geometry in the presence of an electric field ?
To put in an electric field, make the following changes :
I. Change the number of files in SYMBOL to 2 as below :
2 Number of symbolic files
ISYMGEN = INPUT
EFLDVIB = EFIELD
N+3 Number of symbols in list
N Number of nucleiIII. Also make the following changes in SYMBOL:
ELECTRONS = Nup Ndn
EFIELD = x y zwhere Nup and Ndn are the number of spin up and down electrons. x,y,and z are the components of the electric field vector. If optimizing a geometry in the presence of an electric field it may be advisable to pin one of the atoms by setting occurrences of "1 1 1" in the SYMBOL file to "0 0 0". This would be absolutely necessary for a charged molecule or cluster in an electric field. (See item 22)
3. Create a file called EFLDVIB and write EFIELD at the top of the file.
- Can I use the same potential for same nuclear positions but with a different symmetry group ?
Yes, it is possible to use the same potential using different symmetry provided the nuclear positions are fixed. Copy the VMOLD and POTOLD files from the old calculation and change the NCALC in RUNS to the value in the old calculation. Make sure NCALC in RUNS exists in the SYMBOL file.
- Can I fix the spin moment in my calculation ?
To fix the spin magnetic moment, put the following line at the top of the EVALUES file :
FIXM
- Can I specify the occupancies ?
Yes, it is possible to fix the occupancies of the energy levels in the following way :
I. Write the top line of EVALUES as
OCCU
II. Write the number of the spin up states in the first representation and write the occupancies in the following line.
III. Repeat step II for all representation.
IV. Repeat step II and III for spin down states.
Example :
OCCU
9
1.000000 1.000000 1.000000 1.000000 1.000000
1.000000 1.000000 1.000000 1.000000
10
1.000000 1.000000 1.000000 1.000000 1.000000
1.000000 1.000000 1.000000 1.000000 0.0000000E+00
8
1.000000 1.000000 1.000000 1.000000 1.000000
1.000000 1.000000 1.000000
11
1.000000 1.000000 1.000000 1.000000 1.000000
1.000000 1.000000 1.000000 1.000000 1.000000In the above example, there are two representations with 9 and 10 states for spin up and 8 and 11 states for spin down.
Caution: If there is an accidental or near degeneracy for a HOMO and LUMO in the same representation there can be charge oscillations that are akin to those observed when the Fermi Dirac distribution function is used. Sometimes this can be fixed by fractional occupations of the HOMO and LUMO states.
- I have started a calculation with occupied mode for methane. I have given the correct representations and occupancies from the EVALUES file. But my program crashes. Why ?
EVALUES lists only the representations which have at least one states. The correct number of representations can be found in the screen file. If a representation has 0 number of states belonging to it, it still must be specified in the EVALUES file when you run in OCCU mode.
- How do I plot the density ?
Create a file named POTGRID and rerun NRLMOL. IF the calculation is spin polarized then total density is written out in RHOTOT and spin density is written out in RHOSPN. These files have the Gaussian cubic format. Any plotting package that recognises the cubic format can be used to plot these files.
- Can I plot the HOMO density ?
Create a file named WFGRID and rerun NRLMOL. By default the HOMO density is written out in WFHOMO00 file. The last three lines of WFGRID specifies the orbital to be plotted. The states above HOMO are indexed by positive numbers starting with +1 for LUMO and so on. The states below HOMO are indexed by negative numbers starting with -1 for HOMO-1 and so on. To plot LUMO density, open WFGRID and change the index number from 0 to 1. The output file names will be WFHOMO01 for unoccupied higher states and WFHOMO-1 for the occupied lower states. You can write out upto 10 states at a time.
- How can I calculate the vibrational frequencies ?
Later. - Is it possible to fix some atoms in the geometry optimization ?
Yes, you can fix some atoms during the geometry optimization. The SYMBOL file contains the following lines 1.0 1 1 1 corresponding to each atom. The numbers 1(0) 1(0) 1(0) allow the atom to move(remain fixed) along x,y, and z directions according to the forces.
- After convergence, when I looked at the RUNS file, its first line was "0 N+1" where N is the last geometry in the SYMBOL file. Does this mean that, it will first look at N+1st position info in SYMBOl, and if it cannot find it, it will do a new position calculation?
If the geometry convergence criteria is reached in the Nth geometry, the program does not update the nuclear positions in SYMBOL file but the NCALC index in RUNS is changed by +1. If the code is run again, it will look for geometry number N+1 in SYMBOL file and when it does not find it , the program will stop.
- Is it possible to force antiferromagntic spin ordering in some atoms ?
Yes, you can specify the some atoms as spin up and some others as spin down in the SYMBOL file.
Here is an example :ALL-MAN001 = 2.137000 -0.350000 0.003300 SUP
ALL-MAN001 = 3.456900 2.109300 -0.093100 SDN
ALL-OXY003 = -1.271400 0.035600 -0.012700 SUP
ALL-CAR001 = 2.720100 -2.490500 0.106000 UPOThe SUP, SDN and UPO signify spin up , spin down and unpolarized atoms.
The default choice is UPO. - My screen file for the converged geometry got overwritten. Where can I find the spin moment ?
The number of spin up and down electrons are recorded in the VIBINP file.
- Where is the dipole moment written ?
The dipole moment is written in the DIPOLE, VIBINP and the FRCOUT file.
The second line from bottom in FRCOUT and VIBINP gives the x,y , and z components of the dipole vector. - How do I restart my SCF calculation ?
The SCF calculation can be restarted from potential, wavefunction or from an earlier stored Hamiltonian. The user has to make the corresponding choice in the RUNS file. Before restarting a calculation, make sure that the NCALC parameter in RUNS corresponds to the correct geometry in SYMBOL file.
- How do I tell NRLMOL to calculate the molecular properties ?
Some of the molecular properties are calculated automatically once the forces have dropped below the convergence limit. One can still tell the program to calculate properties for any geometry. Most often the user has to create a file with a particular name. Sometimes this file is empty, sometime the user has to write something in it. The following list shows which calculations are triggered by which file.
Property filename initial values
Charges in atomic sphere ATOMSPH
Charge density for plotting RHOGRID
Orbital density for plotting WFGRID
Density of states DOSOCCU
Joint density of states/ DOSJNT
optical properties
Magnetic anisotropy energy SPNORB T F
Electron localization function ELF