/*
* Project: MoleCuilder
* Description: creates and alters molecular systems
* Copyright (C) 2012 University of Bonn. 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 .
*/
/*
* SaturationPotential.cpp
*
* Created on: Oct 11, 2012
* Author: heber
*/
// include config.h
#ifdef HAVE_CONFIG_H
#include
#endif
#include "CodePatterns/MemDebug.hpp"
#include "SaturationPotential.hpp"
#include
#include // for 'map_list_of()'
#include
#include
#include
#include "CodePatterns/Assert.hpp"
#include "CodePatterns/Log.hpp"
#include "FunctionApproximation/Extractors.hpp"
#include "FunctionApproximation/TrainingData.hpp"
#include "Potentials/helpers.hpp"
#include "Potentials/ParticleTypeCheckers.hpp"
class Fragment;
using namespace boost::assign;
// static definitions
const SaturationPotential::ParameterNames_t
SaturationPotential::ParameterNames =
boost::assign::list_of
("all_energy_offset")
("")
("")
("")
("")
("")
;
const std::string SaturationPotential::potential_token("saturation");
SaturationPotential::SaturationPotential(
const ParticleTypes_t &_ParticleTypes) :
EmpiricalPotential(_ParticleTypes),
morse(_ParticleTypes),
angle(addSaturationType(_ParticleTypes)),
energy_offset(0.)
{
// have some decent defaults for parameter_derivative checking
// Morse and Angle have their own defaults, offset is set
ASSERT( _ParticleTypes.size() == (size_t)2,
"SaturationPotential::SaturationPotential() - exactly two types must be given.");
// ASSERT( _ParticleTypes[1] == 1,
// "SaturationPotential::SaturationPotential() - second type must be hydrogen.");
}
SaturationPotential::SaturationPotential(
const ParticleTypes_t &_ParticleTypes,
const double _all_energy_offset,
const double _morse_spring_constant,
const double _morse_equilibrium_distance,
const double _morse_dissociation_energy,
const double _angle_spring_constant,
const double _angle_equilibrium_distance) :
EmpiricalPotential(_ParticleTypes),
morse(_ParticleTypes),
angle(addSaturationType(_ParticleTypes)),
energy_offset(_all_energy_offset)
{
ASSERT( _ParticleTypes.size() == (size_t)2,
"SaturationPotential::SaturationPotential() - exactly two types must be given.");
// ASSERT( _ParticleTypes[1] == 1,
// "SaturationPotential::SaturationPotential() - second type must be hydrogen.");
parameters_t morse_params(morse.getParameterDimension());
morse_params[PairPotential_Morse::spring_constant] = _morse_spring_constant;
morse_params[PairPotential_Morse::equilibrium_distance] = _morse_equilibrium_distance;
morse_params[PairPotential_Morse::dissociation_energy] = _morse_dissociation_energy;
morse_params[PairPotential_Morse::energy_offset] = 0.;
morse.setParameters(morse_params);
parameters_t angle_params(angle.getParameterDimension());
angle_params[PairPotential_Angle::spring_constant] = _angle_spring_constant;
angle_params[PairPotential_Angle::equilibrium_distance] = _angle_equilibrium_distance;
angle_params[PairPotential_Angle::energy_offset] = 0.;
angle.setParameters(angle_params);
}
void SaturationPotential::setParameters(const parameters_t &_params)
{
const size_t paramsDim = _params.size();
ASSERT( paramsDim <= getParameterDimension(),
"SaturationPotential::setParameters() - we need not more than "
+toString(getParameterDimension())+" parameters.");
// LOG(1, "INFO: Setting new SaturationPotential params: " << _params);
// offsets
if (paramsDim > all_energy_offset)
energy_offset = _params[all_energy_offset];
// Morse
{
parameters_t morse_params(morse.getParameters());
if (paramsDim > morse_spring_constant)
morse_params[PairPotential_Morse::spring_constant] = _params[morse_spring_constant];
if (paramsDim > morse_equilibrium_distance)
morse_params[PairPotential_Morse::equilibrium_distance] = _params[morse_equilibrium_distance];
if (paramsDim > morse_dissociation_energy)
morse_params[PairPotential_Morse::dissociation_energy] = _params[morse_dissociation_energy];
morse_params[PairPotential_Morse::energy_offset] = 0.;
morse.setParameters(morse_params);
}
// Angle
{
parameters_t angle_params(angle.getParameters());
if (paramsDim > angle_spring_constant)
angle_params[PairPotential_Angle::spring_constant] = _params[angle_spring_constant];
if (paramsDim > angle_equilibrium_distance)
angle_params[PairPotential_Angle::equilibrium_distance] = _params[angle_equilibrium_distance];
angle_params[PairPotential_Angle::energy_offset] = 0.;
angle.setParameters(angle_params);
}
#ifndef NDEBUG
parameters_t check_params(getParameters());
check_params.resize(paramsDim); // truncate to same size
ASSERT( check_params == _params,
"SaturationPotential::setParameters() - failed, mismatch in to be set "
+toString(_params)+" and set "+toString(check_params)+" params.");
#endif
}
SaturationPotential::parameters_t SaturationPotential::getParameters() const
{
parameters_t params(getParameterDimension());
const parameters_t morse_params = morse.getParameters();
const parameters_t angle_params = angle.getParameters();
params[all_energy_offset] = energy_offset;
params[morse_spring_constant] = morse_params[PairPotential_Morse::spring_constant];
params[morse_equilibrium_distance] = morse_params[PairPotential_Morse::equilibrium_distance];
params[morse_dissociation_energy] = morse_params[PairPotential_Morse::dissociation_energy];
params[angle_spring_constant] = angle_params[PairPotential_Angle::spring_constant];
params[angle_equilibrium_distance] = angle_params[PairPotential_Angle::equilibrium_distance];
return params;
}
void SaturationPotential::stream_to(std::ostream &ost) const
{
morse.stream_to(ost);
ost << std::endl;
angle.stream_to(ost);
}
void SaturationPotential::stream_from(std::istream &ist)
{
morse.stream_from(ist);
ist >> ws;
angle.stream_from(ist);
}
std::vector
triplefunction(
const argument_t &argument,
const FunctionModel::arguments_t& args)
{
const size_t firstindex = argument.indices.first;
const size_t secondindex = argument.indices.second;
// LOG(2, "DEBUG: first index is " << firstindex << ", second index is " << secondindex << ".");
// place all arguments that share either index into a lookup map
typedef std::map< size_t, FunctionModel::arguments_t::const_iterator > IndexLookup_t;
IndexLookup_t LookuptoFirst;
IndexLookup_t LookuptoSecond;
for (FunctionModel::arguments_t::const_iterator iter = args.begin();
iter != args.end();
++iter) {
if (((*iter).indices.first == argument.indices.first)
&& ((*iter).indices.second == argument.indices.second))
continue;
if (firstindex == (*iter).indices.first) {
LookuptoFirst.insert( std::make_pair( (*iter).indices.second, iter) );
}
else if (firstindex == (*iter).indices.second) {
LookuptoFirst.insert( std::make_pair( (*iter).indices.first, iter) );
}
if (secondindex == (*iter).indices.first) {
LookuptoSecond.insert( std::make_pair( (*iter).indices.second, iter) );
}
else if (secondindex == (*iter).indices.second) {
LookuptoSecond.insert( std::make_pair((*iter).indices.first, iter) );
}
}
// {
// std::stringstream lookupstream;
// for (IndexLookup_t::const_iterator iter = LookuptoFirst.begin();
// iter != LookuptoFirst.end();
// ++iter) {
// lookupstream << "(" << iter->first << "," << *(iter->second) << ") ";
// }
// LOG(2, "DEBUG: LookupToFirst is " << lookupstream.str() << ".");
// }
// {
// std::stringstream lookupstream;
// for (IndexLookup_t::const_iterator iter = LookuptoSecond.begin();
// iter != LookuptoSecond.end();
// ++iter) {
// lookupstream << "(" << iter->first << "," << *(iter->second) << ") ";
// }
// LOG(2, "DEBUG: LookuptoSecond is " << lookupstream.str() << ".");
// }
// now go through the first lookup as the second argument and pick the
// corresponding third argument by the matching index
std::vector results;
for (IndexLookup_t::const_iterator iter = LookuptoFirst.begin();
iter != LookuptoFirst.end();
++iter) {
IndexLookup_t::const_iterator otheriter = LookuptoSecond.find(iter->first);
ASSERT( otheriter != LookuptoSecond.end(),
"triplefunction() - cannot find index "+toString(iter->first)
+" in LookupToSecond");
FunctionModel::arguments_t result(1, argument);
result.reserve(3);
result.push_back(*(iter->second));
result.push_back(*(otheriter->second));
results.push_back(result);
}
return results;
}
SaturationPotential::results_t
SaturationPotential::operator()(
const arguments_t &arguments
) const
{
double result = 0.;
const ParticleTypes_t &morse_types = morse.getParticleTypes();
const ParticleTypes_t &angle_types = angle.getParticleTypes();
double multiplicity = 1.;
if ((angle_types[0] == angle_types[1]) && (angle_types[1] == angle_types[2]))
multiplicity = 1./6.;
else if ((angle_types[0] == angle_types[1])
|| (angle_types[1] == angle_types[2])
|| (angle_types[0] == angle_types[2]))
multiplicity = .5;
for(arguments_t::const_iterator argiter = arguments.begin();
argiter != arguments.end();
++argiter) {
const argument_t &r_ij = *argiter;
if (((r_ij.types.first == morse_types[0]) && (r_ij.types.second == morse_types[1]))
|| ((r_ij.types.first == morse_types[1]) && (r_ij.types.second == morse_types[0]))) {
arguments_t args(1, r_ij);
// Morse contribution
const double tmp = morse(args)[0];
// LOG(3, "DEBUG: Morse yields " << tmp << " for << " << r_ij << ".");
result += tmp;
if (result != result)
ELOG(1, "result is NAN.");
}
if (((r_ij.types.first == angle_types[0]) && (r_ij.types.second == angle_types[1]))
|| ((r_ij.types.first == angle_types[1]) && (r_ij.types.second == angle_types[0]))) {
// Angle contribution
{
typedef std::vector tripleargs_t;
tripleargs_t tripleargs =
triplefunction(r_ij, arguments);
for (tripleargs_t::const_iterator iter = tripleargs.begin();
iter != tripleargs.end();
++iter) {
FunctionModel::arguments_t tempargs =
Extractors::reorderArgumentsByParticleTypes(*iter, angle.getParticleTypes());
// We get both angles, e.g. 0-4-1 and 1-4-0, hence multiply with 0.5
const double tmp = multiplicity*angle(tempargs)[0]; // as we have all distances we get both jk and kj
// LOG(3, "DEBUG: angle yields " << tmp << " for << " << tempargs << ".");
result += tmp;
if (result != result)
ELOG(1, "result is NAN.");
}
}
}
}
return std::vector(1, energy_offset + result);
}
SaturationPotential::derivative_components_t
SaturationPotential::derivative(
const arguments_t &arguments
) const
{
ASSERT( 0,
"SaturationPotential::operator() - not implemented.");
derivative_components_t result;
return result;
}
SaturationPotential::results_t
SaturationPotential::parameter_derivative(
const arguments_t &arguments,
const size_t index
) const
{
double result = 0.;
switch (index) {
case all_energy_offset:
result = 1.;
break;
case morse_spring_constant:
case morse_equilibrium_distance:
case morse_dissociation_energy:
{
const ParticleTypes_t &morse_types = morse.getParticleTypes();
for(arguments_t::const_iterator argiter = arguments.begin();
argiter != arguments.end();
++argiter) {
const argument_t &r_ij = *argiter;
if (((r_ij.types.first == morse_types[0]) && (r_ij.types.second == morse_types[1]))
|| ((r_ij.types.first == morse_types[1]) && (r_ij.types.second == morse_types[0]))) {
arguments_t args(1, r_ij);
double tmp = 0.;
switch (index) {
case morse_spring_constant:
tmp += morse.parameter_derivative(args, PairPotential_Morse::spring_constant)[0];
break;
case morse_equilibrium_distance:
tmp += morse.parameter_derivative(args, PairPotential_Morse::equilibrium_distance)[0];
break;
case morse_dissociation_energy:
tmp += morse.parameter_derivative(args, PairPotential_Morse::dissociation_energy)[0];
break;
default:
ASSERT(0, "SaturationPotential::parameter_derivative() - We cannot get here.");
break;
}
// LOG(2, "DEBUG: morse yields " << tmp << " for << " << args << ".");
result += tmp;
}
}
}
break;
case angle_spring_constant:
case angle_equilibrium_distance:
{
const ParticleTypes_t &angle_types = angle.getParticleTypes();
double multiplicity = 1.;
if ((angle_types[0] == angle_types[1]) && (angle_types[1] == angle_types[2]))
multiplicity = 1./6.;
else if ((angle_types[0] == angle_types[1])
|| (angle_types[1] == angle_types[2])
|| (angle_types[0] == angle_types[2]))
multiplicity = .5;
for(arguments_t::const_iterator argiter = arguments.begin();
argiter != arguments.end();
++argiter) {
const argument_t &r_ij = *argiter;
if (((r_ij.types.first == angle_types[0]) && (r_ij.types.second == angle_types[1]))
|| ((r_ij.types.first == angle_types[1]) && (r_ij.types.second == angle_types[0]))) {
typedef std::vector tripleargs_t;
tripleargs_t tripleargs =
triplefunction(r_ij, arguments);
for (tripleargs_t::const_iterator iter = tripleargs.begin();
iter != tripleargs.end();
++iter) {
FunctionModel::arguments_t tempargs =
Extractors::reorderArgumentsByParticleTypes(*iter, angle.getParticleTypes());
// We get both angles, e.g. 0-4-1 and 1-4-0, hence multiply with 0.5
double tmp = 0.;
if (index == angle_spring_constant)
tmp += multiplicity*angle.parameter_derivative(tempargs, PairPotential_Angle::spring_constant)[0];
else if (index == angle_equilibrium_distance)
tmp += multiplicity*angle.parameter_derivative(tempargs, PairPotential_Angle::equilibrium_distance)[0];
// LOG(2, "DEBUG: angle yields " << tmp << " for << " << tempargs << ".");
result += tmp;
if (result != result)
ELOG(1, "result is NAN.");
}
}
}
}
break;
default:
ASSERT( 0, "SaturationPotential::parameter_derivative() - impossible to get here.");
break;
}
return SaturationPotential::results_t(1, result);
}
const SaturationPotential::ParticleTypes_t
SaturationPotential::symmetrizeTypes(const ParticleTypes_t &_ParticleTypes)
{
ASSERT( _ParticleTypes.size() == (size_t)2,
"SaturationPotential::symmetrizeTypes() - require initial _ParticleTypes with two elements.");
// // insert before couple
// ParticleTypes_t types(1, _ParticleTypes[1]);
// types.insert(types.end(), _ParticleTypes.begin(), _ParticleTypes.end());
// insert after the couple
ParticleTypes_t types(_ParticleTypes);
types.push_back( _ParticleTypes.back() );
ASSERT( types.size() == (size_t)3,
"SaturationPotential::symmetrizeTypes() - failed to generate three types for angle.");
return types;
}
const SaturationPotential::ParticleTypes_t
SaturationPotential::addSaturationType(const ParticleTypes_t &_ParticleTypes)
{
ParticleTypes_t types(_ParticleTypes);
types.push_back( ParticleType_t(1) );
return types;
}
FunctionModel::extractor_t
SaturationPotential::getFragmentSpecificExtractor() const
{
// Fragment::charges_t charges;
// charges.resize(getParticleTypes().size());
// std::transform(getParticleTypes().begin(), getParticleTypes().end(),
// charges.begin(), boost::lambda::_1);
FunctionModel::extractor_t returnfunction;
// if (charges[0] == charges[1]) {
// // In case both types are equal there is only a single pair of possible
// // type combinations.
// returnfunction =
// boost::bind(&Extractors::gatherAllDistancesFromFragment,
// boost::bind(&Fragment::getPositions, _1),
// boost::bind(&Fragment::getCharges, _1),
// charges, // is only temporarily created, hence copy
// _2);
// } else {
// we have to chain here a rather complex "tree" of functions
// as we only have a couple of ParticleTypes but need to get
// all possible three pairs of the set of the two types.
// Finally, we also need to arrange them in correct order
// (for PairPotentiale_Angle).
// charges_t firstpair(2, boost::cref(charges[0]));
// only that saturation potential never has its middle element twice!
// hence, we skip the firstpair but keep the code for later generalization
// Fragment::charges_t secondpair(2, boost::cref(charges[1]));
// const Fragment::charges_t &thirdpair = charges;
Fragment::charges_t charges_angle;
{
charges_angle.resize(angle.getParticleTypes().size());
std::transform(angle.getParticleTypes().begin(), angle.getParticleTypes().end(),
charges_angle.begin(), boost::lambda::_1);
}
Fragment::charges_t charges_morse;
{
charges_morse.resize(morse.getParticleTypes().size());
std::transform(morse.getParticleTypes().begin(), morse.getParticleTypes().end(),
charges_morse.begin(), boost::lambda::_1);
}
returnfunction =
// boost::bind(&Extractors::reorderArgumentsByParticleTypes,
boost::bind(&Extractors::combineArguments,
// boost::bind(&Extractors::combineArguments,
boost::bind(&Extractors::gatherAllDistancesFromFragment,
boost::bind(&Fragment::getPositions, _1),
boost::bind(&Fragment::getCharges, _1),
charges_angle, // no crefs here as are temporaries!
_2),
boost::bind(&Extractors::gatherAllDistancesFromFragment,
boost::bind(&Fragment::getPositions, _1),
boost::bind(&Fragment::getCharges, _1),
charges_morse, // no crefs here as are temporaries!
_2)
// )
// boost::bind(&Extractors::gatherAllDistancesFromFragment,
// boost::bind(&Fragment::getPositions, _1),
// boost::bind(&Fragment::getCharges, _1),
// boost::cref(thirdpair), // only the last one is no temporary
// _2)
// ),
// boost::bind(&PairPotential_Angle::getParticleTypes, boost::cref(angle))
);
// }
return returnfunction;
}
void
SaturationPotential::setParametersToRandomInitialValues(
const TrainingData &data)
{
energy_offset = data.getTrainingOutputAverage()[0];
morse.setParametersToRandomInitialValues(data);
angle.setParametersToRandomInitialValues(data);
}