Ignore:
Timestamp:
Feb 27, 2013, 12:43:29 PM (12 years ago)
Author:
Frederik Heber <heber@…>
Branches:
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
Children:
8c9049
Parents:
64bdfd
git-author:
Frederik Heber <heber@…> (12/19/12 15:32:48)
git-committer:
Frederik Heber <heber@…> (02/27/13 12:43:29)
Message:

Modified SaturationPotential to allow fitting to e.g. HCC as well.

  • i.e. now not only saturated elements are fitted but also "saturated" C-C bonds.
File:
1 edited

Legend:

Unmodified
Added
Removed

  • src/Potentials/Specifics/SaturationPotential.cpp

    r64bdfd r5c3ad5  
    7474  SerializablePotential(_ParticleTypes),
    7575  morse(_ParticleTypes),
    76   angle(symmetrizeTypes(_ParticleTypes)),
     76  angle(addSaturationType(_ParticleTypes)),
    7777  energy_offset(0.)
    7878{
    … …  
    8181  ASSERT( _ParticleTypes.size() == (size_t)2,
    8282      "SaturationPotential::SaturationPotential() - exactly two types must be given.");
    83   ASSERT( _ParticleTypes[1] == 1,
    84       "SaturationPotential::SaturationPotential() - second type must be hydrogen.");
     83//  ASSERT( _ParticleTypes[1] == 1,
     84//      "SaturationPotential::SaturationPotential() - second type must be hydrogen.");
    8585}
    8686
    … …  
    9595  SerializablePotential(_ParticleTypes),
    9696  morse(_ParticleTypes),
    97   angle(symmetrizeTypes(_ParticleTypes)),
     97  angle(addSaturationType(_ParticleTypes)),
    9898  energy_offset(_all_energy_offset)
    9999{
    100100  ASSERT( _ParticleTypes.size() == (size_t)2,
    101101      "SaturationPotential::SaturationPotential() - exactly two types must be given.");
    102   ASSERT( _ParticleTypes[1] == 1,
    103       "SaturationPotential::SaturationPotential() - second type must be hydrogen.");
     102//  ASSERT( _ParticleTypes[1] == 1,
     103//      "SaturationPotential::SaturationPotential() - second type must be hydrogen.");
    104104  parameters_t morse_params(morse.getParameterDimension());
    105105  morse_params[PairPotential_Morse::spring_constant] = _morse_spring_constant;
    … …  
    269269  double result = 0.;
    270270  const ParticleTypes_t &morse_types = morse.getParticleTypes();
     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;
    271279  for(arguments_t::const_iterator argiter = arguments.begin();
    272280      argiter != arguments.end();
    … …  
    283291      if (result != result)
    284292        ELOG(1, "result is NAN.");
    285 
     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]))) {
    286296      // Angle contribution
    287297      {
    … …  
    295305              Extractors::reorderArgumentsByParticleTypes(*iter, angle.getParticleTypes());
    296306          // We get both angles, e.g. 0-4-1 and 1-4-0, hence multiply with 0.5
    297           const double tmp = .5*angle(tempargs)[0];  // as we have all distances we get both jk and kj
     307          const double tmp = multiplicity*angle(tempargs)[0];  // as we have all distances we get both jk and kj
    298308//          LOG(3, "DEBUG: angle yields " << tmp << " for << " << tempargs << ".");
    299309          result += tmp;
    … …  
    325335{
    326336  double result = 0.;
    327   if (index == all_energy_offset) {
     337  switch (index) {
     338  case all_energy_offset:
    328339    result = 1.;
    329   } else {
    330     const ParticleTypes_t &morse_types = morse.getParticleTypes();
    331     for(arguments_t::const_iterator argiter = arguments.begin();
    332         argiter != arguments.end();
    333         ++argiter) {
    334       const argument_t &r_ij = *argiter;
    335       if (((r_ij.types.first == morse_types[0]) && (r_ij.types.second == morse_types[1]))
    336           || ((r_ij.types.first == morse_types[1]) && (r_ij.types.second == morse_types[0]))) {
    337         arguments_t args(1, r_ij);
    338 
    339         switch (index) {
    340           case morse_spring_constant:
    341             result += morse.parameter_derivative(args, PairPotential_Morse::spring_constant)[0];
    342             break;
    343           case morse_equilibrium_distance:
    344             result += morse.parameter_derivative(args, PairPotential_Morse::equilibrium_distance)[0];
    345             break;
    346           case morse_dissociation_energy:
    347             result += morse.parameter_derivative(args, PairPotential_Morse::dissociation_energy)[0];
    348             break;
    349           case angle_spring_constant:
    350           case angle_equilibrium_distance:
    351           {
    352             typedef std::vector<FunctionModel::arguments_t> tripleargs_t;
    353             tripleargs_t tripleargs =
    354                 triplefunction(r_ij, arguments);
    355             for (tripleargs_t::const_iterator iter = tripleargs.begin();
    356                 iter != tripleargs.end();
    357                 ++iter) {
    358               FunctionModel::arguments_t tempargs =
    359                   Extractors::reorderArgumentsByParticleTypes(*iter, angle.getParticleTypes());
    360               // We get both angles, e.g. 0-4-1 and 1-4-0, hence multiply with 0.5
    361               if (index == angle_spring_constant)
    362                 result += .5*angle.parameter_derivative(tempargs, PairPotential_Angle::spring_constant)[0];
    363               else if (index == angle_equilibrium_distance)
    364                 result += .5*angle.parameter_derivative(tempargs, PairPotential_Angle::equilibrium_distance)[0];
    365 //              LOG(2, "DEBUG: angle yields " << tmp << " for << " << arguments << ".");
    366 //              result += tmp;
    367               if (result != result)
    368                 ELOG(1, "result is NAN.");
    369             }
     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;
    370368          }
    371             break;
    372           default:
    373             ASSERT( 0, "SaturationPotential::parameter_derivative() - impossible to get here.");
    374             break;
     369//          LOG(2, "DEBUG: morse yields " << tmp << " for << " << args << ".");
     370          result += tmp;
    375371        }
    376372      }
    377373    }
     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.");
     410          }
     411        }
     412      }
     413    }
     414    break;
     415  default:
     416    ASSERT( 0, "SaturationPotential::parameter_derivative() - impossible to get here.");
     417    break;
    378418  }
    379419  return SaturationPotential::results_t(1, result);
    … …  
    396436}
    397437
     438const SaturationPotential::ParticleTypes_t
     439SaturationPotential::addSaturationType(const ParticleTypes_t &_ParticleTypes)
     440{
     441  ParticleTypes_t types(_ParticleTypes);
     442  types.push_back( ParticleType_t(1) );
     443  return types;
     444}
     445
    398446FunctionModel::extractor_t
    399447SaturationPotential::getFragmentSpecificExtractor() const
    400448{
    401   Fragment::charges_t charges;
    402   charges.resize(getParticleTypes().size());
    403   std::transform(getParticleTypes().begin(), getParticleTypes().end(),
    404       charges.begin(), boost::lambda::_1);       
     449//  Fragment::charges_t charges;
     450//  charges.resize(getParticleTypes().size());
     451//  std::transform(getParticleTypes().begin(), getParticleTypes().end(),
     452//      charges.begin(), boost::lambda::_1);     
    405453  FunctionModel::extractor_t returnfunction;
    406   if (charges[0] == charges[1]) {
    407     // In case both types are equal there is only a single pair of possible
    408     // type combinations.
    409      returnfunction =
    410         boost::bind(&Extractors::gatherAllDistancesFromFragment,
    411             boost::bind(&Fragment::getPositions, _1),
    412             boost::bind(&Fragment::getCharges, _1),
    413             charges, // is only temporarily created, hence copy
    414             _2);
    415   } else {
     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 {
    416464    // we have to chain here a rather complex "tree" of functions
    417465    // as we only have a couple of ParticleTypes but need to get
    … …  
    419467    // Finally, we also need to arrange them in correct order
    420468    // (for PairPotentiale_Angle).
    421     //charges_t firstpair(2, boost::cref(charges[0]));
     469//    charges_t firstpair(2, boost::cref(charges[0]));
    422470    // only that saturation potential never has its middle element twice!
    423471    // hence, we skip the firstpair but keep the code for later generalization
    424     Fragment::charges_t secondpair(2, boost::cref(charges[1]));
    425     const Fragment::charges_t &thirdpair = charges;
     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    }
    426486    returnfunction =
    427487//        boost::bind(&Extractors::reorderArgumentsByParticleTypes,
    428488          boost::bind(&Extractors::combineArguments,
    429489//            boost::bind(&Extractors::combineArguments,
    430 //              boost::bind(&Extractors::gatherAllDistancesFromFragment,
    431 //                  boost::bind(&Fragment::getPositions, _1),
    432 //                  boost::bind(&Fragment::getCharges, _1),
    433 //                  firstpair,  // no crefs here as are temporaries!
    434 //                  _2),
    435490              boost::bind(&Extractors::gatherAllDistancesFromFragment,
    436491                  boost::bind(&Fragment::getPositions, _1),
    437492                  boost::bind(&Fragment::getCharges, _1),
    438                   secondpair,  // no crefs here as are temporaries!
     493                  charges_angle,  // no crefs here as are temporaries!
    439494                  _2),
    440 //            ),
    441             boost::bind(&Extractors::gatherAllDistancesFromFragment,
    442                 boost::bind(&Fragment::getPositions, _1),
    443                 boost::bind(&Fragment::getCharges, _1),
    444                 boost::cref(thirdpair), // only the last one is no temporary
    445                 _2)
    446           );
    447 //          boost::cref(angle.getParticleTypes())
    448 //        );
    449 }
     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//  }
    450510  return returnfunction;
    451511}
Note: See TracChangeset for help on using the changeset viewer.