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, |