[4ffbb7] | 1 | /*
|
---|
| 2 | * Project: MoleCuilder
|
---|
| 3 | * Description: creates and alters molecular systems
|
---|
| 4 | * Copyright (C) 2012 University of Bonn. All rights reserved.
|
---|
| 5 | * Please see the COPYING file or "Copyright notice" in builder.cpp for details.
|
---|
| 6 | *
|
---|
| 7 | *
|
---|
| 8 | * This file is part of MoleCuilder.
|
---|
| 9 | *
|
---|
| 10 | * MoleCuilder is free software: you can redistribute it and/or modify
|
---|
| 11 | * it under the terms of the GNU General Public License as published by
|
---|
| 12 | * the Free Software Foundation, either version 2 of the License, or
|
---|
| 13 | * (at your option) any later version.
|
---|
| 14 | *
|
---|
| 15 | * MoleCuilder is distributed in the hope that it will be useful,
|
---|
| 16 | * but WITHOUT ANY WARRANTY; without even the implied warranty of
|
---|
| 17 | * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
|
---|
| 18 | * GNU General Public License for more details.
|
---|
| 19 | *
|
---|
| 20 | * You should have received a copy of the GNU General Public License
|
---|
| 21 | * along with MoleCuilder. If not, see <http://www.gnu.org/licenses/>.
|
---|
| 22 | */
|
---|
| 23 |
|
---|
| 24 | /*
|
---|
| 25 | * SaturationPotential.cpp
|
---|
| 26 | *
|
---|
| 27 | * Created on: Oct 11, 2012
|
---|
| 28 | * Author: heber
|
---|
| 29 | */
|
---|
| 30 |
|
---|
| 31 |
|
---|
| 32 | // include config.h
|
---|
| 33 | #ifdef HAVE_CONFIG_H
|
---|
| 34 | #include <config.h>
|
---|
| 35 | #endif
|
---|
| 36 |
|
---|
| 37 | #include "CodePatterns/MemDebug.hpp"
|
---|
| 38 |
|
---|
| 39 | #include "SaturationPotential.hpp"
|
---|
| 40 |
|
---|
[93e908] | 41 | #include <boost/assign.hpp>
|
---|
[ed2551] | 42 | #include <boost/assign/list_of.hpp> // for 'map_list_of()'
|
---|
[da2d5c] | 43 | #include <boost/lambda/lambda.hpp>
|
---|
[ed2551] | 44 | #include <iostream>
|
---|
| 45 | #include <string>
|
---|
| 46 |
|
---|
[4ffbb7] | 47 | #include "CodePatterns/Assert.hpp"
|
---|
[94f567] | 48 | #include "CodePatterns/Log.hpp"
|
---|
[4ffbb7] | 49 |
|
---|
[7b019a] | 50 | #include "FunctionApproximation/Extractors.hpp"
|
---|
[d52819] | 51 | #include "FunctionApproximation/TrainingData.hpp"
|
---|
[4ffbb7] | 52 | #include "Potentials/helpers.hpp"
|
---|
[b760bc3] | 53 | #include "Potentials/ParticleTypeCheckers.hpp"
|
---|
[4ffbb7] | 54 |
|
---|
[7b019a] | 55 | class Fragment;
|
---|
| 56 |
|
---|
[93e908] | 57 | using namespace boost::assign;
|
---|
| 58 |
|
---|
[ed2551] | 59 | // static definitions
|
---|
| 60 | const SaturationPotential::ParameterNames_t
|
---|
| 61 | SaturationPotential::ParameterNames =
|
---|
| 62 | boost::assign::list_of<std::string>
|
---|
| 63 | ("all_energy_offset")
|
---|
| 64 | ("")
|
---|
| 65 | ("")
|
---|
| 66 | ("")
|
---|
| 67 | ("")
|
---|
| 68 | ("")
|
---|
| 69 | ;
|
---|
| 70 | const std::string SaturationPotential::potential_token("saturation");
|
---|
| 71 |
|
---|
[4ffbb7] | 72 | SaturationPotential::SaturationPotential(
|
---|
[775dd1a] | 73 | const ParticleTypes_t &_ParticleTypes) :
|
---|
[fdd23a] | 74 | EmpiricalPotential(_ParticleTypes),
|
---|
[ed2551] | 75 | morse(_ParticleTypes),
|
---|
[5c3ad5] | 76 | angle(addSaturationType(_ParticleTypes)),
|
---|
[775dd1a] | 77 | energy_offset(0.)
|
---|
[7b019a] | 78 | {
|
---|
| 79 | // have some decent defaults for parameter_derivative checking
|
---|
| 80 | // Morse and Angle have their own defaults, offset is set
|
---|
| 81 | ASSERT( _ParticleTypes.size() == (size_t)2,
|
---|
| 82 | "SaturationPotential::SaturationPotential() - exactly two types must be given.");
|
---|
[5c3ad5] | 83 | // ASSERT( _ParticleTypes[1] == 1,
|
---|
| 84 | // "SaturationPotential::SaturationPotential() - second type must be hydrogen.");
|
---|
[7b019a] | 85 | }
|
---|
[4ffbb7] | 86 |
|
---|
| 87 | SaturationPotential::SaturationPotential(
|
---|
[ed2551] | 88 | const ParticleTypes_t &_ParticleTypes,
|
---|
[b3eabc] | 89 | const double _all_energy_offset,
|
---|
[4ffbb7] | 90 | const double _morse_spring_constant,
|
---|
| 91 | const double _morse_equilibrium_distance,
|
---|
| 92 | const double _morse_dissociation_energy,
|
---|
| 93 | const double _angle_spring_constant,
|
---|
[775dd1a] | 94 | const double _angle_equilibrium_distance) :
|
---|
[fdd23a] | 95 | EmpiricalPotential(_ParticleTypes),
|
---|
[ed2551] | 96 | morse(_ParticleTypes),
|
---|
[5c3ad5] | 97 | angle(addSaturationType(_ParticleTypes)),
|
---|
[775dd1a] | 98 | energy_offset(_all_energy_offset)
|
---|
[4ffbb7] | 99 | {
|
---|
[7b019a] | 100 | ASSERT( _ParticleTypes.size() == (size_t)2,
|
---|
| 101 | "SaturationPotential::SaturationPotential() - exactly two types must be given.");
|
---|
[5c3ad5] | 102 | // ASSERT( _ParticleTypes[1] == 1,
|
---|
| 103 | // "SaturationPotential::SaturationPotential() - second type must be hydrogen.");
|
---|
[4ffbb7] | 104 | parameters_t morse_params(morse.getParameterDimension());
|
---|
| 105 | morse_params[PairPotential_Morse::spring_constant] = _morse_spring_constant;
|
---|
| 106 | morse_params[PairPotential_Morse::equilibrium_distance] = _morse_equilibrium_distance;
|
---|
| 107 | morse_params[PairPotential_Morse::dissociation_energy] = _morse_dissociation_energy;
|
---|
| 108 | morse_params[PairPotential_Morse::energy_offset] = 0.;
|
---|
| 109 | morse.setParameters(morse_params);
|
---|
| 110 | parameters_t angle_params(angle.getParameterDimension());
|
---|
| 111 | angle_params[PairPotential_Angle::spring_constant] = _angle_spring_constant;
|
---|
| 112 | angle_params[PairPotential_Angle::equilibrium_distance] = _angle_equilibrium_distance;
|
---|
| 113 | angle_params[PairPotential_Angle::energy_offset] = 0.;
|
---|
| 114 | angle.setParameters(angle_params);
|
---|
| 115 | }
|
---|
| 116 |
|
---|
| 117 | void SaturationPotential::setParameters(const parameters_t &_params)
|
---|
| 118 | {
|
---|
| 119 | const size_t paramsDim = _params.size();
|
---|
| 120 | ASSERT( paramsDim <= getParameterDimension(),
|
---|
| 121 | "SaturationPotential::setParameters() - we need not more than "
|
---|
| 122 | +toString(getParameterDimension())+" parameters.");
|
---|
| 123 | // LOG(1, "INFO: Setting new SaturationPotential params: " << _params);
|
---|
| 124 |
|
---|
| 125 |
|
---|
| 126 | // offsets
|
---|
| 127 | if (paramsDim > all_energy_offset)
|
---|
| 128 | energy_offset = _params[all_energy_offset];
|
---|
| 129 |
|
---|
| 130 | // Morse
|
---|
| 131 | {
|
---|
| 132 | parameters_t morse_params(morse.getParameters());
|
---|
| 133 | if (paramsDim > morse_spring_constant)
|
---|
| 134 | morse_params[PairPotential_Morse::spring_constant] = _params[morse_spring_constant];
|
---|
| 135 | if (paramsDim > morse_equilibrium_distance)
|
---|
| 136 | morse_params[PairPotential_Morse::equilibrium_distance] = _params[morse_equilibrium_distance];
|
---|
| 137 | if (paramsDim > morse_dissociation_energy)
|
---|
| 138 | morse_params[PairPotential_Morse::dissociation_energy] = _params[morse_dissociation_energy];
|
---|
| 139 | morse_params[PairPotential_Morse::energy_offset] = 0.;
|
---|
| 140 | morse.setParameters(morse_params);
|
---|
| 141 | }
|
---|
| 142 |
|
---|
| 143 | // Angle
|
---|
| 144 | {
|
---|
| 145 | parameters_t angle_params(angle.getParameters());
|
---|
| 146 | if (paramsDim > angle_spring_constant)
|
---|
| 147 | angle_params[PairPotential_Angle::spring_constant] = _params[angle_spring_constant];
|
---|
| 148 | if (paramsDim > angle_equilibrium_distance)
|
---|
| 149 | angle_params[PairPotential_Angle::equilibrium_distance] = _params[angle_equilibrium_distance];
|
---|
| 150 | angle_params[PairPotential_Angle::energy_offset] = 0.;
|
---|
| 151 | angle.setParameters(angle_params);
|
---|
| 152 | }
|
---|
| 153 | #ifndef NDEBUG
|
---|
| 154 | parameters_t check_params(getParameters());
|
---|
| 155 | check_params.resize(paramsDim); // truncate to same size
|
---|
| 156 | ASSERT( check_params == _params,
|
---|
| 157 | "SaturationPotential::setParameters() - failed, mismatch in to be set "
|
---|
| 158 | +toString(_params)+" and set "+toString(check_params)+" params.");
|
---|
| 159 | #endif
|
---|
| 160 | }
|
---|
| 161 |
|
---|
| 162 | SaturationPotential::parameters_t SaturationPotential::getParameters() const
|
---|
| 163 | {
|
---|
| 164 | parameters_t params(getParameterDimension());
|
---|
| 165 | const parameters_t morse_params = morse.getParameters();
|
---|
| 166 | const parameters_t angle_params = angle.getParameters();
|
---|
| 167 |
|
---|
| 168 | params[all_energy_offset] = energy_offset;
|
---|
| 169 |
|
---|
| 170 | params[morse_spring_constant] = morse_params[PairPotential_Morse::spring_constant];
|
---|
| 171 | params[morse_equilibrium_distance] = morse_params[PairPotential_Morse::equilibrium_distance];
|
---|
| 172 | params[morse_dissociation_energy] = morse_params[PairPotential_Morse::dissociation_energy];
|
---|
| 173 |
|
---|
| 174 | params[angle_spring_constant] = angle_params[PairPotential_Angle::spring_constant];
|
---|
| 175 | params[angle_equilibrium_distance] = angle_params[PairPotential_Angle::equilibrium_distance];
|
---|
| 176 | return params;
|
---|
| 177 | }
|
---|
| 178 |
|
---|
[3d2559] | 179 | void SaturationPotential::stream_to(std::ostream &ost) const
|
---|
| 180 | {
|
---|
| 181 | morse.stream_to(ost);
|
---|
| 182 | ost << std::endl;
|
---|
| 183 | angle.stream_to(ost);
|
---|
| 184 | }
|
---|
| 185 |
|
---|
| 186 | void SaturationPotential::stream_from(std::istream &ist)
|
---|
| 187 | {
|
---|
| 188 | morse.stream_from(ist);
|
---|
| 189 | ist >> ws;
|
---|
| 190 | angle.stream_from(ist);
|
---|
| 191 | }
|
---|
| 192 |
|
---|
[0dad5f] | 193 | std::vector<FunctionModel::arguments_t>
|
---|
| 194 | triplefunction(
|
---|
| 195 | const argument_t &argument,
|
---|
| 196 | const FunctionModel::arguments_t& args)
|
---|
| 197 | {
|
---|
| 198 | const size_t firstindex = argument.indices.first;
|
---|
| 199 | const size_t secondindex = argument.indices.second;
|
---|
| 200 | // LOG(2, "DEBUG: first index is " << firstindex << ", second index is " << secondindex << ".");
|
---|
| 201 |
|
---|
| 202 | // place all arguments that share either index into a lookup map
|
---|
| 203 | typedef std::map< size_t, FunctionModel::arguments_t::const_iterator > IndexLookup_t;
|
---|
| 204 | IndexLookup_t LookuptoFirst;
|
---|
| 205 | IndexLookup_t LookuptoSecond;
|
---|
| 206 | for (FunctionModel::arguments_t::const_iterator iter = args.begin();
|
---|
| 207 | iter != args.end();
|
---|
| 208 | ++iter) {
|
---|
| 209 | if (((*iter).indices.first == argument.indices.first)
|
---|
| 210 | && ((*iter).indices.second == argument.indices.second))
|
---|
| 211 | continue;
|
---|
| 212 | if (firstindex == (*iter).indices.first) {
|
---|
| 213 | LookuptoFirst.insert( std::make_pair( (*iter).indices.second, iter) );
|
---|
| 214 | }
|
---|
| 215 | else if (firstindex == (*iter).indices.second) {
|
---|
| 216 | LookuptoFirst.insert( std::make_pair( (*iter).indices.first, iter) );
|
---|
| 217 | }
|
---|
| 218 | if (secondindex == (*iter).indices.first) {
|
---|
| 219 | LookuptoSecond.insert( std::make_pair( (*iter).indices.second, iter) );
|
---|
| 220 | }
|
---|
| 221 | else if (secondindex == (*iter).indices.second) {
|
---|
| 222 | LookuptoSecond.insert( std::make_pair((*iter).indices.first, iter) );
|
---|
| 223 | }
|
---|
| 224 | }
|
---|
| 225 | // {
|
---|
| 226 | // std::stringstream lookupstream;
|
---|
| 227 | // for (IndexLookup_t::const_iterator iter = LookuptoFirst.begin();
|
---|
| 228 | // iter != LookuptoFirst.end();
|
---|
| 229 | // ++iter) {
|
---|
| 230 | // lookupstream << "(" << iter->first << "," << *(iter->second) << ") ";
|
---|
| 231 | // }
|
---|
| 232 | // LOG(2, "DEBUG: LookupToFirst is " << lookupstream.str() << ".");
|
---|
| 233 | // }
|
---|
| 234 | // {
|
---|
| 235 | // std::stringstream lookupstream;
|
---|
| 236 | // for (IndexLookup_t::const_iterator iter = LookuptoSecond.begin();
|
---|
| 237 | // iter != LookuptoSecond.end();
|
---|
| 238 | // ++iter) {
|
---|
| 239 | // lookupstream << "(" << iter->first << "," << *(iter->second) << ") ";
|
---|
| 240 | // }
|
---|
| 241 | // LOG(2, "DEBUG: LookuptoSecond is " << lookupstream.str() << ".");
|
---|
| 242 | // }
|
---|
| 243 |
|
---|
| 244 | // now go through the first lookup as the second argument and pick the
|
---|
| 245 | // corresponding third argument by the matching index
|
---|
| 246 | std::vector<FunctionModel::arguments_t> results;
|
---|
| 247 | for (IndexLookup_t::const_iterator iter = LookuptoFirst.begin();
|
---|
| 248 | iter != LookuptoFirst.end();
|
---|
| 249 | ++iter) {
|
---|
| 250 | IndexLookup_t::const_iterator otheriter = LookuptoSecond.find(iter->first);
|
---|
| 251 | ASSERT( otheriter != LookuptoSecond.end(),
|
---|
| 252 | "triplefunction() - cannot find index "+toString(iter->first)
|
---|
| 253 | +" in LookupToSecond");
|
---|
| 254 | FunctionModel::arguments_t result(1, argument);
|
---|
| 255 | result.reserve(3);
|
---|
| 256 | result.push_back(*(iter->second));
|
---|
| 257 | result.push_back(*(otheriter->second));
|
---|
| 258 | results.push_back(result);
|
---|
| 259 | }
|
---|
| 260 |
|
---|
| 261 | return results;
|
---|
| 262 | }
|
---|
| 263 |
|
---|
[4ffbb7] | 264 | SaturationPotential::results_t
|
---|
| 265 | SaturationPotential::operator()(
|
---|
| 266 | const arguments_t &arguments
|
---|
| 267 | ) const
|
---|
| 268 | {
|
---|
| 269 | double result = 0.;
|
---|
[b760bc3] | 270 | const ParticleTypes_t &morse_types = morse.getParticleTypes();
|
---|
[5c3ad5] | 271 | const ParticleTypes_t &angle_types = angle.getParticleTypes();
|
---|
| 272 | double multiplicity = 1.;
|
---|
| 273 | if ((angle_types[0] == angle_types[1]) && (angle_types[1] == angle_types[2]))
|
---|
| 274 | multiplicity = 1./6.;
|
---|
| 275 | else if ((angle_types[0] == angle_types[1])
|
---|
| 276 | || (angle_types[1] == angle_types[2])
|
---|
| 277 | || (angle_types[0] == angle_types[2]))
|
---|
| 278 | multiplicity = .5;
|
---|
[4ffbb7] | 279 | for(arguments_t::const_iterator argiter = arguments.begin();
|
---|
| 280 | argiter != arguments.end();
|
---|
| 281 | ++argiter) {
|
---|
| 282 | const argument_t &r_ij = *argiter;
|
---|
[b760bc3] | 283 | if (((r_ij.types.first == morse_types[0]) && (r_ij.types.second == morse_types[1]))
|
---|
| 284 | || ((r_ij.types.first == morse_types[1]) && (r_ij.types.second == morse_types[0]))) {
|
---|
[2ba2ed] | 285 | arguments_t args(1, r_ij);
|
---|
| 286 |
|
---|
| 287 | // Morse contribution
|
---|
[b760bc3] | 288 | const double tmp = morse(args)[0];
|
---|
[a82127] | 289 | // LOG(3, "DEBUG: Morse yields " << tmp << " for << " << r_ij << ".");
|
---|
[b760bc3] | 290 | result += tmp;
|
---|
[94f567] | 291 | if (result != result)
|
---|
| 292 | ELOG(1, "result is NAN.");
|
---|
[5c3ad5] | 293 | }
|
---|
| 294 | if (((r_ij.types.first == angle_types[0]) && (r_ij.types.second == angle_types[1]))
|
---|
| 295 | || ((r_ij.types.first == angle_types[1]) && (r_ij.types.second == angle_types[0]))) {
|
---|
[0dad5f] | 296 | // Angle contribution
|
---|
| 297 | {
|
---|
| 298 | typedef std::vector<FunctionModel::arguments_t> tripleargs_t;
|
---|
| 299 | tripleargs_t tripleargs =
|
---|
| 300 | triplefunction(r_ij, arguments);
|
---|
| 301 | for (tripleargs_t::const_iterator iter = tripleargs.begin();
|
---|
| 302 | iter != tripleargs.end();
|
---|
| 303 | ++iter) {
|
---|
| 304 | FunctionModel::arguments_t tempargs =
|
---|
| 305 | Extractors::reorderArgumentsByParticleTypes(*iter, angle.getParticleTypes());
|
---|
[a82127] | 306 | // We get both angles, e.g. 0-4-1 and 1-4-0, hence multiply with 0.5
|
---|
[5c3ad5] | 307 | const double tmp = multiplicity*angle(tempargs)[0]; // as we have all distances we get both jk and kj
|
---|
[a82127] | 308 | // LOG(3, "DEBUG: angle yields " << tmp << " for << " << tempargs << ".");
|
---|
[0dad5f] | 309 | result += tmp;
|
---|
| 310 | if (result != result)
|
---|
| 311 | ELOG(1, "result is NAN.");
|
---|
| 312 | }
|
---|
| 313 | }
|
---|
[4ffbb7] | 314 | }
|
---|
| 315 | }
|
---|
| 316 | return std::vector<result_t>(1, energy_offset + result);
|
---|
| 317 | }
|
---|
| 318 |
|
---|
| 319 | SaturationPotential::derivative_components_t
|
---|
| 320 | SaturationPotential::derivative(
|
---|
| 321 | const arguments_t &arguments
|
---|
| 322 | ) const
|
---|
| 323 | {
|
---|
| 324 | ASSERT( 0,
|
---|
| 325 | "SaturationPotential::operator() - not implemented.");
|
---|
| 326 | derivative_components_t result;
|
---|
| 327 | return result;
|
---|
| 328 | }
|
---|
| 329 |
|
---|
| 330 | SaturationPotential::results_t
|
---|
| 331 | SaturationPotential::parameter_derivative(
|
---|
| 332 | const arguments_t &arguments,
|
---|
| 333 | const size_t index
|
---|
| 334 | ) const
|
---|
| 335 | {
|
---|
| 336 | double result = 0.;
|
---|
[5c3ad5] | 337 | switch (index) {
|
---|
| 338 | case all_energy_offset:
|
---|
[0dad5f] | 339 | result = 1.;
|
---|
[5c3ad5] | 340 | break;
|
---|
| 341 | case morse_spring_constant:
|
---|
| 342 | case morse_equilibrium_distance:
|
---|
| 343 | case morse_dissociation_energy:
|
---|
| 344 | {
|
---|
| 345 | const ParticleTypes_t &morse_types = morse.getParticleTypes();
|
---|
| 346 | for(arguments_t::const_iterator argiter = arguments.begin();
|
---|
| 347 | argiter != arguments.end();
|
---|
| 348 | ++argiter) {
|
---|
| 349 | const argument_t &r_ij = *argiter;
|
---|
| 350 | if (((r_ij.types.first == morse_types[0]) && (r_ij.types.second == morse_types[1]))
|
---|
| 351 | || ((r_ij.types.first == morse_types[1]) && (r_ij.types.second == morse_types[0]))) {
|
---|
| 352 | arguments_t args(1, r_ij);
|
---|
| 353 |
|
---|
| 354 | double tmp = 0.;
|
---|
| 355 | switch (index) {
|
---|
| 356 | case morse_spring_constant:
|
---|
| 357 | tmp += morse.parameter_derivative(args, PairPotential_Morse::spring_constant)[0];
|
---|
| 358 | break;
|
---|
| 359 | case morse_equilibrium_distance:
|
---|
| 360 | tmp += morse.parameter_derivative(args, PairPotential_Morse::equilibrium_distance)[0];
|
---|
| 361 | break;
|
---|
| 362 | case morse_dissociation_energy:
|
---|
| 363 | tmp += morse.parameter_derivative(args, PairPotential_Morse::dissociation_energy)[0];
|
---|
| 364 | break;
|
---|
| 365 | default:
|
---|
| 366 | ASSERT(0, "SaturationPotential::parameter_derivative() - We cannot get here.");
|
---|
| 367 | break;
|
---|
| 368 | }
|
---|
| 369 | // LOG(2, "DEBUG: morse yields " << tmp << " for << " << args << ".");
|
---|
| 370 | result += tmp;
|
---|
| 371 | }
|
---|
| 372 | }
|
---|
| 373 | }
|
---|
| 374 | break;
|
---|
| 375 | case angle_spring_constant:
|
---|
| 376 | case angle_equilibrium_distance:
|
---|
| 377 | {
|
---|
| 378 | const ParticleTypes_t &angle_types = angle.getParticleTypes();
|
---|
| 379 | double multiplicity = 1.;
|
---|
| 380 | if ((angle_types[0] == angle_types[1]) && (angle_types[1] == angle_types[2]))
|
---|
| 381 | multiplicity = 1./6.;
|
---|
| 382 | else if ((angle_types[0] == angle_types[1])
|
---|
| 383 | || (angle_types[1] == angle_types[2])
|
---|
| 384 | || (angle_types[0] == angle_types[2]))
|
---|
| 385 | multiplicity = .5;
|
---|
| 386 | for(arguments_t::const_iterator argiter = arguments.begin();
|
---|
| 387 | argiter != arguments.end();
|
---|
| 388 | ++argiter) {
|
---|
| 389 | const argument_t &r_ij = *argiter;
|
---|
| 390 | if (((r_ij.types.first == angle_types[0]) && (r_ij.types.second == angle_types[1]))
|
---|
| 391 | || ((r_ij.types.first == angle_types[1]) && (r_ij.types.second == angle_types[0]))) {
|
---|
| 392 | typedef std::vector<FunctionModel::arguments_t> tripleargs_t;
|
---|
| 393 | tripleargs_t tripleargs =
|
---|
| 394 | triplefunction(r_ij, arguments);
|
---|
| 395 | for (tripleargs_t::const_iterator iter = tripleargs.begin();
|
---|
| 396 | iter != tripleargs.end();
|
---|
| 397 | ++iter) {
|
---|
| 398 | FunctionModel::arguments_t tempargs =
|
---|
| 399 | Extractors::reorderArgumentsByParticleTypes(*iter, angle.getParticleTypes());
|
---|
| 400 | // We get both angles, e.g. 0-4-1 and 1-4-0, hence multiply with 0.5
|
---|
| 401 | double tmp = 0.;
|
---|
| 402 | if (index == angle_spring_constant)
|
---|
| 403 | tmp += multiplicity*angle.parameter_derivative(tempargs, PairPotential_Angle::spring_constant)[0];
|
---|
| 404 | else if (index == angle_equilibrium_distance)
|
---|
| 405 | tmp += multiplicity*angle.parameter_derivative(tempargs, PairPotential_Angle::equilibrium_distance)[0];
|
---|
| 406 | // LOG(2, "DEBUG: angle yields " << tmp << " for << " << tempargs << ".");
|
---|
| 407 | result += tmp;
|
---|
| 408 | if (result != result)
|
---|
| 409 | ELOG(1, "result is NAN.");
|
---|
[4ffbb7] | 410 | }
|
---|
| 411 | }
|
---|
| 412 | }
|
---|
[b3eabc] | 413 | }
|
---|
[5c3ad5] | 414 | break;
|
---|
| 415 | default:
|
---|
| 416 | ASSERT( 0, "SaturationPotential::parameter_derivative() - impossible to get here.");
|
---|
| 417 | break;
|
---|
[4ffbb7] | 418 | }
|
---|
| 419 | return SaturationPotential::results_t(1, result);
|
---|
| 420 | }
|
---|
[ed2551] | 421 |
|
---|
[93e908] | 422 | const SaturationPotential::ParticleTypes_t
|
---|
| 423 | SaturationPotential::symmetrizeTypes(const ParticleTypes_t &_ParticleTypes)
|
---|
| 424 | {
|
---|
| 425 | ASSERT( _ParticleTypes.size() == (size_t)2,
|
---|
| 426 | "SaturationPotential::symmetrizeTypes() - require initial _ParticleTypes with two elements.");
|
---|
| 427 | // // insert before couple
|
---|
| 428 | // ParticleTypes_t types(1, _ParticleTypes[1]);
|
---|
| 429 | // types.insert(types.end(), _ParticleTypes.begin(), _ParticleTypes.end());
|
---|
| 430 | // insert after the couple
|
---|
| 431 | ParticleTypes_t types(_ParticleTypes);
|
---|
| 432 | types.push_back( _ParticleTypes.back() );
|
---|
| 433 | ASSERT( types.size() == (size_t)3,
|
---|
| 434 | "SaturationPotential::symmetrizeTypes() - failed to generate three types for angle.");
|
---|
| 435 | return types;
|
---|
| 436 | }
|
---|
| 437 |
|
---|
[5c3ad5] | 438 | const SaturationPotential::ParticleTypes_t
|
---|
| 439 | SaturationPotential::addSaturationType(const ParticleTypes_t &_ParticleTypes)
|
---|
| 440 | {
|
---|
| 441 | ParticleTypes_t types(_ParticleTypes);
|
---|
| 442 | types.push_back( ParticleType_t(1) );
|
---|
| 443 | return types;
|
---|
| 444 | }
|
---|
| 445 |
|
---|
[7b019a] | 446 | FunctionModel::extractor_t
|
---|
[da2d5c] | 447 | SaturationPotential::getFragmentSpecificExtractor() const
|
---|
[7b019a] | 448 | {
|
---|
[5c3ad5] | 449 | // Fragment::charges_t charges;
|
---|
| 450 | // charges.resize(getParticleTypes().size());
|
---|
| 451 | // std::transform(getParticleTypes().begin(), getParticleTypes().end(),
|
---|
| 452 | // charges.begin(), boost::lambda::_1);
|
---|
[7b019a] | 453 | FunctionModel::extractor_t returnfunction;
|
---|
[5c3ad5] | 454 | // if (charges[0] == charges[1]) {
|
---|
| 455 | // // In case both types are equal there is only a single pair of possible
|
---|
| 456 | // // type combinations.
|
---|
| 457 | // returnfunction =
|
---|
| 458 | // boost::bind(&Extractors::gatherAllDistancesFromFragment,
|
---|
| 459 | // boost::bind(&Fragment::getPositions, _1),
|
---|
| 460 | // boost::bind(&Fragment::getCharges, _1),
|
---|
| 461 | // charges, // is only temporarily created, hence copy
|
---|
| 462 | // _2);
|
---|
| 463 | // } else {
|
---|
[7b019a] | 464 | // we have to chain here a rather complex "tree" of functions
|
---|
| 465 | // as we only have a couple of ParticleTypes but need to get
|
---|
| 466 | // all possible three pairs of the set of the two types.
|
---|
| 467 | // Finally, we also need to arrange them in correct order
|
---|
| 468 | // (for PairPotentiale_Angle).
|
---|
[5c3ad5] | 469 | // charges_t firstpair(2, boost::cref(charges[0]));
|
---|
[0dad5f] | 470 | // only that saturation potential never has its middle element twice!
|
---|
| 471 | // hence, we skip the firstpair but keep the code for later generalization
|
---|
[5c3ad5] | 472 | // Fragment::charges_t secondpair(2, boost::cref(charges[1]));
|
---|
| 473 | // const Fragment::charges_t &thirdpair = charges;
|
---|
| 474 | Fragment::charges_t charges_angle;
|
---|
| 475 | {
|
---|
| 476 | charges_angle.resize(angle.getParticleTypes().size());
|
---|
| 477 | std::transform(angle.getParticleTypes().begin(), angle.getParticleTypes().end(),
|
---|
| 478 | charges_angle.begin(), boost::lambda::_1);
|
---|
| 479 | }
|
---|
| 480 | Fragment::charges_t charges_morse;
|
---|
| 481 | {
|
---|
| 482 | charges_morse.resize(morse.getParticleTypes().size());
|
---|
| 483 | std::transform(morse.getParticleTypes().begin(), morse.getParticleTypes().end(),
|
---|
| 484 | charges_morse.begin(), boost::lambda::_1);
|
---|
| 485 | }
|
---|
[7b019a] | 486 | returnfunction =
|
---|
[0dad5f] | 487 | // boost::bind(&Extractors::reorderArgumentsByParticleTypes,
|
---|
[7b019a] | 488 | boost::bind(&Extractors::combineArguments,
|
---|
[0dad5f] | 489 | // boost::bind(&Extractors::combineArguments,
|
---|
[7b019a] | 490 | boost::bind(&Extractors::gatherAllDistancesFromFragment,
|
---|
| 491 | boost::bind(&Fragment::getPositions, _1),
|
---|
| 492 | boost::bind(&Fragment::getCharges, _1),
|
---|
[5c3ad5] | 493 | charges_angle, // no crefs here as are temporaries!
|
---|
[0dad5f] | 494 | _2),
|
---|
[5c3ad5] | 495 | boost::bind(&Extractors::gatherAllDistancesFromFragment,
|
---|
| 496 | boost::bind(&Fragment::getPositions, _1),
|
---|
| 497 | boost::bind(&Fragment::getCharges, _1),
|
---|
| 498 | charges_morse, // no crefs here as are temporaries!
|
---|
| 499 | _2)
|
---|
| 500 | // )
|
---|
| 501 | // boost::bind(&Extractors::gatherAllDistancesFromFragment,
|
---|
| 502 | // boost::bind(&Fragment::getPositions, _1),
|
---|
| 503 | // boost::bind(&Fragment::getCharges, _1),
|
---|
| 504 | // boost::cref(thirdpair), // only the last one is no temporary
|
---|
| 505 | // _2)
|
---|
| 506 | // ),
|
---|
| 507 | // boost::bind(&PairPotential_Angle::getParticleTypes, boost::cref(angle))
|
---|
| 508 | );
|
---|
| 509 | // }
|
---|
[7b019a] | 510 | return returnfunction;
|
---|
| 511 | }
|
---|
[d52819] | 512 |
|
---|
| 513 | void
|
---|
| 514 | SaturationPotential::setParametersToRandomInitialValues(
|
---|
| 515 | const TrainingData &data)
|
---|
| 516 | {
|
---|
| 517 | energy_offset = data.getTrainingOutputAverage()[0];
|
---|
| 518 | morse.setParametersToRandomInitialValues(data);
|
---|
| 519 | angle.setParametersToRandomInitialValues(data);
|
---|
| 520 | }
|
---|