Function for charge density difference

[vihardabest]

Dear all,
I just want to ask you if it is possible to be create a script for charge redistrubution.


What I mean:
The transfered charged and the induced change in the electrostatic potential are obtained  by taking the difference between the self-consistent charge/potential  distribution in the equilibrium metal-molecule-metal junction and the charge/potential distribution in the isolated molecule plus the bare bimetallic junction.(Phys.Rev. B 69,085403 (2004))

I know that here we have a function written by ipsecog which calculated density difference.  I successfully used this function to calculated the density difference in different bias.
I trued to make the same for the charged redistribution for the : (metal-molecule-metal) - (molecule) and the error is :

File "elec_density-diffrence.py", line 22, in ?
    diff_density = calculateDensityDifference(density1,density2)
  File "elec_density-diffrence.py", line 15, in calculateDensityDifference
    density._ElectronDensity__electron_density_data = (n1-n2)
ValueError: shape mismatch: objects cannot be broadcast to a single shape


I know that is something wrong with the formats of the input files and I suspect that the reason is that for the metal-molecule-metal file I used Atk.TwoProbe method to calculate while for the molecule I used  ATK.KohnSham and according to me the size of the two files doesn't match. After this I made them text files but they are massive and impossible to handle.

I will really appreciate your help and ideas how to get the charge density difference.


----------------------------- This is my script-----------------------------------

from ATK.TwoProbe import *


def calculateDensityDifference (d1,d2):

    import copy

    n1 = d1.toArray()
    n2 = d2.toArray()
    density = copy.copy(d1)
    density._ElectronDensity__electron_density_data = (n1-n2)
    return density

scf1 = restoreSelfConsistentCalculation ('metal-molecule-metal.nc')
scf2 = restoreSelfConsistentCalculation ('molecule-scf.nc')
density1 = calculateElectronDensity(scf1)   # Two-probe
density2 = calculateElectronDensity(scf2)   # molecule 0.0 V bias
diff_density = calculateDensityDifference(density1,density2)

f = VNLFile('elec_density-difference.vnl')
f.addToSample(atomic_configuration,'elec_density-difference')
f.addToSample(density1,'elec_density-difference','Two-probe')
f.addToSample(density2,'elec_density-difference','Molecule')
f.addToSample(diff_density,'elec_density-difference','Difference density')

[Anders Blom]

Without being able to look closer at the data, I'd guess the two arrays have different sizes, i.e. a different mesh cut-off and/or a different size unit cell was used in the two calculations.

To verify this, you can "print n1.shape, n2.shape" after the toArray statements.

To compute the difference in this way, you'd have to use identical grids in the two calculations.

[vihardabest]

Thank you very much Anders for the fast replay.
I will try to understand what is going on.

[vihardabest]

Hi Anders, it is me again.
I printed the shapes of the two files

"print n1.shape, n2.shape"
(2, 113, 113, 468)
(2, 113, 113, 468)


for me they look the same. But it doesn't work again. 
Is it possible to send you the two *.py input files and have a fast look at them when you have time?
Thank you very much in advance.