use of ffx.potential.bonded.Atom in project ffx by mjschnie.
the class ForceFieldEnergyOpenMM method addTorsionTorsionForce.
private void addTorsionTorsionForce() {
TorsionTorsion[] torsionTorsions = super.getTorsionTorsions();
if (torsionTorsions == null || torsionTorsions.length < 1) {
return;
}
/**
* Load the torsion-torsions.
*/
int nTypes = 0;
LinkedHashMap<String, TorsionTorsionType> torTorTypes = new LinkedHashMap<>();
int nTorsionTorsions = torsionTorsions.length;
amoebaTorsionTorsionForce = OpenMM_AmoebaTorsionTorsionForce_create();
for (int i = 0; i < nTorsionTorsions; i++) {
TorsionTorsion torsionTorsion = torsionTorsions[i];
int ia = torsionTorsion.getAtom(0).getXyzIndex() - 1;
int ib = torsionTorsion.getAtom(1).getXyzIndex() - 1;
int ic = torsionTorsion.getAtom(2).getXyzIndex() - 1;
int id = torsionTorsion.getAtom(3).getXyzIndex() - 1;
int ie = torsionTorsion.getAtom(4).getXyzIndex() - 1;
TorsionTorsionType torsionTorsionType = torsionTorsion.torsionTorsionType;
String key = torsionTorsionType.getKey();
/**
* Check if the TorTor parameters have already been added to the
* Hash.
*/
int gridIndex = 0;
if (torTorTypes.containsKey(key)) {
/**
* If the TorTor has been added, get its (ordered) index in the
* Hash.
*/
int index = 0;
for (String entry : torTorTypes.keySet()) {
if (entry.equalsIgnoreCase(key)) {
gridIndex = index;
break;
} else {
index++;
}
}
} else {
/**
* Add the new TorTor.
*/
torTorTypes.put(key, torsionTorsionType);
gridIndex = nTypes;
nTypes++;
}
Atom atom = torsionTorsion.getChiralAtom();
int iChiral = -1;
if (atom != null) {
iChiral = atom.getXyzIndex() - 1;
}
OpenMM_AmoebaTorsionTorsionForce_addTorsionTorsion(amoebaTorsionTorsionForce, ia, ib, ic, id, ie, iChiral, gridIndex);
}
/**
* Load the Torsion-Torsion parameters.
*/
PointerByReference values = OpenMM_DoubleArray_create(6);
int gridIndex = 0;
for (String key : torTorTypes.keySet()) {
TorsionTorsionType torTorType = torTorTypes.get(key);
int nx = torTorType.nx;
int ny = torTorType.ny;
double[] tx = torTorType.tx;
double[] ty = torTorType.ty;
double[] f = torTorType.energy;
double[] dx = torTorType.dx;
double[] dy = torTorType.dy;
double[] dxy = torTorType.dxy;
/**
* Create the 3D grid.
*/
PointerByReference grid3D = OpenMM_3D_DoubleArray_create(nx, ny, 6);
int xIndex = 0;
int yIndex = 0;
for (int j = 0; j < nx * ny; j++) {
int addIndex = 0;
OpenMM_DoubleArray_set(values, addIndex++, tx[xIndex]);
OpenMM_DoubleArray_set(values, addIndex++, ty[yIndex]);
OpenMM_DoubleArray_set(values, addIndex++, OpenMM_KJPerKcal * f[j]);
OpenMM_DoubleArray_set(values, addIndex++, OpenMM_KJPerKcal * dx[j]);
OpenMM_DoubleArray_set(values, addIndex++, OpenMM_KJPerKcal * dy[j]);
OpenMM_DoubleArray_set(values, addIndex++, OpenMM_KJPerKcal * dxy[j]);
OpenMM_3D_DoubleArray_set(grid3D, yIndex, xIndex, values);
xIndex++;
if (xIndex == nx) {
xIndex = 0;
yIndex++;
}
}
OpenMM_AmoebaTorsionTorsionForce_setTorsionTorsionGrid(amoebaTorsionTorsionForce, gridIndex++, grid3D);
OpenMM_3D_DoubleArray_destroy(grid3D);
}
OpenMM_DoubleArray_destroy(values);
OpenMM_System_addForce(system, amoebaTorsionTorsionForce);
logger.log(Level.INFO, " Added Torsion-Torsions ({0})", nTorsionTorsions);
}
use of ffx.potential.bonded.Atom in project ffx by mjschnie.
the class ForceFieldEnergyOpenMM method createVirtualHydrogenSites.
/**
* Experimental. Virtual hydrogen sites require creation of new particles,
* which then need to be handled (ignored?) for the multiple force.
*/
private void createVirtualHydrogenSites() {
VanDerWaals vdW = super.getVdwNode();
if (vdW == null) {
return;
}
int[] ired = vdW.getReductionIndex();
Atom[] atoms = molecularAssembly.getAtomArray();
int nAtoms = atoms.length;
for (int i = 0; i < nAtoms; i++) {
Atom atom = atoms[i];
VDWType vdwType = atom.getVDWType();
if (vdwType.reductionFactor < 1.0) {
double factor = vdwType.reductionFactor;
// Create the virtual site.
PointerByReference virtualSite = OpenMM_TwoParticleAverageSite_create(i, ired[i], factor, 1.0 - factor);
// Create a massless particle for the hydrogen vdW site.
int id = OpenMM_System_addParticle(system, 0.0);
// Denote the massless particle is a virtual site
OpenMM_System_setVirtualSite(system, id, virtualSite);
}
}
}
use of ffx.potential.bonded.Atom in project ffx by mjschnie.
the class ForceFieldEnergyOpenMM method addCustomGBForce.
private void addCustomGBForce() {
GeneralizedKirkwood gk = super.getGK();
if (gk == null) {
return;
}
customGBForce = OpenMM_CustomGBForce_create();
OpenMM_CustomGBForce_addPerParticleParameter(customGBForce, "q");
OpenMM_CustomGBForce_addPerParticleParameter(customGBForce, "radius");
OpenMM_CustomGBForce_addPerParticleParameter(customGBForce, "scale");
OpenMM_CustomGBForce_addGlobalParameter(customGBForce, "solventDielectric", 78.3);
OpenMM_CustomGBForce_addGlobalParameter(customGBForce, "soluteDielectric", 1.0);
// Factor of 0.1 for Ang to nm.
OpenMM_CustomGBForce_addGlobalParameter(customGBForce, "dOffset", gk.getDielecOffset() * OpenMM_NmPerAngstrom);
OpenMM_CustomGBForce_addComputedValue(customGBForce, "I", // "step(r+sr2-or1)*0.5*((1/L^3-1/U^3)/3+(1/U^4-1/L^4)/8*(r-sr2*sr2/r)+0.25*(1/U^2-1/L^2)/r+C);"
"0.5*((1/L^3-1/U^3)/3.0+(1/U^4-1/L^4)/8.0*(r-sr2*sr2/r)+0.25*(1/U^2-1/L^2)/r+C);" + "U=r+sr2;" + // + "C=2*(1/or1-1/L)*step(sr2-r-or1);"
"C=2/3*(1/or1^3-1/L^3)*step(sr2-r-or1);" + // + "D=step(r-sr2)*(r-sr2) + (1-step(r-sr2))*(sr2-r);"
"L = step(sr2 - r1r)*sr2mr + (1 - step(sr2 - r1r))*L;" + "sr2mr = sr2 - r;" + "r1r = radius1 + r;" + "L = step(r1sr2 - r)*radius1 + (1 - step(r1sr2 - r))*L;" + "r1sr2 = radius1 + sr2;" + "L = r - sr2;" + "sr2 = scale2 * radius2;" + "or1 = radius1; or2 = radius2", OpenMM_CustomGBForce_ParticlePairNoExclusions);
OpenMM_CustomGBForce_addComputedValue(customGBForce, "B", // "psi=I*or; or=radius-0.009"
"step(BB-radius)*BB + (1 - step(BB-radius))*radius;" + "BB = 1 / ( (3.0*III)^(1.0/3.0) );" + "III = step(II)*II + (1 - step(II))*1.0e-9/3.0;" + "II = maxI - I;" + "maxI = 1/(3.0*radius^3)", OpenMM_CustomGBForce_SingleParticle);
double sTens = gk.getSurfaceTension();
logger.info(String.format(" FFX surface tension: %9.5g kcal/mol/Ang^2", sTens));
sTens *= OpenMM_KJPerKcal;
// 100 square Angstroms per square nanometer.
sTens *= 100.0;
logger.info(String.format(" OpenMM surface tension: %9.5g kJ/mol/nm^2", sTens));
String surfaceTension = Double.toString(sTens);
OpenMM_CustomGBForce_addEnergyTerm(customGBForce, surfaceTension + "*(radius+0.14+dOffset)^2*((radius+dOffset)/B)^6/6-0.5*138.935456*(1/soluteDielectric-1/solventDielectric)*q^2/B", OpenMM_CustomGBForce_SingleParticle);
/**
* Particle pair term is the generalized Born cross term.
*/
OpenMM_CustomGBForce_addEnergyTerm(customGBForce, "-138.935456*(1/soluteDielectric-1/solventDielectric)*q1*q2/f;" + "f=sqrt(r^2+B1*B2*exp(-r^2/(2.455*B1*B2)))", OpenMM_CustomGBForce_ParticlePair);
double[] baseRadii = gk.getBaseRadii();
double[] overlapScale = gk.getOverlapScale();
Atom[] atoms = molecularAssembly.getAtomArray();
int nAtoms = atoms.length;
PointerByReference doubleArray = OpenMM_DoubleArray_create(0);
for (int i = 0; i < nAtoms; i++) {
MultipoleType multipoleType = atoms[i].getMultipoleType();
OpenMM_DoubleArray_append(doubleArray, multipoleType.charge);
OpenMM_DoubleArray_append(doubleArray, OpenMM_NmPerAngstrom * baseRadii[i]);
OpenMM_DoubleArray_append(doubleArray, overlapScale[i]);
OpenMM_CustomGBForce_addParticle(customGBForce, doubleArray);
OpenMM_DoubleArray_resize(doubleArray, 0);
}
OpenMM_DoubleArray_destroy(doubleArray);
double cut = gk.getCutoff();
OpenMM_CustomGBForce_setCutoffDistance(customGBForce, cut);
OpenMM_Force_setForceGroup(customGBForce, 1);
OpenMM_System_addForce(system, customGBForce);
logger.log(Level.INFO, " Added generalized Born force");
}
use of ffx.potential.bonded.Atom in project ffx by mjschnie.
the class MolecularAssembly method getMolecule.
private Atom getMolecule(Atom atom, boolean create) {
String resName = atom.getResidueName();
int resNum = atom.getResidueNumber();
Character chainID = atom.getChainID();
String segID = atom.getSegID();
ArrayList<MSNode> list = ions.getChildList();
for (MSNode node : list) {
Molecule m = (Molecule) node;
if (m.getSegID().equalsIgnoreCase(segID) && m.getResidueName().equalsIgnoreCase(resName) && m.getResidueNumber() == resNum) {
return (Atom) m.addMSNode(atom);
}
}
list = water.getChildList();
for (MSNode node : list) {
Molecule m = (Molecule) node;
if (m.getSegID().equalsIgnoreCase(segID) && m.getResidueName().equalsIgnoreCase(resName) && m.getResidueNumber() == resNum) {
return (Atom) m.addMSNode(atom);
}
}
list = molecules.getChildList();
for (MSNode node : list) {
Molecule m = (Molecule) node;
if (m.getSegID().equalsIgnoreCase(segID) && m.getResidueName().equalsIgnoreCase(resName) && m.getResidueNumber() == resNum) {
return (Atom) m.addMSNode(atom);
}
}
if (create) {
Molecule m = new Molecule(resName, resNum, chainID, segID);
m.addMSNode(atom);
if (resName.equalsIgnoreCase("DOD") || resName.equalsIgnoreCase("HOH") || resName.equalsIgnoreCase("WAT")) {
water.add(m);
// NA, K, MG, MG2, CA, CA2, CL
} else if (resName.equalsIgnoreCase("NA") || resName.equalsIgnoreCase("K") || resName.equalsIgnoreCase("MG") || resName.equalsIgnoreCase("MG2") || resName.equalsIgnoreCase("CA") || resName.equalsIgnoreCase("CA2") || resName.equalsIgnoreCase("CL") || resName.equalsIgnoreCase("BR") || resName.equalsIgnoreCase("ZN") || resName.equalsIgnoreCase("ZN2")) {
ions.add(m);
} else {
molecules.add(m);
}
return atom;
} else {
return null;
}
}
use of ffx.potential.bonded.Atom in project ffx by mjschnie.
the class MolecularAssembly method findAtom.
/**
* <p>
* findAtom</p>
*
* @param atom a {@link ffx.potential.bonded.Atom} object.
* @return a {@link ffx.potential.bonded.Atom} object.
*/
public Atom findAtom(Atom atom) {
if (!atom.isHetero() || atom.isModRes()) {
Polymer polymer = getPolymer(atom.getChainID(), atom.getSegID(), false);
if (polymer != null) {
Residue res = polymer.getResidue(atom.getResidueName(), atom.getResidueNumber(), false);
if (res != null) {
MSNode node = res.getAtomNode();
Atom root = (Atom) node.contains(atom);
return root;
}
}
return null;
} else {
ArrayList<Molecule> list = getMolecules();
for (MSNode node : list) {
Molecule m = (Molecule) node;
if (m.getSegID().equalsIgnoreCase(atom.getSegID()) && m.getResidueName().equalsIgnoreCase(atom.getResidueName()) && m.getResidueNumber() == atom.getResidueNumber()) {
Atom root = (Atom) node.contains(atom);
return root;
}
}
ArrayList<MSNode> ionList = getIons();
for (MSNode node : ionList) {
Molecule m = (Molecule) node;
if (m.getSegID().equalsIgnoreCase(atom.getSegID()) && m.getResidueName().equalsIgnoreCase(atom.getResidueName()) && m.getResidueNumber() == atom.getResidueNumber()) {
Atom root = (Atom) node.contains(atom);
return root;
}
}
ArrayList<MSNode> waterList = getWaters();
for (MSNode node : waterList) {
Molecule m = (Molecule) node;
if (m.getSegID().equalsIgnoreCase(atom.getSegID()) && m.getResidueName().equalsIgnoreCase(atom.getResidueName()) && m.getResidueNumber() == atom.getResidueNumber()) {
Atom root = (Atom) node.contains(atom);
return root;
}
}
return null;
}
}
Aggregations