如何正确地将透视添加到表示分子的 Gnuplot 3D 连接点云？
我不了解你们，但我喜欢一些 Gnuplot。如果使用得当，该软件会生成美丽的图像，其简单性和清晰度令人着迷，我非常喜欢。
没有什么特别的原因，有一天，我发现自己在想，如果我能创作出如此卡通魅力和充满活力的清晰度的图片来进入我的论文和个人科学日记该多好。因此，我首先进入了一个糟糕的项目，编写一个基于 gnuplot 的分子可视化工具。
到目前为止，它是为我的特定类型的分子量身定制的。基本上是共价键合的原子形成配体，配体本身通过配位键与一些中心金属离子相互作用。我已经成功地得出了一个非常好的工作概念，如下图所示。
![[插入用 GNUPLOT 制作的复杂图片]](https://i.stack.imgur.com/FGnOa.png)
其中，虚线表示与铕金属离子的配位键，颜色为浅青色，实线是原子之间的共价键。红色是氧，蓝色是氮，白色是氢，灰色是碳。到目前为止一切都很好，看起来相当可靠并且非常符合我想要的。
那么我该怎么做呢，我听到你问？嗯，实际上这很简单。我一次只策划一件事情。首先，我绘制虚线的连接模式，如下所示：
![[虚线图像]](https://i.stack.imgur.com/nx8w4.png)
然后我画出共价键：
![[共价键]](https://i.stack.imgur.com/qYsMK.png)
每个步骤都需要一个或多个单独的文件。每个配体的连接性存储在单独的“bondfile”中，并且点状连接性模式本身位于一个文件中。原子的位置及其颜色被放置在另一个文件中。一个用于每个配体，一个用于中心金属。
然后我有一个单独的文件，用于记录金属和每个配体的原子，我在其中说明它们的颜色。事实上，原子放置在黑点上，这使得点周围具有迷人的黑色布局，否则它们就没有轮廓线。
问题
当您想要旋转复合体以获得更好的角度来保存到图片中时，就会出现问题。为了说明问题，我将用单个配体的图片来展示它的实际情况。我们以联吡啶为例（带有氮原子的，有两个）
所以这是我认为最佳角度的联吡啶：
![[插入联吡啶静止不动的图片]](https://i.stack.imgur.com/IwJB5.png)
现在假设我们沿着下图所示的轴转动联吡啶。
![[联吡啶与旋转轴的图片]](https://i.stack.imgur.com/tkryy.png)
现在问题就显现出来了。因为一些应该位于平面后面的原子实际上位于整个物体的前面，这表明 gnuplot 实际上没有透视图。或者，至少，它确实有，但我使用不正确。
![[吡啶旋转图片]](https://i.stack.imgur.com/jd0RW.png)
到目前为止，一切都很好。我没想到它会自动透视，因为这不是它最初设计的目的。然而，这意味着 gnuplot “splot” 会进行某种程度虚假的 3D 绘图，并且空间中点的实际相对位置并不重要。
所以我的问题是，对于所有 gnuplot/perspective 专家来说：有没有办法巧妙地规避这个限制？
我对任何方法都感兴趣，无论它有多复杂，只要它在 gnuplot 本身的限制内可行即可。
前段时间，我尝试过类似的事情。显然，点和线不遵守 3D 顺序。然而，如果你用曲面来绘制，即原子=球体和键=圆柱体，它就会起作用。
编辑：这是一个完全修改的版本。我知道有专门的程序可以可视化分子。这只是为了好玩并演示 gnuplot 的可行性。我认为随着原子数量的增加，这个脚本会变得相当慢。
可以直接读取结构数据文件（SDF）文件。它包含原子位置和键信息（键的连接性和类型）。原子显示为球体，键显示为圆柱体。因此，数据块$Sphere和$Cylinders包含球体和圆柱体原型的数据点。有关原子的附加信息存储在数据块中$Elements，即原子序数、元素名称、原子大小和颜色。可以将更多元素添加到此列表中。球体只是根据其位置以偏移量绘制。键还需要适当旋转，这需要键向量的旋转。因此，以下基本向量和矩阵运算被实现为函数：
- 矢量长度(V)
- 叉积(a,b)
- 矩阵向量乘法(M,V)
- 向量归一化(V)
该方法可能不是最有效的方法，因为向量和矩阵被作为字符串处理（具有 3 个和 9 个标记）。
作为说明性示例，咖啡因分子的数据取自此处。
数据： Caffeine.sdf
2519
-OEChem-08062013263D
24 25 0 0 0 0 0 0 0999 V2000
0.4700 2.5688 0.0006 O 0 0 0 0 0 0 0 0 0 0 0 0
-3.1271 -0.4436 -0.0003 O 0 0 0 0 0 0 0 0 0 0 0 0
-0.9686 -1.3125 0.0000 N 0 0 0 0 0 0 0 0 0 0 0 0
2.2182 0.1412 -0.0003 N 0 0 0 0 0 0 0 0 0 0 0 0
-1.3477 1.0797 -0.0001 N 0 0 0 0 0 0 0 0 0 0 0 0
1.4119 -1.9372 0.0002 N 0 0 0 0 0 0 0 0 0 0 0 0
0.8579 0.2592 -0.0008 C 0 0 0 0 0 0 0 0 0 0 0 0
0.3897 -1.0264 -0.0004 C 0 0 0 0 0 0 0 0 0 0 0 0
0.0307 1.4220 -0.0006 C 0 0 0 0 0 0 0 0 0 0 0 0
-1.9061 -0.2495 -0.0004 C 0 0 0 0 0 0 0 0 0 0 0 0
2.5032 -1.1998 0.0003 C 0 0 0 0 0 0 0 0 0 0 0 0
-1.4276 -2.6960 0.0008 C 0 0 0 0 0 0 0 0 0 0 0 0
3.1926 1.2061 0.0003 C 0 0 0 0 0 0 0 0 0 0 0 0
-2.2969 2.1881 0.0007 C 0 0 0 0 0 0 0 0 0 0 0 0
3.5163 -1.5787 0.0008 H 0 0 0 0 0 0 0 0 0 0 0 0
-1.0451 -3.1973 -0.8937 H 0 0 0 0 0 0 0 0 0 0 0 0
-2.5186 -2.7596 0.0011 H 0 0 0 0 0 0 0 0 0 0 0 0
-1.0447 -3.1963 0.8957 H 0 0 0 0 0 0 0 0 0 0 0 0
4.1992 0.7801 0.0002 H 0 0 0 0 0 0 0 0 0 0 0 0
3.0468 1.8092 -0.8992 H 0 0 0 0 0 0 0 0 0 0 0 0
3.0466 1.8083 0.9004 H 0 0 0 0 0 0 0 0 0 0 0 0
-1.8087 3.1651 -0.0003 H 0 0 0 0 0 0 0 0 0 0 0 0
-2.9322 2.1027 0.8881 H 0 0 0 0 0 0 0 0 0 0 0 0
-2.9346 2.1021 -0.8849 H 0 0 0 0 0 0 0 0 0 0 0 0
1 9 2 0 0 0 0
2 10 2 0 0 0 0
3 8 1 0 0 0 0
3 10 1 0 0 0 0
3 12 1 0 0 0 0
4 7 1 0 0 0 0
4 11 1 0 0 0 0
4 13 1 0 0 0 0
5 9 1 0 0 0 0
5 10 1 0 0 0 0
5 14 1 0 0 0 0
6 8 1 0 0 0 0
6 11 2 0 0 0 0
7 8 2 0 0 0 0
7 9 1 0 0 0 0
11 15 1 0 0 0 0
12 16 1 0 0 0 0
12 17 1 0 0 0 0
12 18 1 0 0 0 0
13 19 1 0 0 0 0
13 20 1 0 0 0 0
13 21 1 0 0 0 0
14 22 1 0 0 0 0
14 23 1 0 0 0 0
14 24 1 0 0 0 0
M END
> <PUBCHEM_COMPOUND_CID>
2519
> <PUBCHEM_CONFORMER_RMSD>
0.4
> <PUBCHEM_CONFORMER_DIVERSEORDER>
1
> <PUBCHEM_MMFF94_PARTIAL_CHARGES>
15
1 -0.57
10 0.69
11 0.04
12 0.3
13 0.26
14 0.3
15 0.15
2 -0.57
3 -0.42
4 0.05
5 -0.42
6 -0.57
7 -0.24
8 0.29
9 0.71
> <PUBCHEM_EFFECTIVE_ROTOR_COUNT>
0
> <PUBCHEM_PHARMACOPHORE_FEATURES>
5
1 1 acceptor
1 2 acceptor
3 4 6 11 cation
5 4 6 7 8 11 rings
6 3 5 7 8 9 10 rings
> <PUBCHEM_HEAVY_ATOM_COUNT>
14
> <PUBCHEM_ATOM_DEF_STEREO_COUNT>
0
> <PUBCHEM_ATOM_UDEF_STEREO_COUNT>
0
> <PUBCHEM_BOND_DEF_STEREO_COUNT>
0
> <PUBCHEM_BOND_UDEF_STEREO_COUNT>
0
> <PUBCHEM_ISOTOPIC_ATOM_COUNT>
0
> <PUBCHEM_COMPONENT_COUNT>
1
> <PUBCHEM_CACTVS_TAUTO_COUNT>
1
> <PUBCHEM_CONFORMER_ID>
000009D700000001
> <PUBCHEM_MMFF94_ENERGY>
22.901
> <PUBCHEM_FEATURE_SELFOVERLAP>
25.487
> <PUBCHEM_SHAPE_FINGERPRINT>
10967382 1 18338799025773621285
11132069 177 18339075025094499008
12524768 44 18342463625094026902
13140716 1 17978511158789908153
16945 1 18338517550775811621
193761 8 15816500986559935910
20588541 1 18339082691204868851
21501502 16 18338796715286957384
22802520 49 18128840606503503494
2334 1 18338516344016692929
23402539 116 18270382932679789735
23552423 10 18262240993325675966
23559900 14 18199193898169584358
241688 4 18266458702623303353
2748010 2 18266180539182415717
5084963 1 17698433339235542986
528886 8 18267580380709240570
53812653 166 18198902694142226312
66348 1 18339079396917369615
> <PUBCHEM_SHAPE_MULTIPOLES>
256.45
4.01
2.83
0.58
0.71
0.08
0
-0.48
0
-0.81
0
0.01
0
0
> <PUBCHEM_SHAPE_SELFOVERLAP>
550.88
> <PUBCHEM_SHAPE_VOLUME>
143.9
> <PUBCHEM_COORDINATE_TYPE>
2
5
10
$$$$
代码：
### plot a molecule from an SDF file
reset session
FILE = 'Caffeine.sdf'
DATA = '$Molecule'
# get datafile 1:1 into datablock
if (GPVAL_SYSNAME[:7] eq "Windows") { load '< echo '.DATA.' ^<^<EOD & type "'.FILE.'"' } # Windows
if (GPVAL_SYSNAME eq "Linux") { load '< echo "\'.DATA.' << EOD" & cat "'.FILE.'"' } # Linux
if (GPVAL_SYSNAME eq "Darwin") { load '< echo "\'.DATA.' << EOD" & cat "'.FILE.'"' } # MacOS
AtomCount = word($Molecule[4],1) # get number of atoms in molecule
BondCount = word($Molecule[4],2) # get number of bonds in molecule
# put atom data into a datablock
# X, Y, Z, Element
set print $Atoms
do for [i=5:4+AtomCount] { print $Molecule[i] }
set print
# put bond data into a datablock
# Atom1, Atom2, BondType
set print $Bonds
do for [i=5+AtomCount:4+AtomCount+BondCount] { print $Molecule[i] }
set print
# create sphere datapoints (=atom prototype)
set parametric
set isosamples 17
set samples 17
epsilon=1e-8
set urange [epsilon-pi/2:pi/2+epsilon]
set vrange [0:2*pi]
Radius = 1
set table $Sphere
splot Radius*cos(u)*cos(v), Radius*cos(u)*sin(v), Radius*sin(u)
unset table
# create cylinders (=single, double, triple bond prototype)
set isosamples 2
set samples 12
set urange [-pi:pi]
set vrange [0.2:1]
BondRadius = 0.075
set table $Cylinders # single, double, triple bonds
do for [Offset in "0 -1.25 1.25 -2.5 0 2.5"] {
splot BondRadius*(cos(u)+Offset), BondRadius*sin(u), v
}
unset table
unset parametric
# Lookup table for elements
# AtomicNo ElementSymbol Radius Color
$Elements <<EOD
1 H 1.5 #ffffff
6 C 2.5 #888888
7 N 3.0 #0000ff
8 O 2.5 #ff0000
EOD
# lookup function: search for string s in column c1. If found return value in column c2
LookupElement(s,c1,c2) = (tmp = '', sum [iii=1:|$Elements|] (word($Elements[iii],c1) eq s ? \
(tmp=word($Elements[iii],c2),0) : 0), tmp)
Element(n) = word($Atoms[n],4) # get element of nth atom
ElementNo(n) = int(LookupElement(Element(n),2,1)) # lookup atomic number by nth atom
AtomSize(e) = LookupElement(e,2,3) # lookup atom size by element
AtomSizeScaling = 0.2
AtomPos(n,axis) = word($Atoms[n],axis) # get x=1,y=2,z=3 coordinates of nth atom
AtomPoint(n,axis) = AtomPos(n,axis) + (column(axis)*AtomSize(Element(n))*AtomSizeScaling)
# create atom color palette
AtomPalette = "( -1 '#cccccc'"
do for [i=1:|$Elements|] {
AtomPalette = AtomPalette.sprintf(", %s '%s'",word($Elements[i],1),word($Elements[i],4))
}
AtomPalette = AtomPalette.')'
set palette defined @AtomPalette
# functions for vector and marix operations
VectorLength(V) = sqrt(word(V,1)**2 + word(V,2)**2 + word(V,3)**2)
VectorNormalize(V) = sprintf("%g %g %g", \
word(V,1)/VectorLength(V), word(V,2)/VectorLength(V), word(V,3)/VectorLength(V))
# Cross vector product
CrossProduct(a,b) = sprintf("%g %g %g", \
word(a,2)*word(b,3) - word(a,3)*word(b,2), \
word(a,3)*word(b,1) - word(a,1)*word(b,3), \
word(a,1)*word(b,2) - word(a,2)*word(b,1))
# Rotation matrix: Input vector (normalized) and angle
RotationMatrix(Vn,a) = sprintf("%g %g %g %g %g %g %g %g %g", \
word(Vn,1)*word(Vn,1)*(1-cos(a))+cos(a), \
word(Vn,1)*word(Vn,2)*(1-cos(a))-word(Vn,3)*sin(a), \
word(Vn,1)*word(Vn,3)*(1-cos(a))+word(Vn,2)*sin(a), \
word(Vn,2)*word(Vn,1)*(1-cos(a))+word(Vn,3)*sin(a), \
word(Vn,2)*word(Vn,2)*(1-cos(a))+cos(a), \
word(Vn,2)*word(Vn,3)*(1-cos(a))-word(Vn,1)*sin(a), \
word(Vn,3)*word(Vn,1)*(1-cos(a))-word(Vn,2)*sin(a), \
word(Vn,3)*word(Vn,2)*(1-cos(a))+word(Vn,1)*sin(a), \
word(Vn,3)*word(Vn,3)*(1-cos(a))+cos(a))
# define matrix/vector multiplication (Matrix 3x3, Vector 3x1)
MatrixVectorMultiplication(M,V) = sprintf("%g %g %g", \
word(M,1)*word(V,1) + word(M,2)*word(V,2) + word(M,3)*word(V,3), \
word(M,4)*word(V,1) + word(M,5)*word(V,2) + word(M,6)*word(V,3), \
word(M,7)*word(V,1) + word(M,8)*word(V,2) + word(M,9)*word(V,3))
# Rotation of points
RotatedVector(n) = MatrixVectorMultiplication(RotationMatrix(RotationVector(n),RotationAngle(n)), \
sprintf("%g %g %g", column(1),column(2),column(3)))
# Bond start & end
BondStart(i) = int(word($Bonds[i],1))
BondEnd(i) = int(word($Bonds[i],2))
BondVector(n) = sprintf("%g %g %g", \
AtomPos(BondEnd(n),1) - AtomPos(BondStart(n),1), \
AtomPos(BondEnd(n),2) - AtomPos(BondStart(n),2), \
AtomPos(BondEnd(n),3) - AtomPos(BondStart(n),3))
BondLength(n) = VectorLength(BondVector(n))
BondType(i) = int(word($Bonds[i],3)) # get bond type: single, double, triple
BondTypeStart(n) = BondType(n)==3 ? 3 : BondType(n)==2 ? 1 : 0
BondTypeEnd(n) = BondType(n)==3 ? 5 : BondType(n)==2 ? 2 : 0
# rotation axis vector normalized, (cross-product of BondVector and z-axis)
RotationVector(n) = VectorNormalize(CrossProduct(BondVector(n),"0 0 1"))
# rotation angle (between V and z-axis)
RotationAngle(n) = -acos(word(BondVector(n),3)/VectorLength(BondVector(n)))
BondPoint(n,m) = word(RotatedVector(n),m) + AtomPos(BondStart(n),m)
# plot settings
set cbrange [-1:8]
set view equal xyz
unset border
unset tics
unset colorbox
unset key
set style fill solid 1.0 noborder
set pm3d depthorder noborder
set pm3d lighting specular 0.5
set view 26, 329, 2
splot \
for [i=1:|$Bonds|] $Cylinders u \
(BondPoint(i,1)):(BondPoint(i,2)):(BondPoint(i,3)):(-1) \
index BondTypeStart(i):BondTypeEnd(i) w pm3d, \
for [i=1:|$Atoms|] $Sphere u (AtomPoint(i,1)):(AtomPoint(i,2)):(AtomPoint(i,3)):(ElementNo(i)) w pm3d
### end of code
结果：（Windows7下的wxt终端，gnuplot 5.2.8）
可以使用 来完成动画terminal gif animate，但是，我通过使用 ScreenToGif 软件创建 PNG 并将terminal pngcairo它们组合成动画 gif，获得了更好看的结果。
动画片：
