project-gemmi / gemmi

macromolecular crystallography library and utilities
https://project-gemmi.github.io/
Mozilla Public License 2.0
205 stars 42 forks source link

crash with set_extent #307

Closed JoshRackers closed 2 months ago

JoshRackers commented 3 months ago

I am new to Gemmi, but I'm trying to cut out a chunk of electron density around a ligand with code like this:

import gemmi

mtz = gemmi.read_mtz_file("map.mtz")
m = gemmi.Ccp4Map()
m.grid = mtz.transform_f_phi_to_map('FWT','PHWT',sample_rate=3)
m.update_ccp4_header()
ligand = gemmi.read_structure("ligand.pdb")
m.set_extent(ligand.calculate_fractional_box(margin=1))
m.write_ccp4_map("masked.ccp4")

However, the script hangs and then crashes on the set_extent step. Am I doing something wrong or is this a bug?

I have tried versions 0.6.5, 0.6.4 and 0.5.8.

wojdyr commented 3 months ago

What is the span of the box returned from ligand.calculate_fractional_box(margin=1)? (it has properties .minimum and .maximum)

JoshRackers commented 2 months ago
min <gemmi.Fractional(-3.428, 9.159, -13.189)>
max <gemmi.Fractional(27.649, 40.633, 23.763)>
wojdyr commented 2 months ago

What probably happens here is that ligand.pdb has no CRYST1 line, or no unit cell dimensions there. Then ligand.calculate_fractional_box(margin=1) returns a box in Cartesian coordinates, the map extent gets set to a huge size, and the script runs out of memory. I'll add a check to prevent it. A possible workaround (untested) is to set ligand.cell = mtz.get_cell().

JoshRackers commented 2 months ago

I think you are right. When I add ligand.cell = mtz.get_cell(), the min and max change too:

min <gemmi.Fractional(-0.0586284, 0.105958, -0.287907)>
max <gemmi.Fractional(0.472875, 0.470072, 0.51873)>

and the script runs to completion. Thanks for the help!