QuantumATK Forum

QuantumATK => General Questions and Answers => Topic started by: leslie on September 8, 2013, 17:38

Title: potential averaged over the B direction for a fixed A coordinate
Post by: leslie on September 8, 2013, 17:38

Dear ATK developers and all,

May I know how to extract potential averaged over the B direction for a fixed A coordinate (A = 18 Angstrom) as you ploted in FIG. 1. and FIG. 5. (and also the Mulliken charge) in PHYSICAL REVIEW B 85, 165442 (2012). I would appreciate your help.

Best,
leslie

Title: Re: potential averaged over the B direction for a fixed A coordinate
Post by: Anders Blom on September 8, 2013, 19:00

To begin with, I will here assume that the grid cell is orthogonal. If not, you need to do things a bit differently.

The first thing you want to do is figure out which grid index corresponds to the A-coordinate you want.

Code: python

# Read the data
data = nlread("file.nc")[0]

# Extract the X coordinates of the grid
x_axis_length = data.unitCell()[0][0].inUnitsOf(Angstrom)
Nx = data.shape()[0]
x = numpy.array(range(Nx))*x_axis_length/float(Nx)

# Define desired X-coordinate
x0 = 14    # Angstrom
# Find the index in the grid which is closest BELOW (or at) x0
xi = int(numpy.where(x<=x0)[0][-1])
print xi, x[xi]

Now you take this grid index, average of Y, and plot as function of Z.

Code: python

# Extract the Z coordinates
z_axis_length = data.unitCell()[2][2].inUnitsOf(Angstrom)
Nz = data.shape()[2]
z = numpy.array(range(Nz))*z_axis_length/float(Nz)

# numpy.array removes the unit
averaged_data = numpy.array([numpy.average(data[xi,:,zi]) for zi in range(Nz)])

print "z (Ang)  data (%s)" % data.unit()
for zi in range(Nz):
    print z[zi],averaged_data[zi]