0

我正在尝试将 (x,y,z) 坐标包裹在莫比乌斯带上（通过扭曲一次并连接带的末端而获得的拓扑结构）。

结构（我可以为其创建 xyz 坐标）如下（来自RG）使用以下代码。

import numpy as np
import matplotlib.pyplot as plt

bLength=1.6
numPoints=10
radius = bLength*numPoints / (2 * np.pi)
theta = np.linspace(0,2*np.pi,numPoints,endpoint=False)
dtheta=theta[1]-theta[0]
x0,y0=(radius * np.cos(theta)), (radius * np.sin(theta))
x1,y1=(radius * np.cos(theta+dtheta/2)) , (radius * 
np.sin(theta+dtheta/2))
#plt.plot(x0,y0)
#plt.show()
#plt.plot(x1,y1)
#plt.show()
cons0=np.ones(x0.shape)*0
cons1=np.ones(x1.shape)*2
np.savetxt('cooRing1.csv',np.c_[x0,y0,cons0],delimiter=' ')
np.savetxt('cooRing2.csv',np.c_[x1,y1,cons1],delimiter=' ')
cat cooRing1.csv cooRing2.csv > coordinates.csv

[![在此处输入图像描述][3]][3]

我想在莫比乌斯带上映射 xyz 坐标。

以下示例已在网站中给出了不同的条带（代码也使用了一个不公开的模块）

import Twister as Twister    
import math

def displacement(x, width, wrapping_angle):
    """
       Function for converting a nanosheet coordinate into a partly wrapped nanotube
       @param x : Coordinates of nanosheet atom
       @param width : Width of the nano-sheet
       @param wrapping_angle : maximum wrapping angle of the nanotube in radians
    """
    # calculate the average radius of the incomplete wrapped tube
    radius = width/wrapping_angle
    # find the angle of the current atom
    angle = (x[2]-width/2.)/radius
    # calculate the radius of the current atom
    atom_radius = radius+x[1]

    # return atom position of the wrapped atom
    return numpy.array([x[0], atom_radius*math.cos(angle),atom_radius*math.sin(angle)])

def configuration(n, m, repetition):
    """
       Function for generating a moebius molecule
       @param n : Chiral vector index
       @param m : Chiral vector index
       @param repetition : Repetition along z
    """

    # build n,m ribbon
    ribbon = NanoRibbon(n,m)
    ribbon = ribbon.repeat(1,1,repetition)

    # get properties of the ribbon
    lattice = ribbon.bravaisLattice()
    elements = ribbon.elements()
    cartesian_coordinates=ribbon.cartesianCoordinates().inUnitsOf(Angstrom)

    # calculate the length of the 1-d structure
    z_length = numpy.linalg.norm(lattice.primitiveVectors()[2].inUnitsOf(Angstrom)) 

    # calculate twist parameters
    rotation_angle_per_z =  math.pi /z_length
    rotation_axis = numpy.array([0,0,1])
    rotation_axis_center = numpy.sum(cartesian_coordinates,axis=0)/len(cartesian_coordinates)

    # define a function of one variable, f(c),  for displacing the atoms
    f = lambda c : Twister.displacement(c, rotation_angle_per_z, rotation_axis, 
                                rotation_axis_center, 0.,z_length)
    # apply the function to find new displaced coordinates
    cartesian_coordinates = numpy.apply_along_axis(f, 1, cartesian_coordinates)
    cartesian_center = numpy.sum(cartesian_coordinates,axis=0)/len(cartesian_coordinates)
    cartesian_coordinates = cartesian_coordinates - cartesian_center


    # define a function of one variable, f(c),  for displacing the atoms
    f = lambda c : displacement(c, z_length,2.0*math.pi) 
    # apply the function to find new displaced coordinates
    cartesian_coordinates = numpy.apply_along_axis(f, 1, cartesian_coordinates)

    return MoleculeConfiguration(
            elements=elements,
            cartesian_coordinates=cartesian_coordinates * Angstrom
            )

# Instantiate the builder object and choose our title
builder = Builder()
builder.title('Moebius ribbon')

# Set the configuration generator
builder.setConfigurationGenerator(configuration)

# Tube properties group
builder.newGroup('Tube parameters')
builder.integer( 'n', 6, 'n', min=1, max=1000)
builder.integer( 'm', 0, 'm', min=0, max=1000)
builder.integer( 'repetition', 40, 'C-direction', min=1, max=1000)

Python中是否有任何类似的模块，以便我可以构建我想要的结构并创建xyz坐标？或者我如何才能进一步完成代码？

4

1 回答 1

0

无法将您的 2xN 环之一的坐标转换为莫比乌斯带。那是因为一个带有半扭曲的交替带必须有奇数个原子，而你当前的环总是有一个偶数。您需要切掉一个原子或添加一个原子，以使扭曲起作用。也许您可以通过将第一个和最后一个原子彼此堆叠来进行转换，但我认为它会非常难看（无论是数学上还是在绘图中）。

虽然这可能使您对 3D 变换的要求变得不切实际，但您可以改为从头开始创建具有所需扭曲的坐标。只需为带的单个边缘生成一个点数组（只有一个边缘是使莫比乌斯带奇怪的原因之一）。边缘在圆周围形成两个环，在条带的每一侧都有一个，扭曲发生在主旋转速度的一半。

这是我的尝试：

bLength=1.6
n = 10
numPoints = 2*n + 1
radius = bLength*numPoints / (4*np.pi)   # 4*pi because the edge goes around the circle twice
theta = np.linspace(0, 4*np.pi, numPoints, endpoint=False)   # here too
x = (radius + np.cos(theta/2)) * np.cos(theta)
y = (radius + np.cos(theta/2)) * np.sin(theta)
z = np.sin(theta/2)

plt.plot(x, y) # I don't know how to graph in 3d, so only x and y here
plt.show()

绘图的输出

这使得一条宽度为 2 的条带，就像您当前的代码一样。如果这不是正确的宽度，您应该将计算中添加到半径的项乘以所需cos宽度x的一半。ysinz

请注意，虽然此代码可能会产生您想要的内容，但在物理描述实际原子环时，坐标可能是无稽之谈。坐标之间的一些距离将与其他距离非常不同（例如，条带的内侧与几乎平坦时的外侧），并且由原子制成的真正的莫比乌斯条带（如果可以制造这样的东西的话）会可能以某种方式折叠而不是以均匀的速度扭曲（纸制成的莫比乌斯带也不会以这种方式扭曲）。找出真正的原子是如何排列的将是一个更加困难的问题（对于物理学家而不是程序员来说也是一个问题）。

于 2018-09-05T22:30:14.603 回答