0

The code above currently produces a cube, not a sphere. How can I correct this? 

x = list(range(-5,5))
y = list(range(-5,5))
z = list(range(-5,5))
r = [1,2,3,4,5,6,7,10]
xcoords = []
ycoords = []
zcoords = []
"""functions used"""
def square(number):
return number**2
"""determines which values satisfies the equation for a sphere"""
for itemx in x:
    for itemy in y:
        for itemz in z:
            for itemr in r:
                if abs(itemx) == abs(itemy) and abs(itemy) == abs(itemz) and square(itemx) + square(itemy) +square(itemz) == itemr:
                xcoords.append(itemx)
                ycoords.append(itemy)
                zcoords.append(itemz)
"""determines the number of atoms in the system"""
natoms = len(xcoords)
"""writes coords onto txt file"""
out = open("sphere.xyz", "a")
print (natoms)
out.write("%s \n \n" %(natoms))
out.close()
for item in zip(xcoords, ycoords,zcoords):
   out = open("sphere.xyz", "a")
   print (item)
   out.write( "Ar"  " " " " "%2.5s %2.5s %2.5s \n" %(item))
   out.close()

I've considered using spherical coordinates to define the parameters for a sphere, but I do not know how to set this up using python.

  • Can you show me an example of 'sphere point' where all coordinates (x, y, z) and r are integers from ranges defined by you? – Łukasz Rogalski Apr 01 '15 at 19:21
  • 1
    This reminds me a bit of an old answer of mine for [Drawing Sphere in OpenGL without using gluSphere()?](http://stackoverflow.com/q/7687148/953482). Maybe it will be useful to you. – Kevin Apr 01 '15 at 19:25
  • I think the biggest problem is your apparent assumption that you will find points on a sphere that have coordinates on a particular grid. Given a centroid on the grid, and a radius that is a multiple of the grid spacing, the vast majority of spheres will only have 8 such points - one in each direction on the radius aligned to each axis through the centroid - the centers of the faces of a grid-aligned bounding cube... – twalberg Apr 01 '15 at 19:35

1 Answers1

1

As I mentioned in comment, I think you are a victim of misunderstanding conditions. My approach would be different - do not iterate over fixed amount of coordinatess and try to hit the spot, but just generate points which are correct. I used common-known spherical to cartesian formulas and just generated fixed samples for both angles. 

from __future__ import division
from math import sin, cos, pi

SAMPLES = 20

def to_x(r, theta, phi):
    return r*sin(theta)*cos(phi)

def to_y(r, theta, phi):
    return r*sin(theta)*sin(phi)

def to_z(r, theta, phi):
    return r*cos(theta)

samples_r = [r for r in xrange(1, 10)]
samples_theta = [theta * pi / SAMPLES for theta in xrange(SAMPLES)]
samples_phi = [phi * 2*pi / SAMPLES for phi in xrange(SAMPLES)]

coords = []

for r in samples_r:
    for theta in samples_theta:
        for phi in samples_phi:
            coords.append((to_x(r, theta, phi), 
                           to_y(r, theta, phi),
                           to_z(r, theta, phi)))

for coord in coords:
   out = open("sphere.xyz", "a")
   print (coord)
   out.write( "Ar"  " " " " "%2.5s %2.5s %2.5s \n" %(coord))
   out.close()

To improve memory performance you may consider using generators and itertools module (especially itertools.product seems to be designed for similar problems).

Łukasz Rogalski
  • 22,092
  • 8
  • 59
  • 93