Message Boards Message Boards

Interactively query bond information

Posted 8 years ago

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,

enter image description here

or compare the bond angles in five-memebered and six-membered rings in buckminsterfullerene,

enter image description here

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
       ]
    ]
];
POSTED BY: Jason Biggs
5 Replies
Posted 8 years ago

Hi thanks for you reply. I am using 11

POSTED BY: Marie Sylla
Posted 8 years ago

Hi, so I am trying to do the same thing for benzene but it would not work.

POSTED BY: Marie Sylla

I'm glad to help troubleshoot - what kind of errors are you getting? What version of Mathematica are you using?

POSTED BY: Jason Biggs
Posted 8 years ago

I am only seeing the box of bond length and stuff. Nothing is showing up

POSTED BY: Marie Sylla

enter image description here - another post of yours has been selected for the Staff Picks group, congratulations! We are happy to see you at the top of the "Featured Contributor" board. Thank you for your wonderful contributions, and please keep them coming!

POSTED BY: EDITORIAL BOARD
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract