I am working on an illustration for my PhD thesis that will show all interactions commonly included in the molecular quantum Hamiltonian operator.
My idea is as follows:
- Add all nuclei as ico spheres
- Add electrons around the nuclei as tiny spheres
- Add all nuclei-nuclei (NN) interactions as cylinders
- Add all electron-electron (ee) interactions as cylinders
- Add all nuclei-electron (Ne) interactions as cylinders
I separate steps 3-5 so that I easily can keep track of the nature of the interactions and color them differently.
The molecule I want to show is a 54-atom transition metal complex, with 234 electrons.
The number of unique pair interactions $N$ for a set of $n$ particles is
$$ N = \frac{n(n-1)}{2} $$
This amounts to 1431 unique NN interactions and 27261 ee interactions, and the number Ne interactions which I have not calculated.
The first couple of thousand cylinders are relatively efficient, then it becomes much slower. At around 10000 the script adds about 1 cylinder per second, and so I stopped it there. However, I would like to include all the interactions.
Is there a way to achieve this in an efficient way? Could I use an array modifier, and still be able to orient the cylinders correctly? If it is not feasible, then I can always choose a smaller molecule.
Here is the script I use:
import bpy
import numpy as np
import math
from itertools import combinations
settings = {
'H': 1,
'C': 6,
'N': 7,
'O': 8,
'Ni': 28
}
def fiboSphere(n=100, origo=(0, 0, 0), r=1):
goldenRatio = (1 + 5**0.5)/2
grid = np.arange(0, n)
theta = np.pi * grid / goldenRatio
phi = np.arccos(1 - 2 * (grid) / n)
x = r * np.cos(theta) * np.sin(phi) + origo[0]
y = r * np.sin(theta) * np.sin(phi) + origo[1]
z = r * np.cos(phi) + origo[2]
return zip(x, y, z)
def readCoords(f):
with open(f) as file:
lines = file.readlines()
labels = [line.split()[0] for line in lines[2:]]
coords = np.asarray([list(map(float, line.split()[1:])) for line in lines[2:]])
masses = [settings[lab] for lab in labels]
return labels, coords - np.average(coords, axis=0, weights=masses)
def addSphere(x, y, z, r=1):
bpy.ops.mesh.primitive_uv_sphere_add(radius=r, location=(x, y, z))
#bpy.ops.outliner.item_rename('Nuclei')
def addIcoSphere(x, y, z, r=1):
bpy.ops.mesh.primitive_ico_sphere_add(subdivisions=1, radius=r, location=(x, y, z))
def cylinder_between(x1, y1, z1, x2, y2, z2, r):
dx = x2 - x1
dy = y2 - y1
dz = z2 - z1
dist = math.sqrt(dx2 + dy2 + dz**2)
ob = bpy.ops.mesh.primitive_cylinder_add(
radius = r,
depth = dist,
location = (dx/2 + x1, dy/2 + y1, dz/2 + z1)
)
phi = math.atan2(dy, dx)
theta = math.acos(dz/dist)
bpy.context.object.rotation_euler[1] = theta
bpy.context.object.rotation_euler[2] = phi
if name == 'main':
SPHERES = 0
LINE_NUCS = 0
LINE_ELECS = 1
molecule = 'Ni-NHC2.xyz'
r_elec = 0.05
r_nuc = r_elec * 3
r_atom = r_nuc * 3
r_cyl = 0.0001
atoms, coords = readCoords(molecule)
if SPHERES == 1:
for a, c in zip(atoms, coords):
if a == 'Ni':
addIcoSphere(*c, r=r_nuc * np.sqrt(settings[a]/4))
else:
addIcoSphere(*c, r=r_nuc * np.sqrt(settings[a]/2))
for elec in fiboSphere(n=settings[a], origo=c, r=r_atom):
addSphere(*elec, r=r_elec)
if LINE_NUCS == 1:
natoms = len(coords)
combos = list(combinations(np.arange(natoms), 2))
for i, combo in enumerate(combos):
print(i, len(combos))
i, j = combo
xA, yA, zA = coords[i]
xB, yB, zB = coords[j]
cylinder_between(xA, yA, zA, xB, yB, zB, r_cyl)
if LINE_ELECS == 1:
nelecs = 0
for a in atoms:
nelecs += settings[a]
combos = list(combinations(np.arange(nelecs), 2))
all_elecs = []
for a, c in zip(atoms, coords):
for elec in fiboSphere(n=settings[a], origo=c, r=r_atom):
all_elecs.append(elec)
curr = 3000
for i, combo in enumerate(combos[curr:curr+500]):
print(i+curr, len(combos))
i, j = combo
xA, yA, zA = all_elecs[i]
xB, yB, zB = all_elecs[j]
cylinder_between(xA, yA, zA, xB, yB, zB, r_cyl)
and here are the coordinates for the molecule:
54
Ni 0.53508 2.94170 -1.19095
C 0.89172 2.78025 0.55921
C 0.13343 1.34854 -1.93501
C -0.80330 4.11901 -1.38999
O -0.19137 0.31584 -2.35779
O 1.09923 2.70358 1.70102
O -1.66995 4.88913 -1.48201
H 0.51175 2.95147 -6.22064
C 1.10271 3.13388 -5.30374
H 0.84068 0.57894 -6.28354
H 1.83440 3.94105 -5.51296
H 0.41920 3.51974 -4.51700
H 0.02746 7.62618 -2.63086
H 1.69364 7.30829 -3.21356
C 0.88136 6.92729 -2.55858
C 1.51452 0.62494 -5.41305
H 1.50838 -1.80537 -6.60596
H 0.56578 5.95106 -2.97683
C 1.77035 1.87351 -4.82202
C 1.73376 -1.89724 -5.52533
H 0.13563 8.40984 -0.37790
C 2.07482 -0.56412 -4.90989
C 1.34936 6.78856 -1.13348
H 0.83330 -2.32912 -5.03778
C 2.62737 1.90396 -3.70207
C 0.87912 7.63910 -0.11886
H 2.55568 -2.62944 -5.40126
C 3.95328 4.00413 -3.41477
N 2.89355 3.16567 -3.07152
H 4.66590 3.74492 -4.20292
C 2.12853 3.69055 -2.05863
C 3.86116 5.09719 -2.59552
N 2.74966 4.88583 -1.78318
H 4.47675 5.99767 -2.51503
C 2.28833 5.79755 -0.77436
C 2.91681 -0.48577 -3.78475
H 0.54777 9.43656 1.90824
C 3.20418 0.73782 -3.15710
C 1.32300 7.52594 1.21247
C 0.75840 8.41761 2.28950
H -0.20145 8.01023 2.67381
H 3.35077 -1.40920 -3.36894
C 2.78466 5.66588 0.53957
C 2.28179 6.54227 1.51721
C 4.03233 0.81322 -1.90217
H 3.39997 1.15882 -1.05618
H 1.44677 8.50453 3.15264
H 4.86553 1.53987 -1.99398
C 3.82449 4.62798 0.87242
H 4.80908 4.88876 0.42895
H 4.45676 -0.17386 -1.63994
H 3.54684 3.63357 0.47078
H 2.64748 6.44572 2.55195
H 3.95682 4.53377 1.96626
And a rendered image of some of the NN interactions:
Update
I was able to use bmesh to add the cylinders, and to align them correctly, in this way:
(with credit to this answer for the original idea)
def cylinder_between(x1, y1, z1, x2, y2, z2, r, id=''):
dx = x2 - x1
dy = y2 - y1
dz = z2 - z1
dist = math.sqrt(dx**2 + dy**2 + dz**2)
loc = np.asarray((x1 + dx/2, y1 + dy/2, z1 + dz/2))
name = f'Cylinder_{id}'
mesh = bpy.data.meshes.new(name)
obj = bpy.data.objects.new(name, mesh)
bpy.context.collection.objects.link(obj)
bm = bmesh.new()
cone = bmesh.ops.create_cone(bm,
cap_ends=True,
cap_tris=False,
segments=3,
diameter1=r_cyl,
diameter2=r_cyl,
depth=dist)
obj.location.x = loc[0]
obj.location.y = loc[1]
obj.location.z = loc[2]
phi = math.atan2(dy, dx)
theta = math.acos(dz/dist)
obj.rotation_euler[1] = theta
obj.rotation_euler[2] = phi
bm.to_mesh(mesh)
bm.free()
bpy.ops.meshcalls. I believe operators induce a scene update, which walks over all scene objects, and as a result doing soNtimes inducesN^2run-time problems – NeverConvex Oct 20 '21 at 19:42Atomic blenderis able to add tens of thousands of atoms in a few seconds, with bonds and everything. So it should be doable, I feel. I just don't know enough about the workflow and the API to know where to start. Perhaps my bonds don't need to be cylinders, but curves with a thickness.. – Yoda Oct 20 '21 at 19:49