Hi Daniel,
this is the class I implemented to use the Morse potential as described in Marc's thesis (see http://sgu.bioinfo.cipf.es/martirenom/Thesis/pdf/Thesis.pdf , Chapter 4 page 86):
mscore = De * pow( ( 1.0 - exp( - sqrt( alpha * Kb /De ) * ( feature - Harmonic::get_mean() )/2.0 ) ), 2 );
where: feature = value of feature being tested. alpha = smoothly constant (amplitude of the curve). De =bond energy dissociation (depth of the curve). Kb = K bond Extracted from MacKerell et al., 1998
Davide
On Thu, Jan 8, 2009 at 7:58 PM, Daniel Russel 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) > } > > 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). > > > 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> 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 >>> 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, 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.es web 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.de >>>> >>>> www.tomotronic.org >>>> >>>> >>>> >>>> >>>> >>>> >>>> _______________________________________________ >>>> 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 >>> >>> >>> _______________________________________________ >>> IMP-dev mailing list >>> IMP-dev@salilab.org >>> 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.de >> >> www.tomotronic.org >> >> >> >> >> >> >> _______________________________________________ >> IMP-dev mailing list >> IMP-dev@salilab.org >> 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.org > https://salilab.org/mailman/listinfo/imp-dev > > > > _______________________________________________ > IMP-dev mailing list > IMP-dev@salilab.org > https://salilab.org/mailman/listinfo/imp-dev > >