moorhen-coot / Moorhen

A JavaScript molecular graphics program
https://moorhen.org
BSD 3-Clause "New" or "Revised" License
7 stars 3 forks source link

`molecules_container::get_diff_diff_map_peaks` currently not working #281

Open FilomenoSanchez opened 7 months ago

FilomenoSanchez commented 7 months ago

We seem to be getting an empty vector back from molecules_container::get_diff_diff_map_peaks. Here's the test that is currently failing:

test.skip("Test get_diff_diff_map_peaks", () => {
        const molecules_container = new cootModule.molecules_container_js(false)
        const coordMolNo = molecules_container.read_pdb('./5a3h.pdb')
        const mapMolNo = molecules_container.read_mtz('./5a3h_sigmaa.mtz', 'FWT', 'PHWT', 'FOM', false, false)
        const diffMapMolNo = molecules_container.read_mtz('./5a3h_sigmaa.mtz', 'DELFWT', 'PHDELWT', 'FOM', false, true)

        molecules_container.associate_data_mtz_file_with_map(mapMolNo, './5a3h_sigmaa.mtz', 'FP', 'SIGFP', 'FREE')
        molecules_container.connect_updating_maps(coordMolNo, mapMolNo, mapMolNo, diffMapMolNo)
        molecules_container.sfcalc_genmaps_using_bulk_solvent(coordMolNo, mapMolNo, diffMapMolNo, mapMolNo)

        molecules_container.get_r_factor_stats()
        molecules_container.get_map_contours_mesh(mapMolNo, 77.501, 45.049, 22.663, 13, 0.48)

        molecules_container.delete_using_cid(coordMolNo, "/1/A/300/*", "LITERAL")
        molecules_container.get_map_contours_mesh(mapMolNo, 77.501, 45.049, 22.663, 13, 0.48)
        molecules_container.get_r_factor_stats()

        const diff_diff_map_peaks = molecules_container.get_diff_diff_map_peaks(diffMapMolNo, 77.501, 45.049, 22.663)
        expect(diff_diff_map_peaks.size()).toBeGreaterThan(0)
    })
pemsley commented 7 months ago

5a3h_sigmaa.mtz doesn't contain an R-free column hence the function is working as expected - no recalculation of the map and no differences in the difference map.

Capture the return value of sfcalc_genmaps_using_bulk_solvent() and test for r_factor = -1. That will tell you that the structure factor calculation failed.

I will let you close this.

FilomenoSanchez commented 6 months ago

Hi @pemsley, the return value of sfcalc_genmaps_using_bulk_solvent seems to suggest there's nothing wrong with the structure factor calculation (at least as far as I can tell based on the r_factor field):

    {
      r_factor: 0.10371405631303787,
      free_r_factor: 0.15522709488868713,
      bulk_solvent_volume: 0.25201278924942017,
      bulk_correction: 0.2925781309604645,
      n_splines: 16
    }

Also 5a3h_sigmaa.mtz does have a Rfree column (albeit mtzdmp prints weird asterisks that I haven't seen before):

$ mtzdmp ./5a3h_sigmaa.mtz

[ ... ]

 Col Sort    Min    Max    Num      %     Mean     Mean   Resolution   Type Column
 num order               Missing complete          abs.   Low    High       label 

   1 ASC      0      30      0  100.00     11.2     11.2  14.90   1.82   H  H
   2 NONE     0      38      0  100.00     14.4     14.4  14.90   1.82   H  K
   3 NONE     0      42      0  100.00     15.7     15.7  14.90   1.82   H  L
   4 NONE      0.*********     0  100.00******************  14.90   1.82 I  FREE
   5 NONE   11.9  1539.1     0  100.00   184.69   184.69  14.90   1.82   F  FP
   6 NONE    2.1    68.9     0  100.00    12.83    12.83  14.90   1.82   Q  SIGFP
   7 NONE    0.1  2864.1     0  100.00   206.89   206.89  14.90   1.82   F  FC
   8 NONE    0.0   360.0     0  100.00   176.87   176.87  14.90   1.82   P  PHIC
   9 NONE    0.0  1451.5     0  100.00   173.86   173.86  14.90   1.82   F  FC_ALL
  10 NONE    0.0   360.0     0  100.00   178.31   178.31  14.90   1.82   P  PHIC_ALL
  11 NONE    0.0  1626.7     0  100.00   176.99   176.99  14.90   1.82   F  FWT
  12 NONE    0.0   360.0     0  100.00   178.33   178.33  14.90   1.82   P  PHWT
  13 NONE    0.0   858.1     0  100.00    54.44    54.44  14.90   1.82   F  DELFWT
  14 NONE    0.0   360.0     0  100.00   178.24   178.24  14.90   1.82   P  PHDELWT
  15 NONE  0.001   1.000     0  100.00    0.882    0.882  14.90   1.82   W  FOM
  16 NONE    0.1  2062.3     0  100.00   203.85   203.85  14.90   1.82   F  FC_ALL_LS
  17 NONE    0.0   360.0     0  100.00   178.38   178.38  14.90   1.82   P  PHIC_ALL_LS

And finally, if instead of 5a3h.pdb and 5a3h_sigmaa.mtz I use moorhen-tutorial-structure-number-1.pdb and moorhen-tutorial-map-number-1.mtz I still see the same behaviour (test does not pass).

pemsley commented 6 months ago

I have a different version of 5a3h_sigmaa.mtz to you, it seems. OK.