Fork me on GitHub

CPD Results

The following document contains the results of PMD's CPD 6.29.0.

Duplications

File Project Line
ffx/potential/nonbonded/implicit/GKEnergyRegion.java Potential 998
ffx/potential/nonbonded/implicit/GKEnergyRegionQI.java Potential 932
return e + ei;
    }

    private void gradientTensors() {
      // Born radii gradients of unweighted reaction potential tensor.
      gc[21] = b[0][0];
      gux[21] = xr * b[1][0];
      guy[21] = yr * b[1][0];
      guz[21] = zr * b[1][0];
      gqxx[21] = xr2 * b[2][0];
      gqyy[21] = yr2 * b[2][0];
      gqzz[21] = zr2 * b[2][0];
      gqxy[21] = xr * yr * b[2][0];
      gqxz[21] = xr * zr * b[2][0];
      gqyz[21] = yr * zr * b[2][0];

      // Born gradients of the unweighted reaction potential gradient tensor
      gc[22] = xr * b[0][1];
      gc[23] = yr * b[0][1];
      gc[24] = zr * b[0][1];
      gux[22] = b[1][0] + xr2 * b[1][1];
      gux[23] = xr * yr * b[1][1];
      gux[24] = xr * zr * b[1][1];
      guy[22] = gux[23];
      guy[23] = b[1][0] + yr2 * b[1][1];
      guy[24] = yr * zr * b[1][1];
      guz[22] = gux[24];
      guz[23] = guy[24];
      guz[24] = b[1][0] + zr2 * b[1][1];
      gqxx[22] = xr * (2.0 * b[2][0] + xr2 * b[2][1]);
      gqxx[23] = yr * xr2 * b[2][1];
      gqxx[24] = zr * xr2 * b[2][1];
      gqyy[22] = xr * yr2 * b[2][1];
      gqyy[23] = yr * (2.0 * b[2][0] + yr2 * b[2][1]);
      gqyy[24] = zr * yr2 * b[2][1];
      gqzz[22] = xr * zr2 * b[2][1];
      gqzz[23] = yr * zr2 * b[2][1];
      gqzz[24] = zr * (2.0 * b[2][0] + zr2 * b[2][1]);
      gqxy[22] = yr * (b[2][0] + xr2 * b[2][1]);
      gqxy[23] = xr * (b[2][0] + yr2 * b[2][1]);
      gqxy[24] = zr * xr * yr * b[2][1];
      gqxz[22] = zr * (b[2][0] + xr2 * b[2][1]);
      gqxz[23] = gqxy[24];
      gqxz[24] = xr * (b[2][0] + zr2 * b[2][1]);
      gqyz[22] = gqxy[24];
      gqyz[23] = zr * (b[2][0] + yr2 * b[2][1]);
      gqyz[24] = yr * (b[2][0] + zr2 * b[2][1]);

      // Born radii derivatives of the unweighted 2nd reaction potential gradient tensor.
      gc[25] = b[0][1] + xr2 * b[0][2];
      gc[26] = xr * yr * b[0][2];
      gc[27] = xr * zr * b[0][2];
      gc[28] = b[0][1] + yr2 * b[0][2];
      gc[29] = yr * zr * b[0][2];
      gc[30] = b[0][1] + zr2 * b[0][2];
      gux[25] = xr * (3.0 * b[1][1] + xr2 * b[1][2]);
      gux[26] = yr * (b[1][1] + xr2 * b[1][2]);
      gux[27] = zr * (b[1][1] + xr2 * b[1][2]);
      gux[28] = xr * (b[1][1] + yr2 * b[1][2]);
      gux[29] = zr * xr * yr * b[1][2];
      gux[30] = xr * (b[1][1] + zr2 * b[1][2]);
      guy[25] = yr * (b[1][1] + xr2 * b[1][2]);
      guy[26] = xr * (b[1][1] + yr2 * b[1][2]);
      guy[27] = gux[29];
      guy[28] = yr * (3.0 * b[1][1] + yr2 * b[1][2]);
      guy[29] = zr * (b[1][1] + yr2 * b[1][2]);
      guy[30] = yr * (b[1][1] + zr2 * b[1][2]);
      guz[25] = zr * (b[1][1] + xr2 * b[1][2]);
      guz[26] = gux[29];
      guz[27] = xr * (b[1][1] + zr2 * b[1][2]);
      guz[28] = zr * (b[1][1] + yr2 * b[1][2]);
      guz[29] = yr * (b[1][1] + zr2 * b[1][2]);
      guz[30] = zr * (3.0 * b[1][1] + zr2 * b[1][2]);
      gqxx[25] = 2.0 * b[2][0] + xr2 * (5.0 * b[2][1] + xr2 * b[2][2]);
      gqxx[26] = yr * xr * (2.0 * b[2][1] + xr2 * b[2][2]);
      gqxx[27] = zr * xr * (2.0 * b[2][1] + xr2 * b[2][2]);
      gqxx[28] = xr2 * (b[2][1] + yr2 * b[2][2]);
      gqxx[29] = zr * yr * xr2 * b[2][2];
      gqxx[30] = xr2 * (b[2][1] + zr2 * b[2][2]);
      gqyy[25] = yr2 * (b[2][1] + xr2 * b[2][2]);
      gqyy[26] = xr * yr * (2.0 * b[2][1] + yr2 * b[2][2]);
      gqyy[27] = xr * zr * yr2 * b[2][2];
      gqyy[28] = 2.0 * b[2][0] + yr2 * (5.0 * b[2][1] + yr2 * b[2][2]);
      gqyy[29] = yr * zr * (2.0 * b[2][1] + yr2 * b[2][2]);
      gqyy[30] = yr2 * (b[2][1] + zr2 * b[2][2]);
      gqzz[25] = zr2 * (b[2][1] + xr2 * b[2][2]);
      gqzz[26] = xr * yr * zr2 * b[2][2];
      gqzz[27] = xr * zr * (2.0 * b[2][1] + zr2 * b[2][2]);
      gqzz[28] = zr2 * (b[2][1] + yr2 * b[2][2]);
      gqzz[29] = yr * zr * (2.0 * b[2][1] + zr2 * b[2][2]);
      gqzz[30] = 2.0 * b[2][0] + zr2 * (5.0 * b[2][1] + zr2 * b[2][2]);
      gqxy[25] = xr * yr * (3.0 * b[2][1] + xr2 * b[2][2]);
      gqxy[26] = b[2][0] + (xr2 + yr2) * b[2][1] + xr2 * yr2 * b[2][2];
      gqxy[27] = zr * yr * (b[2][1] + xr2 * b[2][2]);
      gqxy[28] = xr * yr * (3.0 * b[2][1] + yr2 * b[2][2]);
      gqxy[29] = zr * xr * (b[2][1] + yr2 * b[2][2]);
      gqxy[30] = xr * yr * (b[2][1] + zr2 * b[2][2]);
      gqxz[25] = xr * zr * (3.0 * b[2][1] + xr2 * b[2][2]);
      gqxz[26] = yr * zr * (b[2][1] + xr2 * b[2][2]);
      gqxz[27] = b[2][0] + (xr2 + zr2) * b[2][1] + xr2 * zr2 * b[2][2];
      gqxz[28] = xr * zr * (b[2][1] + yr2 * b[2][2]);
      gqxz[29] = xr * yr * (b[2][1] + zr2 * b[2][2]);
      gqxz[30] = xr * zr * (3.0 * b[2][1] + zr2 * b[2][2]);
      gqyz[25] = zr * yr * (b[2][1] + xr2 * b[2][2]);
      gqyz[26] = xr * zr * (b[2][1] + yr2 * b[2][2]);
      gqyz[27] = xr * yr * (b[2][1] + zr2 * b[2][2]);
      gqyz[28] = yr * zr * (3.0 * b[2][1] + yr2 * b[2][2]);
      gqyz[29] = b[2][0] + (yr2 + zr2) * b[2][1] + yr2 * zr2 * b[2][2];
      gqyz[30] = yr * zr * (3.0 * b[2][1] + zr2 * b[2][2]);

      // Unweighted 3rd reaction potential gradient tensor.
      gc[11] = xr * (3.0 * a[0][2] + xr2 * a[0][3]);
      gc[12] = yr * (a[0][2] + xr2 * a[0][3]);
      gc[13] = zr * (a[0][2] + xr2 * a[0][3]);
      gc[14] = xr * (a[0][2] + yr2 * a[0][3]);
      gc[15] = xr * yr * zr * a[0][3];
      gc[16] = xr * (a[0][2] + zr2 * a[0][3]);
      gc[17] = yr * (3.0 * a[0][2] + yr2 * a[0][3]);
      gc[18] = zr * (a[0][2] + yr2 * a[0][3]);
      gc[19] = yr * (a[0][2] + zr2 * a[0][3]);
      gc[20] = zr * (3.0 * a[0][2] + zr2 * a[0][3]);

      gux[11] = 3.0 * a[1][1] + xr2 * (6.0 * a[1][2] + xr2 * a[1][3]);
      gux[12] = xr * yr * (3.0 * a[1][2] + xr2 * a[1][3]);
      gux[13] = xr * zr * (3.0 * a[1][2] + xr2 * a[1][3]);
      gux[14] = a[1][1] + (xr2 + yr2) * a[1][2] + xr2 * yr2 * a[1][3];
      gux[15] = yr * zr * (a[1][2] + xr2 * a[1][3]);
      gux[16] = a[1][1] + (xr2 + zr2) * a[1][2] + xr2 * zr2 * a[1][3];
      gux[17] = xr * yr * (3.0 * a[1][2] + yr2 * a[1][3]);
      gux[18] = xr * zr * (a[1][2] + yr2 * a[1][3]);
      gux[19] = xr * yr * (a[1][2] + zr2 * a[1][3]);
      gux[20] = xr * zr * (3.0 * a[1][2] + zr2 * a[1][3]);

      guy[11] = gux[12];
      guy[12] = gux[14];
      guy[13] = gux[15];
      guy[14] = gux[17];
      guy[15] = gux[18];
      guy[16] = gux[19];
      guy[17] = 3.0 * a[1][1] + yr2 * (6.0 * a[1][2] + yr2 * a[1][3]);
      guy[18] = yr * zr * (3.0 * a[1][2] + yr2 * a[1][3]);
      guy[19] = a[1][1] + (yr2 + zr2) * a[1][2] + yr2 * zr2 * a[1][3];
      guy[20] = yr * zr * (3.0 * a[1][2] + zr2 * a[1][3]);

      guz[11] = gux[13];
      guz[12] = gux[15];
      guz[13] = gux[16];
      guz[14] = gux[18];
      guz[15] = gux[19];
      guz[16] = gux[20];
      guz[17] = guy[18];
      guz[18] = guy[19];
      guz[19] = guy[20];
      guz[20] = 3.0 * a[1][1] + zr2 * (6.0 * a[1][2] + zr2 * a[1][3]);

      gqxx[11] = xr * (12.0 * a[2][1] + xr2 * (9.0 * a[2][2] + xr2 * a[2][3]));
      gqxx[12] = yr * (2.0 * a[2][1] + xr2 * (5.0 * a[2][2] + xr2 * a[2][3]));
      gqxx[13] = zr * (2.0 * a[2][1] + xr2 * (5.0 * a[2][2] + xr2 * a[2][3]));
      gqxx[14] = xr * (2.0 * a[2][1] + yr2 * 2.0 * a[2][2] + xr2 * (a[2][2] + yr2 * a[2][3]));
      gqxx[15] = xr * yr * zr * (2.0 * a[2][2] + xr2 * a[2][3]);
      gqxx[16] = xr * (2.0 * a[2][1] + zr2 * 2.0 * a[2][2] + xr2 * (a[2][2] + zr2 * a[2][3]));
      gqxx[17] = yr * xr2 * (3.0 * a[2][2] + yr2 * a[2][3]);
      gqxx[18] = zr * xr2 * (a[2][2] + yr2 * a[2][3]);
      gqxx[19] = yr * xr2 * (a[2][2] + zr2 * a[2][3]);
      gqxx[20] = zr * xr2 * (3.0 * a[2][2] + zr2 * a[2][3]);

      gqxy[11] = yr * (3.0 * a[2][1] + xr2 * (6.0 * a[2][2] + xr2 * a[2][3]));
      gqxy[12] = xr * (3.0 * (a[2][1] + yr2 * a[2][2]) + xr2 * (a[2][2] + yr2 * a[2][3]));
      gqxy[13] = xr * yr * zr * (3.0 * a[2][2] + xr2 * a[2][3]);
      gqxy[14] = yr * (3.0 * (a[2][1] + xr2 * a[2][2]) + yr2 * (a[2][2] + xr2 * a[2][3]));
      gqxy[15] = zr * (a[2][1] + (yr2 + xr2) * a[2][2] + yr2 * xr2 * a[2][3]);
      gqxy[16] = yr * (a[2][1] + (xr2 + zr2) * a[2][2] + xr2 * zr2 * a[2][3]);
      gqxy[17] = xr * (3.0 * (a[2][1] + yr2 * a[2][2]) + yr2 * (3.0 * a[2][2] + yr2 * a[2][3]));
      gqxy[18] = xr * yr * zr * (3.0 * a[2][2] + yr2 * a[2][3]);
      gqxy[19] = xr * (a[2][1] + (yr2 + zr2) * a[2][2] + yr2 * zr2 * a[2][3]);
      gqxy[20] = xr * yr * zr * (3.0 * a[2][2] + zr2 * a[2][3]);

      gqxz[11] = zr * (3.0 * a[2][1] + xr2 * (6.0 * a[2][2] + xr2 * a[2][3]));
      gqxz[12] = xr * yr * zr * (3.0 * a[2][2] + xr2 * a[2][3]);
      gqxz[13] = xr * (3.0 * (a[2][1] + zr2 * a[2][2]) + xr2 * (a[2][2] + zr2 * a[2][3]));
      gqxz[14] = zr * (a[2][1] + (xr2 + yr2) * a[2][2] + xr2 * yr2 * a[2][3]);
      gqxz[15] = yr * (a[2][1] + (xr2 + zr2) * a[2][2] + zr2 * xr2 * a[2][3]);
      gqxz[16] = zr * (3.0 * (a[2][1] + xr2 * a[2][2]) + zr2 * (a[2][2] + xr2 * a[2][3]));
      gqxz[17] = xr * yr * zr * (3.0 * a[2][2] + yr2 * a[2][3]);
      gqxz[18] = xr * (a[2][1] + (zr2 + yr2) * a[2][2] + zr2 * yr2 * a[2][3]);
      gqxz[19] = xr * yr * zr * (3.0 * a[2][2] + zr2 * a[2][3]);
      gqxz[20] = xr * (3.0 * a[2][1] + zr2 * (6.0 * a[2][2] + zr2 * a[2][3]));

      gqyy[11] = xr * yr2 * (3.0 * a[2][2] + xr2 * a[2][3]);
      gqyy[12] = yr * (2.0 * a[2][1] + xr2 * 2.0 * a[2][2] + yr2 * (a[2][2] + xr2 * a[2][3]));
      gqyy[13] = zr * yr2 * (a[2][2] + xr2 * a[2][3]);
      gqyy[14] = xr * (2.0 * a[2][1] + yr2 * (5.0 * a[2][2] + yr2 * a[2][3]));
      gqyy[15] = xr * yr * zr * (2.0 * a[2][2] + yr2 * a[2][3]);
      gqyy[16] = xr * yr2 * (a[2][2] + zr2 * a[2][3]);
      gqyy[17] = yr * (12.0 * a[2][1] + yr2 * (9.0 * a[2][2] + yr2 * a[2][3]));
      gqyy[18] = zr * (2.0 * a[2][1] + yr2 * (5.0 * a[2][2] + yr2 * a[2][3]));
      gqyy[19] = yr * (2.0 * a[2][1] + zr2 * 2.0 * a[2][2] + yr2 * (a[2][2] + zr2 * a[2][3]));
      gqyy[20] = zr * yr2 * (3.0 * a[2][2] + zr2 * a[2][3]);

      gqyz[11] = xr * yr * zr * (3.0 * a[2][2] + xr2 * a[2][3]);
      gqyz[12] = zr * (a[2][1] + (xr2 + yr2) * a[2][2] + xr2 * yr2 * a[2][3]);
      gqyz[13] = yr * (a[2][1] + (xr2 + zr2) * a[2][2] + xr2 * zr2 * a[2][3]);
      gqyz[14] = xr * yr * zr * (3.0 * a[2][2] + yr2 * a[2][3]);
      gqyz[15] = xr * (a[2][1] + (yr2 + zr2) * a[2][2] + yr2 * zr2 * a[2][3]);
      gqyz[16] = xr * yr * zr * (3.0 * a[2][2] + zr2 * a[2][3]);
      gqyz[17] = zr * (3.0 * a[2][1] + yr2 * (6.0 * a[2][2] + yr2 * a[2][3]));
      gqyz[18] = yr * (3.0 * (a[2][1] + zr2 * a[2][2]) + yr2 * (a[2][2] + zr2 * a[2][3]));
      gqyz[19] = zr * (3.0 * (a[2][1] + yr2 * a[2][2]) + zr2 * (a[2][2] + yr2 * a[2][3]));
      gqyz[20] = yr * (3.0 * a[2][1] + zr2 * (6.0 * a[2][2] + zr2 * a[2][3]));

      gqzz[11] = xr * zr2 * (3.0 * a[2][2] + xr2 * a[2][3]);
      gqzz[12] = yr * (zr2 * a[2][2] + xr2 * (zr2 * a[2][3]));
      gqzz[13] = zr * (2.0 * a[2][1] + xr2 * 2.0 * a[2][2] + zr2 * (a[2][2] + xr2 * a[2][3]));
      gqzz[14] = xr * zr2 * (a[2][2] + yr2 * a[2][3]);
      gqzz[15] = xr * yr * zr * (2.0 * a[2][2] + zr2 * a[2][3]);
      gqzz[16] = xr * (2.0 * a[2][1] + zr2 * (5.0 * a[2][2] + zr2 * a[2][3]));
      gqzz[17] = yr * zr2 * (3.0 * a[2][2] + yr2 * a[2][3]);
      gqzz[18] = zr * (2.0 * a[2][1] + yr2 * 2.0 * a[2][2] + zr2 * (a[2][2] + yr2 * a[2][3]));
      gqzz[19] = yr * (2.0 * a[2][1] + zr2 * (5.0 * a[2][2] + zr2 * a[2][3]));
      gqzz[20] = zr * (12.0 * a[2][1] + zr2 * (9.0 * a[2][2] + zr2 * a[2][3]));
    }

    private void permanentEnergyGradient(int i, int k) {
      final double desymdr =
          ci * ck * gc[21]
              - (uxi * (uxk * gux[22] + uyk * guy[22] + uzk * guz[22])
              + uyi * (uxk * gux[23] + uyk * guy[23] + uzk * guz[23])
              + uzi * (uxk * gux[24] + uyk * guy[24] + uzk * guz[24]));
      final double dewidr =
          ci * (uxk * gc[22] + uyk * gc[23] + uzk * gc[24])
              - ck * (uxi * gux[21] + uyi * guy[21] + uzi * guz[21])
              + ci
              * (qxxk * gc[25]
              + qyyk * gc[28]
              + qzzk * gc[30]
              + 2.0 * (qxyk * gc[26] + qxzk * gc[27] + qyzk * gc[29]))
              + ck
              * (qxxi * gqxx[21]
              + qyyi * gqyy[21]
              + qzzi * gqzz[21]
              + 2.0 * (qxyi * gqxy[21] + qxzi * gqxz[21] + qyzi * gqyz[21]))
              - uxi
              * (qxxk * gux[25]
              + qyyk * gux[28]
              + qzzk * gux[30]
              + 2.0 * (qxyk * gux[26] + qxzk * gux[27] + qyzk * gux[29]))
              - uyi
              * (qxxk * guy[25]
              + qyyk * guy[28]
              + qzzk * guy[30]
              + 2.0 * (qxyk * guy[26] + qxzk * guy[27] + qyzk * guy[29]))
              - uzi
              * (qxxk * guz[25]
              + qyyk * guz[28]
              + qzzk * guz[30]
              + 2.0 * (qxyk * guz[26] + qxzk * guz[27] + qyzk * guz[29]))
              + uxk
              * (qxxi * gqxx[22]
              + qyyi * gqyy[22]
              + qzzi * gqzz[22]
              + 2.0 * (qxyi * gqxy[22] + qxzi * gqxz[22] + qyzi * gqyz[22]))
              + uyk
              * (qxxi * gqxx[23]
              + qyyi * gqyy[23]
              + qzzi * gqzz[23]
              + 2.0 * (qxyi * gqxy[23] + qxzi * gqxz[23] + qyzi * gqyz[23]))
              + uzk
              * (qxxi * gqxx[24]
              + qyyi * gqyy[24]
              + qzzi * gqzz[24]
              + 2.0 * (qxyi * gqxy[24] + qxzi * gqxz[24] + qyzi * gqyz[24]))
              + qxxi
              * (qxxk * gqxx[25]
              + qyyk * gqxx[28]
              + qzzk * gqxx[30]
              + 2.0 * (qxyk * gqxx[26] + qxzk * gqxx[27] + qyzk * gqxx[29]))
              + qyyi
              * (qxxk * gqyy[25]
              + qyyk * gqyy[28]
              + qzzk * gqyy[30]
              + 2.0 * (qxyk * gqyy[26] + qxzk * gqyy[27] + qyzk * gqyy[29]))
              + qzzi
              * (qxxk * gqzz[25]
              + qyyk * gqzz[28]
              + qzzk * gqzz[30]
              + 2.0 * (qxyk * gqzz[26] + qxzk * gqzz[27] + qyzk * gqzz[29]))
              + 2.0
              * (qxyi
              * (qxxk * gqxy[25]
              + qyyk * gqxy[28]
              + qzzk * gqxy[30]
              + 2.0 * (qxyk * gqxy[26] + qxzk * gqxy[27] + qyzk * gqxy[29]))
              + qxzi
              * (qxxk * gqxz[25]
              + qyyk * gqxz[28]
              + qzzk * gqxz[30]
              + 2.0 * (qxyk * gqxz[26] + qxzk * gqxz[27] + qyzk * gqxz[29]))
              + qyzi
              * (qxxk * gqyz[25]
              + qyyk * gqyz[28]
              + qzzk * gqyz[30]
              + 2.0 * (qxyk * gqyz[26] + qxzk * gqyz[27] + qyzk * gqyz[29])));
      final double dewkdr =
          ci * (uxk * gux[21] + uyk * guy[21] + uzk * guz[21])
              - ck * (uxi * gc[22] + uyi * gc[23] + uzi * gc[24])
              + ci
              * (qxxk * gqxx[21]
              + qyyk * gqyy[21]
              + qzzk * gqzz[21]
              + 2.0 * (qxyk * gqxy[21] + qxzk * gqxz[21] + qyzk * gqyz[21]))
              + ck
              * (qxxi * gc[25]
              + qyyi * gc[28]
              + qzzi * gc[30]
              + 2.0 * (qxyi * gc[26] + qxzi * gc[27] + qyzi * gc[29]))
              - uxi
              * (qxxk * gqxx[22]
              + qyyk * gqyy[22]
              + qzzk * gqzz[22]
              + 2.0 * (qxyk * gqxy[22] + qxzk * gqxz[22] + qyzk * gqyz[22]))
              - uyi
              * (qxxk * gqxx[23]
              + qyyk * gqyy[23]
              + qzzk * gqzz[23]
              + 2.0 * (qxyk * gqxy[23] + qxzk * gqxz[23] + qyzk * gqyz[23]))
              - uzi
              * (qxxk * gqxx[24]
              + qyyk * gqyy[24]
              + qzzk * gqzz[24]
              + 2.0 * (qxyk * gqxy[24] + qxzk * gqxz[24] + qyzk * gqyz[24]))
              + uxk
              * (qxxi * gux[25]
              + qyyi * gux[28]
              + qzzi * gux[30]
              + 2.0 * (qxyi * gux[26] + qxzi * gux[27] + qyzi * gux[29]))
              + uyk
              * (qxxi * guy[25]
              + qyyi * guy[28]
              + qzzi * guy[30]
              + 2.0 * (qxyi * guy[26] + qxzi * guy[27] + qyzi * guy[29]))
              + uzk
              * (qxxi * guz[25]
              + qyyi * guz[28]
              + qzzi * guz[30]
              + 2.0 * (qxyi * guz[26] + qxzi * guz[27] + qyzi * guz[29]))
              + qxxi
              * (qxxk * gqxx[25]
              + qyyk * gqyy[25]
              + qzzk * gqzz[25]
              + 2.0 * (qxyk * gqxy[25] + qxzk * gqxz[25] + qyzk * gqyz[25]))
              + qyyi
              * (qxxk * gqxx[28]
              + qyyk * gqyy[28]
              + qzzk * gqzz[28]
              + 2.0 * (qxyk * gqxy[28] + qxzk * gqxz[28] + qyzk * gqyz[28]))
              + qzzi
              * (qxxk * gqxx[30]
              + qyyk * gqyy[30]
              + qzzk * gqzz[30]
              + 2.0 * (qxyk * gqxy[30] + qxzk * gqxz[30] + qyzk * gqyz[30]))
              + 2.0
              * (qxyi
              * (qxxk * gqxx[26]
              + qyyk * gqyy[26]
              + qzzk * gqzz[26]
              + 2.0 * (qxyk * gqxy[26] + qxzk * gqxz[26] + qyzk * gqyz[26]))
              + qxzi
              * (qxxk * gqxx[27]
              + qyyk * gqyy[27]
              + qzzk * gqzz[27]
              + 2.0 * (qxyk * gqxy[27] + qxzk * gqxz[27] + qyzk * gqyz[27]))
              + qyzi
              * (qxxk * gqxx[29]
              + qyyk * gqyy[29]
              + qzzk * gqzz[29]
              + 2.0 * (qxyk * gqxy[29] + qxzk * gqxz[29] + qyzk * gqyz[29])));
      final double dsumdr = desymdr + 0.5 * (dewidr + dewkdr);
      final double drbi = rbk * dsumdr;
      final double drbk = rbi * dsumdr;

      double selfScale = 1.0;
      if (i == k) {
        if (iSymm == 0) {
          sharedBornGrad.add(threadID, i, drbi);
          return;
        } else {
          selfScale = 0.5;
        }
      }

      // Increment the gradients and Born chain rule term.
      final double dedx = selfScale * dEdX();
      final double dedy = selfScale * dEdY();
      final double dedz = selfScale * dEdZ();
      dedxi -= dedx;
      dedyi -= dedy;
      dedzi -= dedz;
      dborni += selfScale * drbi;

      final double dedxk = dedx * transOp[0][0] + dedy * transOp[1][0] + dedz * transOp[2][0];
      final double dedyk = dedx * transOp[0][1] + dedy * transOp[1][1] + dedz * transOp[2][1];
      final double dedzk = dedx * transOp[0][2] + dedy * transOp[1][2] + dedz * transOp[2][2];

      grad.add(threadID, k, dedxk, dedyk, dedzk);
      sharedBornGrad.add(threadID, k, selfScale * drbk);
      permanentEnergyTorque(i, k);
    }

    private double dEdZ() {
      final double desymdz =
          ci * ck * gc[4]
              - (uxi * (uxk * gux[7] + uyk * guy[7] + uzk * guz[7])
              + uyi * (uxk * gux[9] + uyk * guy[9] + uzk * guz[9])
              + uzi * (uxk * gux[10] + uyk * guy[10] + uzk * guz[10]));
      final double dewidz =
          ci * (uxk * gc[7] + uyk * gc[9] + uzk * gc[10])
              - ck * (uxi * gux[4] + uyi * guy[4] + uzi * guz[4])
              + ci
              * (qxxk * gc[13]
              + qyyk * gc[18]
              + qzzk * gc[20]
              + 2.0 * (qxyk * gc[15] + qxzk * gc[16] + qyzk * gc[19]))
              + ck
              * (qxxi * gqxx[4]
              + qyyi * gqyy[4]
              + qzzi * gqzz[4]
              + 2.0 * (qxyi * gqxy[4] + qxzi * gqxz[4] + qyzi * gqyz[4]))
              - uxi
              * (qxxk * gux[13]
              + qyyk * gux[18]
              + qzzk * gux[20]
              + 2.0 * (qxyk * gux[15] + qxzk * gux[16] + qyzk * gux[19]))
              - uyi
              * (qxxk * guy[13]
              + qyyk * guy[18]
              + qzzk * guy[20]
              + 2.0 * (qxyk * guy[15] + qxzk * guy[16] + qyzk * guy[19]))
              - uzi
              * (qxxk * guz[13]
              + qyyk * guz[18]
              + qzzk * guz[20]
              + 2.0 * (qxyk * guz[15] + qxzk * guz[16] + qyzk * guz[19]))
              + uxk
              * (qxxi * gqxx[7]
              + qyyi * gqyy[7]
              + qzzi * gqzz[7]
              + 2.0 * (qxyi * gqxy[7] + qxzi * gqxz[7] + qyzi * gqyz[7]))
              + uyk
              * (qxxi * gqxx[9]
              + qyyi * gqyy[9]
              + qzzi * gqzz[9]
              + 2.0 * (qxyi * gqxy[9] + qxzi * gqxz[9] + qyzi * gqyz[9]))
              + uzk
              * (qxxi * gqxx[10]
              + qyyi * gqyy[10]
              + qzzi * gqzz[10]
              + 2.0 * (qxyi * gqxy[10] + qxzi * gqxz[10] + qyzi * gqyz[10]))
              + qxxi
              * (qxxk * gqxx[13]
              + qyyk * gqxx[18]
              + qzzk * gqxx[20]
              + 2.0 * (qxyk * gqxx[15] + qxzk * gqxx[16] + qyzk * gqxx[19]))
              + qyyi
              * (qxxk * gqyy[13]
              + qyyk * gqyy[18]
              + qzzk * gqyy[20]
              + 2.0 * (qxyk * gqyy[15] + qxzk * gqyy[16] + qyzk * gqyy[19]))
              + qzzi
              * (qxxk * gqzz[13]
              + qyyk * gqzz[18]
              + qzzk * gqzz[20]
              + 2.0 * (qxyk * gqzz[15] + qxzk * gqzz[16] + qyzk * gqzz[19]))
              + 2.0
              * (qxyi
              * (qxxk * gqxy[13]
              + qyyk * gqxy[18]
              + qzzk * gqxy[20]
              + 2.0 * (qxyk * gqxy[15] + qxzk * gqxy[16] + qyzk * gqxy[19]))
              + qxzi
              * (qxxk * gqxz[13]
              + qyyk * gqxz[18]
              + qzzk * gqxz[20]
              + 2.0 * (qxyk * gqxz[15] + qxzk * gqxz[16] + qyzk * gqxz[19]))
              + qyzi
              * (qxxk * gqyz[13]
              + qyyk * gqyz[18]
              + qzzk * gqyz[20]
              + 2.0 * (qxyk * gqyz[15] + qxzk * gqyz[16] + qyzk * gqyz[19])));
      final double dewkdz =
          ci * (uxk * gux[4] + uyk * guy[4] + uzk * guz[4])
              - ck * (uxi * gc[7] + uyi * gc[9] + uzi * gc[10])
              + ci
              * (qxxk * gqxx[4]
              + qyyk * gqyy[4]
              + qzzk * gqzz[4]
              + 2.0 * (qxyk * gqxy[4] + qxzk * gqxz[4] + qyzk * gqyz[4]))
              + ck
              * (qxxi * gc[13]
              + qyyi * gc[18]
              + qzzi * gc[20]
              + 2.0 * (qxyi * gc[15] + qxzi * gc[16] + qyzi * gc[19]))
              - uxi
              * (qxxk * gqxx[7]
              + qyyk * gqyy[7]
              + qzzk * gqzz[7]
              + 2.0 * (qxyk * gqxy[7] + qxzk * gqxz[7] + qyzk * gqyz[7]))
              - uyi
              * (qxxk * gqxx[9]
              + qyyk * gqyy[9]
              + qzzk * gqzz[9]
              + 2.0 * (qxyk * gqxy[9] + qxzk * gqxz[9] + qyzk * gqyz[9]))
              - uzi
              * (qxxk * gqxx[10]
              + qyyk * gqyy[10]
              + qzzk * gqzz[10]
              + 2.0 * (qxyk * gqxy[10] + qxzk * gqxz[10] + qyzk * gqyz[10]))
              + uxk
              * (qxxi * gux[13]
              + qyyi * gux[18]
              + qzzi * gux[20]
              + 2.0 * (qxyi * gux[15] + qxzi * gux[16] + qyzi * gux[19]))
              + uyk
              * (qxxi * guy[13]
              + qyyi * guy[18]
              + qzzi * guy[20]
              + 2.0 * (qxyi * guy[15] + qxzi * guy[16] + qyzi * guy[19]))
              + uzk
              * (qxxi * guz[13]
              + qyyi * guz[18]
              + qzzi * guz[20]
              + 2.0 * (qxyi * guz[15] + qxzi * guz[16] + qyzi * guz[19]))
              + qxxi
              * (qxxk * gqxx[13]
              + qyyk * gqyy[13]
              + qzzk * gqzz[13]
              + 2.0 * (qxyk * gqxy[13] + qxzk * gqxz[13] + qyzk * gqyz[13]))
              + qyyi
              * (qxxk * gqxx[18]
              + qyyk * gqyy[18]
              + qzzk * gqzz[18]
              + 2.0 * (qxyk * gqxy[18] + qxzk * gqxz[18] + qyzk * gqyz[18]))
              + qzzi
              * (qxxk * gqxx[20]
              + qyyk * gqyy[20]
              + qzzk * gqzz[20]
              + 2.0 * (qxyk * gqxy[20] + qxzk * gqxz[20] + qyzk * gqyz[20]))
              + 2.0
              * (qxyi
              * (qxxk * gqxx[15]
              + qyyk * gqyy[15]
              + qzzk * gqzz[15]
              + 2.0 * (qxyk * gqxy[15] + qxzk * gqxz[15] + qyzk * gqyz[15]))
              + qxzi
              * (qxxk * gqxx[16]
              + qyyk * gqyy[16]
              + qzzk * gqzz[16]
              + 2.0 * (qxyk * gqxy[16] + qxzk * gqxz[16] + qyzk * gqyz[16]))
              + qyzi
              * (qxxk * gqxx[19]
              + qyyk * gqyy[19]
              + qzzk * gqzz[19]
              + 2.0 * (qxyk * gqxy[19] + qxzk * gqxz[19] + qyzk * gqyz[19])));
      return desymdz + 0.5 * (dewidz + dewkdz);
    }

    private double dEdY() {
      final double desymdy =
          ci * ck * gc[3]
              - (uxi * (uxk * gux[6] + uyk * guy[6] + uzk * guz[6])
              + uyi * (uxk * gux[8] + uyk * guy[8] + uzk * guz[8])
              + uzi * (uxk * gux[9] + uyk * guy[9] + uzk * guz[9]));
      final double dewidy =
          ci * (uxk * gc[6] + uyk * gc[8] + uzk * gc[9])
              - ck * (uxi * gux[3] + uyi * guy[3] + uzi * guz[3])
              + ci
              * (qxxk * gc[12]
              + qyyk * gc[17]
              + qzzk * gc[19]
              + 2.0 * (qxyk * gc[14] + qxzk * gc[15] + qyzk * gc[18]))
              + ck
              * (qxxi * gqxx[3]
              + qyyi * gqyy[3]
              + qzzi * gqzz[3]
              + 2.0 * (qxyi * gqxy[3] + qxzi * gqxz[3] + qyzi * gqyz[3]))
              - uxi
              * (qxxk * gux[12]
              + qyyk * gux[17]
              + qzzk * gux[19]
              + 2.0 * (qxyk * gux[14] + qxzk * gux[15] + qyzk * gux[18]))
              - uyi
              * (qxxk * guy[12]
              + qyyk * guy[17]
              + qzzk * guy[19]
              + 2.0 * (qxyk * guy[14] + qxzk * guy[15] + qyzk * guy[18]))
              - uzi
              * (qxxk * guz[12]
              + qyyk * guz[17]
              + qzzk * guz[19]
              + 2.0 * (qxyk * guz[14] + qxzk * guz[15] + qyzk * guz[18]))
              + uxk
              * (qxxi * gqxx[6]
              + qyyi * gqyy[6]
              + qzzi * gqzz[6]
              + 2.0 * (qxyi * gqxy[6] + qxzi * gqxz[6] + qyzi * gqyz[6]))
              + uyk
              * (qxxi * gqxx[8]
              + qyyi * gqyy[8]
              + qzzi * gqzz[8]
              + 2.0 * (qxyi * gqxy[8] + qxzi * gqxz[8] + qyzi * gqyz[8]))
              + uzk
              * (qxxi * gqxx[9]
              + qyyi * gqyy[9]
              + qzzi * gqzz[9]
              + 2.0 * (qxyi * gqxy[9] + qxzi * gqxz[9] + qyzi * gqyz[9]))
              + qxxi
              * (qxxk * gqxx[12]
              + qyyk * gqxx[17]
              + qzzk * gqxx[19]
              + 2.0 * (qxyk * gqxx[14] + qxzk * gqxx[15] + qyzk * gqxx[18]))
              + qyyi
              * (qxxk * gqyy[12]
              + qyyk * gqyy[17]
              + qzzk * gqyy[19]
              + 2.0 * (qxyk * gqyy[14] + qxzk * gqyy[15] + qyzk * gqyy[18]))
              + qzzi
              * (qxxk * gqzz[12]
              + qyyk * gqzz[17]
              + qzzk * gqzz[19]
              + 2.0 * (qxyk * gqzz[14] + qxzk * gqzz[15] + qyzk * gqzz[18]))
              + 2.0
              * (qxyi
              * (qxxk * gqxy[12]
              + qyyk * gqxy[17]
              + qzzk * gqxy[19]
              + 2.0 * (qxyk * gqxy[14] + qxzk * gqxy[15] + qyzk * gqxy[18]))
              + qxzi
              * (qxxk * gqxz[12]
              + qyyk * gqxz[17]
              + qzzk * gqxz[19]
              + 2.0 * (qxyk * gqxz[14] + qxzk * gqxz[15] + qyzk * gqxz[18]))
              + qyzi
              * (qxxk * gqyz[12]
              + qyyk * gqyz[17]
              + qzzk * gqyz[19]
              + 2.0 * (qxyk * gqyz[14] + qxzk * gqyz[15] + qyzk * gqyz[18])));
      final double dewkdy =
          ci * (uxk * gux[3] + uyk * guy[3] + uzk * guz[3])
              - ck * (uxi * gc[6] + uyi * gc[8] + uzi * gc[9])
              + ci
              * (qxxk * gqxx[3]
              + qyyk * gqyy[3]
              + qzzk * gqzz[3]
              + 2.0 * (qxyk * gqxy[3] + qxzk * gqxz[3] + qyzk * gqyz[3]))
              + ck
              * (qxxi * gc[12]
              + qyyi * gc[17]
              + qzzi * gc[19]
              + 2.0 * (qxyi * gc[14] + qxzi * gc[15] + qyzi * gc[18]))
              - uxi
              * (qxxk * gqxx[6]
              + qyyk * gqyy[6]
              + qzzk * gqzz[6]
              + 2.0 * (qxyk * gqxy[6] + qxzk * gqxz[6] + qyzk * gqyz[6]))
              - uyi
              * (qxxk * gqxx[8]
              + qyyk * gqyy[8]
              + qzzk * gqzz[8]
              + 2.0 * (qxyk * gqxy[8] + qxzk * gqxz[8] + qyzk * gqyz[8]))
              - uzi
              * (qxxk * gqxx[9]
              + qyyk * gqyy[9]
              + qzzk * gqzz[9]
              + 2.0 * (qxyk * gqxy[9] + qxzk * gqxz[9] + qyzk * gqyz[9]))
              + uxk
              * (qxxi * gux[12]
              + qyyi * gux[17]
              + qzzi * gux[19]
              + 2.0 * (qxyi * gux[14] + qxzi * gux[15] + qyzi * gux[18]))
              + uyk
              * (qxxi * guy[12]
              + qyyi * guy[17]
              + qzzi * guy[19]
              + 2.0 * (qxyi * guy[14] + qxzi * guy[15] + qyzi * guy[18]))
              + uzk
              * (qxxi * guz[12]
              + qyyi * guz[17]
              + qzzi * guz[19]
              + 2.0 * (qxyi * guz[14] + qxzi * guz[15] + qyzi * guz[18]))
              + qxxi
              * (qxxk * gqxx[12]
              + qyyk * gqyy[12]
              + qzzk * gqzz[12]
              + 2.0 * (qxyk * gqxy[12] + qxzk * gqxz[12] + qyzk * gqyz[12]))
              + qyyi
              * (qxxk * gqxx[17]
              + qyyk * gqyy[17]
              + qzzk * gqzz[17]
              + 2.0 * (qxyk * gqxy[17] + qxzk * gqxz[17] + qyzk * gqyz[17]))
              + qzzi
              * (qxxk * gqxx[19]
              + qyyk * gqyy[19]
              + qzzk * gqzz[19]
              + 2.0 * (qxyk * gqxy[19] + qxzk * gqxz[19] + qyzk * gqyz[19]))
              + 2.0
              * (qxyi
              * (qxxk * gqxx[14]
              + qyyk * gqyy[14]
              + qzzk * gqzz[14]
              + 2.0 * (qxyk * gqxy[14] + qxzk * gqxz[14] + qyzk * gqyz[14]))
              + qxzi
              * (qxxk * gqxx[15]
              + qyyk * gqyy[15]
              + qzzk * gqzz[15]
              + 2.0 * (qxyk * gqxy[15] + qxzk * gqxz[15] + qyzk * gqyz[15]))
              + qyzi
              * (qxxk * gqxx[18]
              + qyyk * gqyy[18]
              + qzzk * gqzz[18]
              + 2.0 * (qxyk * gqxy[18] + qxzk * gqxz[18] + qyzk * gqyz[18])));

      return desymdy + 0.5 * (dewidy + dewkdy);
    }

    private double dEdX() {
      final double desymdx =
          ci * ck * gc[2]
              - (uxi * (uxk * gux[5] + uyk * guy[5] + uzk * guz[5])
              + uyi * (uxk * gux[6] + uyk * guy[6] + uzk * guz[6])
              + uzi * (uxk * gux[7] + uyk * guy[7] + uzk * guz[7]));
      final double dewidx =
          ci * (uxk * gc[5] + uyk * gc[6] + uzk * gc[7])
              - ck * (uxi * gux[2] + uyi * guy[2] + uzi * guz[2])
              + ci
              * (qxxk * gc[11]
              + qyyk * gc[14]
              + qzzk * gc[16]
              + 2.0 * (qxyk * gc[12] + qxzk * gc[13] + qyzk * gc[15]))
              + ck
              * (qxxi * gqxx[2]
              + qyyi * gqyy[2]
              + qzzi * gqzz[2]
              + 2.0 * (qxyi * gqxy[2] + qxzi * gqxz[2] + qyzi * gqyz[2]))
              - uxi
              * (qxxk * gux[11]
              + qyyk * gux[14]
              + qzzk * gux[16]
              + 2.0 * (qxyk * gux[12] + qxzk * gux[13] + qyzk * gux[15]))
              - uyi
              * (qxxk * guy[11]
              + qyyk * guy[14]
              + qzzk * guy[16]
              + 2.0 * (qxyk * guy[12] + qxzk * guy[13] + qyzk * guy[15]))
              - uzi
              * (qxxk * guz[11]
              + qyyk * guz[14]
              + qzzk * guz[16]
              + 2.0 * (qxyk * guz[12] + qxzk * guz[13] + qyzk * guz[15]))
              + uxk
              * (qxxi * gqxx[5]
              + qyyi * gqyy[5]
              + qzzi * gqzz[5]
              + 2.0 * (qxyi * gqxy[5] + qxzi * gqxz[5] + qyzi * gqyz[5]))
              + uyk
              * (qxxi * gqxx[6]
              + qyyi * gqyy[6]
              + qzzi * gqzz[6]
              + 2.0 * (qxyi * gqxy[6] + qxzi * gqxz[6] + qyzi * gqyz[6]))
              + uzk
              * (qxxi * gqxx[7]
              + qyyi * gqyy[7]
              + qzzi * gqzz[7]
              + 2.0 * (qxyi * gqxy[7] + qxzi * gqxz[7] + qyzi * gqyz[7]))
              + qxxi
              * (qxxk * gqxx[11]
              + qyyk * gqxx[14]
              + qzzk * gqxx[16]
              + 2.0 * (qxyk * gqxx[12] + qxzk * gqxx[13] + qyzk * gqxx[15]))
              + qyyi
              * (qxxk * gqyy[11]
              + qyyk * gqyy[14]
              + qzzk * gqyy[16]
              + 2.0 * (qxyk * gqyy[12] + qxzk * gqyy[13] + qyzk * gqyy[15]))
              + qzzi
              * (qxxk * gqzz[11]
              + qyyk * gqzz[14]
              + qzzk * gqzz[16]
              + 2.0 * (qxyk * gqzz[12] + qxzk * gqzz[13] + qyzk * gqzz[15]))
              + 2.0
              * (qxyi
              * (qxxk * gqxy[11]
              + qyyk * gqxy[14]
              + qzzk * gqxy[16]
              + 2.0 * (qxyk * gqxy[12] + qxzk * gqxy[13] + qyzk * gqxy[15]))
              + qxzi
              * (qxxk * gqxz[11]
              + qyyk * gqxz[14]
              + qzzk * gqxz[16]
              + 2.0 * (qxyk * gqxz[12] + qxzk * gqxz[13] + qyzk * gqxz[15]))
              + qyzi
              * (qxxk * gqyz[11]
              + qyyk * gqyz[14]
              + qzzk * gqyz[16]
              + 2.0 * (qxyk * gqyz[12] + qxzk * gqyz[13] + qyzk * gqyz[15])));
      final double dewkdx =
          ci * (uxk * gux[2] + uyk * guy[2] + uzk * guz[2])
              - ck * (uxi * gc[5] + uyi * gc[6] + uzi * gc[7])
              + ci
              * (qxxk * gqxx[2]
              + qyyk * gqyy[2]
              + qzzk * gqzz[2]
              + 2.0 * (qxyk * gqxy[2] + qxzk * gqxz[2] + qyzk * gqyz[2]))
              + ck
              * (qxxi * gc[11]
              + qyyi * gc[14]
              + qzzi * gc[16]
              + 2.0 * (qxyi * gc[12] + qxzi * gc[13] + qyzi * gc[15]))
              - uxi
              * (qxxk * gqxx[5]
              + qyyk * gqyy[5]
              + qzzk * gqzz[5]
              + 2.0 * (qxyk * gqxy[5] + qxzk * gqxz[5] + qyzk * gqyz[5]))
              - uyi
              * (qxxk * gqxx[6]
              + qyyk * gqyy[6]
              + qzzk * gqzz[6]
              + 2.0 * (qxyk * gqxy[6] + qxzk * gqxz[6] + qyzk * gqyz[6]))
              - uzi
              * (qxxk * gqxx[7]
              + qyyk * gqyy[7]
              + qzzk * gqzz[7]
              + 2.0 * (qxyk * gqxy[7] + qxzk * gqxz[7] + qyzk * gqyz[7]))
              + uxk
              * (qxxi * gux[11]
              + qyyi * gux[14]
              + qzzi * gux[16]
              + 2.0 * (qxyi * gux[12] + qxzi * gux[13] + qyzi * gux[15]))
              + uyk
              * (qxxi * guy[11]
              + qyyi * guy[14]
              + qzzi * guy[16]
              + 2.0 * (qxyi * guy[12] + qxzi * guy[13] + qyzi * guy[15]))
              + uzk
              * (qxxi * guz[11]
              + qyyi * guz[14]
              + qzzi * guz[16]
              + 2.0 * (qxyi * guz[12] + qxzi * guz[13] + qyzi * guz[15]))
              + qxxi
              * (qxxk * gqxx[11]
              + qyyk * gqyy[11]
              + qzzk * gqzz[11]
              + 2.0 * (qxyk * gqxy[11] + qxzk * gqxz[11] + qyzk * gqyz[11]))
              + qyyi
              * (qxxk * gqxx[14]
              + qyyk * gqyy[14]
              + qzzk * gqzz[14]
              + 2.0 * (qxyk * gqxy[14] + qxzk * gqxz[14] + qyzk * gqyz[14]))
              + qzzi
              * (qxxk * gqxx[16]
              + qyyk * gqyy[16]
              + qzzk * gqzz[16]
              + 2.0 * (qxyk * gqxy[16] + qxzk * gqxz[16] + qyzk * gqyz[16]))
              + 2.0
              * (qxyi
              * (qxxk * gqxx[12]
              + qyyk * gqyy[12]
              + qzzk * gqzz[12]
              + 2.0 * (qxyk * gqxy[12] + qxzk * gqxz[12] + qyzk * gqyz[12]))
              + qxzi
              * (qxxk * gqxx[13]
              + qyyk * gqyy[13]
              + qzzk * gqzz[13]
              + 2.0 * (qxyk * gqxy[13] + qxzk * gqxz[13] + qyzk * gqyz[13]))
              + qyzi
              * (qxxk * gqxx[15]
              + qyyk * gqyy[15]
              + qzzk * gqzz[15]
              + 2.0 * (qxyk * gqxy[15] + qxzk * gqxz[15] + qyzk * gqyz[15])));

      return desymdx + 0.5 * (dewidx + dewkdx);
    }

    private void permanentEnergyTorque(int i, int k) {

      // Torque on permanent dipoles due to permanent reaction field.
      final double ix =
          uxk * gux[2]
              + uyk * gux[3]
              + uzk * gux[4]
              + 0.5
              * (ck * gux[1]
              + qxxk * gux[5]
              + qyyk * gux[8]
              + qzzk * gux[10]
              + 2.0 * (qxyk * gux[6] + qxzk * gux[7] + qyzk * gux[9])
              + ck * gc[2]
              + qxxk * gqxx[2]
              + qyyk * gqyy[2]
              + qzzk * gqzz[2]
              + 2.0 * (qxyk * gqxy[2] + qxzk * gqxz[2] + qyzk * gqyz[2]));
      final double iy =
          uxk * guy[2]
              + uyk * guy[3]
              + uzk * guy[4]
              + 0.5
              * (ck * guy[1]
              + qxxk * guy[5]
              + qyyk * guy[8]
              + qzzk * guy[10]
              + 2.0 * (qxyk * guy[6] + qxzk * guy[7] + qyzk * guy[9])
              + ck * gc[3]
              + qxxk * gqxx[3]
              + qyyk * gqyy[3]
              + qzzk * gqzz[3]
              + 2.0 * (qxyk * gqxy[3] + qxzk * gqxz[3] + qyzk * gqyz[3]));
      final double iz =
          uxk * guz[2]
              + uyk * guz[3]
              + uzk * guz[4]
              + 0.5
              * (ck * guz[1]
              + qxxk * guz[5]
              + qyyk * guz[8]
              + qzzk * guz[10]
              + 2.0 * (qxyk * guz[6] + qxzk * guz[7] + qyzk * guz[9])
              + ck * gc[4]
              + qxxk * gqxx[4]
              + qyyk * gqyy[4]
              + qzzk * gqzz[4]
              + 2.0 * (qxyk * gqxy[4] + qxzk * gqxz[4] + qyzk * gqyz[4]));
      final double kx =
          uxi * gux[2]
              + uyi * gux[3]
              + uzi * gux[4]
              - 0.5
              * (ci * gux[1]
              + qxxi * gux[5]
              + qyyi * gux[8]
              + qzzi * gux[10]
              + 2.0 * (qxyi * gux[6] + qxzi * gux[7] + qyzi * gux[9])
              + ci * gc[2]
              + qxxi * gqxx[2]
              + qyyi * gqyy[2]
              + qzzi * gqzz[2]
              + 2.0 * (qxyi * gqxy[2] + qxzi * gqxz[2] + qyzi * gqyz[2]));
      final double ky =
          uxi * guy[2]
              + uyi * guy[3]
              + uzi * guy[4]
              - 0.5
              * (ci * guy[1]
              + qxxi * guy[5]
              + qyyi * guy[8]
              + qzzi * guy[10]
              + 2.0 * (qxyi * guy[6] + qxzi * guy[7] + qyzi * guy[9])
              + ci * gc[3]
              + qxxi * gqxx[3]
              + qyyi * gqyy[3]
              + qzzi * gqzz[3]
              + 2.0 * (qxyi * gqxy[3] + qxzi * gqxz[3] + qyzi * gqyz[3]));
      final double kz =
          uxi * guz[2]
              + uyi * guz[3]
              + uzi * guz[4]
              - 0.5
              * (ci * guz[1]
              + qxxi * guz[5]
              + qyyi * guz[8]
              + qzzi * guz[10]
              + 2.0 * (qxyi * guz[6] + qxzi * guz[7] + qyzi * guz[9])
              + ci * gc[4]
              + qxxi * gqxx[4]
              + qyyi * gqyy[4]
              + qzzi * gqzz[4]
              + 2.0 * (qxyi * gqxy[4] + qxzi * gqxz[4] + qyzi * gqyz[4]));
      double tix = uyi * iz - uzi * iy;
      double tiy = uzi * ix - uxi * iz;
      double tiz = uxi * iy - uyi * ix;
      double tkx = uyk * kz - uzk * ky;
      double tky = uzk * kx - uxk * kz;
      double tkz = uxk * ky - uyk * kx;

      // Torque on quadrupoles due to permanent reaction field gradient.
      final double ixx =
          -0.5
              * (ck * gqxx[1]
              + uxk * gqxx[2]
              + uyk * gqxx[3]
              + uzk * gqxx[4]
              + qxxk * gqxx[5]
              + qyyk * gqxx[8]
              + qzzk * gqxx[10]
              + 2.0 * (qxyk * gqxx[6] + qxzk * gqxx[7] + qyzk * gqxx[9])
              + ck * gc[5]
              + uxk * gux[5]
              + uyk * guy[5]
              + uzk * guz[5]
              + qxxk * gqxx[5]
              + qyyk * gqyy[5]
              + qzzk * gqzz[5]
              + 2.0 * (qxyk * gqxy[5] + qxzk * gqxz[5] + qyzk * gqyz[5]));
      final double ixy =
          -0.5
              * (ck * gqxy[1]
              + uxk * gqxy[2]
              + uyk * gqxy[3]
              + uzk * gqxy[4]
              + qxxk * gqxy[5]
              + qyyk * gqxy[8]
              + qzzk * gqxy[10]
              + 2.0 * (qxyk * gqxy[6] + qxzk * gqxy[7] + qyzk * gqxy[9])
              + ck * gc[6]
              + uxk * gux[6]
              + uyk * guy[6]
              + uzk * guz[6]
              + qxxk * gqxx[6]
              + qyyk * gqyy[6]
              + qzzk * gqzz[6]
              + 2.0 * (qxyk * gqxy[6] + qxzk * gqxz[6] + qyzk * gqyz[6]));
      final double ixz =
          -0.5
              * (ck * gqxz[1]
              + uxk * gqxz[2]
              + uyk * gqxz[3]
              + uzk * gqxz[4]
              + qxxk * gqxz[5]
              + qyyk * gqxz[8]
              + qzzk * gqxz[10]
              + 2.0 * (qxyk * gqxz[6] + qxzk * gqxz[7] + qyzk * gqxz[9])
              + ck * gc[7]
              + uxk * gux[7]
              + uyk * guy[7]
              + uzk * guz[7]
              + qxxk * gqxx[7]
              + qyyk * gqyy[7]
              + qzzk * gqzz[7]
              + 2.0 * (qxyk * gqxy[7] + qxzk * gqxz[7] + qyzk * gqyz[7]));
      final double iyy =
          -0.5
              * (ck * gqyy[1]
              + uxk * gqyy[2]
              + uyk * gqyy[3]
              + uzk * gqyy[4]
              + qxxk * gqyy[5]
              + qyyk * gqyy[8]
              + qzzk * gqyy[10]
              + 2.0 * (qxyk * gqyy[6] + qxzk * gqyy[7] + qyzk * gqyy[9])
              + ck * gc[8]
              + uxk * gux[8]
              + uyk * guy[8]
              + uzk * guz[8]
              + qxxk * gqxx[8]
              + qyyk * gqyy[8]
              + qzzk * gqzz[8]
              + 2.0 * (qxyk * gqxy[8] + qxzk * gqxz[8] + qyzk * gqyz[8]));
      final double iyz =
          -0.5
              * (ck * gqyz[1]
              + uxk * gqyz[2]
              + uyk * gqyz[3]
              + uzk * gqyz[4]
              + qxxk * gqyz[5]
              + qyyk * gqyz[8]
              + qzzk * gqyz[10]
              + 2.0 * (qxyk * gqyz[6] + qxzk * gqyz[7] + qyzk * gqyz[9])
              + ck * gc[9]
              + uxk * gux[9]
              + uyk * guy[9]
              + uzk * guz[9]
              + qxxk * gqxx[9]
              + qyyk * gqyy[9]
              + qzzk * gqzz[9]
              + 2.0 * (qxyk * gqxy[9] + qxzk * gqxz[9] + qyzk * gqyz[9]));
      final double izz =
          -0.5
              * (ck * gqzz[1]
              + uxk * gqzz[2]
              + uyk * gqzz[3]
              + uzk * gqzz[4]
              + qxxk * gqzz[5]
              + qyyk * gqzz[8]
              + qzzk * gqzz[10]
              + 2.0 * (qxyk * gqzz[6] + qxzk * gqzz[7] + qyzk * gqzz[9])
              + ck * gc[10]
              + uxk * gux[10]
              + uyk * guy[10]
              + uzk * guz[10]
              + qxxk * gqxx[10]
              + qyyk * gqyy[10]
              + qzzk * gqzz[10]
              + 2.0 * (qxyk * gqxy[10] + qxzk * gqxz[10] + qyzk * gqyz[10]));
      final double iyx = ixy;
      final double izx = ixz;
      final double izy = iyz;
      final double kxx =
          -0.5
              * (ci * gqxx[1]
              - uxi * gqxx[2]
              - uyi * gqxx[3]
              - uzi * gqxx[4]
              + qxxi * gqxx[5]
              + qyyi * gqxx[8]
              + qzzi * gqxx[10]
              + 2.0 * (qxyi * gqxx[6] + qxzi * gqxx[7] + qyzi * gqxx[9])
              + ci * gc[5]
              - uxi * gux[5]
              - uyi * guy[5]
              - uzi * guz[5]
              + qxxi * gqxx[5]
              + qyyi * gqyy[5]
              + qzzi * gqzz[5]
              + 2.0 * (qxyi * gqxy[5] + qxzi * gqxz[5] + qyzi * gqyz[5]));
      double kxy =
          -0.5
              * (ci * gqxy[1]
              - uxi * gqxy[2]
              - uyi * gqxy[3]
              - uzi * gqxy[4]
              + qxxi * gqxy[5]
              + qyyi * gqxy[8]
              + qzzi * gqxy[10]
              + 2.0 * (qxyi * gqxy[6] + qxzi * gqxy[7] + qyzi * gqxy[9])
              + ci * gc[6]
              - uxi * gux[6]
              - uyi * guy[6]
              - uzi * guz[6]
              + qxxi * gqxx[6]
              + qyyi * gqyy[6]
              + qzzi * gqzz[6]
              + 2.0 * (qxyi * gqxy[6] + qxzi * gqxz[6] + qyzi * gqyz[6]));
      double kxz =
          -0.5
              * (ci * gqxz[1]
              - uxi * gqxz[2]
              - uyi * gqxz[3]
              - uzi * gqxz[4]
              + qxxi * gqxz[5]
              + qyyi * gqxz[8]
              + qzzi * gqxz[10]
              + 2.0 * (qxyi * gqxz[6] + qxzi * gqxz[7] + qyzi * gqxz[9])
              + ci * gc[7]
              - uxi * gux[7]
              - uyi * guy[7]
              - uzi * guz[7]
              + qxxi * gqxx[7]
              + qyyi * gqyy[7]
              + qzzi * gqzz[7]
              + 2.0 * (qxyi * gqxy[7] + qxzi * gqxz[7] + qyzi * gqyz[7]));
      double kyy =
          -0.5
              * (ci * gqyy[1]
              - uxi * gqyy[2]
              - uyi * gqyy[3]
              - uzi * gqyy[4]
              + qxxi * gqyy[5]
              + qyyi * gqyy[8]
              + qzzi * gqyy[10]
              + 2.0 * (qxyi * gqyy[6] + qxzi * gqyy[7] + qyzi * gqyy[9])
              + ci * gc[8]
              - uxi * gux[8]
              - uyi * guy[8]
              - uzi * guz[8]
              + qxxi * gqxx[8]
              + qyyi * gqyy[8]
              + qzzi * gqzz[8]
              + 2.0 * (qxyi * gqxy[8] + qxzi * gqxz[8] + qyzi * gqyz[8]));
      double kyz =
          -0.5
              * (ci * gqyz[1]
              - uxi * gqyz[2]
              - uyi * gqyz[3]
              - uzi * gqyz[4]
              + qxxi * gqyz[5]
              + qyyi * gqyz[8]
              + qzzi * gqyz[10]
              + 2.0 * (qxyi * gqyz[6] + qxzi * gqyz[7] + qyzi * gqyz[9])
              + ci * gc[9]
              - uxi * gux[9]
              - uyi * guy[9]
              - uzi * guz[9]
              + qxxi * gqxx[9]
              + qyyi * gqyy[9]
              + qzzi * gqzz[9]
              + 2.0 * (qxyi * gqxy[9] + qxzi * gqxz[9] + qyzi * gqyz[9]));
      double kzz =
          -0.5
              * (ci * gqzz[1]
              - uxi * gqzz[2]
              - uyi * gqzz[3]
              - uzi * gqzz[4]
              + qxxi * gqzz[5]
              + qyyi * gqzz[8]
              + qzzi * gqzz[10]
              + 2.0 * (qxyi * gqzz[6] + qxzi * gqzz[7] + qyzi * gqzz[9])
              + ci * gc[10]
              - uxi * gux[10]
              - uyi * guy[10]
              - uzi * guz[10]
              + qxxi * gqxx[10]
              + qyyi * gqyy[10]
              + qzzi * gqzz[10]
              + 2.0 * (qxyi * gqxy[10] + qxzi * gqxz[10] + qyzi * gqyz[10]));
      double kyx = kxy;
      double kzx = kxz;
      double kzy = kyz;
      tix += 2.0 * (qxyi * ixz + qyyi * iyz + qyzi * izz - qxzi * ixy - qyzi * iyy - qzzi * izy);
      tiy += 2.0 * (qxzi * ixx + qyzi * iyx + qzzi * izx - qxxi * ixz - qxyi * iyz - qxzi * izz);
      tiz += 2.0 * (qxxi * ixy + qxyi * iyy + qxzi * izy - qxyi * ixx - qyyi * iyx - qyzi * izx);
      tkx += 2.0 * (qxyk * kxz + qyyk * kyz + qyzk * kzz - qxzk * kxy - qyzk * kyy - qzzk * kzy);
      tky += 2.0 * (qxzk * kxx + qyzk * kyx + qzzk * kzx - qxxk * kxz - qxyk * kyz - qxzk * kzz);
      tkz += 2.0 * (qxxk * kxy + qxyk * kyy + qxzk * kzy - qxyk * kxx - qyyk * kyx - qyzk * kzx);
      if (i == k) {
        double selfScale = 0.5;
        tix *= selfScale;
        tiy *= selfScale;
        tiz *= selfScale;
        tkx *= selfScale;
        tky *= selfScale;
        tkz *= selfScale;
      }
      trqxi += tix;
      trqyi += tiy;
      trqzi += tiz;

      final double rtkx = tkx * transOp[0][0] + tky * transOp[1][0] + tkz * transOp[2][0];
      final double rtky = tkx * transOp[0][1] + tky * transOp[1][1] + tkz * transOp[2][1];
      final double rtkz = tkx * transOp[0][2] + tky * transOp[1][2] + tkz * transOp[2][2];
      torque.add(threadID, k, rtkx, rtky, rtkz);
    }

    private void polarizationEnergyGradient(int i, int k) {
      // Electrostatic solvation free energy gradient of the permanent
      // multipoles in the reaction potential of the induced dipoles.
      final double dpsymdx =
          -uxi * (sxk * gux[5] + syk * guy[5] + szk * guz[5])
              - uyi * (sxk * gux[6] + syk * guy[6] + szk * guz[6])
              - uzi * (sxk * gux[7] + syk * guy[7] + szk * guz[7])
              - uxk * (sxi * gux[5] + syi * guy[5] + szi * guz[5])
              - uyk * (sxi * gux[6] + syi * guy[6] + szi * guz[6])
              - uzk * (sxi * gux[7] + syi * guy[7] + szi * guz[7]);
      final double dpwidx =
          ci * (sxk * gc[5] + syk * gc[6] + szk * gc[7])
              - ck * (sxi * gux[2] + syi * guy[2] + szi * guz[2])
              - sxi
              * (qxxk * gux[11]
              + qyyk * gux[14]
              + qzzk * gux[16]
              + 2.0 * (qxyk * gux[12] + qxzk * gux[13] + qyzk * gux[15]))
              - syi
              * (qxxk * guy[11]
              + qyyk * guy[14]
              + qzzk * guy[16]
              + 2.0 * (qxyk * guy[12] + qxzk * guy[13] + qyzk * guy[15]))
              - szi
              * (qxxk * guz[11]
              + qyyk * guz[14]
              + qzzk * guz[16]
              + 2.0 * (qxyk * guz[12] + qxzk * guz[13] + qyzk * guz[15]))
              + sxk
              * (qxxi * gqxx[5]
              + qyyi * gqyy[5]
              + qzzi * gqzz[5]
              + 2.0 * (qxyi * gqxy[5] + qxzi * gqxz[5] + qyzi * gqyz[5]))
              + syk
              * (qxxi * gqxx[6]
              + qyyi * gqyy[6]
              + qzzi * gqzz[6]
              + 2.0 * (qxyi * gqxy[6] + qxzi * gqxz[6] + qyzi * gqyz[6]))
              + szk
              * (qxxi * gqxx[7]
              + qyyi * gqyy[7]
              + qzzi * gqzz[7]
              + 2.0 * (qxyi * gqxy[7] + qxzi * gqxz[7] + qyzi * gqyz[7]));
      final double dpwkdx =
          ci * (sxk * gux[2] + syk * guy[2] + szk * guz[2])
              - ck * (sxi * gc[5] + syi * gc[6] + szi * gc[7])
              - sxi
              * (qxxk * gqxx[5]
              + qyyk * gqyy[5]
              + qzzk * gqzz[5]
              + 2.0 * (qxyk * gqxy[5] + qxzk * gqxz[5] + qyzk * gqyz[5]))
              - syi
              * (qxxk * gqxx[6]
              + qyyk * gqyy[6]
              + qzzk * gqzz[6]
              + 2.0 * (qxyk * gqxy[6] + qxzk * gqxz[6] + qyzk * gqyz[6]))
              - szi
              * (qxxk * gqxx[7]
              + qyyk * gqyy[7]
              + qzzk * gqzz[7]
              + 2.0 * (qxyk * gqxy[7] + qxzk * gqxz[7] + qyzk * gqyz[7]))
              + sxk
              * (qxxi * gux[11]
              + qyyi * gux[14]
              + qzzi * gux[16]
              + 2.0 * (qxyi * gux[12] + qxzi * gux[13] + qyzi * gux[15]))
              + syk
              * (qxxi * guy[11]
              + qyyi * guy[14]
              + qzzi * guy[16]
              + 2.0 * (qxyi * guy[12] + qxzi * guy[13] + qyzi * guy[15]))
              + szk
              * (qxxi * guz[11]
              + qyyi * guz[14]
              + qzzi * guz[16]
              + 2.0 * (qxyi * guz[12] + qxzi * guz[13] + qyzi * guz[15]));
      final double dpsymdy =
          -uxi * (sxk * gux[6] + syk * guy[6] + szk * guz[6])
              - uyi * (sxk * gux[8] + syk * guy[8] + szk * guz[8])
              - uzi * (sxk * gux[9] + syk * guy[9] + szk * guz[9])
              - uxk * (sxi * gux[6] + syi * guy[6] + szi * guz[6])
              - uyk * (sxi * gux[8] + syi * guy[8] + szi * guz[8])
              - uzk * (sxi * gux[9] + syi * guy[9] + szi * guz[9]);
      final double dpwidy =
          ci * (sxk * gc[6] + syk * gc[8] + szk * gc[9])
              - ck * (sxi * gux[3] + syi * guy[3] + szi * guz[3])
              - sxi
              * (qxxk * gux[12]
              + qyyk * gux[17]
              + qzzk * gux[19]
              + 2.0 * (qxyk * gux[14] + qxzk * gux[15] + qyzk * gux[18]))
              - syi
              * (qxxk * guy[12]
              + qyyk * guy[17]
              + qzzk * guy[19]
              + 2.0 * (qxyk * guy[14] + qxzk * guy[15] + qyzk * guy[18]))
              - szi
              * (qxxk * guz[12]
              + qyyk * guz[17]
              + qzzk * guz[19]
              + 2.0 * (qxyk * guz[14] + qxzk * guz[15] + qyzk * guz[18]))
              + sxk
              * (qxxi * gqxx[6]
              + qyyi * gqyy[6]
              + qzzi * gqzz[6]
              + 2.0 * (qxyi * gqxy[6] + qxzi * gqxz[6] + qyzi * gqyz[6]))
              + syk
              * (qxxi * gqxx[8]
              + qyyi * gqyy[8]
              + qzzi * gqzz[8]
              + 2.0 * (qxyi * gqxy[8] + qxzi * gqxz[8] + qyzi * gqyz[8]))
              + szk
              * (qxxi * gqxx[9]
              + qyyi * gqyy[9]
              + qzzi * gqzz[9]
              + 2.0 * (qxyi * gqxy[9] + qxzi * gqxz[9] + qyzi * gqyz[9]));
      final double dpwkdy =
          ci * (sxk * gux[3] + syk * guy[3] + szk * guz[3])
              - ck * (sxi * gc[6] + syi * gc[8] + szi * gc[9])
              - sxi
              * (qxxk * gqxx[6]
              + qyyk * gqyy[6]
              + qzzk * gqzz[6]
              + 2.0 * (qxyk * gqxy[6] + qxzk * gqxz[6] + qyzk * gqyz[6]))
              - syi
              * (qxxk * gqxx[8]
              + qyyk * gqyy[8]
              + qzzk * gqzz[8]
              + 2.0 * (qxyk * gqxy[8] + qxzk * gqxz[8] + qyzk * gqyz[8]))
              - szi
              * (qxxk * gqxx[9]
              + qyyk * gqyy[9]
              + qzzk * gqzz[9]
              + 2.0 * (qxyk * gqxy[9] + qxzk * gqxz[9] + qyzk * gqyz[9]))
              + sxk
              * (qxxi * gux[12]
              + qyyi * gux[17]
              + qzzi * gux[19]
              + 2.0 * (qxyi * gux[14] + qxzi * gux[15] + qyzi * gux[18]))
              + syk
              * (qxxi * guy[12]
              + qyyi * guy[17]
              + qzzi * guy[19]
              + 2.0 * (qxyi * guy[14] + qxzi * guy[15] + qyzi * guy[18]))
              + szk
              * (qxxi * guz[12]
              + qyyi * guz[17]
              + qzzi * guz[19]
              + 2.0 * (qxyi * guz[14] + qxzi * guz[15] + qyzi * guz[18]));
      final double dpsymdz =
          -uxi * (sxk * gux[7] + syk * guy[7] + szk * guz[7])
              - uyi * (sxk * gux[9] + syk * guy[9] + szk * guz[9])
              - uzi * (sxk * gux[10] + syk * guy[10] + szk * guz[10])
              - uxk * (sxi * gux[7] + syi * guy[7] + szi * guz[7])
              - uyk * (sxi * gux[9] + syi * guy[9] + szi * guz[9])
              - uzk * (sxi * gux[10] + syi * guy[10] + szi * guz[10]);
      final double dpwidz =
          ci * (sxk * gc[7] + syk * gc[9] + szk * gc[10])
              - ck * (sxi * gux[4] + syi * guy[4] + szi * guz[4])
              - sxi
              * (qxxk * gux[13]
              + qyyk * gux[18]
              + qzzk * gux[20]
              + 2.0 * (qxyk * gux[15] + qxzk * gux[16] + qyzk * gux[19]))
              - syi
              * (qxxk * guy[13]
              + qyyk * guy[18]
              + qzzk * guy[20]
              + 2.0 * (qxyk * guy[15] + qxzk * guy[16] + qyzk * guy[19]))
              - szi
              * (qxxk * guz[13]
              + qyyk * guz[18]
              + qzzk * guz[20]
              + 2.0 * (qxyk * guz[15] + qxzk * guz[16] + qyzk * guz[19]))
              + sxk
              * (qxxi * gqxx[7]
              + qyyi * gqyy[7]
              + qzzi * gqzz[7]
              + 2.0 * (qxyi * gqxy[7] + qxzi * gqxz[7] + qyzi * gqyz[7]))
              + syk
              * (qxxi * gqxx[9]
              + qyyi * gqyy[9]
              + qzzi * gqzz[9]
              + 2.0 * (qxyi * gqxy[9] + qxzi * gqxz[9] + qyzi * gqyz[9]))
              + szk
              * (qxxi * gqxx[10]
              + qyyi * gqyy[10]
              + qzzi * gqzz[10]
              + 2.0 * (qxyi * gqxy[10] + qxzi * gqxz[10] + qyzi * gqyz[10]));
      final double dpwkdz =
          ci * (sxk * gux[4] + syk * guy[4] + szk * guz[4])
              - ck * (sxi * gc[7] + syi * gc[9] + szi * gc[10])
              - sxi
              * (qxxk * gqxx[7]
              + qyyk * gqyy[7]
              + qzzk * gqzz[7]
              + 2.0 * (qxyk * gqxy[7] + qxzk * gqxz[7] + qyzk * gqyz[7]))
              - syi
              * (qxxk * gqxx[9]
              + qyyk * gqyy[9]
              + qzzk * gqzz[9]
              + 2.0 * (qxyk * gqxy[9] + qxzk * gqxz[9] + qyzk * gqyz[9]))
              - szi
              * (qxxk * gqxx[10]
              + qyyk * gqyy[10]
              + qzzk * gqzz[10]
              + 2.0 * (qxyk * gqxy[10] + qxzk * gqxz[10] + qyzk * gqyz[10]))
              + sxk
              * (qxxi * gux[13]
              + qyyi * gux[18]
              + qzzi * gux[20]
              + 2.0 * (qxyi * gux[15] + qxzi * gux[16] + qyzi * gux[19]))
              + syk
              * (qxxi * guy[13]
              + qyyi * guy[18]
              + qzzi * guy[20]
              + 2.0 * (qxyi * guy[15] + qxzi * guy[16] + qyzi * guy[19]))
              + szk
              * (qxxi * guz[13]
              + qyyi * guz[18]
              + qzzi * guz[20]
              + 2.0 * (qxyi * guz[15] + qxzi * guz[16] + qyzi * guz[19]));

      // Effective radii chain rule terms for the electrostatic solvation free energy
      // gradient of the permanent multipoles in the reaction potential of the induced dipoles.
      final double dsymdr =
          -uxi * (sxk * gux[22] + syk * guy[22] + szk * guz[22])
              - uyi * (sxk * gux[23] + syk * guy[23] + szk * guz[23])
              - uzi * (sxk * gux[24] + syk * guy[24] + szk * guz[24])
              - uxk * (sxi * gux[22] + syi * guy[22] + szi * guz[22])
              - uyk * (sxi * gux[23] + syi * guy[23] + szi * guz[23])
              - uzk * (sxi * gux[24] + syi * guy[24] + szi * guz[24]);
      final double dwipdr =
          ci * (sxk * gc[22] + syk * gc[23] + szk * gc[24])
              - ck * (sxi * gux[21] + syi * guy[21] + szi * guz[21])
              - sxi
              * (qxxk * gux[25]
              + qyyk * gux[28]
              + qzzk * gux[30]
              + 2.0 * (qxyk * gux[26] + qxzk * gux[27] + qyzk * gux[29]))
              - syi
              * (qxxk * guy[25]
              + qyyk * guy[28]
              + qzzk * guy[30]
              + 2.0 * (qxyk * guy[26] + qxzk * guy[27] + qyzk * guy[29]))
              - szi
              * (qxxk * guz[25]
              + qyyk * guz[28]
              + qzzk * guz[30]
              + 2.0 * (qxyk * guz[26] + qxzk * guz[27] + qyzk * guz[29]))
              + sxk
              * (qxxi * gqxx[22]
              + qyyi * gqyy[22]
              + qzzi * gqzz[22]
              + 2.0 * (qxyi * gqxy[22] + qxzi * gqxz[22] + qyzi * gqyz[22]))
              + syk
              * (qxxi * gqxx[23]
              + qyyi * gqyy[23]
              + qzzi * gqzz[23]
              + 2.0 * (qxyi * gqxy[23] + qxzi * gqxz[23] + qyzi * gqyz[23]))
              + szk
              * (qxxi * gqxx[24]
              + qyyi * gqyy[24]
              + qzzi * gqzz[24]
              + 2.0 * (qxyi * gqxy[24] + qxzi * gqxz[24] + qyzi * gqyz[24]));
      final double dwkpdr =
          ci * (sxk * gux[21] + syk * guy[21] + szk * guz[21])
              - ck * (sxi * gc[22] + syi * gc[23] + szi * gc[24])
              - sxi
              * (qxxk * gqxx[22]
              + qyyk * gqyy[22]
              + qzzk * gqzz[22]
              + 2.0 * (qxyk * gqxy[22] + qxzk * gqxz[22] + qyzk * gqyz[22]))
              - syi
              * (qxxk * gqxx[23]
              + qyyk * gqyy[23]
              + qzzk * gqzz[23]
              + 2.0 * (qxyk * gqxy[23] + qxzk * gqxz[23] + qyzk * gqyz[23]))
              - szi
              * (qxxk * gqxx[24]
              + qyyk * gqyy[24]
              + qzzk * gqzz[24]
              + 2.0 * (qxyk * gqxy[24] + qxzk * gqxz[24] + qyzk * gqyz[24]))
              + sxk
              * (qxxi * gux[25]
              + qyyi * gux[28]
              + qzzi * gux[30]
              + 2.0 * (qxyi * gux[26] + qxzi * gux[27] + qyzi * gux[29]))
              + syk
              * (qxxi * guy[25]
              + qyyi * guy[28]
              + qzzi * guy[30]
              + 2.0 * (qxyi * guy[26] + qxzi * guy[27] + qyzi * guy[29]))
              + szk
              * (qxxi * guz[25]
              + qyyi * guz[28]
              + qzzi * guz[30]
              + 2.0 * (qxyi * guz[26] + qxzi * guz[27] + qyzi * guz[29]));
      double dpdx = 0.5 * (dpsymdx + 0.5 * (dpwidx + dpwkdx));
      double dpdy = 0.5 * (dpsymdy + 0.5 * (dpwidy + dpwkdy));
      double dpdz = 0.5 * (dpsymdz + 0.5 * (dpwidz + dpwkdz));
      double dsumdri = dsymdr + 0.5 * (dwipdr + dwkpdr);
      double dbi = 0.5 * rbk * dsumdri;
      double dbk = 0.5 * rbi * dsumdri;
      if (polarization == ParticleMeshEwald.Polarization.MUTUAL) {
        dpdx -=
            0.5
                * (dxi * (pxk * gux[5] + pyk * gux[6] + pzk * gux[7])
                + dyi * (pxk * guy[5] + pyk * guy[6] + pzk * guy[7])
                + dzi * (pxk * guz[5] + pyk * guz[6] + pzk * guz[7])
                + dxk * (pxi * gux[5] + pyi * gux[6] + pzi * gux[7])
                + dyk * (pxi * guy[5] + pyi * guy[6] + pzi * guy[7])
                + dzk * (pxi * guz[5] + pyi * guz[6] + pzi * guz[7]));
        dpdy -=
            0.5
                * (dxi * (pxk * gux[6] + pyk * gux[8] + pzk * gux[9])
                + dyi * (pxk * guy[6] + pyk * guy[8] + pzk * guy[9])
                + dzi * (pxk * guz[6] + pyk * guz[8] + pzk * guz[9])
                + dxk * (pxi * gux[6] + pyi * gux[8] + pzi * gux[9])
                + dyk * (pxi * guy[6] + pyi * guy[8] + pzi * guy[9])
                + dzk * (pxi * guz[6] + pyi * guz[8] + pzi * guz[9]));
        dpdz -=
            0.5
                * (dxi * (pxk * gux[7] + pyk * gux[9] + pzk * gux[10])
                + dyi * (pxk * guy[7] + pyk * guy[9] + pzk * guy[10])
                + dzi * (pxk * guz[7] + pyk * guz[9] + pzk * guz[10])
                + dxk * (pxi * gux[7] + pyi * gux[9] + pzi * gux[10])
                + dyk * (pxi * guy[7] + pyi * guy[9] + pzi * guy[10])
                + dzk * (pxi * guz[7] + pyi * guz[9] + pzi * guz[10]));
        final double duvdr =
            dxi * (pxk * gux[22] + pyk * gux[23] + pzk * gux[24])
                + dyi * (pxk * guy[22] + pyk * guy[23] + pzk * guy[24])
                + dzi * (pxk * guz[22] + pyk * guz[23] + pzk * guz[24])
                + dxk * (pxi * gux[22] + pyi * gux[23] + pzi * gux[24])
                + dyk * (pxi * guy[22] + pyi * guy[23] + pzi * guy[24])
                + dzk * (pxi * guz[22] + pyi * guz[23] + pzi * guz[24]);
        dbi -= 0.5 * rbk * duvdr;
        dbk -= 0.5 * rbi * duvdr;
      }

      // Increment the gradients and Born chain rule term.
      if (i == k && iSymm == 0) {
        dborni += dbi;
      } else {
        if (i == k) {
          dpdx *= 0.5;
          dpdy *= 0.5;
          dpdz *= 0.5;
          dbi *= 0.5;
          dbk *= 0.5;
        }
        dedxi -= dpdx;
        dedyi -= dpdy;
        dedzi -= dpdz;
        dborni += dbi;

        final double rdpdx = dpdx * transOp[0][0] + dpdy * transOp[1][0] + dpdz * transOp[2][0];
        final double rdpdy = dpdx * transOp[0][1] + dpdy * transOp[1][1] + dpdz * transOp[2][1];
        final double rdpdz = dpdx * transOp[0][2] + dpdy * transOp[1][2] + dpdz * transOp[2][2];
        grad.add(threadID, k, rdpdx, rdpdy, rdpdz);
        sharedBornGrad.add(threadID, k, dbk);
      }
      polarizationEnergyTorque(i, k);
    }

    private void polarizationEnergyTorque(int i, int k) {
      double fix = 0.5 * (sxk * gux[2] + syk * guy[2] + szk * guz[2]);
      double fiy = 0.5 * (sxk * gux[3] + syk * guy[3] + szk * guz[3]);
      double fiz = 0.5 * (sxk * gux[4] + syk * guy[4] + szk * guz[4]);
      double fkx = 0.5 * (sxi * gux[2] + syi * guy[2] + szi * guz[2]);
      double fky = 0.5 * (sxi * gux[3] + syi * guy[3] + szi * guz[3]);
      double fkz = 0.5 * (sxi * gux[4] + syi * guy[4] + szi * guz[4]);
      double fixx =
          -0.25
              * ((sxk * gqxx[2] + syk * gqxx[3] + szk * gqxx[4])
              + (sxk * gux[5] + syk * guy[5] + szk * guz[5]));
      double fixy =
          -0.25
              * ((sxk * gqxy[2] + syk * gqxy[3] + szk * gqxy[4])
              + (sxk * gux[6] + syk * guy[6] + szk * guz[6]));
      double fixz =
          -0.25
              * ((sxk * gqxz[2] + syk * gqxz[3] + szk * gqxz[4])
              + (sxk * gux[7] + syk * guy[7] + szk * guz[7]));
      double fiyy =
          -0.25
              * ((sxk * gqyy[2] + syk * gqyy[3] + szk * gqyy[4])
              + (sxk * gux[8] + syk * guy[8] + szk * guz[8]));
      double fiyz =
          -0.25
              * ((sxk * gqyz[2] + syk * gqyz[3] + szk * gqyz[4])
              + (sxk * gux[9] + syk * guy[9] + szk * guz[9]));
      double fizz =
          -0.25
              * ((sxk * gqzz[2] + syk * gqzz[3] + szk * gqzz[4])
              + (sxk * gux[10] + syk * guy[10] + szk * guz[10]));
      double fiyx = fixy;
      double fizx = fixz;
      double fizy = fiyz;
      double fkxx =
          0.25
              * ((sxi * gqxx[2] + syi * gqxx[3] + szi * gqxx[4])
              + (sxi * gux[5] + syi * guy[5] + szi * guz[5]));
      double fkxy =
          0.25
              * ((sxi * gqxy[2] + syi * gqxy[3] + szi * gqxy[4])
              + (sxi * gux[6] + syi * guy[6] + szi * guz[6]));
      double fkxz =
          0.25
              * ((sxi * gqxz[2] + syi * gqxz[3] + szi * gqxz[4])
              + (sxi * gux[7] + syi * guy[7] + szi * guz[7]));
      double fkyy =
          0.25
              * ((sxi * gqyy[2] + syi * gqyy[3] + szi * gqyy[4])
              + (sxi * gux[8] + syi * guy[8] + szi * guz[8]));
      double fkyz =
          0.25
              * ((sxi * gqyz[2] + syi * gqyz[3] + szi * gqyz[4])
              + (sxi * gux[9] + syi * guy[9] + szi * guz[9]));
      double fkzz =
          0.25
              * ((sxi * gqzz[2] + syi * gqzz[3] + szi * gqzz[4])
              + (sxi * gux[10] + syi * guy[10] + szi * guz[10]));
      double fkyx = fkxy;
      double fkzx = fkxz;
      double fkzy = fkyz;
      if (i == k) {
        fix *= 0.5;
        fiy *= 0.5;
        fiz *= 0.5;
        fkx *= 0.5;
        fky *= 0.5;
        fkz *= 0.5;
        fixx *= 0.5;
        fixy *= 0.5;
        fixz *= 0.5;
        fiyx *= 0.5;
        fiyy *= 0.5;
        fiyz *= 0.5;
        fizx *= 0.5;
        fizy *= 0.5;
        fizz *= 0.5;
        fkxx *= 0.5;
        fkxy *= 0.5;
        fkxz *= 0.5;
        fkyx *= 0.5;
        fkyy *= 0.5;
        fkyz *= 0.5;
        fkzx *= 0.5;
        fkzy *= 0.5;
        fkzz *= 0.5;
      }

      // Torque due to induced reaction field on permanent dipoles.
      double tix = uyi * fiz - uzi * fiy;
      double tiy = uzi * fix - uxi * fiz;
      double tiz = uxi * fiy - uyi * fix;
      double tkx = uyk * fkz - uzk * fky;
      double tky = uzk * fkx - uxk * fkz;
      double tkz = uxk * fky - uyk * fkx;

      // Torque due to induced reaction field gradient on quadrupoles.
      tix +=
          2.0 * (qxyi * fixz + qyyi * fiyz + qyzi * fizz - qxzi * fixy - qyzi * fiyy - qzzi * fizy);
      tiy +=
          2.0 * (qxzi * fixx + qyzi * fiyx + qzzi * fizx - qxxi * fixz - qxyi * fiyz - qxzi * fizz);
      tiz +=
          2.0 * (qxxi * fixy + qxyi * fiyy + qxzi * fizy - qxyi * fixx - qyyi * fiyx - qyzi * fizx);
      tkx +=
          2.0 * (qxyk * fkxz + qyyk * fkyz + qyzk * fkzz - qxzk * fkxy - qyzk * fkyy - qzzk * fkzy);
      tky +=
          2.0 * (qxzk * fkxx + qyzk * fkyx + qzzk * fkzx - qxxk * fkxz - qxyk * fkyz - qxzk * fkzz);
      tkz +=
          2.0 * (qxxk * fkxy + qxyk * fkyy + qxzk * fkzy - qxyk * fkxx - qyyk * fkyx - qyzk * fkzx);
      trqxi += tix;
      trqyi += tiy;
      trqzi += tiz;

      final double rx = tkx;
      final double ry = tky;
      final double rz = tkz;
      tkx = rx * transOp[0][0] + ry * transOp[1][0] + rz * transOp[2][0];
      tky = rx * transOp[0][1] + ry * transOp[1][1] + rz * transOp[2][1];
      tkz = rx * transOp[0][2] + ry * transOp[1][2] + rz * transOp[2][2];
      torque.add(threadID, k, tkx, tky, tkz);
    }
  }
}
File Project Line
ffx/potential/nonbonded/implicit/GKEnergyRegion.java Potential 484
ffx/potential/nonbonded/implicit/GKEnergyRegionQI.java Potential 426
final double[] multipolek = globalMultipole[iSymm][k];
      ck = multipolek[t000];
      uxk = multipolek[t100];
      uyk = multipolek[t010];
      uzk = multipolek[t001];
      qxxk = multipolek[t200] * oneThird;
      qxyk = multipolek[t110] * oneThird;
      qxzk = multipolek[t101] * oneThird;
      qyyk = multipolek[t020] * oneThird;
      qyzk = multipolek[t011] * oneThird;
      qzzk = multipolek[t002] * oneThird;
      dxk = inducedDipole[iSymm][k][0];
      dyk = inducedDipole[iSymm][k][1];
      dzk = inducedDipole[iSymm][k][2];
      pxk = inducedDipoleCR[iSymm][k][0];
      pyk = inducedDipoleCR[iSymm][k][1];
      pzk = inducedDipoleCR[iSymm][k][2];
      sxk = dxk + pxk;
      syk = dyk + pyk;
      szk = dzk + pzk;
      final double rb2 = rbi * rbk;
      final double expterm = exp(-r2 / (gkc * rb2));
      final double expc = expterm / gkc;
      final double expc1 = 1.0 - expc;
      final double expcr = r2 * expterm / (gkc * gkc * rb2 * rb2);
      final double dexpc = -2.0 / (gkc * rb2);
      double expcdexpc = -expc * dexpc;
      final double dexpcr = 2.0 / (gkc * rb2 * rb2);
      final double dgfdr = 0.5 * expterm * (1.0 + r2 / (rb2 * gkc));
      final double gf2 = 1.0 / (r2 + rb2 * expterm);
      final double gf = sqrt(gf2);
      final double gf3 = gf2 * gf;
      final double gf5 = gf3 * gf2;
      final double gf7 = gf5 * gf2;
      final double gf9 = gf7 * gf2;
      final double gf11 = gf9 * gf2;

      // Reaction potential auxiliary terms.
      a[0][0] = gf;
      a[1][0] = -gf3;
      a[2][0] = 3.0 * gf5;
      a[3][0] = -15.0 * gf7;
      a[4][0] = 105.0 * gf9;
      a[5][0] = -945.0 * gf11;

      // Reaction potential gradient auxiliary terms.
      a[0][1] = expc1 * a[1][0];
      a[1][1] = expc1 * a[2][0];
      a[2][1] = expc1 * a[3][0];
      a[3][1] = expc1 * a[4][0];
      a[4][1] = expc1 * a[5][0];

      // 2nd reaction potential gradient auxiliary terms.
      a[0][2] = expc1 * a[1][1] + expcdexpc * a[1][0];
      a[1][2] = expc1 * a[2][1] + expcdexpc * a[2][0];
      a[2][2] = expc1 * a[3][1] + expcdexpc * a[3][0];
      a[3][2] = expc1 * a[4][1] + expcdexpc * a[4][0];

      if (gradient) {

        // 3rd reaction potential gradient auxiliary terms.
        expcdexpc = 2.0 * expcdexpc;
        a[0][3] = expc1 * a[1][2] + expcdexpc * a[1][1];
        a[1][3] = expc1 * a[2][2] + expcdexpc * a[2][1];
        a[2][3] = expc1 * a[3][2] + expcdexpc * a[3][1];
        expcdexpc = -expc * dexpc * dexpc;
        a[0][3] = a[0][3] + expcdexpc * a[1][0];
        a[1][3] = a[1][3] + expcdexpc * a[2][0];
        a[2][3] = a[2][3] + expcdexpc * a[3][0];

        // Born radii derivatives of reaction potential auxiliary terms.
        b[0][0] = dgfdr * a[1][0];
        b[1][0] = dgfdr * a[2][0];
        b[2][0] = dgfdr * a[3][0];
        b[3][0] = dgfdr * a[4][0];
        b[4][0] = dgfdr * a[5][0];

        // Born radii gradients of reaction potential gradient auxiliary terms.
        b[0][1] = b[1][0] - expcr * a[1][0] - expc * b[1][0];
        b[1][1] = b[2][0] - expcr * a[2][0] - expc * b[2][0];
        b[2][1] = b[3][0] - expcr * a[3][0] - expc * b[3][0];
        b[3][1] = b[4][0] - expcr * a[4][0] - expc * b[4][0];

        // Born radii derivatives of the 2nd reaction potential gradient auxiliary terms.
        b[0][2] =
            b[1][1]
                - (expcr * (a[1][1] + dexpc * a[1][0])
                + expc * (b[1][1] + dexpcr * a[1][0] + dexpc * b[1][0]));
        b[1][2] =
            b[2][1]
                - (expcr * (a[2][1] + dexpc * a[2][0])
                + expc * (b[2][1] + dexpcr * a[2][0] + dexpc * b[2][0]));
        b[2][2] =
            b[3][1]
                - (expcr * (a[3][1] + dexpc * a[3][0])
                + expc * (b[3][1] + dexpcr * a[3][0] + dexpc * b[3][0]));

        // Multiply the Born radii auxiliary terms by their dielectric functions.
        b[0][0] = electric * fc * b[0][0];
        b[0][1] = electric * fc * b[0][1];
        b[0][2] = electric * fc * b[0][2];
        b[1][0] = electric * fd * b[1][0];
        b[1][1] = electric * fd * b[1][1];
        b[1][2] = electric * fd * b[1][2];
        b[2][0] = electric * fq * b[2][0];
        b[2][1] = electric * fq * b[2][1];
        b[2][2] = electric * fq * b[2][2];
      }

      // Multiply the potential auxiliary terms by their dielectric functions.
      a[0][0] = electric * fc * a[0][0];
      a[0][1] = electric * fc * a[0][1];
      a[0][2] = electric * fc * a[0][2];
      a[0][3] = electric * fc * a[0][3];
      a[1][0] = electric * fd * a[1][0];
      a[1][1] = electric * fd * a[1][1];
      a[1][2] = electric * fd * a[1][2];
      a[1][3] = electric * fd * a[1][3];
      a[2][0] = electric * fq * a[2][0];
      a[2][1] = electric * fq * a[2][1];
      a[2][2] = electric * fq * a[2][2];
      a[2][3] = electric * fq * a[2][3];

      // Compute the GK tensors required to compute the energy.
      energyTensors();

      // Compute the GK interaction energy.
      double eik = energy(i, k);

      // Test the GK Tensor class.
//     if (i == k) {
//        PolarizableMultipole polarizableMultipole = new PolarizableMultipole(multipolek,
//            inducedDipole[iSymm][k], inducedDipoleCR[iSymm][k]);
//        GeneralizedKirkwoodTensor generalizedKirkwoodTensor = new GeneralizedKirkwoodTensor(0, 2,
//            gkc, 1.0, 78.3);
//        double energy = generalizedKirkwoodTensor.selfEnergy(polarizableMultipole, rbk);
//        logger.info(format(" GK Tensor: %16.8f Code: %16.8f Rad: %16.8f", electric * energy, eik, rbk));
//      }

      gkEnergy += eik;
      count++;

      if (gradient) {
        // Compute the additional GK tensors required to compute the energy gradient.
        gradientTensors();

        // Compute the permanent GK energy gradient.
        permanentEnergyGradient(i, k);
        if (polarization != ParticleMeshEwald.Polarization.NONE) {
          // Compute the induced GK energy gradient.
          polarizationEnergyGradient(i, k);
        }
      }
    }

    private void energyTensors() {
      // Unweighted reaction potential tensor.
      gc[1] = a[0][0];
      gux[1] = xr * a[1][0];
      guy[1] = yr * a[1][0];
      guz[1] = zr * a[1][0];
      gqxx[1] = xr2 * a[2][0];
      gqyy[1] = yr2 * a[2][0];
      gqzz[1] = zr2 * a[2][0];
      gqxy[1] = xr * yr * a[2][0];
      gqxz[1] = xr * zr * a[2][0];
      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]);

      // Unweighted 2nd reaction potential gradient tensor.
      gc[5] = a[0][1] + xr2 * a[0][2];
      gc[6] = xr * yr * a[0][2];
      gc[7] = xr * zr * a[0][2];
      gc[8] = a[0][1] + yr2 * a[0][2];
      gc[9] = yr * zr * a[0][2];
      gc[10] = a[0][1] + zr2 * a[0][2];
      gux[5] = xr * (3.0 * a[1][1] + xr2 * a[1][2]);
      gux[6] = yr * (a[1][1] + xr2 * a[1][2]);
      gux[7] = zr * (a[1][1] + xr2 * a[1][2]);
      gux[8] = xr * (a[1][1] + yr2 * a[1][2]);
      gux[9] = zr * xr * yr * a[1][2];
      gux[10] = xr * (a[1][1] + zr2 * a[1][2]);
      guy[5] = yr * (a[1][1] + xr2 * a[1][2]);
      guy[6] = xr * (a[1][1] + yr2 * a[1][2]);
      guy[7] = gux[9];
      guy[8] = yr * (3.0 * a[1][1] + yr2 * a[1][2]);
      guy[9] = zr * (a[1][1] + yr2 * a[1][2]);
      guy[10] = yr * (a[1][1] + zr2 * a[1][2]);
      guz[5] = zr * (a[1][1] + xr2 * a[1][2]);
      guz[6] = gux[9];
      guz[7] = xr * (a[1][1] + zr2 * a[1][2]);
      guz[8] = zr * (a[1][1] + yr2 * a[1][2]);
      guz[9] = yr * (a[1][1] + zr2 * a[1][2]);
      guz[10] = zr * (3.0 * a[1][1] + zr2 * a[1][2]);
      gqxx[5] = 2.0 * a[2][0] + xr2 * (5.0 * a[2][1] + xr2 * a[2][2]);
      gqxx[6] = yr * xr * (2.0 * a[2][1] + xr2 * a[2][2]);
      gqxx[7] = zr * xr * (2.0 * a[2][1] + xr2 * a[2][2]);
      gqxx[8] = xr2 * (a[2][1] + yr2 * a[2][2]);
      gqxx[9] = zr * yr * xr2 * a[2][2];
      gqxx[10] = xr2 * (a[2][1] + zr2 * a[2][2]);
      gqyy[5] = yr2 * (a[2][1] + xr2 * a[2][2]);
      gqyy[6] = xr * yr * (2.0 * a[2][1] + yr2 * a[2][2]);
      gqyy[7] = xr * zr * yr2 * a[2][2];
      gqyy[8] = 2.0 * a[2][0] + yr2 * (5.0 * a[2][1] + yr2 * a[2][2]);
      gqyy[9] = yr * zr * (2.0 * a[2][1] + yr2 * a[2][2]);
      gqyy[10] = yr2 * (a[2][1] + zr2 * a[2][2]);
      gqzz[5] = zr2 * (a[2][1] + xr2 * a[2][2]);
      gqzz[6] = xr * yr * zr2 * a[2][2];
      gqzz[7] = xr * zr * (2.0 * a[2][1] + zr2 * a[2][2]);
      gqzz[8] = zr2 * (a[2][1] + yr2 * a[2][2]);
      gqzz[9] = yr * zr * (2.0 * a[2][1] + zr2 * a[2][2]);
      gqzz[10] = 2.0 * a[2][0] + zr2 * (5.0 * a[2][1] + zr2 * a[2][2]);
      gqxy[5] = xr * yr * (3.0 * a[2][1] + xr2 * a[2][2]);
      gqxy[6] = a[2][0] + (xr2 + yr2) * a[2][1] + xr2 * yr2 * a[2][2];
      gqxy[7] = zr * yr * (a[2][1] + xr2 * a[2][2]);
      gqxy[8] = xr * yr * (3.0 * a[2][1] + yr2 * a[2][2]);
      gqxy[9] = zr * xr * (a[2][1] + yr2 * a[2][2]);
      gqxy[10] = xr * yr * (a[2][1] + zr2 * a[2][2]);
      gqxz[5] = xr * zr * (3.0 * a[2][1] + xr2 * a[2][2]);
      gqxz[6] = yr * zr * (a[2][1] + xr2 * a[2][2]);
      gqxz[7] = a[2][0] + (xr2 + zr2) * a[2][1] + xr2 * zr2 * a[2][2];
      gqxz[8] = xr * zr * (a[2][1] + yr2 * a[2][2]);
      gqxz[9] = xr * yr * (a[2][1] + zr2 * a[2][2]);
      gqxz[10] = xr * zr * (3.0 * a[2][1] + zr2 * a[2][2]);
      gqyz[5] = zr * yr * (a[2][1] + xr2 * a[2][2]);
      gqyz[6] = xr * zr * (a[2][1] + yr2 * a[2][2]);
      gqyz[7] = xr * yr * (a[2][1] + zr2 * a[2][2]);
      gqyz[8] = yr * zr * (3.0 * a[2][1] + yr2 * a[2][2]);
      gqyz[9] = a[2][0] + (yr2 + zr2) * a[2][1] + yr2 * zr2 * a[2][2];
      gqyz[10] = yr * zr * (3.0 * a[2][1] + zr2 * a[2][2]);
    }

    private double energy(int i, int k) {
      // Electrostatic solvation energy of the permanent multipoles in their own GK reaction
      // potential.
      double esym =
          ci * ck * gc[1]
              - (uxi * (uxk * gux[2] + uyk * guy[2] + uzk * guz[2])
              + uyi * (uxk * gux[3] + uyk * guy[3] + uzk * guz[3])
              + uzi * (uxk * gux[4] + uyk * guy[4] + uzk * guz[4]));
      double ewi =
          ci * (uxk * gc[2] + uyk * gc[3] + uzk * gc[4])
              - ck * (uxi * gux[1] + uyi * guy[1] + uzi * guz[1])
              + ci
              * (qxxk * gc[5]
              + qyyk * gc[8]
              + qzzk * gc[10]
              + 2.0 * (qxyk * gc[6] + qxzk * gc[7] + qyzk * gc[9]))
              + ck
              * (qxxi * gqxx[1]
              + qyyi * gqyy[1]
              + qzzi * gqzz[1]
              + 2.0 * (qxyi * gqxy[1] + qxzi * gqxz[1] + qyzi * gqyz[1]))
              - uxi
              * (qxxk * gux[5]
              + qyyk * gux[8]
              + qzzk * gux[10]
              + 2.0 * (qxyk * gux[6] + qxzk * gux[7] + qyzk * gux[9]))
              - uyi
              * (qxxk * guy[5]
              + qyyk * guy[8]
              + qzzk * guy[10]
              + 2.0 * (qxyk * guy[6] + qxzk * guy[7] + qyzk * guy[9]))
              - uzi
              * (qxxk * guz[5]
              + qyyk * guz[8]
              + qzzk * guz[10]
              + 2.0 * (qxyk * guz[6] + qxzk * guz[7] + qyzk * guz[9]))
              + uxk
              * (qxxi * gqxx[2]
              + qyyi * gqyy[2]
              + qzzi * gqzz[2]
              + 2.0 * (qxyi * gqxy[2] + qxzi * gqxz[2] + qyzi * gqyz[2]))
              + uyk
              * (qxxi * gqxx[3]
              + qyyi * gqyy[3]
              + qzzi * gqzz[3]
              + 2.0 * (qxyi * gqxy[3] + qxzi * gqxz[3] + qyzi * gqyz[3]))
              + uzk
              * (qxxi * gqxx[4]
              + qyyi * gqyy[4]
              + qzzi * gqzz[4]
              + 2.0 * (qxyi * gqxy[4] + qxzi * gqxz[4] + qyzi * gqyz[4]))
              + qxxi
              * (qxxk * gqxx[5]
              + qyyk * gqxx[8]
              + qzzk * gqxx[10]
              + 2.0 * (qxyk * gqxx[6] + qxzk * gqxx[7] + qyzk * gqxx[9]))
              + qyyi
              * (qxxk * gqyy[5]
              + qyyk * gqyy[8]
              + qzzk * gqyy[10]
              + 2.0 * (qxyk * gqyy[6] + qxzk * gqyy[7] + qyzk * gqyy[9]))
              + qzzi
              * (qxxk * gqzz[5]
              + qyyk * gqzz[8]
              + qzzk * gqzz[10]
              + 2.0 * (qxyk * gqzz[6] + qxzk * gqzz[7] + qyzk * gqzz[9]))
              + 2.0
              * (qxyi
              * (qxxk * gqxy[5]
              + qyyk * gqxy[8]
              + qzzk * gqxy[10]
              + 2.0 * (qxyk * gqxy[6] + qxzk * gqxy[7] + qyzk * gqxy[9]))
              + qxzi
              * (qxxk * gqxz[5]
              + qyyk * gqxz[8]
              + qzzk * gqxz[10]
              + 2.0 * (qxyk * gqxz[6] + qxzk * gqxz[7] + qyzk * gqxz[9]))
              + qyzi
              * (qxxk * gqyz[5]
              + qyyk * gqyz[8]
              + qzzk * gqyz[10]
              + 2.0 * (qxyk * gqyz[6] + qxzk * gqyz[7] + qyzk * gqyz[9])));
      double ewk =
          ci * (uxk * gux[1] + uyk * guy[1] + uzk * guz[1])
              - ck * (uxi * gc[2] + uyi * gc[3] + uzi * gc[4])
              + ci
              * (qxxk * gqxx[1]
              + qyyk * gqyy[1]
              + qzzk * gqzz[1]
              + 2.0 * (qxyk * gqxy[1] + qxzk * gqxz[1] + qyzk * gqyz[1]))
              + ck
              * (qxxi * gc[5]
              + qyyi * gc[8]
              + qzzi * gc[10]
              + 2.0 * (qxyi * gc[6] + qxzi * gc[7] + qyzi * gc[9]))
              - uxi
              * (qxxk * gqxx[2]
              + qyyk * gqyy[2]
              + qzzk * gqzz[2]
              + 2.0 * (qxyk * gqxy[2] + qxzk * gqxz[2] + qyzk * gqyz[2]))
              - uyi
              * (qxxk * gqxx[3]
              + qyyk * gqyy[3]
              + qzzk * gqzz[3]
              + 2.0 * (qxyk * gqxy[3] + qxzk * gqxz[3] + qyzk * gqyz[3]))
              - uzi
              * (qxxk * gqxx[4]
              + qyyk * gqyy[4]
              + qzzk * gqzz[4]
              + 2.0 * (qxyk * gqxy[4] + qxzk * gqxz[4] + qyzk * gqyz[4]))
              + uxk
              * (qxxi * gux[5]
              + qyyi * gux[8]
              + qzzi * gux[10]
              + 2.0 * (qxyi * gux[6] + qxzi * gux[7] + qyzi * gux[9]))
              + uyk
              * (qxxi * guy[5]
              + qyyi * guy[8]
              + qzzi * guy[10]
              + 2.0 * (qxyi * guy[6] + qxzi * guy[7] + qyzi * guy[9]))
              + uzk
              * (qxxi * guz[5]
              + qyyi * guz[8]
              + qzzi * guz[10]
              + 2.0 * (qxyi * guz[6] + qxzi * guz[7] + qyzi * guz[9]))
              + qxxi
              * (qxxk * gqxx[5]
              + qyyk * gqyy[5]
              + qzzk * gqzz[5]
              + 2.0 * (qxyk * gqxy[5] + qxzk * gqxz[5] + qyzk * gqyz[5]))
              + qyyi
              * (qxxk * gqxx[8]
              + qyyk * gqyy[8]
              + qzzk * gqzz[8]
              + 2.0 * (qxyk * gqxy[8] + qxzk * gqxz[8] + qyzk * gqyz[8]))
              + qzzi
              * (qxxk * gqxx[10]
              + qyyk * gqyy[10]
              + qzzk * gqzz[10]
              + 2.0 * (qxyk * gqxy[10] + qxzk * gqxz[10] + qyzk * gqyz[10]))
              + 2.0
              * (qxyi
              * (qxxk * gqxx[6]
              + qyyk * gqyy[6]
              + qzzk * gqzz[6]
              + 2.0 * (qxyk * gqxy[6] + qxzk * gqxz[6] + qyzk * gqyz[6]))
              + qxzi
              * (qxxk * gqxx[7]
              + qyyk * gqyy[7]
              + qzzk * gqzz[7]
              + 2.0 * (qxyk * gqxy[7] + qxzk * gqxz[7] + qyzk * gqyz[7]))
              + qyzi
              * (qxxk * gqxx[9]
              + qyyk * gqyy[9]
              + qzzk * gqzz[9]
              + 2.0 * (qxyk * gqxy[9] + qxzk * gqxz[9] + qyzk * gqyz[9])));
      double e = esym + 0.5 * (ewi + ewk);
      double ei = 0.0;

      // Electrostatic solvation energy of the permanent multipoles in the
      // GK reaction potential of the induced dipoles.
      if (polarization != ParticleMeshEwald.Polarization.NONE) {
        double esymi =
            -uxi * (dxk * gux[2] + dyk * guy[2] + dzk * guz[2])
                - uyi * (dxk * gux[3] + dyk * guy[3] + dzk * guz[3])
                - uzi * (dxk * gux[4] + dyk * guy[4] + dzk * guz[4])
                - uxk * (dxi * gux[2] + dyi * guy[2] + dzi * guz[2])
                - uyk * (dxi * gux[3] + dyi * guy[3] + dzi * guz[3])
                - uzk * (dxi * gux[4] + dyi * guy[4] + dzi * guz[4]);
        double ewii =
            ci * (dxk * gc[2] + dyk * gc[3] + dzk * gc[4])
                - ck * (dxi * gux[1] + dyi * guy[1] + dzi * guz[1])
                - dxi
                * (qxxk * gux[5]
                + qyyk * gux[8]
                + qzzk * gux[10]
                + 2.0 * (qxyk * gux[6] + qxzk * gux[7] + qyzk * gux[9]))
                - dyi
                * (qxxk * guy[5]
                + qyyk * guy[8]
                + qzzk * guy[10]
                + 2.0 * (qxyk * guy[6] + qxzk * guy[7] + qyzk * guy[9]))
                - dzi
                * (qxxk * guz[5]
                + qyyk * guz[8]
                + qzzk * guz[10]
                + 2.0 * (qxyk * guz[6] + qxzk * guz[7] + qyzk * guz[9]))
                + dxk
                * (qxxi * gqxx[2]
                + qyyi * gqyy[2]
                + qzzi * gqzz[2]
                + 2.0 * (qxyi * gqxy[2] + qxzi * gqxz[2] + qyzi * gqyz[2]))
                + dyk
                * (qxxi * gqxx[3]
                + qyyi * gqyy[3]
                + qzzi * gqzz[3]
                + 2.0 * (qxyi * gqxy[3] + qxzi * gqxz[3] + qyzi * gqyz[3]))
                + dzk
                * (qxxi * gqxx[4]
                + qyyi * gqyy[4]
                + qzzi * gqzz[4]
                + 2.0 * (qxyi * gqxy[4] + qxzi * gqxz[4] + qyzi * gqyz[4]));
        double ewki =
            ci * (dxk * gux[1] + dyk * guy[1] + dzk * guz[1])
                - ck * (dxi * gc[2] + dyi * gc[3] + dzi * gc[4])
                - dxi
                * (qxxk * gqxx[2]
                + qyyk * gqyy[2]
                + qzzk * gqzz[2]
                + 2.0 * (qxyk * gqxy[2] + qxzk * gqxz[2] + qyzk * gqyz[2]))
                - dyi
                * (qxxk * gqxx[3]
                + qyyk * gqyy[3]
                + qzzk * gqzz[3]
                + 2.0 * (qxyk * gqxy[3] + qxzk * gqxz[3] + qyzk * gqyz[3]))
                - dzi
                * (qxxk * gqxx[4]
                + qyyk * gqyy[4]
                + qzzk * gqzz[4]
                + 2.0 * (qxyk * gqxy[4] + qxzk * gqxz[4] + qyzk * gqyz[4]))
                + dxk
                * (qxxi * gux[5]
                + qyyi * gux[8]
                + qzzi * gux[10]
                + 2.0 * (qxyi * gux[6] + qxzi * gux[7] + qyzi * gux[9]))
                + dyk
                * (qxxi * guy[5]
                + qyyi * guy[8]
                + qzzi * guy[10]
                + 2.0 * (qxyi * guy[6] + qxzi * guy[7] + qyzi * guy[9]))
                + dzk
                * (qxxi * guz[5]
                + qyyi * guz[8]
                + qzzi * guz[10]
                + 2.0 * (qxyi * guz[6] + qxzi * guz[7] + qyzi * guz[9]));
        ei = 0.5 * (esymi + 0.5 * (ewii + ewki));
      }

      if (i == k) {
        e *= 0.5;
        ei *= 0.5;
File Project Line
ffx/numerics/multipole/GKTensorGlobal.java Numerics 255
ffx/numerics/multipole/GKTensorQI.java Numerics 253
case QUADRUPOLE:
        quadrupoleIPotentialAtK(mI, 2);
        eK = multipoleEnergy(mK);
        quadrupoleKPotentialAtI(mK, 2);
        eI = multipoleEnergy(mI);
        return 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) {
    switch (multipoleOrder) {
      default:
      case MONOPOLE:
        return monopoleEnergyAndGradient(mI, mK, Gi, Gk, Ti, Tk);
      case DIPOLE:
        return dipoleEnergyAndGradient(mI, mK, Gi, Gk, Ti, Tk);
      case QUADRUPOLE:
        return 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);

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

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

    return 0.5 * (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 site I multipole.
    // Only torque on the site K dipole.
    multipoleIPotentialAtK(mI, 1);
    dipoleTorque(mK, Tk);

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

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

    return 0.5 * (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 pole due to site K multipole.
    // Only torque on the site I quadrupole.
    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 pole due to site I multipole.
    // Only torque on the site K quadrupole.
    multipoleIPotentialAtK(mI, 2);
    quadrupoleTorque(mK, Tk);

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

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

    return 0.5 * (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) {
    GK_TENSOR_MODE currentMode = mode;
    setMode(GK_TENSOR_MODE.BORN);
    generateTensor();
    double db = multipoleEnergy(mI, mK);
    setMode(currentMode);
    return db;
  }

  /**
   * 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) {
    switch (multipoleOrder) {
      default:
      case MONOPOLE:
        // 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);
        return 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.
        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.
        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);
        return 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.
        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.
        eI = polarizationEnergy(mI);
        return 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) {
    switch (multipoleOrder) {
      default:
      case MONOPOLE:
        return monopolePolarizationEnergyAndGradient(mI, mK, Gi, Ti, Tk);
      case DIPOLE:
        return dipolePolarizationEnergyAndGradient(mI, mK, mutualMask, Gi, Ti, Tk);
      case QUADRUPOLE:
        return 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.
   * @param Ti an array of {@link double} objects.
   * @param Tk an array of {@link double} objects.
   * @return a double.
   */
  public double monopolePolarizationEnergyAndGradient(
      PolarizableMultipole mI, PolarizableMultipole mK, double[] Gi, double[] Ti, double[] Tk) {

    // 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);

    Gi[0] *= 0.5;
    Gi[1] *= 0.5;
    Gi[2] *= 0.5;

    // Total polarization energy.
    return 0.5 * (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);
    // Find the potential at K due to the averaged induced dipole at i.
    dipoleIPotentialAtK(mI.sx, mI.sy, mI.sz, 2);
    dipoleTorque(mK, Tk);

    // Find the GK induced dipole potential of site I at site K.
    dipoleIPotentialAtK(mI.ux, mI.uy, mI.uz, 3);
    // Energy of permanent multipole K in the field of induced dipole I.
    eK += 0.5 * multipoleEnergy(mK);
    double[] G = new double[3];
    multipoleGradient(mK, G);
    Gi[0] -= G[0];
    Gi[1] -= G[1];
    Gi[2] -= G[2];
    multipoleTorque(mK, Tk);

    // Find the permanent multipole potential and derivatives at site i.
    dipoleKPotentialAtI(mK.dx, mK.dy, mK.dz, 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);
    // Find the potential at I due to the averaged induced dipole at k.
    dipoleKPotentialAtI(mK.sx, mK.sy, mK.sz, 2);
    dipoleTorque(mI, Ti);

    // Find the GK induced dipole potential of site K at site I.
    dipoleKPotentialAtI(mK.ux, mK.uy, mK.uz, 3);
    // Energy of permanent multipole I in the field of induced dipole K.
    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 energy = 0.5 * (eI + eK);
    Gi[0] *= 0.5;
    Gi[1] *= 0.5;
    Gi[2] *= 0.5;
    Ti[0] *= 0.5;
    Ti[1] *= 0.5;
    Ti[2] *= 0.5;
    Tk[0] *= 0.5;
    Tk[1] *= 0.5;
    Tk[2] *= 0.5;

    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);

    Gi[0] *= 0.5;
    Gi[1] *= 0.5;
    Gi[2] *= 0.5;

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

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

    // Total polarization energy.
    return 0.5 * (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) {
    GK_TENSOR_MODE currentMode = mode;
    setMode(GK_TENSOR_MODE.BORN);
    generateTensor();
    double db = 2.0 * polarizationEnergy(mI, mK);
    setMode(currentMode);
    return db;
  }

  /**
   * 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 db;
  }
File Project Line
ffx/potential/parsers/XPHFilter.java Potential 106
ffx/potential/parsers/XYZFilter.java Potential 119
}

  /**
   * 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 xyzFile = activeMolecularAssembly.getFile();
    int nAtoms = activeMolecularAssembly.getAtomArray().length;
    Pattern crystInfoPattern =
        Pattern.compile(
            "^ *(?:[0-9]+\\.[0-9]+ +){3}(?:-?[0-9]+\\.[0-9]+ +){2}(?:-?[0-9]+\\.[0-9]+) *$");

    try (BufferedReader br = new BufferedReader(new FileReader(xyzFile))) {
      String line = br.readLine();
      int nSnaps = 0;
      // For each header line, will read either X or X+1 lines, where X is the number of atoms.
      while (line != null) {
        assert parseInt(line.trim().split("\\s+")[0]) == nAtoms;
        // Read either the crystal information *or* the first line of the snapshot.
        line = br.readLine();
        Matcher m = crystInfoPattern.matcher(line);
        if (m.matches()) {
          // If this is crystal information, move onto the first line of the snapshot.
          br.readLine();
        }
        // Read lines 2-X of the XYZ.
        for (int i = 1; i < nAtoms; i++) {
          br.readLine();
        }

        ++nSnaps;
        line = br.readLine();
      }
      return nSnaps;
    } catch (Exception ex) {
      logger.log(
          Level.WARNING, String.format(" Exception reading trajectory file %s: %s", xyzFile, ex));
      return 1;
    }
  }

  /** {@inheritDoc} */
  @Override
  public OptionalDouble getLastReadLambda() {
    String[] toks = remarkLine.split("\\s+");
    int nToks = toks.length;
    for (int i = 0; i < (nToks - 1); i++) {
      if (toks[i].equals("Lambda:")) {
        return OptionalDouble.of(Double.parseDouble(toks[i + 1]));
      }
    }
    return OptionalDouble.empty();
  }

  @Override
  public String[] getRemarkLines() {
    return new String[] {remarkLine};
  }

  @Override
  public int getSnapshot() {
    return snapShot;
  }

  /**
   * {@inheritDoc}
   *
   * <p>Parse the XYZ File
   */
  @Override
  public boolean readFile() {
    File xyzFile = activeMolecularAssembly.getFile();

    if (forceField == null) {
      logger.warning(format(" No force field is associated with %s.", xyzFile.toString()));
      return false;
    }
    try (BufferedReader br = new BufferedReader(new FileReader(xyzFile))) {
      String data = br.readLine();
      // Read blank lines at the top of the file
      while (data != null && data.trim().equals("")) {
        data = br.readLine();
      }
      if (data == null) {
        return false;
      }
      String[] tokens = data.trim().split(" +", 2);
      int numberOfAtoms = parseInt(tokens[0]);
      if (numberOfAtoms < 1) {
        return false;
      }
      if (tokens.length == 2) {
        getActiveMolecularSystem().setName(tokens[1]);
      }
      logger.info(format(" Opening %s with %d atoms\n", xyzFile.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<>();
      int[] c = new int[2];
      for (int a1 = 1; a1 <= numberOfAtoms; a1++) {
        int j = -1;
        while (j < 7 && bonds[a1 - 1][++j] > 0) {
          int a2 = bonds[a1 - 1][j];
          if (a1 < a2) {
            if (a2 > numberOfAtoms) {
              logger.warning(
                  format(
                      " Check the bond between %d and %d in %s.",
                      a1, a2, activeMolecularAssembly.getFile().getName()));
              return false;
            }
            // Check for bidirectional connection
            boolean bidirectional = false;
            int k = -1;
            while (k < 7 && bonds[a2 - 1][++k] > 0) {
              int a3 = bonds[a2 - 1][k];
              if (a3 == a1) {
                bidirectional = true;
                break;
              }
            }
            if (!bidirectional) {
              logger.warning(
                  format(
                      " Check the bond between %d and %d in %s.",
                      a1, a2, activeMolecularAssembly.getFile().getName()));
              return false;
            }
            Atom atom1 = atomList.get(a1 - 1);
            Atom atom2 = atomList.get(a2 - 1);
            if (atom1 == null || atom2 == null) {
              logger.warning(
                  format(
                      " Check the bond between %d and %d in %s.",
                      a1, a2, activeMolecularAssembly.getFile().getName()));
              return false;
            }
            Bond bond = new Bond(atom1, atom2);
            BondType bondType = forceField.getBondType(atom1.getAtomType(), atom2.getAtomType());
            if (bondType == null) {
              logNoBondType(atom1, atom2, forceField);
            } else {
              bond.setBondType(bondType);
            }
            bondList.add(bond);
          }
        }
      }
      return true;
    } catch (IOException e) {
      logger.severe(e.toString());
    }
    return false;
  }

  /** {@inheritDoc} */
  @Override
  public boolean readNext() {
    return readNext(false);
  }

  /**
   * {@inheritDoc}
   *
   * <p>Reads the next snap-shot of an archive into the activeMolecularAssembly. After calling this
   * function, a BufferedReader will remain open until the <code>close</code> method is called.
   */
  @Override
  public boolean readNext(boolean resetPosition) {
    return readNext(resetPosition, true);
  }

  /**
   * {@inheritDoc}
   *
   * <p>Reads the next snap-shot of an archive into the activeMolecularAssembly. After calling this
   * function, a BufferedReader will remain open until the <code>close</code> method is called.
   */
  @Override
  public boolean readNext(boolean resetPosition, boolean print) {
    return readNext(resetPosition, print, true);
  }

  /**
   * Reads the next snap-shot of an archive into the activeMolecularAssembly. After calling this
   * function, a BufferedReader will remain open until the <code>close</code> method is called.
   */
  public boolean readNext(boolean resetPosition, boolean print, boolean parse) {
    try {
      String data;
      Atom[] atoms = activeMolecularAssembly.getAtomArray();
      int nSystem = atoms.length;

      if (bufferedReader == null && !resetPosition) {
        bufferedReader = new BufferedReader(new FileReader(currentFile));
        // Read past the first N + 1 lines that begin with an integer.
        for (int i = 0; i < nSystem + 1; i++) {
          data = bufferedReader.readLine();
          while (!firstTokenIsInteger(data)) {
            data = bufferedReader.readLine();
          }
        }
        snapShot = 1;
      } else if (resetPosition) {
        // Reset the reader to the beginning of the file. Do not skip reading the first entry if resetPostion is true.
        bufferedReader = new BufferedReader(new FileReader(currentFile));
        snapShot = 0;
      }

      snapShot++;

      data = bufferedReader.readLine();
      // Read past blank lines
      while (data != null && data.trim().equals("")) {
        data = bufferedReader.readLine();
      }
      if (data == null) {
        return false;
      }

      // Read Past ESV
      if(data.contains("ESV")) {
File Project Line
ffx/potential/utils/Loop.java Potential 2617
ffx/potential/utils/Loop.java Potential 2750
double[] hz3 = coordsArray[HZ3.getIndex() - 1];
            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];
            arraycopy(determineIntxyz(ca, 1.54, n, 109.5, c, 107.8, 1), 0, cb, 0, 3);
            coordsArray = fillCoordsArray(CB, coordsArray, cb);

            arraycopy(
                determineIntxyz(cb, dCG_CB, ca, dCG_CB_CA, n, rotScale * 180.0, 0),
                0,
                cg,
                0,
                3); // CG
            coordsArray = fillCoordsArray(CG, coordsArray, cg);

            arraycopy(determineIntxyz(cg, dCD_CG, cb, dCD_CG_CB, ca, 180.0, 0), 0, cd, 0, 3); // CD
            coordsArray = fillCoordsArray(CD, coordsArray, cd);

            arraycopy(determineIntxyz(cd, dCE_CD, cg, dCE_CD_CG, cb, 180.0, 0), 0, ce, 0, 3); // CE
            coordsArray = fillCoordsArray(CE, coordsArray, ce);

            arraycopy(determineIntxyz(ce, dNZ_CE, cd, dNZ_CE_CD, cg, 180.0, 0), 0, nz, 0, 3); // NZ
            coordsArray = fillCoordsArray(NZ, coordsArray, nz);

            arraycopy(
                determineIntxyz(cb, dHB_CB, ca, dHB_CB_CA, cg, 109.4, 1), 0, hb2, 0, 3); // HB2
            coordsArray = fillCoordsArray(HB2, coordsArray, hb2);

            arraycopy(
                determineIntxyz(cb, dHB_CB, ca, dHB_CB_CA, cg, 109.4, -1), 0, hb3, 0, 3); // HB3
            coordsArray = fillCoordsArray(HB3, coordsArray, hb3);

            arraycopy(
                determineIntxyz(cg, dHG_CG, cb, dHG_CG_CB, cd, 109.4, 1), 0, hg2, 0, 3); // HG2
            coordsArray = fillCoordsArray(HG2, coordsArray, hg2);

            arraycopy(
                determineIntxyz(cg, dHG_CG, cb, dHG_CG_CB, cd, 109.4, -1), 0, hg3, 0, 3); // HG3
            coordsArray = fillCoordsArray(HG3, coordsArray, hg3);

            arraycopy(
                determineIntxyz(cd, dHD_CD, cg, dHD_CD_CG, ce, 109.4, 1), 0, hd2, 0, 3); // HD2
            coordsArray = fillCoordsArray(HD2, coordsArray, hd2);

            arraycopy(
                determineIntxyz(cd, dHD_CD, cg, dHD_CD_CG, ce, 109.4, -1), 0, hd3, 0, 3); // HD3
            coordsArray = fillCoordsArray(HD3, coordsArray, hd3);

            arraycopy(
                determineIntxyz(ce, dHE_CE, cd, dHE_CE_CD, nz, 108.8, 1), 0, he2, 0, 3); // HE2
            coordsArray = fillCoordsArray(HE2, coordsArray, he2);

            arraycopy(
                determineIntxyz(ce, dHE_CE, cd, dHE_CE_CD, nz, 108.8, -1), 0, he3, 0, 3); // HE3
            coordsArray = fillCoordsArray(HE3, coordsArray, he3);

            arraycopy(
                determineIntxyz(nz, dHZ_NZ, ce, dHZ_NZ_CE, cd, 180.0, 0), 0, hz1, 0, 3); // HZ1
            coordsArray = fillCoordsArray(HZ1, coordsArray, hz1);

            arraycopy(
                determineIntxyz(nz, dHZ_NZ, ce, dHZ_NZ_CE, hz1, 109.5, 1), 0, hz2, 0, 3); // HZ2
            coordsArray = fillCoordsArray(HZ2, coordsArray, hz2);
File Project Line
ffx/potential/nonbonded/pme/InducedDipoleFieldRegion.java Potential 374
ffx/potential/nonbonded/pme/PCGSolver.java Potential 978
final double bn2 = (3.0 * bn1 + an1 * exp2a) * rr2;
            double scale3 = 1.0;
            double scale5 = 1.0;
            double damp = pdi * pdk;
            final double pgamma = min(pti, ptk);
            final double rdamp = r * damp;
            damp = -pgamma * rdamp * rdamp * rdamp;
            if (damp > -50.0) {
              final double expdamp = exp(damp);
              scale3 = 1.0 - expdamp;
              scale5 = 1.0 - expdamp * (1.0 - damp);
            }
            double rr3 = rr1 * rr2;
            double rr5 = 3.0 * rr3 * rr2;
            rr3 *= (1.0 - scale3);
            rr5 *= (1.0 - scale5);
            final double xr = dx[0];
            final double yr = dx[1];
            final double zr = dx[2];
            final double[] dipolek = ind[k];
            final double ukx = dipolek[0];
            final double uky = dipolek[1];
            final double ukz = dipolek[2];
            final double ukr = ukx * xr + uky * yr + ukz * zr;
            final double bn2ukr = bn2 * ukr;
            final double fimx = -bn1 * ukx + bn2ukr * xr;
            final double fimy = -bn1 * uky + bn2ukr * yr;
            final double fimz = -bn1 * ukz + bn2ukr * zr;
            final double rr5ukr = rr5 * ukr;
            final double fidx = -rr3 * ukx + rr5ukr * xr;
            final double fidy = -rr3 * uky + rr5ukr * yr;
            final double fidz = -rr3 * ukz + rr5ukr * zr;
            fx += (fimx - fidx);
            fy += (fimy - fidy);
            fz += (fimz - fidz);
            final double[] dipolepk = indCR[k];
            final double pkx = dipolepk[0];
            final double pky = dipolepk[1];
            final double pkz = dipolepk[2];
            final double pkr = pkx * xr + pky * yr + pkz * zr;
            final double bn2pkr = bn2 * pkr;
            final double pimx = -bn1 * pkx + bn2pkr * xr;
            final double pimy = -bn1 * pky + bn2pkr * yr;
            final double pimz = -bn1 * pkz + bn2pkr * zr;
            final double rr5pkr = rr5 * pkr;
            final double pidx = -rr3 * pkx + rr5pkr * xr;
            final double pidy = -rr3 * pky + rr5pkr * yr;
            final double pidz = -rr3 * pkz + rr5pkr * zr;
            px += (pimx - pidx);
            py += (pimy - pidy);
            pz += (pimz - pidz);
            final double uir = uix * xr + uiy * yr + uiz * zr;
            final double bn2uir = bn2 * uir;
            final double fkmx = -bn1 * uix + bn2uir * xr;
            final double fkmy = -bn1 * uiy + bn2uir * yr;
            final double fkmz = -bn1 * uiz + bn2uir * zr;
            final double rr5uir = rr5 * uir;
            final double fkdx = -rr3 * uix + rr5uir * xr;
            final double fkdy = -rr3 * uiy + rr5uir * yr;
            final double fkdz = -rr3 * uiz + rr5uir * zr;
            field.add(threadID, k, fkmx - fkdx, fkmy - fkdy, fkmz - fkdz);
            final double pir = pix * xr + piy * yr + piz * zr;
            final double bn2pir = bn2 * pir;
            final double pkmx = -bn1 * pix + bn2pir * xr;
            final double pkmy = -bn1 * piy + bn2pir * yr;
            final double pkmz = -bn1 * piz + bn2pir * zr;
            final double rr5pir = rr5 * pir;
            final double pkdx = -rr3 * pix + rr5pir * xr;
            final double pkdy = -rr3 * piy + rr5pir * yr;
            final double pkdz = -rr3 * piz + rr5pir * zr;
            fieldCR.add(threadID, k, pkmx - pkdx, pkmy - pkdy, pkmz - pkdz);
          }
          field.add(threadID, i, fx, fy, fz);
          fieldCR.add(threadID, i, px, py, pz);
        }

        // Loop over symmetry mates.
        List<SymOp> symOps = crystal.spaceGroup.symOps;
        int nSymm = symOps.size();
        for (int iSymm = 1; iSymm < nSymm; iSymm++) {
          SymOp symOp = crystal.spaceGroup.getSymOp(iSymm);
          crystal.getTransformationOperator(symOp, transOp);
          lists = realSpaceLists[iSymm];
File Project Line
ffx/potential/bonded/RestraintTorsion.java Potential 37
ffx/potential/bonded/Torsion.java Potential 285
}

    @Override
    public double energy(boolean gradient, int threadID, AtomicDoubleArray3D grad, AtomicDoubleArray3D lambdaGrad) {
        energy = 0.0;
        value = 0.0;
        dEdL = 0.0;
        var atomA = atoms[0];
        var atomB = atoms[1];
        var atomC = atoms[2];
        var atomD = atoms[3];
        var va = atomA.getXYZ();
        var vb = atomB.getXYZ();
        var vc = atomC.getXYZ();
        var vd = atomD.getXYZ();
        var vba = vb.sub(va);
        var vcb = vc.sub(vb);
        var vdc = vd.sub(vc);
        var vt = vba.X(vcb);
        var vu = vcb.X(vdc);
        var rt2 = vt.length2();
        var ru2 = vu.length2();
        var rtru2 = rt2 * ru2;
        if (rtru2 != 0.0) {
            var rr = sqrt(rtru2);
            var rcb = vcb.length();
            var cosine = vt.dot(vu) / rr;
            var sine = vcb.dot(vt.X(vu)) / (rcb * rr);
            value = toDegrees(acos(cosine));
            if (sine < 0.0) {
                value = -value;
            }
            var amp = torsionType.amplitude;
            var tsin = torsionType.sine;
            var tcos = torsionType.cosine;
            energy = amp[0] * (1.0 + cosine * tcos[0] + sine * tsin[0]);
            var dedphi = amp[0] * (cosine * tsin[0] - sine * tcos[0]);
            var cosprev = cosine;
            var sinprev = sine;
            var n = torsionType.terms;
            for (int i = 1; i < n; i++) {
                var cosn = cosine * cosprev - sine * sinprev;
                var sinn = sine * cosprev + cosine * sinprev;
                var phi = 1.0 + cosn * tcos[i] + sinn * tsin[i];
                var dphi = (1.0 + i) * (cosn * tsin[i] - sinn * tcos[i]);
                energy = energy + amp[i] * phi;
                dedphi = dedphi + amp[i] * dphi;
                cosprev = cosn;
                sinprev = sinn;
            }
            energy = units * energy * lambda;
            dEdL = units * energy;
            if (gradient || lambdaTerm) {
                dedphi = units * dedphi;
                var vca = vc.sub(va);
                var vdb = vd.sub(vb);
                var dedt = vt.X(vcb).scaleI(dedphi / (rt2 * rcb));
                var dedu = vu.X(vcb).scaleI(-dedphi / (ru2 * rcb));
                var ga = dedt.X(vcb);
                var gb = vca.X(dedt).addI(dedu.X(vdc));
                var gc = dedt.X(vba).addI(vdb.X(dedu));
                var gd = dedu.X(vcb);
                int ia = atomA.getIndex() - 1;
                int ib = atomB.getIndex() - 1;
                int ic = atomC.getIndex() - 1;
                int id = atomD.getIndex() - 1;
                if (lambdaTerm) {
                    lambdaGrad.add(threadID, ia, ga);
                    lambdaGrad.add(threadID, ib, gb);
                    lambdaGrad.add(threadID, ic, gc);
                    lambdaGrad.add(threadID, id, gd);
                }
                if (gradient) {
                    grad.add(threadID, ia, ga.scaleI(lambda));
                    grad.add(threadID, ib, gb.scaleI(lambda));
                    grad.add(threadID, ic, gc.scaleI(lambda));
                    grad.add(threadID, id, gd.scaleI(lambda));
                }
            }
        }

        return energy;
    }
File Project Line
ffx/numerics/multipole/GKTensorGlobal.java Numerics 785
ffx/numerics/multipole/GKTensorQI.java Numerics 783
protected void setR(double dx, double dy, double dz, double ai, double aj) {
    setR(dx, dy, dz);
    this.ai = ai;
    this.aj = aj;
    expTerm = exp(-r2 / (gc * ai * aj));
    f = sqrt(r2 + ai * aj * expTerm);
  }

  /**
   * Set the "mode" for the tensor (either POTENTIAL or BORN).
   *
   * @param mode The mode for tensor generation.
   */
  protected void setMode(GK_TENSOR_MODE mode) {
    this.mode = mode;
  }

  /**
   * Generate source terms for the Kirkwood version of the Challacombe et al. recursion.
   */
  protected void source(double[] work) {
    int multipoleOrder = this.multipoleOrder.getOrder();

    if (mode == GK_TENSOR_MODE.POTENTIAL) {
      // Prepare the GK Potential tensor.
      for (int n = 0; n <= order; n++) {
        an0[n] = an0(n);
        fn[n] = fn(n);
      }
      for (int n = 0; n <= order; n++) {
        if (n < multipoleOrder) {
          work[n] = 0.0;
        } else {
          work[n] = anm(multipoleOrder, n - multipoleOrder);
        }
      }
    } else {
      // Prepare the GK Born-chain rule tensor.
      for (int n = 0; n <= order; n++) {
        an0[n] = an0(n);
        fn[n] = fn(n);
        bn[n] = bn(n);
      }
      // Only up to order - 1.
      for (int n = 0; n <= order - 1; n++) {
        if (n < multipoleOrder) {
          work[n] = 0.0;
        } else {
          work[n] = bnm(multipoleOrder, n - multipoleOrder);
        }
      }
    }
  }

  /**
   * Compute the potential auxiliary function for a multipole of order n.
   *
   * @param n Multipole order.
   * @return The potential auxiliary function for a multipole of order n.
   */
  protected double an0(int n) {
    return kirkwoodSource[n] / pow(f, 2 * n + 1);
  }

  /**
   * Compute the mth potential gradient auxiliary function for a multipole of order n.
   *
   * @param n Multipole order.
   * @param m Mth potential gradient auxiliary function.
   * @return Returns the mth potential gradient auxiliary function for a multipole of order n.
   */
  protected double anm(int n, int m) {
    if (m == 0) {
      return an0[n];
    }
    var ret = 0.0;
    var coef = anmc[m];
    for (int i = 1; i <= m; i++) {
      ret += coef[i - 1] * fn[i] * anm(n + 1, m - i);
    }
    return ret;
  }

  /**
   * Compute the derivative with respect to a Born radius of the mth potential gradient auxiliary
   * function for a multipole of order n.
   *
   * @param n Multipole order.
   * @param m Mth potential gradient auxiliary function.
   * @return Returns the derivative with respect to a Born radius of the mth potential gradient
   *     auxiliary function for a multipole of order n.
   */
  protected double bnm(int n, int m) {
    if (m == 0) {
      // return bn(0) * an0(n + 1);
      return bn[0] * an0[n + 1];
    }
    var ret = 0.0;
    var coef = anmc[m];
    for (int i = 1; i <= m; i++) {
      ret += coef[i - 1] * bn[i] * anm(n + 1, m - i);
      ret += coef[i - 1] * fn[i] * bnm(n + 1, m - i);
    }
    return ret;
  }

  /**
   * Returns nth value of the function f, which are chain rule terms from differentiating zeroth
   * order auxiliary functions (an0) with respect to x, y or z.
   *
   * @param n Multipole order.
   * @return Returns the nth value of the function f.
   */
  protected double fn(int n) {
    switch (n) {
      case 0:
        return f;
      case 1:
        return 1.0 - expTerm / gc;
      default:
        var gcAiAj = gc * ai * aj;
        var f2 = 2.0 * expTerm / (gc * gcAiAj);
        var fr = -2.0 / gcAiAj;
        return pow(fr, n - 2) * f2;
    }
  }

  /**
   * Returns nth value of the function b, which are chain rule terms from differentiating zeroth
   * order auxiliary functions (an0) with respect to Ai or Aj.
   *
   * @param n Multipole order.
   * @return Returns the nth value of the function f.
   */
  protected double bn(int n) {
    var gcAiAj = gc * ai * aj;
    var ratio = -r2 / gcAiAj;
    switch (n) {
      case 0:
        return 0.5 * expTerm * (1.0 - ratio);
      case 1:
        return -r2 * expTerm / (gcAiAj * gcAiAj);
      default:
        var b2 = 2.0 * expTerm / (gcAiAj * gcAiAj) * (-ratio - 1.0);
        var br = 2.0 / (gcAiAj * ai * aj);
        var f2 = 2.0 / (gc * gcAiAj) * expTerm;
        var fr = -2.0 / (gcAiAj);
        return (n - 2) * pow(fr, n - 3) * br * f2 + pow(fr, n - 2) * b2;
    }
  }

  /**
   * Return coefficients needed when taking derivatives of auxiliary functions.
   *
   * @param n Multipole order.
   * @return Returns coefficients needed when taking derivatives of auxiliary functions.
   */
  protected static double[] anmc(int n) {
File Project Line
ffx/potential/parsers/XPHFilter.java Potential 799
ffx/potential/parsers/XYZFilter.java Potential 731
bw.write("\n");
        }
      } catch (IOException e) {
        String message =
            format(
                " There was an unexpected error writing to %s.",
                getActiveMolecularSystem().toString());
        logger.log(Level.WARNING, message, e);
        return false;
      }
    } catch (IOException e) {
      String message =
          format(
              " There was an unexpected error writing to %s.",
              getActiveMolecularSystem().toString());
      logger.log(Level.WARNING, message, e);
      return false;
    }
    return true;
  }

  /**
   * writeFileAsP1
   *
   * @param saveFile a {@link File} object.
   * @param append a boolean.
   * @param crystal a {@link Crystal} object.
   * @return a boolean.
   */
  public boolean writeFileAsP1(File saveFile, boolean append, Crystal crystal) {
    return writeFileAsP1(saveFile, append, crystal, null);
  }

  /**
   * writeFileAsP1
   *
   * @param saveFile a {@link File} object.
   * @param append a boolean.
   * @param crystal a {@link Crystal} object.
   * @param extraLines Additional lines to print in the header.
   * @return a boolean.
   */
  public boolean writeFileAsP1(
      File saveFile, boolean append, Crystal crystal, String[] extraLines) {
    if (saveFile == null) {
      return false;
    }

    File newFile = saveFile;
    if (!append) {
      newFile = version(saveFile);
    }
    activeMolecularAssembly.setFile(newFile);
    activeMolecularAssembly.setName(newFile.getName());

    try (FileWriter fw = new FileWriter(newFile, append && newFile.exists());
        BufferedWriter bw = new BufferedWriter(fw)) {
      int nSymm = crystal.spaceGroup.symOps.size();
      // XYZ File First Line
      int numberOfAtoms = activeMolecularAssembly.getAtomList().size() * nSymm;
      StringBuilder sb =
          new StringBuilder(format("%7d  %s", numberOfAtoms, activeMolecularAssembly.toString()));
      if (extraLines != null) {
        for (String line : extraLines) {
          line = line.replaceAll("\n", " ");
          sb.append(" ").append(line);
        }
      }
      String output = sb.append("\n").toString();
      bw.write(output);

      if (!crystal.aperiodic()) {
        Crystal uc = crystal.getUnitCell();
        String params = format("%14.8f%14.8f%14.8f%14.8f%14.8f%14.8f\n",
            uc.a, uc.b, uc.c, uc.alpha, uc.beta, uc.gamma);
        bw.write(params);
      }

      Atom a2;
      StringBuilder line;
      StringBuilder[] lines = new StringBuilder[numberOfAtoms];
      // XYZ File Atom Lines
      Atom[] atoms = activeMolecularAssembly.getAtomArray();
      double[] xyz = new double[3];
      for (int iSym = 0; iSym < nSymm; iSym++) {
        SymOp symOp = crystal.spaceGroup.getSymOp(iSym);
        int indexOffset = iSym * atoms.length;
        for (Atom a : atoms) {
          int index = a.getIndex() + indexOffset;
          String id = a.getAtomType().name;
          if (vdwH) {
            a.getRedXYZ(xyz);
          } else {
            xyz[0] = a.getX();
            xyz[1] = a.getY();
            xyz[2] = a.getZ();
          }
          crystal.applySymOp(xyz, xyz, symOp);
          int type = a.getType();
          line =
              new StringBuilder(
                  format("%7d %3s%14.8f%14.8f%14.8f%6d", index, id, xyz[0], xyz[1], xyz[2], type));
          for (Bond b : a.getBonds()) {
            a2 = b.get1_2(a);
            line.append(format("%8d", a2.getIndex() + indexOffset));
          }
          lines[index - 1] = line.append("\n");
        }
      }
      try {
        for (int i = 0; i < numberOfAtoms; i++) {
          bw.write(lines[i].toString());
        }
      } catch (IOException e) {
        String message =
            format(
                " There was an unexpected error writing to %s.",
                getActiveMolecularSystem().toString());
        logger.log(Level.WARNING, message, e);
        return false;
      }
    } catch (IOException e) {
      String message =
          format(
              " There was an unexpected error writing to %s.",
              getActiveMolecularSystem().toString());
      logger.log(Level.WARNING, message, e);
      return false;
    }
    return true;
  }
}
File Project Line
ffx/potential/bonded/RotamerLibrary.java Potential 2265
ffx/potential/bonded/RotamerLibrary.java Potential 2336
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);
File Project Line
edu/rit/pj/WorkerLongForLoop.java Parallel Java 298
edu/rit/pj/WorkerLongStrideForLoop.java Parallel Java 301
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());
File Project Line
ffx/crystal/SpaceGroupDefinitions.java Crystal 6706
ffx/crystal/SpaceGroupDefinitions.java Crystal 7036
"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));
File Project Line
ffx/crystal/SpaceGroupDefinitions.java Crystal 6834
ffx/crystal/SpaceGroupDefinitions.java Crystal 7658
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),
File Project Line
ffx/potential/nonbonded/pme/InducedDipoleFieldRegion.java Potential 516
ffx/potential/nonbonded/pme/PCGSolver.java Potential 1122
final double bn2 = (3.0 * bn1 + an1 * exp2a) * rr2;
              double scale3 = 1.0;
              double scale5 = 1.0;
              double damp = pdi * pdk;
              final double pgamma = min(pti, ptk);
              final double rdamp = r * damp;
              damp = -pgamma * rdamp * rdamp * rdamp;
              if (damp > -50.0) {
                final double expdamp = exp(damp);
                scale3 = 1.0 - expdamp;
                scale5 = 1.0 - expdamp * (1.0 - damp);
              }
              double rr3 = rr1 * rr2;
              double rr5 = 3.0 * rr3 * rr2;
              rr3 *= (1.0 - scale3);
              rr5 *= (1.0 - scale5);
              final double xr = dx[0];
              final double yr = dx[1];
              final double zr = dx[2];
              final double[] dipolek = inds[k];
              final double ukx = dipolek[0];
              final double uky = dipolek[1];
              final double ukz = dipolek[2];
              final double[] dipolepk = indCRs[k];
              final double pkx = dipolepk[0];
              final double pky = dipolepk[1];
              final double pkz = dipolepk[2];
              final double ukr = ukx * xr + uky * yr + ukz * zr;
              final double bn2ukr = bn2 * ukr;
              final double fimx = -bn1 * ukx + bn2ukr * xr;
              final double fimy = -bn1 * uky + bn2ukr * yr;
              final double fimz = -bn1 * ukz + bn2ukr * zr;
              final double rr5ukr = rr5 * ukr;
              final double fidx = -rr3 * ukx + rr5ukr * xr;
              final double fidy = -rr3 * uky + rr5ukr * yr;
              final double fidz = -rr3 * ukz + rr5ukr * zr;
              fx += selfScale * (fimx - fidx);
              fy += selfScale * (fimy - fidy);
              fz += selfScale * (fimz - fidz);
              final double pkr = pkx * xr + pky * yr + pkz * zr;
              final double bn2pkr = bn2 * pkr;
              final double pimx = -bn1 * pkx + bn2pkr * xr;
              final double pimy = -bn1 * pky + bn2pkr * yr;
              final double pimz = -bn1 * pkz + bn2pkr * zr;
              final double rr5pkr = rr5 * pkr;
              final double pidx = -rr3 * pkx + rr5pkr * xr;
              final double pidy = -rr3 * pky + rr5pkr * yr;
              final double pidz = -rr3 * pkz + rr5pkr * zr;
              px += selfScale * (pimx - pidx);
              py += selfScale * (pimy - pidy);
              pz += selfScale * (pimz - pidz);
              final double uir = uix * xr + uiy * yr + uiz * zr;
              final double bn2uir = bn2 * uir;
              final double fkmx = -bn1 * uix + bn2uir * xr;
              final double fkmy = -bn1 * uiy + bn2uir * yr;
              final double fkmz = -bn1 * uiz + bn2uir * zr;
              final double rr5uir = rr5 * uir;
              final double fkdx = -rr3 * uix + rr5uir * xr;
              final double fkdy = -rr3 * uiy + rr5uir * yr;
              final double fkdz = -rr3 * uiz + rr5uir * zr;
              double xc = selfScale * (fkmx - fkdx);
              double yc = selfScale * (fkmy - fkdy);
              double zc = selfScale * (fkmz - fkdz);
              double kx = (xc * transOp[0][0] + yc * transOp[1][0] + zc * transOp[2][0]);
File Project Line
ffx/potential/nonbonded/ReciprocalSpace.java Potential 1606
ffx/potential/nonbonded/ReciprocalSpace.java Potential 1916
if (lbZ <= k && k <= ubZ) {
          atomContributes = true;
          break;
        }
      }

      if (!atomContributes) {
        return;
      }

      splineCount[threadIndex]++;

      // Add atom n to the list for the current symmOp and thread.
      int index = gridAtomCount[iSymm][threadIndex];
      gridAtomList[iSymm][threadIndex][index] = iAtom;
      gridAtomCount[iSymm][threadIndex]++;

      // Convert Cartesian multipoles in the global frame to fractional multipoles.
      final double[] gm = globalMultipoles[iSymm][iAtom];
      final double[] fm = fracMPole;

      // Charge.
      fm[0] = gm[0];

      // Dipole.
      for (int j = 1; j < 4; j++) {
        fm[j] = 0.0;
        for (int k = 1; k < 4; k++) {
          fm[j] = fm[j] + transformMultipoleMatrix[j][k] * gm[k];
        }
      }

      // Quadrupole.
      for (int j = 4; j < 10; j++) {
        fm[j] = 0.0;
        for (int k = 4; k < 7; k++) {
          fm[j] = fm[j] + transformMultipoleMatrix[j][k] * gm[k];
        }
        for (int k = 7; k < 10; k++) {
          fm[j] = fm[j] + transformMultipoleMatrix[j][k] * 2.0 * gm[k];
        }
        // Fractional quadrupole components are pre-multiplied by a
        // factor of 1/3 that arises in their potential.
        fm[j] = fm[j] / 3.0;
      }

      arraycopy(fm, 0, fracMultipoles[iSymm][iAtom], 0, 10);

      // Some atoms are not used during Lambda dynamics.
      if (use != null && !use[iAtom]) {
        return;
      }

      final double[][] splx = bSplines.splineX[iSymm][iAtom];
      final double[][] sply = bSplines.splineY[iSymm][iAtom];
      final double[][] splz = bSplines.splineZ[iSymm][iAtom];
      final int igrd0 = bSplines.initGrid[iSymm][iAtom][0];
      final int jgrd0 = bSplines.initGrid[iSymm][iAtom][1];
      k0 = bSplines.initGrid[iSymm][iAtom][2];
      final double c = fm[t000];
      final double dx = fm[t100];
      final double dy = fm[t010];
      final double dz = fm[t001];
      final double qxx = fm[t200];
      final double qyy = fm[t020];
      final double qzz = fm[t002];
      final double qxy = fm[t110];
      final double qxz = fm[t101];
      final double qyz = fm[t011];
      for (int ith3 = 0; ith3 < bSplineOrder; ith3++) {
        final double[] splzi = splz[ith3];
        final double v0 = splzi[0];
        final double v1 = splzi[1];
        final double v2 = splzi[2];
        final double c0 = c * v0;
        final double dx0 = dx * v0;
        final double dy0 = dy * v0;
        final double dz1 = dz * v1;
        final double qxx0 = qxx * v0;
        final double qyy0 = qyy * v0;
        final double qzz2 = qzz * v2;
        final double qxy0 = qxy * v0;
        final double qxz1 = qxz * v1;
        final double qyz1 = qyz * v1;
        final int k = mod(++k0, fftZ);
        if (k < lbZ || k > ubZ) {
File Project Line
ffx/crystal/SpaceGroupDefinitions.java Crystal 6586
ffx/crystal/SpaceGroupDefinitions.java Crystal 6834
ffx/crystal/SpaceGroupDefinitions.java Crystal 7658
new double[] {f12, f12, f12},
                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));
File Project Line
ffx/xray/CrystalReciprocalSpace.java X-Ray Refinement 1990
ffx/xray/CrystalReciprocalSpace.java X-Ray Refinement 2225
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() {
File Project Line
ffx/potential/nonbonded/implicit/GKEnergyRegion.java Potential 650
ffx/potential/nonbonded/implicit/GKEnergyRegionQI.java Potential 590
ffx/potential/nonbonded/implicit/PermanentGKFieldRegion.java Potential 305
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]);
File Project Line
ffx/potential/nonbonded/implicit/GKEnergyRegion.java Potential 361
ffx/potential/nonbonded/implicit/GKEnergyRegionQI.java Potential 299
sharedPolarizationGKEnergy.addAndGet(gkPolarizationEnergy);
    }

    @Override
    public void run(int lb, int ub) {

      double[] x = sXYZ[0][0];
      double[] y = sXYZ[0][1];
      double[] z = sXYZ[0][2];

      int nSymm = crystal.spaceGroup.symOps.size();
      List<SymOp> symOps = crystal.spaceGroup.symOps;

      for (iSymm = 0; iSymm < nSymm; iSymm++) {
        SymOp symOp = symOps.get(iSymm);
        crystal.getTransformationOperator(symOp, transOp);
        for (int i = lb; i <= ub; i++) {
          if (!use[i]) {
            continue;
          }

          // Zero out force accumulation for atom i.
          dedxi = 0.0;
          dedyi = 0.0;
          dedzi = 0.0;
          dborni = 0.0;
          trqxi = 0.0;
          trqyi = 0.0;
          trqzi = 0.0;

          xi = x[i];
          yi = y[i];
          zi = z[i];
          final double[] multipolei = globalMultipole[0][i];
          ci = multipolei[t000];
          uxi = multipolei[t100];
          uyi = multipolei[t010];
          uzi = multipolei[t001];
          qxxi = multipolei[t200] * oneThird;
          qxyi = multipolei[t110] * oneThird;
          qxzi = multipolei[t101] * oneThird;
          qyyi = multipolei[t020] * oneThird;
          qyzi = multipolei[t011] * oneThird;
          qzzi = multipolei[t002] * oneThird;
          dxi = inducedDipole[0][i][0];
          dyi = inducedDipole[0][i][1];
          dzi = inducedDipole[0][i][2];
          pxi = inducedDipoleCR[0][i][0];
          pyi = inducedDipoleCR[0][i][1];
          pzi = inducedDipoleCR[0][i][2];
          sxi = dxi + pxi;
          syi = dyi + pyi;
          szi = dzi + pzi;
          rbi = born[i];
          int[] list = neighborLists[iSymm][i];
          for (int k : list) {
            if (!use[k]) {
              continue;
            }
            interaction(i, k);
          }
          if (iSymm == 0) {
            // Include self-interactions for the asymmetric unit atoms.
            interaction(i, i);
            /*
             Formula for Born energy approximation for cavitation energy is:
             e = surfaceTension / 6 * (ri + probe)^2 * (ri/rb)^6.
             ri is the base atomic radius the atom.
             rb is Born radius of the atom.
            */
            switch (nonPolar) {
              case BORN_SOLV:
              case BORN_CAV_DISP:
                double r = baseRadius[i] + dOffset + probe;
                double ratio = (baseRadius[i] + dOffset) / born[i];
                ratio *= ratio;
                ratio *= (ratio * ratio);
                double saTerm = surfaceTension * r * r * ratio / 6.0;
                gkEnergy += saTerm;
                sharedBornGrad.sub(threadID, i, 6.0 * saTerm / born[i]);
                break;
              default:
                break;
            }
          }
          if (gradient) {
            grad.add(threadID, i, dedxi, dedyi, dedzi);
            torque.add(threadID, i, trqxi, trqyi, trqzi);
            sharedBornGrad.add(threadID, i, dborni);
          }
        }
      }
    }

    public void setGradient(boolean gradient) {
      this.gradient = gradient;
    }

    @Override
    public void start() {
      gkEnergy = 0.0;
File Project Line
ffx/potential/parsers/XPHFilter.java Potential 648
ffx/potential/parsers/XYZFilter.java Potential 641
}

      return true;

    } catch (FileNotFoundException e) {
      String message = format("Exception opening file %s.", currentFile);
      logger.log(Level.WARNING, message, e);
    } catch (IOException e) {
      String message = format("Exception reading from file %s.", currentFile);
      logger.log(Level.WARNING, message, e);
    }
    return false;
  }

  /** {@inheritDoc} */
  @Override
  public boolean writeFile(File saveFile, boolean append, String[] extraLines) {
    if (saveFile == null) {
      return false;
    }

    File newFile = saveFile;
    if (!append) {
      newFile = version(saveFile);
    }
    activeMolecularAssembly.setFile(newFile);
    if (activeMolecularAssembly.getName() == null) {
      activeMolecularAssembly.setName(newFile.getName());
    }

    try (FileWriter fw = new FileWriter(newFile, append && newFile.exists());
        BufferedWriter bw = new BufferedWriter(fw)) {
      // XYZ File First Line
      int numberOfAtoms = activeMolecularAssembly.getAtomList().size();
      StringBuilder sb =
          new StringBuilder(format("%7d  %s", numberOfAtoms, activeMolecularAssembly.getName()));
      if (extraLines != null) {
        for (String line : extraLines) {
          line = line.replaceAll("\n", " ");
          sb.append(" ").append(line);
        }
      }
      String output = sb.append("\n").toString();
      bw.write(output);

      Crystal crystal = activeMolecularAssembly.getCrystal();
      if (crystal!=null && !crystal.aperiodic()) {
        Crystal uc = crystal.getUnitCell();
        String params =
            format("%14.8f%14.8f%14.8f%14.8f%14.8f%14.8f\n",
                uc.a, uc.b, uc.c, uc.alpha, uc.beta, uc.gamma);
        bw.write(params);
      }

      Atom a2;
      StringBuilder line;
      StringBuilder[] lines = new StringBuilder[numberOfAtoms];
      // XYZ File Atom Lines
      List<Atom> atoms = activeMolecularAssembly.getAtomList();
      Vector3d offset = activeMolecularAssembly.getOffset();
      for (Atom a : atoms) {
        if (vdwH) {
          line =
              new StringBuilder(
                  format(
                      "%7d %3s%14.8f%14.8f%14.8f%6d",
                      a.getIndex(),
                      a.getAtomType().name,
                      a.getRedX() - offset.x,
                      a.getRedY() - offset.y,
                      a.getRedZ() - offset.z,
                      a.getType()));
        } else {
          line =
              new StringBuilder(
                  format(
                      "%7d %3s%14.8f%14.8f%14.8f%6d",
                      a.getIndex(),
                      a.getAtomType().name,
                      a.getX() - offset.x,
                      a.getY() - offset.y,
                      a.getZ() - offset.z,
                      a.getType()));
        }
        for (Bond b : a.getBonds()) {
          a2 = b.get1_2(a);
          line.append(format("%8d", a2.getIndex()));
        }
        lines[a.getIndex() - 1] = line.append("\n");
      }
File Project Line
ffx/potential/utils/Loop.java Potential 1244
ffx/potential/utils/Loop.java Potential 1366
arraycopy(determineIntxyz(ca, dCB_CA, n, 109.5, c, 107.8, 1), 0, cb, 0, 3);
            coordsArray = fillCoordsArray(CB, coordsArray, cb);

            arraycopy(
                determineIntxyz(cb, dCG_CB, ca, dCG_CB_CA, n, rotScale * 180.0, 0),
                0,
                cg,
                0,
                3); // CG
            coordsArray = fillCoordsArray(CG, coordsArray, cg);

            arraycopy(determineIntxyz(cg, dCD_CG, cb, dCD_CG_CB, ca, 90.0, 0), 0, cd1, 0, 3); // CD1
            coordsArray = fillCoordsArray(CD1, coordsArray, cd1);

            arraycopy(
                determineIntxyz(cg, dCD_CG, cb, dCD_CG_CB, cd1, 120.0, 1), 0, cd2, 0, 3); // CD2
            coordsArray = fillCoordsArray(CD2, coordsArray, cd2);

            arraycopy(determineIntxyz(cd1, dCE_CD, cg, dCE_CD_CG, cb, 180, 0), 0, ce1, 0, 3); // CE1
            coordsArray = fillCoordsArray(CE1, coordsArray, ce1);

            arraycopy(determineIntxyz(cd2, dCE_CD, cg, dCE_CD_CG, cb, 180, 0), 0, ce2, 0, 3); // CE2
            coordsArray = fillCoordsArray(CE2, coordsArray, ce2);

            arraycopy(
                determineIntxyz(ce1, dCZ_CE1, cd1, dCZ_CE1_CD1, cg, 0.0, 0), 0, cz, 0, 3); // CZ
            coordsArray = fillCoordsArray(CZ, coordsArray, cz);

            arraycopy(
                determineIntxyz(cz, dOH_CZ, ce2, dOH_CZ_CE2, ce1, 120.0, 1), 0, oh, 0, 3); // OH
            coordsArray = fillCoordsArray(OH, coordsArray, oh);

            arraycopy(
                determineIntxyz(cb, dHB_CB, ca, dHB_CB_CA, cg, 109.4, 1), 0, hb2, 0, 3); // HB2
            coordsArray = fillCoordsArray(HB2, coordsArray, hb2);

            arraycopy(
                determineIntxyz(cb, dHB_CB, ca, dHB_CB_CA, cg, 109.4, -1), 0, hb3, 0, 3); // HB3
            coordsArray = fillCoordsArray(HB3, coordsArray, hb3);

            arraycopy(
                determineIntxyz(cd1, dHD_CD, cg, dHD_CD_CG, ce1, 120.0, 1), 0, hd1, 0, 3); // HD1
            coordsArray = fillCoordsArray(HD1, coordsArray, hd1);

            arraycopy(
                determineIntxyz(cd2, dHD_CD, cg, dHD_CD_CG, ce2, 120.0, 1), 0, hd2, 0, 3); // HD2
            coordsArray = fillCoordsArray(HD2, coordsArray, hd2);

            arraycopy(
                determineIntxyz(ce1, dHE_CE, cd1, dHE_CE_CD, cz, 120.0, 1), 0, he1, 0, 3); // HE1
            coordsArray = fillCoordsArray(HE1, coordsArray, he1);

            arraycopy(
                determineIntxyz(ce2, dHE_CE, cd2, dHE_CE_CD, cz, 120.0, 1), 0, he2, 0, 3); // HE2
            coordsArray = fillCoordsArray(HE2, coordsArray, he2);
File Project Line
edu/rit/pj/WorkerIntegerForLoop.java Parallel Java 297
edu/rit/pj/WorkerIntegerStrideForLoop.java Parallel Java 300
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,