Index: pcp/src/perturbed.c
===================================================================
--- pcp/src/perturbed.c	(revision 51af4a4097ffdd6f0f9ecf4568016f5443005be1)
+++ pcp/src/perturbed.c	(revision c7639399feaf5dc339b10ce8df7efe1d5b741622)
@@ -2939,4 +2939,5 @@
   
   if (P->Call.out[ValueOut]) fprintf(stderr,"(%i) magnetic susceptibility tensor \\Chi_ij = \n",P->Par.me);
+  if (P->Call.out[ReadOut]) fprintf(stderr,"\n");
   for (in=0; in<NDIM; in++) { // index i of integrand vector component
     for(dex=0;dex<4;dex++)  // initialise cross lookup
@@ -2967,5 +2968,5 @@
       I->I[0].chi[in+dex*NDIM] = Chi[in+dex*NDIM];
       Chi[in+dex*NDIM] *= Lat->Volume*loschmidt_constant; // factor for _molar_ susceptibility
-      if (P->Call.out[ValueOut]) {
+      if (P->Call.out[ReadOut]) {
         fprintf(stderr,"%e\t", Chi[in+dex*NDIM]);
         if (dex == NDIM-1) fprintf(stderr,"\n");
@@ -3012,9 +3013,16 @@
     A += pow(-1,i)*pow(0.5*(Chi[in+dex*NDIM]-Chi[dex+in*NDIM]),2);
   }
-  if (P->Call.out[ValueOut]) {
+  if (P->Call.out[ReadOut]) {
     fprintf(stderr,"(%i) converted to Principal Axis System\n==================\nDiagonal entries:", P->Par.me);
     for (i=0;i<NDIM;i++)
       fprintf(stderr,"\t%lg",gsl_vector_get(eval,i));
-    fprintf(stderr,"\nsusceptib. : %e\n", iso);
+  }
+  if (P->Call.out[ValueOut]) {
+    if (P->Call.out[ReadOut])
+      fprintf(stderr,"\nsusceptib. : %e\n", iso);
+    else
+      fprintf(stderr,"%e\n", iso);
+  }
+  if (P->Call.out[ReadOut]) {
     fprintf(stderr,"anisotropy : %e\n", delta_chi);
     fprintf(stderr,"asymmetry  : %e\n", eta);
@@ -3119,5 +3127,5 @@
   for (in=0;in<NDIM*NDIM;in++) {
     CalculateOneDensityC(Lat, R->LevS, Dens0, CurrentDensity[in], CurrentDensityC[in], factorDC);
-    TestReciprocalCurrent(P, CurrentDensityC[in], GArray, in);
+    //TestReciprocalCurrent(P, CurrentDensityC[in], GArray, in);
   }
 
@@ -3167,8 +3175,4 @@
         }
       }
-      // fabs() all sigma values, as we need them as a positive density: OutputVis plots them in logarithmic scale and
-      // thus cannot deal with negative values!
-      for (i=0; i< Dens0->LocalSizeR; i++)
-        sigma_real[in+dex*NDIM][i] = fabs(sigma_real[in+dex*NDIM][i]);
     }
   }
@@ -3194,11 +3198,12 @@
   for (it=0; it < I->Max_Types; it++) {  // integration over all types
     for (ion=0; ion < I->I[it].Max_IonsOfType; ion++) { // and each ion of type
-      if (P->Call.out[ValueOut]) fprintf(stderr,"(%i) Shielding Tensor for Ion %i of element %s \\sigma_ij = \n",P->Par.me, ion, I->I[it].Name);
+      if (P->Call.out[ValueOut]) fprintf(stderr,"(%i) Shielding Tensor for Ion %i of element %s \\sigma_ij = ",P->Par.me, ion, I->I[it].Name);
+      if (P->Call.out[ReadOut]) fprintf(stderr,"\n");
       for (in=0; in<NDIM; in++) { // index i of vector component in integrand
         for (dex=0; dex<NDIM; dex++) {// index j of B component derivation in current density tensor
           gsl_matrix_complex_set(H,in,dex,gsl_complex_rect((I->I[it].sigma_rezi[ion][in+dex*NDIM]+I->I[it].sigma_rezi[ion][dex+in*NDIM])/2.,0));
-          if (P->Call.out[ValueOut]) fprintf(stderr,"%e\t", I->I[it].sigma_rezi[ion][in+dex*NDIM]);
+          if (P->Call.out[ReadOut]) fprintf(stderr,"%e\t", I->I[it].sigma_rezi[ion][in+dex*NDIM]);
         }
-        if (P->Call.out[ValueOut]) fprintf(stderr,"\n");
+        if (P->Call.out[ReadOut]) fprintf(stderr,"\n");
       }
       // output tensor to file
@@ -3238,9 +3243,16 @@
         A += pow(-1,i)*pow(0.5*(I->I[it].sigma_rezi[ion][in+dex*NDIM]-I->I[it].sigma_rezi[ion][dex+in*NDIM]),2);
       }
-      if (P->Call.out[ValueOut]) {
+      if (P->Call.out[ReadOut]) {
         fprintf(stderr,"(%i) converted to Principal Axis System\n==================\nDiagonal entries:", P->Par.me);
         for (i=0;i<NDIM;i++)
           fprintf(stderr,"\t%lg",gsl_vector_get(eval,i));
-        fprintf(stderr,"\nshielding  : %e\n", iso);
+      }
+      if (P->Call.out[ValueOut]) {
+        if (P->Call.out[ReadOut])
+          fprintf(stderr,"\nshielding  : %e\n", iso);
+        else
+          fprintf(stderr,"%e\n", iso);
+      }
+      if (P->Call.out[ReadOut]) {
         fprintf(stderr,"anisotropy : %e\n", delta_sigma);
         fprintf(stderr,"asymmetry  : %e\n", eta);
@@ -3284,5 +3296,13 @@
     }
   }
-   
+
+  // fabs() all sigma values, as we need them as a positive density: OutputVis plots them in logarithmic scale and
+  // thus cannot deal with negative values!
+  for (in=0; in<NDIM; in++) {// index i of vector component in integrand
+    for (dex=0; dex<NDIM; dex++) { // index j of B component derivation in current density tensor
+      for (i=0; i< Dens0->LocalSizeR; i++)
+        sigma_real[in+dex*NDIM][i] = fabs(sigma_real[in+dex*NDIM][i]);
+    }
+  }   
   if (Lev0->LevelNo == 0) {
     if (!P->Par.me && P->Call.out[NormalOut]) fprintf(stderr,"(%i)Output of NICS map ...\n", P->Par.me);
@@ -3371,5 +3391,6 @@
   for (it=0; it < I->Max_Types; it++) {  // integration over all types
     for (ion=0; ion < I->I[it].Max_IonsOfType; ion++) { // and each ion of type
-      if (P->Call.out[ValueOut]) fprintf(stderr,"(%i) Magnetic dipole moment Tensor for Ion %i of element %s \\moment_ij = \n",P->Par.me, ion, I->I[it].Name);
+      if (P->Call.out[ValueOut]) fprintf(stderr,"(%i) Magnetic dipole moment Tensor for Ion %i of element %s \\moment_ij = ",P->Par.me, ion, I->I[it].Name);
+      if (P->Call.out[ReadOut]) fprintf(stderr,"\n");
       for (in=0; in<NDIM; in++) {// index i of vector component in integrand
         for (dex=0; dex<NDIM; dex++) { // index j of B component derivation in current density tensor
@@ -3394,7 +3415,7 @@
           MPI_Allreduce ( &moment[in+dex*NDIM], &Moment[in+dex*NDIM], 1, MPI_DOUBLE, MPI_SUM, P->Par.comm_ST_Psi);  // sum "LocalSize to TotalSize"
           I->I[it].moment[ion][in+dex*NDIM] = Moment[in+dex*NDIM];          
-          if (P->Call.out[ValueOut]) fprintf(stderr," %e", Moment[in+dex*NDIM]);
+          if (P->Call.out[ReadOut]) fprintf(stderr," %e", Moment[in+dex*NDIM]);
         }
-        if (P->Call.out[ValueOut]) fprintf(stderr,"\n");
+        if (P->Call.out[ReadOut]) fprintf(stderr,"\n");
       }
       // store symmetrized matrix
@@ -3441,9 +3462,16 @@
         A += pow(-1,i)*pow(0.5*(Moment[in+dex*NDIM]-Moment[dex+in*NDIM]),2);
       }
-      if (P->Call.out[ValueOut]) {
+      if (P->Call.out[ReadOut]) {
         fprintf(stderr,"(%i) converted to Principal Axis System\n==================\nDiagonal entries:", P->Par.me);
         for (i=0;i<NDIM;i++)
           fprintf(stderr,"\t%lg",gsl_vector_get(eval,i));
-        fprintf(stderr,"\nshielding  : %e\n", iso);
+      }
+      if (P->Call.out[ValueOut]) {
+        if (P->Call.out[ReadOut])
+          fprintf(stderr,"moment  : %e\n", iso);
+        else
+          fprintf(stderr,"%e\n", iso);
+      }
+      if (P->Call.out[ReadOut]) {
         fprintf(stderr,"anisotropy : %e\n", delta_moment);
         fprintf(stderr,"asymmetry  : %e\n", eta);
@@ -3451,6 +3479,5 @@
         fprintf(stderr,"A          : %e\n", A);
         fprintf(stderr,"==================\n");
-        
-      }
+    }
       gsl_vector_free(eval);
     }
