use of ffx.potential.parameters.AtomType in project ffx by mjschnie.
the class ParticleMeshEwaldCart method assignMultipole.
private boolean assignMultipole(int i) {
Atom atom = atoms[i];
AtomType atomType = atoms[i].getAtomType();
if (atomType == null) {
String message = " Multipoles can only be assigned to atoms that have been typed.";
logger.severe(message);
return false;
}
PolarizeType polarizeType = forceField.getPolarizeType(atomType.getKey());
if (polarizeType != null) {
atom.setPolarizeType(polarizeType);
} else {
String message = " No polarization type was found for " + atom.toString();
logger.fine(message);
double polarizability = 0.0;
double thole = 0.0;
int[] polarizationGroup = null;
polarizeType = new PolarizeType(atomType.type, polarizability, thole, polarizationGroup);
forceField.addForceFieldType(polarizeType);
atom.setPolarizeType(polarizeType);
}
String key;
// No reference atoms.
key = atomType.getKey() + " 0 0";
MultipoleType multipoleType = forceField.getMultipoleType(key);
if (multipoleType != null) {
atom.setMultipoleType(multipoleType);
localMultipole[i][t000] = multipoleType.getCharge();
localMultipole[i][t100] = multipoleType.getDipole()[0];
localMultipole[i][t010] = multipoleType.getDipole()[1];
localMultipole[i][t001] = multipoleType.getDipole()[2];
localMultipole[i][t200] = multipoleType.getQuadrupole()[0][0];
localMultipole[i][t020] = multipoleType.getQuadrupole()[1][1];
localMultipole[i][t002] = multipoleType.getQuadrupole()[2][2];
localMultipole[i][t110] = multipoleType.getQuadrupole()[0][1];
localMultipole[i][t101] = multipoleType.getQuadrupole()[0][2];
localMultipole[i][t011] = multipoleType.getQuadrupole()[1][2];
axisAtom[i] = null;
frame[i] = multipoleType.frameDefinition;
return true;
}
// No bonds.
List<Bond> bonds = atom.getBonds();
if (bonds == null || bonds.size() < 1) {
String message = "Multipoles can only be assigned after bonded relationships are defined.\n";
logger.severe(message);
}
// 1 reference atom.
for (Bond b : bonds) {
Atom atom2 = b.get1_2(atom);
key = atomType.getKey() + " " + atom2.getAtomType().getKey() + " 0";
multipoleType = multipoleType = forceField.getMultipoleType(key);
if (multipoleType != null) {
int[] multipoleReferenceAtoms = new int[1];
multipoleReferenceAtoms[0] = atom2.getIndex() - 1;
atom.setMultipoleType(multipoleType);
localMultipole[i][t000] = multipoleType.getCharge();
localMultipole[i][t100] = multipoleType.getDipole()[0];
localMultipole[i][t010] = multipoleType.getDipole()[1];
localMultipole[i][t001] = multipoleType.getDipole()[2];
localMultipole[i][t200] = multipoleType.getQuadrupole()[0][0];
localMultipole[i][t020] = multipoleType.getQuadrupole()[1][1];
localMultipole[i][t002] = multipoleType.getQuadrupole()[2][2];
localMultipole[i][t110] = multipoleType.getQuadrupole()[0][1];
localMultipole[i][t101] = multipoleType.getQuadrupole()[0][2];
localMultipole[i][t011] = multipoleType.getQuadrupole()[1][2];
axisAtom[i] = multipoleReferenceAtoms;
frame[i] = multipoleType.frameDefinition;
return true;
}
}
// 2 reference atoms.
for (Bond b : bonds) {
Atom atom2 = b.get1_2(atom);
String key2 = atom2.getAtomType().getKey();
for (Bond b2 : bonds) {
if (b == b2) {
continue;
}
Atom atom3 = b2.get1_2(atom);
String key3 = atom3.getAtomType().getKey();
key = atomType.getKey() + " " + key2 + " " + key3;
multipoleType = forceField.getMultipoleType(key);
if (multipoleType != null) {
int[] multipoleReferenceAtoms = new int[2];
multipoleReferenceAtoms[0] = atom2.getIndex() - 1;
multipoleReferenceAtoms[1] = atom3.getIndex() - 1;
atom.setMultipoleType(multipoleType);
atom.setMultipoleType(multipoleType);
localMultipole[i][t000] = multipoleType.getCharge();
localMultipole[i][t100] = multipoleType.getDipole()[0];
localMultipole[i][t010] = multipoleType.getDipole()[1];
localMultipole[i][t001] = multipoleType.getDipole()[2];
localMultipole[i][t200] = multipoleType.getQuadrupole()[0][0];
localMultipole[i][t020] = multipoleType.getQuadrupole()[1][1];
localMultipole[i][t002] = multipoleType.getQuadrupole()[2][2];
localMultipole[i][t110] = multipoleType.getQuadrupole()[0][1];
localMultipole[i][t101] = multipoleType.getQuadrupole()[0][2];
localMultipole[i][t011] = multipoleType.getQuadrupole()[1][2];
axisAtom[i] = multipoleReferenceAtoms;
frame[i] = multipoleType.frameDefinition;
return true;
}
}
}
/**
* 3 reference atoms.
*/
for (Bond b : bonds) {
Atom atom2 = b.get1_2(atom);
String key2 = atom2.getAtomType().getKey();
for (Bond b2 : bonds) {
if (b == b2) {
continue;
}
Atom atom3 = b2.get1_2(atom);
String key3 = atom3.getAtomType().getKey();
for (Bond b3 : bonds) {
if (b == b3 || b2 == b3) {
continue;
}
Atom atom4 = b3.get1_2(atom);
String key4 = atom4.getAtomType().getKey();
key = atomType.getKey() + " " + key2 + " " + key3 + " " + key4;
multipoleType = forceField.getMultipoleType(key);
if (multipoleType != null) {
int[] multipoleReferenceAtoms = new int[3];
multipoleReferenceAtoms[0] = atom2.getIndex() - 1;
multipoleReferenceAtoms[1] = atom3.getIndex() - 1;
multipoleReferenceAtoms[2] = atom4.getIndex() - 1;
atom.setMultipoleType(multipoleType);
localMultipole[i][t000] = multipoleType.getCharge();
localMultipole[i][t100] = multipoleType.getDipole()[0];
localMultipole[i][t010] = multipoleType.getDipole()[1];
localMultipole[i][t001] = multipoleType.getDipole()[2];
localMultipole[i][t200] = multipoleType.getQuadrupole()[0][0];
localMultipole[i][t020] = multipoleType.getQuadrupole()[1][1];
localMultipole[i][t002] = multipoleType.getQuadrupole()[2][2];
localMultipole[i][t110] = multipoleType.getQuadrupole()[0][1];
localMultipole[i][t101] = multipoleType.getQuadrupole()[0][2];
localMultipole[i][t011] = multipoleType.getQuadrupole()[1][2];
axisAtom[i] = multipoleReferenceAtoms;
frame[i] = multipoleType.frameDefinition;
return true;
}
}
List<Angle> angles = atom.getAngles();
for (Angle angle : angles) {
Atom atom4 = angle.get1_3(atom);
if (atom4 != null) {
String key4 = atom4.getAtomType().getKey();
key = atomType.getKey() + " " + key2 + " " + key3 + " " + key4;
multipoleType = forceField.getMultipoleType(key);
if (multipoleType != null) {
int[] multipoleReferenceAtoms = new int[3];
multipoleReferenceAtoms[0] = atom2.getIndex() - 1;
multipoleReferenceAtoms[1] = atom3.getIndex() - 1;
multipoleReferenceAtoms[2] = atom4.getIndex() - 1;
atom.setMultipoleType(multipoleType);
localMultipole[i][t000] = multipoleType.getCharge();
localMultipole[i][t100] = multipoleType.getDipole()[0];
localMultipole[i][t010] = multipoleType.getDipole()[1];
localMultipole[i][t001] = multipoleType.getDipole()[2];
localMultipole[i][t200] = multipoleType.getQuadrupole()[0][0];
localMultipole[i][t020] = multipoleType.getQuadrupole()[1][1];
localMultipole[i][t002] = multipoleType.getQuadrupole()[2][2];
localMultipole[i][t110] = multipoleType.getQuadrupole()[0][1];
localMultipole[i][t101] = multipoleType.getQuadrupole()[0][2];
localMultipole[i][t011] = multipoleType.getQuadrupole()[1][2];
axisAtom[i] = multipoleReferenceAtoms;
frame[i] = multipoleType.frameDefinition;
return true;
}
}
}
}
}
/**
* Revert to a 2 reference atom definition that may include a 1-3 site.
* For example a hydrogen on water.
*/
for (Bond b : bonds) {
Atom atom2 = b.get1_2(atom);
String key2 = atom2.getAtomType().getKey();
List<Angle> angles = atom.getAngles();
for (Angle angle : angles) {
Atom atom3 = angle.get1_3(atom);
if (atom3 != null) {
String key3 = atom3.getAtomType().getKey();
key = atomType.getKey() + " " + key2 + " " + key3;
multipoleType = forceField.getMultipoleType(key);
if (multipoleType != null) {
int[] multipoleReferenceAtoms = new int[2];
multipoleReferenceAtoms[0] = atom2.getIndex() - 1;
multipoleReferenceAtoms[1] = atom3.getIndex() - 1;
atom.setMultipoleType(multipoleType);
localMultipole[i][t000] = multipoleType.getCharge();
localMultipole[i][t100] = multipoleType.getDipole()[0];
localMultipole[i][t010] = multipoleType.getDipole()[1];
localMultipole[i][t001] = multipoleType.getDipole()[2];
localMultipole[i][t200] = multipoleType.getQuadrupole()[0][0];
localMultipole[i][t020] = multipoleType.getQuadrupole()[1][1];
localMultipole[i][t002] = multipoleType.getQuadrupole()[2][2];
localMultipole[i][t110] = multipoleType.getQuadrupole()[0][1];
localMultipole[i][t101] = multipoleType.getQuadrupole()[0][2];
localMultipole[i][t011] = multipoleType.getQuadrupole()[1][2];
axisAtom[i] = multipoleReferenceAtoms;
frame[i] = multipoleType.frameDefinition;
return true;
}
for (Angle angle2 : angles) {
Atom atom4 = angle2.get1_3(atom);
if (atom4 != null && atom4 != atom3) {
String key4 = atom4.getAtomType().getKey();
key = atomType.getKey() + " " + key2 + " " + key3 + " " + key4;
multipoleType = forceField.getMultipoleType(key);
if (multipoleType != null) {
int[] multipoleReferenceAtoms = new int[3];
multipoleReferenceAtoms[0] = atom2.getIndex() - 1;
multipoleReferenceAtoms[1] = atom3.getIndex() - 1;
multipoleReferenceAtoms[2] = atom4.getIndex() - 1;
atom.setMultipoleType(multipoleType);
localMultipole[i][t000] = multipoleType.getCharge();
localMultipole[i][t100] = multipoleType.getDipole()[0];
localMultipole[i][t010] = multipoleType.getDipole()[1];
localMultipole[i][t001] = multipoleType.getDipole()[2];
localMultipole[i][t200] = multipoleType.getQuadrupole()[0][0];
localMultipole[i][t020] = multipoleType.getQuadrupole()[1][1];
localMultipole[i][t002] = multipoleType.getQuadrupole()[2][2];
localMultipole[i][t110] = multipoleType.getQuadrupole()[0][1];
localMultipole[i][t101] = multipoleType.getQuadrupole()[0][2];
localMultipole[i][t011] = multipoleType.getQuadrupole()[1][2];
axisAtom[i] = multipoleReferenceAtoms;
frame[i] = multipoleType.frameDefinition;
return true;
}
}
}
}
}
}
return false;
}
use of ffx.potential.parameters.AtomType in project ffx by mjschnie.
the class NucleicAcidUtils method assignNucleicAcidAtomTypes.
/**
* Assign atom types for a nucleic acid polymer.
*
* @param residues A list of residues that form the nucleic acid polymer.
* @param forceField The ForceField in use.
* @param bondList A list of created bonds.
* @throws MissingHeavyAtomException
* @throws MissingAtomTypeException
*/
public static void assignNucleicAcidAtomTypes(List<Residue> residues, ForceField forceField, ArrayList<Bond> bondList) throws MissingHeavyAtomException, MissingAtomTypeException {
/**
* A reference to the O3* atom of the previous base.
*/
Atom pSugarO3 = null;
/**
* Loop over residues.
*/
int numberOfResidues = residues.size();
for (int residueNumber = 0; residueNumber < numberOfResidues; residueNumber++) {
/**
* Match the residue name to a known nucleic acid residue.
*/
Residue residue = residues.get(residueNumber);
String residueName = residue.getName().toUpperCase();
NucleicAcid3 nucleicAcid = null;
int naNumber = -1;
for (NucleicAcid3 nucleic : nucleicAcidList) {
naNumber++;
String nuc3 = nucleic.toString();
nuc3 = nuc3.substring(nuc3.length() - 3);
if (nuc3.equalsIgnoreCase(residueName)) {
nucleicAcid = nucleic;
break;
}
}
/**
* Do atom name conversions.
*/
List<Atom> resAtoms = residue.getAtomList();
int natoms = resAtoms.size();
for (int i = 0; i < natoms; i++) {
Atom atom = resAtoms.get(i);
String name = atom.getName();
name = name.replace('*', '\'');
// name = name.replace('D', 'H');
atom.setName(name);
}
/**
* Check if the sugar is deoxyribose and change the residue name if
* necessary.
*/
boolean isDNA = false;
Atom sugarO2 = (Atom) residue.getAtomNode("O2\'");
if (sugarO2 == null) {
/**
* Assume deoxyribose (DNA) since there is an O2* atom.
*/
isDNA = true;
if (!residueName.startsWith("D")) {
switch(nucleicAcid) {
case ADE:
nucleicAcid = NucleicAcid3.DAD;
residueName = "DAD";
residue.setName(residueName);
break;
case CYT:
nucleicAcid = NucleicAcid3.DCY;
residueName = "DCY";
residue.setName(residueName);
break;
case GUA:
nucleicAcid = NucleicAcid3.DGU;
residueName = "DGU";
residue.setName(residueName);
break;
case THY:
nucleicAcid = NucleicAcid3.DTY;
residueName = "DTY";
residue.setName(residueName);
break;
default:
}
}
} else /**
* Assume ribose (RNA) since there is an O2* atom.
*/
{
if (residueName.startsWith("D")) {
switch(nucleicAcid) {
case DAD:
nucleicAcid = NucleicAcid3.ADE;
residueName = "ADE";
residue.setName(residueName);
break;
case DCY:
nucleicAcid = NucleicAcid3.CYT;
residueName = "CYT";
residue.setName(residueName);
break;
case DGU:
nucleicAcid = NucleicAcid3.GUA;
residueName = "GUA";
residue.setName(residueName);
break;
case DTY:
nucleicAcid = NucleicAcid3.THY;
residueName = "THY";
residue.setName(residueName);
break;
default:
}
}
}
/**
* Set a position flag.
*/
Residue.ResiduePosition position = MIDDLE_RESIDUE;
if (residueNumber == 0) {
position = FIRST_RESIDUE;
} else if (residueNumber == numberOfResidues - 1) {
position = LAST_RESIDUE;
}
/**
* Build the phosphate atoms of the current residue.
*/
Atom phosphate = null;
Atom sugarO5 = null;
if (position == FIRST_RESIDUE) {
/**
* The 5' O5' oxygen of the nucleic acid is generally terminated
* by 1.) A phosphate group PO3 (-3). 2.) A hydrogen.
*
* If the base has phosphate atom we will assume a PO3 group.
*/
phosphate = (Atom) residue.getAtomNode("P");
if (phosphate != null) {
if (isDNA) {
phosphate = buildHeavy(residue, "P", null, 1247, forceField, bondList);
buildHeavy(residue, "OP1", phosphate, 1248, forceField, bondList);
buildHeavy(residue, "OP2", phosphate, 1248, forceField, bondList);
buildHeavy(residue, "OP3", phosphate, 1248, forceField, bondList);
sugarO5 = buildHeavy(residue, "O5\'", phosphate, 1246, forceField, bondList);
} else {
phosphate = buildHeavy(residue, "P", null, 1235, forceField, bondList);
buildHeavy(residue, "OP1", phosphate, 1236, forceField, bondList);
buildHeavy(residue, "OP2", phosphate, 1236, forceField, bondList);
buildHeavy(residue, "OP3", phosphate, 1236, forceField, bondList);
sugarO5 = buildHeavy(residue, "O5\'", phosphate, 1234, forceField, bondList);
}
} else if (isDNA) {
sugarO5 = buildHeavy(residue, "O5\'", phosphate, 1244, forceField, bondList);
} else {
sugarO5 = buildHeavy(residue, "O5\'", phosphate, 1232, forceField, bondList);
}
} else {
phosphate = buildHeavy(residue, "P", pSugarO3, NA_P[naNumber], forceField, bondList);
buildHeavy(residue, "OP1", phosphate, NA_OP[naNumber], forceField, bondList);
buildHeavy(residue, "OP2", phosphate, NA_OP[naNumber], forceField, bondList);
sugarO5 = buildHeavy(residue, "O5\'", phosphate, NA_O5[naNumber], forceField, bondList);
}
/**
* Build the ribose sugar atoms of the current base.
*/
Atom sugarC5 = buildHeavy(residue, "C5\'", sugarO5, NA_C5[naNumber], forceField, bondList);
Atom sugarC4 = buildHeavy(residue, "C4\'", sugarC5, NA_C4[naNumber], forceField, bondList);
Atom sugarO4 = buildHeavy(residue, "O4\'", sugarC4, NA_O4[naNumber], forceField, bondList);
Atom sugarC1 = buildHeavy(residue, "C1\'", sugarO4, NA_C1[naNumber], forceField, bondList);
Atom sugarC3 = buildHeavy(residue, "C3\'", sugarC4, NA_C3[naNumber], forceField, bondList);
Atom sugarC2 = buildHeavy(residue, "C2\'", sugarC3, NA_C2[naNumber], forceField, bondList);
buildBond(sugarC2, sugarC1, forceField, bondList);
Atom sugarO3 = null;
if (position == LAST_RESIDUE || numberOfResidues == 1) {
if (isDNA) {
sugarO3 = buildHeavy(residue, "O3\'", sugarC3, 1249, forceField, bondList);
} else {
sugarO3 = buildHeavy(residue, "O3\'", sugarC3, 1237, forceField, bondList);
}
} else {
sugarO3 = buildHeavy(residue, "O3\'", sugarC3, NA_O3[naNumber], forceField, bondList);
}
if (!isDNA) {
sugarO2 = buildHeavy(residue, "O2\'", sugarC2, NA_O2[naNumber], forceField, bondList);
}
/**
* Build the backbone hydrogen atoms.
*/
if (position == FIRST_RESIDUE && NA_P == null) {
buildHydrogen(residue, "H5T", sugarO5, 1.00e0, sugarC5, 109.5e0, sugarC4, 180.0e0, 0, NA_H5T[naNumber], forceField, bondList);
}
buildHydrogen(residue, "H5\'1", sugarC5, 1.09e0, sugarO5, 109.5e0, sugarC4, 109.5e0, 1, NA_H51[naNumber], forceField, bondList);
buildHydrogen(residue, "H5\'2", sugarC5, 1.09e0, sugarO5, 109.5e0, sugarC4, 109.5e0, -1, NA_H52[naNumber], forceField, bondList);
buildHydrogen(residue, "H4\'", sugarC4, 1.09e0, sugarC5, 109.5e0, sugarC3, 109.5e0, -1, NA_H4[naNumber], forceField, bondList);
buildHydrogen(residue, "H3\'", sugarC3, 1.09e0, sugarC4, 109.5e0, sugarC2, 109.5e0, -1, NA_H3[naNumber], forceField, bondList);
if (isDNA) {
buildHydrogen(residue, "H2\'1", sugarC2, 1.09e0, sugarC3, 109.5e0, sugarC1, 109.5e0, -1, NA_H21[naNumber], forceField, bondList);
buildHydrogen(residue, "H2\'2", sugarC2, 1.09e0, sugarC3, 109.5e0, sugarC1, 109.5e0, 1, NA_H22[naNumber], forceField, bondList);
} else {
buildHydrogen(residue, "H2\'", sugarC2, 1.09e0, sugarC3, 109.5e0, sugarC1, 109.5e0, -1, NA_H21[naNumber], forceField, bondList);
// Add the NA_O2' Methyl for OMC and OMG
if (nucleicAcid == NucleicAcid3.OMC || nucleicAcid == NucleicAcid3.OMG) {
Atom CM2 = buildHeavy(residue, "CM2", sugarO2, 1427, forceField, bondList);
Atom HM21 = buildHydrogen(residue, "HM21", CM2, 1.08e0, sugarO2, 109.5e0, sugarC2, 0.0e0, 0, 1428, forceField, bondList);
buildHydrogen(residue, "HM22", CM2, 1.08e0, sugarO2, 109.5e0, HM21, 109.5e0, 1, 1429, forceField, bondList);
buildHydrogen(residue, "HM23", CM2, 1.08e0, sugarO2, 109.5e0, HM21, 109.5e0, -1, 1430, forceField, bondList);
} else {
buildHydrogen(residue, "HO\'", sugarO2, 1.00e0, sugarC2, 109.5e0, sugarC3, 180.0e0, 0, NA_H22[naNumber], forceField, bondList);
}
}
buildHydrogen(residue, "H1\'", sugarC1, 1.09e0, sugarO4, 109.5e0, sugarC2, 109.5e0, -1, NA_H1[naNumber], forceField, bondList);
if (position == LAST_RESIDUE || numberOfResidues == 1) {
buildHydrogen(residue, "H3T", sugarO3, 1.00e0, sugarC3, 109.5e0, sugarC4, 180.0e0, 0, NA_H3T[naNumber], forceField, bondList);
// Else, if it is terminated by a 3' phosphate cap:
// Will need to see how PDB would label a 3' phosphate cap.
}
/**
* Build the nucleic acid base.
*/
try {
assignNucleicAcidBaseAtomTypes(nucleicAcid, residue, sugarC1, sugarO4, sugarC2, forceField, bondList);
} catch (MissingHeavyAtomException missingHeavyAtomException) {
throw missingHeavyAtomException;
}
/**
* Do some checks on the current base to make sure all atoms have
* been assigned an atom type.
*/
resAtoms = residue.getAtomList();
for (Atom atom : resAtoms) {
AtomType atomType = atom.getAtomType();
if (atomType == null) {
MissingAtomTypeException missingAtomTypeException = new MissingAtomTypeException(residue, atom);
throw missingAtomTypeException;
}
int numberOfBonds = atom.getNumBonds();
if (numberOfBonds != atomType.valence) {
if (atom == sugarO3 && numberOfBonds == atomType.valence - 1 && position != LAST_RESIDUE && numberOfResidues != 1) {
continue;
}
logger.log(Level.WARNING, format(" An atom for residue %s has the wrong number of bonds:\n %s", residueName, atom.toString()));
logger.log(Level.WARNING, format(" Expected: %d Actual: %d.", atomType.valence, numberOfBonds));
}
}
/**
* Save a reference to the current O3* oxygen.
*/
pSugarO3 = sugarO3;
}
}
use of ffx.potential.parameters.AtomType in project ffx by mjschnie.
the class GeneralizedKirkwood method initAtomArrays.
private void initAtomArrays() {
if (fixedRadii) {
fixedRadii = false;
}
x = particleMeshEwald.coordinates[0][0];
y = particleMeshEwald.coordinates[0][1];
z = particleMeshEwald.coordinates[0][2];
globalMultipole = particleMeshEwald.globalMultipole[0];
inducedDipole = particleMeshEwald.inducedDipole[0];
inducedDipoleCR = particleMeshEwald.inducedDipoleCR[0];
neighborLists = particleMeshEwald.neighborLists;
grad = particleMeshEwald.getGradient();
torque = particleMeshEwald.getTorque();
lambdaGrad = particleMeshEwald.getLambdaGradient();
lambdaTorque = particleMeshEwald.getLambdaTorque();
if (sharedGKField[0] == null || sharedGKField[0].length() < nAtoms) {
sharedGKField[0] = new SharedDoubleArray(nAtoms);
sharedGKField[1] = new SharedDoubleArray(nAtoms);
sharedGKField[2] = new SharedDoubleArray(nAtoms);
sharedGKFieldCR[0] = new SharedDoubleArray(nAtoms);
sharedGKFieldCR[1] = new SharedDoubleArray(nAtoms);
sharedGKFieldCR[2] = new SharedDoubleArray(nAtoms);
sharedBornGrad = new SharedDoubleArray(nAtoms);
baseRadius = new double[nAtoms];
baseRadiusWithBondi = new double[nAtoms];
overlapScale = new double[nAtoms];
rDisp = new double[nAtoms];
born = new double[nAtoms];
use = new boolean[nAtoms];
}
fill(use, true);
for (int i = 0; i < nAtoms; i++) {
baseRadius[i] = 2.0;
// Original value based on small molecule parameterization.
overlapScale[i] = DEFAULT_OVERLAP_SCALE;
// overlapScale[i] = 0.60; // New default value based on 2016 amino acid GK parameterization.
if (useFittedRadii) {
overlapScale[i] = solventRadii.getOverlapScale();
}
int atomicNumber = atoms[i].getAtomicNumber();
AtomType atomType = atoms[i].getAtomType();
switch(atomicNumber) {
case 0:
baseRadius[i] = 0.0;
break;
case 1:
baseRadius[i] = 1.2;
overlapScale[i] = hydrogenOverlapScale;
break;
case 2:
baseRadius[i] = 1.4;
break;
case 5:
baseRadius[i] = 1.8;
break;
case 6:
baseRadius[i] = 1.7;
break;
case 7:
baseRadius[i] = 1.55;
break;
case 8:
baseRadius[i] = 1.52;
break;
case 9:
baseRadius[i] = 1.47;
break;
case 10:
baseRadius[i] = 1.54;
break;
case 14:
baseRadius[i] = 2.1;
break;
case 15:
baseRadius[i] = 1.8;
break;
case 16:
baseRadius[i] = 1.8;
break;
case 17:
baseRadius[i] = 1.75;
break;
case 18:
baseRadius[i] = 1.88;
break;
case 34:
baseRadius[i] = 1.9;
break;
case 35:
baseRadius[i] = 1.85;
break;
case 36:
baseRadius[i] = 2.02;
break;
case 53:
baseRadius[i] = 1.98;
break;
case 54:
baseRadius[i] = 2.16;
break;
default:
baseRadius[i] = 2.00;
}
double bondiFactor = bondiScale;
int atomNumber = atoms[i].getIndex() + 1;
if (useFittedRadii) {
// First check to see if this atom is in the hardcoded maps.
switch(radiiMapType) {
default:
case ATOMTYPE:
// Check for hard-coded AtomType bondi factor.
if (solventRadii.getAtomBondiMap().containsKey(atomType.type)) {
bondiFactor = solventRadii.getAtomBondiMap().get(atomType.type);
if (verboseRadii) {
logger.info(String.format(" (GK) TypeToBondi: Atom %3s-%-4s (%d) with AtomType %d to %.2f (bondi factor %.4f)", atoms[i].getResidueName(), atoms[i].getName(), atomNumber, atomType.type, baseRadius[i] * bondiFactor, bondiFactor));
}
}
// Many aminos use types 8,12 for their CA,HA so these can't be specified in the map.
if ((atoms[i].getAtomType().type == 8 || atoms[i].getAtomType().type == 12)) {
if (atoms[i].getResidueName().equals("ALA")) {
bondiFactor = 1.60;
if (verboseRadii) {
logger.info(String.format(" (GK) TypeToBondi: Atom %3s-%-4s (%d) with AtomType %d to %.2f (bondi factor %.4f)", atoms[i].getResidueName(), atoms[i].getName(), atomNumber, atomType.type, baseRadius[i] * bondiFactor, bondiFactor));
}
}
}
// Currently, we only want S+HS on CYS. So the CB,HB (types 43,44) are disabled in the table.
if ((atoms[i].getAtomType().type == 43 || atoms[i].getAtomType().type == 44)) {
if (atoms[i].getResidueName().equals("CYD")) {
bondiFactor = 1.02;
if (verboseRadii) {
logger.info(String.format(" (GK) TypeToBondi: Atom %3s-%-4s (%d) with AtomType %d to %.2f (bondi factor %.4f)", atoms[i].getResidueName(), atoms[i].getName(), atomNumber, atomType.type, baseRadius[i] * bondiFactor, bondiFactor));
}
}
}
break;
case BIOTYPE:
Map<String, BioType> bioTypes = forceField.getBioTypeMap();
BioType bioType = null;
for (BioType one : bioTypes.values()) {
if (one.atomType == atomType.type) {
bioType = one;
break;
}
}
if (bioType == null) {
logger.severe(String.format("Couldn't find biotype for %s", atomType, toString()));
}
// Check for hard-coded BioType bondi factor.
if (solventRadii.getBioBondiMap().containsKey(bioType.index)) {
double factor = solventRadii.getBioBondiMap().get(bioType.index);
bondiFactor = factor;
if (verboseRadii) {
logger.info(String.format(" (GK) BiotypeToBondi: Atom %3s-%-4s (%d) with BioType %d to %.2f (bondi factor %.4f)", atoms[i].getResidueName(), atoms[i].getName(), atomNumber, bioType.index, baseRadius[i] * bondiFactor, bondiFactor));
}
}
// Many aminos use types 8,12 for their CA,HA so these can't be specified in the map.
if (bioType.index == 8 || bioType.index == 12) {
if (atoms[i].getResidueName().equals("ALA")) {
bondiFactor = 1.60;
if (verboseRadii) {
logger.info(String.format(" (GK) BiotypeToBondi: Atom %3s-%-4s (%d) with BioType %d to %.2f (bondi factor %.4f)", atoms[i].getResidueName(), atoms[i].getName(), atomNumber, bioType.index, baseRadius[i] * bondiFactor, bondiFactor));
}
}
}
// Currently, we only want S+HS on CYS. So the CB,HB (types 43,44) are disabled in the table.
if (bioType.index == 83 || bioType.index == 84) {
if (atoms[i].getResidueName().equals("CYD")) {
bondiFactor = 1.02;
if (verboseRadii) {
logger.info(String.format(" (GK) BiotypeToBondi: Atom %3s-%-4s (%d) with BioType %d to %.2f (bondi factor %.4f)", atoms[i].getResidueName(), atoms[i].getName(), atomNumber, bioType.index, baseRadius[i] * bondiFactor, bondiFactor));
}
}
}
break;
}
}
// Next, check if this Atom has an ISolvRad forcefield parameter.
// if (atoms[i].getISolvRadType() != null) {
// bondiFactor = atoms[i].getISolvRadType().radiusScale;
// logger.info(String.format(" (GK) ISolvRad parameter for Atom %3s-%-4s with AtomType %d to %.2f (bondi factor %.2f)",
// atoms[i].getResidueName(), atoms[i].getName(), atomType.type, baseRadius[i] * bondiFactor, bondiFactor));
// }
ISolvRadType iSolvRadType = forceField.getISolvRadType(Integer.toString(atomType.type));
if (iSolvRadType != null) {
bondiFactor = iSolvRadType.radiusScale;
if (verboseRadii) {
logger.info(String.format(" (GK) ISolvRad parameter for Atom %3s-%-4s with AtomType %d to %.2f (bondi factor %.2f)", atoms[i].getResidueName(), atoms[i].getName(), atomType.type, baseRadius[i] * bondiFactor, bondiFactor));
}
}
// Finally, check for command-line bondi factor override.
if (radiiOverride.containsKey(atomType.type) && !radiiByNumberMap.containsKey(atomNumber)) {
bondiFactor = radiiOverride.get(atomType.type);
logger.info(String.format(" (GK) Scaling Atom %3s-%-4s with AtomType %d to %.2f (bondi factor %.2f)", atoms[i].getResidueName(), atoms[i].getName(), atomType.type, baseRadius[i] * bondiFactor, bondiFactor));
}
if (radiiByNumberMap.containsKey(atomNumber)) {
bondiFactor = radiiByNumberMap.get(atomNumber);
logger.info(String.format(" (GK) Scaling Atom number %d, %3s-%-4s, with factor %.2f", atomNumber, atoms[i].getResidueName(), atoms[i].getName(), bondiFactor));
}
if (!radiiOverride.containsKey(atomType.type)) {
bondiFactor *= globalRadiiScale;
} else {
logger.fine(String.format(" Not scaling atom type %d", atomType.type));
}
baseRadiusWithBondi[i] = baseRadius[i] * bondiFactor;
}
if (dispersionRegion != null) {
dispersionRegion.init();
}
if (cavitationRegion != null) {
cavitationRegion.init();
}
}
use of ffx.potential.parameters.AtomType in project ffx by mjschnie.
the class PDBFilter method assignAtomTypes.
/**
* Assign force field atoms types to common chemistries using "biotype"
* records.
*/
private void assignAtomTypes() {
/**
* Create a list to store bonds defined by PDB atom names.
*/
bondList = new ArrayList<>();
/**
* To Do: Look for cyclic peptides and disulfides.
*/
Polymer[] polymers = activeMolecularAssembly.getChains();
/**
* Loop over chains.
*/
if (polymers != null) {
logger.info(format("\n Assigning atom types for %d chains.", polymers.length));
for (Polymer polymer : polymers) {
List<Residue> residues = polymer.getResidues();
int numberOfResidues = residues.size();
/**
* Check if all residues are known amino acids.
*/
boolean isProtein = true;
for (int residueNumber = 0; residueNumber < numberOfResidues; residueNumber++) {
Residue residue = residues.get(residueNumber);
String name = residue.getName().toUpperCase();
boolean aa = false;
for (AminoAcid3 amino : aminoAcidList) {
if (amino.toString().equalsIgnoreCase(name)) {
aa = true;
checkHydrogenAtomNames(residue);
break;
}
}
// Check for a patch.
if (!aa) {
logger.info(" Checking for non-standard amino acid patch " + name);
HashMap<String, AtomType> types = forceField.getAtomTypes(name);
if (types.isEmpty()) {
isProtein = false;
break;
} else {
logger.info(" Patch found for non-standard amino acid " + name);
}
}
}
/**
* If all the residues in this chain have known amino acids
* names, then attempt to assign atom types.
*/
if (isProtein) {
try {
logger.info(format(" Amino acid chain %s", polymer.getName()));
double dist = properties.getDouble("chainbreak", 3.0);
// Detect main chain breaks!
List<List<Residue>> subChains = findChainBreaks(residues, dist);
for (List<Residue> subChain : subChains) {
assignAminoAcidAtomTypes(subChain);
}
} catch (MissingHeavyAtomException missingHeavyAtomException) {
logger.severe(missingHeavyAtomException.toString());
} catch (MissingAtomTypeException missingAtomTypeException) {
logger.severe(missingAtomTypeException.toString());
}
continue;
}
/**
* Check if all residues have known nucleic acids names.
*/
boolean isNucleicAcid = true;
for (int residueNumber = 0; residueNumber < numberOfResidues; residueNumber++) {
Residue residue = residues.get(residueNumber);
String name = residue.getName().toUpperCase();
/**
* Convert 1 and 2-character nucleic acid names to
* 3-character names.
*/
if (name.length() == 1) {
if (name.equals("A")) {
name = NucleicAcid3.ADE.toString();
} else if (name.equals("C")) {
name = NucleicAcid3.CYT.toString();
} else if (name.equals("G")) {
name = NucleicAcid3.GUA.toString();
} else if (name.equals("T")) {
name = NucleicAcid3.THY.toString();
} else if (name.equals("U")) {
name = NucleicAcid3.URI.toString();
}
} else if (name.length() == 2) {
if (name.equals("YG")) {
name = NucleicAcid3.YYG.toString();
}
}
residue.setName(name);
NucleicAcid3 nucleicAcid = null;
for (NucleicAcid3 nucleic : nucleicAcidList) {
String nuc3 = nucleic.toString();
nuc3 = nuc3.substring(nuc3.length() - 3);
if (nuc3.equalsIgnoreCase(name)) {
nucleicAcid = nucleic;
break;
}
}
if (nucleicAcid == null) {
logger.info(format("Nucleic acid was not recognized %s.", name));
isNucleicAcid = false;
break;
}
}
/**
* If all the residues in this chain have known nucleic acids
* names, then attempt to assign atom types.
*/
if (isNucleicAcid) {
try {
logger.info(format(" Nucleic acid chain %s", polymer.getName()));
assignNucleicAcidAtomTypes(residues, forceField, bondList);
} catch (MissingHeavyAtomException | MissingAtomTypeException e) {
logger.severe(e.toString());
}
}
}
}
// Assign ion atom types.
ArrayList<MSNode> ions = activeMolecularAssembly.getIons();
if (ions != null && ions.size() > 0) {
logger.info(format(" Assigning atom types for %d ions.", ions.size()));
for (MSNode m : ions) {
Molecule ion = (Molecule) m;
String name = ion.getResidueName().toUpperCase();
HetAtoms hetatm = HetAtoms.valueOf(name);
Atom atom = ion.getAtomList().get(0);
if (ion.getAtomList().size() != 1) {
logger.severe(format(" Check residue %s of chain %s.", ion.toString(), ion.getChainID()));
}
try {
switch(hetatm) {
case NA:
atom.setAtomType(findAtomType(2003));
break;
case K:
atom.setAtomType(findAtomType(2004));
break;
case MG:
case MG2:
atom.setAtomType(findAtomType(2005));
break;
case CA:
case CA2:
atom.setAtomType(findAtomType(2006));
break;
case CL:
atom.setAtomType(findAtomType(2007));
break;
case ZN:
case ZN2:
atom.setAtomType(findAtomType(2008));
break;
case BR:
atom.setAtomType(findAtomType(2009));
break;
default:
logger.severe(format(" Check residue %s of chain %s.", ion.toString(), ion.getChainID()));
}
} catch (Exception e) {
String message = "Error assigning atom types.";
logger.log(Level.SEVERE, message, e);
}
}
}
// Assign water atom types.
ArrayList<MSNode> water = activeMolecularAssembly.getWaters();
if (water != null && water.size() > 0) {
logger.info(format(" Assigning atom types for %d waters.", water.size()));
for (MSNode m : water) {
Molecule wat = (Molecule) m;
try {
Atom O = buildHeavy(wat, "O", null, 2001);
Atom H1 = buildHydrogen(wat, "H1", O, 0.96e0, null, 109.5e0, null, 120.0e0, 0, 2002);
H1.setHetero(true);
Atom H2 = buildHydrogen(wat, "H2", O, 0.96e0, H1, 109.5e0, null, 120.0e0, 0, 2002);
H2.setHetero(true);
} catch (Exception e) {
String message = "Error assigning atom types to a water.";
logger.log(Level.SEVERE, message, e);
}
}
}
// Assign small molecule atom types.
ArrayList<Molecule> molecules = activeMolecularAssembly.getMolecules();
for (MSNode m : molecules) {
Molecule molecule = (Molecule) m;
String moleculeName = molecule.getResidueName();
logger.info(" Attempting to patch " + moleculeName);
ArrayList<Atom> moleculeAtoms = molecule.getAtomList();
boolean patched = true;
HashMap<String, AtomType> types = forceField.getAtomTypes(moleculeName);
/**
* Assign atom types for all known atoms.
*/
for (Atom atom : moleculeAtoms) {
String atomName = atom.getName().toUpperCase();
AtomType atomType = types.get(atomName);
if (atomType == null) {
logger.info(" No atom type was found for " + atomName + " of " + moleculeName + ".");
patched = false;
break;
} else {
logger.fine(" " + atom.toString() + " -> " + atomType.toString());
atom.setAtomType(atomType);
types.remove(atomName);
}
}
/**
* Create missing hydrogen atoms. Check for missing heavy atoms.
*/
if (patched && !types.isEmpty()) {
for (AtomType type : types.values()) {
if (type.atomicNumber != 1) {
logger.info(" Missing heavy atom " + type.name);
patched = false;
break;
}
}
}
// Create bonds between known atoms.
if (patched) {
for (Atom atom : moleculeAtoms) {
String atomName = atom.getName();
String[] bonds = forceField.getBonds(moleculeName, atomName);
if (bonds != null) {
for (String name : bonds) {
Atom atom2 = molecule.getAtom(name);
if (atom2 != null && !atom.isBonded(atom2)) {
buildBond(atom, atom2);
}
}
}
}
}
// Create missing hydrogen atoms.
if (patched && !types.isEmpty()) {
// Create a hashmap of the molecule's atoms
HashMap<String, Atom> atomMap = new HashMap<String, Atom>();
for (Atom atom : moleculeAtoms) {
atomMap.put(atom.getName().toUpperCase(), atom);
}
for (String atomName : types.keySet()) {
AtomType type = types.get(atomName);
String[] bonds = forceField.getBonds(moleculeName, atomName.toUpperCase());
if (bonds == null || bonds.length != 1) {
patched = false;
logger.info(" Check biotype for hydrogen " + type.name + ".");
break;
}
// Get the heavy atom the hydrogen is bonded to.
Atom ia = atomMap.get(bonds[0].toUpperCase());
Atom hydrogen = new Atom(0, atomName, ia.getAltLoc(), new double[3], ia.getResidueName(), ia.getResidueNumber(), ia.getChainID(), ia.getOccupancy(), ia.getTempFactor(), ia.getSegID());
logger.fine(" Created hydrogen " + atomName + ".");
hydrogen.setAtomType(type);
hydrogen.setHetero(true);
molecule.addMSNode(hydrogen);
int valence = ia.getAtomType().valence;
List<Bond> aBonds = ia.getBonds();
int numBonds = aBonds.size();
/**
* Try to find the following configuration: ib-ia-ic
*/
Atom ib = null;
Atom ic = null;
Atom id = null;
if (numBonds > 0) {
Bond bond = aBonds.get(0);
ib = bond.get1_2(ia);
}
if (numBonds > 1) {
Bond bond = aBonds.get(1);
ic = bond.get1_2(ia);
}
if (numBonds > 2) {
Bond bond = aBonds.get(2);
id = bond.get1_2(ia);
}
/**
* Building the hydrogens depends on hybridization and the
* locations of other bonded atoms.
*/
logger.fine(" Bonding " + atomName + " to " + ia.getName() + " (" + numBonds + " of " + valence + ").");
switch(valence) {
case 4:
switch(numBonds) {
case 3:
// Find the average coordinates of atoms ib, ic and id.
double[] b = ib.getXYZ(null);
double[] c = ib.getXYZ(null);
double[] d = ib.getXYZ(null);
double[] a = new double[3];
a[0] = (b[0] + c[0] + d[0]) / 3.0;
a[1] = (b[1] + c[1] + d[1]) / 3.0;
a[2] = (b[2] + c[2] + d[2]) / 3.0;
// Place the hydrogen at chiral position #1.
intxyz(hydrogen, ia, 1.0, ib, 109.5, ic, 109.5, 1);
double[] e1 = new double[3];
hydrogen.getXYZ(e1);
double[] ret = new double[3];
diff(a, e1, ret);
double l1 = r(ret);
// Place the hydrogen at chiral position #2.
intxyz(hydrogen, ia, 1.0, ib, 109.5, ic, 109.5, -1);
double[] e2 = new double[3];
hydrogen.getXYZ(e2);
diff(a, e2, ret);
double l2 = r(ret);
// Revert to #1 if it is farther from the average.
if (l1 > l2) {
hydrogen.setXYZ(e1);
}
break;
case 2:
intxyz(hydrogen, ia, 1.0, ib, 109.5, ic, 109.5, 0);
break;
case 1:
intxyz(hydrogen, ia, 1.0, ib, 109.5, null, 0.0, 0);
break;
case 0:
intxyz(hydrogen, ia, 1.0, null, 0.0, null, 0.0, 0);
break;
default:
logger.info(" Check biotype for hydrogen " + atomName + ".");
patched = false;
}
break;
case 3:
switch(numBonds) {
case 2:
intxyz(hydrogen, ia, 1.0, ib, 120.0, ic, 0.0, 0);
break;
case 1:
intxyz(hydrogen, ia, 1.0, ib, 120.0, null, 0.0, 0);
break;
case 0:
intxyz(hydrogen, ia, 1.0, null, 0.0, null, 0.0, 0);
break;
default:
logger.info(" Check biotype for hydrogen " + atomName + ".");
patched = false;
}
break;
case 2:
switch(numBonds) {
case 1:
intxyz(hydrogen, ia, 1.0, ib, 120.0, null, 0.0, 0);
break;
case 0:
intxyz(hydrogen, ia, 1.0, null, 0.0, null, 0.0, 0);
break;
default:
logger.info(" Check biotype for hydrogen " + atomName + ".");
patched = false;
}
break;
case 1:
switch(numBonds) {
case 0:
intxyz(hydrogen, ia, 1.0, null, 0.0, null, 0.0, 0);
break;
default:
logger.info(" Check biotype for hydrogen " + atomName + ".");
patched = false;
}
break;
default:
logger.info(" Check biotype for hydrogen " + atomName + ".");
patched = false;
}
if (!patched) {
break;
} else {
buildBond(ia, hydrogen);
}
}
}
if (!patched) {
logger.log(Level.WARNING, format(" Deleting unrecognized molecule %s.", m.toString()));
activeMolecularAssembly.deleteMolecule((Molecule) m);
} else {
logger.info(" Patch for " + moleculeName + " succeeded.");
}
}
resolvePolymerLinks(molecules);
}
use of ffx.potential.parameters.AtomType in project ffx by mjschnie.
the class BiojavaFilter method assignNucleicAcidAtomTypes.
/**
* Assign atom types for a nucleic acid polymer.
*
* @param residues
* @throws ffx.potential.parsers.PDBFilter.MissingHeavyAtomException
*/
private void assignNucleicAcidAtomTypes(List<Residue> residues) throws MissingHeavyAtomException, MissingAtomTypeException {
/**
* A reference to the O3* atom of the previous base.
*/
Atom pO3s = null;
/**
* Loop over residues.
*/
int numberOfResidues = residues.size();
for (int residueNumber = 0; residueNumber < numberOfResidues; residueNumber++) {
/**
* Match the residue name to a known nucleic acid residue.
*/
Residue residue = residues.get(residueNumber);
String residueName = residue.getName().toUpperCase();
NucleicAcid3 nucleicAcid = null;
int naNumber = -1;
for (NucleicAcid3 nucleic : nucleicAcidList) {
naNumber++;
String nuc3 = nucleic.toString();
nuc3 = nuc3.substring(nuc3.length() - 3);
if (nuc3.equalsIgnoreCase(residueName)) {
nucleicAcid = nucleic;
break;
}
}
/**
* Do atom name conversions.
*/
List<Atom> resAtoms = residue.getAtomList();
int natoms = resAtoms.size();
for (int i = 0; i < natoms; i++) {
Atom atom = resAtoms.get(i);
String name = atom.getName();
name = name.replace('*', '\'');
// name = name.replace('D', 'H');
atom.setName(name);
}
/**
* Check if this is a 3' phosphate being listed as its own residue.
*/
/*if (residue.getAtomList().size() == 1) {
Atom P3s = (Atom) residue.getAtomNode("NA_P");
if (P3s != null) {
Residue prevResidue = residue.getPreviousResidue();
if (prevResidue != null) {
Atom O2sPrev = (Atom) prevResidue.getAtomNode("NA_O2\'");
if (O2sPrev == null) {
P3s = buildHeavy(prevResidue, "P3s", null, 1247);
} else {
P3s = buildHeavy(prevResidue, "P3s", null, 1235);
}
} else {
return;
}
} else {
return;
}
}*/
/**
* Check if the sugar is deoxyribose and change the residue name if
* necessary.
*/
boolean isDNA = false;
Atom O2s = (Atom) residue.getAtomNode("O2\'");
if (O2s == null) {
/**
* Assume deoxyribose (DNA) since there is an O2* atom.
*/
isDNA = true;
if (!residueName.startsWith("D")) {
switch(nucleicAcid) {
case ADE:
nucleicAcid = NucleicAcid3.DAD;
residueName = "DAD";
residue.setName(residueName);
break;
case CYT:
nucleicAcid = NucleicAcid3.DCY;
residueName = "DCY";
residue.setName(residueName);
break;
case GUA:
nucleicAcid = NucleicAcid3.DGU;
residueName = "DGU";
residue.setName(residueName);
break;
case THY:
nucleicAcid = NucleicAcid3.DTY;
residueName = "DTY";
residue.setName(residueName);
break;
default:
}
}
} else /**
* Assume ribose (RNA) since there is an O2* atom.
*/
if (residueName.startsWith("D")) {
switch(nucleicAcid) {
case DAD:
nucleicAcid = NucleicAcid3.ADE;
residueName = "ADE";
residue.setName(residueName);
break;
case DCY:
nucleicAcid = NucleicAcid3.CYT;
residueName = "CYT";
residue.setName(residueName);
break;
case DGU:
nucleicAcid = NucleicAcid3.GUA;
residueName = "GUA";
residue.setName(residueName);
break;
case DTY:
nucleicAcid = NucleicAcid3.THY;
residueName = "THY";
residue.setName(residueName);
break;
default:
}
}
/**
* Set a position flag.
*/
ResiduePosition position = MIDDLE_RESIDUE;
if (residueNumber == 0) {
position = FIRST_RESIDUE;
} else if (residueNumber == numberOfResidues - 1) {
position = LAST_RESIDUE;
}
/**
* Build the phosphate atoms of the current residue.
*/
Atom phosphate = null;
Atom O5s = null;
if (position == FIRST_RESIDUE) {
/**
* The 5' O5' oxygen of the nucleic acid is generally terminated
* by 1.) A phosphate group PO3 (-3). 2.) A hydrogen.
*
* If the base has phosphate atom we will assume a PO3 group.
*/
phosphate = (Atom) residue.getAtomNode("P");
if (phosphate != null) {
if (isDNA) {
phosphate = buildHeavy(residue, "P", null, 1247);
buildHeavy(residue, "OP1", phosphate, 1248);
buildHeavy(residue, "OP2", phosphate, 1248);
buildHeavy(residue, "OP3", phosphate, 1248);
O5s = buildHeavy(residue, "O5\'", phosphate, 1246);
} else {
phosphate = buildHeavy(residue, "P", null, 1235);
buildHeavy(residue, "OP1", phosphate, 1236);
buildHeavy(residue, "OP2", phosphate, 1236);
buildHeavy(residue, "OP3", phosphate, 1236);
O5s = buildHeavy(residue, "O5\'", phosphate, 1234);
}
} else if (isDNA) {
O5s = buildHeavy(residue, "O5\'", phosphate, 1244);
} else {
O5s = buildHeavy(residue, "O5\'", phosphate, 1232);
}
} else {
phosphate = buildHeavy(residue, "P", pO3s, NA_P[naNumber]);
buildHeavy(residue, "OP1", phosphate, NA_OP[naNumber]);
buildHeavy(residue, "OP2", phosphate, NA_OP[naNumber]);
O5s = buildHeavy(residue, "O5\'", phosphate, NA_O5[naNumber]);
}
/**
* Build the ribose sugar atoms of the current base.
*/
Atom C5s = buildHeavy(residue, "C5\'", O5s, NA_C5[naNumber]);
Atom C4s = buildHeavy(residue, "C4\'", C5s, NA_C4[naNumber]);
Atom O4s = buildHeavy(residue, "O4\'", C4s, NA_O4[naNumber]);
Atom C1s = buildHeavy(residue, "C1\'", O4s, NA_C1[naNumber]);
Atom C3s = buildHeavy(residue, "C3\'", C4s, NA_C3[naNumber]);
Atom C2s = buildHeavy(residue, "C2\'", C3s, NA_C2[naNumber]);
buildBond(C2s, C1s);
Atom O3s = null;
if (position == LAST_RESIDUE || numberOfResidues == 1) {
if (isDNA) {
O3s = buildHeavy(residue, "O3\'", C3s, 1249);
} else {
O3s = buildHeavy(residue, "O3\'", C3s, 1237);
}
} else {
O3s = buildHeavy(residue, "O3\'", C3s, NA_O3[naNumber]);
}
if (!isDNA) {
O2s = buildHeavy(residue, "O2\'", C2s, NA_O2[naNumber]);
}
/**
* Build the backbone hydrogen atoms.
*/
if (position == FIRST_RESIDUE && phosphate == null) {
buildHydrogen(residue, "H5T", O5s, 1.00e0, C5s, 109.5e0, C4s, 180.0e0, 0, NA_H5T[naNumber]);
}
buildHydrogen(residue, "H5\'1", C5s, 1.09e0, O5s, 109.5e0, C4s, 109.5e0, 1, NA_H51[naNumber]);
buildHydrogen(residue, "H5\'2", C5s, 1.09e0, O5s, 109.5e0, C4s, 109.5e0, -1, NA_H52[naNumber]);
buildHydrogen(residue, "H4\'", C4s, 1.09e0, C5s, 109.5e0, C3s, 109.5e0, -1, NA_H4[naNumber]);
buildHydrogen(residue, "H3\'", C3s, 1.09e0, C4s, 109.5e0, C2s, 109.5e0, -1, NA_H3[naNumber]);
if (isDNA) {
buildHydrogen(residue, "H2\'1", C2s, 1.09e0, C3s, 109.5e0, C1s, 109.5e0, -1, NA_H21[naNumber]);
buildHydrogen(residue, "H2\'2", C2s, 1.09e0, C3s, 109.5e0, C1s, 109.5e0, 1, NA_H22[naNumber]);
} else {
buildHydrogen(residue, "H2\'", C2s, 1.09e0, C3s, 109.5e0, C1s, 109.5e0, -1, NA_H21[naNumber]);
// Add the NA_O2' Methyl for OMC and OMG
if (nucleicAcid == NucleicAcid3.OMC || nucleicAcid == NucleicAcid3.OMG) {
Atom CM2 = buildHeavy(residue, "CM2", O2s, 1427);
Atom HM21 = buildHydrogen(residue, "HM21", CM2, 1.08e0, O2s, 109.5e0, C2s, 0.0e0, 0, 1428);
buildHydrogen(residue, "HM22", CM2, 1.08e0, O2s, 109.5e0, HM21, 109.5e0, 1, 1429);
buildHydrogen(residue, "HM23", CM2, 1.08e0, O2s, 109.5e0, HM21, 109.5e0, -1, 1430);
} else {
buildHydrogen(residue, "HO\'", O2s, 1.00e0, C2s, 109.5e0, C3s, 180.0e0, 0, NA_H22[naNumber]);
}
}
buildHydrogen(residue, "H1\'", C1s, 1.09e0, O4s, 109.5e0, C2s, 109.5e0, -1, NA_H1[naNumber]);
if (position == LAST_RESIDUE || numberOfResidues == 1) {
buildHydrogen(residue, "H3T", O3s, 1.00e0, C3s, 109.5e0, C4s, 180.0e0, 0, NA_H3T[naNumber]);
// Else, if it is terminated by a 3' phosphate cap:
// Will need to see how PDB would label a 3' phosphate cap.
}
/**
* Build the nucleic acid base.
*/
try {
assignNucleicAcidBaseAtomTypes(nucleicAcid, residue, C1s, O4s, C2s);
} catch (MissingHeavyAtomException missingHeavyAtomException) {
logger.throwing(PDBFilter.class.getName(), "assignNucleicAcidAtomTypes", missingHeavyAtomException);
throw missingHeavyAtomException;
}
/**
* Do some checks on the current base to make sure all atoms have
* been assigned an atom type.
*/
resAtoms = residue.getAtomList();
for (Atom atom : resAtoms) {
AtomType atomType = atom.getAtomType();
if (atomType == null) {
MissingAtomTypeException missingAtomTypeException = new MissingAtomTypeException(residue, atom);
logger.throwing(PDBFilter.class.getName(), "assignNucleicAcidAtomTypes", missingAtomTypeException);
throw missingAtomTypeException;
}
int numberOfBonds = atom.getNumBonds();
if (numberOfBonds != atomType.valence) {
if (atom == O3s && numberOfBonds == atomType.valence - 1 && position != LAST_RESIDUE && numberOfResidues != 1) {
continue;
}
logger.log(Level.WARNING, format(" An atom for residue %s has the wrong number of bonds:\n %s", residueName, atom.toString()));
logger.log(Level.WARNING, format(" Expected: %d Actual: %d.", atomType.valence, numberOfBonds));
}
}
/**
* Save a reference to the current O3* oxygen.
*/
pO3s = O3s;
}
}
Aggregations