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 :confused:

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