# How pbc_diff works

Hello. I am relatively new to pymatgen. I am trying to compare two atoms belonging to two different verisons of the same structure. The coordinates of these sites are as follows:

atom1 -> PeriodicSite: Co (17.8204, 24.9218, 54.8488) [2.9850, 3.5300, 7.4060]
atom2 -> PeriodicSite: Co (17.8204, 24.9218, 0.0000) [2.9850, 3.5300, 0.0000]

and my lattice parameters are:

abc: 5.97 7.06 7.406

I extract the site information regarding these atoms from corresponding structures like this:

atom1_fc = str_1.sites[4].frac_coords
atom2_fc = str_2.sites[1].frac_coords

I calculate the pbc_dist between these atoms:
pbc_diff(atom1_fc, atom2_fc)
and this gives:
array([0. , 0. , 0.406001])

I donâ€™t understand why this function does not give an array of zeros. As you see the fractional coords of atom_1 indicate that the c coordinate is the lattice length in that direction. I would expect the pbc_diff of these atoms to be 0, meaning that they are located at the same point on the lattice. I visualise these atoms with VMD and indeed thatâ€™s the case.

I checked the source code of this function and it seems this function does not take the lattice params as input, therefore it really shouldnâ€™t be able to give me an array of all zerosâ€¦

I would like to understand how this function works because it is used in is_periodic_image, which therefore gives me Falseâ€¦

So my question is:
how does pbc_diff() work? maybe it is different than what I am expecting it to be.

Hello Beliz,

It appears that your PeriodicSite coordinates were entered wrongly. Fractional coordinates are typically between 0 and 1. It seems like you entered Cartesian coordinates but treated them as fractional coordinates. So, 7.4060 would correspond to 7 times the unit cell in the z-direction plus another 0.4060 times the unit cell. Hence, considering PBC, the distance (in fractional coords) between the two sites is indeed 0.4060. Therefore, the function `pbc_diff` gives the correct answer.
Just make sure to be consistent in the way you specify the coordinates (either Cartesian or fractional) and you should get the correct answer. Btw, the Pymatgen sites module already has functions that do basically the same: `distance_and_image` and `distance_and_image_from_frac_coords`. But you can also stick to `pbc_diff`, either way.

Welcome to the world of Pymatgen and good luck!

Best,
Peter

Hello Peter,

Thank you for your reply. I see the problem now, however I get these fractional coordinates with pymatgen. I read a XYZ file with IMolecule and then convert to Structure. Maybe this method is problematicâ€¦ I need to do the analysis of these files considering they have pbc but the files are in XYZ file format and Structure cannot read those files. Can you take a quick look at how I do the conversion to get the PeriodicSites and let me know if I am missing something?

``````import pymatgen
from pymatgen import Lattice, Structure, IMolecule

mol_str_1 = IMolecule.from_file('str_1.xyz')
mol_str_2 = IMolecule.from_file('str_2.xyz')

# define lattice
lattice_str = Lattice.from_parameters(a=5.970,b=7.060,c=7.406,alpha=90,beta=90,gamma=90)

# get coordinates of each molecule
mol_str_1_coords = []
for i in range(len( mol_str_1)):
mol_str_1_coords.append(list( mol_str_1.sites[i].coords))
mol_str_2_coords = []
for i in range(len( mol_str_2)):
mol_str_2_coords.append(list( mol_str_2.sites[i].coords))

# convert to Structure (periodic)
str_1 = Structure(lattice=lattice_str, species=mol_str_1.species, coords=mol_str_1_coords)
str_2 = Structure(lattice=lattice_str, species=mol_str_2.species, coords=mol_str_2_coords)
``````

And this gives me fractional coords between 0 and cell lengthsâ€¦

Is there another way I can read XYZ and construct periodic structures?

# Edit:

One important remark about `str_1` is that it is taken from a CIF file and a simulation is performed. You might also know that structures are generated from CIF files in terms of point groups and in this case, some atoms are out of the unit cell but when you consider pbc the unit cell has the correct number of atoms. Maybe pymatgen cannot understand it. If you check the `atom_1` I mentioned in the question, you can see that the atom is actually out of the unit cell but sitting on the boundary.

Little update:

I specified `coords_are_cartesian=True` inside `Structure` and everything is solved. Thanks for all your help.

1 Like

Good job, you found the solution!
Just a bit streamlined: You donâ€™t actually need to save the coordinates in a list first

``````str_1 = Structure(lattice=lattice_str, species=mol_str_1.species,
coords=mol_str_1.coords, coords_are_cartesian=True)
``````

Also, for future use, you can access the fractional coordinates directly with `frac_coords`.

Thanks!

I cannot get the coords with `mol_str_1.coords`, it raises error `'IMolecule' object has no attribute 'coords'`. I can only call `.coords` on sites, thatâ€™s why I was using the loop. The same goes for `.frac_coords`.

Oh, I see, that seems to work only for structures (I rarely work with molecules). Would `mol_str_1.sites.coords` work? If not, then I guess you already have the shortest solution!

No, it only works for one site at a time, like the ones in the for-loop

Hello Peter, I have yet another question for you. Setting `cartesian_coords = True` worked for most of my structures but unfortunately I have some more, which it does not work. Let me show you the problem.

I am again importing as `IMolecule` and then converting to `Structure` with the same code above. Here are the 15th, 16th and 17th atoms in the structure `struct`:

``````PeriodicSite: Fe (-3.0939, -4.2442, 10.4433) [-0.4678, -0.5078, 1.0859]
PeriodicSite: P (1.8263, -5.3201, 13.0661) [0.2762, -0.6365, 1.3586]
PeriodicSite: O (-1.1531, -4.1889, 13.4766) [-0.1744, -0.5012, 1.4013]
``````

As you see, I again have fractional coordinates larger than 1. While I was playing around to fix this, I noticed the functions `get_fractional_coords` and `get_cartesian_coords` within the module `Lattice`.

My lattice is:

``````lattice = Lattice.from_parameters(a=6.6131,b=8.3585,c=9.6175,alpha=90,beta=90,gamma=90)
``````

Then for the site of Fe above, I do the following:

``````>>> lattice.get_fractional_coords(struct.sites[15].coords),struct.frac_coords[15]
(array([ 0.27616452, -0.63649462,  1.35857193]), array([ 0.27616452, -0.63649462,  1.35857193]))
``````

Even when I calculate the fractional coordinates from the lattice, they are larger than 1. I donâ€™t know how to fix thisâ€¦

Hi Beliz,
The fractional coordinates are larger than 1 because your molecule is larger than your unit cell. Your molecule doesnâ€™t fit inside a single unit cell.
Depending on what you are planning to do with these structures this may or may not be a problem. Why are you using this specific lattice?
-Peter

I figured how to access the coordinates of a molecule. `mol_str_1.cart_coords` does the trick. Itâ€™s been some time but I wanted to share anyways.

1 Like