8
$\begingroup$

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

$\endgroup$
2
  • 2
    $\begingroup$ Not an Octopus expert, but it looks like the calculation didn't converge. It only took 3 iterations, which sounds implausibly small - unless you already have a guess that is extremely close to convergence. The log file also says that some of the states aren't converged, and errors of 0.3 seem pretty damn high to me. Maybe there's a problem with a too high convergence threshold, or the density isn't updating at all? $\endgroup$ Feb 9 at 19:53
  • 1
    $\begingroup$ Was Susi's comment helpful? Did you end up figuring this problem out? $\endgroup$ Aug 14 at 16:57

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.