/* * 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)); }