orthorhombic maps lose periodic boundary conditions in P1

[dennisbrookner]

Hey gemmi team,

I've encountered something puzzling in my work with real-space maps (gemmi.FloatGrid and gemmi.Ccp4Map objects) which seem to sometimes lose periodic boundary conditions. To reproduce this behavior, I am starting with two mtz files from the PDB, 1rx2_phases.mtz and 7mm1_phases.mtz. Importantly for these purposes, 1rx2 is in P 21 21 21 and 7mm1 is in P 32 2 1. These files are in the zip archive below, or can be downloaded from the PDB.

input mtzs

The behavior I observe is as follows:

• Writing 1rx2 to a P 21 21 21 map: everything works as expected. Map has periodic boundary conditions.
• Writing 1rx2 to a P 1 map: when I open this map in coot, electron density only appears inside the unit cell, and abruptly disappears outside the unit cell. E.g., this map no longer has periodic boundary conditions.
• Writing 7mm1 to a P 32 2 1 map: I get an error when I attempt to open this map in coot. Upon inspecting this map in python, I'm unable to find anything wrong with it, but something isn't working.
• Writing 7mm1 to a P 1 map: everything works as expected. Map has periodic boundary conditions.

Below is the code I used to create these four map files from the input mtzs. The code uses a mixture of reciprocalspaceship and gemmi functionality. I have linked to the reciprocalspaceship source code where relevant. I have observed this behavior with both gemmi 0.5.5 and gemmi 0.5.7. Hopefully, this code is runnable without issue with a simple pip install reciprocalspaceship and then python thisfile.py *.mtz.

I'm not sure whether this behavior is 1) a bug, 2) an error in my usage, or 3) the intended behavior for a reason that I'm not understanding. I appreciate the time, and I would love to know how to write a map in either P1 or the original spacegroup and maintain the desired periodic boundary conditions. Thanks so much! Let me know if there's any more information I should include!

Here is the code I used:

import numpy as np
import gemmi
import reciprocalspaceship as rs
import sys

def spacing_to_gridsize(spacing, cell):
    """

    Compute the optimal gridsize based on unit cell size and desired spacing

    Parameters
    ----------
    spacing : float
        Desired (approximate) grid spacing in Angstroms
    cell : gemmi.UnitCell
        Or anything similar with attributes a, b, and c
    """
    gridsize = []

    for dim in [cell.a, cell.b, cell.c]:
        gridsize.append(int(dim // spacing))

    return gridsize

def make_floatgrid(mtz, gridsize, F, Phi, spacegroup=None):
    """
    Boilerplate code for making a gemmi.FloatGrid object from an rs.DataSet object

    Parameters
    ----------
    mtz : rs.DataSet
        mtz data to be transformed into real space
    gridsize : list-like of length 3
        The size of the 3D voxel array representing your output map. Can be computed via `spacing_to_gridsize`
    F : str
        Column in mtz containing structure factor amplitudes to use for calculation
    Phi : str
        Column in mtz containing phases to be used for calculation
    spacegroup : str, optional
        Spacegroup for the output FloatGrid. If None (default), FloatGrid will inherit the spacegroup of mtz.

    Returns
    -------
    float_grid : gemmi.FloatGrid
        Fourier transform of mtz
    """

    # drop NAs in either of the specified columns
    # this has the secondary purpose of not silently modifying the input mtz
    new_mtz = mtz[~mtz[F].isnull()]
    new_mtz = new_mtz[~new_mtz[Phi].isnull()]

    # perform FFT using the desired amplitudes and phases
    new_mtz["Fcomplex"] = new_mtz.to_structurefactor(F, Phi)
    # to_structurefactor source code at https://github.com/rs-station/reciprocalspaceship/blob/e848740aabc2a128a226646b7f749479b8dd715a/reciprocalspaceship/dataset.py#L465
    reciprocal_grid = new_mtz.to_reciprocal_grid("Fcomplex", grid_size=gridsize) 
    # to_reciprocal_grid source code at https://github.com/rs-station/reciprocalspaceship/blob/e848740aabc2a128a226646b7f749479b8dd715a/reciprocalspaceship/dataset.py#L1446
    real_grid = np.real(np.fft.fftn(reciprocal_grid))

    # declare gemmi.FloatGrid object
    float_grid = gemmi.FloatGrid(*gridsize)
    float_grid.set_unit_cell(new_mtz.cell)

    if spacegroup is not None:
        float_grid.spacegroup = gemmi.find_spacegroup_by_name(spacegroup)
    else:
        float_grid.spacegroup = new_mtz.spacegroup

    # write real_grid into float_grid via buffer protocol
    temp = np.array(float_grid, copy=False)
    temp[:, :, :] = real_grid[:, :, :]

    # Enforce that mean=0, stdev=1 for voxel intensities
    float_grid.normalize()

    return float_grid

def mtz2map(filename):
    """
    Read in an mtz file and write out a .map file

    Parameters
    ----------
    filename : str
        Filename of mtz file to read in. Output .map file will get the same name
    """

    mtz = rs.read_mtz(filename)
    gridsize = spacing_to_gridsize(0.25, mtz.cell)

    # make map in original spacegroup
    fg = make_floatgrid(mtz, gridsize, F="FP", Phi="PHWT")
    ccp4 = gemmi.Ccp4Map()
    ccp4.grid = fg
    ccp4.update_ccp4_header(2, True)
    ccp4.write_ccp4_map(f"{filename[:-4]}.map")

    # make map in P1
    fg_p1 = make_floatgrid(mtz, gridsize, F="FP", Phi="PHWT", spacegroup="P1")
    ccp4_p1 = gemmi.Ccp4Map()
    ccp4_p1.grid = fg_p1
    ccp4_p1.update_ccp4_header(2, True)
    ccp4_p1.write_ccp4_map(f"{filename[:-4]}_p1.map")

    return

def main():

    args = sys.argv

    for arg in args[1:]:
        mtz2map(arg)

    return

if __name__ == "__main__":
    main()

[keitaroyam]

I think this is because Coot thinks your map is "EM" map. If alpha=beta=gamma=90 and P1, Coot considers it an EM map and doesn't show symmetry. https://github.com/pemsley/coot/blob/523de831e062e10af6c4c06b442da5ab7c40e507/src/molecule-class-info-maps.cc#L2155

[dennisbrookner]

That makes sense, thanks!

Do you know if there's a way to override that when you open the map in Coot?

[keitaroyam]

I don't think we can. It is great that you opened an issue on coot. Actually I also wanted that feature.

[dennisbrookner]

Closing this issue now that I've got a response from coot. Thanks!

https://github.com/pemsley/coot/issues/67#issuecomment-1336489684