Minimum Distance in Atom Position Refinement

[Julian Sarro]

Dear PDF Community,

I am currently working on a project where we model a disordered silver cluster using PDF analysis.
In one of our approaches, we refine the atomic positions. The variables in our script are defined as follows:


    atoms = pdfgenerator.phase.getScatterers()
    for atom in atoms:
        if atom.element == 'Ag':
            recipe.constrain(atom.Biso, AgBiso)
            recipe.restrain(atom.Biso, lb=0, ub=1, sig=0.001)
            recipe.addVar(atom.x, name = atom.name+'_x', tags = [atom.name, 'Ag_all'])
            recipe.restrain(atom.x, lb=atom.x.value-1, ub=atom.x.value+1, sig=0.001)
            recipe.addVar(atom.y, name = atom.name+'_y', tags = [atom.name, 'Ag_all'])
            recipe.restrain(atom.y, lb=atom.y.value-1, ub=atom.y.value+1, sig=0.001)
            recipe.addVar(atom.z, name = atom.name+'_z', tags = [atom.name, 'Ag_all'])
            recipe.restrain(atom.z, lb=atom.z.value-1, ub=atom.z.value+1, sig=0.001)
 
Due to termination ripples in the data, we occasionally find that the refinement produces structures with unphysically short Ag–Ag distances.

Is there a way to set a minimum distance between atoms to prevent them from being too close to each other?

Any insights or suggestions would be greatly appreciated!

Cheers,
Julian

[Tomasz Stawski]

Hi Julian,

yes there is. It follows the following idea.

 # Create new variable for the bond length
            fit.newVar(var_name)
           
            # Constrain the bond length variable to the expression
            fit.constrain(var_name, expression)
                   
            # Adjust the weight of the restraint
            fit.restrain(var_name, lb=0.9*relative_length, ub=1.1*relative_length, scaled = False, sig=sig)

Here the expression could be adapted from this function:

#----------------------------------------------------------------------------
def create_constrain_expressions_for_bonds(bond_pairs, phase):
    constrain_dict = {}
    for bond in bond_pairs:
        central_label = bond['central_label']
        oxygen_label = bond['oxygen_label']
        relative_length = bond['relative_length']

        var_name = f'bond_length_{central_label}_{oxygen_label}_{phase}'
        bond_expr = f"((x_{central_label}_{phase} - x_{oxygen_label}_{phase})**2 + " \
                    f"(y_{central_label}_{phase} - y_{oxygen_label}_{phase})**2 + " \
                    f"(z_{central_label}_{phase} - z_{oxygen_label}_{phase})**2)**0.5"
        constrain_dict[var_name] = {
            'atoms': (central_label, oxygen_label),
            'expression': bond_expr,
            'relative_length': relative_length
        }

    return constrain_dict

I use "relative_length" derived from the atom positions in a unit cell. An important aspect is that specific atom positions need to be already added as variables before you can constrain them with a specific expression. My code is a bit extensive because I have a massive unit cell, but for small cases you may add them one-by-one. 

I hope that my specific example can be adjusted to your needs. Actually, I have spent last several months on writing custom rigid body restrictions for DiffPy and I;ve developed also a few ideas how to work with angle constraints, and how to change space groups semi-dynamically, which can be useful for bigger problems. 

Cheers,

Tomasz

--
You received this message because you are subscribed to the Google Groups "diffpy-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to diffpy-users...@googlegroups.com.
To view this discussion visit https://groups.google.com/d/msgid/diffpy-users/0e9ea7c4-47e5-454c-8d83-b4dcd259d97dn%40googlegroups.com.

[Simon Billinge]

Thanks Tomasz, this is great!

S

To view this discussion visit https://groups.google.com/d/msgid/diffpy-users/CABuBjThx1TAU_C2Fw2tdjqw7pzRi%3DwBQQS-rAzA5H4xS78Yppw%40mail.gmail.com.

--

Simon Billinge

Professor, 

Department of Applied Physics and Applied Mathematics

Columbia University

[Julian Sarro]

Hi Tomasz,

Thank you for your reply! Since I’m fairly new to Diffpy, I’m having some difficulty implementing your code into my fit recipe.
I’m also unsure whether your approach is specifically designed for crystal structures, as I’m working with discrete modeling using the DebyePDFGenerator.

Perhaps I could also simplify things, as my model only includes one type of atom (Silver).

Best regards,
Julian

[Tomasz Stawski]

Hi Julian,
yes it works with crystal structures. My code is an excerpt from a longer script You should follow these steps to generalise it:
1. add all the atom positions as variables. it depends on your approach, depending on whether you use sgpars or getScatterers().
2. Double check that you identify correctly all the names of atoms because for instance sgpars will generate different names for different space groups for the same fundamental structures.

3. You can create constrain expressions only for the atoms which are already added as variable to your fir recipe.

4. for each bond you wish to constrain create a new variable e.g. fit.newVar(bond1....
5. for the expression use a vector expression such as:

        bond_expr = f"((x_{central_label}_{phase} - x_{oxygen_label}_{phase})**2 + " \
                    f"(y_{central_label}_{phase} - y_{oxygen_label}_{phase})**2 + " \
                    f"(z_{central_label}_{phase} - z_{oxygen_label}_{phase})**2)**0.5"

6. labels in the expression need to match exactly the names of your variable in the fit. You may use dir(fit) to figure it out.
7. For a bigger cell you may want automate this procedure (as I did).

hope that helps.
Cheers,
Tomasz

To view this discussion visit https://groups.google.com/d/msgid/diffpy-users/65f946ce-50df-42db-aa6c-b9fe6b9e626an%40googlegroups.com.

[Julian Sarro]

Hi Tomasz,

Thank you so much for your help! I managed to implement it, though it’s probably not the most elegant or computationally efficient solution. Still, it gets the job done.

I might refine the code a bit further, but here’s my implementation in case it helps anyone facing a similar issue:

    atoms = pdfgenerator.phase.getScatterers()

    for atom in atoms:

        if atom.element == 'Ag':

            recipe.addVar(atom.x, name = atom.name+'_x', tags = [atom.name, 'Ag_all'])

            recipe.restrain(atom.x, lb=atom.x.value-1, ub=atom.x.value+1, sig=0.001)

            recipe.addVar(atom.y, name = atom.name+'_y', tags = [atom.name, 'Ag_all'])

            recipe.restrain(atom.y, lb=atom.y.value-1, ub=atom.y.value+1, sig=0.001)

            recipe.addVar(atom.z, name = atom.name+'_z', tags = [atom.name, 'Ag_all'])

            recipe.restrain(atom.z, lb=atom.z.value-1, ub=atom.z.value+1, sig=0.001)

    for atom in atoms:

        for atom2 in atoms:

            if atom2.name > atom.name: # Each pair distace once

                bond_expr = f"(({atom.name}_x - {atom2.name}_x)**2 + " \

                    f"({atom.name}_y - {atom2.name}_y)**2 + " \

                    f"({atom.name}_z - {atom2.name}_z)**2)**0.5"

               

                bond_dist = (

                    (atom.x.value - atom2.x.value)**2 +

                    (atom.y.value - atom2.y.value)**2 +

                    (atom.z.value - atom2.z.value)**2

                )**0.5

                if bond_dist < 4: # Only constrain when bond dist is already close

                    bond_name = f"d_{atom.name}_{atom2.name}"

                    recipe.newVar(bond_name)

                    recipe.constrain(bond_name, bond_expr)

                    recipe.restrain(bond_name, lb=2.7, sig=0.001)

Best regards,

Julian

[Tomasz Stawski]

Hi Julian,

It looks good, just  some comments.
1. check your bond distance values. I think that if you work with atom coordinates from getScatterers(), they are probably expresses as relative values (to unit cell lattice constants, hence taking values from 0 to 1). Therefore, you should calculate a relative bond length value, and use it as a boundary.  So it is probably not 2.7 but 0.0something.

DiffPy implements constrains as penalties to a residual. It helps a lot to play with a value of sig.

The "default" minimiser leastsq is IMHO quite crappy. I find it a lot more efficient to use something from the minimize library. You can use it like so:

from scipy.optimize import minimize

optimizer = 'L-BFGS-B'
minimize(fit.scalarResidual, fit.values, method=optimizer)

(so the trick is to use scalarResidual, instead of fit.residual)

Actually, with this small modification you can start playing with other libraries such as PyTorch, which allow for the use of GPUs. I am currently exploring how this works.

Also DiffPy has a very good support for parallel processing. There are some examples around how to use it. 

Cheers,

Tomasz

To view this discussion visit https://groups.google.com/d/msgid/diffpy-users/98057a93-2e71-49a0-a4a2-ca283bf48afen%40googlegroups.com.

[Julian Sarro]

Hi Tomasz,

Thank you for your feedback! Since I am working with discrete modeling (without a unit cell), the bond distance values are expressed in Angstroms.

I’ve also experimented with the sig parameter, and it indeed makes a noticeable difference—not only in the results but also in computation time in my case.

For most of my refinements, I rely on the least_squares method from scipy.optimize. Currently, it meets my needs, and everything seems to be working well. But, I’ll keep the minimize function in mind for future use.

Best
Julian