Why can not calculate the ElectronDensity of Bulkconfiguration using ATK-SE?

[esp]

>> Even if your system is 2D, the electrons are not, the electron cloud has a finite extension beyond the graphene plane, perpendicular to it. The difference you compute is for the entire simulation cell, and so it is normalized by the cell volume. That is, if you ran the simulation with twice as much vacuum around the graphene sheet (which obviously should not affect the result, provided there is "enough" vacuum in both cases), the density values drop by a factor 2. Hence it is appropriate to multiplying by dimension in which you have vacuum to get a 2D density. I'm not sure that's what you do - I would probably just multiply by the cell length perpendicular to the sheet. All makes sense, except the last part ... I am not sure what is the cell length? >> For the integral over dx, I think you have just misread the formulas. It's not a piecewise integration over x, it's an integral over dydz for each value of x. So the result is a function of x still. What you write d(x) = int { n(x,y,z)dydzdx } makes no mathematical sense. In fact, int { n(x,y,z)dydzdx } is a pure number, viz. the number of electrons in the system. And that is indeed equal to the .sum() you used, provided of course you also multiply by the volume element. Ok, above, great, i understand that .. but then why is this: s1 = diffDensity1[:, :, :].sum() Not equal to the result from this loop (sum1)?

Code

        dz1 = dZ1.norm()
	sum1 = 0 * Units.Bohr**-2
	for i in range(shape1[2]):
		#print dz1*i, n_z1[i]
		sum1+=n_z1[i]*dz1
	print 'total density1=', sum1

In this loop above, n_z1 is the already calculated integral of ElectronDifferenceDensity over X and Y (in my adaptation of the example)  ... so just like your example on the ElectronDifferenceDensity manual page, this is n(z) = int { n(x,y,z)dxdy } ... Then this loop above, integrates over z right? ...  so again my confusion why this is not equal to s1 = diffDensity1[:, :, :].sum()?

[Anders Blom]

The unit cell for the system has a length in all three directions, that's what I meant.

An integral is not only a sum of the values of f(x), it's a sum of f(x)*dx.

[esp]

ok i am almost there but let me go line by line to make sure i understand .. first of all let me say that i realize for what i want, i can use sum(), that is good enough ... but i still want to understand this calculation because it could be useful .... please stop me where i am making a mistake 1) ok read in DD at bias1: diffDensity1 = nlread(fName, ElectronDifferenceDensity, object_id="diffDensity[vg-0.25000][vds0.50000]")[0] 2) this gets dimensions of array? shape1 = diffDensity1.shape() 3) this just separates components? dX1, dY1, dZ1 = diffDensity1.volumeElement() 4) this would get unit area in XY plane? dAYX1 = numpy.linalg.norm( numpy.cross(dY1,dX1) ) 5) this calculates the integral n(z) = int { d(x,y,z)dydx } for each value of z, and stores them in the n_z1 array, indexed by the z index ... so this is the density, per unit area in XY, at each value of z right? n_z1 = [ diffDensity1[:,:,i].sum() * dAYX1 for i in range(shape1[2]) ] 6) this is where my confusion is i think: here we now take each value in z, and we sum them all up multiplying by dz ... so sum1 = int{ n(z)dz } ...  and if it is, then why is this not equivalent to: int{ n(z)dz } = int{ int { d(x,y,z)dydx } dz } = int{ d(x,y,z)dydxdz } ? 

Code

	dz1 = dZ1.norm()
	sum1 = 0 * Units.Bohr**-2
	for i in range(shape1[2]):
		#print dz1*i, n_z1[i]
		sum1+=n_z1[i]*dz1
	print 'total density1=', sum1

If all above are true, then the value of sum1 should be the same as: s1 = diffDensity1[:, :, :].sum() I am sorry, I am not sure where I am going wrong here ,,.. but i do get two different answers .. please tell which step in my explanation is wrong

[Anders Blom]

All is correct, almost

3) Not quite, but by 4) you are on the right track again.

6) Indeed, this is the small misconception. .sum() is different from the integral by a factor equal to the volume element. To check this, remove dAXY1 from the integral over XY, and remove dz1 from the integral over Z, then you get the same values in both methods. Or, conversely, multiply the .sum() by numpy.dot(dZ1,cross(dX1,dY1)).

[Anders Blom]

For a very trivial example:

Code: python

# Estimate the integral of 3x^2 from 0 to 1 (should be 1):
dx = 0.0001
# Generate x-values from 0 to 1 (end-points included, via dx/2)
x = numpy.arange(0,1+dx/2.0,dx)
y = 3*x**2
print sum(y)*dx

Depending on dx, the value is close to 1. sum(y) is however not 1 - far from it, it becomes larger and larger with smaller dx. This is your sum(), which thus is not the same as the integral.

[esp]

i get it now .. sum() is a density sum, not the number of electrons ... i missed that part ..

So then i think this is correct:

   s1 = diffDensity1[:, :, :].sum()
   s2 = diffDensity2[:, :, :].sum()
   dX1, dY1, dZ1 = diffDensity1.volumeElement()
   dX2, dY2, dZ2 = diffDensity2.volumeElement()
   s1 = numpy.dot(dZ1,numpy.cross(dX1,dY1))*s1.inUnitsOf(Bohr**-3)
   s2 = numpy.dot(dZ2,numpy.cross(dX2,dY2))*s2.inUnitsOf(Bohr**-3)
        print "Total # electrons: ", s1, s2   
        numElecDiff = (s2-s1)
   print "numElecDiff", numElecDiff

Where the end answer is the number of electrons added from the increase in voltage, since densities 1 and 2 were at lower and higher voltages, respectively .. but something is not right here, because the # of electrons is less than 1 over the entire volume ... or is that per unit volume? if per unit volume, how do i get the number of electrons ... ?

[Anders Blom]

Yes, it's per unit volume. The number of electrons is obtained by multiplying with the volume of the unit cell.

[esp]

how do i get the entire volume? remember i need charge in the channel not charge in a unit volume

Ok, I did the following ... if i have the density per unit volume, and it seems the unit length is about 0.7 Bohr in the x direction, perpendicular to the 2d plane, some very rough estimate of channel charge (for 2D graphene only) could be to multiply this by the channel length, and its width .... this would assume that the 0.7 Bohr in the x-direction is a large enough "thickness" to encapsulate the charge of the Carbon atoms, which it is not ... but in the absence of knowing how to get the actual volume, ..  very roughly ...

i took my density per unit volume, and multiplied by 160 angstrom (16 nm channel length), and 12.9 angstrom (1.29 nm channel width) .. and i redid this for -0.25V (totaly off) and 0.0 volts (halfway on), and i get

Total # electrons:  before: 1756.83419532 after: 3626105.77915
numElecDiff 3624348.94495
capacitance in pF:  2.31958332477

Well now this sounds more reasonable .. but i know it is not actually correct ... if i could know how to get the actual volume of the ribbon that includes all the electrons, i think i would be set ... now another question is too .. when you run ElectronDifferenceDensity, is it only the scattering region, or the whole device?

[esp]

also i would not be using off and half on as the points, but thats all i have for now ... i just want to get the calculation right first as i wait for results

[Anders Blom]

The electron density is normalized by the actual volume of the simulation cell, the volume of which you can compute by hand by just looking at the configuration (you know the Bravais lattice or UnitCell), or you can do

Code: python

cA, cB, cC = density_object.unitCell()
cell_volume = numpy.dot(cC, numpy.cross(cA,cB))*Bohr**3

[Anders Blom]

Cells? There is only one cell. What's important to understand is that if you integrate the density (not the difference) over the entire volume of the simulation cell you get the number of electrons. That gives you the normalization. I may be wrong regarding the need to multiply your result with the volume, however. Just check that you can integrate the densities (not their difference) to give the number of electrons, then you know it's all normalized.

[esp]

Ok let's try one more time, two difference densities at two different voltages:

Code

	diffDensity1 = nlread(fName, ElectronDifferenceDensity, object_id="diffDensity[vg-0.25000][vds0.50000]")[0]
	diffDensity2 = nlread(fName, ElectronDifferenceDensity, object_id="diffDensity[vg0.00000][vds0.50000]")[0]
	
        s1 = diffDensity1[:, :, :].sum()
	s2 = diffDensity2[:, :, :].sum()
	
        cA1, cB1, cC1 = diffDensity1.unitCell()
	dV1 = numpy.dot(cC1, numpy.cross(cA1,cB1))
	cA2, cB2, cC2 = diffDensity2.unitCell()
	dV2 = numpy.dot(cC2, numpy.cross(cA2,cB2))

	numE1b = s1.inUnitsOf(Angstrom**-3)*dV1
	numE2b = s2.inUnitsOf(Angstrom**-3)*dV2
	print "Total # electrons: ", numE1b, numE2b

        numElecDiff = (numE2b-numE1b)
 	print "numElecDiff", numElecDiff

is numElecDiff the number of electrons that were injected because of raising the voltage?

[Anders Blom]

Start by checking the densities individually, not the difference (since we don't know what to expect there).

But, no, I would use the previous scheme with the volume element, that makes more sense for the normalization, to get the number of electrons.

The reason to involve the unit cell will be separate; after all, you are probably not too interested in the absolute number of electrons, but probably the number per unit area (or unit length, perhaps).

[esp]

number of electrons per unit length in transprt direction would be great ... if only i could decipher how to do it

[Anders Blom]

Are you sure it should be in the transport direction? Is it a ribbon or a sheet you are working with?