Force Field X

public GKTensorGlobal(GK_MULTIPOLE_ORDER multipoleOrder, int order, GKSource gkSource, double Eh,
      double Es) {
    super(order);
    this.multipoleOrder = multipoleOrder;
    this.gkSource = gkSource;

    // Load the dielectric function
    c = GKSource.cn(multipoleOrder.getOrder(), Eh, Es);
  }

  /**
   * GK Permanent multipole energy.
   *
   * @param mI PolarizableMultipole at site I.
   * @param mK PolarizableMultipole at site K.
   * @return a double.
   */
  @Override
  public double multipoleEnergy(PolarizableMultipole mI, PolarizableMultipole mK) {
    return switch (multipoleOrder) {
      default -> {
        chargeIPotentialAtK(mI, 2);
        double eK = multipoleEnergy(mK);
        chargeKPotentialAtI(mK, 2);
        double eI = multipoleEnergy(mI);
        yield c * 0.5 * (eK + eI);
      }
      case DIPOLE -> {
        dipoleIPotentialAtK(mI.dx, mI.dy, mI.dz, 2);
        double eK = multipoleEnergy(mK);
        dipoleKPotentialAtI(mK.dx, mK.dy, mK.dz, 2);
        double eI = multipoleEnergy(mI);
        yield c * 0.5 * (eK + eI);
      }
      case QUADRUPOLE -> {
        quadrupoleIPotentialAtK(mI, 2);
        double eK = multipoleEnergy(mK);
        quadrupoleKPotentialAtI(mK, 2);
        double eI = multipoleEnergy(mI);
        yield c * 0.5 * (eK + eI);
      }
    };
  }

  /**
   * GK Permanent multipole energy and gradient.
   *
   * @param mI PolarizableMultipole at site I.
   * @param mK PolarizableMultipole at site K.
   * @param Gi Coordinate gradient at site I.
   * @param Gk Coordinate gradient at site K.
   * @param Ti Torque at site I.
   * @param Tk Torque at site K.
   * @return the permanent multipole GK energy.
   */
  @Override
  public double multipoleEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
      double[] Gi, double[] Gk, double[] Ti, double[] Tk) {
    return switch (multipoleOrder) {
      default -> monopoleEnergyAndGradient(mI, mK, Gi, Gk, Ti, Tk);
      case DIPOLE -> dipoleEnergyAndGradient(mI, mK, Gi, Gk, Ti, Tk);
      case QUADRUPOLE -> quadrupoleEnergyAndGradient(mI, mK, Gi, Gk, Ti, Tk);
    };
  }

  /**
   * Permanent multipole energy and gradient using the GK monopole tensor.
   *
   * @param mI PolarizableMultipole at site I.
   * @param mK PolarizableMultipole at site K.
   * @param Gi Coordinate gradient at site I.
   * @param Gk Coordinate gradient at site K.
   * @param Ti Torque at site I.
   * @param Tk Torque at site K.
   * @return the permanent multipole GK energy.
   */
  protected double monopoleEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
      double[] Gi, double[] Gk, double[] Ti, double[] Tk) {

    // Compute the potential due to a multipole component at site I.
    chargeIPotentialAtK(mI, 3);
    double eK = multipoleEnergy(mK);
    multipoleGradient(mK, Gk);
    multipoleTorque(mK, Tk);

    // Compute the potential due to a multipole component at site K.
    chargeKPotentialAtI(mK, 3);
    double eI = multipoleEnergy(mI);
    multipoleGradient(mI, Gi);
    multipoleTorque(mI, Ti);

    double scale = c * 0.5;
    Gi[0] = scale * (Gi[0] - Gk[0]);
    Gi[1] = scale * (Gi[1] - Gk[1]);
    Gi[2] = scale * (Gi[2] - Gk[2]);
    Gk[0] = -Gi[0];
    Gk[1] = -Gi[1];
    Gk[2] = -Gi[2];

    Ti[0] = scale * Ti[0];
    Ti[1] = scale * Ti[1];
    Ti[2] = scale * Ti[2];
    Tk[0] = scale * Tk[0];
    Tk[1] = scale * Tk[1];
    Tk[2] = scale * Tk[2];

    return scale * (eK + eI);
  }

  /**
   * Permanent multipole energy and gradient using the GK dipole tensor.
   *
   * @param mI PolarizableMultipole at site I.
   * @param mK PolarizableMultipole at site K.
   * @param Gi Coordinate gradient at site I.
   * @param Gk Coordinate gradient at site K.
   * @param Ti Torque at site I.
   * @param Tk Torque at site K.
   * @return the permanent multipole GK energy.
   */
  protected double dipoleEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
      double[] Gi, double[] Gk, double[] Ti, double[] Tk) {

    // Compute the potential due to a multipole component at site I.
    dipoleIPotentialAtK(mI.dx, mI.dy, mI.dz, 3);
    double eK = multipoleEnergy(mK);
    multipoleGradient(mK, Gk);
    multipoleTorque(mK, Tk);

    // Need the torque on site I pole due to site K multipole.
    // Only torque on the site I dipole.
    multipoleKPotentialAtI(mK, 1);
    dipoleTorque(mI, Ti);

    // Compute the potential due to a multipole component at site K.
    dipoleKPotentialAtI(mK.dx, mK.dy, mK.dz, 3);
    double eI = multipoleEnergy(mI);
    multipoleGradient(mI, Gi);
    multipoleTorque(mI, Ti);

    // Need the torque on site K pole due to multipole on site I.
    // Only torque on the site K dipole.
    multipoleIPotentialAtK(mI, 1);
    dipoleTorque(mK, Tk);

    double scale = c * 0.5;
    Gi[0] = scale * (Gi[0] - Gk[0]);
    Gi[1] = scale * (Gi[1] - Gk[1]);
    Gi[2] = scale * (Gi[2] - Gk[2]);
    Gk[0] = -Gi[0];
    Gk[1] = -Gi[1];
    Gk[2] = -Gi[2];

    Ti[0] = scale * Ti[0];
    Ti[1] = scale * Ti[1];
    Ti[2] = scale * Ti[2];
    Tk[0] = scale * Tk[0];
    Tk[1] = scale * Tk[1];
    Tk[2] = scale * Tk[2];

    return scale * (eK + eI);
  }

  /**
   * Permanent multipole energy and gradient using the GK quadrupole tensor.
   *
   * @param mI PolarizableMultipole at site I.
   * @param mK PolarizableMultipole at site K.
   * @param Gi Coordinate gradient at site I.
   * @param Gk Coordinate gradient at site K.
   * @param Ti Torque at site I.
   * @param Tk Torque at site K.
   * @return the permanent multipole GK energy.
   */
  protected double quadrupoleEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
      double[] Gi, double[] Gk, double[] Ti, double[] Tk) {

    // Compute the potential due to a multipole component at site I.
    quadrupoleIPotentialAtK(mI, 3);
    double eK = multipoleEnergy(mK);
    multipoleGradient(mK, Gk);
    multipoleTorque(mK, Tk);

    // Need the torque on site I quadrupole due to site K multipole.
    multipoleKPotentialAtI(mK, 2);
    quadrupoleTorque(mI, Ti);

    // Compute the potential due to a multipole component at site K.
    quadrupoleKPotentialAtI(mK, 3);
    double eI = multipoleEnergy(mI);
    multipoleGradient(mI, Gi);
    multipoleTorque(mI, Ti);

    // Need the torque on site K quadrupole due to site I multipole.
    multipoleIPotentialAtK(mI, 2);
    quadrupoleTorque(mK, Tk);

    double scale = c * 0.5;
    Gi[0] = scale * (Gi[0] - Gk[0]);
    Gi[1] = scale * (Gi[1] - Gk[1]);
    Gi[2] = scale * (Gi[2] - Gk[2]);
    Gk[0] = -Gi[0];
    Gk[1] = -Gi[1];
    Gk[2] = -Gi[2];

    Ti[0] = scale * Ti[0];
    Ti[1] = scale * Ti[1];
    Ti[2] = scale * Ti[2];
    Tk[0] = scale * Tk[0];
    Tk[1] = scale * Tk[1];
    Tk[2] = scale * Tk[2];

    return scale * (eK + eI);
  }

  /**
   * GK Permanent multipole Born grad.
   *
   * @param mI PolarizableMultipole at site I.
   * @param mK PolarizableMultipole at site K.
   * @return a double.
   */
  public double multipoleEnergyBornGrad(PolarizableMultipole mI, PolarizableMultipole mK) {
    generateTensor();
    return multipoleEnergy(mI, mK);
  }

  /**
   * GK Polarization Energy.
   *
   * @param mI PolarizableMultipole at site I.
   * @param mK PolarizableMultipole at site K.
   * @param scaleEnergy This is ignored, since masking/scaling is not applied to GK interactions
   *     (everything is intermolecular).
   * @return a double.
   */
  @Override
  public double polarizationEnergy(PolarizableMultipole mI, PolarizableMultipole mK,
      double scaleEnergy) {
    return polarizationEnergy(mI, mK);
  }

  /**
   * GK Polarization Energy.
   *
   * @param mI PolarizableMultipole at site I.
   * @param mK PolarizableMultipole at site K.
   * @return a double.
   */
  public double polarizationEnergy(PolarizableMultipole mI, PolarizableMultipole mK) {
    return switch (multipoleOrder) {
      default -> {
        // Find the GK charge potential of site I at site K.
        chargeIPotentialAtK(mI, 1);
        // Energy of induced dipole K in the field of permanent charge I.
        double eK = polarizationEnergy(mK);
        // Find the GK charge potential of site K at site I.
        chargeKPotentialAtI(mK, 1);
        // Energy of induced dipole I in the field of permanent charge K.
        double eI = polarizationEnergy(mI);
        yield c * 0.5 * (eK + eI);
      }
      case DIPOLE -> {
        // Find the GK dipole potential of site I at site K.
        dipoleIPotentialAtK(mI.dx, mI.dy, mI.dz, 1);
        // Energy of induced dipole K in the field of permanent dipole I.
        double eK = polarizationEnergy(mK);
        // Find the GK induced dipole potential of site I at site K.
        dipoleIPotentialAtK(mI.ux, mI.uy, mI.uz, 2);
        // Energy of permanent multipole K in the field of induced dipole I.
        eK += 0.5 * multipoleEnergy(mK);
        // Find the GK dipole potential of site K at site I.
        dipoleKPotentialAtI(mK.dx, mK.dy, mK.dz, 1);
        // Energy of induced dipole I in the field of permanent dipole K.
        double eI = polarizationEnergy(mI);
        // Find the GK induced dipole potential of site K at site I.
        dipoleKPotentialAtI(mK.ux, mK.uy, mK.uz, 2);
        // Energy of permanent multipole I in the field of induced dipole K.
        eI += 0.5 * multipoleEnergy(mI);
        yield c * 0.5 * (eK + eI);
      }
      case QUADRUPOLE -> {
        // Find the GK quadrupole potential of site I at site K.
        quadrupoleIPotentialAtK(mI, 1);
        // Energy of induced dipole K in the field of permanent quadrupole I.
        double eK = polarizationEnergy(mK);
        // Find the GK quadrupole potential of site K at site I.
        quadrupoleKPotentialAtI(mK, 1);
        // Energy of induced dipole I in the field of permanent quadrupole K.
        double eI = polarizationEnergy(mI);
        yield c * 0.5 * (eK + eI);
      }
    };
  }

  /**
   * GK Polarization Energy.
   *
   * @param mI PolarizableMultipole at site I.
   * @param mK PolarizableMultipole at site K.
   * @return a double.
   */
  public double polarizationEnergyBorn(PolarizableMultipole mI, PolarizableMultipole mK) {
    return switch (multipoleOrder) {
      default -> {
        // Find the GK charge potential of site I at site K.
        chargeIPotentialAtK(mI, 1);
        // Energy of induced dipole K in the field of permanent charge I.
        double eK = polarizationEnergyS(mK);
        // Find the GK charge potential of site K at site I.
        chargeKPotentialAtI(mK, 1);
        // Energy of induced dipole I in the field of permanent charge K.
        double eI = polarizationEnergyS(mI);
        yield c * 0.5 * (eK + eI);
      }
      case DIPOLE -> {
        // Find the GK dipole potential of site I at site K.
        dipoleIPotentialAtK(mI.dx, mI.dy, mI.dz, 1);
        // Energy of induced dipole K in the field of permanent dipole I.
        double eK = polarizationEnergyS(mK);
        // Find the GK induced dipole potential of site I at site K.
        dipoleIPotentialAtK(mI.sx, mI.sy, mI.sz, 2);
        // Energy of permanent multipole K in the field of induced dipole I.
        eK += 0.5 * multipoleEnergy(mK);
        // Find the GK dipole potential of site K at site I.
        dipoleKPotentialAtI(mK.dx, mK.dy, mK.dz, 1);
        // Energy of induced dipole I in the field of permanent dipole K.
        double eI = polarizationEnergyS(mI);
        // Find the GK induced dipole potential of site K at site I.
        dipoleKPotentialAtI(mK.sx, mK.sy, mK.sz, 2);
        // Energy of permanent multipole I in the field of induced dipole K.
        eI += 0.5 * multipoleEnergy(mI);
        yield c * 0.5 * (eK + eI);
      }
      case QUADRUPOLE -> {
        // Find the GK quadrupole potential of site I at site K.
        quadrupoleIPotentialAtK(mI, 1);
        // Energy of induced dipole K in the field of permanent quadrupole I.
        double eK = polarizationEnergyS(mK);
        // Find the GK quadrupole potential of site K at site I.
        quadrupoleKPotentialAtI(mK, 1);
        // Energy of induced dipole I in the field of permanent quadrupole K.
        double eI = polarizationEnergyS(mI);
        yield c * 0.5 * (eK + eI);
      }
    };
  }

  /**
   * Polarization Energy and Gradient.
   *
   * @param mI PolarizableMultipole at site I.
   * @param mK PolarizableMultipole at site K.
   * @param inductionMask This is ignored, since masking/scaling is not applied to GK
   *     interactions (everything is intermolecular).
   * @param energyMask This is ignored, since masking/scaling is not applied to GK interactions
   *     (everything is intermolecular).
   * @param mutualMask This should be set to zero for direction polarization.
   * @param Gi an array of {@link double} objects.
   * @param Ti an array of {@link double} objects.
   * @param Tk an array of {@link double} objects.
   * @return a double.
   */
  @Override
  public double polarizationEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
      double inductionMask, double energyMask, double mutualMask, double[] Gi, double[] Ti,
      double[] Tk) {
    return switch (multipoleOrder) {
      default -> monopolePolarizationEnergyAndGradient(mI, mK, Gi);
      case DIPOLE -> dipolePolarizationEnergyAndGradient(mI, mK, mutualMask, Gi, Ti, Tk);
      case QUADRUPOLE -> quadrupolePolarizationEnergyAndGradient(mI, mK, Gi, Ti, Tk);
    };
  }

  /**
   * Monopole Polarization Energy and Gradient.
   *
   * @param mI PolarizableMultipole at site I.
   * @param mK PolarizableMultipole at site K.
   * @param Gi an array of {@link double} objects.
   * @return a double.
   */
  public double monopolePolarizationEnergyAndGradient(PolarizableMultipole mI,
      PolarizableMultipole mK, double[] Gi) {

    // Find the permanent multipole potential at site k.
    chargeIPotentialAtK(mI, 2);
    // Energy of induced dipole k in the field of multipole i.
    double eK = polarizationEnergy(mK);
    // Derivative with respect to moving atom k.
    Gi[0] = -(mK.sx * E200 + mK.sy * E110 + mK.sz * E101);
    Gi[1] = -(mK.sx * E110 + mK.sy * E020 + mK.sz * E011);
    Gi[2] = -(mK.sx * E101 + mK.sy * E011 + mK.sz * E002);

    // Find the permanent multipole potential and derivatives at site i.
    chargeKPotentialAtI(mK, 2);
    // Energy of induced dipole i in the field of multipole k.
    double eI = polarizationEnergy(mI);
    // Derivative with respect to moving atom i.
    Gi[0] += (mI.sx * E200 + mI.sy * E110 + mI.sz * E101);
    Gi[1] += (mI.sx * E110 + mI.sy * E020 + mI.sz * E011);
    Gi[2] += (mI.sx * E101 + mI.sy * E011 + mI.sz * E002);

    double scale = c * 0.5;
    Gi[0] *= scale;
    Gi[1] *= scale;
    Gi[2] *= scale;

    // Total polarization energy.
    return scale * (eI + eK);
  }

  /**
   * Dipole Polarization Energy and Gradient.
   *
   * @param mI PolarizableMultipole at site I.
   * @param mK PolarizableMultipole at site K.
   * @param mutualMask This should be set to zero for direction polarization.
   * @param Gi an array of {@link double} objects.
   * @param Ti an array of {@link double} objects.
   * @param Tk an array of {@link double} objects.
   * @return a double.
   */
  public double dipolePolarizationEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
      double mutualMask, double[] Gi, double[] Ti, double[] Tk) {

    // Find the permanent multipole potential and derivatives at site k.
    dipoleIPotentialAtK(mI.dx, mI.dy, mI.dz, 2);
    // Energy of induced dipole k in the field of multipole i.
    double eK = polarizationEnergy(mK);
    // Derivative with respect to moving atom k.
    Gi[0] = -(mK.sx * E200 + mK.sy * E110 + mK.sz * E101);
    Gi[1] = -(mK.sx * E110 + mK.sy * E020 + mK.sz * E011);
    Gi[2] = -(mK.sx * E101 + mK.sy * E011 + mK.sz * E002);

gkQuadrupole = new GKTensorGlobal(QUADRUPOLE, quadrupoleOrder, gkSource, 1.0, epsilon);
  }

  public void initPotential(double[] r, double r2, double rbi, double rbk) {
    gkSource.generateSource(POTENTIAL, QUADRUPOLE, r2, rbi, rbk);
    gkMonopole.setR(r);
    gkDipole.setR(r);
    gkQuadrupole.setR(r);
    gkMonopole.generateTensor();
    gkDipole.generateTensor();
    gkQuadrupole.generateTensor();
  }

  public void initBorn(double[] r, double r2, double rbi, double rbk) {
    gkSource.generateSource(BORN, QUADRUPOLE, r2, rbi, rbk);
    gkMonopole.setR(r);
    gkDipole.setR(r);
    gkQuadrupole.setR(r);
    gkMonopole.generateTensor();
    gkDipole.generateTensor();
    gkQuadrupole.generateTensor();
  }

  public double multipoleEnergy(PolarizableMultipole mI, PolarizableMultipole mK) {
    double em = gkMonopole.multipoleEnergy(mI, mK);
    double ed = gkDipole.multipoleEnergy(mI, mK);
    double eq = gkQuadrupole.multipoleEnergy(mI, mK);
    return em + ed + eq;
  }

  public double polarizationEnergy(PolarizableMultipole mI, PolarizableMultipole mK) {
    double emp = gkMonopole.polarizationEnergy(mI, mK);
    double edp = gkDipole.polarizationEnergy(mI, mK);
    double eqp = gkQuadrupole.polarizationEnergy(mI, mK);
    return emp + edp + eqp;
  }

  public double multipoleEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
      double[] gradI, double[] torqueI, double[] torqueK) {
    double[] gI = new double[3];
    double[] gK = new double[3];
    double[] tI = new double[3];
    double[] tK = new double[3];
    double em = gkMonopole.multipoleEnergyAndGradient(mI, mK, gI, gK, tI, tK);
    for (int j = 0; j < 3; j++) {
      gradI[j] = gI[j];
      torqueI[j] = tI[j];
      torqueK[j] = tK[j];
    }
    fill(gI, 0.0);
    fill(gK, 0.0);
    fill(tI, 0.0);
    fill(tK, 0.0);
    double ed = gkDipole.multipoleEnergyAndGradient(mI, mK, gI, gK, tI, tK);
    for (int j = 0; j < 3; j++) {
      gradI[j] += gI[j];
      torqueI[j] += tI[j];
      torqueK[j] += tK[j];
    }
    fill(gI, 0.0);
    fill(gK, 0.0);
    fill(tI, 0.0);
    fill(tK, 0.0);
    double eq = gkQuadrupole.multipoleEnergyAndGradient(mI, mK, gI, gK, tI, tK);
    for (int j = 0; j < 3; j++) {
      gradI[j] += gI[j];
      torqueI[j] += tI[j];
      torqueK[j] += tK[j];
    }

    return em + ed + eq;
  }

  public double polarizationEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
      double mutualMask,
      double[] gradI, double[] torqueI, double[] torqueK) {
    double[] gI = new double[3];
    double[] tI = new double[3];
    double[] tK = new double[3];
    double emp = gkMonopole.polarizationEnergyAndGradient(mI, mK, 1.0, 1.0, mutualMask, gI, tI, tK);
    for (int j = 0; j < 3; j++) {
      gradI[j] = gI[j];
      torqueI[j] = tI[j];
      torqueK[j] = tK[j];
    }
    fill(gI, 0.0);
    fill(tI, 0.0);
    fill(tK, 0.0);
    double edp = gkDipole.polarizationEnergyAndGradient(mI, mK, 1.0, 1.0, mutualMask, gI, tI, tK);
    for (int j = 0; j < 3; j++) {
      gradI[j] += gI[j];
      torqueI[j] += tI[j];
      torqueK[j] += tK[j];
    }
    fill(gI, 0.0);
    fill(tI, 0.0);
    fill(tK, 0.0);
    double eqp = gkQuadrupole.polarizationEnergyAndGradient(mI, mK, 1.0, 1.0, mutualMask, gI, tI,
        tK);
    for (int j = 0; j < 3; j++) {
      gradI[j] += gI[j];
      torqueI[j] += tI[j];
      torqueK[j] += tK[j];
    }

    // Sum the GK polarization interaction energy.
    return emp + edp + eqp;
  }

  public double multipoleEnergyBornGrad(PolarizableMultipole mI, PolarizableMultipole mK) {
    double db = gkMonopole.multipoleEnergyBornGrad(mI, mK);
    db += gkDipole.multipoleEnergyBornGrad(mI, mK);
    db += gkQuadrupole.multipoleEnergyBornGrad(mI, mK);
    return db;
  }

  public double polarizationEnergyBornGrad(PolarizableMultipole mI, PolarizableMultipole mK,
      boolean mutual) {
    // Compute the GK polarization Born chain-rule term.
    double db = gkMonopole.polarizationEnergyBornGrad(mI, mK);
    db += gkDipole.polarizationEnergyBornGrad(mI, mK);
    db += gkQuadrupole.polarizationEnergyBornGrad(mI, mK);
    // Add the mutual polarization contribution to Born chain-rule term.
    if (mutual) {
      db += gkDipole.mutualPolarizationEnergyBornGrad(mI, mK);
    }
    return db;
  }

}

eI += 0.5 * multipoleEnergy(mI);
    G = new double[3];
    multipoleGradient(mI, G);
    Gi[0] += G[0];
    Gi[1] += G[1];
    Gi[2] += G[2];
    multipoleTorque(mI, Ti);

    // Get the induced-induced portion of the force (Ud . dC/dX . Up).
    // This contribution does not exist for direct polarization (mutualMask == 0.0).
    if (mutualMask != 0.0) {
      // Find the potential and its derivatives at k due to induced dipole i.
      dipoleIPotentialAtK(mI.ux, mI.uy, mI.uz, 2);
      Gi[0] -= mutualMask * (mK.px * E200 + mK.py * E110 + mK.pz * E101);
      Gi[1] -= mutualMask * (mK.px * E110 + mK.py * E020 + mK.pz * E011);
      Gi[2] -= mutualMask * (mK.px * E101 + mK.py * E011 + mK.pz * E002);

      // Find the potential and its derivatives at i due to induced dipole k.
      dipoleKPotentialAtI(mK.ux, mK.uy, mK.uz, 2);
      Gi[0] += mutualMask * (mI.px * E200 + mI.py * E110 + mI.pz * E101);
      Gi[1] += mutualMask * (mI.px * E110 + mI.py * E020 + mI.pz * E011);
      Gi[2] += mutualMask * (mI.px * E101 + mI.py * E011 + mI.pz * E002);
    }

    // Total polarization energy.
    double scale = c * 0.5;
    double energy = scale * (eI + eK);
    Gi[0] *= scale;
    Gi[1] *= scale;
    Gi[2] *= scale;
    Ti[0] *= scale;
    Ti[1] *= scale;
    Ti[2] *= scale;
    Tk[0] *= scale;
    Tk[1] *= scale;
    Tk[2] *= scale;

    return energy;
  }

  /**
   * Quadrupole Polarization Energy and Gradient.
   *
   * @param mI PolarizableMultipole at site I.
   * @param mK PolarizableMultipole at site K.
   * @param Gi an array of {@link double} objects.
   * @param Ti an array of {@link double} objects.
   * @param Tk an array of {@link double} objects.
   * @return a double.
   */
  public double quadrupolePolarizationEnergyAndGradient(PolarizableMultipole mI,
      PolarizableMultipole mK, double[] Gi, double[] Ti, double[] Tk) {

    // Find the permanent multipole potential and derivatives at site k.
    quadrupoleIPotentialAtK(mI, 2);
    // Energy of induced dipole k in the field of multipole i.
    double eK = polarizationEnergy(mK);
    // Derivative with respect to moving atom k.
    Gi[0] = -(mK.sx * E200 + mK.sy * E110 + mK.sz * E101);
    Gi[1] = -(mK.sx * E110 + mK.sy * E020 + mK.sz * E011);
    Gi[2] = -(mK.sx * E101 + mK.sy * E011 + mK.sz * E002);

    // Find the permanent multipole potential and derivatives at site i.
    quadrupoleKPotentialAtI(mK, 2);
    // Energy of induced dipole i in the field of multipole k.
    double eI = polarizationEnergy(mI);
    // Derivative with respect to moving atom i.
    Gi[0] += (mI.sx * E200 + mI.sy * E110 + mI.sz * E101);
    Gi[1] += (mI.sx * E110 + mI.sy * E020 + mI.sz * E011);
    Gi[2] += (mI.sx * E101 + mI.sy * E011 + mI.sz * E002);

    double scale = c * 0.5;
    Gi[0] *= scale;
    Gi[1] *= scale;
    Gi[2] *= scale;

    // Find the potential and its derivatives at K due to the averaged induced dipole at site i.
    dipoleIPotentialAtK(scale * mI.sx, scale * mI.sy, scale * mI.sz, 2);
    quadrupoleTorque(mK, Tk);

    // Find the potential and its derivatives at I due to the averaged induced dipole at k.
    dipoleKPotentialAtI(scale * mK.sx, scale * mK.sy, scale * mK.sz, 2);
    quadrupoleTorque(mI, Ti);

    // Total polarization energy.
    return scale * (eI + eK);
  }

  /**
   * GK Direct Polarization Born grad.
   *
   * @param mI PolarizableMultipole at site I.
   * @param mK PolarizableMultipole at site K.
   * @return Partial derivative of the Polarization energy with respect to a Born grad.
   */
  public double polarizationEnergyBornGrad(PolarizableMultipole mI, PolarizableMultipole mK) {
    generateTensor();
    return 2.0 * polarizationEnergyBorn(mI, mK);
  }

  /**
   * GK Mutual Polarization Contribution to the Born grad.
   *
   * @param mI PolarizableMultipole at site I.
   * @param mK PolarizableMultipole at site K.
   * @return Mutual Polarization contribution to the partial derivative with respect to a Born grad.
   */
  public double mutualPolarizationEnergyBornGrad(PolarizableMultipole mI, PolarizableMultipole mK) {
    double db = 0.0;
    if (multipoleOrder == GK_MULTIPOLE_ORDER.DIPOLE) {
      // Find the potential and its derivatives at k due to induced dipole i.
      dipoleIPotentialAtK(mI.ux, mI.uy, mI.uz, 2);
      db = 0.5 * (mK.px * E100 + mK.py * E010 + mK.pz * E001);

      // Find the potential and its derivatives at i due to induced dipole k.
      dipoleKPotentialAtI(mK.ux, mK.uy, mK.uz, 2);
      db += 0.5 * (mI.px * E100 + mI.py * E010 + mI.pz * E001);
    }
    return c * db;
  }

  /**
   * Generate source terms for the Kirkwood version of the Challacombe et al. recursion.
   */
  @Override
  protected void source(double[] work) {
    gkSource.source(work, multipoleOrder);
  }

}

logger.info(format(" Opening %s with %d atoms\n", xphFile.getName(), numberOfAtoms));
      remarkLine = data.trim();

      // The header line is reasonable. Check for periodic box dimensions.
      br.mark(10000);
      data = br.readLine();
      if (!readPBC(data, activeMolecularAssembly)) {
        br.reset();
      }

      // Prepare to parse atom lines.
      HashMap<Integer, Integer> labelHash = new HashMap<>();
      int[] label = new int[numberOfAtoms];
      int[][] bonds = new int[numberOfAtoms][8];
      double[][] d = new double[numberOfAtoms][3];
      boolean renumber = false;
      atomList = new ArrayList<>();
      // Loop over the expected number of atoms.
      for (int i = 0; i < numberOfAtoms; i++) {
        if (!br.ready()) {
          return false;
        }
        data = br.readLine();
        if (data == null) {
          logger.warning(
              format(" Check atom %d in %s.", (i + 1), activeMolecularAssembly.getFile().getName()));
          return false;
        }
        tokens = data.trim().split(" +");
        if (tokens.length < 6) {
          logger.warning(
              format(" Check atom %d in %s.", (i + 1), activeMolecularAssembly.getFile().getName()));
          return false;
        }
        // Valid number of tokens, so try to parse this line.
        label[i] = parseInt(tokens[0]);
        // Check for valid atom numbering, or flag for re-numbering.
        if (label[i] != i + 1) {
          renumber = true;
        }
        String atomName = tokens[1];
        d[i][0] = parseDouble(tokens[2]);
        d[i][1] = parseDouble(tokens[3]);
        d[i][2] = parseDouble(tokens[4]);
        int type = parseInt(tokens[5]);
        AtomType atomType = forceField.getAtomType(Integer.toString(type));
        if (atomType == null) {
          StringBuilder message = new StringBuilder("Check atom type ");
          message.append(type).append(" for Atom ").append(i + 1);
          message.append(" in ").append(activeMolecularAssembly.getFile().getName());
          logger.warning(message.toString());
          return false;
        }
        Atom a = new Atom(i + 1, atomName, atomType, d[i]);
        atomList.add(a);
        // Bond Data
        int numberOfBonds = tokens.length - 6;
        for (int b = 0; b < 8; b++) {
          if (b < numberOfBonds) {
            int bond = parseInt(tokens[6 + b]);
            bonds[i][b] = bond;
          } else {
            bonds[i][b] = 0;
          }
        }
      }

      // Check if this is an archive.
      if (br.ready()) {
        // Read past blank lines between archive files
        data = br.readLine().trim();
        while (data.equals("") && br.ready()) {
          data = br.readLine().trim();
        }
        tokens = data.split(" +", 2);
        if (tokens.length > 0) {
          try {
            int archiveNumberOfAtoms = parseInt(tokens[0]);
            if (archiveNumberOfAtoms == numberOfAtoms) {
              setType(FileType.ARC);
            }
          } catch (NumberFormatException e) {
            //
          }
        }
      }
      // Try to renumber
      if (renumber) {
        for (int i = 0; i < numberOfAtoms; i++) {
          if (labelHash.containsKey(label[i])) {
            logger.warning(format(" Two atoms have the same index: %d.", label[i]));
            return false;
          }
          labelHash.put(label[i], i + 1);
        }
        for (int i = 0; i < numberOfAtoms; i++) {
          int j = -1;
          while (j < 3 && bonds[i][++j] > 0) {
            bonds[i][j] = labelHash.get(bonds[i][j]);
          }
        }
      }
      bondList = new ArrayList<>();

}

  /**
   * readOnto
   *
   * @param newFile a {@link File} object.
   * @param oldSystem a {@link MolecularAssembly} object.
   * @return a boolean.
   */
  public static boolean readOnto(File newFile, MolecularAssembly oldSystem) {
    if (newFile == null || !newFile.exists() || oldSystem == null) {
      return false;
    }
    try (BufferedReader br = new BufferedReader(new FileReader(newFile))) {
      String data = br.readLine();
      if (data == null) {
        return false;
      }
      String[] tokens = data.trim().split(" +");
      int num_atoms = parseInt(tokens[0]);
      if (num_atoms != oldSystem.getAtomList().size()) {
        return false;
      }

      br.mark(10000);
      data = br.readLine();
      if (!readPBC(data, oldSystem)) {
        br.reset();
      }

      double[][] d = new double[num_atoms][3];
      for (int i = 0; i < num_atoms; i++) {
        if (!br.ready()) {
          return false;
        }
        data = br.readLine();
        if (data == null) {
          logger.warning(format(" Check atom %d.", (i + 1)));
          return false;
        }
        tokens = data.trim().split(" +");
        if (tokens.length < 6) {
          logger.warning(format(" Check atom %d.", (i + 1)));
          return false;
        }
        d[i][0] = parseDouble(tokens[2]);
        d[i][1] = parseDouble(tokens[3]);
        d[i][2] = parseDouble(tokens[4]);
      }
      List<Atom> atoms = oldSystem.getAtomList();
      for (Atom a : atoms) {
        int index = a.getIndex() - 1;
        a.setXYZ(d[index]);
      }
      oldSystem.center();
      oldSystem.setFile(newFile);
      return true;
    } catch (Exception e) {
      return false;
    }
  }

  private static boolean firstTokenIsInteger(String data) {
    if (data == null) {
      return false;
    }

    // Check for a blank line.
    data = data.trim();
    if (data.equals("")) {
      return false;
    }

    // Check if the first token in an integer.
    try {
      String[] tokens = data.split(" +");
      parseInt(tokens[0]);
      return true;
    } catch (NumberFormatException e) {
      return false;
    }
  }

  /**
   * Attempt to parse the String as unit cell parameters.
   *
   * @param data The String to parse.
   * @return false if the first token in the String is an integer and true otherwise.
   */
  private static boolean readPBC(String data, MolecularAssembly activeMolecularAssembly) {
    if (firstTokenIsInteger(data)) {
      return false;
    }

    String[] tokens = data.trim().split(" +");
    if (tokens.length == 6) {
      CompositeConfiguration config = activeMolecularAssembly.getProperties();
      double a = parseDouble(tokens[0]);
      double b = parseDouble(tokens[1]);
      double c = parseDouble(tokens[2]);
      double alpha = parseDouble(tokens[3]);
      double beta = parseDouble(tokens[4]);
      double gamma = parseDouble(tokens[5]);
      config.setProperty("a-axis", a);
      config.setProperty("b-axis", b);
      config.setProperty("c-axis", c);
      config.setProperty("alpha", alpha);
      config.setProperty("beta", beta);
      config.setProperty("gamma", gamma);

      Crystal crystal = activeMolecularAssembly.getCrystal();
      if (crystal != null) {
        crystal.changeUnitCellParameters(a, b, c, alpha, beta, gamma);
      }
    }
    return true;
  }

  /** {@inheritDoc} */
  @Override
  public void closeReader() {
    if (bufferedReader != null) {
      try {
        bufferedReader.close();
      } catch (IOException ex) {
        logger.warning(format(" Exception in closing XYZ filter: %s", ex));
      }
    }
  }

  @Override
  public int countNumModels() {
    File xphFile = activeMolecularAssembly.getFile();

Atom HZ3 = (Atom) residue.getAtomNode("HZ3");
        Bond CG_CB = CG.getBond(CB);
        Bond CD_CG = CD.getBond(CG);
        Bond CE_CD = CE.getBond(CD);
        Bond NZ_CE = NZ.getBond(CE);
        Bond HB_CB = HB2.getBond(CB);
        Bond HG_CG = HG2.getBond(CG);
        Bond HD_CD = HD2.getBond(CD);
        Bond HE_CE = HE2.getBond(CE);
        Bond HZ_NZ = HZ1.getBond(NZ);
        double dCG_CB = CG_CB.bondType.distance;
        double dCD_CG = CD_CG.bondType.distance;
        double dCE_CD = CE_CD.bondType.distance;
        double dNZ_CE = NZ_CE.bondType.distance;
        double dHB_CB = HB_CB.bondType.distance;
        double dHG_CG = HG_CG.bondType.distance;
        double dHD_CD = HD_CD.bondType.distance;
        double dHE_CE = HE_CE.bondType.distance;
        double dHZ_NZ = HZ_NZ.bondType.distance;
        Angle CG_CB_CA = CG.getAngle(CB, CA);
        Angle CD_CG_CB = CD.getAngle(CG, CB);
        Angle CE_CD_CG = CE.getAngle(CD, CG);
        Angle NZ_CE_CD = NZ.getAngle(CE, CD);
        Angle HB_CB_CA = HB2.getAngle(CB, CA);
        Angle HG_CG_CB = HG2.getAngle(CG, CB);
        Angle HD_CD_CG = HD2.getAngle(CD, CG);
        Angle HE_CE_CD = HE2.getAngle(CE, CD);
        Angle HZ_NZ_CE = HZ1.getAngle(NZ, CE);
        double dCG_CB_CA = CG_CB_CA.angleType.angle[CG_CB_CA.nh];
        double dCD_CG_CB = CD_CG_CB.angleType.angle[CD_CG_CB.nh];
        double dCE_CD_CG = CE_CD_CG.angleType.angle[CE_CD_CG.nh];
        double dNZ_CE_CD = NZ_CE_CD.angleType.angle[NZ_CE_CD.nh];
        double dHB_CB_CA = HB_CB_CA.angleType.angle[HB_CB_CA.nh];
        double dHG_CG_CB = HG_CG_CB.angleType.angle[HG_CG_CB.nh];
        double dHD_CD_CG = HD_CD_CG.angleType.angle[HD_CD_CG.nh];
        double dHE_CE_CD = HE_CE_CD.angleType.angle[HE_CE_CD.nh];
        double dHZ_NZ_CE = HZ_NZ_CE.angleType.angle[HZ_NZ_CE.nh];
        intxyz(CG, CB, dCG_CB, CA, dCG_CB_CA, N, rotamer.chi1, 0);
        intxyz(CD, CG, dCD_CG, CB, dCD_CG_CB, CA, rotamer.chi2, 0);
        intxyz(CE, CD, dCE_CD, CG, dCE_CD_CG, CB, rotamer.chi3, 0);
        intxyz(NZ, CE, dNZ_CE, CD, dNZ_CE_CD, CG, rotamer.chi4, 0);
        intxyz(HB2, CB, dHB_CB, CA, dHB_CB_CA, CG, 109.4, 1);
        intxyz(HB3, CB, dHB_CB, CA, dHB_CB_CA, CG, 109.4, -1);
        intxyz(HG2, CG, dHG_CG, CB, dHG_CG_CB, CD, 109.4, 1);
        intxyz(HG3, CG, dHG_CG, CB, dHG_CG_CB, CD, 109.4, -1);
        intxyz(HD2, CD, dHD_CD, CG, dHD_CD_CG, CE, 109.4, 1);
        intxyz(HD3, CD, dHD_CD, CG, dHD_CD_CG, CE, 109.4, -1);
        intxyz(HE2, CE, dHE_CE, CD, dHE_CE_CD, NZ, 108.8, 1);
        intxyz(HE3, CE, dHE_CE, CD, dHE_CE_CD, NZ, 108.8, -1);
        intxyz(HZ1, NZ, dHZ_NZ, CE, dHZ_NZ_CE, CD, 180.0, 0);
        intxyz(HZ2, NZ, dHZ_NZ, CE, dHZ_NZ_CE, HZ1, 109.5, 1);

long last)
            throws Exception;

    /**
     * Send additional output data associated with a task. Called by a worker
     * thread. The task is denoted by the given chunk of loop iterations. The
     * output data must be sent using the given communicator, to the given
     * master process rank, with the given message tag.
     * <P>
     * The <code>sendTaskOutput()</code> method may be overridden in a subclass. If
     * not overridden, the <code>sendTaskOutput()</code> method does nothing.
     *
     * @param range Chunk of loop iterations.
     * @param comm Communicator.
     * @param mRank Master process rank.
     * @param tag Message tag.
     * @exception IOException Thrown if an I/O error occurred.
     * @throws java.io.IOException if any.
     */
    public void sendTaskOutput(LongRange range,
            Comm comm,
            int mRank,
            int tag)
            throws IOException {
    }

    /**
     * Receive additional output data associated with a task. Called by the
     * master thread. The task is denoted by the given chunk of loop iterations.
     * The output data must be received using the given communicator, from the
     * given worker process rank, with the given message tag.
     * <P>
     * The <code>receiveTaskOutput()</code> method may be overridden in a subclass.
     * If not overridden, the <code>receiveTaskOutput()</code> method does nothing.
     *
     * @param range Chunk of loop iterations.
     * @param comm Communicator.
     * @param wRank Worker process rank.
     * @param tag Message tag.
     * @exception IOException Thrown if an I/O error occurred.
     * @throws java.io.IOException if any.
     */
    public void receiveTaskOutput(LongRange range,
            Comm comm,
            int wRank,
            int tag)
            throws IOException {
    }

    /**
     * Perform per-thread finalization actions after finishing the loop
     * iterations. Called by a worker thread.
     * <P>
     * The <code>finish()</code> method may be overridden in a subclass. If not
     * overridden, the <code>finish()</code> method does nothing.
     *
     * @exception Exception The <code>finish()</code> method may throw any
     * exception.
     * @throws java.lang.Exception if any.
     */
    public void finish()
            throws Exception {
    }

    /**
     * Returns the tag offset for this worker for loop. Each message between the
     * master and worker threads is sent with a message tag equal to
     * <I>W</I>+<I>T</I>, where <I>W</I> is the worker index and <I>T</I> is the
     * tag offset.
     * <P>
     * The <code>tagOffset()</code> method may be overridden in a subclass. If not
     * overridden, the <code>tagOffset()</code> returns a default tag offset of
     * <code>Integer.MIN_VALUE</code>.
     *
     * @return Tag offset.
     */
    public int tagOffset() {
        return Integer.MIN_VALUE;
    }

// Hidden operations.
    /**
     * Execute this worker for loop in the master thread.
     *
     * @param range Loop index range.
     *
     * @exception IOException Thrown if an I/O error occurred.
     */
    void masterExecute(LongRange range)
            throws IOException {
        LongSchedule sch = schedule();
        if (sch.isFixedSchedule()) {
            masterExecuteFixed(range, sch);
        } else {
            masterExecuteNonFixed(range, sch);
        }
    }

    /**
     * Execute this worker for loop in the master thread with a fixed schedule.
     *
     * @param range Loop index range.
     * @param sch Schedule.
     *
     * @exception IOException Thrown if an I/O error occurred.
     */
    void masterExecuteFixed(LongRange range,
            LongSchedule sch)
            throws IOException {
        int count = myTeam.count;
        Comm comm = myTeam.comm;

        // Send additional task input to each worker.
        sch.start(count, range);
        for (int w = 0; w < count; ++w) {
            LongRange chunk = sch.next(w);
            if (chunk != null) {
                sendTaskInput(chunk, comm, myTeam.workerRank(w), tagFor(w));
            }
        }

        // Receive additional task output from each worker.
        sch.start(count, range);
        for (int w = 0; w < count; ++w) {
            LongRange chunk = sch.next(w);
            if (chunk != null) {
                receiveTaskOutput(chunk, comm, myTeam.workerRank(w), tagFor(w));
            }
        }
    }

    /**
     * Execute this worker for loop in the master thread with a non-fixed
     * schedule.
     *
     * @param range Loop index range.
     * @param sch Schedule.
     *
     * @exception IOException Thrown if an I/O error occurred.
     */
    void masterExecuteNonFixed(LongRange range,
            LongSchedule sch)
            throws IOException {
        int count = myTeam.count;
        sch.start(count, range);
        int remaining = count;
        ObjectItemBuf<LongRange> buf = ObjectBuf.buffer();
        Range tagRange = new Range(tagFor(0), tagFor(count - 1));
        Comm comm = myTeam.comm;

        // Send initial task to each worker.
        for (int w = 0; w < count; ++w) {
            LongRange chunk = sch.next(w);
            buf.item = chunk;
            buf.reset();
            int r = myTeam.workerRank(w);
            int tag = tagFor(w);
            comm.send(r, tag, buf);
            if (chunk == null) {
                --remaining;
            } else {
                sendTaskInput(chunk, comm, r, tag);
            }
        }

        // Repeatedly receive a response from a worker and send next task to
        // that worker.
        while (remaining > 0) {
            CommStatus status = comm.receive(null, tagRange, buf);
            LongRange chunk = buf.item;
            int r = status.fromRank;
            int tag = status.tag;
            int w = workerFor(tag);
            receiveTaskOutput(chunk, comm, r, tag);
            chunk = sch.next(w);
            buf.item = chunk;
            buf.reset();
            comm.send(r, tag, buf);
            if (chunk == null) {
                --remaining;
            } else {
                sendTaskInput(chunk, comm, r, tag);
            }
        }
    }

    /**
     * Execute this worker for loop in a worker thread.
     *
     * @param w Worker index.
     * @param range Loop index range.
     *
     * @exception Exception This method may throw any exception.
     */
    void workerExecute(int w,
            LongRange range)
            throws Exception {
        LongSchedule sch = schedule();
        if (sch.isFixedSchedule()) {
            sch.start(myTeam.count, range);
            workerExecuteFixed(sch.next(w), w);
        } else {
            workerExecuteNonFixed(w);
        }
    }

    /**
     * Execute this worker for loop in a worker thread using a fixed schedule.
     *
     * @param range Chunk of loop iterations.
     * @param w Worker index.
     *
     * @exception Exception This method may throw any exception.
     */
    void workerExecuteFixed(LongRange range,
            int w)
            throws Exception {
        start();
        if (range != null) {
            Comm comm = myTeam.comm;
            int r = myTeam.masterRank();
            int tag = tagFor(w);
            receiveTaskInput(range, comm, r, tag);
            run(range.lb(), range.ub());

"P 42/m -3 2/n",
              CUBIC, CUBIC_LATTICE,
              LM3M,
              new ASULimit[] {ASULimit.LT, ASULimit.LT, ASULimit.LT},
              new double[] {-1.0, -1.0, -1.0},
              new SymOp(SymOp.Rot_X_Y_Z, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_mX_mY_Z, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_mX_Y_mZ, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_X_mY_mZ, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_Z_X_Y, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_Z_mX_mY, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_mZ_mX_Y, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_mZ_X_mY, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_Y_Z_X, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_mY_Z_mX, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_Y_mZ_mX, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_mY_mZ_X, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_Y_X_mZ, SymOp.Tr_12_12_12),
              new SymOp(SymOp.Rot_mY_mX_mZ, SymOp.Tr_12_12_12),
              new SymOp(SymOp.Rot_Y_mX_Z, SymOp.Tr_12_12_12),
              new SymOp(SymOp.Rot_mY_X_Z, SymOp.Tr_12_12_12),
              new SymOp(SymOp.Rot_X_Z_mY, SymOp.Tr_12_12_12),
              new SymOp(SymOp.Rot_mX_Z_Y, SymOp.Tr_12_12_12),
              new SymOp(SymOp.Rot_mX_mZ_mY, SymOp.Tr_12_12_12),
              new SymOp(SymOp.Rot_X_mZ_Y, SymOp.Tr_12_12_12),
              new SymOp(SymOp.Rot_Z_Y_mX, SymOp.Tr_12_12_12),
              new SymOp(SymOp.Rot_Z_mY_X, SymOp.Tr_12_12_12),
              new SymOp(SymOp.Rot_mZ_Y_X, SymOp.Tr_12_12_12),
              new SymOp(SymOp.Rot_mZ_mY_mX, SymOp.Tr_12_12_12),
              new SymOp(SymOp.Rot_mX_mY_mZ, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_X_Y_mZ, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_X_mY_Z, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_mX_Y_Z, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_mZ_mX_mY, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_mZ_X_Y, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_Z_X_mY, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_Z_mX_Y, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_mY_mZ_mX, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_Y_mZ_X, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_mY_Z_X, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_Y_Z_mX, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_mY_mX_Z, SymOp.Tr_12_12_12),
              new SymOp(SymOp.Rot_Y_X_Z, SymOp.Tr_12_12_12),
              new SymOp(SymOp.Rot_mY_X_mZ, SymOp.Tr_12_12_12),
              new SymOp(SymOp.Rot_Y_mX_mZ, SymOp.Tr_12_12_12),
              new SymOp(SymOp.Rot_mX_mZ_Y, SymOp.Tr_12_12_12),
              new SymOp(SymOp.Rot_X_mZ_mY, SymOp.Tr_12_12_12),
              new SymOp(SymOp.Rot_X_Z_Y, SymOp.Tr_12_12_12),
              new SymOp(SymOp.Rot_mX_Z_mY, SymOp.Tr_12_12_12),
              new SymOp(SymOp.Rot_mZ_mY_X, SymOp.Tr_12_12_12),
              new SymOp(SymOp.Rot_mZ_Y_mX, SymOp.Tr_12_12_12),
              new SymOp(SymOp.Rot_Z_mY_mX, SymOp.Tr_12_12_12),
              new SymOp(SymOp.Rot_Z_Y_X, SymOp.Tr_12_12_12));

new double[] {f12, f14, f14},
              new SymOp(SymOp.Rot_X_Y_Z, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_mX_mY_Z, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_mX_Y_mZ, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_X_mY_mZ, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_Z_X_Y, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_Z_mX_mY, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_mZ_mX_Y, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_mZ_X_mY, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_Y_Z_X, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_mY_Z_mX, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_Y_mZ_mX, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_mY_mZ_X, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_Y_X_mZ, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_mY_mX_mZ, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_Y_mX_Z, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_mY_X_Z, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_X_Z_mY, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_mX_Z_Y, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_mX_mZ_mY, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_X_mZ_Y, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_Z_Y_mX, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_Z_mY_X, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_mZ_Y_X, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_mZ_mY_mX, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_mX_mY_mZ, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_X_Y_mZ, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_X_mY_Z, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_mX_Y_Z, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_mZ_mX_mY, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_mZ_X_Y, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_Z_X_mY, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_Z_mX_Y, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_mY_mZ_mX, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_Y_mZ_X, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_mY_Z_X, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_Y_Z_mX, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_mY_mX_Z, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_Y_X_Z, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_mY_X_mZ, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_Y_mX_mZ, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_mX_mZ_Y, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_X_mZ_mY, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_X_Z_Y, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_mX_Z_mY, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_mZ_mY_X, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_mZ_Y_mX, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_Z_mY_mX, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_Z_Y_X, SymOp.Tr_0_0_0),
              new SymOp(SymOp.Rot_X_Y_Z, SymOp.Tr_0_12_12),

AtomicRowLoop(RowRegion region) {
      super(region.getNatoms(), region.getNsymm(), region);
      grid = region.getGrid();
      optLocal = new int[fftZ * fftY];
    }

    public void buildList(int iSymm, int iAtom, int lb, int ub) {
      if (!atoms[iAtom].getUse() && !nativeEnvironmentApproximation) {
        return;
      }
      xyz[0] = coordinates[iSymm][0][iAtom];
      xyz[1] = coordinates[iSymm][1][iAtom];
      xyz[2] = coordinates[iSymm][2][iAtom];
      crystal.toFractionalCoordinates(xyz, uvw);
      final int frad =
          min(aRadGrid, (int) floor(atoms[iAtom].getFormFactorWidth() * fftX / crystal.a) + 1);

      final double frz = fftZ * uvw[2];
      final int ifrz = (int) frz;
      final int ifrzu = ifrz + frad;
      final int ifrzl = ifrz - frad;

      final double fry = fftY * uvw[1];
      final int ifry = (int) fry;
      final int ifryu = ifry + frad;
      final int ifryl = ifry - frad;

      // Loop over allowed z coordinates for this Loop
      // Check if the current atom is close enough
      // If so, add to list.
      int buff = bufferSize;

      int lbZ = rowRegion.zFromRowIndex(lb);
      int ubZ = rowRegion.zFromRowIndex(ub);

      for (int iz = ifrzl - buff; iz <= ifrzu + buff; iz++) {
        int giz = mod(iz, fftZ);
        if (lbZ > giz || giz > ubZ) {
          continue;
        }
        int rowLB = rowRegion.rowIndexForYZ(mod(ifryl - buff, fftY), giz);
        int rowUB = rowRegion.rowIndexForYZ(mod(ifryu + buff, fftY), giz);
        if (lb >= rowLB || rowUB <= ub) {
          buildListA.add(iAtom);
          buildListS.add(iSymm);
          break;
        }
      }
    }

    @Override
    public boolean checkList(int[][][] zyAtListBuild, int buff) {
      for (int iSymm = 0; iSymm < nSymm; iSymm++) {
        for (int iAtom = 0; iAtom < nAtoms; iAtom++) {
          if (rowRegion.select[iSymm][iAtom]) {
            if (!atoms[iAtom].getUse() && !nativeEnvironmentApproximation) {
              continue;
            }
            xyz[0] = coordinates[iSymm][0][iAtom];
            xyz[1] = coordinates[iSymm][1][iAtom];
            xyz[2] = coordinates[iSymm][2][iAtom];
            crystal.toFractionalCoordinates(xyz, uvw);
            final double frz = fftZ * uvw[2];
            final int ifrz = (int) frz;
            final int previousZ = zyAtListBuild[iSymm][iAtom][0];

            final double fry = fftY * uvw[1];
            final int ifry = (int) fry;
            final int previousY = zyAtListBuild[iSymm][iAtom][1];
            if (abs(ifrz - previousZ) >= buff || abs(ifry - previousY) >= buff) {
              return true;
            }
          }
        }
      }
      return false;
    }

    @Override
    public void finish() {

gqyz[1] = yr * zr * a[2][0];

      // Unweighted reaction potential gradient tensor.
      gc[2] = xr * a[0][1];
      gc[3] = yr * a[0][1];
      gc[4] = zr * a[0][1];
      gux[2] = a[1][0] + xr2 * a[1][1];
      gux[3] = xr * yr * a[1][1];
      gux[4] = xr * zr * a[1][1];
      guy[2] = gux[3];
      guy[3] = a[1][0] + yr2 * a[1][1];
      guy[4] = yr * zr * a[1][1];
      guz[2] = gux[4];
      guz[3] = guy[4];
      guz[4] = a[1][0] + zr2 * a[1][1];
      gqxx[2] = xr * (2.0 * a[2][0] + xr2 * a[2][1]);
      gqxx[3] = yr * xr2 * a[2][1];
      gqxx[4] = zr * xr2 * a[2][1];
      gqyy[2] = xr * yr2 * a[2][1];
      gqyy[3] = yr * (2.0 * a[2][0] + yr2 * a[2][1]);
      gqyy[4] = zr * yr2 * a[2][1];
      gqzz[2] = xr * zr2 * a[2][1];
      gqzz[3] = yr * zr2 * a[2][1];
      gqzz[4] = zr * (2.0 * a[2][0] + zr2 * a[2][1]);
      gqxy[2] = yr * (a[2][0] + xr2 * a[2][1]);
      gqxy[3] = xr * (a[2][0] + yr2 * a[2][1]);
      gqxy[4] = zr * xr * yr * a[2][1];
      gqxz[2] = zr * (a[2][0] + xr2 * a[2][1]);
      gqxz[3] = gqxy[4];
      gqxz[4] = xr * (a[2][0] + zr2 * a[2][1]);
      gqyz[2] = gqxy[4];
      gqyz[3] = zr * (a[2][0] + yr2 * a[2][1]);
      gqyz[4] = yr * (a[2][0] + zr2 * a[2][1]);

int last)
            throws Exception;

    /**
     * Send additional output data associated with a task. Called by a worker
     * thread. The task is denoted by the given chunk of loop iterations. The
     * output data must be sent using the given communicator, to the given
     * master process rank, with the given message tag.
     * <P>
     * The <code>sendTaskOutput()</code> method may be overridden in a subclass. If
     * not overridden, the <code>sendTaskOutput()</code> method does nothing.
     *
     * @param range Chunk of loop iterations.
     * @param comm Communicator.
     * @param mRank Master process rank.
     * @param tag Message tag.
     * @exception IOException Thrown if an I/O error occurred.
     * @throws java.io.IOException if any.
     */
    public void sendTaskOutput(Range range,
            Comm comm,
            int mRank,
            int tag)
            throws IOException {
    }

    /**
     * Receive additional output data associated with a task. Called by the
     * master thread. The task is denoted by the given chunk of loop iterations.
     * The output data must be received using the given communicator, from the
     * given worker process rank, with the given message tag.
     * <P>
     * The <code>receiveTaskOutput()</code> method may be overridden in a subclass.
     * If not overridden, the <code>receiveTaskOutput()</code> method does nothing.
     *
     * @param range Chunk of loop iterations.
     * @param comm Communicator.
     * @param wRank Worker process rank.
     * @param tag Message tag.
     * @exception IOException Thrown if an I/O error occurred.
     * @throws java.io.IOException if any.
     */
    public void receiveTaskOutput(Range range,
            Comm comm,
            int wRank,
            int tag)
            throws IOException {
    }

    /**
     * Perform per-thread finalization actions after finishing the loop
     * iterations. Called by a worker thread.
     * <P>
     * The <code>finish()</code> method may be overridden in a subclass. If not
     * overridden, the <code>finish()</code> method does nothing.
     *
     * @exception Exception The <code>finish()</code> method may throw any
     * exception.
     * @throws java.lang.Exception if any.
     */
    public void finish()
            throws Exception {
    }

    /**
     * Returns the tag offset for this worker for loop. Each message between the
     * master and worker threads is sent with a message tag equal to
     * <I>W</I>+<I>T</I>, where <I>W</I> is the worker index and <I>T</I> is the
     * tag offset.
     * <P>
     * The <code>tagOffset()</code> method may be overridden in a subclass. If not
     * overridden, the <code>tagOffset()</code> returns a default tag offset of
     * <code>Integer.MIN_VALUE</code>.
     *
     * @return Tag offset.
     */
    public int tagOffset() {
        return Integer.MIN_VALUE;
    }

// Hidden operations.
    /**
     * Execute this worker for loop in the master thread.
     *
     * @param range Loop index range.
     *
     * @exception IOException Thrown if an I/O error occurred.
     */
    void masterExecute(Range range)
            throws IOException {
        IntegerSchedule sch = schedule();
        if (sch.isFixedSchedule()) {
            masterExecuteFixed(range, sch);
        } else {
            masterExecuteNonFixed(range, sch);
        }
    }

    /**
     * Execute this worker for loop in the master thread with a fixed schedule.
     *
     * @param range Loop index range.
     * @param sch Schedule.
     *
     * @exception IOException Thrown if an I/O error occurred.
     */
    void masterExecuteFixed(Range range,
            IntegerSchedule sch)
            throws IOException {
        int count = myTeam.count;
        Comm comm = myTeam.comm;

        // Send additional task input to each worker.
        sch.start(count, range);
        for (int w = 0; w < count; ++w) {
            Range chunk = sch.next(w);
            if (chunk != null) {
                sendTaskInput(chunk, comm, myTeam.workerRank(w), tagFor(w));
            }
        }

        // Receive additional task output from each worker.
        sch.start(count, range);
        for (int w = 0; w < count; ++w) {
            Range chunk = sch.next(w);
            if (chunk != null) {
                receiveTaskOutput(chunk, comm, myTeam.workerRank(w), tagFor(w));
            }
        }
    }

    /**
     * Execute this worker for loop in the master thread with a non-fixed
     * schedule.
     *
     * @param range Loop index range.
     * @param sch Schedule.
     *
     * @exception IOException Thrown if an I/O error occurred.
     */
    void masterExecuteNonFixed(Range range,
            IntegerSchedule sch)
            throws IOException {
        int count = myTeam.count;
        sch.start(count, range);
        int remaining = count;
        ObjectItemBuf<Range> buf = ObjectBuf.buffer();
        Range tagRange = new Range(tagFor(0), tagFor(count - 1));
        Comm comm = myTeam.comm;

        // Send initial task to each worker.
        for (int w = 0; w < count; ++w) {
            Range chunk = sch.next(w);
            buf.item = chunk;
            buf.reset();
            int r = myTeam.workerRank(w);
            int tag = tagFor(w);
            comm.send(r, tag, buf);
            if (chunk == null) {
                --remaining;
            } else {
                sendTaskInput(chunk, comm, r, tag);
            }
        }

        // Repeatedly receive a response from a worker and send next task to
        // that worker.
        while (remaining > 0) {
            CommStatus status = comm.receive(null, tagRange, buf);
            Range chunk = buf.item;
            int r = status.fromRank;
            int tag = status.tag;
            int w = workerFor(tag);
            receiveTaskOutput(chunk, comm, r, tag);
            chunk = sch.next(w);
            buf.item = chunk;
            buf.reset();
            comm.send(r, tag, buf);
            if (chunk == null) {
                --remaining;
            } else {
                sendTaskInput(chunk, comm, r, tag);
            }
        }
    }

    /**
     * Execute this worker for loop in a worker thread.
     *
     * @param w Worker index.
     * @param range Loop index range.
     *
     * @exception Exception This method may throw any exception.
     */
    void workerExecute(int w,