Re: [ccp4bb]: bulk solvent corrected fo-fc map?

[Kevin Cowtan <cowtan@ysbl.york.ac.uk>]

Not sure how to do this in CCP4 if refmac doesn't give you what you want (have 
you fiddled with the bulk solvent options?).

But I can do it in Clipper. Here is a little (largely untested) app I threw 
together on reading your message, using a technique Paul Emsley tried out 
yesterday. It uses a 6-point spline scaling curve to acheive effects of both 
scaling and bulk solvent correction. I can send you an x86 linux binary if 
you want to try it.

Usage:
  diffmap my.mtz '/*/*/[FP,SIGFP]' '/*/*/[FCAL,PHICAL]'
Outputs:
  diff.map

On Friday 07 March 2003 1:54 am, Mark DePristo wrote:
> I've been struggling to create reasonable fo-fc maps with CCP4.  The
> fo-fc maps created directly with xfit in xtalview are very nice,
> clearly highlighting poorly placed sidechain atoms on the protein
> surface.  However, the fo-fc maps calculated with CCP4's fft from the
> equivalent mtz FP, FC, and PHIC are poor and impossible to interpret.
> Interestingly, the xfit fo-fc maps are very similar to the CCP4 maps
> when xfit's fft the bulk solvent correction is disabled.

-- 
Structural Biology Laboratory, University of York, York, YO10 5YW, UK
// Clipper test application
/* (C) 2003 Kevin Cowtan */


#include <clipper/clipper.h>
#include <clipper/clipper-mtz.h>


using clipper::data32::F_sigF;
using clipper::data32::F_phi;


int main(int argc, char** argv)
{
  const int n_scl_param = 6;

  clipper::MTZfile mtzfile;
  clipper::MTZcrystal xtal;
  clipper::MTZdataset dset;

  // make data objects
  clipper::HKL_info wrk_hkl;
  clipper::HKL_data<F_sigF> wrk_fo( wrk_hkl );
  clipper::HKL_data<F_phi> wrk_fc( wrk_hkl );

  // read data
  mtzfile.open_read( argv[1] );
  mtzfile.import_hkl_info( wrk_hkl );
  mtzfile.import_hkl_data( wrk_fo, dset, xtal, argv[2] );
  mtzfile.import_hkl_data( wrk_fc, dset, xtal, argv[3] );
  mtzfile.close_read();

  // scale and difference data
  std::vector<clipper::ftype> params( n_scl_param, 1.0 );
  clipper::BasisFn_spline wrk_basis( wrk_hkl, n_scl_param, 2.0 );
  clipper::TargetFn_scaleF1F2<F_phi,F_sigF> wrk_target( wrk_fc, wrk_fo );
  clipper::ResolutionFn wrk_scale( wrk_hkl, wrk_basis, wrk_target, params );
  clipper::HKL_info::HKL_reference_index ih;
  for ( ih = wrk_hkl.first(); !ih.last(); ih.next() )
    if ( !wrk_fc[ih].missing() ) {
      wrk_fc[ih].scale( sqrt( wrk_scale.f(ih) ) );  // scale
      wrk_fc[ih].f() -= wrk_fo[ih].f();             // difference
    }

  // calculate map
  clipper::Grid_sampling wrk_grid( wrk_hkl.spacegroup(), wrk_hkl.cell(),
				   wrk_hkl.resolution() );
  clipper::Xmap<float> wrk_map( wrk_hkl.spacegroup(), wrk_hkl.cell(),
				wrk_grid );
  wrk_map.fft_from( wrk_fc );

  // output the map
  clipper::MAPfile mapfile;
  mapfile.open_write( "diff.map" );
  mapfile.export_xmap( wrk_map );
  mapfile.close_write();
}