FREQUENTLY ASKED QUESTIONS ABOUT NRLMOLΒΆ
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, vibrational 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-Lundquist
The 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 ?
See the discussion of files at
In the self-consistency cycle, the energy oscillates. How can I reach self-consistency ?
This can for a number of reasons. One of them is that your structure is not correct. Visualize it to ensure that your geometry (and units) are correct. It can also happen 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 ?
Please refer to the discussion of MESHDAT file.
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 :
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 nuclei
Also make the following changes in SYMBOL:
ELECTRONS = Nup Ndn EFIELD = x y z where 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)
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 :
Write the top line of EVALUES as
OCCU
Write the number of the spin up states in the first representation and write the occupancies in the following line.
Repeat step II for all representation.
Repeat step II and III for spin down states.
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 UPO The 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
About this document ... next up previous Back to HomePage. Next: About this document ... Rajendra Zope 2007-11-06