Index: pcp/src/ions.c
===================================================================
--- pcp/src/ions.c	(revision 9275a309437b2d8525beb57bacd29efdecf8278a)
+++ pcp/src/ions.c	(revision 32de28ca6726b18cb1f6631dcc610e6c1a791143)
@@ -61,6 +61,6 @@
   //struct RunStruct *Run = &P->R;
   int i,it,j,IMT,BorAng, d,me, count;
+  struct Lattice *L = &P->Lat;
   int k;
-  struct Lattice *L = &P->Lat;
   int relative; // whether read-in coordinate are relative (1) or absolute  
   double Bohr = ANGSTROEMINBOHRRADIUS; /* Angstroem */
@@ -183,9 +183,50 @@
     }
   }
-  Free(free_name);
+  Free(free_name, "IonsInitRead: free_name");
   I->Max_Max_IonsOfType = 0;
   for (i=0; i < I->Max_Types; i++) 
     I->Max_Max_IonsOfType = MAX(I->Max_Max_IonsOfType, I->I[i].Max_IonsOfType);
   I->FTemp = (double *) Malloc(sizeof(double) * NDIM * I->Max_Max_IonsOfType, "IonsInitRead:");
+}
+
+/** Allocates memory for a IonConstrained item and adds to end of list.
+ * \param *I IonType structure
+ * \param ion which ion of the type
+ * \return pointer to created item
+ */
+struct IonConstrained * AddConstraintItem(struct IonType *I, int ion)
+{
+  struct IonConstrained **ptr = &(I->Constraints[ion]);
+  while (*ptr != NULL) { // advance to end of list
+    ptr = &((*ptr)->next);
+  }
+  // allocated memory for structure and items within
+  *ptr = (struct IonConstrained *) Malloc(sizeof(struct IonConstrained), "CreateConstraintItem: IonConstrained");
+  (*ptr)->R = (double *) Malloc(sizeof(double)*NDIM, "AddConstraintItem: *R");
+  (*ptr)->U = (double *) Malloc(sizeof(double)*NDIM, "AddConstraintItem: *U");
+  (*ptr)->next = NULL;
+  return *ptr;
+}
+
+/** Removes first of item of IonConstrained list for given \a ion.
+ * Frees memory of items within and then
+ * \param *I IonType structure
+ * \param ion which ion of the type
+ * \return 1 - success, 0 - failure (list is already empty, start pointer was NULL)
+ */ 
+int RemoveConstraintItem(struct IonType *I, int ion)
+{
+  struct IonConstrained **ptr = &(I->Constraints[ion]); 
+  struct IonConstrained *next_ptr;
+  if (*ptr != NULL) {
+    next_ptr = (*ptr)->next;
+    Free((*ptr)->R, "RemoveConstraintItem: R");
+    Free((*ptr)->U, "RemoveConstraintItem: U");
+    Free((*ptr), "RemoveConstraintItem: IonConstrained structure");
+    *ptr = next_ptr; // set start of list to next item (automatically is null if ptr was last item)
+    return 1;
+  } else {
+    return 0;
+  }
 }
 
@@ -386,61 +427,39 @@
   int i, it;
   for (i=0; i < I->Max_Types; i++) {
-    //fprintf(stderr,"Name\n");
-    Free(I->I[i].Name);
-    //fprintf(stderr,"Abbreviation\n");
-    Free(I->I[i].Symbol); 
+    Free(I->I[i].Name, "RemoveIonsRead: I->I[i].Name");
+    Free(I->I[i].Symbol, "RemoveIonsRead: I->I[i].Symbol"); 
     for (it=0;it<I->I[i].Max_IonsOfType;it++) {
-    //fprintf(stderr,"sigma[it]\n");
-      Free(I->I[i].sigma[it]);
-    //fprintf(stderr,"sigma_rezi[it]\n");
-      Free(I->I[i].sigma_rezi[it]);
-    //fprintf(stderr,"sigma_PAS[it]\n");
-      Free(I->I[i].sigma_PAS[it]);
-    //fprintf(stderr,"sigma_rezi_PAS[it]\n");
-      Free(I->I[i].sigma_rezi_PAS[it]);
-    }
-    //fprintf(stderr,"sigma\n");
-    Free(I->I[i].sigma);
-    //fprintf(stderr,"sigma_rezi\n");
-    Free(I->I[i].sigma_rezi);
-    //fprintf(stderr,"sigma_PAS\n");
-    Free(I->I[i].sigma_PAS);
-    //fprintf(stderr,"sigma_rezi_PAS\n");
-    Free(I->I[i].sigma_rezi_PAS);
-    //fprintf(stderr,"R\n");
-    Free(I->I[i].R);
-    //fprintf(stderr,"R_old\n");
-    Free(I->I[i].R_old);
-    //fprintf(stderr,"R_old_old\n");
-    Free(I->I[i].R_old_old);
-    //fprintf(stderr,"FIon\n");
-    Free(I->I[i].FIon);
-    //fprintf(stderr,"FIon_old\n");
-    Free(I->I[i].FIon_old);
-    //fprintf(stderr,"SearchDir\n");
-    Free(I->I[i].SearchDir);
-    //fprintf(stderr,"GammaA\n");
-    Free(I->I[i].GammaA);
-    //fprintf(stderr,"FIonL\n");
-    Free(I->I[i].FIonL);
-    //fprintf(stderr,"FIonNL\n");
-    Free(I->I[i].FIonNL);
-    //fprintf(stderr,"U\n");
-    Free(I->I[i].U);
-    //fprintf(stderr,"SFactor\n");
-    Free(I->I[i].SFactor);
-    //fprintf(stderr,"IMT\n");
-    Free(I->I[i].IMT);
-    //fprintf(stderr,"FEwald\n");
-    Free(I->I[i].FEwald);
-    Free(I->I[i].alpha);
-  }
-  //fprintf(stderr,"RLatticeVec\n");
+      Free(I->I[i].moment[it], "RemoveIonsRead: I->I[i].moment[it]");
+      Free(I->I[i].sigma[it], "RemoveIonsRead: I->I[i].sigma[it]");
+      Free(I->I[i].sigma_rezi[it], "RemoveIonsRead: I->I[i].sigma_rezi[it]");
+      Free(I->I[i].sigma_PAS[it], "RemoveIonsRead: I->I[i].sigma_PAS[it]");
+      Free(I->I[i].sigma_rezi_PAS[it], "RemoveIonsRead: I->I[i].sigma_rezi_PAS[it]");
+      while (RemoveConstraintItem(&I->I[i], it)); // empty the constrained item list
+    }
+    Free(I->I[i].sigma, "RemoveIonsRead: I->I[i].sigma");
+    Free(I->I[i].sigma_rezi, "RemoveIonsRead: I->I[i].sigma_rezi");
+    Free(I->I[i].sigma_PAS, "RemoveIonsRead: I->I[i].sigma_PAS");
+    Free(I->I[i].sigma_rezi_PAS, "RemoveIonsRead: I->I[i].sigma_rezi_PAS");
+    Free(I->I[i].R, "RemoveIonsRead: I->I[i].R");
+    Free(I->I[i].R_old, "RemoveIonsRead: I->I[i].R_old");
+    Free(I->I[i].R_old_old, "RemoveIonsRead: I->I[i].R_old_old");
+    Free(I->I[i].Constraints, "RemoveIonsRead: I->I[i].Constraints");
+    Free(I->I[i].FIon, "RemoveIonsRead: I->I[i].FIon");
+    Free(I->I[i].FIon_old, "RemoveIonsRead: I->I[i].FIon_old");
+    Free(I->I[i].SearchDir, "RemoveIonsRead: I->I[i].SearchDir");
+    Free(I->I[i].GammaA, "RemoveIonsRead: I->I[i].GammaA");
+    Free(I->I[i].FIonL, "RemoveIonsRead: I->I[i].FIonL");
+    Free(I->I[i].FIonNL, "RemoveIonsRead: I->I[i].FIonNL");
+    Free(I->I[i].FMagnetic, "RemoveIonsRead: I->I[i].FMagnetic");
+    Free(I->I[i].U, "RemoveIonsRead: I->I[i].U");
+    Free(I->I[i].SFactor, "RemoveIonsRead: I->I[i].SFactor");
+    Free(I->I[i].IMT, "RemoveIonsRead: I->I[i].IMT");
+    Free(I->I[i].FEwald, "RemoveIonsRead: I->I[i].FEwald");
+    Free(I->I[i].alpha, "RemoveIonsRead: I->I[i].alpha");
+  }
   if (I->R_cut == 0.0)
-    Free(I->RLatticeVec);
-  //fprintf(stderr,"FTemp\n");
-  Free(I->FTemp);
-  //fprintf(stderr,"I\n");
-  Free(I->I);
+    Free(I->RLatticeVec, "RemoveIonsRead: I->RLatticeVec");
+  Free(I->FTemp, "RemoveIonsRead: I->FTemp");
+  Free(I->I, "RemoveIonsRead: I->I");
 }
 
Index: pcp/src/ions.h
===================================================================
--- pcp/src/ions.h	(revision 9275a309437b2d8525beb57bacd29efdecf8278a)
+++ pcp/src/ions.h	(revision 32de28ca6726b18cb1f6631dcc610e6c1a791143)
@@ -38,4 +38,16 @@
 enum thermostats { None, Woodcock, Gaussian, Langevin, Berendsen, NoseHoover };   //!< Thermostat names for output
 
+/** Structure containing coordinate and velocity of one ion for a single MD step.
+ * This is an item in a forward-referenced list, such that at each step the current
+ * constrained can be extracted and the item subsequently removed from the list.
+ */
+struct IonConstrained {
+  int step;           //!< constrained motion for MD step ...
+  double *R;          //!< ion coordinates (null if not specified)
+  double *U;          //!< ion velocity (null if not specified)
+  enum IonMoveType IMT;  //!< Ion MoveType for this step
+  struct IonConstrained *next;  //!< pointer to next item in forward-referenced list
+};
+
 /** Structure containing data for a single ion type.
  * Contains Mass, atomic number, number of types, coordinates (current and the
@@ -48,13 +60,15 @@
   int Max_IonsOfType;	//!< Maximum number of ions of this type (needed for memory allocation)
   double *R;					//!< coordinate (x,y,z) array for all ions of this type R[i][j] (i=x,y,z;j=1,..)
-  double *R_old;			//!< (?) Remembers last position (two-dim array, see IonType::R)
-  double *R_old_old;	//!< (?) Remembers position before the last one (two-dim array, see IonType::R)
-  double *FIon;				//!< (?) Overall force acting on the ion (two-dim array, see IonType::R)
-  double *FIonL;			//!< (?) Linear force acting on the ion (two-dim array, see IonType::R)
-  double *FIonNL;			//!< (?) Nonlinear force acting on the ion (two-dim array, see IonType::R)
+  double *R_old;			//!< Remembers last position (two-dim array, see IonType::R)
+  double *R_old_old;	//!< Remembers position before the last one (two-dim array, see IonType::R)
+  struct IonConstrained **Constraints; //!< pointer to each Ion of this IonType with a pointer to the current start item of a list of constrained values
+  double *FIon;				//!< Overall force acting on the ion (two-dim array, see IonType::R)
+  double *FIonL;			//!< Linear force acting on the ion (two-dim array, see IonType::R)
+  double *FIonNL;			//!< Nonlinear force acting on the ion (two-dim array, see IonType::R)
+  double *FEwald;     //!< Ewald force (two-dim array, see IonType::R)
+  double *FMagnetic;  //!< Magnetic force acting on the ion (two-dim array, see IonType::R)
+  double *FConstraint;//!< Constraint force acting on the ion, needed for Thermostat() (two-dim array, see IonType::R)
   double *U; 					//!< Velocity of the ion (two-dim array, see IonType::R)
-  double *FIon_old;		//!< (?) Remembers prior force acting on the ion (two-dim array, see IonType::R)
-  double *FEwald;			//!< (?) Ewald force (two-dim array, see IonType::R)
-  double *FConstraint;//!< Constraint force acting on the ion, needed for Thermostat() (two-dim array, see IonType::R)  
+  double *FIon_old;		//!< Remembers prior force acting on the ion (two-dim array, see IonType::R)
   double IonFac;      //!< mass-like factor in CG structure optimization to determinte step width
   double *GammaA;     //!< 
@@ -70,4 +84,5 @@
   char *Name;         //!< name of ion's element
   char *Symbol; //!< short form of ion's element
+  double **moment;     //!< magnetic moment
   double **sigma;     //!< shielding tensor
   double **sigma_rezi;      //!< shielding tensor for reciprocal calculation
@@ -106,4 +121,5 @@
 /* Functions */
 void IonsInitRead(struct Problem *P, FILE *source);
+void InitThermostats(struct Problem *P, FILE *source);
 void CalculateEwald(struct Problem *P, int first);
 void RemoveIonsRead(struct Ions *I);
@@ -111,4 +127,5 @@
 void CalculateIonForce(struct Problem *P);
 void OutputIonForce(struct Problem *P);
+void ParseIonForce(struct Problem *P);
 void OutputIonCoordinates(struct Problem *P);
 void UpdateIons(struct Problem *P);
