Changeset 68cfd1 for molecuilder/src/molecules.cpp

Timestamp:

Sep 28, 2008, 5:55:04 PM (17 years ago)

Author:

Frederik Heber <heber@…>

Children:

a7b3b8

Parents:

556157

Message:

corrected molecule::VerletForceIntegration() and analyzer creates data and plot file for shielding values

File:

• molecuilder/src/molecules.cpp (modified) (4 diffs)

molecuilder/src/molecules.cpp

@@ -1493 +1493 @@
   string token;
   stringstream item;
-  double a, IonMass, Vector[NDIM], ConstrainedPotentialEnergy, ActualTemp = 0.;
+  double a, IonMass, Vector[NDIM], ConstrainedPotentialEnergy, ActualTemp;
   ForceMatrix Force;
 
@@ -1535 +1535 @@
     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
+      //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 (Trajectories[walker].R.size() <= (unsigned int)(MDSteps)) {
@@ -1546 +1546 @@
       // Update R (and F)
       for (int d=0; d<NDIM; d++) {
-        Trajectories[walker].F.at(MDSteps).x[d] = -Force.Matrix[0][AtomNo][d+5]*(configuration.GetIsAngstroem() ? AtomicLengthToAngstroem : 1.);
+        Trajectories[walker].F.at(MDSteps).x[d] = -Force.Matrix[0][walker->nr][d+5]*(configuration.GetIsAngstroem() ? AtomicLengthToAngstroem : 1.);
         Trajectories[walker].R.at(MDSteps).x[d] = Trajectories[walker].R.at(MDSteps-1).x[d];
-        Trajectories[walker].R.at(MDSteps).x[d] += configuration.Deltat*(Trajectories[walker].U.at(MDSteps-1).x[d]);
-        Trajectories[walker].R.at(MDSteps).x[d] += configuration.Deltat*a*(Trajectories[walker].F.at(MDSteps).x[d]);     // F = m * a and s = 0.5 * F/m * t^2 = F * a * t
+        Trajectories[walker].R.at(MDSteps).x[d] += configuration.Deltat*(Trajectories[walker].U.at(MDSteps-1).x[d]);     // s(t) = s(0) + v * deltat + 1/2 a * deltat^2 
+        Trajectories[walker].R.at(MDSteps).x[d] += 0.5*configuration.Deltat*configuration.Deltat*(Trajectories[walker].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++) {
         Trajectories[walker].U.at(MDSteps).x[d] = Trajectories[walker].U.at(MDSteps-1).x[d]; 
-        Trajectories[walker].U.at(MDSteps).x[d] += a * (Trajectories[walker].F.at(MDSteps).x[d]+Trajectories[walker].F.at(MDSteps-1).x[d]);
+        Trajectories[walker].U.at(MDSteps).x[d] += configuration.Deltat * (Trajectories[walker].F.at(MDSteps).x[d]+Trajectories[walker].F.at(MDSteps-1).x[d]/walker->type->mass); // v = F/m * t
       }
 //      out << "Integrated position&velocity of step " << (MDSteps) << ": ("; 
@@ -1578 +1578 @@
     }
   }
+  for(int d=0;d<NDIM;d++)
+    Vector[d] /= 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++) {
-      Trajectories[walker].U.at(MDSteps).x[d] -= Vector[d]*walker->type->mass/IonMass;
+      Trajectories[walker].U.at(MDSteps).x[d] -= Vector[d];
       ActualTemp += 0.5 * walker->type->mass * Trajectories[walker].U.at(MDSteps).x[d] * Trajectories[walker].U.at(MDSteps).x[d];
     }