private void addDataD(int nGaussians) {
    //for now just assumes 6 orbitals in the order XX YY ZZ XY XZ YZ
    if (atomCoordBohr[atomIndex] == null) {
      moCoeff += 6;
      return;
    }
    if (doDebug)
      dumpInfo(nGaussians, "XXYYZZXYXZYZ");
    setMinMax(nGaussians);
    int moCoeff0 = moCoeff;
    for (int ig = 0; ig < nGaussians; ig++) {
      moCoeff = moCoeff0;
      float alpha = gaussians[gaussianPtr + ig][0];
      float c1 = gaussians[gaussianPtr + ig][1];
      // xx|yy|zz: (2048 alpha^7/9pi^3)^0.25 [xx|yy|zz]exp(-alpha r^2)
      // xy|xz|yz: (2048 alpha^7/pi^3)^0.25 [xy|xz|yz]exp(-alpha r^2)
      float a = c1 * (float) Math.pow(alpha, 1.75) * 2.8508219178923f;
      float axx = a / ROOT3 * moCoefficients[moCoeff++];
      float ayy = a / ROOT3 * moCoefficients[moCoeff++];
      float azz = a / ROOT3 * moCoefficients[moCoeff++];
      float axy = a * moCoefficients[moCoeff++];
      float axz = a * moCoefficients[moCoeff++];
      float ayz = a * moCoefficients[moCoeff++];
      setCE(alpha, 0, axx, ayy, azz);

      for (int i = xMax; --i >= xMin;) {
        DXY[i] = axy * X[i];
        DXZ[i] = axz * X[i];
      }
      for (int i = yMax; --i >= yMin;) {
        DYZ[i] = ayz * Y[i];
      }
      for (int ix = xMax; --ix >= xMin;) {
        float axx_x2 = CX[ix] * X[ix];
        float axy_x = DXY[ix];
        float axz_x = DXZ[ix];
        for (int iy = yMax; --iy >= yMin;) {
          float axx_x2__ayy_y2__axy_xy = axx_x2 + (CY[iy] + axy_x) * Y[iy];
          float axz_x__ayz_y = axz_x + DYZ[iy];
          for (int iz = zMax; --iz >= zMin;)
            voxelData[ix][iy][iz] += (axx_x2__ayy_y2__axy_xy + (CZ[iz] + axz_x__ayz_y)
                * Z[iz])
                * EX[ix] * EY[iy] * EZ[iz];
          // giving (axx_x2 + ayy_y2 + azz_z2 + axy_xy + axz_xz + ayz_yz)e^-br2; 
        }
      }
    }
  }