source: src/LevMartester.cpp@ ed2551

Action_Thermostats Add_AtomRandomPerturbation Add_FitFragmentPartialChargesAction Add_RotateAroundBondAction Add_SelectAtomByNameAction Added_ParseSaveFragmentResults AddingActions_SaveParseParticleParameters Adding_Graph_to_ChangeBondActions Adding_MD_integration_tests Adding_ParticleName_to_Atom Adding_StructOpt_integration_tests AtomFragments Automaking_mpqc_open AutomationFragmentation_failures Candidate_v1.5.4 Candidate_v1.6.0 Candidate_v1.6.1 ChangeBugEmailaddress ChangingTestPorts ChemicalSpaceEvaluator CombiningParticlePotentialParsing Combining_Subpackages Debian_Package_split Debian_package_split_molecuildergui_only Disabling_MemDebug Docu_Python_wait EmpiricalPotential_contain_HomologyGraph EmpiricalPotential_contain_HomologyGraph_documentation Enable_parallel_make_install Enhance_userguide Enhanced_StructuralOptimization Enhanced_StructuralOptimization_continued Example_ManyWaysToTranslateAtom Exclude_Hydrogens_annealWithBondGraph FitPartialCharges_GlobalError Fix_BoundInBox_CenterInBox_MoleculeActions Fix_ChargeSampling_PBC Fix_ChronosMutex Fix_FitPartialCharges Fix_FitPotential_needs_atomicnumbers Fix_ForceAnnealing Fix_IndependentFragmentGrids Fix_ParseParticles Fix_ParseParticles_split_forward_backward_Actions Fix_PopActions Fix_QtFragmentList_sorted_selection Fix_Restrictedkeyset_FragmentMolecule Fix_StatusMsg Fix_StepWorldTime_single_argument Fix_Verbose_Codepatterns Fix_fitting_potentials Fixes ForceAnnealing_goodresults ForceAnnealing_oldresults ForceAnnealing_tocheck ForceAnnealing_with_BondGraph ForceAnnealing_with_BondGraph_continued ForceAnnealing_with_BondGraph_continued_betteresults ForceAnnealing_with_BondGraph_contraction-expansion FragmentAction_writes_AtomFragments FragmentMolecule_checks_bonddegrees GeometryObjects Gui_Fixes Gui_displays_atomic_force_velocity ImplicitCharges IndependentFragmentGrids IndependentFragmentGrids_IndividualZeroInstances IndependentFragmentGrids_IntegrationTest IndependentFragmentGrids_Sole_NN_Calculation JobMarket_RobustOnKillsSegFaults JobMarket_StableWorkerPool JobMarket_unresolvable_hostname_fix MoreRobust_FragmentAutomation ODR_violation_mpqc_open PartialCharges_OrthogonalSummation PdbParser_setsAtomName PythonUI_with_named_parameters QtGui_reactivate_TimeChanged_changes Recreated_GuiChecks Rewrite_FitPartialCharges RotateToPrincipalAxisSystem_UndoRedo SaturateAtoms_findBestMatching SaturateAtoms_singleDegree StoppableMakroAction Subpackage_CodePatterns Subpackage_JobMarket Subpackage_LinearAlgebra Subpackage_levmar Subpackage_mpqc_open Subpackage_vmg Switchable_LogView ThirdParty_MPQC_rebuilt_buildsystem TrajectoryDependenant_MaxOrder TremoloParser_IncreasedPrecision TremoloParser_MultipleTimesteps TremoloParser_setsAtomName Ubuntu_1604_changes stable
Last change on this file since ed2551 was ed2551, checked in by Frederik Heber <heber@…>, 12 years ago

All Pair_.. and ManyBodyPotentials are now also derived from SerializablePotential.

  • added static ParameterNames, skipping "energy_offset" which should not go to file.
  • overrode operator<<() and ..>>() for SaturationPotential to print both potentials independently.
  • adapted LevMartest and all unit tests.
  • Property mode set to 100644
File size: 15.9 KB
Line 
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 * LevMartester.cpp
26 *
27 * Created on: Sep 27, 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 <boost/archive/text_iarchive.hpp>
38
39#include "CodePatterns/MemDebug.hpp"
40
41#include <boost/assign.hpp>
42#include <boost/assign/list_of.hpp>
43#include <boost/bind.hpp>
44#include <boost/filesystem.hpp>
45#include <boost/function.hpp>
46#include <boost/program_options.hpp>
47
48#include <cstdlib>
49#include <ctime>
50#include <fstream>
51#include <iostream>
52#include <iterator>
53#include <list>
54#include <vector>
55
56#include <levmar.h>
57
58#include "CodePatterns/Assert.hpp"
59#include "CodePatterns/Log.hpp"
60
61#include "LinearAlgebra/Vector.hpp"
62
63#include "Fragmentation/Homology/HomologyContainer.hpp"
64#include "Fragmentation/SetValues/Fragment.hpp"
65#include "FunctionApproximation/Extractors.hpp"
66#include "FunctionApproximation/FunctionApproximation.hpp"
67#include "FunctionApproximation/FunctionModel.hpp"
68#include "FunctionApproximation/TrainingData.hpp"
69#include "FunctionApproximation/writeDistanceEnergyTable.hpp"
70#include "Helpers/defs.hpp"
71#include "Potentials/Specifics/PairPotential_Morse.hpp"
72#include "Potentials/Specifics/PairPotential_Angle.hpp"
73#include "Potentials/Specifics/SaturationPotential.hpp"
74#include "types.hpp"
75
76namespace po = boost::program_options;
77
78using namespace boost::assign;
79
80HomologyGraph getFirstGraphwithTimesSpecificElement(
81 const HomologyContainer &homologies,
82 const size_t _number,
83 const size_t _times)
84{
85 for (HomologyContainer::container_t::const_iterator iter =
86 homologies.begin(); iter != homologies.end(); ++iter) {
87 if (iter->first.hasTimesAtomicNumber(_number,_times))
88 return iter->first;
89 }
90 return HomologyGraph();
91}
92
93/** This function returns the elements of the sum over index "k" for an
94 * argument containing indices "i" and "j"
95 * @param inputs vector of all configuration (containing each a vector of all arguments)
96 * @param arg argument containing indices "i" and "j"
97 * @param cutoff cutoff criterion for sum over k
98 * @return vector of argument pairs (a vector) of ik and jk for at least all k
99 * within distance of \a cutoff to i
100 */
101std::vector<FunctionModel::arguments_t>
102getTripleFromArgument(const FunctionApproximation::inputs_t &inputs, const argument_t &arg, const double cutoff)
103{
104 typedef std::list<argument_t> arg_list_t;
105 typedef std::map<size_t, arg_list_t > k_args_map_t;
106 k_args_map_t tempresult;
107 ASSERT( inputs.size() > arg.globalid,
108 "getTripleFromArgument() - globalid "+toString(arg.globalid)
109 +" is greater than all inputs "+toString(inputs.size())+".");
110 const FunctionModel::arguments_t &listofargs = inputs[arg.globalid];
111 for (FunctionModel::arguments_t::const_iterator argiter = listofargs.begin();
112 argiter != listofargs.end();
113 ++argiter) {
114 // first index must be either i or j but second index not
115 if (((argiter->indices.first == arg.indices.first)
116 || (argiter->indices.first == arg.indices.second))
117 && ((argiter->indices.second != arg.indices.first)
118 && (argiter->indices.second != arg.indices.second))) {
119 // we need arguments ik and jk
120 std::pair< k_args_map_t::iterator, bool> inserter =
121 tempresult.insert( std::make_pair( argiter->indices.second, arg_list_t(1,*argiter)));
122 if (!inserter.second) {
123 // is present one ik or jk, if ik insert jk at back
124 if (inserter.first->second.begin()->indices.first == arg.indices.first)
125 inserter.first->second.push_back(*argiter);
126 else // if jk, insert ik at front
127 inserter.first->second.push_front(*argiter);
128 }
129 }
130// // or second index must be either i or j but first index not
131// else if (((argiter->indices.first != arg.indices.first)
132// && (argiter->indices.first != arg.indices.second))
133// && ((argiter->indices.second == arg.indices.first)
134// || (argiter->indices.second == arg.indices.second))) {
135// // we need arguments ki and kj
136// std::pair< k_args_map_t::iterator, bool> inserter =
137// tempresult.insert( std::make_pair( argiter->indices.first, arg_list_t(1,*argiter)));
138// if (!inserter.second) {
139// // is present one ki or kj, if ki insert kj at back
140// if (inserter.first->second.begin()->indices.second == arg.indices.first)
141// inserter.first->second.push_back(*argiter);
142// else // if kj, insert ki at front
143// inserter.first->second.push_front(*argiter);
144// }
145// }
146 }
147 // check that i,j are NOT contained
148 ASSERT( tempresult.count(arg.indices.first) == 0,
149 "getTripleFromArgument() - first index of argument present in k_args_map?");
150 ASSERT( tempresult.count(arg.indices.second) == 0,
151 "getTripleFromArgument() - first index of argument present in k_args_map?");
152
153 // convert
154 std::vector<FunctionModel::arguments_t> result;
155 for (k_args_map_t::const_iterator iter = tempresult.begin();
156 iter != tempresult.end();
157 ++iter) {
158 ASSERT( iter->second.size() == 2,
159 "getTripleFromArgument() - for index "+toString(iter->first)+" we did not find both ik and jk.");
160 result.push_back( FunctionModel::arguments_t(iter->second.begin(), iter->second.end()) );
161 }
162 return result;
163}
164
165int main(int argc, char **argv)
166{
167 std::cout << "Hello to the World from LevMar!" << std::endl;
168
169 // load homology file
170 po::options_description desc("Allowed options");
171 desc.add_options()
172 ("help", "produce help message")
173 ("homology-file", po::value< boost::filesystem::path >(), "homology file to parse")
174 ;
175
176 po::variables_map vm;
177 po::store(po::parse_command_line(argc, argv, desc), vm);
178 po::notify(vm);
179
180 if (vm.count("help")) {
181 std::cout << desc << "\n";
182 return 1;
183 }
184
185 boost::filesystem::path homology_file;
186 if (vm.count("homology-file")) {
187 homology_file = vm["homology-file"].as<boost::filesystem::path>();
188 LOG(1, "INFO: Parsing " << homology_file.string() << ".");
189 } else {
190 LOG(0, "homology-file level was not set.");
191 }
192 HomologyContainer homologies;
193 if (boost::filesystem::exists(homology_file)) {
194 std::ifstream returnstream(homology_file.string().c_str());
195 if (returnstream.good()) {
196 boost::archive::text_iarchive ia(returnstream);
197 ia >> homologies;
198 } else {
199 ELOG(2, "Failed to parse from " << homology_file.string() << ".");
200 }
201 returnstream.close();
202 } else {
203 ELOG(0, homology_file << " does not exist.");
204 }
205
206 // first we try to look into the HomologyContainer
207 LOG(1, "INFO: Listing all present homologies ...");
208 for (HomologyContainer::container_t::const_iterator iter =
209 homologies.begin(); iter != homologies.end(); ++iter) {
210 LOG(1, "INFO: graph " << iter->first << " has Fragment "
211 << iter->second.first << " and associated energy " << iter->second.second << ".");
212 }
213
214 /******************** Angle TRAINING ********************/
215 {
216 // then we ought to pick the right HomologyGraph ...
217 const HomologyGraph graph = getFirstGraphwithTimesSpecificElement(homologies,8,1);
218 if (graph != HomologyGraph()) {
219 LOG(1, "First representative graph containing three saturated carbons is " << graph << ".");
220
221 // Afterwards we go through all of this type and gather the distance and the energy value
222 Fragment::charges_t h2o;
223 h2o += 8,1,1;
224 TrainingData AngleData(
225 boost::bind(&Extractors::gatherDistancesFromFragment,
226 boost::bind(&Fragment::getPositions, _1),
227 boost::bind(&Fragment::getCharges, _1),
228 boost::cref(h2o),
229 _2)
230 );
231 AngleData(homologies.getHomologousGraphs(graph));
232 LOG(1, "INFO: I gathered the following training data: " << AngleData);
233 // NOTICE that distance are in bohrradi as they come from MPQC!
234
235 // now perform the function approximation by optimizing the model function
236 FunctionModel::parameters_t params(PairPotential_Angle::MAXPARAMS, 0.);
237 params[PairPotential_Angle::energy_offset] = -1.;
238 params[PairPotential_Angle::spring_constant] = 1.;
239 params[PairPotential_Angle::equilibrium_distance] = 0.2;
240 PairPotential_Angle::ParticleTypes_t types =
241 boost::assign::list_of<PairPotential_Angle::ParticleType_t>
242 (8)(1)(1)
243 ;
244 PairPotential_Angle angle(types);
245 LOG(0, "INFO: Initial parameters are " << params << ".");
246 angle.setParameters(params);
247 FunctionModel &model = angle;
248 FunctionApproximation approximator(AngleData, model);
249 if (model.isBoxConstraint() && approximator.checkParameterDerivatives())
250 approximator(FunctionApproximation::ParameterDerivative);
251 else {
252 ELOG(0, "We require parameter derivatives for a box constraint minimization.");
253 return 1;
254 }
255 params = model.getParameters();
256
257 LOG(0, "RESULT: Best parameters are " << params << ".");
258 }
259 }
260
261 /******************** MORSE TRAINING ********************/
262 {
263 // then we ought to pick the right HomologyGraph ...
264 const HomologyGraph graph = getFirstGraphwithTimesSpecificElement(homologies,8,1);
265 if (graph != HomologyGraph()) {
266 LOG(1, "First representative graph containing two saturated carbons is " << graph << ".");
267
268 // Afterwards we go through all of this type and gather the distance and the energy value
269 Fragment::charges_t h2o;
270 h2o += 8,1;
271 TrainingData MorseData(
272 boost::bind(&Extractors::gatherDistancesFromFragment,
273 boost::bind(&Fragment::getPositions, _1),
274 boost::bind(&Fragment::getCharges, _1),
275 boost::cref(h2o),
276 _2)
277 );
278 MorseData(homologies.getHomologousGraphs(graph));
279 LOG(1, "INFO: I gathered the following training data: " << MorseData);
280 // NOTICE that distance are in bohrradi as they come from MPQC!
281
282 // now perform the function approximation by optimizing the model function
283 FunctionModel::parameters_t params(PairPotential_Morse::MAXPARAMS, 0.);
284 params[PairPotential_Morse::dissociation_energy] = 0.5;
285 params[PairPotential_Morse::energy_offset] = -1.;
286 params[PairPotential_Morse::spring_constant] = 1.;
287 params[PairPotential_Morse::equilibrium_distance] = 2.9;
288 PairPotential_Morse::ParticleTypes_t types =
289 boost::assign::list_of<PairPotential_Morse::ParticleType_t>
290 (8)(1)
291 ;
292 PairPotential_Morse morse(types);
293 LOG(0, "INFO: Initial parameters are " << params << ".");
294 morse.setParameters(params);
295 FunctionModel &model = morse;
296 FunctionApproximation approximator(MorseData, model); // we only give CC distance, hence 1 input dim
297 if (model.isBoxConstraint() && approximator.checkParameterDerivatives())
298 approximator(FunctionApproximation::ParameterDerivative);
299 else {
300 ELOG(0, "We require parameter derivatives for a box constraint minimization.");
301 return 1;
302 }
303 params = model.getParameters();
304
305 LOG(0, "RESULT: Best parameters are " << params << ".");
306 }
307 }
308
309 /******************* SATURATION TRAINING *******************/
310 FunctionModel::parameters_t params(SaturationPotential::MAXPARAMS, 0.);
311 {
312 // then we ought to pick the right HomologyGraph ...
313 const HomologyGraph graph = getFirstGraphwithTimesSpecificElement(homologies,8,1);
314 if (graph != HomologyGraph()) {
315 LOG(1, "First representative graph containing one saturated carbon is " << graph << ".");
316
317 // Afterwards we go through all of this type and gather the distance and the energy value
318 Fragment::charges_t h2o;
319 h2o += 8,1,1;
320 TrainingData TersoffData(
321 boost::bind(&Extractors::gatherDistancesFromFragment,
322 boost::bind(&Fragment::getPositions, _1),
323 boost::bind(&Fragment::getCharges, _1),
324 boost::cref(h2o),
325 _2)
326 );
327// TrainingData::extractor_t(&Extractors::gatherAllDistances) // gather first carbon pair
328// );
329 TersoffData( homologies.getHomologousGraphs(graph) );
330 LOG(1, "INFO: I gathered the following training data: " << TersoffData);
331 // NOTICE that distance are in bohrradi as they come from MPQC!
332
333 // now perform the function approximation by optimizing the model function
334 boost::function< std::vector<FunctionModel::arguments_t>(const argument_t &, const double)> triplefunction =
335 boost::bind(&getTripleFromArgument, boost::cref(TersoffData.getTrainingInputs()), _1, _2);
336 srand((unsigned)time(0)); // seed with current time
337 params[SaturationPotential::all_energy_offset] =
338 TersoffData.getTrainingOutputAverage()[0]; //1e+1*rand()/(double)RAND_MAX;//3.487900e+00;
339 params[SaturationPotential::morse_spring_constant] = 1e+1*rand()/(double)RAND_MAX;//1.393600e+03;
340 params[SaturationPotential::morse_equilibrium_distance] = 1e+1*rand()/(double)RAND_MAX;//3.467000e+02;
341 params[SaturationPotential::morse_dissociation_energy] = 1e+1*rand()/(double)RAND_MAX;//3.487900e+00;
342 params[SaturationPotential::angle_spring_constant] = 1e+1*rand()/(double)RAND_MAX;//3.487900e+00;
343 params[SaturationPotential::angle_equilibrium_distance] = 2e+0*rand()/(double)RAND_MAX - 1.;//3.487900e+00;
344 LOG(0, "INFO: Initial parameters are " << params << ".");
345 SaturationPotential::ParticleTypes_t types =
346 boost::assign::list_of<SaturationPotential::ParticleType_t>
347 (8)(1)(1)
348 ;
349
350 SaturationPotential saturation(types, 2.5, triplefunction);
351 saturation.setParameters(params);
352 FunctionModel &model = saturation;
353 FunctionApproximation approximator(TersoffData, model); // CH4 has 5 atoms, hence 5*4/2 distances
354 if (model.isBoxConstraint() && approximator.checkParameterDerivatives())
355 approximator(FunctionApproximation::ParameterDerivative);
356 else {
357 ELOG(0, "We require parameter derivatives for a box constraint minimization.");
358 return 1;
359 }
360
361 params = model.getParameters();
362 LOG(0, "RESULT: Best parameters are " << params << ".");
363
364 // std::cout << "\tsaturationparticle:";
365 // std::cout << "\tparticle_type=C,";
366 // std::cout << "\tA=" << params[SaturationPotential::A] << ",";
367 // std::cout << "\tB=" << params[SaturationPotential::B] << ",";
368 // std::cout << "\tlambda=" << params[SaturationPotential::lambda] << ",";
369 // std::cout << "\tmu=" << params[SaturationPotential::mu] << ",";
370 // std::cout << "\tbeta=" << params[SaturationPotential::beta] << ",";
371 // std::cout << "\tn=" << params[SaturationPotential::n] << ",";
372 // std::cout << "\tc=" << params[SaturationPotential::c] << ",";
373 // std::cout << "\td=" << params[SaturationPotential::d] << ",";
374 // std::cout << "\th=" << params[SaturationPotential::h] << ",";
375 //// std::cout << "\toffset=" << params[SaturationPotential::offset] << ",";
376 // std::cout << "\tR=" << saturation.R << ",";
377 // std::cout << "\tS=" << saturation.S << ";";
378 // std::cout << std::endl;
379
380 // check L2 and Lmax error against training set
381 LOG(1, "INFO: L2sum = " << TersoffData.getL2Error(model)
382 << ", LMax = " << TersoffData.getLMaxError(model) << ".");
383 }
384
385 }
386
387 return 0;
388}
Note: See TracBrowser for help on using the repository browser.