FormFactorTable.h

Go to the documentation of this file.

/**
 *  \file IMP/saxs/FormFactorTable.h   \brief A class for computation of
 * atomic and residue level form factors for SAXS calculations
 *
 *  Copyright 2007-2014 IMP Inventors. All rights reserved.
 *
 */

#ifndef IMPSAXS_FORM_FACTOR_TABLE_H
#define IMPSAXS_FORM_FACTOR_TABLE_H

#include <IMP/saxs/saxs_config.h>

#include <IMP/kernel/Particle.h>
#include <IMP/base_types.h>
#include <IMP/atom/Residue.h>
#include <IMP/atom/Atom.h>
#include <IMP/atom/element.h>
#include <IMP/algebra/utility.h>

#include <iostream>
#include <vector>

IMPSAXS_BEGIN_NAMESPACE

//! type of the form factors for profile calculations
/*
 ALL_ATOMS - all atoms including hydrogens
 HEAVY_ATOMS - no hydrogens, all other atoms included
 CA_ATOMS - residue level, residue represented by CA
*/
enum FormFactorType {
  ALL_ATOMS,
  HEAVY_ATOMS,
  CA_ATOMS
};

/**
   class that deals with form factor computation
   two form factors are supported:
   (i) zero form factors for faster approximated calculations
   (ii) full form factors for slower accurate calculations

   Each form factor can be divided into two parts: vacuum and dummy.
   dummy is an approximated excluded volume (solvent) form factor.
   The approximation is done using Fraser, MacRae and Suzuki (1978) model.
*/
class IMPSAXSEXPORT FormFactorTable {
 public:
  //! default constructor
  FormFactorTable();

  //! constructor with form factor table file (required for full form factors)
  FormFactorTable(const String& table_name, Float min_q, Float max_q,
                  Float delta_q);

  // 1. Zero form factors

  //! get f(0), ie q=0 for real space profile calculation
  Float get_form_factor(kernel::Particle* p,
                        FormFactorType ff_type = HEAVY_ATOMS) const;

  //! f(0) in vacuum
  Float get_vacuum_form_factor(kernel::Particle* p,
                               FormFactorType ff_type = HEAVY_ATOMS) const;

  //! f(0) for solvent
  Float get_dummy_form_factor(kernel::Particle* p,
                              FormFactorType ff_type = HEAVY_ATOMS) const;

  //! f(0) for water
  Float get_water_form_factor() const { return zero_form_factors_[OH2]; }

  //! f(0) for water in vacuum
  Float get_vacuum_water_form_factor() const {
    return vacuum_zero_form_factors_[OH2];
  }

  //! f(0) for water (solvent)
  Float get_dummy_water_form_factor() const {
    return dummy_zero_form_factors_[OH2];
  }

  // 2. Full form factors

  //! full form factor for reciprocal space profile calculation
  const Floats& get_form_factors(kernel::Particle* p,
                                 FormFactorType ff_type = HEAVY_ATOMS) const;

  //! for reciprocal space profile calculation
  const Floats& get_vacuum_form_factors(kernel::Particle* p,
                                        FormFactorType ff_type =
                                            HEAVY_ATOMS) const;

  //! for reciprocal space profile calculation
  const Floats& get_dummy_form_factors(kernel::Particle* p,
                                       FormFactorType ff_type =
                                           HEAVY_ATOMS) const;

  //! full water form factor
  const Floats& get_water_form_factors() const { return form_factors_[OH2]; }

  //! full water vacuum form factor
  const Floats& get_water_vacuum_form_factors() const {
    return vacuum_form_factors_[OH2];
  }

  //! full water dummy form factor
  const Floats& get_water_dummy_form_factors() const {
    return dummy_form_factors_[OH2];
  }

  //! radius
  Float get_radius(kernel::Particle* p,
                   FormFactorType ff_type = HEAVY_ATOMS) const;

  //! volume
  Float get_volume(kernel::Particle* p,
                   FormFactorType ff_type = HEAVY_ATOMS) const;

  //! print tables
  void show(std::ostream& out = std::cout, std::string prefix = "") const;

  // electron density of solvent - default=0.334 e/A^3 (H2O)
  static Float rho_;

 private:
  // atom types for heavy atoms according to the number of hydrogens
  // connected to them
  // ALL_ATOM_SIZE is number of types needed for all atom representation
  // this indexing is used in form_factors arrays
  enum FormFactorAtomType {
    H,
    He,
    Li,
    Be,
    B,
    C,
    N,
    O,
    F,
    Ne,  // periodic table, lines 1-2 (10)
    Na,
    Mg,
    Al,
    Si,
    P,
    S,
    Cl,
    Ar,  // line 3 (8)
    K,
    Ca,
    Cr,
    Mn,
    Fe,
    Co,
    Ni,
    Cu,
    Zn,
    Se,
    Br,  // line 4 (11)
    I,
    Ir,
    Pt,
    Au,
    Hg,
    ALL_ATOM_SIZE = 34,
    CH = 34,
    CH2 = 35,
    CH3 = 36,
    NH = 37,
    NH2 = 38,
    NH3 = 39,
    OH = 40,
    OH2 = 41,
    SH = 42,
    HEAVY_ATOM_SIZE = 43,
    UNK = 44
  };

  // map between atom element and FormFactorAtomType
  static std::map<atom::Element, FormFactorAtomType> element_ff_type_map_;

  struct FormFactor {
    FormFactor() {}
    FormFactor(double ff, double vacuum_ff, double dummy_ff)
        : ff_(ff), vacuum_ff_(vacuum_ff), dummy_ff_(dummy_ff) {}
    double ff_, vacuum_ff_, dummy_ff_;
  };

  // map between residue type and residue level form factors
  static std::map<atom::ResidueType, FormFactor> residue_type_form_factor_map_;

  // form factors for q=0, the order as in the FormFactorAtomType enum
  static Float zero_form_factors_[];

  static Float vacuum_zero_form_factors_[];
  // those represent excluded volume
  static Float dummy_zero_form_factors_[];

  // a key for storing zero form factor in kernel::Particle as attribute
  static IntKey form_factor_type_key_;

  // class for storing form factors solvation table
  class AtomFactorCoefficients {
   public:
    String atom_type_;
    Float a_[5];
    Float b_[5];
    Float c_;
    Float excl_vol_;
  };

  // read entry
  friend std::istream& operator>>(
      std::istream& s, AtomFactorCoefficients& atom_factor_coefficients);

  // write entry
  friend std::ostream& operator<<(
      std::ostream& s, const AtomFactorCoefficients& atom_factor_coefficients);

 private:
  int read_form_factor_table(const String& table_name);

  void init_element_form_factor_map();

  void init_residue_type_form_factor_map();

  void compute_form_factors_all_atoms();

  void compute_form_factors_heavy_atoms();

  float get_form_factor(atom::ResidueType rt) const;

  float get_vacuum_form_factor(atom::ResidueType rt) const;

  float get_dummy_form_factor(atom::ResidueType rt) const;

  FormFactorAtomType get_form_factor_atom_type(atom::Element e) const;

  FormFactorAtomType get_form_factor_atom_type(kernel::Particle* p,
                                               FormFactorType ff_type) const;

  FormFactorAtomType get_carbon_atom_type(const atom::AtomType& atom_type,
                                          const atom::ResidueType& residue_type)
      const;

  FormFactorAtomType get_nitrogen_atom_type(
      const atom::AtomType& atom_type,
      const atom::ResidueType& residue_type) const;

  FormFactorAtomType get_oxygen_atom_type(const atom::AtomType& atom_type,
                                          const atom::ResidueType& residue_type)
      const;

  FormFactorAtomType get_sulfur_atom_type(const atom::AtomType& atom_type,
                                          const atom::ResidueType& residue_type)
      const;

 private:
  // read from lib file
  std::vector<AtomFactorCoefficients> form_factors_coefficients_;

  // table of full form factors for 14 atom types
  std::vector<Floats> form_factors_;

  // vacuum full form factors for 14 atom types
  std::vector<Floats> vacuum_form_factors_;

  // dummy full form factors for 14 atom types
  std::vector<Floats> dummy_form_factors_;

  // min/max q and sampling resolution for form factor computation
  Float min_q_, max_q_, delta_q_;

  base::WarningContext warn_context_;
};

/** Get the default table.*/
IMPSAXSEXPORT FormFactorTable* default_form_factor_table();

IMPSAXS_END_NAMESPACE

#endif /* IMPSAXS_FORM_FACTOR_TABLE_H */