Changeset 6590194 for pcp/src

Timestamp:

Apr 21, 2008, 2:19:26 PM (17 years ago)

Author:

Frederik Heber <heber@…>

Children:

05fec0

Parents:

4782e78

git-author:

Frederik Heber <heber@…> (04/18/08 16:18:38)

git-committer:

Frederik Heber <heber@…> (04/21/08 14:19:26)

Message:

UpdateIonsR(): included constraint motion, and don't reset Forces afterwards, as they are needed for OutputVis()

File:

• pcp/src/ions.c (modified) (2 diffs)

pcp/src/ions.c ¶

@@ -853 +853 @@
       R_old = &I->I[is].R_old[NDIM*ia];
       R_old_old = &I->I[is].R_old_old[NDIM*ia];
-      if (I->I[is].IMT[ia] == MoveIon) {
+      //if (I->I[is].IMT[ia] == MoveIon) { // even FixedIon moves, only not by other's forces
                                 for (d=0; d<NDIM; d++) {
                                   R_old_old[d] = R_old[d];                              // shift old values
                                   R_old[d] = R[d];
-                                  R[d] += delta_t*(U[d]+a*FIon[d]);                     // F = m * a and s = 0.5 * a * t^2
+                            R[d] += delta_t*(U[d]+a*FIon[d]);                   // F = m * a and s = 0.5 * a * t^2
                                   FIon_old[d] = FIon[d];                                        // store old force
-                                  FIon[d] = 0;                                                                          // zero all as a sign that's been moved
-                                  FIonL[d] = 0;
-                                  FIonNL[d] = 0;
-                                  FEwald[d] = 0;
+//                                FIon[d] = 0.;                                                                         // zero all as a sign that's been moved
+//                                FIonL[d] = 0.;
+//                                FIonNL[d] = 0.;
+//                                FEwald[d] = 0.;
+// don't reset forces here anymore, as we OutputVis() after this and before next CalculateForce() (rendering them null in the cute graphs)
                                 }
+        if (I->I[is].Constraints[ia] != NULL) { // override current resulting R if constrained is given
+          fprintf(stderr, "(%i) Using constrained coordinates R of step %d on ion (%i,%i)\n", P->Par.me, I->I[is].Constraints[ia]->step, is, ia);
+          if ((I->I[is].Constraints[ia]->step != P->R.OuterStep) && (P->Par.me == 0))
+            fprintf(stderr, "(%i) WARNING: constraint step %d does not match Outerstep %d for ion (%d,%d)\n", P->Par.me, I->I[is].Constraints[ia]->step, P->R.OuterStep, is, ia);
+          for (d=0; d<NDIM; d++)
+            R[d] = I->I[is].Constraints[ia]->R[d];
+          I->I[is].IMT[ia] = I->I[is].Constraints[ia]->IMT;
+          //RemoveConstraintItem(&I->I[is],ia); // and remove this first item from the list: // is done now  at UpdateIonsU()
+        }
                                 RMat33Vec3(sR, P->Lat.InvBasis, R);
                                 RMat33Vec3(sRold, P->Lat.InvBasis, R_old);
@@ -882 +892 @@
                                 RMat33Vec3(R_old, P->Lat.RealBasis, sRold);
                                 RMat33Vec3(R_old_old, P->Lat.RealBasis, sRoldold);
-      }
+      //}
     }
   }