如何正确地向代表分子的 Gnuplot 3D 连接点云添加透视图？
How to Correctly Add Perspective to a Gnuplot 3D connected Point Cloud Representing a Molecule?
我不了解你们，但我喜欢一些 Gnuplot。如果使用得当，该软件可以生成精美的图像，其简洁明了的魅力，我非常喜欢。
无缘无故，有一天我发现自己在想，如果我能创作出如此具有卡通魅力和充满活力的清晰图片，并将其用于我的论文和个人科学期刊中，那该有多好。所以首先进入了一个 batshit 项目，编写一个基于 gnuplot 的分子可视化器。
到目前为止，它是为我的特定分子类型量身定做的。基本上是形成配体的共价键合原子，它们本身通过配位键与一些中心金属离子相互作用。我已经设法得出一个非常好的工作概念，如下图所示。
其中，虚线表示与铕金属离子的配位键，浅青色，实线为原子间的共价键。红色是氧，蓝色是氮，白色是氢，灰色是碳。到目前为止一切顺利，看起来很扎实，非常符合我的要求。
那么我该怎么做，我听到你在问？其实这很简单。我一次一个地策划事情。首先，我绘制虚线的连接模式，如下所示：
然后我在共价键上作画：
每个步骤都需要一个或多个单独的文件。每个配体的连通性都存储在一个单独的“键合文件”中，点状连接模式本身就在一个文件中。具有颜色的原子的位置被放置在另一个文件中。每个配体一个，中心金属一个。
然后我有一个单独的文件，用于金属和每个配体的原子，我在其中说明它们是什么颜色。原子被放置在黑点上的事实是点周围迷人的黑色布局，否则它们没有轮廓线。
问题
当您想旋转复合体以获得更好的角度以保存到图片中时，就会出现问题。为了说明这个问题，我将用单个配体的图片来展示它的实际效果。让我们拿联吡啶（带氮的那个，有两个）
所以这是我认为最佳角度的联吡啶：
现在假设我们沿着下图所示的轴转动联吡啶。
现在问题出现了。因为一些本应位于后面平面的原子实际上位于整个平面的前面，这表明 gnuplot 实际上没有透视图。或者，至少，它确实有，但我使用不正确。
到目前为止一切顺利。我没想到它会自动拥有透视图，因为这不是它最初制作的目的。然而，这意味着 gnuplot“splot”做了一个有点假的 3D 绘图，并且 space 中点的实际相对位置无关紧要。
所以我的问题是，对于你们所有的 gnuplot/perspective 学者：有没有办法巧妙地绕过这个限制？
我对任何方法都感兴趣，但只要在gnuplot本身的限制范围内可行，就可以。
呵呵。我自己就是一个分子图形专家，从 1970 年代的研究生时代开始就编写了查看器和可视化工具。你知道吗？我真的不喜欢在分子图形中使用透视法。如此之多，以至于我将 gnuplot 中的缺失称为功能而不是限制。
gnuplot 集合中有显示简单分子图形的演示 molecule.dem。它在 gnuplot (5.3) 的开发版本中工作得更好，您可以在其中为原子使用绘图样式 "with circles" 而不是 "with points"。给你：
set title "GM1 pentasaccharide ball-and-stick representation"
set hidden3d
set border 0
unset tics
unset key
set title offset 0, screen -0.85
set view equal xyz
set view 348, 163, 1.64872, 1.14
set style fill transparent solid 0.9 border -1
atomcolor(name) = name[1:1] eq "O" ? 0xdd2222 : name [1:1] eq "N" ? 0x4444ff : 0x888888
splot 'GM1_sugar.pdb' using 6:7:8:(0.6):(atomcolor(strcol(3))) with circles fc rgb var, \
'GM1_bonds.r3d' using 1:2:3:(-):(-):(-) with vectors nohead lw 3 lc "black"
备注：
- 原子位置直接从 PDB 文件中读取
- 原子颜色根据原子名称生成（灰色代表碳，蓝色代表氮等）
- 这些键是通过 Raster3D 分子图形包中的 "bonds" 实用程序从同一个 PDB 文件生成的
- 后面的原子被前面的原子遮挡由"set hidden3d"
处理
- 键的封闭不太令人满意，因为线段一直绘制到原子中心，而在视觉上终止于原子的投影球面看起来会更好。
- 通过使原子部分透明，进一步有助于深度的视觉印象。
前段时间，我尝试过类似的东西。显然，点和线不遵守 3D 顺序。但是，如果您使用曲面进行绘制，它将起作用，即 atoms=spheres 和 bonds=cylinders。
编辑： 这是一个完全修订的版本。
我知道有专门的程序可以可视化分子。这只是为了好玩，并展示了 gnuplot 的可行性。
我假设随着原子数量的增加，这个脚本会变得很慢。
可以直接读取结构数据文件（SDF）文件。它包含原子位置和键信息（连接性和键类型）。
原子显示为球体，键显示为圆柱体。因此，数据块 $Sphere 和 $Cylinders 包含球体和圆柱体原型的数据点。
有关原子的其他信息存储在数据块 $Elements 中，即原子序数、元素名称、原子大小和颜色。可以将更多元素添加到此列表中。
球体只是根据它们的位置用偏移量绘制的。
键也需要适当旋转，这需要键向量的旋转。
因此，将以下基本向量和矩阵运算实现为函数：
- 向量长度(V)
- 叉积(a,b)
- 矩阵向量乘法(M,V)
- 向量标准化(V)
该方法可能不是最有效的方法，因为向量和矩阵被作为字符串处理（具有 3 和 9 个标记）。
举个例子，咖啡因分子的数据取自here。
数据：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 来完成，但是，我通过使用 terminal pngcairo 创建 PNG 并使用软件 ScreenToGif 将它们组合成动画 gif 获得了更好看的结果。
动画：
我不了解你们，但我喜欢一些 Gnuplot。如果使用得当，该软件可以生成精美的图像，其简洁明了的魅力，我非常喜欢。
无缘无故，有一天我发现自己在想，如果我能创作出如此具有卡通魅力和充满活力的清晰图片，并将其用于我的论文和个人科学期刊中，那该有多好。所以首先进入了一个 batshit 项目，编写一个基于 gnuplot 的分子可视化器。
到目前为止，它是为我的特定分子类型量身定做的。基本上是形成配体的共价键合原子，它们本身通过配位键与一些中心金属离子相互作用。我已经设法得出一个非常好的工作概念，如下图所示。
其中，虚线表示与铕金属离子的配位键，浅青色，实线为原子间的共价键。红色是氧，蓝色是氮，白色是氢，灰色是碳。到目前为止一切顺利，看起来很扎实，非常符合我的要求。
那么我该怎么做，我听到你在问？其实这很简单。我一次一个地策划事情。首先，我绘制虚线的连接模式，如下所示：
然后我在共价键上作画：
每个步骤都需要一个或多个单独的文件。每个配体的连通性都存储在一个单独的“键合文件”中，点状连接模式本身就在一个文件中。具有颜色的原子的位置被放置在另一个文件中。每个配体一个，中心金属一个。
然后我有一个单独的文件，用于金属和每个配体的原子，我在其中说明它们是什么颜色。原子被放置在黑点上的事实是点周围迷人的黑色布局，否则它们没有轮廓线。
问题
当您想旋转复合体以获得更好的角度以保存到图片中时，就会出现问题。为了说明这个问题，我将用单个配体的图片来展示它的实际效果。让我们拿联吡啶（带氮的那个，有两个）
所以这是我认为最佳角度的联吡啶：
现在假设我们沿着下图所示的轴转动联吡啶。
现在问题出现了。因为一些本应位于后面平面的原子实际上位于整个平面的前面，这表明 gnuplot 实际上没有透视图。或者，至少，它确实有，但我使用不正确。
到目前为止一切顺利。我没想到它会自动拥有透视图，因为这不是它最初制作的目的。然而，这意味着 gnuplot“splot”做了一个有点假的 3D 绘图，并且 space 中点的实际相对位置无关紧要。
所以我的问题是，对于你们所有的 gnuplot/perspective 学者：有没有办法巧妙地绕过这个限制？
我对任何方法都感兴趣，但只要在gnuplot本身的限制范围内可行，就可以。
呵呵。我自己就是一个分子图形专家，从 1970 年代的研究生时代开始就编写了查看器和可视化工具。你知道吗？我真的不喜欢在分子图形中使用透视法。如此之多，以至于我将 gnuplot 中的缺失称为功能而不是限制。
gnuplot 集合中有显示简单分子图形的演示 molecule.dem。它在 gnuplot (5.3) 的开发版本中工作得更好，您可以在其中为原子使用绘图样式 "with circles" 而不是 "with points"。给你：
set title "GM1 pentasaccharide ball-and-stick representation"
set hidden3d
set border 0
unset tics
unset key
set title offset 0, screen -0.85
set view equal xyz
set view 348, 163, 1.64872, 1.14
set style fill transparent solid 0.9 border -1
atomcolor(name) = name[1:1] eq "O" ? 0xdd2222 : name [1:1] eq "N" ? 0x4444ff : 0x888888
splot 'GM1_sugar.pdb' using 6:7:8:(0.6):(atomcolor(strcol(3))) with circles fc rgb var, \
'GM1_bonds.r3d' using 1:2:3:(-):(-):(-) with vectors nohead lw 3 lc "black"
备注：
- 原子位置直接从 PDB 文件中读取
- 原子颜色根据原子名称生成（灰色代表碳，蓝色代表氮等）
- 这些键是通过 Raster3D 分子图形包中的 "bonds" 实用程序从同一个 PDB 文件生成的
- 后面的原子被前面的原子遮挡由"set hidden3d" 处理
- 键的封闭不太令人满意，因为线段一直绘制到原子中心，而在视觉上终止于原子的投影球面看起来会更好。
- 通过使原子部分透明，进一步有助于深度的视觉印象。
前段时间，我尝试过类似的东西。显然，点和线不遵守 3D 顺序。但是，如果您使用曲面进行绘制，它将起作用，即 atoms=spheres 和 bonds=cylinders。
编辑： 这是一个完全修订的版本。 我知道有专门的程序可以可视化分子。这只是为了好玩，并展示了 gnuplot 的可行性。 我假设随着原子数量的增加，这个脚本会变得很慢。
可以直接读取结构数据文件（SDF）文件。它包含原子位置和键信息（连接性和键类型）。
原子显示为球体，键显示为圆柱体。因此，数据块 $Sphere 和 $Cylinders 包含球体和圆柱体原型的数据点。
有关原子的其他信息存储在数据块 $Elements 中，即原子序数、元素名称、原子大小和颜色。可以将更多元素添加到此列表中。
球体只是根据它们的位置用偏移量绘制的。
键也需要适当旋转，这需要键向量的旋转。
因此，将以下基本向量和矩阵运算实现为函数：
- 向量长度(V)
- 叉积(a,b)
- 矩阵向量乘法(M,V)
- 向量标准化(V)
该方法可能不是最有效的方法，因为向量和矩阵被作为字符串处理（具有 3 和 9 个标记）。
举个例子，咖啡因分子的数据取自here。
数据：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 来完成，但是，我通过使用 terminal pngcairo 创建 PNG 并使用软件 ScreenToGif 将它们组合成动画 gif 获得了更好看的结果。
动画：