Index: pcp/src/init.c
===================================================================
--- pcp/src/init.c	(revision e00f4732f333881727b6265e9f995246359ea74f)
+++ pcp/src/init.c	(revision 961b75c543d6081b2758726da4c939e33cfc6b28)
@@ -1299,21 +1299,26 @@
  * \param type Type of the Parameter to be read
  * \param value address of the value to be read (must have been allocated)
+ * \param repetition determines, if the keyword appears multiply in the config file, which repetition shall be parsed, i.e. 1 if not multiply
  * \param critical necessity of this keyword being specified (optional, critical)
  * \return 1 - found, 0 - not found
  * \sa ReadMainParameterFile Calling routine
  */
-int ParseForParameter(int verbose, FILE *file, const char *const name, int sequential, int const xth, int const yth, enum value_type type, void *value, enum necessity critical) {
+int ParseForParameter(int verbose, FILE *file, const char *const name, int sequential, int const xth, int const yth, enum value_type type, void *value, int repetition, enum necessity critical) {
 	int i,j;	// loop variables
 	int me, length, maxlength = -1;
 	long file_position = ftell(file);	// mark current position
-	char *dummy1, *dummy, *free_dummy;		// pointers in the line that is read in per step
+	char *dummy2, *dummy1, *dummy, *free_dummy;		// pointers in the line that is read in per step
+
+  if (repetition == 0)
+    Error(SomeError, "ParseForParameter(): argument repetition must not be 0!");
+
+  // allocated memory and rank in multi-process word
 	dummy1 = free_dummy = (char *) Malloc(256 * sizeof(char),"ParseForParameter: MemAlloc for dummy1 failed");
 	MPI_Comm_rank(MPI_COMM_WORLD, &me);
 
 	//fprintf(stderr,"(%i) Parsing for %s\n",me,name);	
-
 	int line = 0;	// marks line where parameter was found
 	int found = (type >= grid) ? 0 : (-yth + 1);	// marks if yth parameter name was found
-	while((found != 1)) {
+	while((found != repetition)) {
 		dummy1 = dummy = free_dummy;
 		do {
@@ -1348,8 +1353,8 @@
 		if (name != NULL) {
 			dummy = strchr(dummy1,'\t');	// set dummy on first tab or space which ever nearer
-      if (dummy == NULL) {
-        dummy = strchr(dummy1, ' ');  // if not found seek for space
-        while ((dummy != NULL) && ((*dummy == '\t') || (*dummy == ' ')))    // skip some more tabs and spaces if necessary
-          dummy++;
+      dummy2 = strchr(dummy1, ' ');  // if not found seek for space
+      if ((dummy != NULL) || (dummy2 != NULL)) {
+        if ((dummy == NULL) || ((dummy2 != NULL) && (dummy2 < dummy))) 
+          dummy = dummy2;
       }
 			if (dummy == NULL) {
@@ -1359,13 +1364,13 @@
 		    //Error(FileOpenParams, NULL);			
 			} else {
-				//fprintf(stderr,"(%i) found tab at %i\n",me,(char *)dummy-(char *)dummy1);
+				//fprintf(stderr,"(%i) found tab at %li\n",me,(long)((char *)dummy-(char *)dummy1));
 			}
 		} else dummy = dummy1;
 		// ... and check if it is the keyword!
-		if ((name == NULL) || ((dummy-dummy1 >= 3) && (strncmp(dummy1, name, strlen(name)) == 0))) {
+		if ((name == NULL) || ((dummy-dummy1 >= 3) && (strncmp(dummy1, name, strlen(name)) == 0) && (dummy-dummy1 == strlen(name)))) {
 			found++; // found the parameter!
 			//fprintf(stderr,"(%i) found %s at line %i\n",me, name, line);		
 			
-			if (found == 1) {
+			if (found == repetition) {
 				for (i=0;i<xth;i++) {	// i = rows
 					if (type >= grid) {
@@ -1413,18 +1418,30 @@
 							// otherwise we must skip all interjacent tabs and spaces and find next value
 							dummy1 = dummy;
-							dummy = strchr(dummy1, '\t');	// seek for tab or space
-              if (dummy == NULL) {
-                dummy = strchr(dummy1, ' ');  // if not found seek for space
-                while ((dummy != NULL) && ((*dummy == '\t') || (*dummy == ' ')))    // skip some more tabs and spaces if necessary
-                  dummy++;
+							dummy = strchr(dummy1, '\t');	// seek for tab or space (space if appears sooner)
+              dummy2 = strchr(dummy1, ' ');
+              //fprintf(stderr,"dummy (%p) is %c\t dummy2 (%p) is %c\n", dummy, (dummy == NULL ? '-' : *dummy), dummy2, (dummy2 == NULL ? '-' : *dummy2));
+              if ((dummy != NULL) || (dummy2 != NULL)) {
+                if ((dummy == NULL) || ((dummy2 != NULL) && (dummy2 < dummy))) 
+                  dummy = dummy2;
+                //while ((dummy != NULL) && ((*dummy == '\t') || (*dummy == ' ')))    // skip some more tabs and spaces if necessary
+                //  dummy++;
               }
 /*							while ((dummy != NULL) && ((*dummy == '\t') || (*dummy == ' ')))		// skip some more tabs and spaces if necessary
 								dummy++;*/
-							if (dummy == NULL) { // if still zero returned ...
+              //fprintf(stderr,"dummy is %c\n", *dummy);
+              dummy2 = strchr(dummy1, '#');
+              if ((dummy == NULL) || ((dummy2 != NULL) && (dummy1 >= dummy2))) { // if still zero returned or in comment area ...
 								dummy = strchr(dummy1, '\n');	 // ... at line end then
-								if ((j < yth-1) && (type < 4)) {	// check if xth value or not yet
-									if (verbose) fprintf(stderr,"(%i)Error: EoL at %i and still missing %i value(s) for parameter %s\n", me, line, yth-j, name);
-                  Free(free_dummy, "ParseForParameter: free_dummy");
-									return 0;
+								if ((i < xth-1) && (type < 4)) {	// check if xth value or not yet
+                  if (critical) {
+                    if (verbose) fprintf(stderr,"(%i)Error: EoL or comment at %i and still missing %i value(s) for parameter %s\n", me, line, yth-j, name);
+                    Free(free_dummy, "ParseForParameter: free_dummy");
+                    Error(InitReading, name);
+                  } else {
+                    //if (!sequential) 
+                    fseek(file, file_position, SEEK_SET);  // rewind to start position
+                    Free(free_dummy, "ParseForParameter: free_dummy");
+                    return 0;
+                  }
 							    //Error(FileOpenParams, NULL);			
 								}
@@ -1482,6 +1499,6 @@
 	}	
   if ((type >= row_int) && (verbose)) fprintf(stderr,"\n");
+	if (!sequential) fseek(file, file_position, SEEK_SET);	// rewind to start position
   Free(free_dummy, "ParseForParameter: free_dummy");
-	if (!sequential) fseek(file, file_position, SEEK_SET);	// rewind to start position
 	//fprintf(stderr, "(%i) End of Parsing\n\n",me);
 	
@@ -1503,5 +1520,5 @@
   FILE *file;
   int i, di, me;
-  double BField[NDIM_NDIM];
+  //double BField[NDIM];
   
   MPI_Comm_rank(MPI_COMM_WORLD, &me);
@@ -1516,13 +1533,13 @@
 
   P->Files.filename = MallocString(MAXDUMMYSTRING,"ReadParameters: filename");
-	ParseForParameter(P->Call.out[ReadOut],file, "mainname", 0, 1, 1, string_type, P->Files.filename, critical);
+	ParseForParameter(P->Call.out[ReadOut],file, "mainname", 0, 1, 1, string_type, P->Files.filename, 1, critical);
   //debug(P,"mainname");
   CreateMainname(P, filename);
   P->Files.default_path = MallocString(MAXDUMMYSTRING,"ReadParameters: default_path");
-  ParseForParameter(P->Call.out[ReadOut],file, "defaultpath", 0, 1, 1, string_type, P->Files.default_path, critical);
+  ParseForParameter(P->Call.out[ReadOut],file, "defaultpath", 0, 1, 1, string_type, P->Files.default_path, 1, critical);
   P->Files.pseudopot_path = MallocString(MAXDUMMYSTRING,"ReadParameters: pseudopot_path");
-  ParseForParameter(P->Call.out[ReadOut],file, "pseudopotpath", 0, 1, 1, string_type, P->Files.pseudopot_path, critical);
-  ParseForParameter(P->Call.out[ReadOut],file,"ProcPEGamma", 0, 1, 1, int_type, &(P->Par.proc[PEGamma]), critical);
-  ParseForParameter(P->Call.out[ReadOut],file,"ProcPEPsi", 0, 1, 1, int_type, &(P->Par.proc[PEPsi]), critical);
+  ParseForParameter(P->Call.out[ReadOut],file, "pseudopotpath", 0, 1, 1, string_type, P->Files.pseudopot_path, 1, critical);
+  ParseForParameter(P->Call.out[ReadOut],file,"ProcPEGamma", 0, 1, 1, int_type, &(P->Par.proc[PEGamma]), 1, critical);
+  ParseForParameter(P->Call.out[ReadOut],file,"ProcPEPsi", 0, 1, 1, int_type, &(P->Par.proc[PEPsi]), 1, critical);
   if (P->Call.proc[PEPsi]) { /* Kommandozeilenoption */
     P->Par.proc[PEPsi]  = P->Call.proc[PEPsi];
@@ -1532,12 +1549,12 @@
     Run
   */
-  if (!ParseForParameter(P->Call.out[ReadOut],file,"Seed", 0, 1, 1, int_type, &(R->Seed), optional))
+  if (!ParseForParameter(P->Call.out[ReadOut],file,"Seed", 0, 1, 1, int_type, &(R->Seed), 1, optional))
     R->Seed = 1;
-  if (!ParseForParameter(P->Call.out[ReadOut],file,"WaveNr", 0, 1, 1, int_type, &(R->WaveNr), optional))
+  if (!ParseForParameter(P->Call.out[ReadOut],file,"WaveNr", 0, 1, 1, int_type, &(R->WaveNr), 1, optional))
     R->WaveNr = 0;
-  if (!ParseForParameter(P->Call.out[ReadOut],file,"Lambda", 0, 1, 1, double_type, &R->Lambda, optional))
+  if (!ParseForParameter(P->Call.out[ReadOut],file,"Lambda", 0, 1, 1, double_type, &R->Lambda, 1, optional))
     R->Lambda = 1.;
 
-  if(!ParseForParameter(P->Call.out[ReadOut],file,"DoOutOrbitals", 0, 1, 1, int_type, &F->DoOutOrbitals, optional)) {
+  if(!ParseForParameter(P->Call.out[ReadOut],file,"DoOutOrbitals", 0, 1, 1, int_type, &F->DoOutOrbitals, 1, optional)) {
     F->DoOutOrbitals = 0;
   } else {
@@ -1545,8 +1562,8 @@
     if (F->DoOutOrbitals > 1) F->DoOutOrbitals = 1;
   }
-  ParseForParameter(P->Call.out[ReadOut],file,"DoOutVis", 0, 1, 1, int_type, &F->DoOutVis, critical);
+  ParseForParameter(P->Call.out[ReadOut],file,"DoOutVis", 0, 1, 1, int_type, &F->DoOutVis, 1, critical);
   if (F->DoOutVis < 0) F->DoOutVis = 0;
-  if (F->DoOutVis > 1) F->DoOutVis = 1;
-  if (!ParseForParameter(P->Call.out[ReadOut],file,"VectorPlane", 0, 1, 1, int_type, &R->VectorPlane, optional))
+  if (F->DoOutVis > 2) F->DoOutVis = 1;
+  if (!ParseForParameter(P->Call.out[ReadOut],file,"VectorPlane", 0, 1, 1, int_type, &R->VectorPlane, 1, optional))
     R->VectorPlane = -1;
   if (R->VectorPlane < -1 || R->VectorPlane > 2) {
@@ -1554,18 +1571,20 @@
     R->VectorPlane = -1;
   }
-  if (!ParseForParameter(P->Call.out[ReadOut],file,"VectorCut", 0, 1, 1, double_type, &R->VectorCut, optional))
+  if (!ParseForParameter(P->Call.out[ReadOut],file,"VectorCut", 0, 1, 1, double_type, &R->VectorCut, 1, optional))
     R->VectorCut = 0.;
-  ParseForParameter(P->Call.out[ReadOut],file,"DoOutMes", 0, 1, 1, int_type, &F->DoOutMes, critical);
+  ParseForParameter(P->Call.out[ReadOut],file,"DoOutMes", 0, 1, 1, int_type, &F->DoOutMes, 1, critical);
   if (F->DoOutMes < 0) F->DoOutMes = 0;
   if (F->DoOutMes > 1) F->DoOutMes = 1;
-  if (!ParseForParameter(P->Call.out[ReadOut],file,"DoOutCurr", 0, 1, 1, int_type, &F->DoOutCurr, optional))
+  if (!ParseForParameter(P->Call.out[ReadOut],file,"DoOutCurr", 0, 1, 1, int_type, &F->DoOutCurr, 1, optional))
     F->DoOutCurr = 0;
+  if (!ParseForParameter(P->Call.out[ReadOut],file,"DoOutNICS", 0, 1, 1, int_type, &F->DoOutNICS, 1, optional))
+    F->DoOutNICS = 0;
   if (F->DoOutCurr < 0) F->DoOutCurr = 0;
   if (F->DoOutCurr > 1) F->DoOutCurr = 1; 
-  ParseForParameter(P->Call.out[ReadOut],file,"AddGramSch", 0, 1, 1, int_type, &R->UseAddGramSch, critical);
+  ParseForParameter(P->Call.out[ReadOut],file,"AddGramSch", 0, 1, 1, int_type, &R->UseAddGramSch, 1, critical);
   if (R->UseAddGramSch < 0) R->UseAddGramSch = 0;
   if (R->UseAddGramSch > 2) R->UseAddGramSch = 2;
   if(P->Call.out[ReadOut]) fprintf(stderr,"(%i)UseAddGramSch = %i\n",me,R->UseAddGramSch);
-  if(!ParseForParameter(P->Call.out[ReadOut],file,"CommonWannier", 0, 1, 1, int_type, &R->CommonWannier, optional)) {
+  if(!ParseForParameter(P->Call.out[ReadOut],file,"CommonWannier", 0, 1, 1, int_type, &R->CommonWannier, 1, optional)) {
     R->CommonWannier = 0;
   } else {
@@ -1573,5 +1592,5 @@
     if (R->CommonWannier > 4) R->CommonWannier = 4;
   }
-  if(!ParseForParameter(P->Call.out[ReadOut],file,"SawtoothStart", 0, 1, 1, double_type, &P->Lat.SawtoothStart, optional)) {
+  if(!ParseForParameter(P->Call.out[ReadOut],file,"SawtoothStart", 0, 1, 1, double_type, &P->Lat.SawtoothStart, 1, optional)) {
     P->Lat.SawtoothStart = 0.01;
   } else {
@@ -1579,5 +1598,5 @@
     if (P->Lat.SawtoothStart > 1.) P->Lat.SawtoothStart = 1.;
   }
-  if (!ParseForParameter(P->Call.out[ReadOut],file,"DoBrent", 0, 1, 1, int_type, &R->DoBrent, optional)) {
+  if (!ParseForParameter(P->Call.out[ReadOut],file,"DoBrent", 0, 1, 1, int_type, &R->DoBrent, 1, optional)) {
     R->DoBrent = 0;
   } else {
@@ -1587,18 +1606,24 @@
    
   
-	ParseForParameter(P->Call.out[ReadOut],file,"MaxOuterStep", 0, 1, 1, int_type, &R->MaxOuterStep, critical);
+  if (!ParseForParameter(P->Call.out[ReadOut],file,"MaxOuterStep", 0, 1, 1, int_type, &R->MaxOuterStep, 1, optional))
+    R->MaxOuterStep = 0;
+  if (!ParseForParameter(P->Call.out[ReadOut],file,"MaxStructOptStep", 0, 1, 1, int_type, &R->MaxStructOptStep, 1, optional)) {
+    if (R->MaxOuterStep != 0) { // if StructOptStep not explicitely specified, then use OuterStep
+      R->MaxStructOptStep = R->MaxOuterStep;
+    } else {
+      R->MaxStructOptStep = 100;
+    }
+  }
   // the following values are MD specific and should be dropped in further revisions
-	if (!ParseForParameter(P->Call.out[ReadOut],file,"Deltat", 0, 1, 1, double_type, &R->delta_t, optional))
-    R->delta_t = 0;
-	if (!ParseForParameter(P->Call.out[ReadOut],file,"OutVisStep", 0, 1, 1, int_type, &R->OutVisStep, optional))
+	if (!ParseForParameter(P->Call.out[ReadOut],file,"Deltat", 0, 1, 1, double_type, &R->delta_t, 1, optional))
+    R->delta_t = 1.;
+	if (!ParseForParameter(P->Call.out[ReadOut],file,"OutVisStep", 0, 1, 1, int_type, &R->OutVisStep, 1, optional))
     R->OutVisStep = 10;
-	if (!ParseForParameter(P->Call.out[ReadOut],file,"OutSrcStep", 0, 1, 1, int_type, &R->OutSrcStep, optional))
+	if (!ParseForParameter(P->Call.out[ReadOut],file,"OutSrcStep", 0, 1, 1, int_type, &R->OutSrcStep, 1, optional))
     R->OutSrcStep = 5;
-	if (!ParseForParameter(P->Call.out[ReadOut],file,"TargetTemp", 0, 1, 1, double_type, &R->TargetTemp, optional))
+	if (!ParseForParameter(P->Call.out[ReadOut],file,"TargetTemp", 0, 1, 1, double_type, &R->TargetTemp, 1, optional))
     R->TargetTemp = 0;
-	if (!ParseForParameter(P->Call.out[ReadOut],file,"ScaleTempStep", 0, 1, 1, int_type, &R->ScaleTempStep, optional))
-    R->ScaleTempStep = 25;
     
-  if (!ParseForParameter(P->Call.out[ReadOut],file,"EpsWannier", 0, 1, 1, double_type, &R->EpsWannier, optional))
+  if (!ParseForParameter(P->Call.out[ReadOut],file,"EpsWannier", 0, 1, 1, double_type, &R->EpsWannier, 1, optional))
     R->EpsWannier = 1e-8;
   
@@ -1607,14 +1632,14 @@
   //if (R->MaxOuterStep <= 0) R->MaxOuterStep = 1;
   if(P->Call.out[NormalOut]) fprintf(stderr,"(%i)MaxOuterStep = %i\n",me,R->MaxOuterStep);
-  if(P->Call.out[NormalOut]) fprintf(stderr,"(%i)t = %g\n",me,R->t);
-  ParseForParameter(P->Call.out[ReadOut],file,"MaxPsiStep", 0, 1, 1, int_type, &R->MaxPsiStep, critical);
+  if(P->Call.out[NormalOut]) fprintf(stderr,"(%i)Deltat = %g\n",me,R->delta_t);
+  ParseForParameter(P->Call.out[ReadOut],file,"MaxPsiStep", 0, 1, 1, int_type, &R->MaxPsiStep, 1, critical);
   if (R->MaxPsiStep <= 0) R->MaxPsiStep = 3;
   if(P->Call.out[NormalOut]) fprintf(stderr,"(%i)MaxPsiStep = %i\n",me,R->MaxPsiStep);
   
-  ParseForParameter(P->Call.out[ReadOut],file,"MaxMinStep", 0, 1, 1, int_type, &R->MaxMinStep, critical);
-	ParseForParameter(P->Call.out[ReadOut],file,"RelEpsTotalE", 0, 1, 1, double_type, &R->RelEpsTotalEnergy, critical);
-	ParseForParameter(P->Call.out[ReadOut],file,"RelEpsKineticE", 0, 1, 1, double_type, &R->RelEpsKineticEnergy, critical);
-  ParseForParameter(P->Call.out[ReadOut],file,"MaxMinStopStep", 0, 1, 1, int_type, &R->MaxMinStopStep, critical);
-  ParseForParameter(P->Call.out[ReadOut],file,"MaxMinGapStopStep", 0, 1, 1, int_type, &R->MaxMinGapStopStep, critical);
+  ParseForParameter(P->Call.out[ReadOut],file,"MaxMinStep", 0, 1, 1, int_type, &R->MaxMinStep, 1, critical);
+	ParseForParameter(P->Call.out[ReadOut],file,"RelEpsTotalE", 0, 1, 1, double_type, &R->RelEpsTotalEnergy, 1, critical);
+	ParseForParameter(P->Call.out[ReadOut],file,"RelEpsKineticE", 0, 1, 1, double_type, &R->RelEpsKineticEnergy, 1, critical);
+  ParseForParameter(P->Call.out[ReadOut],file,"MaxMinStopStep", 0, 1, 1, int_type, &R->MaxMinStopStep, 1, critical);
+  ParseForParameter(P->Call.out[ReadOut],file,"MaxMinGapStopStep", 0, 1, 1, int_type, &R->MaxMinGapStopStep, 1, critical);
   if (R->MaxMinStep <= 0) R->MaxMinStep = R->MaxPsiStep;
   if (R->MaxMinStopStep < 1) R->MaxMinStopStep = 1;
@@ -1622,9 +1647,9 @@
   if(P->Call.out[NormalOut]) fprintf(stderr,"(%i)MaxMinStep = %i\tRelEpsTotalEnergy: %e\tRelEpsKineticEnergy %e\tMaxMinStopStep %i\n",me,R->MaxMinStep,R->RelEpsTotalEnergy,R->RelEpsKineticEnergy,R->MaxMinStopStep);
   
-  ParseForParameter(P->Call.out[ReadOut],file,"MaxInitMinStep", 0, 1, 1, int_type, &R->MaxInitMinStep, critical);
-	ParseForParameter(P->Call.out[ReadOut],file,"InitRelEpsTotalE", 0, 1, 1, double_type, &R->InitRelEpsTotalEnergy, critical);
-	ParseForParameter(P->Call.out[ReadOut],file,"InitRelEpsKineticE", 0, 1, 1, double_type, &R->InitRelEpsKineticEnergy, critical);
-  ParseForParameter(P->Call.out[ReadOut],file,"InitMaxMinStopStep", 0, 1, 1, int_type, &R->InitMaxMinStopStep, critical);
-  ParseForParameter(P->Call.out[ReadOut],file,"InitMaxMinGapStopStep", 0, 1, 1, int_type, &R->InitMaxMinGapStopStep, critical);
+  ParseForParameter(P->Call.out[ReadOut],file,"MaxInitMinStep", 0, 1, 1, int_type, &R->MaxInitMinStep, 1, critical);
+	ParseForParameter(P->Call.out[ReadOut],file,"InitRelEpsTotalE", 0, 1, 1, double_type, &R->InitRelEpsTotalEnergy, 1, critical);
+	ParseForParameter(P->Call.out[ReadOut],file,"InitRelEpsKineticE", 0, 1, 1, double_type, &R->InitRelEpsKineticEnergy, 1, critical);
+  ParseForParameter(P->Call.out[ReadOut],file,"InitMaxMinStopStep", 0, 1, 1, int_type, &R->InitMaxMinStopStep, 1, critical);
+  ParseForParameter(P->Call.out[ReadOut],file,"InitMaxMinGapStopStep", 0, 1, 1, int_type, &R->InitMaxMinGapStopStep, 1, critical);
   if (R->MaxInitMinStep <= 0) R->MaxInitMinStep = R->MaxPsiStep;
   if (R->InitMaxMinStopStep < 1) R->InitMaxMinStopStep = 1;
@@ -1633,12 +1658,13 @@
   
   // Unit cell and magnetic field
-  ParseForParameter(P->Call.out[ReadOut],file, "BoxLength", 0, 3, 3, lower_trigrid, &P->Lat.RealBasis, critical); /* Lattice->RealBasis */
-  ParseForParameter(P->Call.out[ReadOut],file, "DoPerturbation", 0, 1, 1, int_type, &R->DoPerturbation, optional);
+  ParseForParameter(P->Call.out[ReadOut],file, "BoxLength", 0, 3, 3, lower_trigrid, &P->Lat.RealBasis, 1, critical); /* Lattice->RealBasis */
+  ParseForParameter(P->Call.out[ReadOut],file, "DoPerturbation", 0, 1, 1, int_type, &R->DoPerturbation, 1, optional);
   if (!R->DoPerturbation) {
-    if (!ParseForParameter(P->Call.out[ReadOut],file, "BField", 0, 1, 3, grid, &BField, optional)) { // old parameter for perturbation with field direction
+    if (!ParseForParameter(P->Call.out[ReadOut],file, "BField", 0, 1, 3, grid, &R->BField, 1, optional)) { // old parameter for perturbation with field direction
       if (P->Call.out[ReadOut]) fprintf(stderr,"(%i) No perturbation specified.\n", P->Par.me);
-      if (F->DoOutCurr) {
-        if (P->Call.out[ReadOut]) fprintf(stderr,"(%i) DoOutCurr=  1 though no DoPerturbation specified: setting DoOutCurr = 0\n", P->Par.me);
+      if (F->DoOutCurr || F->DoOutNICS) {
+        if (P->Call.out[ReadOut]) fprintf(stderr,"(%i) DoOutCurr/DoOutNICS = 1 though no DoPerturbation specified: setting DoOutCurr = 0\n", P->Par.me);
         F->DoOutCurr = 0;
+        F->DoOutNICS = 0;
         R->DoPerturbation = 0;
       }
@@ -1653,5 +1679,5 @@
   RMatReci3(P->Lat.InvBasis,  P->Lat.RealBasis);
 
-  if (!ParseForParameter(P->Call.out[ReadOut],file,"DoFullCurrent", 0, 1, 1, int_type, &R->DoFullCurrent, optional))
+  if (!ParseForParameter(P->Call.out[ReadOut],file,"DoFullCurrent", 0, 1, 1, int_type, &R->DoFullCurrent, 1, optional))
     R->DoFullCurrent = 0;
   if (R->DoFullCurrent < 0) R->DoFullCurrent = 0;
@@ -1662,13 +1688,13 @@
   }
 
-	ParseForParameter(P->Call.out[ReadOut],file,"ECut", 0, 1, 1, double_type, &P->Lat.ECut, critical);
+	ParseForParameter(P->Call.out[ReadOut],file,"ECut", 0, 1, 1, double_type, &P->Lat.ECut, 1, critical);
   if(P->Call.out[NormalOut]) fprintf(stderr,"(%i) ECut = %e -> %e Rydberg\n", me, P->Lat.ECut, 2.*P->Lat.ECut);
-  ParseForParameter(P->Call.out[ReadOut],file,"MaxLevel", 0, 1, 1, int_type, &P->Lat.MaxLevel, critical);
-  ParseForParameter(P->Call.out[ReadOut],file,"Level0Factor", 0, 1, 1, int_type, &P->Lat.Lev0Factor, critical);
+  ParseForParameter(P->Call.out[ReadOut],file,"MaxLevel", 0, 1, 1, int_type, &P->Lat.MaxLevel, 1, critical);
+  ParseForParameter(P->Call.out[ReadOut],file,"Level0Factor", 0, 1, 1, int_type, &P->Lat.Lev0Factor, 1, critical);
   if (P->Lat.Lev0Factor < 2) {
     if(P->Call.out[ReadOut]) fprintf(stderr,"(%i) Set Lev0Factor to 2!\n",me);
     P->Lat.Lev0Factor = 2;
   }
-  ParseForParameter(P->Call.out[ReadOut],file,"RiemannTensor", 0, 1, 1, int_type, &di, critical);
+  ParseForParameter(P->Call.out[ReadOut],file,"RiemannTensor", 0, 1, 1, int_type, &di, 1, critical);
   if (di >= 0 && di < 2) {
     P->Lat.RT.Use = (enum UseRiemannTensor)di;
@@ -1691,5 +1717,5 @@
 	      P->Lat.MaxLevel = 3;
 	    }
-	    ParseForParameter(P->Call.out[ReadOut],file,"RiemannLevel", 0, 1, 1, int_type, &P->Lat.RT.RiemannLevel, critical);
+	    ParseForParameter(P->Call.out[ReadOut],file,"RiemannLevel", 0, 1, 1, int_type, &P->Lat.RT.RiemannLevel, 1, critical);
 		  if(P->Call.out[ReadOut]) fprintf(stderr,"(%i)!RiemannLevel = %i\n",me,P->Lat.RT.RiemannLevel);
 	    if (P->Lat.RT.RiemannLevel < 2) {
@@ -1701,5 +1727,5 @@
 	      P->Lat.RT.RiemannLevel = P->Lat.MaxLevel-1;
 	    }
-	    ParseForParameter(P->Call.out[ReadOut],file,"LevRFactor", 0, 1, 1, int_type, &P->Lat.LevRFactor, critical);
+	    ParseForParameter(P->Call.out[ReadOut],file,"LevRFactor", 0, 1, 1, int_type, &P->Lat.LevRFactor, 1, critical);
 	    if(P->Call.out[ReadOut]) fprintf(stderr,"(%i)!LevRFactor = %i\n",me,P->Lat.LevRFactor);
 	    if (P->Lat.LevRFactor < 2) {
@@ -1711,5 +1737,5 @@
 	    break;
   }
-  ParseForParameter(P->Call.out[ReadOut],file,"PsiType", 0, 1, 1, int_type, &di, critical);
+  ParseForParameter(P->Call.out[ReadOut],file,"PsiType", 0, 1, 1, int_type, &di, 1, critical);
   if (di >= 0 && di < 2) {
     P->Lat.Psi.Use = (enum UseSpinType)di;
@@ -1722,6 +1748,6 @@
   switch (P->Lat.Psi.Use) {
   case UseSpinDouble:
-  	ParseForParameter(P->Call.out[ReadOut],file,"MaxPsiDouble", 0, 1, 1, int_type, &P->Lat.Psi.GlobalNo[PsiMaxNoDouble], critical);
-    if (!ParseForParameter(P->Call.out[ReadOut],file,"AddPsis", 0, 1, 1, int_type, &P->Lat.Psi.GlobalNo[PsiMaxAdd], optional))
+  	ParseForParameter(P->Call.out[ReadOut],file,"MaxPsiDouble", 0, 1, 1, int_type, &P->Lat.Psi.GlobalNo[PsiMaxNoDouble], 1, critical);
+    if (!ParseForParameter(P->Call.out[ReadOut],file,"AddPsis", 0, 1, 1, int_type, &P->Lat.Psi.GlobalNo[PsiMaxAdd], 1, optional))
       P->Lat.Psi.GlobalNo[PsiMaxAdd] = 0;
     if (P->Lat.Psi.GlobalNo[PsiMaxAdd] != 0)
@@ -1735,7 +1761,7 @@
   case UseSpinUpDown:
     if (P->Par.proc[PEPsi] % 2) Error(SomeError,"UseSpinUpDown & P->Par.proc[PEGamma] % 2");
-    ParseForParameter(P->Call.out[ReadOut],file,"PsiMaxNoUp", 0, 1, 1, int_type, &P->Lat.Psi.GlobalNo[PsiMaxNoUp], critical);
-    ParseForParameter(P->Call.out[ReadOut],file,"PsiMaxNoDown", 0, 1, 1, int_type, &P->Lat.Psi.GlobalNo[PsiMaxNoDown], critical);
-    if (!ParseForParameter(P->Call.out[ReadOut],file,"AddPsis", 0, 1, 1, int_type, &P->Lat.Psi.GlobalNo[PsiMaxAdd], optional))
+    ParseForParameter(P->Call.out[ReadOut],file,"PsiMaxNoUp", 0, 1, 1, int_type, &P->Lat.Psi.GlobalNo[PsiMaxNoUp], 1, critical);
+    ParseForParameter(P->Call.out[ReadOut],file,"PsiMaxNoDown", 0, 1, 1, int_type, &P->Lat.Psi.GlobalNo[PsiMaxNoDown], 1, critical);
+    if (!ParseForParameter(P->Call.out[ReadOut],file,"AddPsis", 0, 1, 1, int_type, &P->Lat.Psi.GlobalNo[PsiMaxAdd], 1, optional))
       P->Lat.Psi.GlobalNo[PsiMaxAdd] = 0;
     if (P->Lat.Psi.GlobalNo[PsiMaxAdd] != 0)
Index: pcp/src/init.h
===================================================================
--- pcp/src/init.h	(revision e00f4732f333881727b6265e9f995246359ea74f)
+++ pcp/src/init.h	(revision 961b75c543d6081b2758726da4c939e33cfc6b28)
@@ -29,5 +29,5 @@
 
 /* Einlesen der Parameterdatei */
-int ParseForParameter(int verbose, FILE *file, const char *const name, int sequential, int const xth, int const yth, enum value_type type, void *value, enum necessity critical);
+int ParseForParameter(int verbose, FILE *file, const char *const name, int sequential, int const xth, int const yth, enum value_type type, void *value, int repetition, enum necessity critical);
 void ReadParameters(struct Problem *const  P, const char *const filename);
 
Index: pcp/src/ions.c
===================================================================
--- pcp/src/ions.c	(revision e00f4732f333881727b6265e9f995246359ea74f)
+++ pcp/src/ions.c	(revision 961b75c543d6081b2758726da4c939e33cfc6b28)
@@ -29,4 +29,5 @@
 #include <string.h>
 #include <unistd.h>
+#include <gsl/gsl_randist.h>
 #include "mymath.h"
 
@@ -50,4 +51,98 @@
 #include "init.h"
 
+/** Parses for a given keyword the ion coordinates and velocites.
+ * \param *P Problem at hand, contains pointer to Ion and IonType structure
+ * \param *source
+ * \param *keyword keyword string
+ * \param repetition which repeated version of the keyword shall be read by ReadParameters()
+ * \param relative whether atom coordinates are relative to unit cell size or absolute
+ * \param Bohr conversion factor to atomiclength (e.g. 1.0 if coordinates in bohr radius, 0.529... if given in Angstrom)
+ * \param sigma variance of maxwell-boltzmann distribution
+ * \param *R where to put the read coordinates
+ * \param *U where to put the read velocities
+ * \param *IMT whether its MoveType#FixedIon or not
+ * \return if 1 - success, 0 - failure
+ */
+int ParseIonsCoordinatesAndVelocities(struct Problem *P, FILE *source, char *keyword, int repetition, int relative, double sigma, double *R, double *U, int *IMT) 
+{
+  //struct RunStruct *Run = &P->R;
+  struct Ions *I = &P->Ion;
+  struct Lattice *L = &P->Lat;
+  int k;
+  double R_trafo[3];
+  gsl_rng * r;
+  const gsl_rng_type * T;
+  double Bohr = (I->IsAngstroem) ? ANGSTROEMINBOHRRADIUS : 1.;
+  
+  gsl_rng_env_setup();
+  T = gsl_rng_default;
+  r = gsl_rng_alloc (T);
+  if (ParseForParameter(P->Call.out[ReadOut],source, keyword, 0, 1, 1, double_type, &R_trafo[0], repetition, optional) &&
+    ParseForParameter(P->Call.out[ReadOut],source, keyword, 0, 2, 1, double_type, &R_trafo[1], repetition, optional) &&
+    ParseForParameter(P->Call.out[ReadOut],source, keyword, 0, 3, 1, double_type, &R_trafo[2], repetition, optional) &&
+    ParseForParameter(P->Call.out[ReadOut],source, keyword, 0, 4, 1, int_type, IMT, repetition, optional)) {
+    // if read successful, then try also parsing velocities
+    if ((ParseForParameter(P->Call.out[ReadOut],source, keyword, 0, 5, 1, double_type, &U[0], repetition, optional)) &&
+      (ParseForParameter(P->Call.out[ReadOut],source, keyword, 0, 6, 1, double_type, &U[1], repetition, optional)) && 
+      (ParseForParameter(P->Call.out[ReadOut],source, keyword, 0, 7, 1, double_type, &U[2], repetition, optional))) {
+      //nothing
+    } else {  // if no velocities specified, set to maxwell-boltzmann distributed
+      if ((I->Thermostat != None)) {
+        if (P->Par.me == 0) fprintf(stderr, "(%i) Missing velocities of ion %s: Maxwell-Boltzmann with %lg\n", P->Par.me, keyword, sigma);
+        for (k=0;k<NDIM;k++)
+          U[k] = gsl_ran_gaussian (r, sigma);
+      } else {
+        for (k=0;k<NDIM;k++)
+          U[k] = 0.;
+      }
+    }
+    // change if coordinates were relative
+    if (relative) {
+      //fprintf(stderr,"(%i)Ion coordinates are relative %i ... \n",P->Par.me, relative);
+      RMat33Vec3(R, L->RealBasis, R_trafo); // multiply with real basis
+    } else {
+      for (k=0;k<NDIM;k++) // otherweise just copy to vector
+        R[k] = R_trafo[k];
+    }
+    // check if given MoveType is valid
+    if (*IMT < 0 || *IMT >= MaxIonMoveType) {
+      fprintf(stderr,"Bad Ion MoveType set to MoveIon for %s\n", keyword);
+      *IMT = 0;
+    }
+    if (!relative) {
+      //fprintf(stderr,"converting with %lg\n", Bohr);
+      SM(R, Bohr, NDIM); // convert angstroem to bohrradius
+      SM(U, Bohr, NDIM); // convert angstroem to bohrradius
+    }
+    // finally periodicly correct coordinates to remain in unit cell
+    RMat33Vec3(R_trafo, L->InvBasis, R);
+    for (k=0; k <NDIM; k++) { // periodic correction until within unit cell
+      while (R_trafo[k] < 0) 
+        R_trafo[k] += 1.0;
+      while (R_trafo[k] >= 1.0) 
+        R_trafo[k] -= 1.0;
+    }
+    RMat33Vec3(R, L->RealBasis, R_trafo);
+    // print coordinates if desired
+    if ((P->Call.out[ReadOut])) {
+      fprintf(stderr,"(%i) coordinates of Ion %s: (x,y,z) = (",P->Par.me, keyword);
+      for(k=0;k<NDIM;k++) {
+        fprintf(stderr,"%lg ",R[k]);
+        if (k != NDIM-1) fprintf(stderr,", ");
+        else fprintf(stderr,")\n");
+      }
+    }
+    gsl_rng_free (r);
+    return 1;
+  } else {
+    for (k=0;k<NDIM;k++) {
+      U[k] = R[k] = 0.;
+    }
+    *IMT = 0; 
+    gsl_rng_free (r);
+    return 0;
+  }
+}
+
 /** Readin of the ion section of parameter file.
  * Among others the following paramaters are read from file:
@@ -60,29 +155,30 @@
 	struct Ions *I = &P->Ion;
   //struct RunStruct *Run = &P->R;
-  int i,it,j,IMT,BorAng, d,me, count;
-  struct Lattice *L = &P->Lat;
-  int k;
+  int i,it,j,IMT,d,me, count;
   int relative; // whether read-in coordinate are relative (1) or absolute  
   double Bohr = ANGSTROEMINBOHRRADIUS; /* Angstroem */
   double sigma, R[3], U[3];
-
+  struct IonConstrained *ptr = NULL;
+     
   I->Max_TotalIons = 0;
   I->TotalZval = 0;
   MPI_Comm_rank(MPI_COMM_WORLD, &me);
-  ParseForParameter(P->Call.out[ReadOut],source,"RCut", 0, 1, 1, double_type, &I->R_cut, critical);
-  ParseForParameter(P->Call.out[ReadOut],source,"IsAngstroem", 0, 1, 1, int_type, &BorAng, critical);
+  ParseForParameter(P->Call.out[ReadOut],source,"RCut", 0, 1, 1, double_type, &I->R_cut, 1, critical);
+  ParseForParameter(P->Call.out[ReadOut],source,"IsAngstroem", 0, 1, 1, int_type, &I->IsAngstroem, 1, critical);
   if (!I->IsAngstroem) {
-  	Bohr = 1.0;
-	} else {  // adapt RealBasis (is in Angstroem as well)
+    Bohr = 1.0;
+  } else {  // adapt RealBasis (is in Angstroem as well)
     SM(P->Lat.RealBasis, Bohr, NDIM_NDIM); 
     RMatReci3(P->Lat.InvBasis,  P->Lat.RealBasis);
   }
-  ParseForParameter(P->Call.out[ReadOut],source,"MaxTypes", 0, 1, 1, int_type, &I->Max_Types, critical);
+  ParseForParameter(P->Call.out[ReadOut],source,"MaxTypes", 0, 1, 1, int_type, &I->Max_Types, 1, critical);
   I->I = (struct IonType *) Malloc(sizeof(struct IonType)*I->Max_Types, "IonsInitRead: IonType");
-  if (!ParseForParameter(P->Call.out[ReadOut],source,"RelativeCoord", 0, 1, 1, int_type, &relative, optional))
+  if (!ParseForParameter(P->Call.out[ReadOut],source,"RelativeCoord", 0, 1, 1, int_type, &relative, 1, optional))
     relative = 0;
-  if (!ParseForParameter(P->Call.out[ReadOut],source,"StructOpt", 0, 1, 1, int_type, &I->StructOpt, optional))
+  if (!ParseForParameter(P->Call.out[ReadOut],source,"StructOpt", 0, 1, 1, int_type, &I->StructOpt, 1, optional))
     I->StructOpt = 0;
-
+  
+  //InitThermostats(P, source);
+  
   /* Ions Data */
   I->RLatticeVec = NULL;
@@ -95,15 +191,16 @@
     I->I[i].Name = MallocString(255, "IonsInitRead: Name");
     I->I[i].Symbol = MallocString(255, "IonsInitRead: Symbol");
-		ParseForParameter(P->Call.out[ReadOut],source, name, 0, 1, 1, int_type, &I->I[i].Max_IonsOfType, critical);
-		ParseForParameter(P->Call.out[ReadOut],source, name, 0, 2, 1, int_type, &I->I[i].Z, critical);
-		ParseForParameter(P->Call.out[ReadOut],source, name, 0, 3, 1, double_type, &I->I[i].rgauss, critical);
-		ParseForParameter(P->Call.out[ReadOut],source, name, 0, 4, 1, int_type, &I->I[i].l_max, critical);
-		ParseForParameter(P->Call.out[ReadOut],source, name, 0, 5, 1, int_type, &I->I[i].l_loc, critical);
-		ParseForParameter(P->Call.out[ReadOut],source, name, 0, 6, 1, double_type, &I->I[i].IonMass, critical);
-    if(!ParseForParameter(P->Call.out[ReadOut],source, name, 0, 7, 1, string_type, &I->I[i].Name[0], optional))
+		ParseForParameter(P->Call.out[ReadOut],source, name, 0, 1, 1, int_type, &I->I[i].Max_IonsOfType, 1, critical);
+		ParseForParameter(P->Call.out[ReadOut],source, name, 0, 2, 1, int_type, &I->I[i].Z, 1, critical);
+		ParseForParameter(P->Call.out[ReadOut],source, name, 0, 3, 1, double_type, &I->I[i].rgauss, 1, critical);
+		ParseForParameter(P->Call.out[ReadOut],source, name, 0, 4, 1, int_type, &I->I[i].l_max, 1, critical);
+		ParseForParameter(P->Call.out[ReadOut],source, name, 0, 5, 1, int_type, &I->I[i].l_loc, 1, critical);
+		ParseForParameter(P->Call.out[ReadOut],source, name, 0, 6, 1, double_type, &I->I[i].IonMass, 1, critical);
+    if(!ParseForParameter(P->Call.out[ReadOut],source, name, 0, 7, 1, string_type, &I->I[i].Name[0], 1, optional))
       sprintf(&I->I[i].Name[0],"type%i",i);
-    if(!ParseForParameter(P->Call.out[ReadOut],source, name, 0, 8, 1, string_type, &I->I[i].Symbol[0], optional))
+    if(!ParseForParameter(P->Call.out[ReadOut],source, name, 0, 8, 1, string_type, &I->I[i].Symbol[0], 1, optional))
       sprintf(&I->I[i].Symbol[0],"t%i",i);
     if (P->Call.out[ReadOut]) fprintf(stderr,"(%i) Element name: %s, symbol: %s\n",P->Par.me, I->I[i].Name, I->I[i].Symbol);
+    I->I[i].moment = (double **) Malloc(sizeof(double *) * I->I[i].Max_IonsOfType, "IonsInitRead: moment");
     I->I[i].sigma = (double **) Malloc(sizeof(double *) * I->I[i].Max_IonsOfType, "IonsInitRead: sigma");
     I->I[i].sigma_rezi = (double **) Malloc(sizeof(double *) * I->I[i].Max_IonsOfType, "IonsInitRead: sigma_rezi");
@@ -111,4 +208,5 @@
     I->I[i].sigma_rezi_PAS = (double **) Malloc(sizeof(double *) * I->I[i].Max_IonsOfType, "IonsInitRead: sigma_rezi_PAS");
     for (it=0;it<I->I[i].Max_IonsOfType;it++) {
+      I->I[i].moment[it] = (double *) Malloc(sizeof(double) * NDIM*NDIM, "IonsInitRead: moment");
       I->I[i].sigma[it] = (double *) Malloc(sizeof(double) * NDIM*NDIM, "IonsInitRead: sigma");
       I->I[i].sigma_rezi[it] = (double *) Malloc(sizeof(double) * NDIM*NDIM, "IonsInitRead: sigma_rezi");
@@ -117,6 +215,9 @@
     }
     //I->I[i].IonFac = I->I[i].IonMass;
-    I->I[i].IonFac = 1/I->I[i].IonMass;
+    I->I[i].IonMass *= Units2Electronmass;  // in config given in amu not in "atomic units" (electron mass, m_e = 1)
+                                            // proportional factor between electron and proton mass
+    I->I[i].IonFac = Units2Electronmass/I->I[i].IonMass;
     I->TotalMass += I->I[i].Max_IonsOfType*I->I[i].IonMass;
+    sigma = sqrt(P->R.TargetTemp/I->I[i].IonMass);
     I->I[i].ZFactor = -I->I[i].Z*4.*PI;
     I->I[i].alpha = (double *) Malloc(sizeof(double) * I->I[i].Max_IonsOfType, "IonsInitRead: alpha");
@@ -124,4 +225,7 @@
     I->I[i].R_old = (double *) Malloc(sizeof(double) *NDIM* I->I[i].Max_IonsOfType, "IonsInitRead: R_old");
     I->I[i].R_old_old = (double *) Malloc(sizeof(double) *NDIM* I->I[i].Max_IonsOfType, "IonsInitRead: R_old_old");
+    I->I[i].Constraints = (struct IonConstrained **) Malloc(sizeof(struct IonConstrained *)*I->I[i].Max_IonsOfType, "IonsInitRead: **Constraints");
+    for (it=0;it<I->I[i].Max_IonsOfType;it++)
+      I->I[i].Constraints[it] = NULL;
     I->I[i].FIon = (double *) Malloc(sizeof(double) * NDIM * I->I[i].Max_IonsOfType, "IonsInitRead: FIon");
     I->I[i].FIon_old = (double *) Malloc(sizeof(double) * NDIM * I->I[i].Max_IonsOfType, "IonsInitRead: FIon_old");
@@ -136,48 +240,42 @@
     I->I[i].FIonL = (double *) Malloc(sizeof(double) * NDIM * I->I[i].Max_IonsOfType, "IonsInitRead: FIonL");
     I->I[i].FIonNL = (double *) Malloc(sizeof(double) * NDIM * I->I[i].Max_IonsOfType, "IonsInitRead: FIonNL");
+    I->I[i].FMagnetic = (double *) Malloc(sizeof(double) * NDIM * I->I[i].Max_IonsOfType, "IonsInitRead: FMagnetic");
+    I->I[i].FConstraint = (double *) Malloc(sizeof(double) * NDIM * I->I[i].Max_IonsOfType, "IonsInitRead: FConstraint");
     I->I[i].FEwald = (double *) Malloc(sizeof(double) * NDIM * I->I[i].Max_IonsOfType, "IonsInitRead: FEwald");
     I->I[i].alpha = (double *) Malloc(sizeof(double) * I->I[i].Max_IonsOfType, "IonsInitRead: alpha");
+    // now parse ion coordination
     for (j=0; j < I->I[i].Max_IonsOfType; j++) {
       I->I[i].alpha[j] = 2.;
     	sprintf(name,"Ion_Type%i_%i",i+1,j+1);
-			ParseForParameter(P->Call.out[ReadOut],source, name, 0, 1, 1, double_type, &R[0], critical);
-			ParseForParameter(P->Call.out[ReadOut],source, name, 0, 2, 1, double_type, &R[1], critical);
-			ParseForParameter(P->Call.out[ReadOut],source, name, 0, 3, 1, double_type, &R[2], critical);
-			ParseForParameter(P->Call.out[ReadOut],source, name, 0, 4, 1, int_type, &IMT, critical);
-      // change if coordinates were relative
-      if (relative) {
-        //fprintf(stderr,"(%i)Ion coordinates are relative %i ... \n",P->Par.me, relative);
-        RMat33Vec3(&I->I[i].R[0+NDIM*j], L->RealBasis, R); // multiply with real basis
+      for (d=0; d <NDIM; d++) // transfor to old and old_old coordinates as well
+        I->I[i].FMagnetic[d+j*NDIM] = 0.;
+      // first ones are starting positions
+      if ((count = ParseIonsCoordinatesAndVelocities(P, source, name, 1, relative, sigma, &I->I[i].R[NDIM*j], &I->I[i].U[NDIM*j], &IMT))) {
+        for (d=0; d <NDIM; d++) // transfor to old and old_old coordinates as well
+          I->I[i].R_old_old[d+NDIM*j] = I->I[i].R_old[d+NDIM*j] = I->I[i].R[d+NDIM*j];
+        I->I[i].IMT[j] = IMT;
+        if (P->Call.out[ReadOut]) fprintf(stderr, "(%i) R[%d,%d] = (%lg, %lg, %lg)\n", P->Par.me, i,j,I->I[i].R[0+NDIM*j],I->I[i].R[1+NDIM*j],I->I[i].R[2+NDIM*j]);
+        if (P->Call.out[ReadOut]) fprintf(stderr, "(%i) U[%d,%d] = (%lg, %lg, %lg)\n", P->Par.me, i,j,I->I[i].U[0+NDIM*j],I->I[i].U[1+NDIM*j],I->I[i].U[2+NDIM*j]);
+        if (P->Call.out[ReadOut]) fprintf(stderr, "(%i) IMT[%d,%d] = %i\n", P->Par.me, i,j,I->I[i].IMT[j]);
       } else {
-        for (k=0;k<NDIM;k++) // and copy to old vector
-          I->I[i].R[k+NDIM*j] = R[k];
-      }
-      if (IMT < 0 || IMT >= MaxIonMoveType) {
-				fprintf(stderr,"Bad Ion MoveType set to MoveIon for Ion (%i,%i)\n",i,j);
-				IMT = 0;
-      }
-        
-      if ((P->Call.out[ReadOut])) {  // rotate coordinates
-        fprintf(stderr,"(%i) coordinates of Ion %i: (x,y,z) = (",P->Par.me,j);
-        for(k=0;k<NDIM;k++) {
-          fprintf(stderr,"%lg ",I->I[i].R[k+NDIM*j]);
-          if (k != NDIM-1) fprintf(stderr,", ");
-          else fprintf(stderr,")\n");
+        Error(SomeError, "Could not find or parse coordinates of an ion");
+      }
+      // following ones are constrained motions
+      while (ParseIonsCoordinatesAndVelocities(P, source, name, ++count, relative, sigma, R, U, &IMT)) {
+        if (P->Call.out[ReadOut]) fprintf(stderr, "(%i) Constrained read %d ...\n", P->Par.me, count-1);
+        ptr = AddConstraintItem(&I->I[i], j);  // allocate memory
+        for (d=0; d< NDIM; d++) { // fill the constraint item
+          ptr->R[d] = R[d];
+          ptr->U[d] = U[d];
         }
-      }
-
-      I->I[i].IMT[j] = (enum IonMoveType)IMT;
-      SM(&I->I[i].R[NDIM*j], 1./Bohr, NDIM);
-      RMat33Vec3(R, L->InvBasis, &I->I[i].R[NDIM*j]);
-      for (d=0; d <NDIM; d++) {
-				while (R[d] < 0) 
-				  R[d] += 1.0;
-				while (R[d] >= 1.0) 
-				  R[d] -= 1.0;
-      }
-      RMat33Vec3(&I->I[i].R[NDIM*j], L->RealBasis, R);
-      for (d=0; d <NDIM; d++) {
-				I->I[i].R_old_old[d+NDIM*j] = I->I[i].R_old[d+NDIM*j] = I->I[i].R[d+NDIM*j];
-      }
+        ptr->IMT = IMT;
+        ptr->step = count-1;
+        if (P->Call.out[ReadOut]) fprintf(stderr, "(%i) R[%d,%d] = (%lg, %lg, %lg)\n", P->Par.me, i,j,ptr->R[0],ptr->R[1],ptr->R[2]);
+        if (P->Call.out[ReadOut]) fprintf(stderr, "(%i) U[%d,%d] = (%lg, %lg, %lg)\n", P->Par.me, i,j,ptr->U[0],ptr->U[1],ptr->U[2]);
+        if (P->Call.out[ReadOut]) fprintf(stderr, "(%i) IMT[%d,%d] = %i\n", P->Par.me, i,j,ptr->IMT);
+      }
+      //if (P->Call.out[ReadOut]) 
+      fprintf(stderr,"(%i) A total of %d additional constrained moves for ion (%d,%d) was parsed.\n", P->Par.me, count-2, i,j);
+
       I->Max_TotalIons++;
     }
Index: pcp/src/ions.h
===================================================================
--- pcp/src/ions.h	(revision e00f4732f333881727b6265e9f995246359ea74f)
+++ pcp/src/ions.h	(revision 961b75c543d6081b2758726da4c939e33cfc6b28)
@@ -121,5 +121,4 @@
 /* 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);
@@ -138,3 +137,5 @@
 void ResetForces(struct Problem *P);
 void CorrectVelocity(struct Problem *P);
+struct IonConstrained * AddConstraintItem(struct IonType *I, int ion);
+int RemoveConstraintItem(struct IonType *I, int ion);
 #endif
Index: pcp/src/output.c
===================================================================
--- pcp/src/output.c	(revision e00f4732f333881727b6265e9f995246359ea74f)
+++ pcp/src/output.c	(revision 961b75c543d6081b2758726da4c939e33cfc6b28)
@@ -606,10 +606,10 @@
     }
     // ... and parse critical ...
-    readnr += ParseForParameter(0,SrcPsiDoc,"Mintype",0,1,1,int_type,(int *)&read_type, test ? optional : critical);
-    readnr += ParseForParameter(0,SrcPsiDoc,"LevelNo",0,1,1,int_type,&LevelNo, test ? optional : critical);
-    readnr += 3*ParseForParameter(0,SrcPsiDoc,"GridNodes",0,3,1,row_int,&N[0], test ? optional : critical);
-    readnr += 2*ParseForParameter(0,SrcPsiDoc,"PsiNo",0,2,1,row_int,&GlobalNo[0], test ? optional : critical);
+    readnr += ParseForParameter(0,SrcPsiDoc,"Mintype",0,1,1,int_type,(int *)&read_type, 1, test ? optional : critical);
+    readnr += ParseForParameter(0,SrcPsiDoc,"LevelNo",0,1,1,int_type,&LevelNo,1,  test ? optional : critical);
+    readnr += 3*ParseForParameter(0,SrcPsiDoc,"GridNodes",0,3,1,row_int,&N[0], 1, test ? optional : critical);
+    readnr += 2*ParseForParameter(0,SrcPsiDoc,"PsiNo",0,2,1,row_int,&GlobalNo[0], 1, test ? optional : critical);
     // and optional items ...
-    if (ParseForParameter(0,SrcPsiDoc,"Epsilon",0,2,1,row_double,&Eps[0],optional))
+    if (ParseForParameter(0,SrcPsiDoc,"Epsilon",0,2,1,row_double,&Eps[0],1,optional))
       if ((P->Call.ReadSrcFiles == 1) && ((Eps[1] < R->RelEpsKineticEnergy) || (Eps[0] < R->RelEpsTotalEnergy))) {
         //fprintf(stderr,"(%i) Eps %lg %lg\tRelEps %lg %lg\n", P->Par.me, Eps[0], Eps[1], R->RelEpsTotalEnergy, R->RelEpsKineticEnergy);
Index: pcp/src/perturbed.c
===================================================================
--- pcp/src/perturbed.c	(revision e00f4732f333881727b6265e9f995246359ea74f)
+++ pcp/src/perturbed.c	(revision 961b75c543d6081b2758726da4c939e33cfc6b28)
@@ -2314,5 +2314,5 @@
   int wished = -1;
   FILE *file =  fopen(P->Call.MainParameterFile,"r");
-  if (!ParseForParameter(1,file,"Orbital",0,1,1,int_type,&wished, optional)) {
+  if (!ParseForParameter(1,file,"Orbital",0,1,1,int_type,&wished, 1, optional)) {
     fprintf(stderr,"Desired Orbital missing!\n");
     wished = -1;
Index: pcp/src/pseudo.c
===================================================================
--- pcp/src/pseudo.c	(revision e00f4732f333881727b6265e9f995246359ea74f)
+++ pcp/src/pseudo.c	(revision 961b75c543d6081b2758726da4c939e33cfc6b28)
@@ -506,13 +506,13 @@
     	if (P->Call.out[ReadOut]) fprintf(stderr,"%s was found\n",cpiInputFileName);
     int zeile = 1;
-    ParseForParameter(0,cpiInputFile, NULL, 0, 1, zeile, double_type, &PP->zval[it], critical);
-    ParseForParameter(0,cpiInputFile, NULL, 1, 2, zeile, int_type, &PP->nang[it], critical);
+    ParseForParameter(0,cpiInputFile, NULL, 0, 1, zeile, double_type, &PP->zval[it], 1, critical);
+    ParseForParameter(0,cpiInputFile, NULL, 1, 2, zeile, int_type, &PP->nang[it], 1, critical);
     I->TotalZval += PP->zval[it]*(I->I[it].Max_IonsOfType);
     PP->core[it] = (double *) Malloc(sizeof(double)*2, "InitPseudoRead: ");
     PP->rc[it] = (double *) Malloc(sizeof(double)*2, "InitPseudoRead: ");
-    ParseForParameter(0,cpiInputFile, NULL, 0, 1, zeile, double_type, &PP->core[it][0], critical);
-    ParseForParameter(0,cpiInputFile, NULL, 0, 2, zeile, double_type, &PP->rc[it][0], critical);
-    ParseForParameter(0,cpiInputFile, NULL, 0, 3, zeile, double_type, &PP->core[it][1], critical);
-    ParseForParameter(0,cpiInputFile, NULL, 1, 4, zeile, double_type, &PP->rc[it][1], critical);
+    ParseForParameter(0,cpiInputFile, NULL, 0, 1, zeile, double_type, &PP->core[it][0], 1, critical);
+    ParseForParameter(0,cpiInputFile, NULL, 0, 2, zeile, double_type, &PP->rc[it][0], 1, critical);
+    ParseForParameter(0,cpiInputFile, NULL, 0, 3, zeile, double_type, &PP->core[it][1], 1, critical);
+    ParseForParameter(0,cpiInputFile, NULL, 1, 4, zeile, double_type, &PP->rc[it][1], 1, critical);
     PP->rcl[it] = (double **) Malloc(sizeof(double*)*3, "InitPseudoRead: ");
     PP->al[it] =  (double **) Malloc(sizeof(double*)*3, "InitPseudoRead: ");
@@ -526,7 +526,7 @@
       PP->bl[it][il] = (double *) Malloc(sizeof(double)*3, "InitPseudoRead: ");
       for (ib = 0; ib < 3; ib++) {
-		    ParseForParameter(0,cpiInputFile, NULL, 0, 1, zeile, double_type, &PP->rcl[it][il][ib], critical);
-		    ParseForParameter(0,cpiInputFile, NULL, 0, 2, zeile, double_type, &PP->al[it][il][ib], critical);
-		    ParseForParameter(0,cpiInputFile, NULL, 1, 3, zeile, double_type, &PP->bl[it][il][ib], critical);
+		    ParseForParameter(0,cpiInputFile, NULL, 0, 1, zeile, double_type, &PP->rcl[it][il][ib], 1, critical);
+		    ParseForParameter(0,cpiInputFile, NULL, 0, 2, zeile, double_type, &PP->al[it][il][ib], 1, critical);
+		    ParseForParameter(0,cpiInputFile, NULL, 1, 3, zeile, double_type, &PP->bl[it][il][ib], 1, critical);
 				if (PP->rcl[it][il][ib] < MYEPSILON)
 				  PP->rcl[it][il][ib] = 0.0;
@@ -535,6 +535,6 @@
       }
     }
-    ParseForParameter(0,cpiInputFile, NULL, 0, 1, zeile, int_type, &PP->mmax[it], critical);
-    ParseForParameter(0,cpiInputFile, NULL, 1, 2, zeile, double_type, &PP->clog[it], critical);
+    ParseForParameter(0,cpiInputFile, NULL, 0, 1, zeile, int_type, &PP->mmax[it], 1, critical);
+    ParseForParameter(0,cpiInputFile, NULL, 1, 2, zeile, double_type, &PP->clog[it], 1, critical);
     if (PP->mmax[it] > PP->Mmax) PP->Mmax = PP->mmax[it];
     PP->clog[it] = log(PP->clog[it]);
@@ -547,12 +547,12 @@
       PP->v_loc[it][il] = (double *) Malloc(sizeof(double)*PP->mmax[it], "InitPseudoRead: ");
       for (j=0;j< PP->mmax[it];j++) {
-		    ParseForParameter(0,cpiInputFile, NULL, 0, 1, zeile, int_type, &m, critical);
-		    ParseForParameter(0,cpiInputFile, NULL, 0, 2, zeile, double_type, &PP->r[it][j], critical);
-		    ParseForParameter(0,cpiInputFile, NULL, 0, 3, zeile, double_type, &PP->R[it][il][j], critical);
-		    ParseForParameter(0,cpiInputFile, NULL, 1, 4, zeile, double_type, &PP->v_loc[it][il][j], critical);
+		    ParseForParameter(0,cpiInputFile, NULL, 0, 1, zeile, int_type, &m, 1, critical);
+		    ParseForParameter(0,cpiInputFile, NULL, 0, 2, zeile, double_type, &PP->r[it][j], 1, critical);
+		    ParseForParameter(0,cpiInputFile, NULL, 0, 3, zeile, double_type, &PP->R[it][il][j], 1, critical);
+		    ParseForParameter(0,cpiInputFile, NULL, 1, 4, zeile, double_type, &PP->v_loc[it][il][j], 1, critical);
       }
       if (il < PP->nang[it]-1) {
-		    ParseForParameter(0,cpiInputFile, NULL, 0, 1, zeile, int_type, &PP->mmax[it], critical);
-		    ParseForParameter(0,cpiInputFile, NULL, 1, 2, zeile, double_type, &PP->clog[it], critical);
+		    ParseForParameter(0,cpiInputFile, NULL, 0, 1, zeile, int_type, &PP->mmax[it], 1, critical);
+		    ParseForParameter(0,cpiInputFile, NULL, 1, 2, zeile, double_type, &PP->clog[it], 1, critical);
 				PP->clog[it] = log(PP->clog[it]);
       }
@@ -560,6 +560,6 @@
     I->I[it].corecorr = NotCoreCorrected;
     count = 0;
-    count += ParseForParameter(0,cpiInputFile, NULL, 0, 1, zeile, double_type, &dummyr, optional);
-    count += ParseForParameter(0,cpiInputFile, NULL, 0, 2, zeile, double_type, &dummycorewave, optional);
+    count += ParseForParameter(0,cpiInputFile, NULL, 0, 1, zeile, double_type, &dummyr, 1, optional);
+    count += ParseForParameter(0,cpiInputFile, NULL, 0, 2, zeile, double_type, &dummycorewave, 1, optional);
     //fprintf(stderr,"(%i) %lg %lg\n",P->Par.me,dummyr,dummycorewave);
     if (count == 2) {
@@ -570,11 +570,11 @@
       PP->corecorr = CoreCorrected;
       PP->corewave[it][0] = dummycorewave/(4.*PI);
-	    ParseForParameter(0,cpiInputFile, NULL, 0, 3, zeile, double_type, &dummyr, critical);
-	    ParseForParameter(0,cpiInputFile, NULL, 1, 4, zeile, double_type, &dummycorewave, critical);
+	    ParseForParameter(0,cpiInputFile, NULL, 0, 3, zeile, double_type, &dummyr, 1, critical);
+	    ParseForParameter(0,cpiInputFile, NULL, 1, 4, zeile, double_type, &dummycorewave, 1, critical);
       for (j=1;j < PP->mmax[it]; j++) {
-		    ParseForParameter(0,cpiInputFile, NULL, 0, 1, zeile, double_type, &PP->r[it][j], critical);
-		    ParseForParameter(0,cpiInputFile, NULL, 0, 2, zeile, double_type, &PP->corewave[it][j], critical);
-		    ParseForParameter(0,cpiInputFile, NULL, 0, 3, zeile, double_type, &dummyr, critical);
-		    ParseForParameter(0,cpiInputFile, NULL, 1, 4, zeile, double_type, &dummycorewave, critical);
+		    ParseForParameter(0,cpiInputFile, NULL, 0, 1, zeile, double_type, &PP->r[it][j], 1, critical);
+		    ParseForParameter(0,cpiInputFile, NULL, 0, 2, zeile, double_type, &PP->corewave[it][j], 1, critical);
+		    ParseForParameter(0,cpiInputFile, NULL, 0, 3, zeile, double_type, &dummyr, 1, critical);
+		    ParseForParameter(0,cpiInputFile, NULL, 1, 4, zeile, double_type, &dummycorewave, 1, critical);
 				PP->corewave[it][j] /= (4.*PI);
       }
Index: pcp/src/wannier.c
===================================================================
--- pcp/src/wannier.c	(revision e00f4732f333881727b6265e9f995246359ea74f)
+++ pcp/src/wannier.c	(revision 961b75c543d6081b2758726da4c939e33cfc6b28)
@@ -1319,5 +1319,5 @@
         sprintf(tagname,"Psi%d_Lev%d",i,R->LevSNo);
         signal = 0;
-        if (!ParseForParameter(0,SpreadFile,tagname,0,3,1,row_double,WannierCentre,optional)) {
+        if (!ParseForParameter(0,SpreadFile,tagname,0,3,1,row_double,WannierCentre,1,optional)) {
           //Error(SomeError,"ParseWannierFile: Parsing WannierCentre failed");
           if (MPI_Bcast(&signal,1,MPI_INT,0,P->Par.comm_ST) != MPI_SUCCESS)
@@ -1325,5 +1325,5 @@
           return 0;
         }
-        if (!ParseForParameter(0,SpreadFile,tagname,0,4,1,double_type,&WannierCentre[NDIM],optional)) {
+        if (!ParseForParameter(0,SpreadFile,tagname,0,4,1,double_type,&WannierCentre[NDIM],1,optional)) {
           //Error(SomeError,"ParseWannierFile: Parsing WannierSpread failed");
           if (MPI_Bcast(&signal,1,MPI_INT,0,P->Par.comm_ST) != MPI_SUCCESS)
