confuse about my results on electron-density-difference and mulliken-population

[fangyongxinxi]

Dear Sir,

I have one question about my resluts on electron-density-difference and mulliken-population for my system.

My system is Nb film on MgO(111) substrate, and the calculation model is built like a sandwich (a few layers of Nb on top of a few layers of MgO, with vacuum layer). I hope to get elctrons redistribution and transformation along the interface, thus,  atomic mulliken population and electron density difference are calculated.

Firstly, we calculate atomic mulliken population for the atoms near interface, taking Nb atoms for example, the difference in atomic mulliken population is
m(Nb) = m(Nb in Nb/MgO structure) - m(Nb in bare Nb slab)
m(Nb in Nb/MgO structure) is the atomic mulliken population of Nb atoms in Nb/MgO structure,
m(Nb in bare Nb slab) is the atomic mulliken population of Nb atoms in bare Nb slab
thus, m(Nb) gives the elctrons transformation for Nb atoms.


secondly, the electron density difference along the interface is calculated
pho(interface) = pho(total) - pho(mgo) - pho(nb)
pho(total) is the total electron density of the whole structure,
pho(mgo) is the electron density of mgo layers (MgO layers as a whole part),
pho(nb) is the electron density of nb layers.
thus, pho(interface) gives electrons redistribution.

the input file, and calculation result are shown in attach files.
(the electron density difference map is viewed by VESTA software).

Results are like follows:
Firstly, for the atomic mulliken population, we see each Nb atoms ( Nb atoms in the first layer in Nb slab) loose 0.25 electrons, and each Oxygen atoms gain 0.24 electrons. It's easy to understand that the Oxygen-terminated MgO(111) gain electrons from Nb slab.

However, the electron density difference map is not agree with mulliken population results,
it shows both nb atoms and oxygen atoms seem loose electrons (more precisely，strong electron redistribution around there atoms) , as the blue color shows, and only nb atoms seem gain electrons (because there is no light red color around oxygen atoms, which means there is no electrons enhance around oxygen atoms).

In my view, oxygen atoms gain electrons, as mulliken population result show, and there should be electrons enhance around oxygen atoms, but the electron density difference result give the opposite result. This is my confuse.

How to understand the difference ?


Thanks.




# the input.txt is like follows:



# Set up lattice
savefile = 'lsda_1.nc'


vector_a = [2.98186929626, 0.0, 0.0]*Angstrom
vector_b = [0.0, 5.16474912266, 0.0]*Angstrom
vector_c = [0.0, 0.0, 38.0247821144]*Angstrom
lattice = UnitCell(vector_a, vector_b, vector_c)

# Define elements
elements = [Magnesium, Magnesium, Oxygen, Oxygen, Magnesium, Magnesium, Oxygen,
            Oxygen, Magnesium, Magnesium, Oxygen, Oxygen, Niobium, Niobium,
            Niobium, Niobium, Niobium, Niobium, Niobium, Niobium, Niobium,
            Niobium, Niobium, Niobium]

# Define coordinates
fractional_coordinates = [[ 0.            ,  0.            ,  0.25611571715 ],
                          [ 0.5           ,  0.5           ,  0.25611571715 ],
                          [ 0.5           ,  0.166666666667,  0.288130181794],
                          [ 0.            ,  0.666666666667,  0.288130181794],
                          [ 0.            ,  0.333333333333,  0.320144646437],
                          [ 0.5           ,  0.833333333333,  0.320144646437],
                          [ 0.            ,  0.            ,  0.352159111081],
                          [ 0.5           ,  0.5           ,  0.352159111081],
                          [ 0.5           ,  0.166666666667,  0.384173575725],
                          [ 0.            ,  0.666666666667,  0.384173575725],
                          [ 0.            ,  0.333333333333,  0.416188040369],
                          [ 0.5           ,  0.833333333333,  0.416188040369], # ############
                          [ 0.            ,  0.            ,  0.455618743112],
                          [ 0.5           ,  0.5           ,  0.455618743112],
                          [ 0.5           ,  0.            ,  0.517272419054],
                          [ 0.            ,  0.5           ,  0.517272419054],
                          [ 0.            ,  0.            ,  0.578926094997],
                          [ 0.5           ,  0.5           ,  0.578926094997],
                          [ 0.5           ,  0.            ,  0.640579770939],
                          [ 0.            ,  0.5           ,  0.640579770939],
                          [ 0.            ,  0.            ,  0.702233446882],
                          [ 0.5           ,  0.5           ,  0.702233446882],
                          [ 0.5           ,  0.            ,  0.763887122824],
                          [ 0.            ,  0.5           ,  0.763887122824]]

# Set up configuration
bulk_configuration = BulkConfiguration(
    bravais_lattice=lattice,
    elements=elements,
    fractional_coordinates=fractional_coordinates
    )



# Set up configuration
bulk_configuration = BulkConfiguration(
    bravais_lattice=lattice,
    elements=elements,
    fractional_coordinates=fractional_coordinates
    )

#----------------------------------------
# Basis Set
#----------------------------------------
basis_set = [
    LDABasis.Oxygen_SingleZetaPolarized,
    LDABasis.Magnesium_SingleZetaPolarized,
    LDABasis.Niobium_SingleZetaPolarized,
    ]

exchange_correlation = LSDA.PZ

numerical_accuracy_parameters = NumericalAccuracyParameters(
    k_point_sampling=(5, 5, 1),
    )

calculator = LCAOCalculator(
    basis_set=basis_set,
    exchange_correlation=exchange_correlation,
    numerical_accuracy_parameters=numerical_accuracy_parameters,
    )


bulk_configuration.setCalculator(calculator)
nlprint(bulk_configuration)
bulk_configuration.update()
nlsave(savefile, bulk_configuration)

# -------------------------------------------------------------
# Total energy
# -------------------------------------------------------------
total_energy = TotalEnergy(bulk_configuration)
nlsave(savefile, total_energy)
nlprint(total_energy)

# -------------------------------------------------------------
# Mulliken population
# -------------------------------------------------------------
mulliken_population = MullikenPopulation(bulk_configuration)
nlsave(savefile, mulliken_population)
nlprint(mulliken_population)

# -------------------------------------------------------------
# Electron density
# -------------------------------------------------------------
total = ElectronDensity(
    configuration=bulk_configuration,
    )
nlsave(savefile, total)




vector_a = [2.98186929626, 0.0, 0.0]*Angstrom
vector_b = [0.0, 5.16474912266, 0.0]*Angstrom
vector_c = [0.0, 0.0, 38.0247821144]*Angstrom
lattice = UnitCell(vector_a, vector_b, vector_c)

# Define elements
elements = [Magnesium, Magnesium, Oxygen, Oxygen, Magnesium, Magnesium, Oxygen,
            Oxygen, Magnesium, Magnesium, Oxygen, Oxygen]

# Define coordinates
fractional_coordinates = [[ 0.            ,  0.            ,  0.25611571715 ],
                          [ 0.5           ,  0.5           ,  0.25611571715 ],
                          [ 0.5           ,  0.166666666667,  0.288130181794],
                          [ 0.            ,  0.666666666667,  0.288130181794],
                          [ 0.            ,  0.333333333333,  0.320144646437],
                          [ 0.5           ,  0.833333333333,  0.320144646437],
                          [ 0.            ,  0.            ,  0.352159111081],
                          [ 0.5           ,  0.5           ,  0.352159111081],
                          [ 0.5           ,  0.166666666667,  0.384173575725],
                          [ 0.            ,  0.666666666667,  0.384173575725],
                          [ 0.            ,  0.333333333333,  0.416188040369],
                          [ 0.5           ,  0.833333333333,  0.416188040369],
]

# Set up configuration
bulk_configuration = BulkConfiguration(
    bravais_lattice=lattice,
    elements=elements,
    fractional_coordinates=fractional_coordinates
    )



# Set up configuration
bulk_configuration = BulkConfiguration(
    bravais_lattice=lattice,
    elements=elements,
    fractional_coordinates=fractional_coordinates
    )

basis_set = [
    LDABasis.Oxygen_SingleZetaPolarized,
    LDABasis.Magnesium_SingleZetaPolarized,
    LDABasis.Niobium_SingleZetaPolarized,
    ]

exchange_correlation = LSDA.PZ

numerical_accuracy_parameters = NumericalAccuracyParameters(
    k_point_sampling=(5, 5, 1),
    )

calculator = LCAOCalculator(
    basis_set=basis_set,
    exchange_correlation=exchange_correlation,
    numerical_accuracy_parameters=numerical_accuracy_parameters,
    )

bulk_configuration.setCalculator(calculator)
nlprint(bulk_configuration)
bulk_configuration.update()
nlsave(savefile, bulk_configuration)

# -------------------------------------------------------------
# Total energy
# -------------------------------------------------------------
total_energy = TotalEnergy(bulk_configuration)
nlsave(savefile, total_energy)
nlprint(total_energy)

# -------------------------------------------------------------
# Mulliken population
# -------------------------------------------------------------
mulliken_population = MullikenPopulation(bulk_configuration)
nlsave(savefile, mulliken_population)
nlprint(mulliken_population)

# -------------------------------------------------------------
# Electron density
# -------------------------------------------------------------
mgo = ElectronDensity(
    configuration=bulk_configuration,
    )
nlsave(savefile, mgo)




vector_a = [2.98186929626, 0.0, 0.0]*Angstrom
vector_b = [0.0, 5.16474912266, 0.0]*Angstrom
vector_c = [0.0, 0.0, 38.0247821144]*Angstrom
lattice = UnitCell(vector_a, vector_b, vector_c)

# Define elements
elements = [Niobium, Niobium,
            Niobium, Niobium, Niobium, Niobium, Niobium, Niobium, Niobium,
            Niobium, Niobium, Niobium]

# Define coordinates
fractional_coordinates = [
                          [ 0.            ,  0.            ,  0.455618743112],
                          [ 0.5           ,  0.5           ,  0.455618743112],
                          [ 0.5           ,  0.            ,  0.517272419054],
                          [ 0.            ,  0.5           ,  0.517272419054],
                          [ 0.            ,  0.            ,  0.578926094997],
                          [ 0.5           ,  0.5           ,  0.578926094997],
                          [ 0.5           ,  0.            ,  0.640579770939],
                          [ 0.            ,  0.5           ,  0.640579770939],
                          [ 0.            ,  0.            ,  0.702233446882],
                          [ 0.5           ,  0.5           ,  0.702233446882],
                          [ 0.5           ,  0.            ,  0.763887122824],
                          [ 0.            ,  0.5           ,  0.763887122824]]

# Set up configuration
bulk_configuration = BulkConfiguration(
    bravais_lattice=lattice,
    elements=elements,
    fractional_coordinates=fractional_coordinates
    )



# Set up configuration
bulk_configuration = BulkConfiguration(
    bravais_lattice=lattice,
    elements=elements,
    fractional_coordinates=fractional_coordinates
    )

#----------------------------------------
# Basis Set
#----------------------------------------
basis_set = [
    LDABasis.Oxygen_SingleZetaPolarized,
    LDABasis.Magnesium_SingleZetaPolarized,
    LDABasis.Niobium_SingleZetaPolarized,
    ]

exchange_correlation = LSDA.PZ

numerical_accuracy_parameters = NumericalAccuracyParameters(
    k_point_sampling=(5, 5, 1),
    )

calculator = LCAOCalculator(
    basis_set=basis_set,
    exchange_correlation=exchange_correlation,
    numerical_accuracy_parameters=numerical_accuracy_parameters,
    )

bulk_configuration.setCalculator(calculator)
nlprint(bulk_configuration)
bulk_configuration.update()
nlsave(savefile, bulk_configuration)

# -------------------------------------------------------------
# Total energy
# -------------------------------------------------------------
total_energy = TotalEnergy(bulk_configuration)
nlsave(savefile, total_energy)
nlprint(total_energy)

# -------------------------------------------------------------
# Mulliken population
# -------------------------------------------------------------
mulliken_population = MullikenPopulation(bulk_configuration)
nlsave(savefile, mulliken_population)
nlprint(mulliken_population)

# -------------------------------------------------------------
# Electron density
# -------------------------------------------------------------
nb = ElectronDensity(
    configuration=bulk_configuration,
    )
nlsave(savefile, nb)

final = total - mgo - nb
nlsave(savefile, final)

[zh]

The discrepancy for the charge of atom analyzed from the basis set  and from the real space division is quite common, especially for the charge estimated from the Mulliken population, in which the charge is strongly dependent on the basis set of atom used in the analysis.

In your case, you used the SZP basis set, which might not give reliable results.

[fangyongxinxi]

ok, thanks.

[fangyongxinxi]

The discrepancy for the charge of atom analyzed from the basis set  and from the real space division is quite common, especially for the charge estimated from the Mulliken population, in which the charge is strongly dependent on the basis set of atom used in the analysis.

In your case, you used the SZP basis set, which might not give reliable results.

Dear,
I try DZP and DZDP, they all show the samilar  result as I used szp. so, I am still confused.
Any other suggestion ?
Thanks.

[fangyongxinxi]

Problem solved,
I change lsda.pz to lda.pz,
and I get the result that I expected.

[zh]

Great.  But you may check the "Spin.Sum" component for Mulliken population and ElectronDensity in the LSDA calculations and see whether it gives the consistent results as you expected.