RFE: "intelligent" atom transform

[Tristan Croll]

Would be great to have a function to apply a rigid-body transformation to a selection of atoms, that also transforms any altlocs and correctly adjusts the aniso_u values (the 3x3 aniso_u matrix U should be transformed by (rot).U.(rot)T). The need pops up fairly regularly (creating an assembly from a single asymmetric unit; rearranging which symmetry copies constitute "the" asymmetric unit; etc.).

[Eric Pettersen]

Ressigning to T.G., since it is undoubtedly his code that would make use of any such call, so at the very least he would be the one to suggest an API for this (which I might then be the one implementing).

[Tom Goddard]

I think a rigid motion of a subset of atoms of a structure is very rare. I don't think I have any code that does that (maybe fitmap can do it, but I'm sure it ignores alt locs and aniso_u). Seems reasonable to have some API that would transform all altlocs and aniso_u. This is a small bit of code in Python but would be slow, having to set the current altloc cycling through all of them. Maybe Tristan needs better performance having it done in C++?

[tic20@…]

I’m on a bit of an optimisation kick at the moment, moving various functions to C++ (including the conceptually somewhat similar task of copying atoms into Clipper’s data structures for map calculations). I can whip an implementation together and pass it on, if you like.
 

 

[Eric Pettersen]

Since you know what API you want and I don't -- sounds like a plan!

[Tristan Croll]

This appears to be working (coordinates including altlocs transform sensibly, and transformed aniso_u values checked against Clipper's implementation). If the following is in molc.cpp:

template <class C, typename T>
void affine_transform(const C& coord, T* tf, C& result)
{
    for (size_t i=0; i<3; ++i)
    {
        result[i] = tf[4*i] * coord[0] + tf[4*i+1] * coord[1] + tf[4*i+2]*coord[2] + tf[4*i+3];
    }
}

template <typename T>
void transform_u_aniso(const std::vector<float>* aup, T* tf, std::vector<float>& result)
{
    // Need to apply only rotation component of transform, as (rot).U.(rot)T
    // aniso_u6 is stored as u11, u12, u13, u22, u23, u33
    const auto& au = *aup;
    std::vector<float> full = {au[0],au[1],au[2],au[1],au[3],au[4],au[2],au[4],au[5]};
    std::vector<float> ir(9);
    for (size_t i=0; i<3; ++i){
        for (size_t j=0; j<3; ++j){
            ir[3*i+j] = tf[4*i] * full[j] + tf[4*i+1] * full[3+j] + tf[4*i+2] * full[6+j];
        }
    }
    for (size_t i=0; i<3; ++i) {
        for (size_t j=0; j<3; ++j){
            result[3*i+j] = ir[3*i] * tf[4*j] + ir[3*i+1] * tf[4*j+1] + ir[3*i+2] * tf[4*j+2];
        }
    }
}

void transform_atom(Atom* atom, double* tf)
{
    Coord transformed;
    affine_transform<Coord, double>(atom->coord(), tf, transformed);
    atom->set_coord(transformed);
    if (atom->has_aniso_u())
    {
        std::vector<float> ua(9);
        transform_u_aniso<double>(atom->aniso_u(), tf, ua);
        atom->set_aniso_u(ua[0],ua[1],ua[2],ua[4],ua[5],ua[8]);
    }
}

extern "C" EXPORT void transform_atoms(void* atom, size_t n, double* tf)
{
    try {
        auto a = static_cast<Atom**>(atom);
        char current_altloc;
        for (size_t i=0; i<n; ++i)
        {
            auto atom = *(a++);
            auto altlocs = atom->alt_locs();
            if (altlocs.size())
            {
                current_altloc = atom->alt_loc();
                for (const auto& altloc: altlocs)
                {
                    atom->set_alt_loc(altloc);
                    transform_atom(atom, tf);
                }
                atom->set_alt_loc(current_altloc);
            } else {
                transform_atom(atom, tf);
            }
        }
    } catch (...) {
        molc_error();
    }
}

... then on the Python side, using the molobject/molc API, where tf is a Place:

def transform_atoms(atoms, tf):
    f = c_function('transform_atoms',
        args=(ctypes.c_void_p, ctypes.c_size_t, ctypes.POINTER(ctypes.c_double)))
    f(atoms._c_pointers, len(atoms), pointer(tf.matrix))

I'd suggest implementing it as a method of the Atoms class (e.g. Atoms.transform(place)).

[Eric Pettersen]

Okay, this patch has been applied to both the develop and release branches, and is available via Atoms.transform(place).