Changeset bcb593


Ignore:
Timestamp:
Jul 11, 2013, 9:43:31 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:
9c1324
Parents:
11f0fa
git-author:
Frederik Heber <heber@…> (06/18/13 14:35:36)
git-committer:
Frederik Heber <heber@…> (07/11/13 21:43:31)
Message:

FIX: VerletIntegrationAction now assumes forces have just been calculated.

  • according to Wikipedia's Velocity_Verlet info, we first integrate position, then calculate forces, then integrate velocity. This assumes that forces based on next time step's position are already known. This is possible when parsed from file but not if they are calculated dynamically via fragmentation. Hence, we now integrate velocity from last time step to current, then integrate position from current time step to next. Then we are in the position to calculate forces and do this cycle again.
  • for this, VelocityVerletUpdate was split into ..X() and ..U() for independent integration of position and velocity.
  • VelocityVerletIntegration::operator() now first corrects forces, then integrates velocites, corrects them too and finally integrates positions according to the above new scheme.
  • removed option MDSteps from VerletIntegrationAction.
  • TESTFIX: Had to change regression rest on VerletIntegration since the cycle sequence has changed. It was nonsense before to have the forces already in some file respective to future position that actually first need to come out of the time integration.
Files:
8 edited

Legend:

Unmodified
Added
Removed

  • src/Actions/MoleculeAction/VerletIntegrationAction.cpp

    r11f0fa rbcb593  
    4242#include "molecule.hpp"
    4343#include "World.hpp"
     44#include "WorldTime.hpp"
    4445
    4546#include <vector>
    … …  
    5960  // TODO: sollte besser stream nutzen, nicht filename direkt (es sei denn, ist prefix), besser fuer unit test
    6061  AtomSetMixin<std::vector<atom *> > set(World::getInstance().getSelectedAtoms());
    61   for (unsigned int step = 0; step < params.MDSteps.get(); ++step) {
    62     VerletForceIntegration<std::vector<atom *> > Verlet(set, params.Deltat.get(), false);
    63     // parse forces into next step
    64     if (!params.forcesfile.get().string().empty()) {
    65       LOG(1, "Parsing forces file.");
    66       if (!Verlet.parseForcesFile(params.forcesfile.get().string().c_str(), step))
    67         LOG(2, "File " << params.forcesfile.get() << " not found.");
    68       else
    69         LOG(2, "File " << params.forcesfile.get() << " found and parsed.");
    70     }
    71     // perform velocity verlet update
    72     LOG(1, "Verlet integrating.");
    73     if (!Verlet(step+1, 1, 0, params.FixedCenterOfMass.get()))
     62  // we always operate relative to current time step
     63  size_t CurrentStep = WorldTime::getInstance().getTime();
     64  VerletForceIntegration<std::vector<atom *> > Verlet(set, params.Deltat.get(), true);
     65  // parse forces into next step
     66  if (!params.forcesfile.get().string().empty()) {
     67    LOG(1, "Parsing forces file.");
     68    if (!Verlet.parseForcesFile(params.forcesfile.get().string().c_str(), CurrentStep))
    7469      LOG(2, "File " << params.forcesfile.get() << " not found.");
    7570    else
    7671      LOG(2, "File " << params.forcesfile.get() << " found and parsed.");
    7772  }
     73  // perform velocity verlet update
     74  LOG(1, "Verlet integrating.");
     75  if (!Verlet(CurrentStep+1, 1, 0, params.FixedCenterOfMass.get()))
     76    ELOG(1, "Velocity verlet failed.");
     77  else
     78    LOG(2, "DEBUG: Successfully performed updates on velocity and position.");
     79  // increment to next time step
     80  World::getInstance().setTime(CurrentStep+1);
    7881
    7982  return Action::success;

  • src/Actions/MoleculeAction/VerletIntegrationAction.def

    r11f0fa rbcb593  
    1717// ValueStorage by the token "Z" -> first column: int, Z, "Z"
    1818// "undefine" if no parameters are required, use (NOPARAM_DEFAULT) for each (undefined) default value
    19 #define paramtypes (boost::filesystem::path)(double)(unsigned int)(bool)
    20 #define paramtokens ("forces-file")("deltat")("MDSteps")("keep-fixed-CenterOfMass")
    21 #define paramdescriptions ("file containing")("time step width")("number of MDSteps to integrate")("whether forces and velocities shall be corrected such that center of mass remains at rest")
    22 #define paramdefaults (PARAM_DEFAULT(""))(NOPARAM_DEFAULT)(NOPARAM_DEFAULT)(NOPARAM_DEFAULT)
    23 #define paramreferences (forcesfile)(Deltat)(MDSteps)(FixedCenterOfMass)
     19#define paramtypes (boost::filesystem::path)(double)(bool)
     20#define paramtokens ("forces-file")("deltat")("keep-fixed-CenterOfMass")
     21#define paramdescriptions ("file containing")("time step width")("whether forces and velocities shall be corrected such that center of mass remains at rest")
     22#define paramdefaults (PARAM_DEFAULT(""))(NOPARAM_DEFAULT)(NOPARAM_DEFAULT)
     23#define paramreferences (forcesfile)(Deltat)(FixedCenterOfMass)
    2424#define paramvalids \
    2525(DummyValidator< boost::filesystem::path >()) \
    2626(PositiveValidator< double >()) \
    27 (NotZeroValidator< unsigned int >()) \
    2827(DummyValidator< bool >())
    2928

  • src/Atom/atom_atominfo.cpp

    r11f0fa rbcb593  
    279279{
    280280  OBSERVE;
    281   NOTIFY(AtomObservable::VelocityChanged);
     281  NOTIFY(AtomObservable::ForceChanged);
    282282  ASSERT(WorldTime::getTime() < AtomicForce.size(),
    283283      "AtomInfo::setAtomicForce() - Access out of range: "
    … …  
    291291  OBSERVE;
    292292  if (WorldTime::getTime() == _step)
    293     NOTIFY(AtomObservable::VelocityChanged);
     293    NOTIFY(AtomObservable::ForceChanged);
    294294  const unsigned int size = AtomicForce.size();
    295295  ASSERT(_step <= size,
    … …  
    579579};
    580580
    581 /** Performs a velocity verlet update of the trajectory.
    582  * Parameters are according to those in configuration class.
     581/** Performs a velocity verlet update of the position at \a NextStep from \a LastStep information only.
     582 *
     583 * We calculate \f$x(t + \delta t) = x(t) + v(t)* \delta t + .5 * \delta t * \delta t * F(t)/m \f$.
     584 *
     585 *
    583586 * \param NextStep index of sequential step to set
    584587 * \param Deltat time step width
    585588 * \param IsAngstroem whether the force's underlying unit of length is angstroem or bohr radii
    586589 */
    587 void AtomInfo::VelocityVerletUpdate(int nr, const unsigned int NextStep, double Deltat, bool IsAngstroem)
     590void AtomInfo::VelocityVerletUpdateX(int nr, const unsigned int NextStep, double Deltat, bool IsAngstroem)
    588591{
    589592  const unsigned int LastStep = NextStep == 0 ? 0 : NextStep-1;
    590593
    591594  LOG(2, "INFO: Particle that currently " << *this);
    592   LOG(2, "INFO: Integrating with mass=" << getMass() << " and Deltat="
     595  LOG(2, "INFO: Integrating position with mass=" << getMass() << " and Deltat="
    593596      << Deltat << " at NextStep=" << NextStep);
    594   // update force
    595   // (F+F_old)/2m = a and thus: v = (F+F_old)/2m * t = (F + F_old) * a
    596   LOG(3, "INFO: Force at step " << NextStep << " is " << getAtomicForceAtStep(NextStep));
    597597
    598598  // update position
    … …  
    603603    LOG(4, "INFO: position with velocity " << getAtomicVelocityAtStep(LastStep) << " from last step " << tempVector);
    604604    tempVector += .5*Deltat*Deltat*(getAtomicForceAtStep(LastStep))*(1./getMass());     // F = m * a and s =
    605     LOG(4, "INFO: position with force " << getAtomicForceAtStep(NextStep) << " from last step " << tempVector);
     605    LOG(4, "INFO: position with force " << getAtomicForceAtStep(LastStep) << " from last step " << tempVector);
    606606    setPositionAtStep(NextStep, tempVector);
    607607    LOG(3, "INFO: Position at step " << NextStep << " set to " << tempVector);
    608608  }
     609};
     610
     611/** Performs a velocity verlet update of the velocity at \a NextStep.
     612 *
     613 * \note forces at NextStep should have been calculated based on position at NextStep prior
     614 * to calling this function.
     615 *
     616 * We calculate \f$v(t) = v(t - \delta t) + \delta _t * .5 * (F(t - \delta t) + F(t))/m \f$.
     617 *
     618 * Parameters are according to those in configuration class.
     619 * \param NextStep index of sequential step to set
     620 * \param Deltat time step width
     621 * \param IsAngstroem whether the force's underlying unit of length is angstroem or bohr radii
     622 */
     623void AtomInfo::VelocityVerletUpdateU(int nr, const unsigned int NextStep, double Deltat, bool IsAngstroem)
     624{
     625  const unsigned int LastStep = NextStep == 0 ? 0 : NextStep-1;
     626
     627  LOG(2, "INFO: Particle that currently " << *this);
     628  LOG(2, "INFO: Integrating velocity with mass=" << getMass() << " and Deltat="
     629      << Deltat << " at NextStep=" << NextStep);
    609630
    610631  // Update U

  • src/Atom/atom_atominfo.hpp

    r11f0fa rbcb593  
    267267  size_t getTrajectorySize() const;
    268268  void CopyStepOnStep(const unsigned int dest, const unsigned int src);
    269   void VelocityVerletUpdate(int nr, const unsigned int NextStep, double Deltat, bool IsAngstroem);
     269  void VelocityVerletUpdateX(int nr, const unsigned int NextStep, double Deltat, bool IsAngstroem);
     270  void VelocityVerletUpdateU(int nr, const unsigned int NextStep, double Deltat, bool IsAngstroem);
    270271  double getKineticEnergy(const unsigned int step) const;
    271272  Vector getMomentum(const unsigned int step) const;

  • src/Dynamics/VerletForceIntegration.hpp

    r11f0fa rbcb593  
    7878
    7979  /** Performs Verlet integration.
    80    * Note that we assume the parsed forces to be in atomic units (hence, if coordinates are in angstroem, we
    81    * have to transform them).
    82    * This adds a new MD step \f$ t + \Delta t \f$ to the config file.
     80   *
     81   * We assume that forces have just been calculated. Then, we perform the velocity
     82   * and the position calculation for \f$ t + \Delta t \f$, such that forces may be
     83   * again calculated with respect to the new position.
     84   *
    8385   * \param NextStep current time step (i.e. \f$ t + \Delta t \f$ in the sense of the velocity verlet)
    8486   * \param offset offset in matrix file to the first force component
    … …  
    100102    // make sum of forces equal zero
    101103    if (FixedCenterOfMass)
    102       correctForceMatrixForFixedCenterOfMass(offset,NextStep);
     104      correctForceMatrixForFixedCenterOfMass(offset,NextStep-1);
    103105
    104106    // solve a constrained potential if we are meant to
    105     if (DoConstrainedMD) {
    106       performConstraintMinimization(DoConstrainedMD,NextStep);
     107    if (DoConstrainedMD)
     108      performConstraintMinimization(DoConstrainedMD,NextStep-1);
     109
     110    if (NextStep > 0) {
     111      for(typename AtomSetMixin<T>::iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
     112        //std::cout << "Id of atom is " << (*iter)->getId() << std::endl;
     113        (*iter)->VelocityVerletUpdateU((*iter)->getId(), NextStep-1, Deltat, IsAngstroem);
     114      }
     115
     116      // make sum of velocities equal zero
     117      if (FixedCenterOfMass)
     118        correctVelocitiesForFixedCenterOfMass(NextStep-1);
     119
     120      // thermostat
     121      performThermostatControl(NextStep-1);
    107122    }
    108123
    … …  
    112127    for(typename AtomSetMixin<T>::iterator iter = atoms.begin(); iter != atoms.end(); ++iter) {
    113128      //std::cout << "Id of atom is " << (*iter)->getId() << std::endl;
    114       (*iter)->VelocityVerletUpdate((*iter)->getId(), NextStep, Deltat, IsAngstroem);
    115     }
    116 
    117     // make sum of velocities equal zero
    118     if (FixedCenterOfMass)
    119       correctVelocitiesForFixedCenterOfMass(NextStep);
    120 
    121     // thermostat
    122     performThermostatControl(NextStep);
     129      (*iter)->VelocityVerletUpdateX((*iter)->getId(), NextStep, Deltat, IsAngstroem);
     130    }
    123131
    124132    // exit

  • tests/Python/AllActions/options.dat

    r11f0fa rbcb593  
    9191load    "test.data"
    9292MaxDistance     "-1"
    93 MDSteps "1"
    9493mesh-offset     "0.5,0.5,0.5"
    9594mesh-size       "10,10,10"

  • tests/regression/Molecules/VerletIntegration/post/test.conf

    r11f0fa rbcb593  
    7171Ion_Type2       3       6       1.0     3       3       12.01100000000  Carbon  C
    7272#Ion_TypeNr._Nr.R[0]    R[1]    R[2]    MoveType (0 MoveIon, 1 FixedIon)
    73 Ion_Type1_1     9.782085945     2.645886050     2.645886050     0       1.000000e-01    0.000000e+00    0.000000e+00     # molecule nr 0
    74 Ion_Type1_2     9.782085945     2.645886050     4.425886024     0       1.000000e-01    0.000000e+00    0.000000e+00     # molecule nr 1
    75 Ion_Type1_3     10.672039608    3.904536878     3.535886037     0       1.000000e-01    0.000000e+00    0.000000e+00     # molecule nr 2
    76 Ion_Type1_4     8.532785963     4.787886018     2.645886050     0       1.000000e-01    0.000000e+00    0.000000e+00     # molecule nr 3
    77 Ion_Type1_5     8.532785963     4.787886018     4.425886024     0       1.000000e-01    0.000000e+00    0.000000e+00     # molecule nr 4
    78 Ion_Type1_6     6.393632318     3.904536877     3.535886037     0       1.000000e-01    0.000000e+00    0.000000e+00     # molecule nr 5
    79 Ion_Type1_7     7.283585982     2.645886050     2.645886050     0       1.000000e-01    0.000000e+00    0.000000e+00     # molecule nr 6
    80 Ion_Type1_8     7.283585982     2.645886050     4.425886024     0       1.000000e-01    0.000000e+00    0.000000e+00     # molecule nr 7
    81 Ion_Type2_1     9.782085945     3.275186040     3.535886037     0       1.000000e-01    0.000000e+00    0.000000e+00     # molecule nr 8
    82 Ion_Type2_2     8.532785963     4.158586027     3.535886037     0       1.000000e-01    0.000000e+00    0.000000e+00     # molecule nr 9
    83 Ion_Type2_3     7.283585982     3.275186040     3.535886037     0       1.000000e-01    0.000000e+00    0.000000e+00     # molecule nr 10
     73Ion_Type1_1     9.782085945     2.645886050     2.645886050     0       1.003242e-01    5.239287e-04    0.000000e+00     # molecule nr 0
     74Ion_Type1_2     9.782085945     2.645886050     4.425886024     0       1.003242e-01    5.239287e-04    0.000000e+00     # molecule nr 1
     75Ion_Type1_3     10.672039608    3.904536878     3.535886037     0       1.003242e-01    5.239287e-04    0.000000e+00     # molecule nr 2
     76Ion_Type1_4     8.532785963     4.787886018     2.645886050     0       1.003242e-01    5.239287e-04    0.000000e+00     # molecule nr 3
     77Ion_Type1_5     8.532785963     4.787886018     4.425886024     0       1.003242e-01    5.239287e-04    0.000000e+00     # molecule nr 4
     78Ion_Type1_6     6.393632318     3.904536877     3.535886037     0       1.003242e-01    5.239287e-04    0.000000e+00     # molecule nr 5
     79Ion_Type1_7     7.283585982     2.645886050     2.645886050     0       1.003242e-01    5.239287e-04    0.000000e+00     # molecule nr 6
     80Ion_Type1_8     7.283585982     2.645886050     4.425886024     0       1.003242e-01    5.239287e-04    0.000000e+00     # molecule nr 7
     81Ion_Type2_1     9.782085945     3.275186040     3.535886037     0       9.984420e-02    4.396970e-05    0.000000e+00     # molecule nr 8
     82Ion_Type2_2     8.532785963     4.158586027     3.535886037     0       9.984420e-02    4.396970e-05    0.000000e+00     # molecule nr 9
     83Ion_Type2_3     7.283585982     3.275186040     3.535886037     0       9.984420e-02    4.396970e-05    0.000000e+00     # molecule nr 10
    8484
    85 Ion_Type1_1     9.783090905     2.645891010     2.645886050     0       1.002953e-01    4.950408e-04    0.000000e+00    # molecule nr 0
    86 Ion_Type1_2     9.783090905     2.645891010     4.425886024     0       1.002953e-01    4.950408e-04    0.000000e+00    # molecule nr 1
    87 Ion_Type1_3     10.673044568    3.904541838     3.535886037     0       1.002953e-01    4.950408e-04    0.000000e+00    # molecule nr 2
    88 Ion_Type1_4     8.533790923     4.787890978     2.645886050     0       1.002953e-01    4.950408e-04    0.000000e+00    # molecule nr 3
    89 Ion_Type1_5     8.533790923     4.787890978     4.425886024     0       1.002953e-01    4.950408e-04    0.000000e+00    # molecule nr 4
    90 Ion_Type1_6     6.394637278     3.904541837     3.535886037     0       1.002953e-01    4.950408e-04    0.000000e+00    # molecule nr 5
    91 Ion_Type1_7     7.284590942     2.645891010     2.645886050     0       1.002953e-01    4.950408e-04    0.000000e+00    # molecule nr 6
    92 Ion_Type1_8     7.284590942     2.645891010     4.425886024     0       1.002953e-01    4.950408e-04    0.000000e+00    # molecule nr 7
    93 Ion_Type2_1     9.783086361     3.275186456     3.535886037     0       9.984178e-02    4.154535e-05    0.000000e+00    # molecule nr 8
    94 Ion_Type2_2     8.533786379     4.158586443     3.535886037     0       9.984178e-02    4.154535e-05    0.000000e+00    # molecule nr 9
    95 Ion_Type2_3     7.284586398     3.275186456     3.535886037     0       9.984178e-02    4.154535e-05    0.000000e+00    # molecule nr 10
     85Ion_Type1_1     9.783091811     2.645893914     2.645886050     0 # molecule nr 0
     86Ion_Type1_2     9.783091811     2.645893914     4.425886024     0 # molecule nr 1
     87Ion_Type1_3     10.673045474    3.904544742     3.535886037     0 # molecule nr 2
     88Ion_Type1_4     8.533791829     4.787893882     2.645886050     0 # molecule nr 3
     89Ion_Type1_5     8.533791829     4.787893882     4.425886024     0 # molecule nr 4
     90Ion_Type1_6     6.394638184     3.904544741     3.535886037     0 # molecule nr 5
     91Ion_Type1_7     7.284591848     2.645893914     2.645886050     0 # molecule nr 6
     92Ion_Type1_8     7.284591848     2.645893914     4.425886024     0 # molecule nr 7
     93Ion_Type2_1     9.783084607     3.275186700     3.535886037     0 # molecule nr 8
     94Ion_Type2_2     8.533784625     4.158586687     3.535886037     0 # molecule nr 9
     95Ion_Type2_3     7.284584644     3.275186700     3.535886037     0 # molecule nr 10

  • tests/regression/Molecules/VerletIntegration/testsuite-molecules-verlet-integration.at

    r11f0fa rbcb593  
    2525AT_CHECK([/bin/cp -f ${abs_top_srcdir}/tests/regression/Molecules/VerletIntegration/pre/test.forces .], 0)
    2626AT_CHECK([chmod u+w $file], 0)
    27 AT_CHECK([../../molecuilder -i $file --select-all-atoms --verlet-integration --forces-file test.forces --MDSteps 1 --deltat 0.01 --keep-fixed-CenterOfMass 0], 0, [stdout], [stderr])
     27AT_CHECK([../../molecuilder -i $file --select-all-atoms --verlet-integration --forces-file test.forces --deltat 0.01 --keep-fixed-CenterOfMass 0], 0, [stdout], [stderr])
    2828AT_CHECK([diff $file ${abs_top_srcdir}/tests/regression/Molecules/VerletIntegration/post/test.conf], 0, [ignore], [ignore])
    2929
Note: See TracChangeset for help on using the changeset viewer.