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: