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