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

[eastnobil]

Why can not calculate the ElectronDensity of Bulk configuration using ATK-SE? But I found the ElectronDifferent Density can be calculated. see attachment as following.
What is the meaning of ElectronDifferentDensity if the Electron Density can not be calculated?

[kstokbro]

You can calculate the electron density with ATK-DFT, however, electron density is not supported fro ATK-SE.

Attached is from the reference manual:
ElectronDifferenceDensity, returns the electron difference density, i.e. the difference between the self-consistent valence charge density and the superposition of atomic valence densities.

[esp]

I need to calculate the "charge" in the channel to try to estimate capacitance for a GNR device ... from some relevant papers, i see the technique many use use is to calculate a small change in charge with respect to gate voltage (for Cgate for example) .. Q = CV, C = dQ/dV ... now, the idea i found in another post was to use the ElecronDensity, then add it up for example over the entire channel to get charge .. I was planning on calculating Q1 at V1, then Q2 and V2, then C = (Q2-Q1)/(V2-V1) .... does this all sound reasonable?

Now, I have been running ElectronDensity all week .. I cannot get it to converge .. my current job is 36 hours running, still on the base calculation .. no points completed yet ... it is taking forever ... i even reduced the k points a bit and the energy range ... it never got to the Elec.Dens. calc yet ...   [Update:  As Dr. Blom noted, it is the base calculation taking a long time to converge, not the E.D. calculation]

So, I am wondering, I believe i tried before and electron difference density was faster ... is this expected?  if so, then wouldn't the method above of calculating difference in charge at two gate voltages give the same answer?  In other words if i am interested in the difference at two different gate voltages and not the absolute value of charge or elec dens, then wouldnt both methods give the same answer, if i used either ElectronDensity, or ElectronDifferenceDensity?

[Anders Blom]

The electron density calculation itself is a simple step after the main calculation is converged. So, one cannot say that the electron density doesn't converge.

So, I would assume you have trouble converging the self-consistent device calculation. But without knowing anything about that, it's impossible to provide any guidance beyond the generic, which I think you know very well already.

But this really has nothing to do with which quantity you attempt to compute, since as you mention, you apparently never even get to this point in the script. To reiterate: ATK first converges the self-consistent calculation, and then, afterwards, computes stuff like band structure, potentials, LDOS, and electron density. But all of that is not taken into consideration for the self-consistent calculation, nor are (typically) the post-processing steps in any way difficult (time-consuming, perhaps, but there is a progress bar).

Finally, indeed the difference between two ElectronDensity objects and two ElectronDifferenceDensity objects is identical.

[esp]

ok thank you very much

[esp]

i would add though that it seems that ElectronDifferenceDendity uses Huckel, while ElectronDensity uses DFT, which in my case seems to be why the base calculation converges much more easily for EDD than ED ... I will retry with Huckel

[Anders Blom]

No, you cannot say it this way. If you have a Huckel calculation, then you can only ask for the difference density, if you have DFT you can ask for both. It's not the other way around (i.e. the density calculators do not "use" any method).

[esp]

i understand, thank you ...

[esp]

Ok so i got some results, ran Electron Difference Density at two different gate voltages ... then i did this:

   diffDensity1 = nlread(fName, ElectronDifferenceDensity, object_id="diffDensity[vg-0.25000][vds0.50000]")[0]
   diffDensity2 = nlread(fName, ElectronDifferenceDensity, object_id="diffDensity[vg-0.20000][vds0.50000]")[0]

   s1 = diffDensity1[:, :, :].sum()
   s2 = diffDensity2[:, :, :].sum()
   print s1, s2, s2-s1

I got this:

EDD at V1:  2.49328980937 1/Bohr**3
EDD at V2:  2.62952644165 1/Bohr**3
EDD difference:  0.136236632277 1/Bohr**3

So does this mean 0.136 electrons per Bohr^3 ? Also I cannot figure out how to convert to 1/Angstrom**3 .. I have tried inUnitsOf() in different forms but it did not work...

also, since this is graphene .. volume does not seem to make sense .. how do i get a 2d calculation, or is it already 2d anyway, since there is only 2d graphene?

[esp]

ok i did this, does it make sense?

Code

	diffDensity1 = nlread(fName, ElectronDifferenceDensity, object_id="diffDensity[vg-0.25000][vds0.50000]")[0]
	diffDensity2 = nlread(fName, ElectronDifferenceDensity, object_id="diffDensity[vg-0.20000][vds0.50000]")[0]

	s1 = diffDensity1[:, :, :].sum()
	s2 = diffDensity2[:, :, :].sum()
	print s1, s2, s2-s1

	shape1 = diffDensity1.shape()
	shape2 = diffDensity2.shape()

	# Find the volume elements.
	dX1, dY1, dZ1 = diffDensity1.volumeElement()
	dX2, dY2, dZ2 = diffDensity2.volumeElement()

	# Calculate the unit area in the y-x plane.
	dAYX1 = numpy.linalg.norm( numpy.cross(dY1,dX1) )
	dAYX2 = numpy.linalg.norm( numpy.cross(dY2,dX2) )

	# calculate density along z integrated over y,x
	n_z1 = [ diffDensity1[:,:,i].sum() * dAYX1 for i in range(shape1[2]) ]
	n_z2 = [ diffDensity2[:,:,i].sum() * dAYX2 for i in range(shape2[2]) ]

	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

	dz2 = dZ2.norm()
	sum2 = 0 * Units.Bohr**-2
	for i in range(shape2[2]):
		print dz2*i, n_z2[i]
		sum2+=n_z2[i]*dz2
	print 'total density2=', sum2

	print s1, s2, s2-s1
	print 'diff', sum2-sum1

I get: 0.0465095588926 1/Bohr**2 now i think this means, in a slice of X-Y plane, this is the density difference (difference) per area, due to the change in voltage ... correct?  and so to translate this over the whole channel in the Z direction .. i would have to do what exactly? 

[esp]

actually the documentation page does not make sense to me ..

http://quantumwise.com/documents/manuals/latest/ReferenceManual/index.html/ref.electrondifferencedensity.html

it says:
Read in the electron difference density from a file and print d(x) = int { n(x,y,z)dydz }

but the calculation looks like it first calculates d(x) as above, then integrates that over dx, so it seems it calculates the integral over all x,y,z isnt that correct?  then why does it say print out d(x) = int { n(x,y,z)dydz } when it is printing out d(x) = int { n(x,y,z)dydzdx } ??

Update: ok i see that is referring to the intermediate print statements before the sum is printed .. but then still the question is, why is sum not equal to a simple sum like this, without doing the whole piece by piece integration?
   s1 = diffDensity1[:, :, :].sum()

[esp]

ok here is a full script, does it makes sense for calcaulting gate capacitance in channel of GNR? 

Code

	diffDensity1 = nlread(fName, ElectronDifferenceDensity, object_id="diffDensity[vg-0.25000][vds0.50000]")[0]
	diffDensity2 = nlread(fName, ElectronDifferenceDensity, object_id="diffDensity[vg-0.20000][vds0.50000]")[0]

	s1 = diffDensity1[:, :, :].sum()
	s2 = diffDensity2[:, :, :].sum()
	print s1, s2, s2-s1

	shape1 = diffDensity1.shape()
	shape2 = diffDensity2.shape()

	# Find the volume elements.
	dX1, dY1, dZ1 = diffDensity1.volumeElement()
	dX2, dY2, dZ2 = diffDensity2.volumeElement()

	# Calculate the unit area in the y-x plane.
	dAYX1 = numpy.linalg.norm( numpy.cross(dY1,dX1) )
	dAYX2 = numpy.linalg.norm( numpy.cross(dY2,dX2) )

	print dAYX1, dAYX2

	# calculate density along z integrated over y,x
	n_z1 = [ diffDensity1[:,:,i].sum() * dAYX1 for i in range(shape1[2]) ]
	n_z2 = [ diffDensity2[:,:,i].sum() * dAYX2 for i in range(shape2[2]) ]

	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

	dz2 = dZ2.norm()
	sum2 = 0 * Units.Bohr**-2
	for i in range(shape2[2]):
		#print dz2*i, n_z2[i]
		sum2+=n_z2[i]*dz2
	print 'total density2=', sum2

	print 'diff', s1, s2, s2-s1
	print 'diff', sum1, sum2, sum2-sum1

	# let's use s2-s1 as difference in # electrons in channel due to voltage
	numElecDiff = (s2-s1).inUnitsOf(Bohr**-3)
	print "numElecDiff", numElecDiff
	# convert to coulombs
	chargeDiff = numElecDiff*1.6e-19
	voltDiff = 0.05
	capacitance = chargeDiff/voltDiff
	print "capacitance in F", capacitance

the voltage difference and voltages i used would need to be chosen carefully, and in fact i am going to average over two cap values found at off end, and at on end ... fyi

[Anders Blom]

You're fast on the trigger with the posts, so I have to answer out of sequence Let's see if I cover all points. First of all, you will find this useful, for converting units:

Code: python

x=(1/Bohr)**3
x.inUnitsOf(Ang**-3)

which prints 6.748342563230854. 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. 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.

[esp]

ok thank you ...

So then why is the integral method not exactly equal to the .sum() method ... ?

[Anders Blom]

You mean d(x) or the point about the volume element?