Index: pcp/src/ions.c
===================================================================
--- pcp/src/ions.c	(revision 3716b2ef4660a374ad11bc3bc01dd9b133c6ce37)
+++ pcp/src/ions.c	(revision c5bdb2362e31e9ba6a11194a01d498a04e79d94e)
@@ -438,4 +438,8 @@
 
 /** Calculated Ewald force.
+ * The basic idea of the ewald method is to add a countercharge - here of gaussian form -
+ * which shields the long-range charges making them effectively short-ranged, while the
+ * countercharges themselves can be analytically (here by (hidden) Fouriertransform: it's
+ * added to the electron density) integrated and withdrawn. 
  * \f[
  *  E_{K-K} = \frac{1}{2} \sum_L \sum_{K=\{(I_s,I_a,J_s,J_a)|R_{I_s,I_a} - R_{J_s,J_a} - L\neq 0\}} \frac{Z_{I_s} Z_{J_s}} {|R_{I_s,I_a} - R_{J_s,J_a} - L|} 
@@ -448,9 +452,9 @@
  * in a circular motion from the current cell to the outside up to Ions::R_cut.
  * \param *P Problem at hand
- * \param first additional calculation beforehand if != 0
+ * \param first if != 0 table mirror cells to take into (local!) account of ewald summation 
  */
 void CalculateEwald(struct Problem *P, int first) 
 {
-  int r,is,is1,is2,ia1,ia2,ir,ir2,i,j,k,ActualVec,MaxVec,MaxCell,Is_Neighb_Cell,LocalActualVec;
+  int r,is1,is2,ia1,ia2,ir,ir2,i,j,k,ActualVec,MaxVec,MaxCell,Is_Neighb_Cell,LocalActualVec; //,is
   int MaxPar=P->Par.procs, ParMe=P->Par.me, MaxLocalVec, StartVec;
   double rcut2,r2,erre2,addesr,addpre,arg,fac,gkl,esr,fxx,repand;
@@ -462,122 +466,132 @@
   struct PseudoPot *PP = &P->PP;
   struct Ions *I = &P->Ion;
-  if (I->R_cut != 0.0) {
-    if (first) {
-      SpeedMeasure(P, InitSimTime, StopTimeDo);
-      rcut2 = I->R_cut*I->R_cut;
-      ActualVec =0;
-      MaxCell = (int)5.*I->R_cut/pow(L->Volume,1./3.);
-      if (MaxCell < 2) MaxCell = 2;
-      for (i = -MaxCell; i <= MaxCell; i++) {
-        for (j = -MaxCell; j <= MaxCell; j++) {
-  				for (k = -MaxCell; k <= MaxCell; k++) {
-  				  r2 = 0.0;
-  				  for (ir=0; ir <3; ir++) {
-  				    t[ir] = i*L->RealBasis[0*NDIM+ir]+j*L->RealBasis[1*NDIM+ir]+k*L->RealBasis[2*NDIM+ir];
-  				    r2 = t[ir]*t[ir];
-  				  }
-  				  Is_Neighb_Cell = ((abs(i)<=1) && (abs(j)<=1) && (abs(k)<=1));
-  				  if ((r2 <= rcut2) || Is_Neighb_Cell) {
-  				    ActualVec++;
-  				  }
-  				}
-        }
-      }
-      MaxVec = ActualVec;
-      I->MaxVec = ActualVec;
-      MaxLocalVec = MaxVec / MaxPar;
-      StartVec = ParMe*MaxLocalVec;
-      r = MaxVec % MaxPar;
-      if (ParMe >= r) {
-        StartVec += r;
-      } else {
-        StartVec += ParMe;
-        MaxLocalVec++;
-      }
-      I->MaxLocalVec = MaxLocalVec;
-      LocalActualVec = ActualVec = 0;
+  if (first) {
+    SpeedMeasure(P, InitSimTime, StopTimeDo);
+    rcut2 = I->R_cut*I->R_cut;
+    ActualVec = 0;  // counts number of cells taken into account
+    MaxCell = (int)(5.*I->R_cut/pow(L->Volume,1./3.));    // the number of cells in each direction is prop. to axis length over cutoff
+    if (MaxCell < 2) MaxCell = 2;     // take at least 2
+    for (i = -MaxCell; i <= MaxCell; i++) {
+      for (j = -MaxCell; j <= MaxCell; j++) {
+				for (k = -MaxCell; k <= MaxCell; k++) {
+				  r2 = 0.;
+				  for (ir=0; ir <NDIM; ir++) {  // t is the offset vector pointing to distant cell (only diagonal entries)
+				    t[ir] = i*L->RealBasis[0*NDIM+ir]+j*L->RealBasis[1*NDIM+ir]+k*L->RealBasis[2*NDIM+ir];
+				    r2 = t[ir]*t[ir];   // is squared distance to other cell 
+				  }
+				  Is_Neighb_Cell = ((abs(i)<=1) && (abs(j)<=1) && (abs(k)<=1));   // whether it's directly adjacent
+				  if ((r2 <= rcut2) || Is_Neighb_Cell) {    // if it's either adjacent or within cutoff, count it in
+				    ActualVec++;
+				  }
+				}
+      }
+    }
+    MaxVec = ActualVec;   // store number of cells
+    I->MaxVec = ActualVec;
+    MaxLocalVec = MaxVec / MaxPar;  // share among processes
+    StartVec = ParMe*MaxLocalVec; // for this process start at its patch
+    r = MaxVec % MaxPar;  // rest not fitting...
+    if (ParMe >= r) {   // ones who don't get extra cells have to shift their start vector
+      StartVec += r;
+    } else {    // others get extra cell and have to shift also by a bit
+      StartVec += ParMe;
+      MaxLocalVec++;
+    }
+    I->MaxLocalVec = MaxLocalVec;
+    LocalActualVec = ActualVec = 0;
+    // now go through the same loop again and note down every offset vector for cells in this process' patch
+    if (MaxLocalVec != 0) {
       I->RLatticeVec = (double *)
-        Realloc(I->RLatticeVec, sizeof(double)*NDIM*MaxLocalVec, "CalculateEwald:");
-      for (i = -MaxCell; i <= MaxCell; i++) {
-        for (j = -MaxCell; j <= MaxCell; j++) {
-  				for (k = -MaxCell; k <= MaxCell; k++) {
-  				  r2 = 0.0;
-  				  for (ir=0; ir <3; ir++) {
-  				    t[ir] = i*L->RealBasis[0*NDIM+ir]+j*L->RealBasis[1*NDIM+ir]+k*L->RealBasis[2*NDIM+ir];
-  				    r2 = t[ir]*t[ir];
-  				  }
-  				  Is_Neighb_Cell = ((abs(i)<=1) && (abs(j)<=1) && (abs(k)<=1));
-  				  if ((r2 <= rcut2) || Is_Neighb_Cell) {
-  				    if (ActualVec >= StartVec && ActualVec < StartVec+MaxLocalVec) {
-  				      I->RLatticeVec[0+NDIM*LocalActualVec] = t[0];
-  				      I->RLatticeVec[1+NDIM*LocalActualVec] = t[1];
-  				      I->RLatticeVec[2+NDIM*LocalActualVec] = t[2];
-  				      LocalActualVec++;
-  				    }
-  				    ActualVec++;
-  				  }
-  				}
-        }
-      }
-      SpeedMeasure(P, InitSimTime, StartTimeDo);
-    }
-    SpeedMeasure(P, EwaldTime, StartTimeDo);
-    esr = 0.0;
-    for (is1=0;is1 < I->Max_Types; is1++) 
-      for (ia1=0;ia1 < I->I[is1].Max_IonsOfType; ia1++)
-        for (i=0; i< NDIM; i++)
-  				I->FTemp[i+NDIM*(ia1)] = 0;
-    for (is1=0;is1 < I->Max_Types; is1++) {
-      for (ia1=0;ia1 < I->I[is1].Max_IonsOfType; ia1++) {
-        for (i=0;i<3;i++) {
-  				R1[i]=I->I[is1].R[i+NDIM*ia1];
-        }
-        for (ir2=0;ir2<I->MaxLocalVec;ir2++) {
-  				for (is2=0;is2 < I->Max_Types; is2++) {
-  				  gkl=1./sqrt(I->I[is1].rgauss*I->I[is1].rgauss + I->I[is2].rgauss*I->I[is2].rgauss);
-  				  for (ia2=0;ia2<I->I[is2].Max_IonsOfType; ia2++) {
-  				    for (i=0;i<3;i++) {
-  				      R2[i]=I->RLatticeVec[i+NDIM*ir2]+I->I[is2].R[i+NDIM*ia2];
-  				    }
-  				    erre2=0.0;
-  				    for (i=0;i<3;i++) {
-  				      R3[i] = R1[i]-R2[i];
-  				      erre2+=R3[i]*R3[i];
-  				    }
-  				    if (erre2 > MYEPSILON) {
-  				      arg=sqrt(erre2);
-  				      fac=PP->zval[is1]*PP->zval[is2]/arg*0.5;
-  				      
-  				      arg *= gkl;
-  				      addesr = derf(arg);
-  				      addesr = (1.0-addesr)*fac;
-  				      esr += addesr;
-  				      addpre=exp(-arg*arg)*gkl;
-  				      addpre=PP->fac1sqrtPI*PP->zval[is1]*PP->zval[is2]*addpre;
-  				      repand=(addesr+addpre)/erre2;
-  				      for (i=0;i<3;i++) {
-  								fxx=repand*R3[i];
-  								/*fprintf(stderr,"%i %i %i %i %i %g\n",is1+1,ia1+1,is2+1,ia2+1,i+1,fxx);*/
-  								I->FTemp[i+NDIM*(ia1)] += fxx;
-  								I->FTemp[i+NDIM*(ia2)] -= fxx;
-  				      }
-  				    }
-  				  }
-  				}
-        }
-      }
-    }
-  } else {
-    esr = 0.0;
-    for (is1=0;is1 < I->Max_Types; is1++) 
-      for (ia1=0;ia1 < I->I[is1].Max_IonsOfType; ia1++)
-        for (i=0; i< NDIM; i++)
-          I->FTemp[i+NDIM*(ia1)] = 0.;
+        Realloc(I->RLatticeVec, sizeof(double)*NDIM*MaxLocalVec, "CalculateEwald: I->RLatticeVec");
+    } else {
+      fprintf(stderr,"(%i) Warning in Ewald summation: Got MaxLocalVec == 0\n", P->Par.me);
+    }
+    for (i = -MaxCell; i <= MaxCell; i++) {
+      for (j = -MaxCell; j <= MaxCell; j++) {
+				for (k = -MaxCell; k <= MaxCell; k++) {
+				  r2 = 0.;
+				  for (ir=0; ir <NDIM; ir++) {  // create offset vector from diagonal entries
+				    t[ir] = i*L->RealBasis[0*NDIM+ir]+j*L->RealBasis[1*NDIM+ir]+k*L->RealBasis[2*NDIM+ir];
+				    r2 = t[ir]*t[ir];
+				  }
+				  Is_Neighb_Cell = ((abs(i)<=1) && (abs(j)<=1) && (abs(k)<=1));
+				  if ((r2 <= rcut2) || Is_Neighb_Cell) {  // if adjacent or within cutoff ...
+				    if (ActualVec >= StartVec && ActualVec < StartVec+MaxLocalVec) { // ... and belongs to patch, store 'em
+				      I->RLatticeVec[0+NDIM*LocalActualVec] = t[0];
+				      I->RLatticeVec[1+NDIM*LocalActualVec] = t[1];
+				      I->RLatticeVec[2+NDIM*LocalActualVec] = t[2];
+				      LocalActualVec++;
+				    }
+				    ActualVec++;
+				  }
+				}
+      }
+    }
+    SpeedMeasure(P, InitSimTime, StartTimeDo);
+  }
+  SpeedMeasure(P, EwaldTime, StartTimeDo);
+
+  // first set everything to zero
+  esr = 0.0;
+  //for (is1=0;is1 < I->Max_Types; is1++) 
+  // then take each ion of each type once 
+  for (is1=0;is1 < I->Max_Types; is1++) {
+    // reset FTemp for IonType
+    for (ia1=0;ia1 < I->I[is1].Max_IonsOfType; ia1++)
+      for (i=0; i< NDIM; i++)
+        I->FTemp[i+NDIM*(ia1)] = 0.;
+    // then sum for each on of the IonType
+    for (ia1=0;ia1 < I->I[is1].Max_IonsOfType; ia1++) {
+      for (i=0;i<NDIM;i++) {
+				R1[i]=I->I[is1].R[i+NDIM*ia1];    // R1 is local coordinate of current first ion
+      }
+      for (ir2=0;ir2<I->MaxLocalVec;ir2++) {  // for every cell of local patch (each with the same atoms)
+				for (is2=0;is2 < I->Max_Types; is2++) { // and for each ion of each type a second time
+          // gkl = 1./(sqrt(r_is1,gauss^2 + r_is2,gauss^2))
+				  gkl=1./sqrt(I->I[is1].rgauss*I->I[is1].rgauss + I->I[is2].rgauss*I->I[is2].rgauss);
+				  for (ia2=0;ia2<I->I[is2].Max_IonsOfType; ia2++) {
+				    for (i=0;i<3;i++) {
+				      R2[i]=I->RLatticeVec[i+NDIM*ir2]+I->I[is2].R[i+NDIM*ia2]; // R2 is coordinate of current second ion in the distant cell!
+				    }
+
+				    erre2=0.0;
+				    for (i=0;i<NDIM;i++) {
+				      R3[i] = R1[i]-R2[i];    // R3 = R_is1,ia1 - R_is2,ia2 - L
+				      erre2+=R3[i]*R3[i];     // erre2 = |R_is1,ia1 - R_is2,ia2 - L|^2
+				    }
+				    if (erre2 > MYEPSILON) {  // appears as one over ... in fac, thus check if zero
+              arg = sqrt(erre2); // arg = |R_is1,ia1 - R_is2,ia2 - L|
+				      fac=PP->zval[is1]*PP->zval[is2]/arg*0.5;    // fac = -1/2 * Z_is1 * Z_is2 / arg
+				      
+				      arg *= gkl;             // arg = |R_is1,ia1 - R_is2,ia2 - L|/(sqrt(r_is1,gauss^2 + r_is2,gauss^2))
+				      addesr = derf(arg);  
+              addesr = (1.-addesr)*fac; // complementary error function: addesr = 1/2*erfc(arg)/|R_is1,ia1 - R_is2,ia2 - L|
+				      esr += addesr;
+              addpre=exp(-arg*arg)*gkl;
+              addpre = PP->fac1sqrtPI*PP->zval[is1]*PP->zval[is2]*addpre;  // addpre = exp(-|R_is1,ia1 - R_is2,ia2 - L|^2/(r_is1,gauss^2 + r_is2,gauss^2))/(sqrt(pi)*sqrt(r_is1,gauss^2 + r_is2,gauss^2))
+              repand = (addesr+addpre)/erre2;
+//              if ((fabs(repand) > MYEPSILON)) {
+//                fprintf(stderr,"(%i) Ewald[is%d,ia%d/is%d,ia%d,(%d)]: Z1 %lg, Z2 %lg\tpre %lg\t esr %lg\t fac %lg\t arg %lg\tR3 (%lg,%lg,%lg)\n", P->Par.me,  is1, ia1, is2, ia2, ir2, PP->zval[is1], PP->zval[is2],addpre,addesr, fac, arg, R3[0], R3[1], R3[2]);
+//              }
+				      for (i=0;i<NDIM;i++) {
+								fxx=repand*R3[i];
+//                if ((fabs(repand) > MYEPSILON)) {
+//								  fprintf(stderr,"%i %i %i %i %i %g\n",is1+1,ia1+1,is2+1,ia2+1,i+1,fxx);
+//                }
+								I->FTemp[i+NDIM*(ia1)] += fxx*2.0;      // actio = reactio?
+								//I->FTemp[i+NDIM*(ia2)] -= fxx;
+				      }
+				    }
+				  }
+				}
+      }
+    }
+    MPI_Allreduce (I->FTemp , I->I[is1].FEwald, NDIM*I->I[is1].Max_IonsOfType, MPI_DOUBLE, MPI_SUM, P->Par.comm);  
   }
   SpeedMeasure(P, EwaldTime, StopTimeDo);
-  for (is=0;is < I->Max_Types; is++) {
-    MPI_Allreduce (I->FTemp , I->I[is].FEwald, NDIM*I->I[is].Max_IonsOfType, MPI_DOUBLE, MPI_SUM, P->Par.comm);
-  }
+  //for (is=0;is < I->Max_Types; is++) {
+  //}
   MPI_Allreduce ( &esr, &L->E->AllTotalIonsEnergy[EwaldEnergy], 1, MPI_DOUBLE, MPI_SUM, P->Par.comm);
+//  fprintf(stderr,"\n");
 }
 
