How to Correctly Add Perspective to a Gnuplot 3D connected Point Cloud Representing a Molecule?

Question

I don't know about you people but I love me some Gnuplot. Properly used, that software produces beautiful images, charming in their simplicity and clarity, that I am very fond of.

For no particular reason, one day I caught myself thinking how good would it be if I could create pictures of such cartoonish charm and vibrant clarity to go into my papers and personal scientific journal. So went head first into a batshit project to code a gnuplot-based molecule visualizer.

So far it is tailor made for my specific type of molecule. Basically covalently bonded atoms that form ligands, which themselves interact with some central metal ions via coordination bonds. I have managed to arrive at a pretty good working concept, pictured down below.

In it, the dotted lines denote a coordination bond with a metallic ion of Europium, colored in light cyan, the solid lines are covalent bonds between atoms. Red is oxygen, blue is nitrogen, white is hydrogen and grey is carbon. So far so good, seems pretty solid and very much in line with what I wanted.

So how do I do it, I hear you asking? Well that is pretty simple, actually. I plot things one at a time. First I plot the connectivity pattern of the dotted lines, like so:

Then I paint in the covalent bonds:

Each of the steps requires one or more separate files. The connectivity of each ligand is stored in a separate "bondfile" and the dotted conectivity pattern is in a file all of itself. The positions of the atoms with the color that they have is placed in yet another file. One for each ligand and one for the central metal.

Then I have a separate file for the atoms of the metal and of each ligand, where I say what color they are. The fact that the atoms are placed over the black dots is what gives the charming black layout around the points, otherwise they have no contour line.

THE ISSUE

The issue arises when you want to rotate the complex to get a better angle to save into the picture. In ordert o illustrate the problem I am going to show it in action with the picture of a single ligand. Let's take the bipyridine (the one with the nitrogens, there are two of them)

So here is a bipyridine in what I think is its optimal angle:

Now suppose we turn the bipyridine along the axis shown in the figure below.

Now the problem shows up. Because some atoms that should be in the behind plane are in fact in front of the entire thing, revealing that gnuplot actually don't have perspective. Or, at the very least, that it indeed have but I am using it incorrectly.

So far so good. I didn't expect it to have perspective automagically since this is not what it was originally made for. However, that means that gnuplot "splot" does a somewhat fake 3D plotting and that the actual relative positions of the points in space matters little.

So my question is, for all of you gnuplot/perspective savants out there:Is there a way to cleverly circumvent this limitation?

I am interested in any method, however involved it may be as long as it is feasible within the limitations of gnuplot itself.

Answer 1

Some time ago, I tried something similar. Apparently, points and lines do not respect the 3D-order. However, it will work if you draw with surfaces, i.e. atoms=spheres and bonds=cylinders.

Edit: This is a completely revised version. I'm aware that there are dedicated programs to visualize molecules. This is just for fun and to demonstrate the feasibility with gnuplot. I assume that this script will get pretty slow with increasing number of atoms.

A structure-data file (SDF) file can be read directly. It contains the atomic positions and bond information (connectivity and type of bond). Atoms are displayed as spheres and bonds as cylinders. Hence, the datablocks $Sphere and $Cylinders contain datapoints of a sphere and cylinder prototype. Additional information about atoms are stored in the datablock $Elements, i.e. atomic number, element name, atom size and color. More elements can be added to this list. Spheres are simply plotted with an offset according to their position. Bonds need also to be rotated appropriately which requires rotation of the bond vectors. Therefore, the following basic vector and matrix operations are implemented as functions:

• VectorLength(V)
• CrossProduct(a,b)
• MatrixVectorMultiplication(M,V)
• VectorNormalize(V)

The approach might not be the most efficient way, since vectors and matrices are handled as strings (with 3 and 9 tokens).

As an illustrative example, the data of a Caffeine molecule is taken from here.

Data: 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

$$$$

Code:

### 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

Result: (wxt terminal under Windows7, gnuplot 5.2.8)

An animation can be done by using terminal gif animate, however, I got better looking results by creating PNGs with terminal pngcairo and putting them together to an animated gif with the software ScreenToGif.

Animation:

Answer 2

Heh. I'm a molecular graphics wonk myself, having written viewers and visualization tools since grad student days in the 1970s. And you know what? I really dislike the use of perspective in molecular graphics. So much so that I'll call the absence in gnuplot a feature rather than a limitation.

There is demo molecule.dem in the gnuplot collection showing simple molecular graphics. It works better in the development version of gnuplot (5.3), where you can use plotting style "with circles" rather than "with points" for the atoms. Here you go:

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:($5-$1):($6-$2):($7-$3) with vectors nohead lw 3 lc "black"

Notes:

• The atom positions are read directly from a PDB file
• Atom coloring is generated from the atom name (gray for carbon, blue for nitrogen, etc)
• The bonds were generated from the same PDB file by the "bonds" utility in the Raster3D molecular graphics package
• Occlusion of the atoms in back by the ones in front is handled by "set hidden3d"
• Occlusion of the bonds is less satisfactory because the line segment is drawn all the way to the atom center whereas visually it would look better to terminate at the projected spherical surface of the atom.
• Visual impression of depth is further helped by making the atoms partially transparent.