CPD Results
The following document contains the results of PMD's CPD 7.0.0.
Duplications
| File | Project | Line |
|---|---|---|
| ffx/numerics/multipole/GKTensorGlobalSIMD.java | Numerics | 74 |
| ffx/numerics/multipole/GKTensorQISIMD.java | Numerics | 74 |
public GKTensorGlobalSIMD(GKMultipoleOrder multipoleOrder, int order, GKSourceSIMD gkSource, double Eh, double Es) {
super(order);
this.multipoleOrder = multipoleOrder;
this.gkSource = gkSource;
// Load the dielectric function
c = GKSource.cn(multipoleOrder.getOrder(), Eh, Es);
}
/**
* GK Permanent multipole energy.
*
* @param mI PolarizableMultipole at site I.
* @param mK PolarizableMultipole at site K.
* @return a double.
*/
@Override
public DoubleVector multipoleEnergy(PolarizableMultipoleSIMD mI, PolarizableMultipoleSIMD mK) {
return switch (multipoleOrder) {
default -> {
chargeIPotentialAtK(mI, 2);
DoubleVector eK = multipoleEnergy(mK);
chargeKPotentialAtI(mK, 2);
DoubleVector eI = multipoleEnergy(mI);
yield eK.add(eI).mul(c * 0.5);
}
case DIPOLE -> {
dipoleIPotentialAtK(mI.dx, mI.dy, mI.dz, 2);
DoubleVector eK = multipoleEnergy(mK);
dipoleKPotentialAtI(mK.dx, mK.dy, mK.dz, 2);
DoubleVector eI = multipoleEnergy(mI);
yield eK.add(eI).mul(c * 0.5);
}
case QUADRUPOLE -> {
quadrupoleIPotentialAtK(mI, 2);
DoubleVector eK = multipoleEnergy(mK);
quadrupoleKPotentialAtI(mK, 2);
DoubleVector eI = multipoleEnergy(mI);
yield eK.add(eI).mul(c * 0.5);
}
};
}
/**
* 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 DoubleVector multipoleEnergyAndGradient(PolarizableMultipoleSIMD mI, PolarizableMultipoleSIMD mK,
DoubleVector[] Gi, DoubleVector[] Gk,
DoubleVector[] Ti, DoubleVector[] Tk) {
return switch (multipoleOrder) {
default -> monopoleEnergyAndGradient(mI, mK, Gi, Gk, Ti, Tk);
case DIPOLE -> dipoleEnergyAndGradient(mI, mK, Gi, Gk, Ti, Tk);
case QUADRUPOLE -> quadrupoleEnergyAndGradient(mI, mK, Gi, Gk, Ti, Tk);
};
}
/**
* Permanent multipole energy and gradient using the GK monopole tensor.
*
* @param mI PolarizableMultipole at site I.
* @param mK PolarizableMultipole at site K.
* @param Gi Coordinate gradient at site I.
* @param Gk Coordinate gradient at site K.
* @param Ti Torque at site I.
* @param Tk Torque at site K.
* @return the permanent multipole GK energy.
*/
protected DoubleVector monopoleEnergyAndGradient(PolarizableMultipoleSIMD mI, PolarizableMultipoleSIMD mK,
DoubleVector[] Gi, DoubleVector[] Gk,
DoubleVector[] Ti, DoubleVector[] Tk) {
// Compute the potential due to a multipole component at site I.
chargeIPotentialAtK(mI, 3);
DoubleVector eK = multipoleEnergy(mK);
multipoleGradient(mK, Gk);
multipoleTorque(mK, Tk);
// Compute the potential due to a multipole component at site K.
chargeKPotentialAtI(mK, 3);
DoubleVector eI = multipoleEnergy(mI);
multipoleGradient(mI, Gi);
multipoleTorque(mI, Ti);
double scale = c * 0.5;
Gi[0] = Gi[0].sub(Gk[0]).mul(scale);
Gi[1] = Gi[1].sub(Gk[1]).mul(scale);
Gi[2] = Gi[2].sub(Gk[2]).mul(scale);
Gk[0] = Gi[0].neg();
Gk[1] = Gi[1].neg();
Gk[2] = Gi[2].neg();
Ti[0] = Ti[0].mul(scale);
Ti[1] = Ti[1].mul(scale);
Ti[2] = Ti[2].mul(scale);
Tk[0] = Tk[0].mul(scale);
Tk[1] = Tk[1].mul(scale);
Tk[2] = Tk[2].mul(scale);
return eK.add(eI).mul(scale);
}
/**
* 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 DoubleVector dipoleEnergyAndGradient(PolarizableMultipoleSIMD mI, PolarizableMultipoleSIMD mK,
DoubleVector[] Gi, DoubleVector[] Gk,
DoubleVector[] Ti, DoubleVector[] Tk) {
// Compute the potential due to a multipole component at site I.
dipoleIPotentialAtK(mI.dx, mI.dy, mI.dz, 3);
DoubleVector 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);
DoubleVector eI = multipoleEnergy(mI);
multipoleGradient(mI, Gi);
multipoleTorque(mI, Ti);
// Need the torque on site K pole due to multipole on site I.
// Only torque on the site K dipole.
multipoleIPotentialAtK(mI, 1);
dipoleTorque(mK, Tk);
double scale = c * 0.5;
Gi[0] = Gi[0].sub(Gk[0]).mul(scale);
Gi[1] = Gi[1].sub(Gk[1]).mul(scale);
Gi[2] = Gi[2].sub(Gk[2]).mul(scale);
Gk[0] = Gi[0].neg();
Gk[1] = Gi[1].neg();
Gk[2] = Gi[2].neg();
Ti[0] = Ti[0].mul(scale);
Ti[1] = Ti[1].mul(scale);
Ti[2] = Ti[2].mul(scale);
Tk[0] = Tk[0].mul(scale);
Tk[1] = Tk[1].mul(scale);
Tk[2] = Tk[2].mul(scale);
return eK.add(eI).mul(scale);
}
/**
* 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 DoubleVector quadrupoleEnergyAndGradient(PolarizableMultipoleSIMD mI, PolarizableMultipoleSIMD mK,
DoubleVector[] Gi, DoubleVector[] Gk,
DoubleVector[] Ti, DoubleVector[] Tk) {
// Compute the potential due to a multipole component at site I.
quadrupoleIPotentialAtK(mI, 3);
DoubleVector eK = multipoleEnergy(mK);
multipoleGradient(mK, Gk);
multipoleTorque(mK, Tk);
// Need the torque on site I quadrupole due to site K multipole.
multipoleKPotentialAtI(mK, 2);
quadrupoleTorque(mI, Ti);
// Compute the potential due to a multipole component at site K.
quadrupoleKPotentialAtI(mK, 3);
DoubleVector eI = multipoleEnergy(mI);
multipoleGradient(mI, Gi);
multipoleTorque(mI, Ti);
// Need the torque on site K quadrupole due to site I multipole.
multipoleIPotentialAtK(mI, 2);
quadrupoleTorque(mK, Tk);
double scale = c * 0.5;
Gi[0] = Gi[0].sub(Gk[0]).mul(scale);
Gi[1] = Gi[1].sub(Gk[1]).mul(scale);
Gi[2] = Gi[2].sub(Gk[2]).mul(scale);
Gk[0] = Gi[0].neg();
Gk[1] = Gi[1].neg();
Gk[2] = Gi[2].neg();
Ti[0] = Ti[0].mul(scale);
Ti[1] = Ti[1].mul(scale);
Ti[2] = Ti[2].mul(scale);
Tk[0] = Tk[0].mul(scale);
Tk[1] = Tk[1].mul(scale);
Tk[2] = Tk[2].mul(scale);
return eK.add(eI).mul(scale);
}
/**
* GK Permanent multipole Born grad.
*
* @param mI PolarizableMultipole at site I.
* @param mK PolarizableMultipole at site K.
* @return a double.
*/
public DoubleVector multipoleEnergyBornGrad(PolarizableMultipoleSIMD mI, PolarizableMultipoleSIMD mK) {
return multipoleEnergy(mI, mK);
}
/**
* GK Polarization Energy.
*
* @param mI PolarizableMultipole at site I.
* @param mK PolarizableMultipole at site K.
* @param scaleEnergy This is ignored, since masking/scaling is not applied to GK interactions
* (everything is intermolecular).
* @return a double.
*/
public DoubleVector polarizationEnergy(PolarizableMultipoleSIMD mI, PolarizableMultipoleSIMD mK, DoubleVector scaleEnergy) {
return polarizationEnergy(mI, mK);
}
/**
* GK Polarization Energy.
*
* @param mI PolarizableMultipole at site I.
* @param mK PolarizableMultipole at site K.
* @return a double.
*/
public DoubleVector polarizationEnergy(PolarizableMultipoleSIMD mI, PolarizableMultipoleSIMD mK) {
return switch (multipoleOrder) {
default -> {
// Find the GK charge potential of site I at site K.
chargeIPotentialAtK(mI, 1);
// Energy of induced dipole K in the field of permanent charge I.
DoubleVector 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.
DoubleVector eI = polarizationEnergy(mI);
yield eK.add(eI).mul(c * 0.5);
}
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.
DoubleVector 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 = multipoleEnergy(mK).mul(0.5).add(eK);
// 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.
DoubleVector 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 = multipoleEnergy(mI).mul(0.5).add(eI);
yield eK.add(eI).mul(c * 0.5);
}
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.
DoubleVector 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.
DoubleVector eI = polarizationEnergy(mI);
yield eK.add(eI).mul(c * 0.5);
}
};
}
/**
* GK Polarization Energy.
*
* @param mI PolarizableMultipole at site I.
* @param mK PolarizableMultipole at site K.
* @return a double.
*/
public DoubleVector polarizationEnergyBorn(PolarizableMultipoleSIMD mI, PolarizableMultipoleSIMD mK) {
return switch (multipoleOrder) {
default -> {
// Find the GK charge potential of site I at site K.
chargeIPotentialAtK(mI, 1);
// Energy of induced dipole K in the field of permanent charge I.
DoubleVector eK = polarizationEnergyS(mK);
// Find the GK charge potential of site K at site I.
chargeKPotentialAtI(mK, 1);
// Energy of induced dipole I in the field of permanent charge K.
DoubleVector eI = polarizationEnergyS(mI);
yield eK.add(eI).mul(c * 0.5);
}
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.
DoubleVector eK = polarizationEnergyS(mK);
// Find the GK induced dipole potential of site I at site K.
dipoleIPotentialAtK(mI.sx, mI.sy, mI.sz, 2);
// Energy of permanent multipole K in the field of induced dipole I.
eK = multipoleEnergy(mK).mul(0.5).add(eK);
// 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.
DoubleVector eI = polarizationEnergyS(mI);
// Find the GK induced dipole potential of site K at site I.
dipoleKPotentialAtI(mK.sx, mK.sy, mK.sz, 2);
// Energy of permanent multipole I in the field of induced dipole K.
eI = multipoleEnergy(mI).mul(0.5).add(eI);
yield eK.add(eI).mul(c * 0.5);
}
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.
DoubleVector eK = polarizationEnergyS(mK);
// Find the GK quadrupole potential of site K at site I.
quadrupoleKPotentialAtI(mK, 1);
// Energy of induced dipole I in the field of permanent quadrupole K.
DoubleVector eI = polarizationEnergyS(mI);
yield eK.add(eI).mul(c * 0.5);
}
};
}
/**
* 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 DoubleVector polarizationEnergyAndGradient(PolarizableMultipoleSIMD mI, PolarizableMultipoleSIMD mK,
DoubleVector inductionMask, DoubleVector energyMask, DoubleVector mutualMask,
DoubleVector[] Gi, DoubleVector[] Ti, DoubleVector[] Tk) {
return switch (multipoleOrder) {
default -> monopolePolarizationEnergyAndGradient(mI, mK, Gi);
case DIPOLE -> dipolePolarizationEnergyAndGradient(mI, mK, mutualMask, Gi, Ti, Tk);
case QUADRUPOLE -> quadrupolePolarizationEnergyAndGradient(mI, mK, Gi, Ti, Tk);
};
}
/**
* Monopole Polarization Energy and Gradient.
*
* @param mI PolarizableMultipole at site I.
* @param mK PolarizableMultipole at site K.
* @param Gi an array of {@link double} objects.
* @return a double.
*/
public DoubleVector monopolePolarizationEnergyAndGradient(PolarizableMultipoleSIMD mI,
PolarizableMultipoleSIMD mK, DoubleVector[] Gi) {
// Find the permanent multipole potential at site k.
chargeIPotentialAtK(mI, 2);
// Energy of induced dipole k in the field of multipole i.
DoubleVector eK = polarizationEnergy(mK);
// Derivative with respect to moving atom k.
Gi[0] = mK.sx.mul(E200).add(mK.sy.mul(E110)).add(mK.sz.mul(E101)).neg();
Gi[1] = mK.sx.mul(E110).add(mK.sy.mul(E020)).add(mK.sz.mul(E011)).neg();
Gi[2] = mK.sx.mul(E101).add(mK.sy.mul(E011)).add(mK.sz.mul(E002)).neg();
// Find the permanent multipole potential and derivatives at site i.
chargeKPotentialAtI(mK, 2);
// Energy of induced dipole i in the field of multipole k.
DoubleVector eI = polarizationEnergy(mI);
// Derivative with respect to moving atom i.
Gi[0] = Gi[0].add(mI.sx.mul(E200)).add(mI.sy.mul(E110)).add(mI.sz.mul(E101));
Gi[1] = Gi[1].add(mI.sx.mul(E110)).add(mI.sy.mul(E020)).add(mI.sz.mul(E011));
Gi[2] = Gi[2].add(mI.sx.mul(E101)).add(mI.sy.mul(E011)).add(mI.sz.mul(E002));
double scale = c * 0.5;
Gi[0] = Gi[0].mul(scale);
Gi[1] = Gi[1].mul(scale);
Gi[2] = Gi[2].mul(scale);
// Total polarization energy.
return eI.add(eK).mul(scale);
}
/**
* 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 DoubleVector dipolePolarizationEnergyAndGradient(PolarizableMultipoleSIMD mI, PolarizableMultipoleSIMD mK,
DoubleVector mutualMask, DoubleVector[] Gi,
DoubleVector[] Ti, DoubleVector[] 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.
DoubleVector eK = polarizationEnergy(mK);
// Derivative with respect to moving atom k.
Gi[0] = mK.sx.mul(E200).add(mK.sy.mul(E110)).add(mK.sz.mul(E101)).neg();
Gi[1] = mK.sx.mul(E110).add(mK.sy.mul(E020)).add(mK.sz.mul(E011)).neg();
Gi[2] = mK.sx.mul(E101).add(mK.sy.mul(E011)).add(mK.sz.mul(E002)).neg();
// Find the potential at K due to the averaged induced dipole at site 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 = eK.add(multipoleEnergy(mK).mul(0.5));
DoubleVector[] G = new DoubleVector[3];
multipoleGradient(mK, G);
Gi[0] = Gi[0].sub(G[0]);
Gi[1] = Gi[1].sub(G[1]);
Gi[2] = Gi[2].sub(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.
DoubleVector eI = polarizationEnergy(mI);
// Derivative with respect to moving atom i.
Gi[0] = Gi[0].add(mI.sx.mul(E200)).add(mI.sy.mul(E110)).add(mI.sz.mul(E101));
Gi[1] = Gi[1].add(mI.sx.mul(E110)).add(mI.sy.mul(E020)).add(mI.sz.mul(E011));
Gi[2] = Gi[2].add(mI.sx.mul(E101)).add(mI.sy.mul(E011)).add(mI.sz.mul(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 = eI.add(multipoleEnergy(mI).mul(0.5));
multipoleGradient(mI, G);
Gi[0] = Gi[0].add(G[0]);
Gi[1] = Gi[1].add(G[1]);
Gi[2] = Gi[2].add(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.eq(zero).allTrue()) {
// Find the potential and its derivatives at k due to induced dipole i.
dipoleIPotentialAtK(mI.ux, mI.uy, mI.uz, 2);
Gi[0] = Gi[0].sub((mK.px.mul(E200).add(mK.py.mul(E110)).add(mK.pz.mul(E101)).mul(mutualMask)));
Gi[1] = Gi[1].sub((mK.px.mul(E110).add(mK.py.mul(E020)).add(mK.pz.mul(E011)).mul(mutualMask)));
Gi[2] = Gi[2].sub((mK.px.mul(E101).add(mK.py.mul(E011)).add(mK.pz.mul(E002)).mul(mutualMask)));
// Find the potential and its derivatives at i due to induced dipole k.
dipoleKPotentialAtI(mK.ux, mK.uy, mK.uz, 2);
Gi[0] = Gi[0].sub((mI.px.mul(E200).add(mI.py.mul(E110)).add(mI.pz.mul(E101)).mul(mutualMask)));
Gi[1] = Gi[1].sub((mI.px.mul(E110).add(mI.py.mul(E020)).add(mI.pz.mul(E011)).mul(mutualMask)));
Gi[2] = Gi[2].sub((mI.px.mul(E101).add(mI.py.mul(E011)).add(mI.pz.mul(E002)).mul(mutualMask)));
}
// Total polarization energy.
double scale = c * 0.5;
DoubleVector energy = eI.add(eK).mul(scale);
Gi[0] = Gi[0].mul(scale);
Gi[1] = Gi[1].mul(scale);
Gi[2] = Gi[2].mul(scale);
Ti[0] = Ti[0].mul(scale);
Ti[1] = Ti[1].mul(scale);
Ti[2] = Ti[2].mul(scale);
Tk[0] = Tk[0].mul(scale);
Tk[1] = Tk[1].mul(scale);
Tk[2] = Tk[2].mul(scale);
return energy;
}
/**
* Quadrupole Polarization Energy and Gradient.
*
* @param mI PolarizableMultipole at site I.
* @param mK PolarizableMultipole at site K.
* @param Gi an array of {@link double} objects.
* @param Ti an array of {@link double} objects.
* @param Tk an array of {@link double} objects.
* @return a double.
*/
public DoubleVector quadrupolePolarizationEnergyAndGradient(PolarizableMultipoleSIMD mI, PolarizableMultipoleSIMD mK,
DoubleVector[] Gi, DoubleVector[] Ti, DoubleVector[] 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.
DoubleVector eK = polarizationEnergy(mK);
// Derivative with respect to moving atom k.
Gi[0] = (mK.sx.mul(E200).add(mK.sy.mul(E110)).add(mK.sz.mul(E101))).neg();
Gi[1] = (mK.sx.mul(E110).add(mK.sy.mul(E020)).add(mK.sz.mul(E011))).neg();
Gi[2] = (mK.sx.mul(E101).add(mK.sy.mul(E011)).add(mK.sz.mul(E002))).neg();
// Find the permanent multipole potential and derivatives at site i.
quadrupoleKPotentialAtI(mK, 2);
// Energy of induced dipole i in the field of multipole k.
DoubleVector eI = polarizationEnergy(mI);
// Derivative with respect to moving atom i.
Gi[0] = Gi[0].add(mI.sx.mul(E200)).add(mI.sy.mul(E110)).add(mI.sz.mul(E101));
Gi[1] = Gi[1].add(mI.sx.mul(E110)).add(mI.sy.mul(E020)).add(mI.sz.mul(E011));
Gi[2] = Gi[2].add(mI.sx.mul(E101)).add(mI.sy.mul(E011)).add(mI.sz.mul(E002));
double scale = c * 0.5;
Gi[0] = Gi[0].mul(scale);
Gi[1] = Gi[1].mul(scale);
Gi[2] = Gi[2].mul(scale);
// Find the potential and its derivatives at K due to the averaged induced dipole at site i.
dipoleIPotentialAtK(mI.sx.mul(scale), mI.sy.mul(scale), mI.sz.mul(scale), 2);
quadrupoleTorque(mK, Tk);
// Find the potential and its derivatives at I due to the averaged induced dipole at k.
dipoleKPotentialAtI(mK.sx.mul(scale), mK.sy.mul(scale), mK.sz.mul(scale), 2);
quadrupoleTorque(mI, Ti);
// Total polarization energy.
return eI.add(eK).mul(scale);
}
/**
* 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 DoubleVector polarizationEnergyBornGrad(PolarizableMultipoleSIMD mI, PolarizableMultipoleSIMD mK) {
return polarizationEnergyBorn(mI, mK).mul(2.0);
}
/**
* 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 DoubleVector mutualPolarizationEnergyBornGrad(PolarizableMultipoleSIMD mI, PolarizableMultipoleSIMD mK) {
DoubleVector db = zero;
if (multipoleOrder == GKMultipoleOrder.DIPOLE) {
// Find the potential and its derivatives at k due to induced dipole i.
dipoleIPotentialAtK(mI.ux, mI.uy, mI.uz, 2);
db = mK.px.mul(E100).add(mK.py.mul(E010)).add(mK.pz.mul(E001)).mul(0.5);
// Find the potential and its derivatives at i due to induced dipole k.
dipoleKPotentialAtI(mK.ux, mK.uy, mK.uz, 2);
db = db.add(mI.px.mul(E100).add(mI.py.mul(E010)).add(mI.pz.mul(E001)).mul(0.5));
}
return db.mul(c);
}
/**
* Generate source terms for the Kirkwood version of the Challacombe et al. recursion.
*/
@Override
protected void source(DoubleVector[] work) {
gkSource.source(work, multipoleOrder);
}
} | ||
| File | Project | Line |
|---|---|---|
| ffx/numerics/multipole/MultipoleTensor.java | Numerics | 298 |
| ffx/numerics/multipole/MultipoleTensorSIMD.java | Numerics | 287 |
T000 = new double[order + 1];
// l + m + n = 0 (1)
t000 = MultipoleUtilities.ti(0, 0, 0, order);
// l + m + n = 1 (3) 4
t100 = MultipoleUtilities.ti(1, 0, 0, order);
t010 = MultipoleUtilities.ti(0, 1, 0, order);
t001 = MultipoleUtilities.ti(0, 0, 1, order);
// l + m + n = 2 (6) 10
t200 = MultipoleUtilities.ti(2, 0, 0, order);
t020 = MultipoleUtilities.ti(0, 2, 0, order);
t002 = MultipoleUtilities.ti(0, 0, 2, order);
t110 = MultipoleUtilities.ti(1, 1, 0, order);
t101 = MultipoleUtilities.ti(1, 0, 1, order);
t011 = MultipoleUtilities.ti(0, 1, 1, order);
// l + m + n = 3 (10) 20
t300 = MultipoleUtilities.ti(3, 0, 0, order);
t030 = MultipoleUtilities.ti(0, 3, 0, order);
t003 = MultipoleUtilities.ti(0, 0, 3, order);
t210 = MultipoleUtilities.ti(2, 1, 0, order);
t201 = MultipoleUtilities.ti(2, 0, 1, order);
t120 = MultipoleUtilities.ti(1, 2, 0, order);
t021 = MultipoleUtilities.ti(0, 2, 1, order);
t102 = MultipoleUtilities.ti(1, 0, 2, order);
t012 = MultipoleUtilities.ti(0, 1, 2, order);
t111 = MultipoleUtilities.ti(1, 1, 1, order);
// l + m + n = 4 (15) 35
t400 = MultipoleUtilities.ti(4, 0, 0, order);
t040 = MultipoleUtilities.ti(0, 4, 0, order);
t004 = MultipoleUtilities.ti(0, 0, 4, order);
t310 = MultipoleUtilities.ti(3, 1, 0, order);
t301 = MultipoleUtilities.ti(3, 0, 1, order);
t130 = MultipoleUtilities.ti(1, 3, 0, order);
t031 = MultipoleUtilities.ti(0, 3, 1, order);
t103 = MultipoleUtilities.ti(1, 0, 3, order);
t013 = MultipoleUtilities.ti(0, 1, 3, order);
t220 = MultipoleUtilities.ti(2, 2, 0, order);
t202 = MultipoleUtilities.ti(2, 0, 2, order);
t022 = MultipoleUtilities.ti(0, 2, 2, order);
t211 = MultipoleUtilities.ti(2, 1, 1, order);
t121 = MultipoleUtilities.ti(1, 2, 1, order);
t112 = MultipoleUtilities.ti(1, 1, 2, order);
// l + m + n = 5 (21) 56
t500 = MultipoleUtilities.ti(5, 0, 0, order);
t050 = MultipoleUtilities.ti(0, 5, 0, order);
t005 = MultipoleUtilities.ti(0, 0, 5, order);
t410 = MultipoleUtilities.ti(4, 1, 0, order);
t401 = MultipoleUtilities.ti(4, 0, 1, order);
t140 = MultipoleUtilities.ti(1, 4, 0, order);
t041 = MultipoleUtilities.ti(0, 4, 1, order);
t104 = MultipoleUtilities.ti(1, 0, 4, order);
t014 = MultipoleUtilities.ti(0, 1, 4, order);
t320 = MultipoleUtilities.ti(3, 2, 0, order);
t302 = MultipoleUtilities.ti(3, 0, 2, order);
t230 = MultipoleUtilities.ti(2, 3, 0, order);
t032 = MultipoleUtilities.ti(0, 3, 2, order);
t203 = MultipoleUtilities.ti(2, 0, 3, order);
t023 = MultipoleUtilities.ti(0, 2, 3, order);
t311 = MultipoleUtilities.ti(3, 1, 1, order);
t131 = MultipoleUtilities.ti(1, 3, 1, order);
t113 = MultipoleUtilities.ti(1, 1, 3, order);
t221 = MultipoleUtilities.ti(2, 2, 1, order);
t212 = MultipoleUtilities.ti(2, 1, 2, order);
t122 = MultipoleUtilities.ti(1, 2, 2, order);
// l + m + n = 6 (28) 84
t600 = MultipoleUtilities.ti(6, 0, 0, order);
t060 = MultipoleUtilities.ti(0, 6, 0, order);
t006 = MultipoleUtilities.ti(0, 0, 6, order);
t510 = MultipoleUtilities.ti(5, 1, 0, order);
t501 = MultipoleUtilities.ti(5, 0, 1, order);
t150 = MultipoleUtilities.ti(1, 5, 0, order);
t051 = MultipoleUtilities.ti(0, 5, 1, order);
t105 = MultipoleUtilities.ti(1, 0, 5, order);
t015 = MultipoleUtilities.ti(0, 1, 5, order);
t420 = MultipoleUtilities.ti(4, 2, 0, order);
t402 = MultipoleUtilities.ti(4, 0, 2, order);
t240 = MultipoleUtilities.ti(2, 4, 0, order);
t042 = MultipoleUtilities.ti(0, 4, 2, order);
t204 = MultipoleUtilities.ti(2, 0, 4, order);
t024 = MultipoleUtilities.ti(0, 2, 4, order);
t411 = MultipoleUtilities.ti(4, 1, 1, order);
t141 = MultipoleUtilities.ti(1, 4, 1, order);
t114 = MultipoleUtilities.ti(1, 1, 4, order);
t330 = MultipoleUtilities.ti(3, 3, 0, order);
t303 = MultipoleUtilities.ti(3, 0, 3, order);
t033 = MultipoleUtilities.ti(0, 3, 3, order);
t321 = MultipoleUtilities.ti(3, 2, 1, order);
t231 = MultipoleUtilities.ti(2, 3, 1, order);
t213 = MultipoleUtilities.ti(2, 1, 3, order);
t312 = MultipoleUtilities.ti(3, 1, 2, order);
t132 = MultipoleUtilities.ti(1, 3, 2, order);
t123 = MultipoleUtilities.ti(1, 2, 3, order);
t222 = MultipoleUtilities.ti(2, 2, 2, order);
}
/**
* getd2EdZ2.
*
* @return a double.
*/
public double getd2EdZ2() { | ||
| File | Project | Line |
|---|---|---|
| ffx/numerics/multipole/GKTensorGlobal.java | Numerics | 72 |
| ffx/numerics/multipole/GKTensorQI.java | Numerics | 70 |
public GKTensorGlobal(GKMultipoleOrder multipoleOrder, int order, GKSource gkSource,
double Eh, double Es) {
super(order);
this.multipoleOrder = multipoleOrder;
this.gkSource = gkSource;
// Load the dielectric function
c = GKSource.cn(multipoleOrder.getOrder(), Eh, Es);
}
/**
* GK Permanent multipole energy.
*
* @param mI PolarizableMultipole at site I.
* @param mK PolarizableMultipole at site K.
* @return a double.
*/
@Override
public double multipoleEnergy(PolarizableMultipole mI, PolarizableMultipole mK) {
return switch (multipoleOrder) {
default -> {
chargeIPotentialAtK(mI, 2);
double eK = multipoleEnergy(mK);
chargeKPotentialAtI(mK, 2);
double eI = multipoleEnergy(mI);
yield c * 0.5 * (eK + eI);
}
case DIPOLE -> {
dipoleIPotentialAtK(mI.dx, mI.dy, mI.dz, 2);
double eK = multipoleEnergy(mK);
dipoleKPotentialAtI(mK.dx, mK.dy, mK.dz, 2);
double eI = multipoleEnergy(mI);
yield c * 0.5 * (eK + eI);
}
case QUADRUPOLE -> {
quadrupoleIPotentialAtK(mI, 2);
double eK = multipoleEnergy(mK);
quadrupoleKPotentialAtI(mK, 2);
double eI = multipoleEnergy(mI);
yield c * 0.5 * (eK + eI);
}
};
}
/**
* GK Permanent multipole energy and gradient.
*
* @param mI PolarizableMultipole at site I.
* @param mK PolarizableMultipole at site K.
* @param Gi Coordinate gradient at site I.
* @param Gk Coordinate gradient at site K.
* @param Ti Torque at site I.
* @param Tk Torque at site K.
* @return the permanent multipole GK energy.
*/
@Override
public double multipoleEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
double[] Gi, double[] Gk, double[] Ti, double[] Tk) {
return switch (multipoleOrder) {
default -> monopoleEnergyAndGradient(mI, mK, Gi, Gk, Ti, Tk);
case DIPOLE -> dipoleEnergyAndGradient(mI, mK, Gi, Gk, Ti, Tk);
case QUADRUPOLE -> quadrupoleEnergyAndGradient(mI, mK, Gi, Gk, Ti, Tk);
};
}
/**
* Permanent multipole energy and gradient using the GK monopole tensor.
*
* @param mI PolarizableMultipole at site I.
* @param mK PolarizableMultipole at site K.
* @param Gi Coordinate gradient at site I.
* @param Gk Coordinate gradient at site K.
* @param Ti Torque at site I.
* @param Tk Torque at site K.
* @return the permanent multipole GK energy.
*/
protected double monopoleEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
double[] Gi, double[] Gk, double[] Ti, double[] Tk) {
// Compute the potential due to a multipole component at site I.
chargeIPotentialAtK(mI, 3);
double eK = multipoleEnergy(mK);
multipoleGradient(mK, Gk);
multipoleTorque(mK, Tk);
// Compute the potential due to a multipole component at site K.
chargeKPotentialAtI(mK, 3);
double eI = multipoleEnergy(mI);
multipoleGradient(mI, Gi);
multipoleTorque(mI, Ti);
double scale = c * 0.5;
Gi[0] = scale * (Gi[0] - Gk[0]);
Gi[1] = scale * (Gi[1] - Gk[1]);
Gi[2] = scale * (Gi[2] - Gk[2]);
Gk[0] = -Gi[0];
Gk[1] = -Gi[1];
Gk[2] = -Gi[2];
Ti[0] = scale * Ti[0];
Ti[1] = scale * Ti[1];
Ti[2] = scale * Ti[2];
Tk[0] = scale * Tk[0];
Tk[1] = scale * Tk[1];
Tk[2] = scale * Tk[2];
return scale * (eK + eI);
}
/**
* Permanent multipole energy and gradient using the GK dipole tensor.
*
* @param mI PolarizableMultipole at site I.
* @param mK PolarizableMultipole at site K.
* @param Gi Coordinate gradient at site I.
* @param Gk Coordinate gradient at site K.
* @param Ti Torque at site I.
* @param Tk Torque at site K.
* @return the permanent multipole GK energy.
*/
protected double dipoleEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
double[] Gi, double[] Gk, double[] Ti, double[] Tk) {
// Compute the potential due to a multipole component at site I.
dipoleIPotentialAtK(mI.dx, mI.dy, mI.dz, 3);
double eK = multipoleEnergy(mK);
multipoleGradient(mK, Gk);
multipoleTorque(mK, Tk);
// Need the torque on site I pole due to site K multipole.
// Only torque on the site I dipole.
multipoleKPotentialAtI(mK, 1);
dipoleTorque(mI, Ti);
// Compute the potential due to a multipole component at site K.
dipoleKPotentialAtI(mK.dx, mK.dy, mK.dz, 3);
double eI = multipoleEnergy(mI);
multipoleGradient(mI, Gi);
multipoleTorque(mI, Ti);
// Need the torque on site K pole due to multipole on site I.
// Only torque on the site K dipole.
multipoleIPotentialAtK(mI, 1);
dipoleTorque(mK, Tk);
double scale = c * 0.5;
Gi[0] = scale * (Gi[0] - Gk[0]);
Gi[1] = scale * (Gi[1] - Gk[1]);
Gi[2] = scale * (Gi[2] - Gk[2]);
Gk[0] = -Gi[0];
Gk[1] = -Gi[1];
Gk[2] = -Gi[2];
Ti[0] = scale * Ti[0];
Ti[1] = scale * Ti[1];
Ti[2] = scale * Ti[2];
Tk[0] = scale * Tk[0];
Tk[1] = scale * Tk[1];
Tk[2] = scale * Tk[2];
return scale * (eK + eI);
}
/**
* Permanent multipole energy and gradient using the GK quadrupole tensor.
*
* @param mI PolarizableMultipole at site I.
* @param mK PolarizableMultipole at site K.
* @param Gi Coordinate gradient at site I.
* @param Gk Coordinate gradient at site K.
* @param Ti Torque at site I.
* @param Tk Torque at site K.
* @return the permanent multipole GK energy.
*/
protected double quadrupoleEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
double[] Gi, double[] Gk, double[] Ti, double[] Tk) {
// Compute the potential due to a multipole component at site I.
quadrupoleIPotentialAtK(mI, 3);
double eK = multipoleEnergy(mK);
multipoleGradient(mK, Gk);
multipoleTorque(mK, Tk);
// Need the torque on site I quadrupole due to site K multipole.
multipoleKPotentialAtI(mK, 2);
quadrupoleTorque(mI, Ti);
// Compute the potential due to a multipole component at site K.
quadrupoleKPotentialAtI(mK, 3);
double eI = multipoleEnergy(mI);
multipoleGradient(mI, Gi);
multipoleTorque(mI, Ti);
// Need the torque on site K quadrupole due to site I multipole.
multipoleIPotentialAtK(mI, 2);
quadrupoleTorque(mK, Tk);
double scale = c * 0.5;
Gi[0] = scale * (Gi[0] - Gk[0]);
Gi[1] = scale * (Gi[1] - Gk[1]);
Gi[2] = scale * (Gi[2] - Gk[2]);
Gk[0] = -Gi[0];
Gk[1] = -Gi[1];
Gk[2] = -Gi[2];
Ti[0] = scale * Ti[0];
Ti[1] = scale * Ti[1];
Ti[2] = scale * Ti[2];
Tk[0] = scale * Tk[0];
Tk[1] = scale * Tk[1];
Tk[2] = scale * Tk[2];
return scale * (eK + eI);
}
/**
* GK Permanent multipole Born grad.
*
* @param mI PolarizableMultipole at site I.
* @param mK PolarizableMultipole at site K.
* @return a double.
*/
public double multipoleEnergyBornGrad(PolarizableMultipole mI, PolarizableMultipole mK) {
return multipoleEnergy(mI, mK);
} | ||
| File | Project | Line |
|---|---|---|
| ffx/numerics/multipole/GKEnergyGlobal.java | Numerics | 76 |
| ffx/numerics/multipole/GKEnergyQI.java | Numerics | 77 |
gkQuadrupole = new GKTensorGlobal(QUADRUPOLE, quadrupoleOrder, gkSource, 1.0, epsilon);
}
/**
* Initialize the potential.
*
* @param r The separation. vector.
* @param r2 The squared separation.
* @param rbi The Born radius of atom i.
* @param rbk The Born radius of atom k.
*/
public void initPotential(double[] r, double r2, double rbi, double rbk) {
gkSource.generateSource(POTENTIAL, QUADRUPOLE, r2, rbi, rbk);
gkMonopole.setR(r);
gkDipole.setR(r);
gkQuadrupole.setR(r);
gkMonopole.generateTensor();
gkDipole.generateTensor();
gkQuadrupole.generateTensor();
}
/**
* Initialize for computing Born chain-rule terms.
*
* @param r The separation vector.
* @param r2 The squared separation.
* @param rbi The Born radius of atom i.
* @param rbk The Born radius of atom k.
*/
public void initBorn(double[] r, double r2, double rbi, double rbk) {
gkSource.generateSource(BORN, QUADRUPOLE, r2, rbi, rbk);
gkMonopole.setR(r);
gkDipole.setR(r);
gkQuadrupole.setR(r);
gkMonopole.generateTensor();
gkDipole.generateTensor();
gkQuadrupole.generateTensor();
}
/**
* Compute the multipole energy.
*
* @param mI The polarizable multipole of atom i.
* @param mK The polarizable multipole of atom k.
* @return The multipole energy.
*/
public double multipoleEnergy(PolarizableMultipole mI, PolarizableMultipole mK) {
double em = gkMonopole.multipoleEnergy(mI, mK);
double ed = gkDipole.multipoleEnergy(mI, mK);
double eq = gkQuadrupole.multipoleEnergy(mI, mK);
return em + ed + eq;
}
/**
* Compute the polarization energy.
*
* @param mI The polarizable multipole of atom i.
* @param mK The polarizable multipole of atom k.
* @return The polarization energy.
*/
public double polarizationEnergy(PolarizableMultipole mI, PolarizableMultipole mK) {
double emp = gkMonopole.polarizationEnergy(mI, mK);
double edp = gkDipole.polarizationEnergy(mI, mK);
double eqp = gkQuadrupole.polarizationEnergy(mI, mK);
return emp + edp + eqp;
}
/**
* Compute the multipole energy and gradient.
*
* @param mI The polarizable multipole of atom i.
* @param mK The polarizable multipole of atom k.
* @param gradI The gradient for atom i.
* @param torqueI The torque on atom i.
* @param torqueK The torque on atom k.
* @return The multipole energy.
*/
public double multipoleEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
double[] gradI, double[] torqueI, double[] torqueK) {
double[] gI = new double[3];
double[] gK = new double[3];
double[] tI = new double[3];
double[] tK = new double[3];
double em = gkMonopole.multipoleEnergyAndGradient(mI, mK, gI, gK, tI, tK);
for (int j = 0; j < 3; j++) {
gradI[j] = gI[j];
torqueI[j] = tI[j];
torqueK[j] = tK[j];
}
fill(gI, 0.0);
fill(gK, 0.0);
fill(tI, 0.0);
fill(tK, 0.0);
double ed = gkDipole.multipoleEnergyAndGradient(mI, mK, gI, gK, tI, tK);
for (int j = 0; j < 3; j++) {
gradI[j] += gI[j];
torqueI[j] += tI[j];
torqueK[j] += tK[j];
}
fill(gI, 0.0);
fill(gK, 0.0);
fill(tI, 0.0);
fill(tK, 0.0);
double eq = gkQuadrupole.multipoleEnergyAndGradient(mI, mK, gI, gK, tI, tK);
for (int j = 0; j < 3; j++) {
gradI[j] += gI[j];
torqueI[j] += tI[j];
torqueK[j] += tK[j];
}
return em + ed + eq;
}
/**
* Compute the polarization energy and gradient.
*
* @param mI The polarizable multipole of atom i.
* @param mK The polarizable multipole of atom k.
* @param mutualMask The mutual polarization mask.
* @param gradI The gradient for atom i.
* @param torqueI The torque on atom i.
* @param torqueK The torque on atom k.
* @return The polarization energy.
*/
public double polarizationEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK, double mutualMask,
double[] gradI, double[] torqueI, double[] torqueK) {
double[] gI = new double[3];
double[] tI = new double[3];
double[] tK = new double[3];
double emp = gkMonopole.polarizationEnergyAndGradient(mI, mK, 1.0, 1.0, mutualMask, gI, tI, tK);
for (int j = 0; j < 3; j++) {
gradI[j] = gI[j];
torqueI[j] = tI[j];
torqueK[j] = tK[j];
}
fill(gI, 0.0);
fill(tI, 0.0);
fill(tK, 0.0);
double edp = gkDipole.polarizationEnergyAndGradient(mI, mK, 1.0, 1.0, mutualMask, gI, tI, tK);
for (int j = 0; j < 3; j++) {
gradI[j] += gI[j];
torqueI[j] += tI[j];
torqueK[j] += tK[j];
}
fill(gI, 0.0);
fill(tI, 0.0);
fill(tK, 0.0);
double eqp = gkQuadrupole.polarizationEnergyAndGradient(mI, mK, 1.0, 1.0, mutualMask, gI, tI, tK);
for (int j = 0; j < 3; j++) {
gradI[j] += gI[j];
torqueI[j] += tI[j];
torqueK[j] += tK[j];
}
// Sum the GK polarization interaction energy.
return emp + edp + eqp;
}
/**
* Compute the Born chain-rule term for the multipole energy.
*
* @param mI The polarizable multipole of atom i.
* @param mK The polarizable multipole of atom k.
* @return The Born chain-rule term.
*/
public double multipoleEnergyBornGrad(PolarizableMultipole mI, PolarizableMultipole mK) {
double db = gkMonopole.multipoleEnergyBornGrad(mI, mK);
db += gkDipole.multipoleEnergyBornGrad(mI, mK);
db += gkQuadrupole.multipoleEnergyBornGrad(mI, mK);
return db;
}
/**
* Compute the Born chain-rule term for the polarization energy.
*
* @param mI The polarizable multipole of atom i.
* @param mK The polarizable multipole of atom k.
* @param mutual If true, compute the mutual polarization contribution.
* @return The Born chain-rule term.
*/
public double polarizationEnergyBornGrad(PolarizableMultipole mI, PolarizableMultipole mK,
boolean mutual) {
// Compute the GK polarization Born chain-rule term.
double db = gkMonopole.polarizationEnergyBornGrad(mI, mK);
db += gkDipole.polarizationEnergyBornGrad(mI, mK);
db += gkQuadrupole.polarizationEnergyBornGrad(mI, mK);
// Add the mutual polarization contribution to Born chain-rule term.
if (mutual) {
db += gkDipole.mutualPolarizationEnergyBornGrad(mI, mK);
}
return db;
}
} | ||
| File | Project | Line |
|---|---|---|
| ffx/numerics/multipole/GKTensorGlobal.java | Numerics | 307 |
| ffx/numerics/multipole/GKTensorQI.java | Numerics | 305 |
public double polarizationEnergy(PolarizableMultipole mI, PolarizableMultipole mK, DoubleVector scaleEnergy) {
return polarizationEnergy(mI, mK);
}
/**
* GK Polarization Energy.
*
* @param mI PolarizableMultipole at site I.
* @param mK PolarizableMultipole at site K.
* @return a double.
*/
public double polarizationEnergy(PolarizableMultipole mI, PolarizableMultipole mK) {
return switch (multipoleOrder) {
default -> {
// Find the GK charge potential of site I at site K.
chargeIPotentialAtK(mI, 1);
// Energy of induced dipole K in the field of permanent charge I.
double eK = polarizationEnergy(mK);
// Find the GK charge potential of site K at site I.
chargeKPotentialAtI(mK, 1);
// Energy of induced dipole I in the field of permanent charge K.
double eI = polarizationEnergy(mI);
yield c * 0.5 * (eK + eI);
}
case DIPOLE -> {
// Find the GK dipole potential of site I at site K.
dipoleIPotentialAtK(mI.dx, mI.dy, mI.dz, 1);
// Energy of induced dipole K in the field of permanent dipole I.
double eK = polarizationEnergy(mK);
// Find the GK induced dipole potential of site I at site K.
dipoleIPotentialAtK(mI.ux, mI.uy, mI.uz, 2);
// Energy of permanent multipole K in the field of induced dipole I.
eK += 0.5 * multipoleEnergy(mK);
// Find the GK dipole potential of site K at site I.
dipoleKPotentialAtI(mK.dx, mK.dy, mK.dz, 1);
// Energy of induced dipole I in the field of permanent dipole K.
double eI = polarizationEnergy(mI);
// Find the GK induced dipole potential of site K at site I.
dipoleKPotentialAtI(mK.ux, mK.uy, mK.uz, 2);
// Energy of permanent multipole I in the field of induced dipole K.
eI += 0.5 * multipoleEnergy(mI);
yield c * 0.5 * (eK + eI);
}
case QUADRUPOLE -> {
// Find the GK quadrupole potential of site I at site K.
quadrupoleIPotentialAtK(mI, 1);
// Energy of induced dipole K in the field of permanent quadrupole I.
double eK = polarizationEnergy(mK);
// Find the GK quadrupole potential of site K at site I.
quadrupoleKPotentialAtI(mK, 1);
// Energy of induced dipole I in the field of permanent quadrupole K.
double eI = polarizationEnergy(mI);
yield c * 0.5 * (eK + eI);
}
};
}
/**
* GK Polarization Energy.
*
* @param mI PolarizableMultipole at site I.
* @param mK PolarizableMultipole at site K.
* @return a double.
*/
public double polarizationEnergyBorn(PolarizableMultipole mI, PolarizableMultipole mK) {
return switch (multipoleOrder) {
default -> {
// Find the GK charge potential of site I at site K.
chargeIPotentialAtK(mI, 1);
// Energy of induced dipole K in the field of permanent charge I.
double eK = polarizationEnergyS(mK);
// Find the GK charge potential of site K at site I.
chargeKPotentialAtI(mK, 1);
// Energy of induced dipole I in the field of permanent charge K.
double eI = polarizationEnergyS(mI);
yield c * 0.5 * (eK + eI);
}
case DIPOLE -> {
// Find the GK dipole potential of site I at site K.
dipoleIPotentialAtK(mI.dx, mI.dy, mI.dz, 1);
// Energy of induced dipole K in the field of permanent dipole I.
double eK = polarizationEnergyS(mK);
// Find the GK induced dipole potential of site I at site K.
dipoleIPotentialAtK(mI.sx, mI.sy, mI.sz, 2);
// Energy of permanent multipole K in the field of induced dipole I.
eK += 0.5 * multipoleEnergy(mK);
// Find the GK dipole potential of site K at site I.
dipoleKPotentialAtI(mK.dx, mK.dy, mK.dz, 1);
// Energy of induced dipole I in the field of permanent dipole K.
double eI = polarizationEnergyS(mI);
// Find the GK induced dipole potential of site K at site I.
dipoleKPotentialAtI(mK.sx, mK.sy, mK.sz, 2);
// Energy of permanent multipole I in the field of induced dipole K.
eI += 0.5 * multipoleEnergy(mI);
yield c * 0.5 * (eK + eI);
}
case QUADRUPOLE -> {
// Find the GK quadrupole potential of site I at site K.
quadrupoleIPotentialAtK(mI, 1);
// Energy of induced dipole K in the field of permanent quadrupole I.
double eK = polarizationEnergyS(mK);
// Find the GK quadrupole potential of site K at site I.
quadrupoleKPotentialAtI(mK, 1);
// Energy of induced dipole I in the field of permanent quadrupole K.
double eI = polarizationEnergyS(mI);
yield c * 0.5 * (eK + eI);
}
};
}
/**
* Polarization Energy and Gradient.
*
* @param mI PolarizableMultipole at site I.
* @param mK PolarizableMultipole at site K.
* @param inductionMask This is ignored, since masking/scaling is not applied to GK
* interactions (everything is intermolecular).
* @param energyMask This is ignored, since masking/scaling is not applied to GK interactions
* (everything is intermolecular).
* @param mutualMask This should be set to zero for direction polarization.
* @param Gi an array of {@link double} objects.
* @param Ti an array of {@link double} objects.
* @param Tk an array of {@link double} objects.
* @return a double.
*/
@Override
public double polarizationEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
double inductionMask, double energyMask, double mutualMask, double[] Gi, double[] Ti,
double[] Tk) {
return switch (multipoleOrder) {
default -> monopolePolarizationEnergyAndGradient(mI, mK, Gi);
case DIPOLE -> dipolePolarizationEnergyAndGradient(mI, mK, mutualMask, Gi, Ti, Tk);
case QUADRUPOLE -> quadrupolePolarizationEnergyAndGradient(mI, mK, Gi, Ti, Tk);
};
}
/**
* Monopole Polarization Energy and Gradient.
*
* @param mI PolarizableMultipole at site I.
* @param mK PolarizableMultipole at site K.
* @param Gi an array of {@link double} objects.
* @return a double.
*/
public double monopolePolarizationEnergyAndGradient(PolarizableMultipole mI,
PolarizableMultipole mK, double[] Gi) {
// Find the permanent multipole potential at site k.
chargeIPotentialAtK(mI, 2);
// Energy of induced dipole k in the field of multipole i.
double eK = polarizationEnergy(mK);
// Derivative with respect to moving atom k.
Gi[0] = -(mK.sx * E200 + mK.sy * E110 + mK.sz * E101);
Gi[1] = -(mK.sx * E110 + mK.sy * E020 + mK.sz * E011);
Gi[2] = -(mK.sx * E101 + mK.sy * E011 + mK.sz * E002);
// Find the permanent multipole potential and derivatives at site i.
chargeKPotentialAtI(mK, 2);
// Energy of induced dipole i in the field of multipole k.
double eI = polarizationEnergy(mI);
// Derivative with respect to moving atom i.
Gi[0] += (mI.sx * E200 + mI.sy * E110 + mI.sz * E101);
Gi[1] += (mI.sx * E110 + mI.sy * E020 + mI.sz * E011);
Gi[2] += (mI.sx * E101 + mI.sy * E011 + mI.sz * E002);
double scale = c * 0.5;
Gi[0] *= scale;
Gi[1] *= scale;
Gi[2] *= scale;
// Total polarization energy.
return scale * (eI + eK);
}
/**
* Dipole Polarization Energy and Gradient.
*
* @param mI PolarizableMultipole at site I.
* @param mK PolarizableMultipole at site K.
* @param mutualMask This should be set to zero for direction polarization.
* @param Gi an array of {@link double} objects.
* @param Ti an array of {@link double} objects.
* @param Tk an array of {@link double} objects.
* @return a double.
*/
public double dipolePolarizationEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
double mutualMask, double[] Gi, double[] Ti, double[] Tk) {
// Find the permanent multipole potential and derivatives at site k.
dipoleIPotentialAtK(mI.dx, mI.dy, mI.dz, 2);
// Energy of induced dipole k in the field of multipole i.
double eK = polarizationEnergy(mK);
// Derivative with respect to moving atom k.
Gi[0] = -(mK.sx * E200 + mK.sy * E110 + mK.sz * E101);
Gi[1] = -(mK.sx * E110 + mK.sy * E020 + mK.sz * E011);
Gi[2] = -(mK.sx * E101 + mK.sy * E011 + mK.sz * E002); | ||
| File | Project | Line |
|---|---|---|
| ffx/numerics/multipole/GKTensorGlobal.java | Numerics | 533 |
| ffx/numerics/multipole/GKTensorQI.java | Numerics | 535 |
eI += 0.5 * multipoleEnergy(mI);
G = new double[3];
multipoleGradient(mI, G);
Gi[0] += G[0];
Gi[1] += G[1];
Gi[2] += G[2];
multipoleTorque(mI, Ti);
// Get the induced-induced portion of the force (Ud . dC/dX . Up).
// This contribution does not exist for direct polarization (mutualMask == 0.0).
if (mutualMask != 0.0) {
// Find the potential and its derivatives at k due to induced dipole i.
dipoleIPotentialAtK(mI.ux, mI.uy, mI.uz, 2);
Gi[0] -= mutualMask * (mK.px * E200 + mK.py * E110 + mK.pz * E101);
Gi[1] -= mutualMask * (mK.px * E110 + mK.py * E020 + mK.pz * E011);
Gi[2] -= mutualMask * (mK.px * E101 + mK.py * E011 + mK.pz * E002);
// Find the potential and its derivatives at i due to induced dipole k.
dipoleKPotentialAtI(mK.ux, mK.uy, mK.uz, 2);
Gi[0] += mutualMask * (mI.px * E200 + mI.py * E110 + mI.pz * E101);
Gi[1] += mutualMask * (mI.px * E110 + mI.py * E020 + mI.pz * E011);
Gi[2] += mutualMask * (mI.px * E101 + mI.py * E011 + mI.pz * E002);
}
// Total polarization energy.
double scale = c * 0.5;
double energy = scale * (eI + eK);
Gi[0] *= scale;
Gi[1] *= scale;
Gi[2] *= scale;
Ti[0] *= scale;
Ti[1] *= scale;
Ti[2] *= scale;
Tk[0] *= scale;
Tk[1] *= scale;
Tk[2] *= scale;
return energy;
}
/**
* Quadrupole Polarization Energy and Gradient.
*
* @param mI PolarizableMultipole at site I.
* @param mK PolarizableMultipole at site K.
* @param Gi an array of {@link double} objects.
* @param Ti an array of {@link double} objects.
* @param Tk an array of {@link double} objects.
* @return a double.
*/
public double quadrupolePolarizationEnergyAndGradient(PolarizableMultipole mI, PolarizableMultipole mK,
double[] Gi, double[] Ti, double[] Tk) {
// Find the permanent multipole potential and derivatives at site k.
quadrupoleIPotentialAtK(mI, 2);
// Energy of induced dipole k in the field of multipole i.
double eK = polarizationEnergy(mK);
// Derivative with respect to moving atom k.
Gi[0] = -(mK.sx * E200 + mK.sy * E110 + mK.sz * E101);
Gi[1] = -(mK.sx * E110 + mK.sy * E020 + mK.sz * E011);
Gi[2] = -(mK.sx * E101 + mK.sy * E011 + mK.sz * E002);
// Find the permanent multipole potential and derivatives at site i.
quadrupoleKPotentialAtI(mK, 2);
// Energy of induced dipole i in the field of multipole k.
double eI = polarizationEnergy(mI);
// Derivative with respect to moving atom i.
Gi[0] += (mI.sx * E200 + mI.sy * E110 + mI.sz * E101);
Gi[1] += (mI.sx * E110 + mI.sy * E020 + mI.sz * E011);
Gi[2] += (mI.sx * E101 + mI.sy * E011 + mI.sz * E002);
double scale = c * 0.5;
Gi[0] *= scale;
Gi[1] *= scale;
Gi[2] *= scale;
// Find the potential and its derivatives at K due to the averaged induced dipole at site i.
dipoleIPotentialAtK(scale * mI.sx, scale * mI.sy, scale * mI.sz, 2);
quadrupoleTorque(mK, Tk);
// Find the potential and its derivatives at I due to the averaged induced dipole at k.
dipoleKPotentialAtI(scale * mK.sx, scale * mK.sy, scale * mK.sz, 2);
quadrupoleTorque(mI, Ti);
// Total polarization energy.
return scale * (eI + eK);
}
/**
* GK Direct Polarization Born grad.
*
* @param mI PolarizableMultipole at site I.
* @param mK PolarizableMultipole at site K.
* @return Partial derivative of the Polarization energy with respect to a Born grad.
*/
public double polarizationEnergyBornGrad(PolarizableMultipole mI, PolarizableMultipole mK) {
return 2.0 * polarizationEnergyBorn(mI, mK);
}
/**
* GK Mutual Polarization Contribution to the Born grad.
*
* @param mI PolarizableMultipole at site I.
* @param mK PolarizableMultipole at site K.
* @return Mutual Polarization contribution to the partial derivative with respect to a Born grad.
*/
public double mutualPolarizationEnergyBornGrad(PolarizableMultipole mI, PolarizableMultipole mK) {
double db = 0.0;
if (multipoleOrder == GKMultipoleOrder.DIPOLE) {
// Find the potential and its derivatives at k due to induced dipole i.
dipoleIPotentialAtK(mI.ux, mI.uy, mI.uz, 2);
db = 0.5 * (mK.px * E100 + mK.py * E010 + mK.pz * E001);
// Find the potential and its derivatives at i due to induced dipole k.
dipoleKPotentialAtI(mK.ux, mK.uy, mK.uz, 2);
db += 0.5 * (mI.px * E100 + mI.py * E010 + mI.pz * E001);
}
return c * db;
}
/**
* Generate source terms for the Kirkwood version of the Challacombe et al. recursion.
*/
@Override
protected void source(double[] work) {
gkSource.source(work, multipoleOrder);
}
} | ||
| File | Project | Line |
|---|---|---|
| ffx/potential/bonded/RotamerLibrary.java | Potential | 2196 |
| ffx/potential/bonded/RotamerLibrary.java | Potential | 2266 |
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 | 6270 |
| ffx/crystal/SpaceGroupDefinitions.java | Crystal | 6594 |
"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/numerics/multipole/GKEnergyGlobalSIMD.java | Numerics | 79 |
| ffx/numerics/multipole/GKEnergyQISIMD.java | Numerics | 80 |
gkQuadrupole = new GKTensorGlobalSIMD(QUADRUPOLE, quadrupoleOrder, gkSource, 1.0, epsilon);
}
/**
* Initialize the potential.
*
* @param r The separation. vector.
* @param r2 The squared separation.
* @param rbi The Born radius of atom i.
* @param rbk The Born radius of atom k.
*/
public void initPotential(DoubleVector[] r, DoubleVector r2, DoubleVector rbi, DoubleVector rbk) {
gkSource.generateSource(POTENTIAL, QUADRUPOLE, r2, rbi, rbk);
gkMonopole.setR(r);
gkDipole.setR(r);
gkQuadrupole.setR(r);
gkMonopole.generateTensor();
gkDipole.generateTensor();
gkQuadrupole.generateTensor();
}
/**
* Initialize for computing Born chain-rule terms.
*
* @param r The separation vector.
* @param r2 The squared separation.
* @param rbi The Born radius of atom i.
* @param rbk The Born radius of atom k.
*/
public void initBorn(DoubleVector[] r, DoubleVector r2, DoubleVector rbi, DoubleVector rbk) {
gkSource.generateSource(BORN, QUADRUPOLE, r2, rbi, rbk);
gkMonopole.setR(r);
gkDipole.setR(r);
gkQuadrupole.setR(r);
gkMonopole.generateTensor();
gkDipole.generateTensor();
gkQuadrupole.generateTensor();
}
/**
* Compute the multipole energy.
*
* @param mI The polarizable multipole of atom i.
* @param mK The polarizable multipole of atom k.
* @return The multipole energy.
*/
public DoubleVector multipoleEnergy(PolarizableMultipoleSIMD mI, PolarizableMultipoleSIMD mK) {
DoubleVector em = gkMonopole.multipoleEnergy(mI, mK);
DoubleVector ed = gkDipole.multipoleEnergy(mI, mK);
DoubleVector eq = gkQuadrupole.multipoleEnergy(mI, mK);
return em.add(ed).add(eq);
}
/**
* Compute the polarization energy.
*
* @param mI The polarizable multipole of atom i.
* @param mK The polarizable multipole of atom k.
* @return The polarization energy.
*/
public DoubleVector polarizationEnergy(PolarizableMultipoleSIMD mI, PolarizableMultipoleSIMD mK) {
DoubleVector emp = gkMonopole.polarizationEnergy(mI, mK);
DoubleVector edp = gkDipole.polarizationEnergy(mI, mK);
DoubleVector eqp = gkQuadrupole.polarizationEnergy(mI, mK);
return emp.add(edp).add(eqp);
}
/**
* Compute the multipole energy and gradient.
*
* @param mI The polarizable multipole of atom i.
* @param mK The polarizable multipole of atom k.
* @param gI The gradient for atom i.
* @param tI The torque on atom i.
* @param tK The torque on atom k.
* @return The multipole energy.
*/
public DoubleVector multipoleEnergyAndGradient(PolarizableMultipoleSIMD mI, PolarizableMultipoleSIMD mK,
DoubleVector[] gI, DoubleVector[] tI, DoubleVector[] tK) {
DoubleVector[] gK = new DoubleVector[3];
DoubleVector em = gkMonopole.multipoleEnergyAndGradient(mI, mK, gI, gK, tI, tK);
DoubleVector ed = gkDipole.multipoleEnergyAndGradient(mI, mK, gI, gK, tI, tK);
DoubleVector eq = gkQuadrupole.multipoleEnergyAndGradient(mI, mK, gI, gK, tI, tK);
return em.add(ed).add(eq);
}
/**
* Compute the polarization energy and gradient.
*
* @param mI The polarizable multipole of atom i.
* @param mK The polarizable multipole of atom k.
* @param mutualMask The mutual polarization mask.
* @param gI The gradient for atom i.
* @param tI The torque on atom i.
* @param tK The torque on atom k.
* @return The polarization energy.
*/
public DoubleVector polarizationEnergyAndGradient(PolarizableMultipoleSIMD mI, PolarizableMultipoleSIMD mK, DoubleVector mutualMask,
DoubleVector[] gI, DoubleVector[] tI, DoubleVector[] tK) {
DoubleVector emp = gkMonopole.polarizationEnergyAndGradient(mI, mK, one, one, mutualMask, gI, tI, tK);
DoubleVector edp = gkDipole.polarizationEnergyAndGradient(mI, mK, one, one, mutualMask, gI, tI, tK);
DoubleVector eqp = gkQuadrupole.polarizationEnergyAndGradient(mI, mK, one, one, mutualMask, gI, tI, tK);
// Sum the GK polarization interaction energy.
return emp.add(edp).add(eqp);
}
/**
* Compute the Born chain-rule term for the multipole energy.
*
* @param mI The polarizable multipole of atom i.
* @param mK The polarizable multipole of atom k.
* @return The Born chain-rule term.
*/
public DoubleVector multipoleEnergyBornGrad(PolarizableMultipoleSIMD mI, PolarizableMultipoleSIMD mK) {
DoubleVector db = gkMonopole.multipoleEnergyBornGrad(mI, mK);
db = db.add(gkDipole.multipoleEnergyBornGrad(mI, mK));
db = db.add(gkQuadrupole.multipoleEnergyBornGrad(mI, mK));
return db;
}
/**
* Compute the Born chain-rule term for the polarization energy.
*
* @param mI The polarizable multipole of atom i.
* @param mK The polarizable multipole of atom k.
* @param mutual If true, compute the mutual polarization contribution.
* @return The Born chain-rule term.
*/
public DoubleVector polarizationEnergyBornGrad(PolarizableMultipoleSIMD mI, PolarizableMultipoleSIMD mK, boolean mutual) {
// Compute the GK polarization Born chain-rule term.
DoubleVector db = gkMonopole.polarizationEnergyBornGrad(mI, mK);
db = db.add(gkDipole.polarizationEnergyBornGrad(mI, mK));
db = db.add(gkQuadrupole.polarizationEnergyBornGrad(mI, mK));
// Add the mutual polarization contribution to Born chain-rule term.
if (mutual) {
db = db.add(gkDipole.mutualPolarizationEnergyBornGrad(mI, mK));
}
return db;
}
} | ||
| File | Project | Line |
|---|---|---|
| ffx/crystal/SpaceGroupDefinitions.java | Crystal | 6394 |
| ffx/crystal/SpaceGroupDefinitions.java | Crystal | 7210 |
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/crystal/SpaceGroupDefinitions.java | Crystal | 6154 |
| ffx/crystal/SpaceGroupDefinitions.java | Crystal | 6394 |
| ffx/crystal/SpaceGroupDefinitions.java | Crystal | 7210 |
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 | Refinement | 1997 |
| ffx/xray/CrystalReciprocalSpace.java | Refinement | 2232 |
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 | 698 |
| ffx/potential/nonbonded/implicit/PermanentGKFieldRegion.java | Potential | 306 |
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/parsers/XPHFilter.java | Potential | 345 |
| ffx/potential/parsers/XYZFilter.java | Potential | 332 |
logger.info(format(" Opening %s with %d atoms\n", xphFile.getName(), numberOfAtoms));
remarkLine = data.trim();
// The header line is reasonable. Check for periodic box dimensions.
br.mark(10000);
data = br.readLine();
if (!readPBC(data, activeMolecularAssembly)) {
br.reset();
}
// Prepare to parse atom lines.
HashMap<Integer, Integer> labelHash = new HashMap<>();
int[] label = new int[numberOfAtoms];
int[][] bonds = new int[numberOfAtoms][8];
double[][] d = new double[numberOfAtoms][3];
boolean renumber = false;
atomList = new ArrayList<>();
// Loop over the expected number of atoms.
for (int i = 0; i < numberOfAtoms; i++) {
if (!br.ready()) {
return false;
}
data = br.readLine();
if (data == null) {
logger.warning(
format(" Check atom %d in %s.", (i + 1), activeMolecularAssembly.getFile().getName()));
return false;
}
tokens = data.trim().split(" +");
if (tokens.length < 6) {
logger.warning(
format(" Check atom %d in %s.", (i + 1), activeMolecularAssembly.getFile().getName()));
return false;
}
// Valid number of tokens, so try to parse this line.
label[i] = parseInt(tokens[0]);
// Check for valid atom numbering, or flag for re-numbering.
if (label[i] != i + 1) {
renumber = true;
}
String atomName = tokens[1];
d[i][0] = parseDouble(tokens[2]);
d[i][1] = parseDouble(tokens[3]);
d[i][2] = parseDouble(tokens[4]);
int type = parseInt(tokens[5]);
AtomType atomType = forceField.getAtomType(Integer.toString(type));
if (atomType == null) {
StringBuilder message = new StringBuilder("Check atom type ");
message.append(type).append(" for Atom ").append(i + 1);
message.append(" in ").append(activeMolecularAssembly.getFile().getName());
logger.warning(message.toString());
return false;
}
Atom a = new Atom(i + 1, atomName, atomType, d[i]);
atomList.add(a);
// Bond Data
int numberOfBonds = tokens.length - 6;
for (int b = 0; b < 8; b++) {
if (b < numberOfBonds) {
int bond = parseInt(tokens[6 + b]);
bonds[i][b] = bond;
} else {
bonds[i][b] = 0;
}
}
}
// Check if this is an archive.
if (br.ready()) {
// Read past blank lines between archive files
data = br.readLine().trim();
while (data.equals("") && br.ready()) { | ||
| 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, | ||