-  for (int ipos = 0; ipos < ipow(length,dimension); ipos++)
-    {
-      for (int idir = 0; idir < 2*dimension; idir++)
-       {
-         if ( idir < dimension ) //pos direction
-           {
-             if ( ( ipos/ipow(length,idir%dimension) ) % length == length-1 )
-               nfield[ipos*2*dimension + idir] = ipos - ipow(length,idir%dimension)*(length-1);
-             else
-               nfield[ipos*2*dimension + idir] = ipos + ipow(length,idir%dimension);
-           }
-         else //negative dir
-           {
-             if ( ( ipos/ipow(length,idir%dimension) ) % length == 0 )
-               nfield[ipos*2*dimension + idir] = ipos + ipow(length,idir%dimension)*(length-1);
-             else
-               nfield[ipos*2*dimension + idir] = ipos - ipow(length,idir%dimension);
-           }
-         
-       }
-    }
+  for( int ipos=0; ipos<totsize; ipos++ )
+    for( int idir=0; idir<2*dimension; idir++) {
+      if( idir < dimension ) { // positive direction
+       if ( ( ipos/dirstep(idir%dimension) ) % len[idir%dimension] == len[idir%dimension]-1 )
+         nfield[ipos*2*dimension + idir] = ipos - dirstep(idir%dimension)*(len[idir%dimension]-1);
+       else
+         nfield[ipos*2*dimension + idir] = ipos + dirstep(idir%dimension);
+      }
+      else { // negative direction
+       if ( ( ipos/dirstep(idir%dimension) ) % len[idir%dimension] == 0 )
+         nfield[ipos*2*dimension + idir] = ipos + dirstep(idir%dimension)*(len[idir%dimension]-1);
+       else
+         nfield[ipos*2*dimension + idir] = ipos - dirstep(idir%dimension);
+      }
+    }