Definition of position delta in compare_structures

[atztogo]

Environment

• aiida-vasp version/branch+commit: Latest develop (c41eef0)

Question

I was looking for a function to compare two structures (e.g. of StructureData). I found a function aiida_vasp.utils.workchains.compare_structures. This is used for checking the convergence of crystal structure in relax workchain. In VASP, atomic position in fractional coordinates may be set in the interval of [0, 1]. So if position of an atom cross the border (0 or 1) under the relaxation or numerical treatment, the value can change by +1 or -1. In the following code, I am wondering if this situation is considered or not.

aiida-vasp/aiida_vasp/utils/workchains.py

Lines 80 to 87 in c41eef0

	pos_a = np.array([site.position for site in structure_a.sites])
	pos_b = np.array([site.position for site in structure_b.sites])
	delta.absolute.pos = pos_a - pos_b
	site_vectors = [delta.absolute.pos[i, :] for i in range(delta.absolute.pos.shape[0])]
	a_lengths = np.linalg.norm(pos_a, axis=1)
	delta.absolute.pos_lengths = np.array([np.linalg.norm(vector) for vector in site_vectors])
	delta.relative.pos_lengths = np.array([np.linalg.norm(vector) for vector in site_vectors]) / a_lengths

Is any trick implemented, for example, in parsevasp or StructureData?

[zhubonan]

I don't think that is has been considered. For cases where the atom can cross the boundary, it is very likely that the workchain will never converge, despite the relaxation is converged. The problem is that the output structure is alway wrapped such that the atoms are inside the units cell.

A fix can be done by computing the true displacement taking minimum image convention into account.

Using the routione in ase, the function can be updated as:

def compare_structures(structure_a, structure_b):
    """Compare two StructreData objects A, B and return a delta (A - B) of the relevant properties."""

    delta = AttributeDict()
    delta.absolute = AttributeDict()
    delta.relative = AttributeDict()
    volume_a = structure_a.get_cell_volume()
    volume_b = structure_b.get_cell_volume()
    delta.absolute.volume = np.absolute(volume_a - volume_b)
    delta.relative.volume = np.absolute(volume_a - volume_b) / volume_a

    # Check the change in positions taking account of the pbc
    atoms_a = structure_a.get_ase()
    atoms_a.set_pbc(True)
    atoms_b = structure_b.get_ase()
    atoms_b.set_pbc(True)
    n_at = len(atoms_a)
    atoms_a.extend(atoms_b)
    pos_change_abs = np.zeros((n_at, 3))
    for isite in range(n_at):
        dist = atoms_a.get_distance(isite, isite + n_at, mic=True, vector=True)
        pos_change_abs[isite] = dist

    pos_a = np.array([site.position for site in structure_a.sites])
    # pos_b = np.array([site.position for site in structure_b.sites])
    delta.absolute.pos = pos_change_abs

    site_vectors = [delta.absolute.pos[i, :] for i in range(delta.absolute.pos.shape[0])]
    a_lengths = np.linalg.norm(pos_a, axis=1)
    delta.absolute.pos_lengths = np.array([np.linalg.norm(vector) for vector in site_vectors])
    delta.relative.pos_lengths = np.array([np.linalg.norm(vector) for vector in site_vectors]) / a_lengths

    cell_lengths_a = np.array(structure_a.cell_lengths)
    delta.absolute.cell_lengths = np.absolute(cell_lengths_a - np.array(structure_b.cell_lengths))
    delta.relative.cell_lengths = np.absolute(cell_lengths_a - np.array(structure_b.cell_lengths)) / cell_lengths_a

    cell_angles_a = np.array(structure_a.cell_angles)
    delta.absolute.cell_angles = np.absolute(cell_angles_a - np.array(structure_b.cell_angles))
    delta.relative.cell_angles = np.absolute(cell_angles_a - np.array(structure_b.cell_angles)) / cell_angles_a

    return delta

However, I think the relative position displacement still does not make any sense, if an atom is at [0.001, 0.001, 0.001] in the first place, then any displacement can be very large! We should disallow any "relative" displacement comparison.

[atztogo]

@zhubonan, thank you for your reply. I am sorry for my late reply.
I prefer not relying on ASE. In our case, we can assume the displacements are small. So in this case, I transform position in Cartesian coordinates into fractional by

frac_position = np.dot(cart_position, np.linalg.inv(cell))
frac_position -= np.rint(frac_position)
distance = np.linalg.norm(np.dot(frac_position, cell))

Here, I use 3x3 matrix of basis vectors that are stored as row vectors. The second line brings the fractional position into the interval of [-0.5, 0.5], and I expect it is near zero for this purpose.

I agree with you that relative position displacement is not a good measure because choice of basis vectors is not unique. Also for the angles, but so we may have to discuss on the measures and their definitions.

[zhubonan]

Good point. I think the assumption of small displacement is sound here. If not the case, the relaxation should restart again anyway.