1.问题背景

  在一个半径为1的球面上，均匀的放置N个点，使得两两之间的最近距离最大，也就是说使得其中任意一点到离它最近的点的距离最大。比如说有两个点，那么放在两极的位置上。
  这样可以用来生成类似C-60的结构，用于构建大型的碳球。

2.构建代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from matplotlib.tri import Triangulation

def generate_uniform_points_on_sphere(num_points):
    points = []
    phi = np.arccos(1 - 2 * np.linspace(0, 1, num_points))
    #np.linspace(0, 1, num_points) 生成一个包含 num_points 个元素的线性间隔数组，范围在 [0, 1] 之间。
    #1 - 2 * np.linspace(0, 1, num_points) 将上述数组映射到 [-1, 1] 的范围。
    #np.arccos(1 - 2 * np.linspace(0, 1, num_points)) 计算反余弦，得到一组角度值 phi。这些角度值用于表示球坐标系中的纬度
    theta = np.pi * (3. - np.sqrt(5.)) * np.arange(num_points)
    #np.arange(num_points) 生成一个包含 num_points 个元素的数组，表示经度。
    #np.pi * (3. - np.sqrt(5.)) 是常数，用于调整生成的角度值。
    #np.pi * (3. - np.sqrt(5.)) * np.arange(num_points) 得到一组角度值 theta。这些角度值用于表示球坐标系中的经度。
    x, y, z = np.cos(theta) * np.sin(phi), np.sin(theta) * np.sin(phi), np.cos(phi)
    # np.cos(theta) * np.sin(phi)、np.sin(theta) * np.sin(phi) 和 np.cos(phi) 分别计算球坐标系中的 x、y 和 z 分量。
    #这些计算转换球坐标系中的角度值（phi 和 theta）为笛卡尔坐标系中的坐标。
    points = np.column_stack((x, y, z))

    return points

def select_uniform_points(base_points, num_selected_points):
    num_base_points = len(base_points)

    # 从已有的点中选择均匀分布的新点
    selected_indices = np.linspace(0, num_base_points - 1, num_selected_points, dtype=int)
    #selected_indices2 = np.linspace(0, num_base_points - 1, num_selected_points, dtype=int)
    #print(selected_indices2)
    selected_points = base_points[selected_indices]

    return selected_points


num_points = 120
base_points = generate_uniform_points_on_sphere(num_points)

num_selected_points = 60
selected_points = select_uniform_points(base_points, num_selected_points)

jg=open("xyz.txt",'w')

# 输出所有点的坐标
print("All points:")
for i, point in enumerate(base_points):
    #print(f"Point {i + 1}: {point}")
    jg.writelines(str(point[0]) + ' ' + str(point[1]) + ' ' + str(point[2]) + '\n')

"""for i, point in enumerate(selected_points): #生成O
    #print(f"Point {i + 1}: {point}")
    jg.writelines(str(point[0]) + ' ' + str(point[1]) + ' ' + str(point[2]) + '\n')
for i, point in enumerate(selected_points): #生成H
    jg.writelines(str(point[0]) + ' ' + str(point[1]) + ' ' + str(point[2]) + '\n')"""

3.转变成VMD可读取的.xyz文件

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
import math
import numpy as np

numC = 120 #C的数目
numO = 0 #O的数目
xyz = open("xyz.txt",'r') #之前生成的文件存放的位置
result = open("result.xyz",'w') 

rC_O = 1.35 #C和O之间的距离，单位为埃
rO_H = 0.97 #O和H之间的距离，单位为埃

numH = numO
d60 = 7.068 #标准C60的半径，单位为埃
rC = d60 /2 * math.sqrt(numC/60)
rO = rC + rC_O
rH = rO + rO_H
num = numC+numO+numH
result.writelines(str(num)+'\nPDB File\n')

#计算H原子坐标
def rotate_point(x, y, z, angle_degrees, distance):
    # 将角度转换为弧度
    angle_radians = np.radians(angle_degrees)

    # 计算原点到给定点的单位向量
    magnitude = np.sqrt(x**2 + y**2 + z**2)
    unit_vector = np.array([x, y, z]) / magnitude

    # 计算旋转后的坐标
    new_x = x + distance * np.cos(angle_radians)
    new_y = y + distance * np.sin(angle_radians)
    new_z = z

    return new_x, new_y, new_z


xh = 0
for line in xyz:
    line = line.strip('\n')
    data = line.split()
    s3=0
    s4=0
    s5=0

    if(xh<numC):
        s3 = str('%.6f' % (float(data[0]) * rC))
        s4 = str('%.6f' % (float(data[1]) * rC))
        s5 = str('%.6f' % (float(data[2]) * rC))
        result.writelines(" C")
    elif(xh<numO+numC):
        s3 = str('%.6f' % (float(data[0]) * rO))
        s4 = str('%.6f' % (float(data[1]) * rO))
        s5 = str('%.6f' % (float(data[2]) * rO))
        result.writelines(" O")
    elif(xh<numH+numO+numC):
        s3 = str('%.6f' % (float(data[0]) * rH))
        s4 = str('%.6f' % (float(data[1]) * rH))
        s5 = str('%.6f' % (float(data[2]) * rH))
        """a, b, c = rotate_point(float(data[0])*rO, float(data[1])*rO, float(data[2])*rO,109,rO_H)
        s3 = str('%.6f' % (a))
        s4 = str('%.6f' % (b))
        s5 = str('%.6f' % (c))"""
        result.writelines(" H")
    xh = xh + 1
    fs3 = f"{s3:>9}"
    fs4 = f"{s4:>9}"
    fs5 = f"{s5:>9}"

    result.writelines("   "+fs3+"   "+fs4+"   "+fs5+"\n")
  如果需要在C球外生成连着的OH或者其他原子，可以改变第一个代码的最后和第二个代码的numO等数据。