# How to generate the remaining sites in a cubic lattice from a given POSCAR/cif in Python?

For example, I intend to generate the remaining sites of the Cr cubic lattice but the POSCAR file provides limited info on the number of atoms present as shown below. Is there a way to implement this in python (maybe using pymatgen)? It shows 8 atoms which are the total number of atoms within a unit cell volumetrically but how do I generate all the sites of Cr, if you visualize the below in VESTA you will notice that there are a total of 21 atomic sites in a lattice.

1.0
4.528916128472692 0.0 -0.0
0.0 4.528916128472692 0.0
-0.0 0.0 4.528916128472692
Cr
8
direct
0.0 0.0 0.0
0.5 0.5 0.5
0.25 0.0 0.5
0.75 0.0 0.5
0.0 0.5 0.25
0.0 0.5 0.75
0.5 0.75 0.0
0.5 0.25 0.0

• If you are using a CIF, all the atoms for the unit cell are added (8 looks ok), without the need to add/complete them manually. Normally, VESTA shows some atoms besides the unit cell, so maybe you are counting 21 because of that.
– Camps
Commented Jul 8, 2023 at 1:34
• @Camps how can I add them manually in python if I need to? I want to create the structure shown in VESTA! Commented Jul 8, 2023 at 4:45

## 1 Answer

I sort of get what OP wants to achieve based on the comment section, but what a strange request! Nonetheless, this can be easily done by both pymatgen and the Atomic Simulation Environment (ASE).

The idea is to double the structure in all 3 dimensions and then only keep the atoms that are inside or on the boundary of the original cell. To do this, you just need to find atoms that have fractional coordinate no larger than 1/2 in all 3 dimensions (this is a general solution for cubic and non-cubic systems). Below are some simple Python codes that can generate the structure you want.

## ASE

from ase.io import read, write
from ase.visualize import view

# Read in the POSCAR file
atoms = read('POSCAR')
# Store the original unit cell
old_cell = atoms.cell.copy()
# Expand the unit cell in all 3 dimensions
atoms *= (2, 2, 2)
# Find the atom indices that are inside or on the boundary of the original cell
indices = [a.index for a in atoms if (a.scaled_position <= 1/2).all()]
# Reduce the structure
atoms = atoms[indices]
# Assign the original cell back
atoms.cell = old_cell
# Visualize
view(atoms)
# Export to a cif file
write('Cr21.cif', atoms)


## Pymatgen

from pymatgen.core import Structure

struct = Structure.from_file('POSCAR')
old_lattice = struct.lattice.copy()
struct.make_supercell([2, 2, 2])
symbols, coords = zip(*[(a.species, a.coords) for a in struct
if (a.frac_coords <= 1/2).all()])
struct = Structure(old_lattice, symbols, coords)
struct.to(filename='Cr21.cif')

• This code is pretty slick. Nice job. Was wondering if you've come across a pymatgen implementation? Commented Jul 10, 2023 at 18:56
• @PranoyRay The two packages are very similar so yeah I have added a pymatgen implementation for you. I have also made the implementation general for cubic and non-cubic systems. I personally prefer ASE because it's easier for atom manipulation (and I work at DTU). Commented Jul 11, 2023 at 2:22
• Thanks for your help! Commented Jul 11, 2023 at 14:53
• This code(pymatgen) is written focused on Cr21, how to make it more generalizable for a multi element system? Basically looking for a more generalized adaptation of the line: struct = Structure(old_lattice, ['Cr']*len(coords), coords) what do I put instead of ['Cr']? assuming I have the struct object available. Commented Jul 19, 2023 at 20:25
• @PranoyRay I have editted the pymatgen code to be general for all systems. Commented Jul 20, 2023 at 6:23