Daniel,
here is the code for the Morse potential.
Thanks a lot! Davide
/** * \file Morse.h \brief Morse potential function. * * */
#ifndef __IMP_MORSE_H #define __IMP_MORSE_H
#include "Harmonic.h" #include <cmath>
namespace IMP {
//! Morse potential function /** \ingroup unaryf */ class Morse : public Harmonic { public: Morse(Float mean, Float k, Float threshold, Float alpha, Float De) : Harmonic(mean, k) {} virtual ~Morse() {}
//! Calculate Morse potential with respect to the given feature. /** If the feature is greater than threshold, the Morse function is applied. \param[in] feature Value of feature being tested. \param[in] threshold Value for feature being tested. \param[in] smoothly constant alpha (amplitude of the curve). \param[in] bond energy dissociation De (depth of the curve). \return mscore */ virtual Float evaluate(Float feature, Float threshold, Float alpha, Float De) const { Float Kb = 173; // Kcal/mol. K bond Extracted from MacKerell et al., 1998 Float mscore = 0.0;
if ( feature <= threshold ) { return evaluate_with_derivative(feature).first; } else { mscore = De * pow( ( 1.0 - exp( - sqrt( alpha * Kb /De ) * ( feature - Harmonic::get_mean() )/2.0 ) ), 2 ); // MORSE return mscore; } }
//! Calculate harmonic score and derivative with respect to the given feature. /** \param[in] feature Value of feature being tested. \return Score */ virtual FloatPair evaluate_with_derivative(Float feature) const { Float e = (feature - mean_); Float deriv = k_ * e; return std::make_pair(0.5 * k_ * e * e, deriv); }
void show(std::ostream &out=std::cout) const { out << "Morse: " << Harmonic::get_mean() << " and " << Harmonic::get_k() << std::endl; } };
} // namespace IMP
#endif /* __IMP_MORSE_H */
On Mon, Jan 12, 2009 at 2:21 PM, Daniel Russel drussel@gmail.com wrote:
> Sure. Send it to me. Or, better yet, commit a patch in the patches die > which adds it to core. > > > > On Jan 12, 2009, at 7:42 AM, "Davide Baú" davide.bau@gmail.com wrote: > > Hi Daniel, > > would it be possible to check in the Morse potential too, so that to test > if it works fine? > > Cheers, > Davide > > > > On Mon, Jan 12, 2009 at 2:14 AM, Daniel Russel < drussel@gmail.com > drussel@gmail.com> wrote: > >> I checked in a truncated harmonic (IMP.core.TruncatedHarmonicBound, >> TruncatedHarmonicLowerBound, TruncatedHarmonicUpperBound) for people to look >> at. It currently uses Frido's function, but can be easily changed to >> something else. I deliberately kept the documentation vague on the details >> of the internal function as I am not sure we want the exact nature to be >> important (at least until we come up with a function which we have a strong >> justification for). Anyway, let me know if you have comments. >> >> >> On Jan 8, 2009, at 11:11 AM, Friedrich Foerster wrote: >> >> >> >> On Thu, Jan 8, 2009 at 7:58 PM, Daniel Russel < drussel@gmail.com >> drussel@gmail.com> wrote: >> >>> Having such a general function in IMP would be very nice (I was using a >>> non-smooth truncated harmonic before). Such a function should not be >>> dependent on the parameter tuning (since the parameters would simply be >>> passes by a user). >>> Davide, if you have a version of your code which allows the user to pass >>> in the threshold and cutoff, we can put a version into core for people to >>> look at as a starting point. Otherwise, it would only take me a few minutes >>> to throw something together. >>> >>> It sounds like we would like something like: >>> >>> class TruncatedHarmonicLowerBound { >>> /** threshold > center, capped_value > .5*k*(center-threshold)^2 */ >>> TruncatedHarmonicLowerBound(Float center, Float k, Float threshold, >>> Float capped_value) >>> } >>> >> that'd be great. for me that'd help a lot >> >>> >>> If we use Frido's smooth function, capped_value could default to >>> k*(threshold-center)^2, i.e. 2x the value at the threshold to give it space >>> to go smoothly to the value. And we could have Upper and symmetric versions >>> too (as that is little more code). >>> >> also fine with me. >> for the specific function: i just came up with 2-x^(-2) spontaneously. i'd >> also be happy with another solution (morse, fermi, ...) if that is >> advantageous. if anyone has a probabilistic reasoning for a specific >> function, that'd be nice. >> >> thanks >> >> frido >> >>> >>> >>> On Jan 8, 2009, at 3:44 AM, Davide Baú wrote: >>> >>> It's the exp() function of the standard cmath. >>> Performance is indeed another thing to be improved. >>> >>> Davide >>> >>> >>> >>> On Thu, Jan 8, 2009 at 12:37 PM, Friedrich Foerster <foerster@biochem.mpg.de >>> foerster@biochem.mpg.de> wrote: >>> >>>> thanks. >>>> >>>> just for info: isn't the exp() a bit slow, or are you using some dodgy >>>> macro? >>>> >>>> best >>>> >>>> frido >>>> >>>> >>>> On Jan 8, 2009, at 12:05 PM, Davide Baú wrote: >>>> >>>> Hi Frido, >>>>> >>>>> these are the methods I implemented so far in a new class Morse.h. >>>>> >>>>> The first one returns a constant score if feature is greater then >>>>> threshold. At the moment the returned score is the Harmonic potential >>>>> calculated when feature >= threshold (i.e. when feature is >= threshold, the >>>>> returned score will always be Harmonic::evaluate(threshold) ). >>>>> The second one implements the Morse potential as described by Marc. >>>>> >>>>> These methods are still in a development phase since there are some >>>>> parameters which need to be better defined, such as the threshold to be >>>>> applied and the parameters that define the Morse potential. >>>>> >>>>> Regards, >>>>> Davide >>>>> >>>>> >>>>> // 1. a function that is a Harmonic for x<std and constant for x>std >>>>> /** If the feature is greater than or equal to the threshold, the >>>>> score is constant. >>>>> \param[in] feature Value of feature being tested. >>>>> \return Score >>>>> */ >>>>> virtual Float evaluate(Float feature) const { >>>>> Float k = Harmonic::get_k(); >>>>> Float mean = Harmonic::get_mean(); >>>>> Float threshold = 1.0 * mean + XYZRDecorator::get_radius(); >>>>> if ( feature <= threshold ) { >>>>> return Harmonic::evaluate(feature); >>>>> } else { >>>>> return evaluate_with_derivative(threshold, k).first; >>>>> } >>>>> } >>>>> >>>>> >>>>> // 2. a function that is a Harmonic for x<std and asymptotically >>>>> approaches a constant value for x-> inf: >>>>> /** If the feature is greater than or equal to the threshold, apply >>>>> the Morse potential. >>>>> \param[in] feature Value of feature being tested. >>>>> \return Score >>>>> */ >>>>> virtual Float evaluate(Float feature) const { >>>>> Float mean = Harmonic::get_mean(); >>>>> Float k = Harmonic::get_k(); >>>>> Float threshold = 1.0 * mean + XYZRDecorator::get_radius(); >>>>> >>>>> // MORSE POTENTIAL PARAMETERS >>>>> Float var = 1.0/k; >>>>> Float alpha = 2.0; >>>>> Float De = 40.0; >>>>> Float beta = 1.0; >>>>> Float swf = 0.0; >>>>> >>>>> if ( feature <= threshold ) { >>>>> return Harmonic::evaluate(feature); >>>>> } else { >>>>> swf = beta * pow( ( 1.0 - exp( - sqrt( alpha/De ) * ( feature - >>>>> threshold )/2.0 ) ), 2 ); // MORSE POTENTIAL >>>>> return swf; >>>>> } >>>>> } >>>>> >>>>> >>>>> >>>>> >>>>> >>>>> On Thu, Jan 8, 2009 at 10:20 AM, Marc A. Marti-Renom <mmarti@cipf.es >>>>> mmarti@cipf.es> wrote: >>>>> Hi Frido, >>>>> >>>>> Indeed we are also in need of similar stuff. Davide Bau (a postdoc in >>>>> my lab and also in the IMP-Dev list) has already implemented a sort of a >>>>> MORSE potential in IMP restraints, which is not what you ask but may be >>>>> similar. I implemented many many many years ago something similar in CHARMM >>>>> ( http://sgu.bioinfo.cipf.es/martirenom/Thesis/pdf/Thesis.pdf >>>>> http://sgu.bioinfo.cipf.es/martirenom/Thesis/pdf/Thesis.pdf, Chapter 4 >>>>> page 86). >>>>> >>>>> Unfortunately, we are NOT yet ready to commit the changes in the code >>>>> to the SVN since we are in the middle of production phase (ie, making >>>>> models) thus would prefer not to update IMP. But I am sure we could do soon. >>>>> >>>>> Davide, any comments? :-) >>>>> >>>>> m >>>>> >>>>> ------- >>>>> Marc A. Marti-Renom, Head of the Structural Genomics Unit >>>>> Bioinformatics & Genomics Department, Prince Felipe Research Center >>>>> Avda. Autopista del Saler, 16 , 46012 Valencia, Spain >>>>> Tel +34 96 328 96 80 (ext. 1015) Fax +34 96 328 97 01 >>>>> email mmarti@cipf.esmmarti@cipf.es web http://sgu.bioinfo.cipf.es >>>>> http://sgu.bioinfo.cipf.es >>>>> >>>>> *** Consider the environment. Print responsibly. *** >>>>> >>>>> On Jan 8, 2009, at 10:00 AM, Friedrich Foerster wrote: >>>>> >>>>> hi all, >>>>>> >>>>>> when dealing with experimental restraints i often feel that harmonic >>>>>> restraints are somewhat unsatisfactory: >>>>>> experimental data can be wrong and thus there is not necessarily a >>>>>> solution that fulfills all restraints (score=0). in those cases, a model >>>>>> with a lower score is not necessarily any more correct than a higher scoring >>>>>> one. thus, i think it'd be good to have an additional scoring function that >>>>>> levels off if the feature exceeds a certain value. for example, i'd suggest >>>>>> the following two functions alternatively to Harmonic: >>>>>> 1. a function that is a Harmonic for x<std and constant for x>std >>>>>> 2. a function that is a Harmonic for x<std and asymptotically >>>>>> approaches a constant value for x-> inf: >>>>>> e.g.: f=x^2 for x<1 and f=2-x^-2 for x>1 (would even be continuous >>>>>> in 1st derivative) >>>>>> >>>>>> is anyone else also interested in something like that? >>>>>> any better suggestions? >>>>>> is any IMP core expert interested and willing to code that (including >>>>>> analogous functions for Upper and Lower Harmonic)? probably i'd be one of >>>>>> the slowest and clumsiest persons in doing that... >>>>>> >>>>>> thanks >>>>>> >>>>>> frido >>>>>> >>>>>> >>>>>> -- >>>>>> >>>>>> Friedrich Foerster >>>>>> Max-Planck Institut fuer Biochemie >>>>>> Am Klopferspitz 18 >>>>>> D-82152 Martinsried >>>>>> >>>>>> Tel: +49 89 8578 2651 >>>>>> Fax: +49 89 8578 2641 >>>>>> >>>>>> foerster@biochem.mpg.defoerster@biochem.mpg.de >>>>>> >>>>>> http://www.tomotronic.orgwww.tomotronic.org >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> >>>>>> _______________________________________________ >>>>>> IMP-dev mailing list >>>>>> IMP-dev@salilab.orgIMP-dev@salilab.org >>>>>> https://salilab.org/mailman/listinfo/imp-dev >>>>>> https://salilab.org/mailman/listinfo/imp-dev >>>>>> >>>>>> >>>>> >>>>> _______________________________________________ >>>>> IMP-dev mailing list >>>>> IMP-dev@salilab.orgIMP-dev@salilab.org >>>>> https://salilab.org/mailman/listinfo/imp-dev >>>>> https://salilab.org/mailman/listinfo/imp-dev >>>>> >>>>> >>>>> _______________________________________________ >>>>> IMP-dev mailing list >>>>> IMP-dev@salilab.orgIMP-dev@salilab.org >>>>> https://salilab.org/mailman/listinfo/imp-dev >>>>> https://salilab.org/mailman/listinfo/imp-dev >>>>> >>>> >>>> -- >>>> >>>> Friedrich Foerster >>>> Max-Planck Institut fuer Biochemie >>>> Am Klopferspitz 18 >>>> D-82152 Martinsried >>>> >>>> Tel: +49 89 8578 2651 >>>> Fax: +49 89 8578 2641 >>>> >>>> foerster@biochem.mpg.defoerster@biochem.mpg.de >>>> >>>> http://www.tomotronic.orgwww.tomotronic.org >>>> >>>> >>>> >>>> >>>> >>>> >>>> _______________________________________________ >>>> IMP-dev mailing list >>>> IMP-dev@salilab.orgIMP-dev@salilab.org >>>> https://salilab.org/mailman/listinfo/imp-dev >>>> https://salilab.org/mailman/listinfo/imp-dev >>>> >>> >>> >>> ___________ >>> *Carpe diem.* >>> >>> Chi non ride non è una persona seria. >>> Who does not laugh is not a serious person. >>> >>> *"Di solito non sono un uomo religioso, ma se tu sei lassu', salvami, >>> Superman!"* >>> > Homer Simpson >>> >>> _______________________________________________ >>> IMP-dev mailing list >>> IMP-dev@salilab.orgIMP-dev@salilab.org >>> https://salilab.org/mailman/listinfo/imp-dev >>> https://salilab.org/mailman/listinfo/imp-dev >>> >>> >>> >>> _______________________________________________ >>> IMP-dev mailing list >>> IMP-dev@salilab.orgIMP-dev@salilab.org >>> https://salilab.org/mailman/listinfo/imp-dev >>> https://salilab.org/mailman/listinfo/imp-dev >>> >>> >> _______________________________________________ >> IMP-dev mailing list >> IMP-dev@salilab.orgIMP-dev@salilab.org >> https://salilab.org/mailman/listinfo/imp-dev >> https://salilab.org/mailman/listinfo/imp-dev >> >> >> >> _______________________________________________ >> IMP-dev mailing list >> IMP-dev@salilab.orgIMP-dev@salilab.org >> https://salilab.org/mailman/listinfo/imp-dev >> https://salilab.org/mailman/listinfo/imp-ZGV2PC9hPjxicj4KPGJyPjwvYmxvY2txdW9... >> https://salilab.org/mailman/listinfo/imp-dev >> https://salilab.org/mailman/listinfo/imp-ZGV2PC9hPjxicj4KPGJyPjwvYmxvY2txdW9... > > _______________________________________________ > IMP-dev mailing list > IMP-dev@salilab.org > https://salilab.org/mailman/listinfo/imp-dev > > > _______________________________________________ > IMP-dev mailing list > IMP-dev@salilab.org > https://salilab.org/mailman/listinfo/imp-dev > >