/*
* Project: MoleCuilder
* Description: creates and alters molecular systems
* Copyright (C) 2012 University of Bonn. All rights reserved.
* Copyright (C) 2013 Frederik Heber. All rights reserved.
* Please see the COPYING file or "Copyright notice" in builder.cpp for details.
*
*
* This file is part of MoleCuilder.
*
* MoleCuilder is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 2 of the License, or
* (at your option) any later version.
*
* MoleCuilder is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with MoleCuilder. If not, see .
*/
/*
* PairPotential_Morse.cpp
*
* Created on: Oct 03, 2012
* Author: heber
*/
// include config.h
#ifdef HAVE_CONFIG_H
#include
#endif
#include "CodePatterns/MemDebug.hpp"
#include "PairPotential_Morse.hpp"
#include // for 'map_list_of()'
#include
#include
#include
#include
#include "CodePatterns/Assert.hpp"
#include "FunctionApproximation/Extractors.hpp"
#include "FunctionApproximation/TrainingData.hpp"
#include "Potentials/helpers.hpp"
#include "Potentials/InternalCoordinates/TwoBody_Length.hpp"
#include "Potentials/ParticleTypeCheckers.hpp"
#include "RandomNumbers/RandomNumberGeneratorFactory.hpp"
#include "RandomNumbers/RandomNumberGenerator.hpp"
class Fragment;
// static definitions
const PairPotential_Morse::ParameterNames_t
PairPotential_Morse::ParameterNames =
boost::assign::list_of
("spring_constant")
("equilibrium_distance")
("dissociation_energy")
;
const std::string PairPotential_Morse::potential_token("morse");
Coordinator::ptr PairPotential_Morse::coordinator(Memory::ignore(new TwoBody_Length()));
PairPotential_Morse::PairPotential_Morse() :
EmpiricalPotential(),
params(parameters_t(MAXPARAMS, 0.))
{
// have some decent defaults for parameter_derivative checking
params[spring_constant] = 1.;
params[equilibrium_distance] = 1.;
params[dissociation_energy] = 0.1;
}
PairPotential_Morse::PairPotential_Morse(
const ParticleTypes_t &_ParticleTypes
) :
EmpiricalPotential(_ParticleTypes),
params(parameters_t(MAXPARAMS, 0.))
{
// have some decent defaults for parameter_derivative checking
params[spring_constant] = 1.;
params[equilibrium_distance] = 1.;
params[dissociation_energy] = 0.1;
}
PairPotential_Morse::PairPotential_Morse(
const ParticleTypes_t &_ParticleTypes,
const double _spring_constant,
const double _equilibrium_distance,
const double _dissociation_energy) :
EmpiricalPotential(_ParticleTypes),
params(parameters_t(MAXPARAMS, 0.))
{
params[spring_constant] = _spring_constant;
params[equilibrium_distance] = _equilibrium_distance;
params[dissociation_energy] = _dissociation_energy;
}
void PairPotential_Morse::setParameters(const parameters_t &_params)
{
const size_t paramsDim = _params.size();
ASSERT( paramsDim <= getParameterDimension(),
"PairPotential_Morse::setParameters() - we need not more than "
+toString(getParameterDimension())+" parameters.");
for(size_t i=0;i(1, result);
}
PairPotential_Morse::derivative_components_t
PairPotential_Morse::derivative(
const list_of_arguments_t &listarguments
) const
{
result_t result = 0.;
for(list_of_arguments_t::const_iterator iter = listarguments.begin();
iter != listarguments.end(); ++iter) {
const arguments_t &arguments = *iter;
ASSERT( arguments.size() == 1,
"PairPotential_Morse::operator() - requires exactly one argument.");
ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
arguments, getParticleTypes()),
"PairPotential_Morse::operator() - types don't match with ones in arguments.");
const argument_t &r_ij = arguments[0];
// Maple result: 2*D*(1-exp(-k*(r-R)))*k*exp(-k*(r-R))
result +=
2. * params[dissociation_energy]
* ( 1. - exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance])))
* (- params[spring_constant]) * exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance]));
}
return derivative_components_t(1, result);
}
PairPotential_Morse::results_t
PairPotential_Morse::parameter_derivative(
const list_of_arguments_t &listarguments,
const size_t index
) const
{
result_t result = 0.;
for(list_of_arguments_t::const_iterator iter = listarguments.begin();
iter != listarguments.end(); ++iter) {
const arguments_t &arguments = *iter;
ASSERT( arguments.size() == 1,
"PairPotential_Morse::parameter_derivative() - requires exactly one argument.");
ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
arguments, getParticleTypes()),
"PairPotential_Morse::operator() - types don't match with ones in arguments.");
const argument_t &r_ij = arguments[0];
switch (index) {
case spring_constant:
{
// Maple result: -2*D*(1-exp(-k*(r-R)))*(-r+R)*exp(-k*(r-R))
result +=
- 2. * params[dissociation_energy]
* ( 1. - exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance])))
* (- r_ij.distance + params[equilibrium_distance])
* exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance]))
;
break;
}
case equilibrium_distance:
{
// Maple result: -2*D*(1-exp(-k*(r-R)))*k*exp(-k*(r-R))
result +=
- 2. * params[dissociation_energy]
* ( 1. - exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance])))
* params[spring_constant] * exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance]))
;
break;
}
case dissociation_energy:
{
// Maple result: (1-exp(-k*(r-R)))^2
result +=
Helpers::pow(1.
- exp( - params[spring_constant] * ( r_ij.distance - params[equilibrium_distance])),2);
break;
}
default:
ASSERT(0, "PairPotential_Morse::parameter_derivative() - derivative to unknown parameter desired.");
break;
}
}
return results_t(1, result);
}
FunctionModel::filter_t PairPotential_Morse::getSpecificFilter() const
{
FunctionModel::filter_t returnfunction =
boost::bind(&Extractors::filterArgumentsByParticleTypes,
_1,
getParticleTypes());
return returnfunction;
}
void
PairPotential_Morse::setParametersToRandomInitialValues(
const TrainingData &data)
{
RandomNumberGenerator &random = RandomNumberGeneratorFactory::getInstance().makeRandomNumberGenerator();
const double rng_min = random.min();
const double rng_max = random.max();
params[PairPotential_Morse::dissociation_energy] = 1e+0*(random()/(rng_max-rng_min));// 0.5;
params[PairPotential_Morse::spring_constant] = 1e+0*(random()/(rng_max-rng_min));// 1.;
params[PairPotential_Morse::equilibrium_distance] = 3e+0*(random()/(rng_max-rng_min));//2.9;
}