I am trying to do a simple ground state energy calculation of methane molecule using the Octopus (TDDFT) code. The example is already available in the octopus manual. Here is the link. As per the manual, the total converged energy for CH4 molecule is -219.01537542 eV. When I am repeating the calculation with the same input file, the output is -138.66600447 eV. The different versions of the code may be one of reasons (I am using octopus-10.0), but the difference between the energies is very high. I am attaching the log file and the output file for your reference.
**LOG FILE**
************************** Calculation Mode **************************
Input: [CalculationMode = gs]
**********************************************************************
Input: [PseudopotentialSet = standard]
Reading Coordinates from Coordinates block
****************************** Species *******************************
Species 'C'
type : pseudopotential
file : '/usr/local/share/octopus/pseudopotentials/PSF/C.psf'
file format : PSF
valence charge : 4.0
atomic number : 6
form on file : semilocal
orbital origin : calculated
lmax : 1
llocal : 0
projectors per l : 1
total projectors : 1
application form : kleinman-bylander
orbitals : 16
bound orbitals : 4
Species 'H'
type : pseudopotential
file : '/usr/local/share/octopus/pseudopotentials/PSF/H.psf'
file format : PSF
valence charge : 1.0
atomic number : 1
form on file : semilocal
orbital origin : calculated
lmax : 0
llocal : 0
projectors per l : 1
total projectors : 0
application form : local
orbitals : 16
bound orbitals : 1
***************************** Symmetries *****************************
Symmetry elements : 4*(C3) 3*(C2) 3*(S4) 6*(sigma)
Symmetry group : Td
**********************************************************************
Input: [SpinComponents = unpolarized]
Input: [SmearingFunction = semiconducting]
Input: [SymmetrizeDensity = no]
******************************* States *******************************
Total electronic charge = 8.000
Number of states = 4
States block-size = 4
**********************************************************************
Input: [CurvMethod = curv_uniform]
Input: [DerivativesStencil = stencil_star]
************************** Parallelization ***************************
Info: Octopus will run in *serial*
**********************************************************************
Info: Generating weights for finite-difference discretization of x-gradient
Info: Generating weights for finite-difference discretization of y-gradient
Info: Generating weights for finite-difference discretization of z-gradient
Info: Generating weights for finite-difference discretization of Laplacian
******************************** Grid ********************************
Simulation Box:
Type = minimum
Radius [A] = 3.500
Octopus will run in 3 dimension(s).
Octopus will treat the system as periodic in 0 dimension(s).
Main mesh:
Spacing [A] = ( 0.220, 0.220, 0.220) volume/point [A^3] = 0.01065
# inner mesh = 31525
# total mesh = 51605
Grid Cutoff [eV] = 776.921749 Grid Cutoff [Ry] = 57.102700
**********************************************************************
Info: states-block size = 1.6 MiB
Input: [StatesPack = yes]
Input: [StatesOrthogonalization = cholesky_serial]
Info: the XCFunctional has been selected to match the pseudopotentials
used in the calculation.
**************************** Theory Level ****************************
Input: [TheoryLevel = dft]
Exchange-correlation:
Exchange
Slater exchange (LDA)
[1] P. A. M. Dirac, Math. Proc. Cambridge Philos. Soc. 26, 376 (1930)
Correlation
Perdew & Zunger (Modified) (LDA)
Input: [SICCorrection = sic_none]
**********************************************************************
****************************** Hartree *******************************
Input: [DressedOrbitals = no]
The chosen Poisson solver is 'interpolating scaling functions'
**********************************************************************
Input: [FilterPotentials = filter_TS]
Info: Pseudopotential for C
Radii for localized parts:
local part = 1.4 A
non-local part = 1.0 A
orbitals = 4.8 A
Info: Pseudopotential for H
Radii for localized parts:
local part = 1.0 A
non-local part = 0.0 A
orbitals = 4.6 A
Input: [RelativisticCorrection = non_relativistic]
Input: [DFTULevel = dft_u_none]
****************** Approximate memory requirements *******************
Mesh
global : 1.0 MiB
local : 1.2 MiB
total : 2.2 MiB
States
real : 1.6 MiB (par_kpoints + par_states + par_domains)
complex : 3.1 MiB (par_kpoints + par_states + par_domains)
**********************************************************************
Info: Generating external potential
done.
Info: Octopus initialization completed.
Info: Starting calculation mode.
Info: Allocating ground state wave-functions
Info: Blocks of states
Block 1 contains 4 states: 1 - 4
Info: Ground-state allocation done.
** Warning:
** Could not find 'restart/gs' directory for restart.
** No restart information will be read.
** Warning:
** Unable to read wavefunctions.
** Starting from scratch!
Input: [MixField = potential] (what to mix during SCF cycles)
Input: [MixingScheme = broyden]
**************************** Eigensolver *****************************
Input: [Eigensolver = cg]
Input: [Preconditioner = pre_filter]
Input: [PreconditionerFilterFactor = 0.5000]
Input: [SubspaceDiagonalization = standard]
**********************************************************************
Input: [EigensolverSkipKpoints = no]
Input: [LCAOStart = lcao_states]
Input: [LCAOScaleFactor = 1.000]
Input: [LCAOMaximumOrbitalRadius = 10.58 A]
Info: Single-precision storage for 4 extra orbitals will be allocated.
Info: Unnormalized total charge = 7.997443
Info: Renormalized total charge = 8.000000
Info: Setting up Hamiltonian.
Info: Performing initial LCAO calculation with 8 orbitals.
Info: Getting Hamiltonian matrix elements.
ETA: .......1......2.......3......4......5.......6......7.......8......9......0
Eigenvalues [eV]
#st Spin Eigenvalue Occupation
1 -- ************ 2.000000
2 -- ************ 2.000000
3 -- ************ 2.000000
4 -- ************ 2.000000
Info: Ground-state restart information will be written to 'restart/gs'.
Info: SCF using real wavefunctions.
Info: Starting SCF iteration.
ETA: .......1......2.......3......4......5.......6......7.......8......9......0
*********************** SCF CYCLE ITER # 1 ************************
etot = -1.52078533E+02 abs_ev = 1.19E+08 rel_ev = 2.66E+06
ediff = -5.59E+00 abs_dens = 3.63E+00 rel_dens = 4.54E-01
Matrix vector products: 8
Converged eigenvectors: 0
# State Eigenvalue [eV] Occupation Error
1 -13.785511 2.000000 (3.4E-01)
2 -8.680159 2.000000 (4.1E-01)
3 -0.798552 2.000000 (3.8E-01)
4 0.777761 2.000000 (3.3E-01)
Density of states:
----------------------------------------------------------------------
----------------------------------------------------------------------
Elapsed time for SCF step 1: 0.09
**********************************************************************
ETA: .......1......2.......3......4......5.......6......7.......8......9......0
*********************** SCF CYCLE ITER # 2 ************************
etot = -1.48054766E+02 abs_ev = 4.02E+00 rel_ev = 9.83E-02
ediff = 1.48E-01 abs_dens = 1.09E-14 rel_dens = 1.36E-15
Matrix vector products: 8
Converged eigenvectors: 0
# State Eigenvalue [eV] Occupation Error
1 -12.755744 2.000000 (3.3E-01)
2 -8.213151 2.000000 (4.2E-01)
3 -0.790464 2.000000 (3.7E-01)
4 1.284781 2.000000 (3.4E-01)
Density of states:
----------------------------------------------------------------------
----------------------------------------------------------------------
ETA: .......1......2.......3......4......5.......6......7.......8......9......0
*********************** SCF CYCLE ITER # 3 ************************
etot = -1.38666004E+02 abs_ev = 9.39E+00 rel_ev = 2.97E-01
ediff = 3.45E-01 abs_dens = 1.88E-14 rel_dens = 2.35E-15
Matrix vector products: 8
Converged eigenvectors: 0
# State Eigenvalue [eV] Occupation Error
1 -10.617759 2.000000 (3.1E-01)
2 -7.350349 2.000000 (4.5E-01)
3 -0.625441 2.000000 (3.8E-01)
4 2.813351 2.000000 (3.8E-01)
Density of states:
----------------------------------------------------------------------
----------------------------------------------------------------------
Elapsed time for SCF step 3: 0.12
*********************************************************************
Info: SCF converged in 3 iterations
** Warning:
** Some of the states are not fully converged!
Info: Finished writing information to 'restart/gs'.
Octopus emitted 3 warnings.
**OUTPUT**
******************************** Grid ********************************
Simulation Box:
Type = minimum
Radius [A] = 3.500
Octopus will run in 3 dimension(s).
Octopus will treat the system as periodic in 0 dimension(s).
Main mesh:
Spacing [A] = ( 0.220, 0.220, 0.220) volume/point [A^3] = 0.01065
# inner mesh = 31525
# total mesh = 51605
Grid Cutoff [eV] = 776.921749 Grid Cutoff [Ry] = 57.102700
**********************************************************************
***************************** Symmetries *****************************
Symmetry elements : 4*(C3) 3*(C2) 3*(S4) 6*(sigma)
Symmetry group : Td
**********************************************************************
**************************** Theory Level ****************************
Input: [TheoryLevel = dft]
Exchange-correlation:
Exchange
Slater exchange (LDA)
[1] P. A. M. Dirac, Math. Proc. Cambridge Philos. Soc. 26, 376 (1930)
Correlation
Perdew & Zunger (Modified) (LDA)
Input: [SICCorrection = sic_none]
**********************************************************************
SCF converged in 3 iterations
Some of the states are not fully converged!
Eigenvalues [eV]
#st Spin Eigenvalue Occupation
1 -- -10.617759 2.000000
2 -- -7.350349 2.000000
3 -- -0.625441 2.000000
4 -- 2.813351 2.000000
Energy [eV]:
Total = -138.66600447
Free = -138.66600447
-----------
Ion-ion = 236.08498119
Eigenvalues = -31.56039516
Hartree = 367.70074449
Int[n*v_xc] = -104.40563731
Exchange = -69.09214325
Correlation = -10.80334008
vanderWaals = 0.00000000
Delta XC = 0.00000000
Entropy = 0.00000000
-TS = -0.00000000
Photon ex. = 0.00000000
Kinetic = 233.69335590
External = -896.24957644
Non-local = -56.82391015
Int[n*v_E] = 0.00000000
Dipole: [A] [Debye]
<x> = -1.28092E+00 -6.15254E+00
<y> = 3.40491E-01 1.63545E+00
<z> = -3.29699E-01 -1.58361E+00
Convergence:
abs_dens = 1.88070103E-14 ( 0.00000000E+00)
rel_dens = 2.35087629E-15 ( 1.00000000E-06)
abs_ev = 3.45030665E-01 ( 0.00000000E+00) [eV]
rel_ev = 2.97485561E-01 ( 0.00000000E+00)
Forces on the ions [eV/A]
Ion x y z
1 C 0.362937E+02 -0.275830E+02 0.116276E+02
2 H 0.170877E+02 0.757380E+01 0.150166E+02
3 H -0.650774E+01 -0.134651E+02 0.129681E+02
4 H 0.147686E+02 -0.112452E+02 -0.725906E+01
5 H -0.142230E+02 0.106085E+02 -0.154793E+02
----------------------------------------------------------
Max abs force 0.362937E+02 0.275830E+02 0.154793E+02
Total force 0.474192E+02 -0.341111E+02 0.168739E+02
Total torque -0.635395E+00 -0.162418E+00 0.317413E+01
Where am I going wrong?
Thank you
