/*
* Project: MoleCuilder
* Description: creates and alters molecular systems
* 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 .
*/
/*
* FourBodyPotential_Torsion.cpp
*
* Created on: Jul 08, 2013
* Author: heber
*/
// include config.h
#ifdef HAVE_CONFIG_H
#include
#endif
#include "CodePatterns/MemDebug.hpp"
#include "FourBodyPotential_Torsion.hpp"
#include // for 'map_list_of()'
#include
#include
#include
#include "CodePatterns/Assert.hpp"
#include "FunctionApproximation/Extractors.hpp"
#include "FunctionApproximation/TrainingData.hpp"
#include "Potentials/helpers.hpp"
#include "Potentials/InternalCoordinates/FourBody_TorsionAngle.hpp"
#include "Potentials/ParticleTypeCheckers.hpp"
#include "RandomNumbers/RandomNumberGeneratorFactory.hpp"
#include "RandomNumbers/RandomNumberGenerator.hpp"
class Fragment;
// static definitions1
const FourBodyPotential_Torsion::ParameterNames_t
FourBodyPotential_Torsion::ParameterNames =
boost::assign::list_of
("spring_constant")
("equilibrium_distance")
;
const std::string FourBodyPotential_Torsion::potential_token("torsion");
Coordinator::ptr FourBodyPotential_Torsion::coordinator(new FourBody_TorsionAngle());
FourBodyPotential_Torsion::FourBodyPotential_Torsion() :
EmpiricalPotential(),
params(parameters_t(MAXPARAMS, 0.))
{
// have some decent defaults for parameter_derivative checking
params[spring_constant] = 1.;
params[equilibrium_distance] = 0.1;
}
FourBodyPotential_Torsion::FourBodyPotential_Torsion(
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] = 0.1;
}
FourBodyPotential_Torsion::FourBodyPotential_Torsion(
const ParticleTypes_t &_ParticleTypes,
const double _spring_constant,
const double _equilibrium_distance) :
EmpiricalPotential(_ParticleTypes),
params(parameters_t(MAXPARAMS, 0.))
{
params[spring_constant] = _spring_constant;
params[equilibrium_distance] = _equilibrium_distance;
}
void FourBodyPotential_Torsion::setParameters(const parameters_t &_params)
{
const size_t paramsDim = _params.size();
ASSERT( paramsDim <= getParameterDimension(),
"FourBodyPotential_Torsion::setParameters() - we need not more than "
+toString(getParameterDimension())+" parameters.");
for(size_t i=0;i(1, result);
}
FourBodyPotential_Torsion::derivative_components_t
FourBodyPotential_Torsion::derivative(
const arguments_t &arguments
) const
{
ASSERT( arguments.size() == getSpecificArgumentCount(),
"FourBodyPotential_Torsion::operator() - requires exactly three arguments.");
ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
arguments, getParticleTypes()),
"FourBodyPotential_Torsion::operator() - types don't match with ones in arguments.");
derivative_components_t result;
const argument_t &r_ij = arguments[0]; // 01
const argument_t &r_ik = arguments[1]; // 02
const argument_t &r_il = arguments[2]; // 03
const argument_t &r_jk = arguments[3]; // 12
const argument_t &r_jl = arguments[4]; // 13
const argument_t &r_kl = arguments[5]; // 23
result.push_back( 2. * params[spring_constant] * ( function_theta(r_ij.distance, r_ik.distance, r_il.distance, r_jk.distance, r_jl.distance, r_kl.distance) - params[equilibrium_distance]) );
ASSERT( result.size() == 1,
"FourBodyPotential_Torsion::operator() - we did not create exactly one component.");
return result;
}
FourBodyPotential_Torsion::results_t
FourBodyPotential_Torsion::parameter_derivative(
const arguments_t &arguments,
const size_t index
) const
{
ASSERT( arguments.size() == getSpecificArgumentCount(),
"FourBodyPotential_Torsion::parameter_derivative() - requires exactly three arguments.");
ASSERT( ParticleTypeChecker::checkArgumentsAgainstParticleTypes(
arguments, getParticleTypes()),
"FourBodyPotential_Torsion::operator() - types don't match with ones in arguments.");
const argument_t &r_ij = arguments[0]; // 01
const argument_t &r_ik = arguments[1]; // 02
const argument_t &r_il = arguments[2]; // 03
const argument_t &r_jk = arguments[3]; // 12
const argument_t &r_jl = arguments[4]; // 13
const argument_t &r_kl = arguments[5]; // 23
switch (index) {
case spring_constant:
{
const result_t result =
Helpers::pow( function_theta(r_ij.distance, r_ik.distance, r_il.distance, r_jk.distance, r_jl.distance, r_kl.distance) - params[equilibrium_distance], 2 );
return std::vector(1, result);
break;
}
case equilibrium_distance:
{
const result_t result =
-2. * params[spring_constant]
* ( function_theta(r_ij.distance, r_ik.distance, r_il.distance, r_jk.distance, r_jl.distance, r_kl.distance) - params[equilibrium_distance]);
return std::vector(1, result);
break;
}
default:
ASSERT(0, "FourBodyPotential_Torsion::parameter_derivative() - derivative to unknown parameter desired.");
break;
}
return std::vector(1);
}
FunctionModel::extractor_t
FourBodyPotential_Torsion::getSpecificExtractor() const
{
Fragment::charges_t charges;
charges.resize(getParticleTypes().size());
std::transform(getParticleTypes().begin(), getParticleTypes().end(),
charges.begin(), boost::lambda::_1);
FunctionModel::extractor_t returnfunction =
boost::bind(&Extractors::gatherDistancesFromFragment,
boost::bind(&Fragment::getPositions, _1),
boost::bind(&Fragment::getCharges, _1),
charges,
_2);
return returnfunction;
}
FunctionModel::filter_t
FourBodyPotential_Torsion::getSpecificFilter() const
{
FunctionModel::filter_t returnfunction =
boost::bind(&Extractors::reorderArgumentsByParticleTypes,
boost::bind(&Extractors::filterArgumentsByParticleTypes,
_1,
getParticleTypes()),
getParticleTypes()
);
return returnfunction;
}
void
FourBodyPotential_Torsion::setParametersToRandomInitialValues(
const TrainingData &data)
{
RandomNumberGenerator &random = RandomNumberGeneratorFactory::getInstance().makeRandomNumberGenerator();
const double rng_min = random.min();
const double rng_max = random.max();
params[FourBodyPotential_Torsion::spring_constant] = 2.*(random()/(rng_max-rng_min));
params[FourBodyPotential_Torsion::equilibrium_distance] = 2.*(random()/(rng_max-rng_min));
}