Index: kernel/include/IMP/score_states/NonbondedListScoreState.h =================================================================== --- kernel/include/IMP/score_states/NonbondedListScoreState.h (revision 593) +++ kernel/include/IMP/score_states/NonbondedListScoreState.h (working copy) @@ -210,7 +210,7 @@ */ void set_slack(float slack) { IMP_check(slack>= 0, "Slack must be nonnegative", - ValueException("Negative slack")); + ValueException); slack_=slack; } @@ -236,13 +236,13 @@ //! Iterates through the pairs of non-bonded particles NonbondedIterator nonbonded_begin() const { IMP_check(get_nbl_is_valid(), "Must call update first", - ValueException("Must call update first")); + ValueException); return NonbondedIterator(boxes_overlap_object(cutoff_), nbl_.begin(), nbl_.end()); } NonbondedIterator nonbonded_end() const { IMP_check(get_nbl_is_valid(), "Must call update first", - ValueException("Must call update first")); + ValueException); return NonbondedIterator(boxes_overlap_object(cutoff_), nbl_.end(), nbl_.end()); } @@ -253,7 +253,7 @@ unsigned int number_of_nonbonded() const { IMP_check(get_nbl_is_valid(), "Must call update first", - ValueException("Must call update first")); + ValueException); return nbl_.size(); } Index: kernel/include/IMP/Restraint.h =================================================================== --- kernel/include/IMP/Restraint.h (revision 593) +++ kernel/include/IMP/Restraint.h (working copy) @@ -41,6 +41,9 @@ constructor and should provide methods so that the set of particles can be modified after construction. + A restraint can be added to the model multiple times or to multiple + restraint sets in the same model. + \note When logging is VERBOSE, restraints should print enough information in evaluate to reproduce the the entire flow of data in evaluate. When logging is TERSE the restraint should print out only a constant number of Index: kernel/include/IMP/optimizers/MonteCarlo.h =================================================================== --- kernel/include/IMP/optimizers/MonteCarlo.h (revision 593) +++ kernel/include/IMP/optimizers/MonteCarlo.h (working copy) @@ -74,7 +74,7 @@ */ void set_move_probability(Float p) { IMP_check(p > 0 && p <= 1, "Not a valid probability", - ValueException("Not a probability")); + ValueException); probability_=p; } Index: kernel/include/IMP/optimizers/movers/NormalMover.h =================================================================== --- kernel/include/IMP/optimizers/movers/NormalMover.h (revision 593) +++ kernel/include/IMP/optimizers/movers/NormalMover.h (working copy) @@ -30,7 +30,7 @@ void set_sigma(Float sigma) { IMP_check(sigma > 0, "Sigma must be positive", - ValueException("Negative sigma")); + ValueException); stddev_=sigma; } Float get_sigma() const { Index: kernel/include/IMP/optimizers/movers/BallMover.h =================================================================== --- kernel/include/IMP/optimizers/movers/BallMover.h (revision 593) +++ kernel/include/IMP/optimizers/movers/BallMover.h (working copy) @@ -35,7 +35,7 @@ //! void set_radius(Float radius) { IMP_check(radius > 0, "The radius must be positive", - ValueException("Negative radius")); + ValueException); radius_=radius; } //! Index: kernel/include/IMP/optimizers/BrownianDynamics.h =================================================================== --- kernel/include/IMP/optimizers/BrownianDynamics.h (revision 593) +++ kernel/include/IMP/optimizers/BrownianDynamics.h (working copy) @@ -40,17 +40,20 @@ BrownianDynamics(FloatKey dkey= FloatKey("D")); virtual ~BrownianDynamics(); - IMP_OPTIMIZER(internal::kernel_version_info) + IMP_OPTIMIZER(internal::kernel_version_info); + //! Simulate until the given time in fs + void simulate(float time_in_fs); + //! Set the max step that should be allowed in angstroms /** The timestep is shrunk if any particles are moved more than this */ void set_max_step(Float a) { - set_max_step(internal::Angstrom(a)); + set_max_step(unit::Angstrom(a)); } //! Set time step in fs void set_time_step(Float t) { - set_time_step(internal::FemtoSecond(t)); + set_time_step(unit::Femtosecond(t)); } //! Time step in fs @@ -64,33 +67,37 @@ } - void set_temperature(Float t) { set_temperature(internal::Kelvin(t)); } + void set_temperature(Float t) { set_temperature(unit::Kelvin(t)); } Float get_temperature() const { return get_temperature_units().get_value(); } FloatKey get_d_key() const { return dkey_; } - //! Estimate the diffusion coeffient from the mass in kD - /** This method does not do anything fancy, but is OK for large - globular proteins. + //! Estimate the radius of a protein from the mass + /** Proteins are assumed to be spherical. The density is estimated + using the formula from + - Note that it depends on temperature and so can't be static. + The formula is: + density= 1.410+ 0.145 exp(-M/13) g/cm^3 */ - Float estimate_radius_from_mass(Float mass_in_kd) const { - return estimate_radius_from_mass_units(mass_in_kd).get_value(); + static Float estimate_radius_from_mass(Float mass_in_kd) { + return + estimate_radius_from_mass_units(unit::Kilodalton(mass_in_kd)).get_value(); } //! Return the expected distance moved for a particle with a given D /** The units on D are cm^2/sec and the return has units of Angstroms. */ Float compute_sigma_from_D(Float D) const { - return compute_sigma_from_D(internal::Cm2PerSecond(D)).get_value(); + unit::SquareCentimeterPerSecond du(D); + return compute_sigma_from_D(du).get_value(); } //! Estimate the diffusion coefficient from the radius in Angstroms /** This depends on the temperature. */ Float estimate_D_from_radius(Float radius_in_angstroms) const { - internal::Angstrom r(radius_in_angstroms); + unit::Angstrom r(radius_in_angstroms); return estimate_D_from_radius(r).get_value(); } @@ -98,7 +105,8 @@ /** This value is in KCal/A/mol */ Float get_force_scale_from_D(Float D) const { - return get_force_scale_from_D(internal::Cm2PerSecond(D)).get_value(); + unit::SquareCentimeterPerSecond du(D); + return get_force_scale_from_D(du).get_value(); } //! Get the current time in femtoseconds @@ -108,7 +116,7 @@ //! Set the current time in femtoseconds void set_current_time(Float fs) { - set_current_time( internal::FemtoSecond(fs)); + set_current_time( unit::Femtosecond(fs)); } //! Return a histogram of timesteps @@ -122,55 +130,62 @@ IMP_LIST(private, Particle, particle, Particle*); protected: + //! Perform a single dynamics step. // \return true if the initial step size was OK bool step(); + void take_step(); + /** Propose a single step, this will be accepted if all moves are small enough. \return true if it should be accepted. */ bool propose_step(std::vector &proposal); - internal::FemtoJoule kt() const; + unit::Femtojoule kt() const; void setup_particles(); - void set_max_step(internal::Angstrom a) { + void set_max_step(unit::Angstrom a) { max_change_= a; } - void set_time_step(internal::FemtoSecond t); + void set_time_step(unit::Femtosecond t); - internal::FemtoSecond get_time_step_units() const {return max_dt_;} + unit::Femtosecond get_time_step_units() const {return max_dt_;} - internal::FemtoSecond get_last_time_step_units() const {return cur_dt_;} + unit::Femtosecond get_last_time_step_units() const {return cur_dt_;} - void set_temperature(internal::Kelvin t) { T_=t;} + void set_temperature(unit::Kelvin t) { T_=t;} - internal::Kelvin get_temperature_units() const { return T_;} + unit::Kelvin get_temperature_units() const { return T_;} - internal::Angstrom estimate_radius_from_mass_units(Float mass_in_kd) const; + static unit::Angstrom + estimate_radius_from_mass_units(unit::Kilodalton mass_in_kd); - internal::Angstrom compute_sigma_from_D(internal::Cm2PerSecond D) const; + unit::Angstrom + compute_sigma_from_D(unit::SquareCentimeterPerSecond D) const; - internal::Cm2PerSecond - estimate_D_from_radius(internal::Angstrom radius) const; + unit::SquareCentimeterPerSecond + estimate_D_from_radius(unit::Angstrom radius) const; - internal::KCalPerAMol get_force_scale_from_D(internal::Cm2PerSecond D) const; + unit::KilocaloriePerAngstromPerMol + get_force_scale_from_D(unit::SquareCentimeterPerSecond D) const; - internal::FemtoSecond get_current_time_units() const { + unit::Femtosecond get_current_time_units() const { return cur_time_; } //! Set the current time in femtoseconds - void set_current_time(internal::FemtoSecond fs) { + void set_current_time(unit::Femtosecond fs) { cur_time_= fs; } - internal::Angstrom max_change_; - internal::FemtoSecond max_dt_, cur_dt_, cur_time_; - internal::Kelvin T_; + unsigned int num_const_dt_; + unit::Angstrom max_change_; + unit::Femtosecond max_dt_, cur_dt_, cur_time_; + unit::Kelvin T_; FloatKey dkey_; Index: kernel/include/IMP/optimizers/Mover.h =================================================================== --- kernel/include/IMP/optimizers/Mover.h (revision 593) +++ kernel/include/IMP/optimizers/Mover.h (working copy) @@ -51,7 +51,7 @@ moves. This number can be either used to scale a continuous move or affect the probability of a discrete move. */ - virtual void propose_move(float size)=0; + virtual void propose_move(float probability)=0; //! set whether the proposed modification is accepted /** \note Accepting should not change the Particles at all. */ Index: kernel/include/IMP/decorators/XYZDecorator.h =================================================================== --- kernel/include/IMP/decorators/XYZDecorator.h (revision 593) +++ kernel/include/IMP/decorators/XYZDecorator.h (working copy) @@ -117,7 +117,7 @@ protected: static FloatKey get_coordinate_key(unsigned int i) { IMP_check(i <3, "Out of range coordinate", - IndexException("Out of range coordinate")); + IndexException); return key_[i]; } }; Index: kernel/include/IMP/decorators/MolecularHierarchyDecorator.h =================================================================== --- kernel/include/IMP/decorators/MolecularHierarchyDecorator.h (revision 593) +++ kernel/include/IMP/decorators/MolecularHierarchyDecorator.h (working copy) @@ -93,12 +93,11 @@ unsigned int add_child(This o) { IMP_check(get_type() > o.get_type(), "Parent type must subsume child type", - InvalidStateException("Type of hierarchy parent must be "\ - "higher than type of child")); + InvalidStateException); IMP_check(get_type() != UNKNOWN, "Parent must have known type", - InvalidStateException("Hierarchy parent must have known type")); + InvalidStateException); IMP_check(o.get_type() != UNKNOWN, "Child must have known type", - InvalidStateException("Hierarchy child must have known type")); + InvalidStateException); return P::add_child(o); } @@ -109,11 +108,11 @@ void add_child_at(This o, unsigned int i) { IMP_check(get_type() > o.get_type(), "Parent type must subsume child type", - InvalidStateException("")); + InvalidStateException); IMP_check(get_type() != UNKNOWN, "Parent must have known type", - InvalidStateException("")); + InvalidStateException); IMP_check(o.get_type() != UNKNOWN, "Child must have known type", - InvalidStateException("")); + InvalidStateException); P::add_child_at(o, i); } Index: kernel/include/IMP/internal/ExponentialNumber.h =================================================================== --- kernel/include/IMP/internal/ExponentialNumber.h (revision 0) +++ kernel/include/IMP/internal/ExponentialNumber.h (revision 0) @@ -0,0 +1,129 @@ +/** + * \file ExponentialNumber.h + * \brief Classes to add compile time exponents to numbers. + * + * Copyright 2007-8 Sali Lab. All rights reserved. + * + */ + +#ifndef __IMP_EXPONENTIAL_NUMBER_H +#define __IMP_EXPONENTIAL_NUMBER_H + +#include "../macros.h" + +#include +#include +#include +#include + +namespace IMP +{ + +namespace internal +{ + +namespace unit +{ + +/** This helper class implements a floating point number + with an extra exponent. The actual value is the number + times 10^EXP. I needed this since I was running out of bits + in the exponent when converting between vastly different units. + It also means that you can have a type for both nanometer + and meter without much pain. Since the exponent is stored + as a template argument, there shouldn't be any runtime overhead. + */ +template +class ExponentialNumber +{ + double v_; + + typedef ExponentialNumber This; + + template + void copy_from(ExponentialNumber o) { + const int diff = OEXP-EXP; + v_=o.v_; + float factor= std::pow(10.0, static_cast(diff)); + v_*=factor; + } + bool is_default() const { + return v_== std::numeric_limits::infinity(); + } + + template + friend class ExponentialNumber; +public: + + + ExponentialNumber(): v_(std::numeric_limits::infinity()) {} + template + ExponentialNumber(ExponentialNumber o) { + copy_from(o); + } + explicit ExponentialNumber(double d): v_(d){ + } + template + This &operator=(ExponentialNumber o) { + copy_from(o); + return *this; + } + + //! Compute the value with the exponent + double get_normalized_value() const { + double ret=v_; + return ret*std::pow(10.0, EXP); + } + //! Get the stored value ignoring the exponent + double get_value() const { + return v_; + } + template + ExponentialNumber operator+(ExponentialNumber o) const { + return This(v_+ This(o).v_); + } + template + ExponentialNumber operator-(ExponentialNumber o) const { + return This(v_- This(o).v_); + } + + template + ExponentialNumber operator/(ExponentialNumber o) const { + return ExponentialNumber(v_/o.v_); + } + + template + ExponentialNumber operator*(ExponentialNumber o) const { + return ExponentialNumber(v_*o.v_); + } + + This operator-() const { + return This(-v_); + } + + void show(std::ostream &out) const { + std::ios::fmtflags of = out.flags(); + out << setiosflags(std::ios::fixed) << v_; + if (EXP != 0) out << "e" << EXP; + out.flags(of); + /** \todo I should restore the io flags. + */ + } + + IMP_COMPARISONS_1(v_); +}; + +template +std::ostream &operator<<(std::ostream &out, ExponentialNumber o) +{ + o.show(out); + return out; +} + +} // namespace unit + +} // namespace internal + +} // namespace IMP + +#endif /* __IMP_EXPONENTIAL_NUMBER_H */ Index: kernel/include/IMP/internal/Vector.h =================================================================== --- kernel/include/IMP/internal/Vector.h (revision 593) +++ kernel/include/IMP/internal/Vector.h (working copy) @@ -48,14 +48,14 @@ const D& operator[](unsigned int i) const { IMP_check(i < P::size(), "Index " << i << " out of range", - IndexException("")); + IndexException); return P::operator[](i); } D& operator[](unsigned int i) { IMP_check(i < P::size(), "Index " << i << " out of range", - IndexException("")); + IndexException); return P::operator[](i); } Index: kernel/include/IMP/internal/units.h =================================================================== --- kernel/include/IMP/internal/units.h (revision 593) +++ kernel/include/IMP/internal/units.h (working copy) @@ -8,11 +8,12 @@ #ifndef __IMP_UNITS_H #define __IMP_UNITS_H -#include "../utility.h" +#include "Unit.h" +#include "../base_types.h" + #include #include -#include #include namespace IMP @@ -21,478 +22,211 @@ namespace internal { -/** This helper class implements a floating point number - with an extra exponent. The actual value is the number - times 10^EXP. I needed this since I was running out of bits - in the exponent when converting between vastly different units. - It also means that you can have a type for both nanometer - and meter without much pain. Since the exponent is stored - as a template argument, there shouldn't be any runtime overhead. - */ -template -class ExponentialNumber -{ - typedef ExponentialNumber This; +// want it for conversions +extern IMPDLLEXPORT const unit::ExponentialNumber<3> JOULES_PER_KILOCALORIE; +extern IMPDLLEXPORT const unit::ExponentialNumber<23> NA; - template - void copy_from(ExponentialNumber o) { - const int diff = OEXP-EXP; - v_=o.v_; - float factor= std::pow(10.0, static_cast(diff)); - v_*=factor; - IMP_assert(std::abs(get_normalized_value() -o.get_normalized_value()) - <= .09 * std::abs(std::max(get_normalized_value(), - o.get_normalized_value())), - "Mismatch after scaling"); - } - bool is_default() const { - return v_== std::numeric_limits::infinity(); - } - // should be private - double v_; - template - friend class ExponentialNumber; -public: +namespace unit +{ - ExponentialNumber(): v_(std::numeric_limits::infinity()) {} - template - ExponentialNumber(ExponentialNumber o) { - copy_from(o); - } - ExponentialNumber(double d): v_(d){ - } - template - This &operator=(ExponentialNumber o) { - copy_from(o); - return *this; - } +namespace internal +{ - //! Compute the value with the exponent - double get_normalized_value() const { - double ret=v_; - return ret*std::pow(10.0, EXP); - } - //! Get the stored value ignoring the exponent - double get_value() const { - return v_; - } - template - ExponentialNumber operator+(ExponentialNumber o) const { - return This(v_+ This(o).v_); - } - template - ExponentialNumber operator-(ExponentialNumber o) const { - return This(v_- This(o).v_); - } +struct AtomsPerMol {}; - template - ExponentialNumber operator/(ExponentialNumber o) const { - return ExponentialNumber(v_/o.v_); - } - template - ExponentialNumber operator*(ExponentialNumber o) const { - return ExponentialNumber(v_*o.v_); - } +struct MDEnergyTag; - This operator-() const { - return This(-v_); - } +template <> +inline std::string get_unit_name(int o) { + std::string os[]= {"Cal/Mol"}; + return os[0]; +} - void show(std::ostream &out) const { - std::ios::fmtflags of = out.flags(); - out << setiosflags(std::ios::fixed) << v_; - if (EXP != 0) out << "e" << EXP; - out.flags(of); - /** \todo I should restore the io flags. - */ - } - IMP_COMPARISONS_1(v_); -}; +struct MDDerivativeTag; -template -std::ostream &operator<<(std::ostream &out, ExponentialNumber o) -{ - o.show(out); - return out; +template <> +inline std::string get_unit_name(int o) { + std::string os[]= {"Cal/(A Mol)"}; + return os[0]; } -// helper function generate unit string names -inline void handle_unit(std::ostringstream &out, std::string name, int m) -{ - if (m>0) { - out << name; - if (m>1) { - out << "^" << m<< " "; - } - } -} +struct MKSTag{}; -//! A class to represent a quantity in the MKS system -/** The correct units are preserved with arithmetic operations - (and only operations with compatible units are allowed). The EXP - parameter is there to allow representation of quantities in cm or - other prefixed units. +template <> +inline std::string get_unit_name(int o) { + std::string os[]= {"kg", "m", "s", "k", "Cal"}; + return os[o]; +} - The template parameters are - - EXP the shift in base 10 - - M the power for meters - - KG the power for kilograms - - K the power for Kelvin - - S the power for seconds - */ -template -class MKSUnit -{ - typedef MKSUnit This; - typedef ExponentialNumber V; - template - friend class MKSUnit; +struct DaltonTag; - // for debugging since the debugger doesn't always show - // template arguments - int get_exponent() const {return EXP;} - int get_meters() const {return M;} - int get_kilograms() const {return KG;} - int get_kelvin() const {return K;} - int get_seconds() const {return S;} +template <> +inline std::string get_unit_name(int o) { + std::string os[]= {"Da"}; + return os[0]; +} +} // ns unit::internal - V v_; -public: - explicit MKSUnit(V v): v_(v) {} - MKSUnit(){} - MKSUnit(double v): v_(v) { - if (v < -std::numeric_limits::max()) { - std::cout << get_exponent() << get_meters() - << get_kilograms() << get_kelvin() << get_seconds(); - } - } - MKSUnit(float v): v_(static_cast(v)) {} - MKSUnit(int v): v_(static_cast(v)) {} - template - MKSUnit(MKSUnit o): v_(o.v_) {} +typedef boost::mpl::vector_c Mass; +typedef boost::mpl::vector_c Length; +typedef boost::mpl::vector_c Time; +typedef boost::mpl::vector_c Temperature; +typedef boost::mpl::vector_c HeatEnergy; - //! MKS units can be compared if their types match - template - bool operator<(MKSUnit o) const { - return v_ < o.v_; - } +typedef boost::mpl::vector_c Energy; +typedef boost::mpl::vector_c Force; +typedef boost::mpl::vector_c Pressure; +typedef boost::mpl::vector_c HeatEnergyDerivative; - //! MKS units can be compared if their types match - template - bool operator<=(MKSUnit o) const { - return v_ <= o.v_; - } - //! MKS units can be compared if their types match - template - bool operator>=(MKSUnit o) const { - return v_ >= o.v_; - } - //! only works with matching types - template - bool operator>(MKSUnit o) const { - return v_ > o.v_; - } - void show(std::ostream &out) const { - out << v_ << get_name(); - } +typedef Unit Meter; +typedef Unit Centimeter; +typedef Unit Kilogram; +typedef Unit Gram; +typedef Unit Second; +typedef Unit Joule; +typedef Unit Kelvin; +typedef Shift::type Angstrom; +typedef Multiply::type SquareAngstrom; +typedef Shift::type Femtojoule; +typedef Shift::type Picojoule; +typedef Divide::type JoulePerKelvin; +typedef Multiply::type SquareMeter; +typedef Divide::type SquareMeterPerSecond; +typedef Shift::type SquareCentimeterPerSecond; +typedef Unit Femtosecond; +typedef Unit Femtonewton; +typedef Unit Piconewton; +typedef Unit Femtogram; +typedef Divide::type +PiconewtonPerAngstrom; +typedef Unit Pascal; +typedef Unit Kilocalorie; +typedef Unit YoctoKilocalorie; +typedef Divide::type KilocaloriePerMeter; +typedef Divide::type KilocaloriePerAngstrom; +// Calorie is ambiguous given our naming convention +typedef Shift::type YoctoKilocaloriePerMeter; +typedef Shift::type YoctoKilocaloriePerAngstrom; +typedef Multiply::type, + Centimeter>::type CubicCentimeter; +typedef Divide::type GramPerCubicCentimeter; - //! Get the current value - /** \note The value returned is the value before it is multiplied - by the appropriate power of 10. That means, 1 Angstrom returns 1, - as does 1 Meter. - */ - double get_value() const { - return v_.get_value(); - } - //! Get the value with the proper exponent - /** 1 Angstrom returns 1e-10. - */ - double get_normalized_value() const { - return v_.get_normalized_value(); - } - V get_exponential_value() const { - return v_; - } - template - MKSUnit - operator*(MKSUnit o) const { - return MKSUnit(v_*o.v_); - } +template +typename Exchange, + Kilocalorie, Joule, 4>::type +convert_Cal_to_J(Unit i) { + typedef typename Exchange, + Kilocalorie, Joule, 4>::type Return; + return Return(i.get_exponential_value() + * JOULES_PER_KILOCALORIE); +} - template - MKSUnit - operator/(MKSUnit o) const { - return MKSUnit(v_/o.v_); - } +template +typename Exchange, + Joule, Kilocalorie, -3>::type +convert_J_to_Cal(Unit i) { + typedef typename Exchange, + Joule, Kilocalorie, -3>::type Return; + std::cout << i.get_exponential_value() + / JOULES_PER_KILOCALORIE << std::endl; + return Return(i.get_exponential_value() + / JOULES_PER_KILOCALORIE); +} - This operator+(This o) const { - return This(v_+o.v_); - } - This operator-(This o) const { - return This(v_-o.v_); - } - This operator-() const { - return This(-v_); - } +//! Marker constant to handle coversions +extern const IMPDLLEXPORT internal::AtomsPerMol ATOMS_PER_MOL; - //! Get a string version of the name of the units - std::string get_name() const { - int m=M; - int kg=KG; - int k=K; - int s=S; - int j=0; - while (kg >=1 && m > 1 && s<=-2) { - m-=2; - kg-=1; - s+=2; - j+=1; - } - std::ostringstream name; - handle_unit(name, "m", m); - handle_unit(name, "kg", kg); - handle_unit(name, "K", k); - handle_unit(name, "s", s); - handle_unit(name, "J", j); - if (m < 0 || kg < 0 || k < 0 || s < 0 || j < 0) { - name << "/"; - handle_unit(name, "m", -m); - handle_unit(name, "kg", -kg); - handle_unit(name, "K", -k); - handle_unit(name, "s", -s); - handle_unit(name, "J", -j); - } - return name.str(); - } -}; +typedef Unit KilocaloriePerMol; +typedef Unit +KilocaloriePerAngstromPerMol; -template -MKSUnit -sqrt(MKSUnit o) -{ - return MKSUnit(std::sqrt(o.get_value())); +inline Unit +operator/(KilocaloriePerMol k, + internal::AtomsPerMol) { + return Unit(k.get_value()*1.661); } -template -MKSUnit -square(MKSUnit o) -{ - return MKSUnit(IMP::square(o.get_value())); +inline Unit +operator/(KilocaloriePerAngstromPerMol k, + internal::AtomsPerMol) { + return Unit(k.get_value()*1.661); } -//! They have to be the same EXP to make the return type make sense -template -MKSUnit -min(MKSUnit a, MKSUnit b) { - if (a +KilocaloriePerMol operator*(Unit k, + internal::AtomsPerMol) { + return KilocaloriePerMol(k.get_exponential_value() + *NA); } -//! They have to be the same EXP to make the return type make sense -template -MKSUnit -max(MKSUnit a, MKSUnit b) { - if (a>b) return a; - else return b; +template +KilocaloriePerAngstromPerMol +operator*(Unit k, + internal::AtomsPerMol) { + return KilocaloriePerAngstromPerMol(k.get_exponential_value() + *NA*ExponentialNumber<-10>(1)); } -typedef MKSUnit<-10, 1, 0, 0, 0> Angstrom; -typedef MKSUnit<-20, 2, 0, 0, 0> SquaredAngstrom; -typedef MKSUnit<0, 2, 1, 0, -2> Joule; -typedef MKSUnit<-15, 2, 1, 0, -2> FemtoJoule; -typedef MKSUnit<-12, 2, 1, 0, -2> PicoJoule; -typedef MKSUnit<0, 0, 0, 1, 0> Kelvin; -typedef MKSUnit<0,2,1,-1,-2> JoulePerKelvin; -typedef MKSUnit<-15, 0,0,0,1> FemtoSecond; -typedef MKSUnit<-4,2,0,0,-1> Cm2PerSecond; -typedef MKSUnit<-15,1,1,0,-2> FemtoNewton; -typedef MKSUnit<-12,1,1,0,-2> PicoNewton; -typedef MKSUnit<-2, 0, 1, 0, -2> PicoNewtonPerAngstrom; -typedef MKSUnit<0,1,0,0,0> Meter; -typedef MKSUnit<0,0,0,0,0> Scalar; - -template -std::ostream &operator<<(std::ostream &out, MKSUnit u) -{ - u.show(out); - return out; +template +KilocaloriePerMol operator*(internal::AtomsPerMol, + Unit k) { + return operator*(k, ATOMS_PER_MOL); } -//! Represent a value with units in kcal and angstroms per mol -/** This class is just for representing energies and derivatives - for MD-type computations and so is a bit more limited than - the MKS unit type. - */ -template -class KCalPerMolUnit -{ - typedef KCalPerMolUnit This; - template - friend class KCalPerMolUnit; - typedef ExponentialNumber V; - V v_; - int get_exponent() const {return EXP;} - int get_kilocalories() const {return KCal;} - int get_angstroms() const {return A;} -public: - explicit KCalPerMolUnit(ExponentialNumber e): v_(e) {} - KCalPerMolUnit(float f): v_(f) { - if (0 < -std::numeric_limits::max()) { - std::cout << get_exponent() << get_angstroms() - << get_kilocalories(); - } - } +template +KilocaloriePerAngstromPerMol +operator*(internal::AtomsPerMol, + Unit + /*typename YoctoKilocaloriePerAngstrom::Units>*/ k) { + return operator*(k, ATOMS_PER_MOL); +} - template - KCalPerMolUnit(KCalPerMolUnit o):v_(o.v_) { - } - template - KCalPerMolUnit(MKSUnit o): v_(o.v_) {} - //! Get the value ignoring the exponent - double get_value() const {return v_.get_value();} - This operator-() const { - return This(-v_); - } +// define Daltons - //! Get the value ignoring the exponent - V get_exponential_value() const {return v_;} +typedef Unit Dalton; +typedef Unit Kilodalton; - void show(std::ostream &out) const { - out << v_ << get_name(); - } - - //! Get a string representing the units - std::string get_name() const { - int exp=EXP; - int kcal=KCal; - int a=A; - - std::ostringstream name; - handle_unit(name, "kCal", kcal); - handle_unit(name, "A", a); - - if (kcal < 0 || a < 0) { - name << "/"; - handle_unit(name, "kCal", -kcal); - handle_unit(name, "A", -a); - } - return name.str(); - } - - template - KCalPerMolUnit - operator*(KCalPerMolUnit o) const { - return KCalPerMolUnit(v_*o.v_); - } - - template - KCalPerMolUnit - operator/(KCalPerMolUnit o) const { - return KCalPerMolUnit(v_/o.v_); - } - - template - This operator+(KCalPerMolUnit o) const { - return This(v_+o.v_); - } - template - This operator-(KCalPerMolUnit o) const { - return This(v_-o.v_); - } - - template - bool operator<(KCalPerMolUnit o) const { - return v_ < o.v_; - } - template - bool operator>(KCalPerMolUnit o) const { - return v_ > o.v_; - } -}; - - - -// want it for conversions -extern IMPDLLEXPORT const ExponentialNumber<23> NA; - -extern IMPDLLEXPORT const ExponentialNumber<3> JOULES_PER_KCAL; -extern IMPDLLEXPORT const ExponentialNumber<-4> KCALS_PER_JOULE; - - -/** \fn convert_to_kcal - \brief Convert from MKS units to kcal units - - \fn convert_to_mks - \brief Convert from kcal units to MKS units. - - These don't yet work for arbitrary units since it is kind of - annoying to get right and there are not equivalents for everything. - */ template -KCalPerMolUnit -convert_to_kcal(MKSUnit j) { - return KCalPerMolUnit(j.get_exponential_value() - *KCALS_PER_JOULE*NA); +Unit +convert_to_mks(Unit d) { + return Unit(d.get_value()/NA.get_value()); } template -KCalPerMolUnit -convert_to_kcal(MKSUnit j) { - internal::ExponentialNumber en(j.get_exponential_value() - *KCALS_PER_JOULE*NA - *internal::ExponentialNumber<-10>(1)); - return KCalPerMolUnit(en); +Unit +convert_to_Dalton(Unit d) { + return Unit(d.get_value()*NA.get_value()); } -template -MKSUnit -convert_to_mks(KCalPerMolUnit o) { - return MKSUnit(o.get_exponential_value() - *JOULES_PER_KCAL/NA); -} +} // namespace unit -template -MKSUnit -convert_to_mks(KCalPerMolUnit o) { - return MKSUnit(o.get_exponential_value() - *JOULES_PER_KCAL/NA - *internal::ExponentialNumber<10>(1)); -} - -//! The unit for MD energies -typedef KCalPerMolUnit<0, 1, 0> KCalPerMol; -//! The unit for MD-compatible derivatives -typedef KCalPerMolUnit<0, 1, -1> KCalPerAMol; - - -template -std::ostream &operator<<(std::ostream &out, KCalPerMolUnit u) -{ - u.show(out); - return out; -} - } // namespace internal +namespace unit= internal::unit; + } // namespace IMP #endif /* __IMP_UNITS_H */ Index: kernel/include/IMP/internal/ObjectContainer.h =================================================================== --- kernel/include/IMP/internal/ObjectContainer.h (revision 593) +++ kernel/include/IMP/internal/ObjectContainer.h (working copy) @@ -46,8 +46,32 @@ unsigned int get_index(II i) const {return i.get_index();} unsigned int get_index(unsigned int i) const {return i;} + void check_unique(O* o) { + IMP_CHECK_OBJECT(o); + for (typename Vector::const_iterator it= data_.begin(); + it != data_.end(); ++it) { + IMP_assert(*it != o, "IMP Containers can only have one copy of " + << " each object"); + } + } + + template + void check_all_unique(It b, It e) { + for (It cc= b; cc != e; ++cc) { + check_unique(*cc); + } + } + + I fill_free(O *o) { + unsigned int i= free_.back(); + free_.pop_back(); + data_[i]= o; + return I(i); + } + // hide it typedef std::vector Vector; + public: ObjectContainer(){} @@ -100,31 +124,22 @@ O* operator[](I i) const { IMP_check(get_index(i) < data_.size(), "Index " << i << " out of range", - IndexException("Out of range")); + IndexException); IMP_assert(data_.operator[](get_index(i)) != NULL, "Attempting to access invalid slot in container"); return data_[get_index(i)]; } I push_back(O* d) { - IMP_CHECK_OBJECT(d); own(d); IMP_IF_CHECK(EXPENSIVE) { - for (typename Vector::const_iterator it= data_.begin(); - it != data_.end(); ++it) { - IMP_assert(*it != d, "IMP Containers can only have one copy of " - << " each object"); - } + check_unique(d); } if (free_.empty()) { data_.push_back(d); - unsigned int idx= data_.size()-1; - return I(idx); + return I(data_.size()-1); } else { - unsigned int i= free_.back(); - free_.pop_back(); - data_[i]= d; - return I(i); + return fill_free(d); } } @@ -132,22 +147,13 @@ void insert(iterator c, It b, It e) { IMP_assert(c== end(), "Insert position is ignored in ObjectContainer"); IMP_IF_CHECK(EXPENSIVE) { - for (It cc= b; cc != e; ++cc) { - IMP_CHECK_OBJECT(*cc); - for (typename Vector::const_iterator it= data_.begin(); - it != data_.end(); ++it) { - IMP_assert(*it != *cc, "IMP Containers can only have one copy of " - << " each object"); - } - } + check_all_unique(b,e); } for (It cc= b; cc != e; ++cc) { own(*cc); } while (!free_.empty()) { - int i= free_.back(); - free_.pop_back(); - data_[i]= *b; + fill_free(*b); ++b; } data_.insert(data_.end(), b, e); Index: kernel/include/IMP/internal/AttributeTable.h =================================================================== --- kernel/include/IMP/internal/AttributeTable.h (revision 593) +++ kernel/include/IMP/internal/AttributeTable.h (working copy) @@ -146,7 +146,7 @@ IMP_check(Traits::get_is_valid(map_[k.get_index()]), "Attribute \"" << k.get_string() << "\" not found in table.", - IndexException("Invalid attribute")); + IndexException); } public: @@ -179,7 +179,7 @@ check_contains(k); IMP_check(Traits::get_is_valid(v), "Cannot set value of attribute to " << v, - ValueException("Invalid value for attribute")); + ValueException); map_[k.get_index()] = v; } @@ -260,7 +260,7 @@ IMP_check(Traits::get_is_valid(v), "Trying to insert invalid value for attribute " << v << " into attribute " << k, - ValueException("Invalid attribute value")); + ValueException); if (size_ <= k.get_index()) { boost::scoped_array old; swap(old, map_); Index: kernel/include/IMP/internal/Grid3D.h =================================================================== --- kernel/include/IMP/internal/Grid3D.h (revision 593) +++ kernel/include/IMP/internal/Grid3D.h (working copy) @@ -156,7 +156,7 @@ friend class GridIndexIterator; GridIndex(int x, int y, int z): VirtualGridIndex(x,y,z) { IMP_check(x>=0 && y>=0 && z>=0, "Bad indexes in grid index", - IndexException("Bad indexes in GridIndex")); + IndexException); } }; Index: kernel/include/IMP/internal/constants.h =================================================================== --- kernel/include/IMP/internal/constants.h (revision 593) +++ kernel/include/IMP/internal/constants.h (working copy) @@ -21,16 +21,20 @@ // make them doubles so we don't have to worry about digits //! Avagadro's number -extern IMPDLLEXPORT const ExponentialNumber<23> NA; +extern IMPDLLEXPORT const unit::ExponentialNumber<23> NA; + //! Boltzmann constant in J/K -extern IMPDLLEXPORT const MKSUnit<-23, 2, 1, -1, -2> KB; +extern IMPDLLEXPORT const +unit::Shift::type, -23>::type KB; //! the default temperature -extern IMPDLLEXPORT const Kelvin T; +extern IMPDLLEXPORT const unit::Kelvin T; //! Pi static const double PI = 3.1415926535897931; +extern IMPDLLEXPORT const unit::ExponentialNumber<3> JOULES_PER_KILOCALORIE; + } // namespace internal } // namespace IMP Index: kernel/include/IMP/internal/Unit.h =================================================================== --- kernel/include/IMP/internal/Unit.h (revision 0) +++ kernel/include/IMP/internal/Unit.h (revision 0) @@ -0,0 +1,408 @@ +/** + * \file Unit.h \brief Classes to help with converting between units. + * + * Copyright 2007-8 Sali Lab. All rights reserved. + * + */ + +#ifndef __IMP_UNIT_H +#define __IMP_UNIT_H + + +#include "ExponentialNumber.h" +#include "../utility.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include + + +namespace IMP +{ + +namespace internal +{ + + +namespace unit +{ + + +namespace internal +{ + +template +std::string get_unit_name(int o) { + bool should_not_be_compiled; + return ""; +} + +// for stupid line length limits +using namespace boost::mpl; +using namespace boost::mpl::placeholders; + + +/* + These classes allow the Units part of a unit to be manipulated + */ +template +struct Divide { + typedef typename transform >::type type; +}; + +template +struct Multiply { + typedef typename transform >::type type; +}; + +template +struct Sqrt { + typedef typename transform > >::type type; +}; + + +template +struct Exponentiate { + typedef typename transform >::type type; +}; + +template +struct DoNormalize { + typedef int type; +}; + +/* + The boost MPI multiply sorts of routines return a slightly different + type that we want so we have to change it into the exact type to get + certain things work work right. + */ +template +struct DoNormalize { + typedef vector_c type; +}; + +template +struct DoNormalize { + typedef vector_c >::type::value> type; +}; + +template +struct DoNormalize { + typedef vector_c >::type::value, + at >::type::value> type; +}; + +template +struct DoNormalize { + typedef vector_c >::type::value, + at >::type::value, + at >::type::value> type; +}; + +template +struct DoNormalize { + typedef vector_c >::type::value, + at >::type::value, + at >::type::value, + at >::type::value> type; +}; + +template +struct DoNormalize { + typedef vector_c >::type::value, + at >::type::value, + at >::type::value, + at >::type::value, + at >::type::value> type; +}; + +template +struct Normalize { + typedef typename DoNormalize::type::value > ::type type; +}; + +// Recursive type to print out the unit names and powers +template +struct PrintUnits { + PrintUnits p_; + + void operator()(std::ostream &out) const { + std::string str= get_unit_name(O); + int e= boost::mpl::at >::type::value; + if (e != 0) { + out << " " << str; + if (e != 1) { + out << "^" << e; + } + } + p_(out); + } +}; + +//! Terminate the recursion +template +struct PrintUnits { + void operator()(std::ostream &out) const{}; +}; + + +//! Specializaton for singleton units +template +struct PrintUnits > { + void operator()(std::ostream &out) const{ + std::string str= get_unit_name(O); + out << " " << str; + } +}; + +// recursive type to check if all the unit powers are zero +template +struct IsNoUnits { + typedef IsNoUnits Next; + static const bool value= Next::value + && !(boost::mpl::at >::type::value); +}; + + +template +struct IsNoUnits { + static const bool value =true; +}; + + + +} // unit::internal + + + + +//! A base class for units +/** + \internal + */ +template +class Unit { + template + friend class Unit; + typedef Unit This; + typedef ExponentialNumber V; + V v_; + + + bool is_default() const { + return v_== V(); + } +public: + typedef TagT Tag; + typedef UnitsT Units; + static const int EXP= EXPT; + + explicit Unit(V v): v_(v){} + explicit Unit(double v): v_(v){} + //Unit(int v): v_(v){} + Unit(){} + template + Unit(Unit o): v_(o.v_) { + BOOST_STATIC_ASSERT(( + boost::mpl::equal::type::value + )); + } + + template + Unit(Unit o): v_(o.v_) {} + + double to_scalar() const { + BOOST_STATIC_ASSERT((internal::IsNoUnits<0, + boost::mpl::size::type::value, Units>::value)); + return get_normalized_value(); + } + + IMP_COMPARISONS_1(v_); + + //! Get the current value + /** \note The value returned is the value before it is multiplied + by the appropriate power of 10. That means, 1 Angstrom returns 1, + as does 1 Meter. + */ + double get_value() const { + return v_.get_value(); + } + //! Get the value with the proper exponent + /** 1 Angstrom returns 1e-10. + */ + double get_normalized_value() const { + return v_.get_normalized_value(); + } + + V get_exponential_value() const { + return v_; + } + + + template + Unit::type > + operator*(Unit o) const { + return Unit::type >(v_*o.v_); + } + + template + Unit::type > + operator/(Unit o) const { + return Unit::type > (v_/o.v_); + } + + This operator+(This o) const { + return This(v_+o.v_); + } + This operator-(This o) const { + return This(v_-o.v_); + } + + This operator-() const { + return This(-v_); + } + + void show(std::ostream &out) const { + out << v_; + internal::PrintUnits::type::value, Units> p; + p(out); + } +}; + +template +std::ostream &operator<<(std::ostream &out, + Unit o) { + o.show(out); + return out; +} + + +typedef boost::mpl::vector_c SingletonUnit; + + + +/** need to be careful of integer division + \internal + */ +template +Unit::type > +sqrt(Unit o) { + return Unit::type > + (ExponentialNumber(std::sqrt(o.get_normalized_value()))); +} + +/** \internal + */ +template +Unit::type > +square(Unit o) { + return Unit::type > + (::IMP::square(o.get_value())); +} + +template +Unit +operator*(Unit o, double d) { + return Unit(d*o.get_value()); +} + +template +Unit +operator*(double d, Unit o) { + return Unit(d*o.get_value()); +} + +template +Unit +operator/(Unit o, double d) { + return Unit(o.get_value()/d); +} + +template +Unit +operator/(double d, Unit o) { + return Unit(d/o.get_value()); +} + +// Multiply and divide Unit instantiations + +template +struct Divide { + BOOST_STATIC_ASSERT((boost::mpl::equal::type::value)); + typedef typename internal::Divide::type raw_units; + typedef typename internal::Normalize::type units; + typedef Unit type; + +}; + +template +struct Multiply { + BOOST_STATIC_ASSERT((boost::mpl::equal::type::value)); + typedef typename internal::Multiply::type raw_units; + typedef typename internal::Normalize::type units; + typedef Unit type; +}; + +template +struct Shift { + typedef Unit type; + +}; + + + +template +struct Exchange { + BOOST_STATIC_ASSERT((boost::mpl::equal::type::value)); + BOOST_STATIC_ASSERT((boost::mpl::equal::type::value)); + typedef typename internal::Divide::type Div; + typedef typename internal::Multiply::type Mul; + typedef typename internal::Normalize::type units; + typedef Unit type; + +}; + +} // unit + +} // internal + + +} // IMP + + +#endif Index: kernel/include/IMP/internal/ArrayOnAttributesHelper.h =================================================================== --- kernel/include/IMP/internal/ArrayOnAttributesHelper.h (revision 593) +++ kernel/include/IMP/internal/ArrayOnAttributesHelper.h (working copy) @@ -39,10 +39,10 @@ Value get_value(const Particle *p, unsigned int i) const { IMP_check(keys_.size() > i, "Out of range attribute in array", - IndexException("out of range")); + IndexException); IMP_check(static_cast(p->get_value(num_key_)) > i, "Out of range attribute in array", - IndexException("out of range")); + IndexException); return p->get_value(keys_[i]); } @@ -50,9 +50,9 @@ unsigned int i, Value v) const { IMP_check(keys_.size() > i, "Out of range attribute in array", - IndexException("out of range")); + IndexException); IMP_check(p->get_value(num_key_) > i, "Out of range attribute in array", - IndexException("out of range")); + IndexException); p->set_value(keys_[i], v); } @@ -70,7 +70,7 @@ Value v) { unsigned int osz= p->get_value(num_key_); IMP_check(loc <= osz, "Attribute array must be contiguous", - IndexException("out of range")); + IndexException); for (unsigned int i=loc; i < osz; ++i) { Key k= get_key(i); Value t= p->get_value(k); @@ -87,7 +87,7 @@ unsigned int loc) const { unsigned int osz= p->get_value(num_key_); IMP_check(loc <= osz, "Can only erase values in array", - IndexException("out of range")); + IndexException); for (unsigned int i=loc+1; i < osz; ++i) { Key k= keys_[i]; Key kl= keys_[i-1]; Index: kernel/include/IMP/internal/RefCountedObject.h =================================================================== --- kernel/include/IMP/internal/RefCountedObject.h (revision 593) +++ kernel/include/IMP/internal/RefCountedObject.h (working copy) @@ -21,8 +21,7 @@ { //! Common base class for ref counted objects. -/** Currently the only ref counted objects are particles. - This class acts as a tag rather than providing any functionality. +/** This class acts as a tag rather than providing any functionality. \internal */ @@ -104,8 +103,6 @@ void unref(O o) { BOOST_STATIC_ASSERT(!boost::is_pointer::value); - /*IMP_LOG(VERBOSE, "NonUnRef called with nonpointer for " - << o << std::endl);*/ } @@ -114,17 +111,12 @@ void ref(O o) { BOOST_STATIC_ASSERT(!boost::is_pointer::value); - /*IMP_LOG(VERBOSE, "NonRef count called with nonpointer for " - << o << std::endl);*/ } //! Can be called on any object and will only unref it if appropriate template void unref(O* o) { - /*IMP_LOG(VERBOSE, "Unref count called with " - << (boost::is_base_of::value) - << " for " << o << std::endl);*/ UnRef<(boost::is_base_of::value)>::eval(o); } @@ -133,9 +125,6 @@ template void ref(O* o) { - /*IMP_LOG(VERBOSE, "Ref called with " - << (boost::is_base_of::value) - << " for " << o << std::endl);*/ Ref<(boost::is_base_of::value)>::eval(o); } @@ -167,7 +156,7 @@ } else { IMP_check(!o->get_has_ref(), "Trying to own already owned but " << "non-reference-counted object: " << *o, - ValueException("Already owned object")); + ValueException); } o->ref(); } Index: kernel/include/IMP/exception.h =================================================================== --- kernel/include/IMP/exception.h (revision 593) +++ kernel/include/IMP/exception.h (working copy) @@ -16,6 +16,7 @@ #include #include #include +#include namespace IMP { @@ -69,7 +70,7 @@ */ struct IMPDLLEXPORT ErrorException: public Exception { - ErrorException(): Exception("Fatal error"){} + ErrorException(const char *msg="Fatal error"): Exception(msg){} }; //! An exception for an invalid model state @@ -87,8 +88,9 @@ class IMPDLLEXPORT InactiveParticleException : public Exception { public: - InactiveParticleException(): - Exception("Attempting to use inactive particle"){} + InactiveParticleException(const char *msg + ="Attempting to use inactive particle"): + Exception(msg){} }; //! An exception for a request for an invalid member of a container @@ -137,13 +139,13 @@ /** Break on exception.cpp:31 to catch assertion failures. \ingroup assert */ -IMPDLLEXPORT void assert_fail(); +IMPDLLEXPORT void assert_fail(const char *msg); //! Here so you can catch check failures more easily in the debugger /** Break on exception.cpp:35 to catch check failures. \ingroup assert */ -IMPDLLEXPORT void check_fail(); +IMPDLLEXPORT void check_fail(const char *msg); } // namespace internal @@ -161,12 +163,15 @@ \param[in] message Write this message if the assertion fails. \ingroup assert */ -#define IMP_assert(expr, message) \ - do { \ - if (IMP::get_check_level() >= IMP::EXPENSIVE && !(expr)) { \ - IMP_ERROR(message); \ - IMP::internal::assert_fail(); \ - } \ +#define IMP_assert(expr, message) \ + do { \ + if (IMP::get_check_level() >= IMP::EXPENSIVE && !(expr)) { \ + std::ostringstream oss; \ + oss << message << std::endl \ + << " File \"" << __FILE__ << "\", line " << __LINE__ \ + << std::endl; \ + IMP::internal::assert_fail(oss.str().c_str()); \ + } \ } while(false) #else #define IMP_assert(expr, message) @@ -175,16 +180,18 @@ //! A runtime check for IMP. /** \param[in] expr The assertion expression. \param[in] message Write this message if the assertion fails. - \param[in] exception Throw the object constructed by this expression. + \param[in] ExceptionType Throw an exception of this type. The exception + must be constructable from a char *. \ingroup assert */ -#define IMP_check(expr, message, exception) \ - do { \ - if (IMP::get_check_level() >= IMP::CHEAP && !(expr)) { \ - IMP_ERROR(message); \ - IMP::internal::check_fail(); \ - throw exception; \ - } \ +#define IMP_check(expr, message, ExceptionType) \ + do { \ + if (IMP::get_check_level() >= IMP::CHEAP && !(expr)) { \ + std::ostringstream oss; \ + oss << message << std::endl; \ + IMP::internal::check_fail(oss.str().c_str()); \ + throw ExceptionType(oss.str().c_str()); \ + } \ } while (false) //! A runtime failure for IMP. Index: kernel/include/IMP/unary_functions/WormLikeChain.h =================================================================== --- kernel/include/IMP/unary_functions/WormLikeChain.h (revision 593) +++ kernel/include/IMP/unary_functions/WormLikeChain.h (working copy) @@ -31,7 +31,11 @@ /** \param[in] l_max maximum length of the chain in angstroms \param[in] lp persistence length in angstroms */ - WormLikeChain(Float l_max, Float lp) : lmax_(l_max), lp_(lp) {} + WormLikeChain(Float l_max, Float lp) : lmax_(l_max), lp_(lp) { + IMP_check(l_max > lp, "The persistance length should be less " + << "than the total length for this model", + ValueException); + } virtual ~WormLikeChain() {} @@ -40,19 +44,20 @@ \return Energy in kcal/mol */ virtual Float evaluate(Float lf) { - static const internal::PicoJoule zero=eval(internal::Angstrom(0)); - internal::Angstrom l(lf); - if (l < internal::Angstrom(0)) l=internal::Angstrom(0); - internal::PicoJoule ret; + static const unit::Picojoule zero=eval(unit::Angstrom(0)); + unit::Angstrom l(lf); + if (l < unit::Angstrom(0)) l=unit::Angstrom(0); + unit::Picojoule ret; if (l < cutoff()) { ret= (eval(l) - zero); } else { - internal::PicoJoule springterm=(l-cutoff())*cderiv(cutoff()); + unit::Picojoule springterm=(l-cutoff())*cderiv(cutoff()); ret= (eval(cutoff())+ springterm -zero); } - /*std::cout << "Return is " << ret <<" " << l << " " << lp_ << " " - << lmax_ << std::endl;*/ - return internal::KCalPerMol(convert_to_kcal(ret)).get_value(); + std::cout << "Return is " << ret <<" " << l << " " << lp_ << " " + << lmax_ << std::endl; + unit::YoctoKilocalorie zc= convert_J_to_Cal(ret); + return (zc*unit::ATOMS_PER_MOL).get_value(); } //! Calculate the WormLikeChain energy given the length @@ -61,9 +66,9 @@ \return Score */ virtual Float evaluate_deriv(Float fl, Float& deriv) { - internal::Angstrom l(fl); - if (l < internal::Angstrom(0)) l=internal::Angstrom(0); - internal::PicoNewton doubled; + unit::Angstrom l(fl); + if (l < unit::Angstrom(0)) l=unit::Angstrom(0); + unit::Piconewton doubled; if (l < cutoff()) { doubled= cderiv(l); } else { @@ -74,8 +79,9 @@ } //std::cout << "Force is " << doubled << " for length " << l << std::endl; // convert from picoNewton - deriv = internal::KCalPerAMol(internal::convert_to_kcal(doubled)) - .get_value(); + unit::YoctoKilocaloriePerAngstrom du= unit::convert_J_to_Cal(doubled); + + deriv = (du*unit::ATOMS_PER_MOL).get_value(); //std::cout << "Which converts to " << d << std::endl; return evaluate(fl); } @@ -86,28 +92,28 @@ protected: //! \note named to avoid clash with 'deriv' argument - internal::PicoNewton cderiv(internal::Angstrom l) const { - internal::PicoNewton pn= internal::KB*internal::T - /lp_*(internal::Scalar(.25)/ square(internal::Scalar(1)-l/lmax_) - -internal::Scalar(.25)+l/lmax_); + unit::Piconewton cderiv(unit::Angstrom l) const { + unit::Piconewton pn= internal::KB*internal::T + /lp_*(.25/ square(1.0-(l/lmax_).get_normalized_value()) + -.25+(l/lmax_).to_scalar()); return pn; } - internal::PicoJoule eval(internal::Angstrom m) const { - internal::PicoJoule J - = internal::KB*internal::T/lp_*(internal::Scalar(.25)*square(lmax_) + unit::Picojoule eval(unit::Angstrom m) const { + unit::Picojoule J + = internal::KB*internal::T/lp_*(.25*square(lmax_) /(lmax_-m) - -m*internal::Scalar(.25) - +internal::Scalar(.5)*square(m) + -m*.25 + +.5*square(m) /lmax_); return J; } - internal::Angstrom cutoff() const { - return internal::Scalar(.99)*lmax_; + unit::Angstrom cutoff() const { + return .99*lmax_; } - internal::Angstrom lmax_, lp_; + unit::Angstrom lmax_, lp_; }; } // namespace IMP Index: kernel/include/IMP/Particle.h =================================================================== --- kernel/include/IMP/Particle.h (revision 593) +++ kernel/include/IMP/Particle.h (working copy) @@ -407,41 +407,41 @@ inline Float Particle::get_value(FloatKey name) const { IMP_check(get_is_active(), "Do not touch inactive particles", - InactiveParticleException()); + InactiveParticleException); return floats_.get_value(name); } inline Float Particle::get_derivative(FloatKey name) const { IMP_check(get_is_active(), "Do not touch inactive particles", - InactiveParticleException()); + InactiveParticleException); return derivatives_.get_value(name); } inline void Particle::set_value(FloatKey name, Float value) { IMP_check(get_is_active(), "Do not touch inactive particles", - InactiveParticleException()); + InactiveParticleException); floats_.set_value(name, value); } inline bool Particle::get_is_optimized(FloatKey name) const { IMP_check(get_is_active(), "Do not touch inactive particles", - InactiveParticleException()); + InactiveParticleException); IMP_check(floats_.contains(name), "get_is_optimized called " << "with invalid attribute" << name, - IndexException("Invalid float attribute")); + IndexException); return optimizeds_.contains(name); } inline void Particle::set_is_optimized(FloatKey name, bool tf) { IMP_check(get_is_active(), "Do not touch inactive particles", - InactiveParticleException()); + InactiveParticleException); IMP_check(floats_.contains(name), "set_is_optimized called " << "with invalid attribute" << name, - IndexException("Invalid float attribute")); + IndexException); if (tf) { optimizeds_.insert_always(name, true); } else { @@ -453,7 +453,7 @@ const DerivativeAccumulator &da) { IMP_check(get_is_active(), "Do not touch inactive particles", - InactiveParticleException()); + InactiveParticleException); IMP_assert(value==value, "Can't add NaN to derivative in particle " << *this); derivatives_.set_value(name, derivatives_.get_value(name)+ da(value)); } @@ -461,7 +461,7 @@ inline bool Particle::has_attribute(IntKey name) const { IMP_check(get_is_active(), "Do not touch inactive particles", - InactiveParticleException()); + InactiveParticleException); return ints_.contains(name); } @@ -470,7 +470,7 @@ inline Int Particle::get_value(IntKey name) const { IMP_check(get_is_active(), "Do not touch inactive particles", - InactiveParticleException()); + InactiveParticleException); return ints_.get_value(name); } @@ -478,14 +478,14 @@ inline void Particle::set_value(IntKey name, Int value) { IMP_check(get_is_active(), "Do not touch inactive particles", - InactiveParticleException()); + InactiveParticleException); ints_.set_value(name, value); } inline bool Particle::has_attribute(StringKey name) const { IMP_check(get_is_active(), "Do not touch inactive particles", - InactiveParticleException()); + InactiveParticleException); return strings_.contains(name); } @@ -494,14 +494,14 @@ inline String Particle::get_value(StringKey name) const { IMP_check(get_is_active(), "Do not touch inactive particles", - InactiveParticleException()); + InactiveParticleException); return strings_.get_value(name); } inline void Particle::set_value(StringKey name, String value) { IMP_check(get_is_active(), "Do not touch inactive particles", - InactiveParticleException()); + InactiveParticleException); strings_.set_value(name, value); } @@ -509,7 +509,7 @@ inline bool Particle::has_attribute(ParticleKey name) const { IMP_check(get_is_active(), "Do not touch inactive particles", - InactiveParticleException()); + InactiveParticleException); return particles_.contains(name); } @@ -518,7 +518,7 @@ inline Particle* Particle::get_value(ParticleKey name) const { IMP_check(get_is_active(), "Do not touch inactive particles", - InactiveParticleException()); + InactiveParticleException); return particles_.get_value(name).get(); } @@ -526,7 +526,7 @@ inline void Particle::set_value(ParticleKey name, Particle* value) { IMP_check(get_is_active(), "Do not touch inactive particles", - InactiveParticleException()); + InactiveParticleException); particles_.set_value(name, internal::ObjectPointer(value)); } Index: kernel/include/IMP/Index.h =================================================================== --- kernel/include/IMP/Index.h (revision 593) +++ kernel/include/IMP/Index.h (working copy) @@ -26,7 +26,7 @@ //! Construct an index from a nonnegative int Index(unsigned int i): i_(i) { IMP_check(i >= 0, "Index initializer must be positive. " << i << " is not.", - ErrorException()); + IndexException); } //! Construct a default index /** This can be used as a sentinal value */ @@ -36,7 +36,7 @@ /** The integer is unique within the container */ unsigned int get_index() const { IMP_check(i_ >= 0, "get_index() called on defaultly constructed Index", - ErrorException()); + ValueException); return i_; } void show(std::ostream &out) const { Index: kernel/src/exception.cpp =================================================================== --- kernel/src/exception.cpp (revision 593) +++ kernel/src/exception.cpp (working copy) @@ -30,13 +30,22 @@ namespace internal { -void assert_fail() +// The error message is already in the exception +bool print_exceptions=false; + +void assert_fail(const char *msg) { - throw ErrorException(); + if (print_exceptions) { + IMP_ERROR(msg); + } + throw ErrorException(msg); } -void check_fail() +void check_fail(const char *msg) { + if (print_exceptions) { + IMP_ERROR(msg); + } } } // namespace internal Index: kernel/src/Particle.cpp =================================================================== --- kernel/src/Particle.cpp (revision 593) +++ kernel/src/Particle.cpp (working copy) @@ -27,12 +27,12 @@ { IMP_check(model_==NULL || md==NULL, "Set_model called for particle already in model", - ValueException("Cannot transfer particles directly")); + ValueException); model_ = md; pi_ = pi; IMP_check(model_==NULL || model_->get_particle(pi_)== this, "Set_model called with inconsistent data", - ValueException("Cannot transfer particles directly")); + ValueException); } void Particle::set_is_active(const bool is_active) Index: kernel/src/restraints/ConnectivityRestraint.cpp =================================================================== --- kernel/src/restraints/ConnectivityRestraint.cpp (revision 593) +++ kernel/src/restraints/ConnectivityRestraint.cpp (working copy) @@ -53,7 +53,7 @@ void ConnectivityRestraint::add_set(const Particles &ps) { IMP_check(!ps.empty(), "Cannot add empty set to ConnectivityRestraint", - InvalidStateException("Empty set")); + InvalidStateException); unsigned int sz= ps.size(); Restraint::add_particles(ps); set_offsets_.push_back(set_offsets_.back()+sz); @@ -92,8 +92,7 @@ << tag_weight << " as a marker, so the distance function " << " should not return it. This can be fixed if you " << " complain.", - ValueException("Tag weight matches distance in"\ - " ConnectivityRestraint")); + ValueException); IMP_LOG(VERBOSE, "ConnectivityRestraint edge between " << get_particle(ii)->get_index() << " and " << get_particle(ij)->get_index() << " with weight " Index: kernel/src/score_states/MaxChangeScoreState.cpp =================================================================== --- kernel/src/score_states/MaxChangeScoreState.cpp (revision 593) +++ kernel/src/score_states/MaxChangeScoreState.cpp (working copy) @@ -33,7 +33,7 @@ IMP_check(obj->has_attribute(keys_[i]), "Particle missing needed attribute " << keys_[i] << obj, - ValueException("Particle missing attribute")); + ValueException); }; for (unsigned int i=0; i< origkeys_.size(); ++i) { if (!obj->has_attribute(origkeys_[i])) { Index: kernel/src/decorators/MolecularHierarchyDecorator.cpp =================================================================== --- kernel/src/decorators/MolecularHierarchyDecorator.cpp (revision 593) +++ kernel/src/decorators/MolecularHierarchyDecorator.cpp (working copy) @@ -113,7 +113,7 @@ || mhd.get_type() == MolecularHierarchyDecorator::CHAIN || mhd.get_type() == MolecularHierarchyDecorator::NUCLEOTIDE, "Invalid type of MolecularHierarchyDecorator passed to get_residue", - ValueException("Bad val")); + ValueException); MatchResidueIndex mi(index); HierarchyDecorator hd= hierarchy_find(mhd, mi); if (hd== HierarchyDecorator()) { @@ -129,14 +129,14 @@ create_fragment(const MolecularHierarchyDecorators &ps) { IMP_check(!ps.empty(), "Need some particles", - ValueException("")); + ValueException); MolecularHierarchyDecorator parent= ps[0].get_parent(); unsigned int index= ps[0].get_parent_index(); IMP_IF_CHECK(CHEAP) { for (unsigned int i=0; i< ps.size(); ++i) { IMP_check(ps[i].get_parent() == parent, "Parents don't match", - ValueException("")); + ValueException); } } Index: kernel/src/decorators/bond_decorators.cpp =================================================================== --- kernel/src/decorators/bond_decorators.cpp (revision 593) +++ kernel/src/decorators/bond_decorators.cpp (working copy) @@ -80,7 +80,7 @@ { IMP_check(a.get_particle() != b.get_particle(), "The endpoints of a bond must be disjoint", - ValueException("The endpoints must be disjoint")); + ValueException); Particle *p= internal::graph_connect(a.get_particle(), b.get_particle(), internal::bond_graph_data_); Index: kernel/src/decorators/XYZDecorator.cpp =================================================================== --- kernel/src/decorators/XYZDecorator.cpp (revision 593) +++ kernel/src/decorators/XYZDecorator.cpp (working copy) @@ -30,7 +30,7 @@ float radius) { IMP_check(radius > 0, "Radius in randomize must be postive", - ValueException("Radius must be positive")); + ValueException); Vector3D min(center[0]-radius, center[1]-radius, center[2]-radius); Vector3D max(center[0]+radius, center[1]+radius, center[2]+radius); float norm; @@ -49,7 +49,7 @@ { for (unsigned int i=0; i< 3; ++i) { IMP_check(min[i] < max[i], "Box for randomize must be non-empty", - ValueException("Box must be non-empty")); + ValueException); ::boost::uniform_real<> rand(min[i], max[i]); set_coordinate(i, rand(random_number_generator)); } Index: kernel/src/optimizers/ConjugateGradients.cpp =================================================================== --- kernel/src/optimizers/ConjugateGradients.cpp (revision 593) +++ kernel/src/optimizers/ConjugateGradients.cpp (working copy) @@ -34,7 +34,7 @@ /* set model state */ for (i = 0; i < opt_var_cnt; i++) { IMP_check(x[i] == x[i], "Got NaN in CG", - ValueException("Got NaN in CG")); + ValueException); if (std::abs(x[i] - get_value(float_indices[i])) > max_change_) { if (x[i] < get_value(float_indices[i])) { x[i] = get_value(float_indices[i]) - max_change_; @@ -54,7 +54,7 @@ IMP_check(dscore[i] == dscore[i] && dscore[i] != std::numeric_limits::infinity() && dscore[i] != - std::numeric_limits::infinity(), - "Bad input to CG", ValueException("Bad input to CG")); + "Bad input to CG", ValueException); } return score; } @@ -230,7 +230,7 @@ { IMP_check(get_model(), "Must set the model on the optimizer before optimizing", - ValueException("Must set the model before optimizing")); + ValueException); std::vector x, dx; int i; //ModelData* model_data = get_model()->get_model_data(); @@ -246,7 +246,7 @@ x[i] = get_value(float_indices[i]); IMP_check(x[i] == x[i] && x[i] != std::numeric_limits::infinity() && x[i] != - std::numeric_limits::infinity(), - "Bad input to CG", ValueException("Bad input to CG")); + "Bad input to CG", ValueException); } // Initialize optimization variables Index: kernel/src/optimizers/MonteCarlo.cpp =================================================================== --- kernel/src/optimizers/MonteCarlo.cpp (revision 593) +++ kernel/src/optimizers/MonteCarlo.cpp (working copy) @@ -44,7 +44,7 @@ IMP_check(cg_->get_model() == get_model(), "The model used by the local optimizer does not match "\ " that used by the montecarlo optimizer", - InvalidStateException("Bad model pointer")); + InvalidStateException); } update_states(); prior_energy_ =get_model()->evaluate(0); Index: kernel/src/optimizers/BrownianDynamics.cpp =================================================================== --- kernel/src/optimizers/BrownianDynamics.cpp (revision 593) +++ kernel/src/optimizers/BrownianDynamics.cpp (working copy) @@ -17,31 +17,34 @@ namespace IMP { -typedef internal::MKSUnit<-3, -1, 1, 0, -1> MilliPascalSeconds; +typedef +unit::Shift::type, + -3>::type MillipascalSecond; -static MilliPascalSeconds eta(internal::Kelvin T) +static MillipascalSecond eta(unit::Kelvin T) { - const std::pair points[] - ={ std::make_pair(internal::Kelvin(273+10.0), - MilliPascalSeconds(1.308)), - std::make_pair(internal::Kelvin(273+20.0), - MilliPascalSeconds(1.003)), - std::make_pair(internal::Kelvin(273+30.0), - MilliPascalSeconds(0.7978)), - std::make_pair(internal::Kelvin(273+40.0), - MilliPascalSeconds(0.6531)), - std::make_pair(internal::Kelvin(273+50.0), - MilliPascalSeconds(0.5471)), - std::make_pair(internal::Kelvin(273+60.0), - MilliPascalSeconds(0.4668)), - std::make_pair(internal::Kelvin(273+70.0), - MilliPascalSeconds(0.4044)), - std::make_pair(internal::Kelvin(273+80.0), - MilliPascalSeconds(0.3550)), - std::make_pair(internal::Kelvin(273+90.0), - MilliPascalSeconds(0.3150)), - std::make_pair(internal::Kelvin(273+100.0), - MilliPascalSeconds(0.2822))}; + const std::pair points[] + ={ std::make_pair(unit::Kelvin(273+10.0), + MillipascalSecond(1.308)), + std::make_pair(unit::Kelvin(273+20.0), + MillipascalSecond(1.003)), + std::make_pair(unit::Kelvin(273+30.0), + MillipascalSecond(0.7978)), + std::make_pair(unit::Kelvin(273+40.0), + MillipascalSecond(0.6531)), + std::make_pair(unit::Kelvin(273+50.0), + MillipascalSecond(0.5471)), + std::make_pair(unit::Kelvin(273+60.0), + MillipascalSecond(0.4668)), + std::make_pair(unit::Kelvin(273+70.0), + MillipascalSecond(0.4044)), + std::make_pair(unit::Kelvin(273+80.0), + MillipascalSecond(0.3550)), + std::make_pair(unit::Kelvin(273+90.0), + MillipascalSecond(0.3150)), + std::make_pair(unit::Kelvin(273+100.0), + MillipascalSecond(0.2822))}; const unsigned int npoints= sizeof(points)/sizeof(std::pair); if (T < points[0].first) { @@ -49,10 +52,11 @@ } else { for (unsigned int i=1; i< npoints; ++i) { if (points[i].first > T) { - internal::Scalar f= (T - points[i-1].first)/(points[i].first - - points[i-1].first); - MilliPascalSeconds ret= - (internal::Scalar(1.0)-f) *points[i-1].second + f*points[i].second; + float f= ((T - points[i-1].first) + /(points[i].first - points[i-1].first)) + .get_normalized_value(); + MillipascalSecond ret= + (1.0-f) *points[i-1].second + f*points[i].second; return ret; } } @@ -66,7 +70,7 @@ { IMP_check(dkey_ != FloatKey(), "BrownianDynamics needs a valid key for the " << "diffusion coefficient", - ValueException("Bad diffusion coef key")); + ValueException); } @@ -74,11 +78,10 @@ { } -IMP_LIST_IMPL(BrownianDynamics, Particle, particle, Particle*, - {if (0) std::cout << obj << index;},); +IMP_LIST_IMPL(BrownianDynamics, Particle, particle, Particle*,,); -void BrownianDynamics::set_time_step(internal::FemtoSecond t) +void BrownianDynamics::set_time_step(unit::Femtosecond t) { time_steps_.clear(); max_dt_ = t; @@ -118,7 +121,7 @@ std::vector new_pos(number_of_particles()); bool noshrink=true; while (!propose_step(new_pos)) { - cur_dt_= cur_dt_/internal::Scalar(2.0); + cur_dt_= cur_dt_/2.0; noshrink=false; } @@ -135,13 +138,6 @@ //! Perform a single dynamics step. bool BrownianDynamics::propose_step(std::vector& new_pos) { - // Assuming score is in kcal/mol, its derivatives in kcal/mol/angstrom, - // and mass is in g/mol, conversion factor necessary to get accelerations - // in angstrom/fs/fs from raw derivatives - // we want - //const double eta= 10;// kg/(m s) from ~1 cg/(cm s) - //const double pi=M_PI; - //const double kTnu= 1.3806504*T_.get_value(); FloatKeys xyzk=XYZDecorator::get_xyz_keys(); IMP_LOG(VERBOSE, "dt is " << cur_dt_ << std::endl); @@ -150,18 +146,41 @@ for (unsigned int i=0; i< number_of_particles(); ++i) { Particle *p= get_particle(i); XYZDecorator d(p); + IMP_IF_CHECK(CHEAP) { + for (unsigned int j=0; j< 3; ++j) { + // GDB 6.6 prints infinity as 0 on 64 bit machines. Grumble. + /*int szf= sizeof(Float); + int szi= sizeof(int); + Float one=1;*/ + Float mx= std::numeric_limits::max(); + Float c= d.get_coordinate(j); + bool ba= (c != c); + bool bb = (c >= mx); + bool bc= -d.get_coordinate(j) >= std::numeric_limits::max(); + if (ba || bb || bc ) { + IMP_WARN("Bad value for coordinate in Brownian Dynamics on " + << "particle " << p->get_index() << std::endl); + throw ValueException("Bad coordinate value"); + } + } + } + //double xi= 6*pi*eta*radius; // kg/s - internal::Cm2PerSecond D(p->get_value(dkey_)); - internal::Angstrom sigma(compute_sigma_from_D(D)); -#ifndef NDEBUG - internal::Angstrom osigma(sqrt(internal::Scalar(2.0)*D*cur_dt_)); -#endif - IMP_assert(sigma - - osigma - <= internal::Scalar(.01)* sigma, - "Sigma computations don't match " << sigma - << " " - << sqrt(internal::Scalar(2.0)*D*cur_dt_)); + unit::SquareCentimeterPerSecond D(p->get_value(dkey_)); + IMP_check(D.get_value() > 0 + && D.get_value() < std::numeric_limits::max(), + "Bad diffusion coefficient on particle " << p->get_index(), + ValueException); + unit::Angstrom sigma(compute_sigma_from_D(D)); + IMP_IF_CHECK(EXPENSIVE) { + unit::Angstrom osigma(sqrt(2.0*D*cur_dt_)); + IMP_check(sigma - osigma + <= .01* sigma, + "Sigma computations don't match " << sigma + << " " + << sqrt(2.0*D*cur_dt_), + ErrorException); + } IMP_LOG(VERBOSE, p->get_index() << ": sigma is " << sigma << std::endl); boost::normal_distribution mrng(0, sigma.get_value()); @@ -171,25 +190,22 @@ //std::cout << p->get_index() << std::endl; - internal::Angstrom delta[3]; + unit::Angstrom delta[3]; for (unsigned j = 0; j < 3; ++j) { - // force originally in kcal/mol/A - // derive* 2.390e-4 gives it in J/(mol A) - // that * e10 gives it in J/(mol m) - // that / NA gives it in J/m - internal::KCalPerAMol force(-d.get_coordinate_derivative(j)); - //*4.1868e+13/NA; old conversion - internal::FemtoNewton nforce= convert_to_mks(force); - internal::Angstrom R(sampler()); - internal::Angstrom force_term=nforce*cur_dt_*D/kt(); + unit::KilocaloriePerAngstromPerMol + force( -d.get_coordinate_derivative(j)); + unit::Femtonewton nforce + = unit::convert_Cal_to_J(force/unit::ATOMS_PER_MOL); + unit::Angstrom R(sampler()); + unit::Angstrom force_term(nforce*cur_dt_*D/kt()); //std::cout << "Force term is " << force_term << " and R is " //<< R << std::endl; - internal::Angstrom dr= force_term + R; //std::sqrt(2*kb*T_*ldt/xi) * + unit::Angstrom dr= force_term + R; //std::sqrt(2*kb*T_*ldt/xi) * // get back from meters delta[j]=dr; } - //internal::Angstrom max_motion= internal::Scalar(4)*sigma; + //unit::Angstrom max_motion= unit::Scalar(4)*sigma; /*std::cout << "delta is " << delta << " mag is " << delta.get_magnitude() << " sigma " << sigma << std::endl;*/ @@ -200,12 +216,12 @@ << ", " << d.get_coordinate_derivative(1) << ", " << d.get_coordinate_derivative(2) << "]" << std::endl); - internal::SquaredAngstrom t= square(delta[0]); + unit::SquareAngstrom t= square(delta[0]); - internal::SquaredAngstrom magnitude2=square(delta[0])+square(delta[1]) + unit::SquareAngstrom magnitude2=square(delta[0])+square(delta[1]) +square(delta[2]); - //internal::SquaredAngstrom max_motion2= square(max_motion); + //unit::SquaredAngstrom max_motion2= square(max_motion); if (magnitude2 > square(max_change_)) { return false; } @@ -223,85 +239,112 @@ */ Float BrownianDynamics::optimize(unsigned int max_steps) { + IMP_check(get_model() != NULL, "Must set model before calling optimize", + ValueException); setup_particles(); - IMP_LOG(SILENT, "Running brownian dynamics on " << get_particles().size() + IMP_LOG(TERSE, "Running brownian dynamics on " << get_particles().size() << " particles with a step of " << cur_dt_ << std::endl); - unsigned int num_const_dt=0; + setup_particles(); for (unsigned int i = 0; i < max_steps; ++i) { - update_states(); - get_model()->evaluate(true); - if (step()) { - ++num_const_dt; - } + take_step(); + } + return get_model()->evaluate(false); +} - cur_time_= cur_time_+cur_dt_; - if (num_const_dt == 10) { - num_const_dt=0; - cur_dt_= std::min(cur_dt_*internal::Scalar(2), max_dt_); + +void BrownianDynamics::take_step() { + update_states(); + get_model()->evaluate(true); + if (step()) ++ num_const_dt_; + + cur_time_= cur_time_+cur_dt_; + if (num_const_dt_ == 10) { + num_const_dt_=0; + cur_dt_= std::min(cur_dt_*2.0, max_dt_); + } + { + unit::Femtosecond lb(max_dt_/2.0); + std::vector::size_type d=0; + while (lb > cur_dt_) { + lb = lb/2.0; + ++d; } - { - internal::FemtoSecond lb(max_dt_/internal::Scalar(2.0)); - std::vector::size_type d=0; - while (lb > cur_dt_) { - lb = lb/internal::Scalar(2.0); - ++d; - } - time_steps_.resize(std::max(time_steps_.size(), d+1), 0); - ++time_steps_[d]; - } + time_steps_.resize(std::max(time_steps_.size(), d+1), 0); + ++time_steps_[d]; } - return get_model()->evaluate(false); } +void BrownianDynamics::simulate(float max_time) { + IMP_check(get_model() != NULL, "Must set model before calling simulate", + ValueException); + setup_particles(); + unit::Femtosecond mt(max_time); + num_const_dt_=0; + while (cur_time_ < mt) { + take_step(); + } +} -internal::FemtoJoule BrownianDynamics::kt() const +unit::Femtojoule BrownianDynamics::kt() const { - return internal::KB*T_; + return unit::Femtojoule(internal::KB*T_); } -internal::Angstrom -BrownianDynamics::estimate_radius_from_mass_units(Float mass_in_kd) const +unit::Angstrom +BrownianDynamics +::estimate_radius_from_mass_units(unit::Kilodalton mass_in_kd) { - /* Data points used: - thyroglobulin 669kD 85A - catalase 232kD 52A - aldolase 158kD 48A - */ + //unit::KG kg= convert_to_mks(mass_in_kd) + unit::Kilodalton kd(mass_in_kd); + unit::GramPerCubicCentimeter p= unit::GramPerCubicCentimeter(1.410) + + unit::GramPerCubicCentimeter(0.145)* exp(-mass_in_kd.get_value()/13); - internal::Angstrom r(std::pow(mass_in_kd, .3333f)*10); - return r; + unit::Kilogram m= convert_to_mks(kd); + // m/(4/3 pi r^3) =p + // r= (m/p)/(4/(3 pi))^(1/3) + typedef unit::Multiply::type SquareAngstrom; + typedef unit::Multiply::type CubicAngstrom; + + CubicAngstrom v((m/p)*(4.0/(3.0*internal::PI))); + + return unit::Angstrom(std::pow(v.get_value(), .3333)); } -internal::Cm2PerSecond -BrownianDynamics::estimate_D_from_radius(internal::Angstrom r) const +unit::SquareCentimeterPerSecond +BrownianDynamics::estimate_D_from_radius(unit::Angstrom r) const { - MilliPascalSeconds e=eta(T_); - internal::MKSUnit<-13, 0, 1, 0, -1> etar( e*r); + MillipascalSecond e=eta(T_); + //unit::MKSUnit<-13, 0, 1, 0, -1> etar( e*r); /*std::cout << e << " " << etar << " " << kt << std::endl; - std::cout << "scalar etar " << (internal::Scalar(6*internal::PI)*etar) + std::cout << "scalar etar " << (unit::Scalar(6*unit::PI)*etar) << std::endl; - std::cout << "ret pre conv " << (kt/(internal::Scalar(6* internal::PI)*etar)) + std::cout << "ret pre conv " << (kt/(unit::Scalar(6* unit::PI)*etar)) << std::endl;*/ - internal::Cm2PerSecond ret( kt()/(internal::Scalar(6* internal::PI)*etar)); + unit::SquareCentimeterPerSecond ret(kt()/(6.0* internal::PI*e*r)); //std::cout << "ret " << ret << std::endl; return ret; } -internal::Angstrom -BrownianDynamics::compute_sigma_from_D(internal::Cm2PerSecond D) const +unit::Angstrom +BrownianDynamics +::compute_sigma_from_D(unit::SquareCentimeterPerSecond D) const { - return sqrt(internal::Scalar(2.0)*D*cur_dt_); //6*xi*kb*T_; + return sqrt(2.0*D*cur_dt_); //6*xi*kb*T_; } -internal::KCalPerAMol BrownianDynamics -::get_force_scale_from_D(internal::Cm2PerSecond D) const +unit::KilocaloriePerAngstromPerMol BrownianDynamics +::get_force_scale_from_D(unit::SquareCentimeterPerSecond D) const { // force motion is f*cur_dt_*D/kT // sigma*kT/(cur_dt_*D) - return convert_to_kcal(internal::Angstrom(compute_sigma_from_D(D))*kt() - /(max_dt_*D)); + unit::Angstrom s=compute_sigma_from_D(D); + unit::Piconewton pn= s*kt()/(max_dt_*D); + unit::YoctoKilocaloriePerAngstrom yc=unit::convert_J_to_Cal(pn); + return unit::operator*(unit::ATOMS_PER_MOL,yc); } Index: kernel/src/pair_scores/SphereDistancePairScore.cpp =================================================================== --- kernel/src/pair_scores/SphereDistancePairScore.cpp (revision 593) +++ kernel/src/pair_scores/SphereDistancePairScore.cpp (working copy) @@ -23,10 +23,10 @@ { IMP_check(a->has_attribute(radius_), "Particle " << a->get_index() << "missing radius in SphereDistancePairScore", - ValueException("Missing radius")); + ValueException); IMP_check(b->has_attribute(radius_), "Particle " << b->get_index() << "missing radius in SphereDistancePairScore", - ValueException("Missing radius")); + ValueException); Float ra = a->get_value(radius_); Float rb = b->get_value(radius_); return internal::evaluate_distance_pair_score(a,b, da, f_.get(), Index: kernel/src/internal/constants.cpp =================================================================== --- kernel/src/internal/constants.cpp (revision 593) +++ kernel/src/internal/constants.cpp (working copy) @@ -13,16 +13,15 @@ namespace internal { -const ExponentialNumber<23> NA(6.02214179); - +const unit::ExponentialNumber<23> NA(6.02214179); // definition of KB -const MKSUnit<-23, 2, 1, -1, -2> KB(1.3806503); +const unit::Shift::type, -23>::type +KB(1.3806503); -const Kelvin T(297.15); +const unit::Kelvin T(297.15); -const ExponentialNumber<3> JOULES_PER_KCAL(4.1868); -const ExponentialNumber<-4> KCALS_PER_JOULE(2.388459); +const unit::ExponentialNumber<3> JOULES_PER_KILOCALORIE(4.1868); } // namespace internal