There was a recent question on the stack exchange was about how to find dihedral angles from an XYZ file, which contains atom types and 3D coordinates for a molecule or group of molecules. A dihedral angle is the angle between two intersecting planes. In chemistry, we define the dihedral among four atoms as the angle between the plane containing the first three atoms and the last three atoms. I got to wishing that I could do it interactively like you can in other chemistry-specific software.
This is my attempt at such a tool, code below. You can query for bond lengths, bond angles, and dihedral angles. You could look at the angles found in our 3D model of ibuprofen,
or compare the bond angles in five-memebered and six-membered rings in buckminsterfullerene,
Here is the code to generate these plots. We use EventHandler
along with MousePosition["Graphics3DBoxIntercepts"]
to correlate mouse-clicks with atom selection, borrowing from this post
interactiveBondTool[chemical_String]:= Module[{plot, coords, atomLabels, bonds},
{plot, coords, atomLabels, bonds} = EntityValue[
Entity["Chemical",chemical],
{"MoleculePlot", "AtomPositions", "VertexTypes", "EdgeRules"}
];
If[
Head /@ {plot, coords, atomLabels, bonds} === {Graphics3D, List, List, List},
interactiveBondTool[ {plot, coords, atomLabels, bonds} ],
Missing["NotAvailable"]
]
];
interactiveBondTool[{plot_, coords_, atomLabels_, bonds_}] := Module[
{dihedralFromVectors, dihedralFromAtomNumbers, bondLength, bondAngle,
findAtomNearestToLine, bondInfoBox},
dihedralFromVectors[{b1_, b2_, b3_}] := Module[{n1, n2},
(*http://math.stackexchange.com/a/47084/210969*)
n1 = Normalize@Cross[b1, b2];
n2 = Normalize@Cross[b2, b3];
ArcTan[n1.n2, Cross[n1, Normalize@b2].n2]
];
dihedralFromAtomNumbers[{a1_,a2_,a3_,a4_}]:=dihedralFromVectors[
(Subtract@@coords[[#]])&/@{{a1,a2},{a2,a3},{a3,a4}}
];
bondLength[{a1_,a2_}]:=EuclideanDistance@@coords[[{a1,a2}]];
bondAngle[{a1_,a2_,a3_}]:=VectorAngle @@ ((Subtract@@coords[[#]]) &/@ {{a2,a1},{a2,a3}});
findAtomNearestToLine[{v1_,v2_},pts_]:=Module[{nearestFunc},
(* adapted from this answer: http://mathematica.stackexchange.com/a/28004/9490 *)
nearestFunc=Function[{u},Norm/@({#/10,u-v1-#}&@Projection[u-v1,v2-v1])];
First@Nearest[(nearestFunc/@pts)->pts,{0,0}]
];
findAtomNearestToLine[None,pts_]:=Nothing;
bondInfoBox[pts_]:=
Grid[{
{"atom (atom number)",Grid@
Thread[{atomLabels[[pts]],pts,{Red,Yellow,Green,Blue}[[;;(Length@pts)]]}]},
{"bond length", If[Length@Union@pts>1,
(bondLength@pts[[;;2]])/100,""]},
{"bond angle", If[Length@Union@pts>2,
(bondAngle@pts[[;;3]])/Degree,""]},
{"dihedral angle", If[Length@Union@pts>3,
(dihedralFromAtomNumbers@pts[[;;4]])/Degree,""]}
},
Frame->All
];
DynamicModule[
{clicked={},atoms={},spheres={},atomlabels={}},
atoms =Dynamic[Flatten[Position[coords,#]&/@clicked]];
atomlabels:=With[
{pos=atoms},
If[pos==={},
{},
atomLabels[[#]]&/@pos
]
];
spheres=Dynamic[
Transpose[{
{Red,Yellow,Green,Blue}[[;;Length@clicked]],
Sphere[#,40]&/@clicked }
]
];
EventHandler[
Row[{
MouseAppearance[
Show[
plot,
Graphics3D[spheres],
ImageSize->500
],
"Arrow"],
Dynamic@bondInfoBox[Setting[atoms]]
}],
{"MouseClicked":>
If[
Length@clicked===4,
clicked={},
AppendTo[clicked,
findAtomNearestToLine[MousePosition["Graphics3DBoxIntercepts"],coords]
]
]
},PassEventsDown->True
]
]
];