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.

4.528916128472692 0.0 -0.0
0.0 4.528916128472692 0.0
-0.0 0.0 4.528916128472692
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 
  • $\begingroup$ 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. $\endgroup$
    – Camps
    Jul 8, 2023 at 1:34
  • $\begingroup$ @Camps how can I add them manually in python if I need to? I want to create the structure shown in VESTA! $\endgroup$
    – Pranoy Ray
    Jul 8, 2023 at 4:45

1 Answer 1


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.


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
# Export to a cif file
write('Cr21.cif', atoms)


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)
  • 1
    $\begingroup$ This code is pretty slick. Nice job. Was wondering if you've come across a pymatgen implementation? $\endgroup$
    – Pranoy Ray
    Jul 10, 2023 at 18:56
  • $\begingroup$ @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). $\endgroup$
    – Shaun Han
    Jul 11, 2023 at 2:22
  • $\begingroup$ Thanks for your help! $\endgroup$
    – Pranoy Ray
    Jul 11, 2023 at 14:53
  • $\begingroup$ 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. $\endgroup$
    – Pranoy Ray
    Jul 19, 2023 at 20:25
  • 1
    $\begingroup$ @PranoyRay I have editted the pymatgen code to be general for all systems. $\endgroup$
    – Shaun Han
    Jul 20, 2023 at 6:23

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .