Changeset 4a7776a for src/molecule_dynamics.cpp

Timestamp:

Oct 12, 2009, 10:30:02 AM (16 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, Candidate_v1.7.0, 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:

5034e1

Parents:

ccd9f5

Message:

Complete refactoring of molecule_dynamics.cpp

• new small functions: SumDistanceOfTrajectories(), PenalizeEqualTargets(), (), CreateInitialLists(), TryNextNearestNeighbourForInjectivePermutation(), MakeInjectivePermutation()
• Thermostat routines (case blocks) have been shifted to small functions in atom class (and soon to be extracted from there into its own class: ThermostatedParticle ?)
• old functions are a lot smaller now: MinimiseConstrainedPotential(), molecule::EvaluateConstrainedForces(), VerletForceIntegration(), also using ActOnAllAtoms
• atom class has a lot of new small functions:
  • trajectory: ResizeTrajectory(), CopyStepOnStep(), VelocityVerletUpdate(), VelocityVerletUpdate()
  • constraint potential: CorrectVelocity(), Thermostat_...()

File:

• src/molecule_dynamics.cpp (modified) (7 diffs)

src/molecule_dynamics.cpp

@@ -497 +497 @@
   else {
     PermutationMap = Malloc<atom *>(AtomCount, "molecule::LinearInterpolationBetweenConfiguration: **PermutationMap");
-    Walker = start;
-    while (Walker->next != end) {
-      Walker = Walker->next;
-      PermutationMap[Walker->nr] = Walker;   // always pick target with the smallest distance
-    }
+    SetIndexedArrayForEachAtomTo( PermutationMap, &atom::nr );
   }
 
   // check whether we have sufficient space in Trajectories for each atom
-  Walker = start;
-  while (Walker->next != end) {
-    Walker = Walker->next;
-    if (Walker->Trajectory.R.size() <= (unsigned int)(MaxSteps)) {
-      //cout << "Increasing size for trajectory array of " << keyword << " to " << (MaxSteps+1) << "." << endl;
-      Walker->Trajectory.R.resize(MaxSteps+1);
-      Walker->Trajectory.U.resize(MaxSteps+1);
-      Walker->Trajectory.F.resize(MaxSteps+1);
-    }
-  }
+  ActOnAllAtoms( &atom::ResizeTrajectory, MaxSteps );
   // push endstep to last one
-  Walker = start;
-  while (Walker->next != end) { // remove the endstep (is now the very last one)
-    Walker = Walker->next;
-    for (int n=NDIM;n--;) {
-      Walker->Trajectory.R.at(MaxSteps).x[n] = Walker->Trajectory.R.at(endstep).x[n];
-      Walker->Trajectory.U.at(MaxSteps).x[n] = Walker->Trajectory.U.at(endstep).x[n];
-      Walker->Trajectory.F.at(MaxSteps).x[n] = Walker->Trajectory.F.at(endstep).x[n];
-    }
-  }
+  ActOnAllAtoms( &atom::CopyStepOnStep, MaxSteps, endstep );
   endstep = MaxSteps;
 
@@ -578 +557 @@
 bool molecule::VerletForceIntegration(ofstream *out, char *file, config &configuration)
 {
-  atom *walker = NULL;
   ifstream input(file);
   string token;
   stringstream item;
-  double IonMass, Vector[NDIM], ConstrainedPotentialEnergy, ActualTemp;
+  double IonMass, ConstrainedPotentialEnergy, ActualTemp;
+  Vector Velocity;
   ForceMatrix Force;
 
@@ -601 +580 @@
     }
     // correct Forces
-    for(int d=0;d<NDIM;d++)
-      Vector[d] = 0.;
+    Velocity.Zero();
     for(int i=0;i<AtomCount;i++)
       for(int d=0;d<NDIM;d++) {
-        Vector[d] += Force.Matrix[0][i][d+5];
+        Velocity.x[d] += Force.Matrix[0][i][d+5];
       }
     for(int i=0;i<AtomCount;i++)
       for(int d=0;d<NDIM;d++) {
-        Force.Matrix[0][i][d+5] -= Vector[d]/(double)AtomCount;
+        Force.Matrix[0][i][d+5] -= Velocity.x[d]/(double)AtomCount;
       }
     // solve a constrained potential if we are meant to
@@ -621 +599 @@
 
     // and perform Verlet integration for each atom with position, velocity and force vector
-    walker = start;
-    while (walker->next != end) { // go through every atom of this element
-      walker = walker->next;
-      //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
-      // check size of vectors
-      if (walker->Trajectory.R.size() <= (unsigned int)(MDSteps)) {
-        //out << "Increasing size for trajectory array of " << *walker << " to " << (size+10) << "." << endl;
-        walker->Trajectory.R.resize(MDSteps+10);
-        walker->Trajectory.U.resize(MDSteps+10);
-        walker->Trajectory.F.resize(MDSteps+10);
-      }
-
-      // Update R (and F)
-      for (int d=0; d<NDIM; d++) {
-        walker->Trajectory.F.at(MDSteps).x[d] = -Force.Matrix[0][walker->nr][d+5]*(configuration.GetIsAngstroem() ? AtomicLengthToAngstroem : 1.);
-        walker->Trajectory.R.at(MDSteps).x[d] = walker->Trajectory.R.at(MDSteps-1).x[d];
-        walker->Trajectory.R.at(MDSteps).x[d] += configuration.Deltat*(walker->Trajectory.U.at(MDSteps-1).x[d]);     // s(t) = s(0) + v * deltat + 1/2 a * deltat^2
-        walker->Trajectory.R.at(MDSteps).x[d] += 0.5*configuration.Deltat*configuration.Deltat*(walker->Trajectory.F.at(MDSteps).x[d]/walker->type->mass);     // F = m * a and s = 0.5 * F/m * t^2 = F * a * t
-      }
-      // Update U
-      for (int d=0; d<NDIM; d++) {
-        walker->Trajectory.U.at(MDSteps).x[d] = walker->Trajectory.U.at(MDSteps-1).x[d];
-        walker->Trajectory.U.at(MDSteps).x[d] += configuration.Deltat * (walker->Trajectory.F.at(MDSteps).x[d]+walker->Trajectory.F.at(MDSteps-1).x[d]/walker->type->mass); // v = F/m * t
-      }
-//      out << "Integrated position&velocity of step " << (MDSteps) << ": (";
-//      for (int d=0;d<NDIM;d++)
-//        out << walker->Trajectory.R.at(MDSteps).x[d] << " ";          // next step
-//      out << ")\t(";
-//      for (int d=0;d<NDIM;d++)
-//        cout << walker->Trajectory.U.at(MDSteps).x[d] << " ";          // next step
-//      out << ")" << endl;
-            // next atom
-    }
+    // check size of vectors
+    ActOnAllAtoms( &atom::ResizeTrajectory, MDSteps+10 );
+
+    ActOnAllAtoms( &atom::VelocityVerletUpdate, MDSteps, &configuration, &Force);
   }
   // correct velocities (rather momenta) so that center of mass remains motionless
-  for(int d=0;d<NDIM;d++)
-    Vector[d] = 0.;
+  Velocity.Zero();
   IonMass = 0.;
-  walker = start;
-  while (walker->next != end) { // go through every atom
-    walker = walker->next;
-    IonMass += walker->type->mass;  // sum up total mass
-    for(int d=0;d<NDIM;d++) {
-      Vector[d] += walker->Trajectory.U.at(MDSteps).x[d]*walker->type->mass;
-    }
-  }
+  ActOnAllAtoms ( &atom::SumUpKineticEnergy, MDSteps, &IonMass, &Velocity );
+
   // correct velocities (rather momenta) so that center of mass remains motionless
-  for(int d=0;d<NDIM;d++)
-    Vector[d] /= IonMass;
+  Velocity.Scale(1./IonMass);
   ActualTemp = 0.;
-  walker = start;
-  while (walker->next != end) { // go through every atom of this element
-    walker = walker->next;
-    for(int d=0;d<NDIM;d++) {
-      walker->Trajectory.U.at(MDSteps).x[d] -= Vector[d];
-      ActualTemp += 0.5 * walker->type->mass * walker->Trajectory.U.at(MDSteps).x[d] * walker->Trajectory.U.at(MDSteps).x[d];
-    }
-  }
+  ActOnAllAtoms ( &atom::CorrectVelocity, &ActualTemp, MDSteps, &Velocity );
   Thermostats(configuration, ActualTemp, Berendsen);
   MDSteps++;
-
 
   // exit
@@ -712 +645 @@
   double delta_alpha = 0.;
   double ScaleTempFactor;
-  double sigma;
-  double IonMass;
-  int d;
   gsl_rng * r;
   const gsl_rng_type * T;
-  double *U = NULL, *F = NULL, FConstraint[NDIM];
-  atom *walker = NULL;
 
   // calculate scale configuration
@@ -731 +659 @@
       if ((configuration.ScaleTempStep > 0) && ((MDSteps-1) % configuration.ScaleTempStep == 0)) {
         cout << Verbose(2) <<  "Applying Woodcock thermostat..." << endl;
-        walker = start;
-        while (walker->next != end) { // go through every atom of this element
-          walker = walker->next;
-          IonMass = walker->type->mass;
-          U = walker->Trajectory.U.at(MDSteps).x;
-          if (walker->FixedIon == 0) // even FixedIon moves, only not by other's forces
-            for (d=0; d<NDIM; d++) {
-              U[d] *= sqrt(ScaleTempFactor);
-              ekin += 0.5*IonMass * U[d]*U[d];
-            }
-        }
+        ActOnAllAtoms( &atom::Thermostat_Woodcock, sqrt(ScaleTempFactor), MDSteps, &ekin );
       }
       break;
      case Gaussian:
       cout << Verbose(2) <<  "Applying Gaussian thermostat..." << endl;
-      walker = start;
-      while (walker->next != end) { // go through every atom of this element
-        walker = walker->next;
-        IonMass = walker->type->mass;
-        U = walker->Trajectory.U.at(MDSteps).x;
-        F = walker->Trajectory.F.at(MDSteps).x;
-        if (walker->FixedIon == 0) // even FixedIon moves, only not by other's forces
-          for (d=0; d<NDIM; d++) {
-            G += U[d] * F[d];
-            E += U[d]*U[d]*IonMass;
-          }
-      }
+      ActOnAllAtoms( &atom::Thermostat_Gaussian_init, MDSteps, &G, &E );
+
       cout << Verbose(1) << "Gaussian Least Constraint constant is " << G/E << "." << endl;
-      walker = start;
-      while (walker->next != end) { // go through every atom of this element
-        walker = walker->next;
-        IonMass = walker->type->mass;
-        U = walker->Trajectory.U.at(MDSteps).x;
-        F = walker->Trajectory.F.at(MDSteps).x;
-        if (walker->FixedIon == 0) // even FixedIon moves, only not by other's forces
-          for (d=0; d<NDIM; d++) {
-            FConstraint[d] = (G/E) * (U[d]*IonMass);
-            U[d] += configuration.Deltat/IonMass * (FConstraint[d]);
-            ekin += IonMass * U[d]*U[d];
-          }
-      }
+      ActOnAllAtoms( &atom::Thermostat_Gaussian_least_constraint, MDSteps, G/E, &ekin, &configuration);
+
       break;
      case Langevin:
@@ -780 +677 @@
       r = gsl_rng_alloc (T);
       // Go through each ion
-      walker = start;
-      while (walker->next != end) { // go through every atom of this element
-        walker = walker->next;
-        IonMass = walker->type->mass;
-        sigma  = sqrt(configuration.TargetTemp/IonMass); // sigma = (k_b T)/m (Hartree/atomicmass = atomiclength/atomictime)
-        U = walker->Trajectory.U.at(MDSteps).x;
-        F = walker->Trajectory.F.at(MDSteps).x;
-        if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
-          // throw a dice to determine whether it gets hit by a heat bath particle
-          if (((((rand()/(double)RAND_MAX))*configuration.TempFrequency) < 1.)) {
-            cout << Verbose(3) << "Particle " << *walker << " was hit (sigma " << sigma << "): " << sqrt(U[0]*U[0]+U[1]*U[1]+U[2]*U[2]) << " -> ";
-            // pick three random numbers from a Boltzmann distribution around the desired temperature T for each momenta axis
-            for (d=0; d<NDIM; d++) {
-              U[d] = gsl_ran_gaussian (r, sigma);
-            }
-            cout << sqrt(U[0]*U[0]+U[1]*U[1]+U[2]*U[2]) << endl;
-          }
-          for (d=0; d<NDIM; d++)
-            ekin += 0.5*IonMass * U[d]*U[d];
-        }
-      }
+      ActOnAllAtoms( &atom::Thermostat_Langevin, MDSteps, r, &ekin, &configuration );
       break;
+
      case Berendsen:
       cout << Verbose(2) <<  "Applying Berendsen-VanGunsteren thermostat..." << endl;
-      walker = start;
-      while (walker->next != end) { // go through every atom of this element
-        walker = walker->next;
-        IonMass = walker->type->mass;
-        U = walker->Trajectory.U.at(MDSteps).x;
-        F = walker->Trajectory.F.at(MDSteps).x;
-        if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
-          for (d=0; d<NDIM; d++) {
-            U[d] *= sqrt(1+(configuration.Deltat/configuration.TempFrequency)*(ScaleTempFactor-1));
-            ekin += 0.5*IonMass * U[d]*U[d];
-          }
-        }
-      }
+      ActOnAllAtoms( &atom::Thermostat_Berendsen, MDSteps, ScaleTempFactor, &ekin, &configuration );
       break;
+
      case NoseHoover:
       cout << Verbose(2) <<  "Applying Nose-Hoover thermostat..." << endl;
       // dynamically evolve alpha (the additional degree of freedom)
       delta_alpha = 0.;
-      walker = start;
-      while (walker->next != end) { // go through every atom of this element
-        walker = walker->next;
-        IonMass = walker->type->mass;
-        U = walker->Trajectory.U.at(MDSteps).x;
-        if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
-          for (d=0; d<NDIM; d++) {
-            delta_alpha += U[d]*U[d]*IonMass;
-          }
-        }
-      }
+      ActOnAllAtoms( &atom::Thermostat_NoseHoover_init, MDSteps, &delta_alpha );
       delta_alpha = (delta_alpha - (3.*AtomCount+1.) * configuration.TargetTemp)/(configuration.HooverMass*Units2Electronmass);
       configuration.alpha += delta_alpha*configuration.Deltat;
       cout << Verbose(3) << "alpha = " << delta_alpha << " * " << configuration.Deltat << " = " << configuration.alpha << "." << endl;
       // apply updated alpha as additional force
-      walker = start;
-      while (walker->next != end) { // go through every atom of this element
-        walker = walker->next;
-        IonMass = walker->type->mass;
-        U = walker->Trajectory.U.at(MDSteps).x;
-        if (walker->FixedIon == 0) { // even FixedIon moves, only not by other's forces
-          for (d=0; d<NDIM; d++) {
-              FConstraint[d] = - configuration.alpha * (U[d] * IonMass);
-              U[d] += configuration.Deltat/IonMass * (FConstraint[d]);
-              ekin += (0.5*IonMass) * U[d]*U[d];
-            }
-        }
-      }
+      ActOnAllAtoms( &atom::Thermostat_NoseHoover_scale, MDSteps, &ekin, &configuration );
       break;
   }