# Substituting Atoms in .XYZ File

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"]];

Do[
If[
structure[[1,i]] == Atom[Cu],
If[RandomReal[] > 0.5, structure[[1,i]] = Atom[Ga], structure[[1,i]] = Atom[Cu]],
Unevaluated[Sequence[]]
],
{i,1,Dimensions[structure[[1]]][[1]]}
]


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:

structure[[1,600]]==Atom[Cu]


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?

Thanks.

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];
OutputForm[m2]

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.

• 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". – BesselFunct Oct 27 at 18:40
• 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. – Jason B. Oct 27 at 18:46
• Excellent. Thanks so much! – BesselFunct Oct 27 at 18:53