Substituting Atoms in .XYZ File

Question

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.

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