How Can I Relax a Molecule in Mathematica?

I have a molecule in Mathematica. It consists of around 70 atoms, but let us for simplicity assume that we only have four atoms:

atoms={"Au","Au","S","C"}


that have the coordinates:

coords={{3,1.1,0.4},{-0.37,1.2,0.4},{0.1,0.06,3.2},{-0.43,0.2,6.1}}


Now, I want to make a change to the molecule, this could be a rotation of the two latter atoms, {"S","C"}, around the x-axis by 60°. I want to relax the structure and find the relaxation path between these two points. I know this is something people usually do in DFT, but I have heard that it should be possible to do a rough estimate in Mathematica. So what I want to find is the positions for the four atoms along the way when this rotation is performed (this could be in steps of 10°). I have searched a lot and found that a potential method is the BFGS method, but I don't understand how to do it in practice.

I have found this introduction, but since I have never done anything like this before I do not understand how to start. Which function am I supposed to minimize? How do I get the Hessian? And how is this related to finding the new coordinates of the molecule? Can anyone please provide me with a brief overview of how to relax a structure in Mathematica or with a link to where this is described?

• I want to relax something, I try rubbing its ears. Commented May 30 at 1:04

Mathematica has a couple of force fields available for relaxing atomic coordinates. For example you can start with ethane in a nearly planar geometry and find a more reasonable geometry via

m1 = Molecule[{"C", "C", "H", "H", "H", "H", "H",
"H"}, {Bond[{1, 2}, "Single"], Bond[{1, 3}, "Single"],
Bond[{1, 4}, "Single"], Bond[{1, 5}, "Single"],
Bond[{2, 6}, "Single"], Bond[{2, 7}, "Single"],
Bond[{2, 8}, "Single"]},
AtomCoordinates -> {{-0.5, 0, -0.01}, {0.5, 0, -0.01}, {-0.5,
0.65, 0.01}, {-0.5, -0.65, -0.01}, {-1.15, 0,
0.01}, {0.5, -0.65, 0.01}, {0.5, 0.65, -0.01}, {1.15, 0,
0.01}}];
m2 = MoleculeModify[m1, "EnergyMinimizeAtomCoordinates"];
MoleculePlot3D /@ {m1, m2}


But these force fields are not applicable to metallic systems generally, so they will not work with the gold-sulfur-carbon complex described in the OP. For this you would need to use an electronic structure package.

• The force field support is cool, I think you could also do something really useful with this to begin on ES support by creating a canonical form for specifying scans / geometry optimizations and the results of such calculations that could then be specialized to automate building jobs to run on a cluster with a proper ES package (like Psi4) compiled and available, even if for licensing reasons you can't provide direct access to a Psi4 binary Commented May 28 at 17:15
• Wow, this example is neat @jason-b, but you are right it doesn't work on my example molecule (or my full molecule). Where can I find these electronic structure packages, and an introduction on how to use them? I found this though it does not seem to work with relaxation of structures.
– flg
Commented May 29 at 9:21
• @flg you'll be better off with this kind of thing psicode.org/psi4manual/master/optking.html and to make the calculations tractable for a system with two gold atoms you'll need an effective core potential psicode.org/psi4manual/master/… Commented May 30 at 17:28
• @b3m2a1 thanks for the very nice input. I think I will try to write my own routine in mathematica using the Nudged Elastic Band (NEB) method.
– flg
Commented May 31 at 8:42