Changeset ef7d30

Timestamp:

Jun 29, 2010, 1:04:40 PM (15 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:

8e0c63

Parents:

52d777

Message:

MEMFIX: ForceMatrix in VelocityVerletIntegration had a stupid offset.

• this offset is now generalized and needed as an additional parameter in TrajectoryParticle::VelocityVerletUpdate() and molecule::VerletForceIntegration()

Location:

src

Files:

• Actions/MoleculeAction/VerletIntegrationAction.cpp (modified) (1 diff)
• atom_trajectoryparticle.cpp (modified) (1 diff)
• atom_trajectoryparticle.hpp (modified) (1 diff)
• molecule.hpp (modified) (1 diff)
• molecule_dynamics.cpp (modified) (5 diffs)

src/Actions/MoleculeAction/VerletIntegrationAction.cpp

@@ -62 +62 @@
     char outputname[MAXSTRINGSIZE];
     strcpy(outputname, filename.c_str());
-    if (!mol->VerletForceIntegration(outputname, *(World::getInstance().getConfig())))
+    if (!mol->VerletForceIntegration(outputname, *(World::getInstance().getConfig()), 0))
       DoLog(2) && (Log() << Verbose(2) << "File not found." << endl);
     else

src/atom_trajectoryparticle.cpp

@@ -106 +106 @@
  * \param *Force matrix with forces
  */
-void TrajectoryParticle::VelocityVerletUpdate(int NextStep, config *configuration, ForceMatrix *Force)
+void TrajectoryParticle::VelocityVerletUpdate(int NextStep, config *configuration, ForceMatrix *Force, const size_t offset)
 {
   //a = configuration.Deltat*0.5/walker->type->mass;        // (F+F_old)/2m = a and thus: v = (F+F_old)/2m * t = (F + F_old) * a
   for (int d=0; d<NDIM; d++) {
-    Trajectory.F.at(NextStep)[d] = -Force->Matrix[0][nr][d+5]*(configuration->GetIsAngstroem() ? AtomicLengthToAngstroem : 1.);
+    Trajectory.F.at(NextStep)[d] = -Force->Matrix[0][nr][d+offset]*(configuration->GetIsAngstroem() ? AtomicLengthToAngstroem : 1.);
     Trajectory.R.at(NextStep)[d] = Trajectory.R.at(NextStep-1)[d];
     Trajectory.R.at(NextStep)[d] += configuration->Deltat*(Trajectory.U.at(NextStep-1)[d]);     // s(t) = s(0) + v * deltat + 1/2 a * deltat^2

src/atom_trajectoryparticle.hpp

@@ -50 +50 @@
   void ResizeTrajectory(int MaxSteps);
   void CopyStepOnStep(int dest, int src);
-  void VelocityVerletUpdate(int MDSteps, config *configuration, ForceMatrix *Force);
+  void VelocityVerletUpdate(int MDSteps, config *configuration, ForceMatrix *Force, const size_t offset);
   void SumUpKineticEnergy( int Step, double *TotalMass, Vector *TotalVelocity ) const;
 

src/molecule.hpp

@@ -279 +279 @@
   void SetBoxDimension(Vector *dim);
   void ScanForPeriodicCorrection();
-  bool VerletForceIntegration(char *file, config &configuration);
+  bool VerletForceIntegration(char *file, config &configuration, const size_t offset);
   void Thermostats(config &configuration, double ActualTemp, int Thermostat);
   void PrincipalAxisSystem(bool DoRotate);

src/molecule_dynamics.cpp

@@ -540 +540 @@
  * have to transform them).
  * This adds a new MD step to the config file.
- * \param *out output stream for debugging
  * \param *file filename
  * \param config structure with config::Deltat, config::IsAngstroem, config::DoConstrained
- * \param delta_t time step width in atomic units
- * \param IsAngstroem whether coordinates are in angstroem (true) or bohrradius (false)
- * \param DoConstrained whether we perform a constrained (>0, target step in molecule::trajectories) or unconstrained (0) molecular dynamics, \sa molecule::MinimiseConstrainedPotential()
+ * \param offset offset in matrix file to the first force component
  * \return true - file found and parsed, false - file not found or imparsable
  * \todo This is not yet checked if it is correctly working with DoConstrained set to true.
  */
-bool molecule::VerletForceIntegration(char *file, config &configuration)
+bool molecule::VerletForceIntegration(char *file, config &configuration, const size_t offset)
 {
   Info FunctionInfo(__func__);
@@ -561 +558 @@
   CountElements();  // make sure ElementsInMolecule is up to date
 
+  const int AtomCount = getAtomCount();
   // check file
   if (input == NULL) {
@@ -571 +569 @@
       return false;
     }
-    if (Force.RowCounter[0] != getAtomCount()) {
+    if (Force.RowCounter[0] != AtomCount) {
       DoeLog(0) && (eLog()<< Verbose(0) << "Mismatch between number of atoms in file " << Force.RowCounter[0] << " and in molecule " << getAtomCount() << "." << endl);
       performCriticalExit();
@@ -578 +576 @@
     // correct Forces
     Velocity.Zero();
-    for(int i=0;i<getAtomCount();i++)
+    for(int i=0;i<AtomCount;i++)
       for(int d=0;d<NDIM;d++) {
-        Velocity[d] += Force.Matrix[0][i][d+5];
+        Velocity[d] += Force.Matrix[0][i][d+offset];
       }
-    for(int i=0;i<getAtomCount();i++)
+    for(int i=0;i<AtomCount;i++)
       for(int d=0;d<NDIM;d++) {
-        Force.Matrix[0][i][d+5] -= Velocity[d]/static_cast<double>(getAtomCount());
+        Force.Matrix[0][i][d+offset] -= Velocity[d]/static_cast<double>(AtomCount);
       }
     // solve a constrained potential if we are meant to
@@ -599 +597 @@
     //ActOnAllAtoms( &atom::ResizeTrajectory, MDSteps+10 );
 
-    ActOnAllAtoms( &atom::VelocityVerletUpdate, MDSteps+1, &configuration, &Force);
+    ActOnAllAtoms( &atom::VelocityVerletUpdate, MDSteps+1, &configuration, &Force, (const size_t) 0);
   }
   // correct velocities (rather momenta) so that center of mass remains motionless