Force Field XCPD Results
The following document contains the results of PMD's CPD 7.17.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 double values.
* @param Ti an array of double values.
* @param Tk an array of double values.
* @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 double values.
* @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 double values.
* @param Ti an array of double values.
* @param Tk an array of double values.
* @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 double values.
* @param Ti an array of double values.
* @param Tk an array of double values.
* @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 | 302 |
| 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/fft/MixedRadixFactor7.java | Numerics | 499 |
| ffx/numerics/fft/MixedRadixFactor7.java | Numerics | 580 |
for (int k1 = 0; k1 < innerLoopLimit; k1 += BLOCK_LOOP, i += LENGTH, j += LENGTH) {
final DoubleVector
z0r = fromArray(DOUBLE_SPECIES, data, i),
z1r = fromArray(DOUBLE_SPECIES, data, i + di),
z2r = fromArray(DOUBLE_SPECIES, data, i + di2),
z3r = fromArray(DOUBLE_SPECIES, data, i + di3),
z4r = fromArray(DOUBLE_SPECIES, data, i + di4),
z5r = fromArray(DOUBLE_SPECIES, data, i + di5),
z6r = fromArray(DOUBLE_SPECIES, data, i + di6),
z0i = fromArray(DOUBLE_SPECIES, data, i + im),
z1i = fromArray(DOUBLE_SPECIES, data, i + di + im),
z2i = fromArray(DOUBLE_SPECIES, data, i + di2 + im),
z3i = fromArray(DOUBLE_SPECIES, data, i + di3 + im),
z4i = fromArray(DOUBLE_SPECIES, data, i + di4 + im),
z5i = fromArray(DOUBLE_SPECIES, data, i + di5 + im),
z6i = fromArray(DOUBLE_SPECIES, data, i + di6 + im);
final DoubleVector
t0r = z1r.add(z6r), t0i = z1i.add(z6i),
t1r = z1r.sub(z6r), t1i = z1i.sub(z6i),
t2r = z2r.add(z5r), t2i = z2i.add(z5i),
t3r = z2r.sub(z5r), t3i = z2i.sub(z5i),
t4r = z4r.add(z3r), t4i = z4i.add(z3i),
t5r = z4r.sub(z3r), t5i = z4i.sub(z3i),
t6r = t2r.add(t0r), t6i = t2i.add(t0i),
t7r = t5r.add(t3r), t7i = t5i.add(t3i);
final DoubleVector
b0r = z0r.add(t6r).add(t4r), b0i = z0i.add(t6i).add(t4i),
b1r = t6r.add(t4r).mul(v1), b1i = t6i.add(t4i).mul(v1),
b2r = t0r.sub(t4r).mul(v2), b2i = t0i.sub(t4i).mul(v2),
b3r = t4r.sub(t2r).mul(v3), b3i = t4i.sub(t2i).mul(v3),
b4r = t2r.sub(t0r).mul(v4), b4i = t2i.sub(t0i).mul(v4),
b5r = t7r.add(t1r).mul(v5), b5i = t7i.add(t1i).mul(v5),
b6r = t1r.sub(t5r).mul(v6), b6i = t1i.sub(t5i).mul(v6),
b7r = t5r.sub(t3r).mul(v7), b7i = t5i.sub(t3i).mul(v7),
b8r = t3r.sub(t1r).mul(v8), b8i = t3i.sub(t1i).mul(v8);
final DoubleVector
u0r = b0r.add(b1r), u0i = b0i.add(b1i),
u1r = b2r.add(b3r), u1i = b2i.add(b3i),
u2r = b4r.sub(b3r), u2i = b4i.sub(b3i),
u3r = b2r.add(b4r).neg(), u3i = b2i.add(b4i).neg(),
u4r = b6r.add(b7r), u4i = b6i.add(b7i),
u5r = b8r.sub(b7r), u5i = b8i.sub(b7i),
u6r = b8r.add(b6r).neg(), u6i = b8i.add(b6i).neg(),
u7r = u0r.add(u1r), u7i = u0i.add(u1i),
u8r = u0r.add(u2r), u8i = u0i.add(u2i),
u9r = u0r.add(u3r), u9i = u0i.add(u3i),
u10r = u4r.add(b5r), u10i = u4i.add(b5i),
u11r = u5r.add(b5r), u11i = u5i.add(b5i),
u12r = u6r.add(b5r), u12i = u6i.add(b5i);
b0r.intoArray(ret, j);
b0i.intoArray(ret, j + im); | ||
| 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 double values.
* @param Ti an array of double values.
* @param Tk an array of double values.
* @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 double values.
* @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 double values.
* @param Ti an array of double values.
* @param Tk an array of double values.
* @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/realspace/commands/test/Alchemical.java | Refinement | 192 |
| ffx/xray/commands/test/Alchemical.java | Refinement | 203 |
File structureFile = new File(FilenameUtils.normalize(modelFilename));
structureFile = new File(structureFile.getAbsolutePath());
String baseFilename = FilenameUtils.removeExtension(structureFile.getName());
File histogramRestart = new File(baseFilename + ".his");
File lambdaRestart = new File(baseFilename + ".lam");
File dyn = new File(baseFilename + ".dyn");
Comm world = Comm.world();
int size = world.size();
int rank = 0;
double[] energyArray = new double[world.size()];
for (int i = 0; i < world.size(); i++) {
energyArray[i] = Double.MAX_VALUE;
}
// For a multi-process job, try to get the restart files from rank sub-directories.
if (size > 1) {
rank = world.rank();
File rankDirectory = new File(
structureFile.getParent() + File.separator + Integer.toString(rank));
if (!rankDirectory.exists()) {
rankDirectory.mkdir();
}
lambdaRestart = new File(rankDirectory.getPath() + File.separator + baseFilename + ".lam");
dyn = new File(rankDirectory.getPath() + File.separator + baseFilename + ".dyn");
structureFile = new File(rankDirectory.getPath() + File.separator + structureFile.getName());
}
if (!dyn.exists()) {
dyn = null;
}
// Set built atoms active/use flags to true (false for other atoms).
Atom[] atoms = activeAssembly.getAtomArray();
// Get a reference to the first system's ForceFieldEnergy.
ForceFieldEnergy forceFieldEnergy = activeAssembly.getPotentialEnergy();
forceFieldEnergy.setPrintOnFailure(false, false);
// Configure all atoms to:
// 1) be used in the potential
// 2) be inactive (i.e. cannot move)
// 3) not be controlled by the lambda state variable.
for (int i = 0; i <= atoms.length; i++) {
Atom ai = atoms[i - 1];
ai.setUse(true);
ai.setActive(false);
ai.setApplyLambda(false);
}
double crystalCharge = activeAssembly.getCharge(true);
logger.info(" Overall crystal charge: " + crystalCharge);
List<MSNode> ions = assemblies[0].getIons();
List<MSNode> water = assemblies[0].getWater();
// Consider the option of creating a composite lambda gradient from vapor phase to crystal phase
if (!onlyWater) {
logger.info("Doing ions.");
if (ions == null || ions.size() == 0) {
logger.info("\n Please add an ion to the PDB file to scan with.");
return this;
}
for (MSNode msNode : ions) {
// Check if ionType is specified and matches this ion's atom names
boolean ionSelected = false;
if (ionType != null) {
logger.info("Selecting ion.");
List<Atom> atomList = msNode.getAtomList();
if (!atomList.isEmpty()) {
String atomName = atomList.get(0).getName();
for (String targetIonType : ionType) {
if (atomName.equals(targetIonType)) {
ionSelected = true;
break;
}
}
}
}
if (ionSelected) {
logger.info("Ion has been selected.");
for (Atom atom : msNode.getAtomList()) {
System.out.println("Activating ions");
atom.setUse(true);
atom.setActive(true);
atom.setApplyLambda(true);
logger.info(" Alchemical atom: " + atom.toString());
}
} else {
logger.info("Ion has not been selected.");
if (neutralize) {
logger.info("Neutralizing crystal.");
double ionCharge = 0;
for (Atom atom : msNode.getAtomList()) {
ionCharge += atom.getMultipoleType().getCharge();
}
logger.info("Ion charge is: " + Double.toString(ionCharge));
int numIons = (int) (-1 * (Math.ceil(crystalCharge / ionCharge)));
if (numIons > 0) {
List<Atom> atomList = msNode.getAtomList();
String atomName = atomList.isEmpty() ? "" : atomList.get(0).getName();
logger.info(numIons + " " + atomName
+ " ions needed to neutralize the crystal.");
ionType = new String[]{atomName};
for (Atom atom : msNode.getAtomList()) {
atom.setUse(true);
atom.setActive(true);
atom.setApplyLambda(true);
logger.info(" Alchemical atom: " + atom.toString());
}
}
} else {
List<Atom> atomList = msNode.getAtomList();
String atomName = atomList.isEmpty() ? "" : atomList.get(0).getName();
ionType = new String[]{atomName};
for (Atom atom : msNode.getAtomList()) {
atom.setUse(true);
atom.setActive(true);
atom.setApplyLambda(true);
logger.info(" Alchemical atom: " + atom.toString());
}
}
}
}
}
// Lambdize water for position optimization
if (!onlyIons) { | ||
| 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 double values.
* @param Ti an array of double values.
* @param Tk an array of double values.
* @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/numerics/fft/MixedRadixFactor7.java | Numerics | 117 |
| ffx/numerics/fft/MixedRadixFactor7.java | Numerics | 223 |
for (int k1 = 0; k1 < innerLoopLimit; k1++, i += ii, j += ii) {
final double z0r = data[i];
final double z1r = data[i + di];
final double z2r = data[i + di2];
final double z3r = data[i + di3];
final double z4r = data[i + di4];
final double z5r = data[i + di5];
final double z6r = data[i + di6];
final double z0i = data[i + im];
final double z1i = data[i + di + im];
final double z2i = data[i + di2 + im];
final double z3i = data[i + di3 + im];
final double z4i = data[i + di4 + im];
final double z5i = data[i + di5 + im];
final double z6i = data[i + di6 + im];
final double t0r = z1r + z6r;
final double t0i = z1i + z6i;
final double t1r = z1r - z6r;
final double t1i = z1i - z6i;
final double t2r = z2r + z5r;
final double t2i = z2i + z5i;
final double t3r = z2r - z5r;
final double t3i = z2i - z5i;
final double t4r = z4r + z3r;
final double t4i = z4i + z3i;
final double t5r = z4r - z3r;
final double t5i = z4i - z3i;
final double t6r = t2r + t0r;
final double t6i = t2i + t0i;
final double t7r = t5r + t3r;
final double t7i = t5i + t3i;
final double b0r = z0r + t6r + t4r;
final double b0i = z0i + t6i + t4i;
final double b1r = v1 * (t6r + t4r);
final double b1i = v1 * (t6i + t4i);
final double b2r = v2 * (t0r - t4r);
final double b2i = v2 * (t0i - t4i);
final double b3r = v3 * (t4r - t2r);
final double b3i = v3 * (t4i - t2i);
final double b4r = v4 * (t2r - t0r);
final double b4i = v4 * (t2i - t0i);
final double b5r = v5 * (t7r + t1r);
final double b5i = v5 * (t7i + t1i);
final double b6r = v6 * (t1r - t5r);
final double b6i = v6 * (t1i - t5i);
final double b7r = v7 * (t5r - t3r);
final double b7i = v7 * (t5i - t3i);
final double b8r = v8 * (t3r - t1r);
final double b8i = v8 * (t3i - t1i);
final double u0r = b0r + b1r;
final double u0i = b0i + b1i;
final double u1r = b2r + b3r;
final double u1i = b2i + b3i;
final double u2r = b4r - b3r;
final double u2i = b4i - b3i;
final double u3r = -b2r - b4r;
final double u3i = -b2i - b4i;
final double u4r = b6r + b7r;
final double u4i = b6i + b7i;
final double u5r = b8r - b7r;
final double u5i = b8i - b7i;
final double u6r = -b8r - b6r;
final double u6i = -b8i - b6i;
final double u7r = u0r + u1r;
final double u7i = u0i + u1i;
final double u8r = u0r + u2r;
final double u8i = u0i + u2i;
final double u9r = u0r + u3r;
final double u9i = u0i + u3i;
final double u10r = u4r + b5r;
final double u10i = u4i + b5i;
final double u11r = u5r + b5r;
final double u11i = u5i + b5i;
final double u12r = u6r + b5r;
final double u12i = u6i + b5i;
ret[j] = b0r;
ret[j + im] = b0i; | ||
| File | Project | Line |
|---|---|---|
| ffx/potential/commands/test/ReorderAtoms.java | Potential | 803 |
| ffx/potential/parsers/CIFFilter.java | Potential | 828 |
for (Atom hydrogen : xyzatoms.get(i)) {
if (hydrogen.isHydrogen()) {
Bond bond0 = hydrogen.getBonds().get(0);
Atom atom1 = bond0.get1_2(hydrogen);
double angle0_2 = 0;
List<Angle> anglesList = hydrogen.getAngles(); // Same as bond
Atom atom2 = null;
switch (anglesList.size()) {
case 0 ->
// H-Cl No angles
// Place hydrogen slightly off center of bonded atom (~1Å away).
hydrogen.moveTo(new double[]{atom1.getX() - 0.6, atom1.getY() - 0.6,
atom1.getZ() - 0.6});
case 1 -> {
// H-O=C Need torsion
// H-O-H
for (Angle angle : anglesList) {
atom2 = angle.get1_3(hydrogen);
if (atom2 != null) {
angle0_2 = angle.angleType.angle[0];
}
if (angle0_2 != 0) {
//if angle0_2 is found then no need to look for another atom.
break;
}
}
assert atom2 != null;
List<Bond> bonds2 = atom2.getBonds();
Atom atom3 =
(bonds2.size() > 1 && atom1 == bonds2.get(0).get1_2(atom2)) ? bonds2.get(
1).get1_2(atom2) : bonds2.get(0).get1_2(atom2);
double diAng = dihedralAngle(hydrogen.getXYZ(null), atom1.getXYZ(null),
atom2.getXYZ(null), atom3.getXYZ(null));
if (atom1 != atom3) {
intxyz(hydrogen, atom1, bond0.bondType.distance, atom2, angle0_2, atom3,
Math.toDegrees(diAng), 0);
} else {
// Likely water as Atom 3 is not unique. Since no hydrogen atoms
// are present, there isn't a third atom...
double[] coord = new double[]{atom2.getX(), atom2.getY(), atom3.getZ()};
double mag = sqrt(
coord[0] * coord[0] + coord[1] * coord[1] + coord[2] * coord[2]);
coord[0] /= mag;
coord[1] /= mag;
coord[2] /= mag;
hydrogen.moveTo(atom1.getX() - coord[0], atom1.getY() - coord[1],
atom1.getZ() - coord[2]);
}
}
default -> {
// H-C(-C)(-C)-C
Atom atom2B = null;
double angle0_2B = 0;
Atom proposedAtom;
double proposedAngle = 0;
int chiral = 1;
for (Angle angle : anglesList) {
proposedAtom = angle.get1_3(hydrogen);
if (proposedAtom != null && !proposedAtom.isHydrogen()) {
proposedAngle = angle.angleType.angle[0];
}
if (proposedAngle != 0) {
// If angle1_3 is found then no need to look for another atom.
if (angle0_2 != 0) {
atom2B = proposedAtom;
angle0_2B = proposedAngle;
} else {
atom2 = proposedAtom;
angle0_2 = proposedAngle;
}
proposedAngle = 0.0;
}
if (angle0_2 != 0 && angle0_2B != 0) {
break;
}
}
if (lastKnownAtom1 == null || lastKnownAtom1 != atom1) {
lastKnownAtom1 = atom1;
knownCount = 0;
} else {
knownCount++;
}
if (angle0_2B == 0.0) {
// Hydrogen position depends on other hydrogen, use generic location
chiral = 0;
assert atom2 != null;
List<Bond> bonds2 = atom2.getBonds();
atom2B =
(atom1 == bonds2.get(0).get1_2(atom2)) ? bonds2.get(1).get1_2(atom2)
: bonds2.get(0).get1_2(atom2);
// Evenly space out hydrogen
angle0_2B = 180.0 - 120.0 * knownCount;
} else if (anglesList.size() == 2) {
// Trigonal hydrogen use equipartition between selected atoms
chiral = 3;
} else if (knownCount == 1) {
// Already created one hydrogen (chiral = 1), use other.
chiral = -1;
}
//TODO discern whether to use chiral = 1 or -1 when angles are to three heavy atoms.
intxyz(hydrogen, atom1, bond0.bondType.distance, atom2, angle0_2, atom2B,
angle0_2B, chiral);
}
}
}
}
} else { | ||
| File | Project | Line |
|---|---|---|
| ffx/potential/bonded/RotamerLibrary.java | Potential | 2198 |
| ffx/potential/bonded/RotamerLibrary.java | Potential | 2268 |
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 | 1998 |
| ffx/xray/CrystalReciprocalSpace.java | Refinement | 2233 |
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 | 696 |
| 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 |
|---|---|---|
| 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, | ||