gnuplot - 如何正确地将透视图添加到代表分子的 Gnuplot 3D 连接点云？

Question

我不了解你们，但我爱我一些 Gnuplot。如果使用得当，该软件会产生漂亮的图像，它们的简单和清晰令人着迷，我非常喜欢。
没有什么特别的原因，有一天我发现自己在想，如果我能创作出具有如此卡通魅力和充满活力的清晰图片，以进入我的论文和个人科学期刊，那该多好。因此，首先进入一个 batshit 项目，编写一个基于 gnuplot 的分子可视化器。
到目前为止，它是为我特定类型的分子量身定制的。基本上共价键合的原子形成配体，配体本身通过配位键与一些中心金属 ionic 相互作用。我已经设法得出了一个非常好的工作概念，如下图所示。

其中，虚线表示与铕金属 ionic 的配位键，以浅青色着色，实线表示原子之间的共价键。红色是氧，蓝色是氮，白色是氢，灰色是碳。到目前为止一切都很好，看起来很扎实，非常符合我想要的。
那么我该怎么做呢，我听到你问了？实际上，这很简单。我一次策划一件事情。首先，我绘制虚线的连接模式，如下所示:

然后我画出共价键:

每个步骤都需要一个或多个单独的文件。每个配体的连接性存储在一个单独的“bondfile”中，而虚线连接性模式本身就在一个文件中。具有其颜色的原子的位置被放置在另一个文件中。每个配体一个，中心金属一个。
然后我有一个单独的文件，用于金属和每个配体的原子，我在其中说明它们是什么颜色。原子被放置在黑点上的事实是在这些点周围形成迷人的黑色布局，否则它们没有轮廓线。
问题
当您想要旋转复合体以获得更好的角度以保存到图片中时，就会出现问题。为了说明这个问题，我将用单个配体的图片来展示它。让我们拿联吡啶(带氮的那个，有两个)
所以这是我认为最佳角度的联吡啶:

现在假设我们沿着下图所示的轴转动联吡啶。

现在问题出现了。因为一些应该在后面的原子实际上在整个事物的前面，这表明 gnuplot 实际上没有透视。或者，至少，它确实有，但我使用不正确。

到现在为止还挺好。我没想到它会自动获得透视，因为这不是它最初的用途。然而，这意味着 gnuplot “splot” 做了一个有点假的 3D 绘图，并且空间中点的实际相对位置无关紧要。
所以我的问题是，对于你们所有的 gnuplot/perspective 专家:有没有办法巧妙地规避这个限制？
我对任何方法都感兴趣，无论涉及多少，只要它在 gnuplot 本身的限制范围内是可行的。

Answer 1

前段时间，我尝试了类似的东西。显然，点和线不遵守 3D 顺序。但是，如果您使用表面绘制，即原子=球体和键=圆柱体，它将起作用。
编辑:这是一个完全修订的版本。
我知道有专门的程序来可视化分子。这只是为了好玩，并展示使用 gnuplot 的可行性。
我假设随着原子数量的增加，这个脚本会变得很慢。
可以直接读取结构数据文件 (SDF) 文件。它包含原子位置和键信息(连接性和键类型)。
原子显示为球体，键显示为圆柱体。因此，数据 block $Sphere和 $Cylinders包含球体和圆柱体原型(prototype)的数据点。
有关原子的附加信息存储在数据 block $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。
动画: