debugging of dipole iterations

examples/amoeba/in.water_box.amoeba

@@ -43,5 +43,5 @@ dump_modify     1 sort id
 
 fix             1 all nve
 
-thermo          10
+thermo          1
-run             100
+run             5

examples/amoeba/in.water_box.hippo

@@ -43,5 +43,5 @@ dump_modify     1 sort id
 
 fix             1 all nve
 
-thermo          10
+thermo          1
-run             100
+run             5

src/AMOEBA/amoeba_induce.cpp

@@ -25,6 +25,10 @@
 #include "memory.h"
 #include "error.h"
 
+// DEBUG
+
+#include "update.h"
+
 #include <cmath>
 #include <cstring>
 
@@ -381,6 +385,10 @@ void PairAmoeba::induce()
       }
     }
 
+    // DEBUG
+
+    printf("CG iteration count = %d\n",iter);
+
     // terminate the calculation if dipoles failed to converge
     // NOTE: could make this an error
 

src/AMOEBA/pair_amoeba.cpp

@@ -1756,7 +1756,10 @@ void PairAmoeba::precond_neigh()
   choose(USOLV);
 
   //double off = sqrt(off2);
+  //if (comm->me == 0) printf("PCN off %g off2 %g\n",off,off2);
   //off2 = (off + neighbor->skin) * (off + neighbor->skin);
+  //if (comm->me == 0) 
+  // printf("PCNnew off %g off2 %g skin %g\n",sqrt(off2),off2,neighbor->skin);
 
   // atoms and neighbor list
 
@@ -1770,6 +1773,13 @@ void PairAmoeba::precond_neigh()
   // store all induce neighs of owned atoms within shorter cutoff
   // scan longer-cutoff neighbor list of I
 
+
+  // DEBUG
+
+  int *ncount;
+  memory->create(ncount,atom->nlocal,"amoeba:ncount");
+  for (int i = 0; i < atom->nlocal; i++) ncount[i] = 0;
+
   ipage_precond->reset();
 
   for (ii = 0; ii < inum; ii++) {
@@ -1793,15 +1803,38 @@ void PairAmoeba::precond_neigh()
       delz = ztmp - x[j][2];
       rsq = delx*delx + dely*dely + delz*delz;
 
-      if (rsq < off2) neighptr[n++] = jlist[jj];
+      if (rsq < off2) {
+        neighptr[n++] = jlist[jj];
+
+        // DEBUG
+
+        int jcount = atom->map(atom->tag[j]);
+        ncount[jcount]++;
+      }
     }
 
     firstneigh_precond[i] = neighptr;
     numneigh_precond[i] = n;
     ipage_precond->vgot(n);
+
+    // DEBUG
+
+    ncount[i] += n;
   }
+
+  // DEBUG
+
+  if (update->ntimestep == 0) {
+    for (int i = 1; i <= atom->nlocal; i++) {
+      int ilocal = atom->map(i);
+      printf("PRECOND id %d ilocal %d numneigh %d\n",i,ilocal,ncount[ilocal]);
+    }
+  }
+
+  memory->destroy(ncount);
 }
 
+
 /* ----------------------------------------------------------------------
    allocate Vdwl arrays
    note that n_amclass = # of classes in Tinker PRM file