I have a crystal structure in which I'd like to randomize the placement of cations from an initial structure. I have a .xyz file (given here https://pastebin.com/d5UcNr2h) which gives atomic identities and coordinates, and the following code:

structure = Import[SystemDialogInput["FileOpen"]];

    structure[[1,i]] == Atom[Cu],
    If[RandomReal[] > 0.5, structure[[1,i]] = Atom[Ga], structure[[1,i]] = Atom[Cu]],

In short, my current approach is to read in the file, and test if the atomic identity at a given position is Copper (Cu). If the atom is copper, then I roll the dice, and if RandomReal is > 0.5, then I want to substitute the atom with Gallium (Ga). I've tried the individual logical test:


which does return "True" when the atomic identity is Cu, but the code doesn't modify the actual atomic identity when I run it. What am I missing?



1 Answer 1


Import the molecule

In[14]:= m = Import["https://pastebin.com/raw/d5UcNr2h", "XYZ"];

In[15]:= OutputForm[m]

Out[15]//OutputForm= Molecule[<Cu513H210S672>]

Get a list of copper atom indices

In[16]:= copperSites = AtomList[m, "Cu", "AtomIndex"]

Out[16]= {483,484,485,486,487,488,489,490,491,492,<<493>>,986,987,988,989,990,991,992,993,994,995}

Now choose some of these at random, however you wish to pick them, and use "ReplaceAtom" to modify the molecule:

In[17]:= toReplace = RandomSample[copperSites, 10]

Out[17]= {919, 613, 543, 675, 587, 867, 715, 700, 527, 775}

In[18]:= m2 = MoleculeModify[m, {"ReplaceAtom", Thread[toReplace -> "Ga"]}, ValenceErrorHandling -> None];

During evaluation of In[18]:= Molecule::valenc: Invalid valence for gallium atom at position 527.

Out[19]//OutputForm= Molecule[<Cu503Ga10H210S672>]

You should avoid using Part when working with objects like Molecule, and instead use accessor functions like AtomList and BondList. In a future release, Molecule will be AtomQ and accessing its parts directly will not be possible.

  • $\begingroup$ Wow, I've never seen any of these techniques. Thanks for this. I tried exporting the resulting data back as a .XYZ and I'm running into issues again. Is there something in this method that is dissociating the atomic coordinates from the individual atoms? The particular error is: "Unable to generate 3D coordinates for molecule" and "There are insufficient elements to export to XYZ format". $\endgroup$ Commented Oct 27, 2020 at 18:40
  • $\begingroup$ That is a bug - MoleculeModify is removing the coordinates. I modified the answer to use ValenceErrorHandling -> None as an option to MoleculeModify and it should work. $\endgroup$
    – Jason B.
    Commented Oct 27, 2020 at 18:46
  • $\begingroup$ Excellent. Thanks so much! $\endgroup$ Commented Oct 27, 2020 at 18:53

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.