2
$\begingroup$

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?

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

1 Answer 1

5
$\begingroup$

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}

enter image description here

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.

$\endgroup$
4
  • $\begingroup$ 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 $\endgroup$
    – b3m2a1
    Commented May 28 at 17:15
  • $\begingroup$ 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. $\endgroup$
    – flg
    Commented May 29 at 9:21
  • $\begingroup$ @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/… $\endgroup$
    – b3m2a1
    Commented May 30 at 17:28
  • $\begingroup$ @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. $\endgroup$
    – flg
    Commented May 31 at 8:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.