Unexpected behavior when updating 3D array
The real problem is that I am trying to take a list of the positions of atoms in a unit cell and build the corresponding supercell given the input of the number of times the cell is repeated in each direction.
The result is a loop structure:
#Create Every Unit Cell in SuperCell
aNum=2
bNum=2
cNum=2
atomPos = copy.deepcopy(atomPositions)
for l in range(len(atomPos)):
index=0
for i in range(cNum):
for j in range(bNum):
for k in range(aNum):
for _ in range(numEachAtom[l]):
atomPositions[l][index][0] = atomPos[l][index][0] + 1*k
atomPositions[l][index][1] = atomPos[l][index][1] + 1*j
atomPositions[l][index][2] = atomPos[l][index][2] + 1*i
print atomPositions[0][0]
index += 1
Where atomPositions
is a 3D array such that:,
atomPositions[atomtype=l][atom=index][atomposition=0] = [x,y,z]
and the print instruction is for diagnostics.
The problem is that from the print statement it seems to atomPositions[0][0]
change more often than the atom type and the index seems to be updating correctly, I just don't get it.
Why atomPositions[0][0]
does it change more often than l
?
I originally ran into the problem of modifying the list you are iterating over, hence deep copy at the beginning. Any comments would be much appreciated!
PS This is my first question asking a question and I am a complete Python newbie, so please feel free to comment on my lack of formatting / clarity / style, etc.
Edit:
Example of starting output for numEachAtom=[4,6]
:
[0.17611251675504244, 0.17611251675504244, 0.17611251675504244]
[0.17611251675504244, 0.17611251675504244, 0.17611251675504244]
[0.17611251675504244, 0.17611251675504244, 0.17611251675504244]
[0.17611251675504244, 0.17611251675504244, 0.17611251675504244]
[1.1761125167550424, 0.17611251675504244, 0.17611251675504244]
[1.1761125167550424, 0.17611251675504244, 0.17611251675504244]
[1.1761125167550424, 0.17611251675504244, 0.17611251675504244]
[1.1761125167550424, 0.17611251675504244, 0.17611251675504244]
[0.17611251675504244, 1.1761125167550424, 0.17611251675504244]
[0.17611251675504244, 1.1761125167550424, 0.17611251675504244]
[0.17611251675504244, 1.1761125167550424, 0.17611251675504244]
and etc.
Edit: How atomic objects are initialized:
#Separate out positions of different atom types
atomPositions = []
counter=0
for i,atomN in enumerate(numEachAtom):
atomPositions.append( origAtomPositions[counter:counter+atomN] )
counter += atomN
Edit: Dependencies, etc.
import sys
import scipy as sp
from scipy import *
import copy
Edit: atomPositions [0] is initialized as
[[ 0.17611252 0.17611252 0.17611252]
[ 0.32388748 0.32388748 0.32388748]
[ 0.42611252 0.42611252 0.42611252]
[ 0.07388748 0.07388748 0.07388748]
[ 0.17611252 0.17611252 0.17611252]
[ 0.32388748 0.32388748 0.32388748]
[ 0.42611252 0.42611252 0.42611252]
[ 0.07388748 0.07388748 0.07388748]
[ 0.17611252 0.17611252 0.17611252]
[ 0.32388748 0.32388748 0.32388748]
[ 0.42611252 0.42611252 0.42611252]
[ 0.07388748 0.07388748 0.07388748]
[ 0.17611252 0.17611252 0.17611252]
[ 0.32388748 0.32388748 0.32388748]
[ 0.42611252 0.42611252 0.42611252]
[ 0.07388748 0.07388748 0.07388748]
[ 0.17611252 0.17611252 0.17611252]
[ 0.32388748 0.32388748 0.32388748]
[ 0.42611252 0.42611252 0.42611252]
[ 0.07388748 0.07388748 0.07388748]
[ 0.17611252 0.17611252 0.17611252]
[ 0.32388748 0.32388748 0.32388748]
[ 0.42611252 0.42611252 0.42611252]
[ 0.07388748 0.07388748 0.07388748]
[ 0.17611252 0.17611252 0.17611252]
[ 0.32388748 0.32388748 0.32388748]
[ 0.42611252 0.42611252 0.42611252]
[ 0.07388748 0.07388748 0.07388748]
[ 0.17611252 0.17611252 0.17611252]
[ 0.32388748 0.32388748 0.32388748]
[ 0.42611252 0.42611252 0.42611252]
[ 0.07388748 0.07388748 0.07388748]]
source to share
I use both itertools
and numpy
here, two very useful libraries that you should familiarize yourself with. Question:
take a list of atomic positions in a unit cell and build the corresponding supercell, given the input of how many times to repeat the cell in each direction.
Below is the general answer for any measurement. For illustrative purposes I am using a 2D array, but it will work fine if you pass a 3D array to it. The square loop you are using in your post is replaced with one loop across all possible cell directions with product
. Divide it and play with it to see how it works.
from numpy import *
from itertools import product
# Pass in the points, plus a vector that indicates the repeats in each direction
def supercell(R, v):
v = array(v) + 1
n,d = R.shape
# Construct the copy directions
CV = list(product(*map(range,v)))
R2 = zeros((len(CV)*n,d))
for i,cell in enumerate(CV):
R2[i*n:(i+1)*n,:] = cell+R
return R2
# Construct some random points within a unit cell to work with
N,dimension = 100, 2
# Contrain them so you can see the supercell
R = random.random((N,dimension))*.6 + .5
R2 = supercell(R,[2,3])
# Plot the results
from pylab import *
scatter(*zip(*R2),color='b')
scatter(*zip(*R),color='r')
show()
In your example with a copy of 2 in each positive direction, the vector CV
would look like this:
[(0, 0, 0), (0, 0, 1), (0, 0, 2), (0, 1, 0), (0, 1, 1), (0, 1, 2), (0, 2, 0), (0, 2, 1), (0, 2, 2), (1, 0, 0), (1, 0, 1), (1, 0, 2), (1, 1, 0), (1, 1, 1), (1, 1, 2), (1, 2, 0), (1, 2, 1), (1, 2, 2), (2, 0, 0), (2, 0, 1), (2, 0, 2), (2, 1, 0), (2, 1, 1), (2, 1, 2), (2, 2, 0), (2, 2, 1), (2, 2, 2)]
this is the magic of the function product
.
source to share
atomPositions[0][0]
will change a couple of times, which does yours for _ in range(numEachAtom[l])
.
What is this cycle? Can't you just write
#Create Every Unit Cell in SuperCell
aNum=2
bNum=2
cNum=2
atomPos = copy.deepcopy(atomPositions)
for l in range(len(atomPos)):
index=0
for i in range(cNum):
for j in range(bNum):
for k in range(aNum):
atomPositions[l][index][0] = atomPos[l][index][0] + numEachAtom[l]*k
atomPositions[l][index][1] = atomPos[l][index][1] + numEachAtom[l]*j
atomPositions[l][index][2] = atomPos[l][index][2] + numEachAtom[l]*i
print atomPositions[0][0]
Edit: op just updated his code where he inserted the index into the loop, so it doesn't make sense anymore.
source to share