Why does changing the arguments in my variable have no effect on my transformed coordinates?
I have written this code to wrap a set of 2d coordinates of a hexagonal connected lattice into a 3d cylinder:
import numpy as np
import matplotlib.pyplot as plt
import mpl_toolkits.mplot3d.axes3d as axes3d
import matplotlib.colors
with open("test_A_aux.dat", "r") as aux:
dim = np.genfromtxt(aux, skip_header =3, max_rows = 1)
n = int(input("Enter the value for n: "))
m = int(input("Enter the value for m: "))
# Open coordinate file and import data to array
with open("test_A_crds.dat", "r") as crds:
array = np.genfromtxt(crds)
# Defines the chiral angle based on a1 and a2 vectors
# This is then used to manipulate x,y coords onto the vector
def chiral_angle(n, m):
# a1 vector is set along x axis
a1 = np.array([np.sqrt(3),0])
a2 = np.array([np.sqrt(3)/2, 3/2])
# chiral vector is the dot product of n * a1 and m * a2
chiral_vector = np.add(np.multiply(a1, n), np.multiply(a2, m))
# magnitude of chiral vector
mag = np.linalg.norm(chiral_vector)
# dot product formula to find cos theta
cos_theta = np.dot(chiral_vector, np.array([1,0])) / mag
theta = np.arccos(cos_theta)
return theta
# This function wraps the 2d coordinates into a cylinder by maniuplating
# the y cooridnates and generating z coordinates
# x coordinates unchaged
def wrap(n, m):
# normalises to max length of the dimensions of the input cell
normalise = matplotlib.colors.Normalize(0, dim[1])
X = array[:,0]
THETA = np.array([2*np.pi * normalise(i) for i in array[:,1]])
r = dim[1] / (2 * np.pi)
Y = np.array([r * np.sin(angle) for angle in THETA])
Z = np.array([r * np.cos(angle) for angle in THETA])
# stacks X and Y into n*2 array
tube = np.vstack((X, Y))
# stacks above n*2 array with Z to make n*3 array
tube = np.vstack((tube, Z))
# transposes array to make 3*n array of 3D coordinates
tube = np.transpose(tube)
return tube
# Plotting function
def plot_tube():
fig = plt.figure()
ax = fig.add_subplot(1,1,1, projection='3d')
plot = ax.scatter(tube[:,0], tube[:,1], tube[:,2])
plt.show()
return
# creates file and writes it to vmd for visualisation
# with open("nanotube.xyz", "w") as vmd:
# vmd.write("{:}\n\n " .format(array.shape[0]))
# for i in range(array.shape[0]):
# vmd.write("{:} {:<10} {:<10} {:<10}\n" .format(i, tube[i,0], tube[i,1], tube[i,2]))
tube = wrap(n, m)
print(tube)
plot_tube()
But changing the n and m variables has no affect on the 3d coordinates of the output. For reference the chiral vector is like so:
With a1 set along the x axis.
How can I make n and m affect the chiral vector and hence the angle theta and the resulting shape of the tube?
The lack of usage of chiral_angle in the wrap function explains why you always get the same result when changing n and m parameters.
However, even if you change that, the is a problem with your formulation of the problem. So here is how to have it more clear:
You want to take a set of points in a plane, select a direction and roll the plane into a cylinder. for that you select a chiral angle, with components (m, n). This vector will work as the axis around which you roll the plane, so, for example, in the zigzag configuration, you roll the plane around the x axis, i.e. if you trace a vertical line in the original coordinates, you will get a perfect circle on the cylinder.
The problem now is that you cannot assume that the X axis will be unchanged for any chiral angle, it will only happen if you roll on a (m, 0) vector, which corresponds to zigzag.
Instead, what happens is that the point coordinates along the direction of the chiral vector will become your new X axis, and the orthogonal coordinates will be rolled into a cylinder (more on how to do this next).
So first we calculate the chiral vector (m, n), returning the normalized vector will make things easier.
Then we compute the point projections of all points to this chiral vector (this will return the "part" of each point that is along the chiral vector.
then we compute the orthogonals (find the orthognal vector to the chiral, if chiral=(x, y), the orthogonal will be (-y, x) to rotate 90° counter-clockwise) (this will return the "part" of each point that is perpendicular to the chiral vector.
Your X coordinates will be the the projections (how far from the origin is this point along the chiral direction)
and for the Y and Z we first need to convert the orthogonals into a corresponding angle in the roll. The orthogonals are how far is it from the chiral direction. Then normalize them so that the most distant points from the chiral vector convert to 2*pi, I use the factor 2pi/(max-min).
Now you can compute Y as the sine of this "angle" and Z as the cossine, and thats how you get all the coordinates.
Bellow is the code I adapted from you:
import numpy as np
import matplotlib.pyplot as plt
import mpl_toolkits.mplot3d.axes3d as axes3d
import matplotlib.colors
# with open("test_A_aux.dat", "r") as aux:
# dim = np.genfromtxt(aux, skip_header =3, max_rows = 1)
dim = np.array([69.282032]*2)
n = int(input("Enter the value for n: "))
m = int(input("Enter the value for m: "))
# Open coordinate file and import data to array
with open("test_A_crds.dat", "r") as crds:
array = np.genfromtxt(crds)
# Defines the chiral vector based on a1 and a2 vectors
# This is then used to manipulate x,y coords onto the vector
def chiral_vector(n, m):
# a1 vector is set along x axis
a1 = np.array([np.sqrt(3),0])
a2 = np.array([np.sqrt(3)/2, 3/2])
# chiral vector is the dot product of n * a1 and m * a2
chiral = np.add(np.multiply(a1, n), np.multiply(a2, m))
return chiral/np.linalg.norm(chiral)
# compute projections (direction* (dot product between direction and vector))
def project_point(point, chiral):
projection = (point.T@chiral)
return projection
# This function wraps the 2d coordinates into a cylinder by maniuplating
# the y cooridnates and generating z coordinates
# x coordinates unchaged
def wrap(n, m):
# normalises to max length of the dimensions of the input cell
cv = chiral_vector(n, m)
# make projections
X = np.array([project_point(x, cv) for x in array])
# make orthogonals
ortho_cv = np.array([-cv[1], cv[0]])
Y_ = np.array([project_point(x, ortho_cv) for x in array])
# normaliye those distances so that the furthest point gets converted to 2*pi
y_length_factor = (2*np.pi)/(np.max(Y_)-np.min(Y_))
# get the corresponding roll angle
THETA = Y_*y_length_factor
# there is probably a way to optimize this r to maintain proportions in the final plot
r = 1
# numpy applies these functions to each element of the input array
Y = r * np.sin(THETA)
Z = r * np.cos(THETA)
# you can stack all of them together
tube = np.vstack((X, Y, Z))
# transposes array to make 3*n array of 3D coordinates
tube = np.transpose(tube)
return tube
# Plotting function
def plot_tube():
fig = plt.figure()
ax = fig.add_subplot(1,1,1, projection='3d')
plot = ax.scatter(tube[:,0], tube[:,1], tube[:,2])
plt.show()
return
# creates file and writes it to vmd for visualisation
# with open("nanotube.xyz", "w") as vmd:
# vmd.write("{:}\n\n " .format(array.shape[0]))
# for i in range(array.shape[0]):
# vmd.write("{:} {:<10} {:<10} {:<10}\n" .format(i, tube[i,0], tube[i,1], tube[i,2]))
tube = wrap(n, m)
print(tube)
plot_tube()