use of ffx.potential.bonded.Bond in project ffx by mjschnie.
the class ParticleMeshEwald method assignPolarizationGroups.
protected void assignPolarizationGroups() {
/**
* Find directly connected group members for each atom.
*/
List<Integer> group = new ArrayList<>();
for (int i = 0; i < nAtoms; i++) {
Atom a = atoms[i];
if (a.getIndex() - 1 != i) {
logger.severe(" Atom indexing is not consistent in PME.");
}
}
for (Atom ai : atoms) {
group.clear();
Integer index = ai.getIndex() - 1;
group.add(index);
PolarizeType polarizeType = ai.getPolarizeType();
if (polarizeType != null) {
if (polarizeType.polarizationGroup != null) {
growGroup(group, ai);
Collections.sort(group);
ip11[index] = new int[group.size()];
int j = 0;
for (int k : group) {
ip11[index][j++] = k;
}
} else {
ip11[index] = new int[group.size()];
int j = 0;
for (int k : group) {
ip11[index][j++] = k;
}
}
} else {
String message = "The polarize keyword was not found for atom " + (index + 1) + " with type " + ai.getType();
logger.severe(message);
}
}
/**
* Find 1-2 group relationships.
*/
int[] mask = new int[nAtoms];
List<Integer> list = new ArrayList<>();
List<Integer> keep = new ArrayList<>();
for (int i = 0; i < nAtoms; i++) {
mask[i] = -1;
}
for (int i = 0; i < nAtoms; i++) {
list.clear();
for (int j : ip11[i]) {
list.add(j);
mask[j] = i;
}
keep.clear();
for (int j : list) {
Atom aj = atoms[j];
ArrayList<Bond> bonds = aj.getBonds();
for (Bond b : bonds) {
Atom ak = b.get1_2(aj);
int k = ak.getIndex() - 1;
if (mask[k] != i) {
keep.add(k);
}
}
}
list.clear();
for (int j : keep) {
for (int k : ip11[j]) {
list.add(k);
}
}
Collections.sort(list);
ip12[i] = new int[list.size()];
int j = 0;
for (int k : list) {
ip12[i][j++] = k;
}
}
/**
* Find 1-3 group relationships.
*/
for (int i = 0; i < nAtoms; i++) {
mask[i] = -1;
}
for (int i = 0; i < nAtoms; i++) {
for (int j : ip11[i]) {
mask[j] = i;
}
for (int j : ip12[i]) {
mask[j] = i;
}
list.clear();
for (int j : ip12[i]) {
for (int k : ip12[j]) {
if (mask[k] != i) {
if (!list.contains(k)) {
list.add(k);
}
}
}
}
ip13[i] = new int[list.size()];
Collections.sort(list);
int j = 0;
for (int k : list) {
ip13[i][j++] = k;
}
}
}
use of ffx.potential.bonded.Bond in project ffx by mjschnie.
the class ParticleMeshEwald method growGroup.
/**
* A recursive method that checks all atoms bonded to the seed atom for
* inclusion in the polarization group. The method is called on each newly
* found group member.
*
* @param group XYZ indeces of current group members.
* @param seed The bonds of the seed atom are queried for inclusion in the
* group.
*/
private void growGroup(List<Integer> group, Atom seed) {
List<Bond> bonds = seed.getBonds();
for (Bond bi : bonds) {
Atom aj = bi.get1_2(seed);
int tj = aj.getType();
boolean added = false;
PolarizeType polarizeType = seed.getPolarizeType();
for (int type : polarizeType.polarizationGroup) {
if (type == tj) {
Integer index = aj.getIndex() - 1;
if (!group.contains(index)) {
group.add(index);
added = true;
break;
}
}
}
if (added) {
growGroup(group, aj);
}
}
}
use of ffx.potential.bonded.Bond in project ffx by mjschnie.
the class PhMD method meltdown.
private void meltdown() {
writeSnapshot(".meltdown-");
ffe.energy(false, true);
if (ffe.getBondEnergy() > 1000) {
for (ROLS rols : previousTarget.getBondList()) {
((Bond) rols).log();
}
}
if (ffe.getAngleEnergy() > 1000) {
for (ROLS rols : previousTarget.getAngleList()) {
((Angle) rols).log();
}
}
if (ffe.getVanDerWaalsEnergy() > 1000) {
for (Atom a1 : previousTarget.getAtomList()) {
for (Atom a2 : mola.getAtomArray()) {
if (a1 == a2 || a1.getBond(a2) != null) {
continue;
}
double dist = sqrt(pow((a1.getX() - a2.getX()), 2) + pow((a1.getY() - a2.getY()), 2) + pow((a1.getZ() - a2.getZ()), 2));
if (dist < 0.8 * (a1.getVDWR() + a2.getVDWR())) {
logger.warning(String.format("Close vdW contact for atoms: \n %s\n %s", a1, a2));
}
}
}
}
logger.severe(String.format("Temperature above critical threshold: %f", thermostat.getCurrentTemperature()));
}
use of ffx.potential.bonded.Bond in project ffx by mjschnie.
the class RosenbluthChiAllMove method createBackBondedMap.
/**
* Maps the back-bonded terms affected by key atoms in an amino acid. Here,
* 'key atom' refers to each new rotamer-torsion-completing atom. e.g. VAL
* has 1 key atom (CG1), ARG has 4 key atoms (CG,CD,NE,CZ). 'Back-bonded'
* means we only map terms that lead toward the backbone.
*/
private HashMap<Integer, BackBondedList> createBackBondedMap(AminoAcid3 name) {
HashMap<Integer, BackBondedList> map = new HashMap<>();
List<Atom> chain = new ArrayList<>();
Atom N = (Atom) target.getAtomNode("N");
Atom CA = (Atom) target.getAtomNode("CA");
Atom CB = (Atom) target.getAtomNode("CB");
List<Atom> keyAtoms = new ArrayList<>();
switch(name) {
case VAL:
{
Atom CG1 = (Atom) target.getAtomNode("CG1");
keyAtoms.add(CG1);
keyAtoms.add(CB);
break;
}
case LEU:
{
Atom CG = (Atom) target.getAtomNode("CG");
Atom CD1 = (Atom) target.getAtomNode("CD1");
keyAtoms.add(CG);
keyAtoms.add(CD1);
break;
}
case ILE:
{
Atom CD1 = (Atom) target.getAtomNode("CD1");
Atom CG1 = (Atom) target.getAtomNode("CG1");
keyAtoms.add(CD1);
keyAtoms.add(CG1);
break;
}
case SER:
{
Atom OG = (Atom) target.getAtomNode("OG");
Atom HG = (Atom) target.getAtomNode("HG");
keyAtoms.add(OG);
keyAtoms.add(HG);
break;
}
case THR:
{
Atom OG1 = (Atom) target.getAtomNode("OG1");
Atom HG1 = (Atom) target.getAtomNode("HG1");
keyAtoms.add(OG1);
keyAtoms.add(HG1);
break;
}
case CYX:
case CYD:
{
Atom SG = (Atom) target.getAtomNode("SG");
keyAtoms.add(SG);
break;
}
case PHE:
{
Atom CD1 = (Atom) target.getAtomNode("CD1");
Atom CG = (Atom) target.getAtomNode("CG");
keyAtoms.add(CG);
break;
}
case PRO:
{
// Not allowed yet.
Atom CD = (Atom) target.getAtomNode("CD");
Atom CG = (Atom) target.getAtomNode("CG");
keyAtoms.add(CG);
keyAtoms.add(CD);
break;
}
case TYR:
{
Atom CD1 = (Atom) target.getAtomNode("CD1");
Atom CE2 = (Atom) target.getAtomNode("CE2");
Atom CG = (Atom) target.getAtomNode("CG");
Atom CZ = (Atom) target.getAtomNode("CZ");
Atom OH = (Atom) target.getAtomNode("OH");
Atom HH = (Atom) target.getAtomNode("HH");
// SPECIAL CASE: have to create map manualy.
Bond b1 = CG.getBond(CB);
Angle a1 = CG.getAngle(CB, CA);
Torsion t1 = CG.getTorsion(CB, CA, N);
Bond b2 = CD1.getBond(CG);
Angle a2 = CD1.getAngle(CG, CB);
Torsion t2 = CD1.getTorsion(CG, CB, CA);
Bond b3 = HH.getBond(OH);
Angle a3 = HH.getAngle(OH, CZ);
Torsion t3 = HH.getTorsion(OH, CZ, CE2);
BackBondedList bbl1 = new BackBondedList(b1, a1, t1);
BackBondedList bbl2 = new BackBondedList(b2, a2, t2);
BackBondedList bbl3 = new BackBondedList(b3, a3, t3);
map.put(0, bbl1);
map.put(1, bbl2);
map.put(2, bbl3);
// Note the return here.
return map;
}
case TYD:
{
Atom CD1 = (Atom) target.getAtomNode("CD1");
Atom CG = (Atom) target.getAtomNode("CG");
keyAtoms.add(CG);
keyAtoms.add(CD1);
break;
}
case TRP:
{
Atom CD1 = (Atom) target.getAtomNode("CD1");
Atom CG = (Atom) target.getAtomNode("CG");
keyAtoms.add(CG);
keyAtoms.add(CD1);
break;
}
case HIS:
case HID:
case HIE:
{
Atom CG = (Atom) target.getAtomNode("CG");
Atom ND1 = (Atom) target.getAtomNode("ND1");
keyAtoms.add(CG);
keyAtoms.add(ND1);
break;
}
case ASP:
{
Atom CG = (Atom) target.getAtomNode("CG");
keyAtoms.add(CG);
break;
}
case ASH:
{
Atom CG = (Atom) target.getAtomNode("CG");
Atom OD1 = (Atom) target.getAtomNode("OD1");
keyAtoms.add(CG);
keyAtoms.add(OD1);
break;
}
case ASN:
{
Atom CG = (Atom) target.getAtomNode("CG");
Atom OD1 = (Atom) target.getAtomNode("OD1");
keyAtoms.add(CG);
keyAtoms.add(OD1);
break;
}
case GLU:
case GLH:
{
Atom CG = (Atom) target.getAtomNode("CG");
Atom CD = (Atom) target.getAtomNode("CD");
Atom OE1 = (Atom) target.getAtomNode("OE1");
keyAtoms.add(CG);
keyAtoms.add(CD);
keyAtoms.add(OE1);
break;
}
case GLN:
{
Atom CG = (Atom) target.getAtomNode("CG");
Atom CD = (Atom) target.getAtomNode("CD");
Atom OE1 = (Atom) target.getAtomNode("OE1");
keyAtoms.add(CG);
keyAtoms.add(CD);
keyAtoms.add(OE1);
break;
}
case MET:
{
Atom CG = (Atom) target.getAtomNode("CG");
Atom CE = (Atom) target.getAtomNode("CE");
Atom SD = (Atom) target.getAtomNode("SD");
keyAtoms.add(CG);
keyAtoms.add(SD);
keyAtoms.add(CE);
break;
}
case LYS:
case LYD:
{
Atom CD = (Atom) target.getAtomNode("CD");
Atom CE = (Atom) target.getAtomNode("CE");
Atom CG = (Atom) target.getAtomNode("CG");
Atom NZ = (Atom) target.getAtomNode("NZ");
keyAtoms.add(CG);
keyAtoms.add(CD);
keyAtoms.add(CE);
keyAtoms.add(NZ);
break;
}
case ARG:
{
Atom CD = (Atom) target.getAtomNode("CD");
Atom CG = (Atom) target.getAtomNode("CG");
Atom CZ = (Atom) target.getAtomNode("CZ");
Atom NE = (Atom) target.getAtomNode("NE");
keyAtoms.add(CG);
keyAtoms.add(CD);
keyAtoms.add(NE);
keyAtoms.add(CZ);
break;
}
default:
logger.severe(String.format("CBMC called on unsupported residue: %s", name.toString()));
}
// Build the chain and assign back-bonded terms.
chain.add(N);
chain.add(CA);
chain.add(CB);
chain.addAll(keyAtoms);
for (int i = 3; i < chain.size(); i++) {
Atom key = chain.get(i);
Bond bond = key.getBond(chain.get(i - 1));
Angle angle = key.getAngle(chain.get(i - 1), chain.get(i - 2));
Torsion torsion = key.getTorsion(chain.get(i - 1), chain.get(i - 2), chain.get(i - 3));
BackBondedList bbl = new BackBondedList(bond, angle, torsion);
map.put(i - 3, bbl);
// report.append(String.format(" BackBondedMap: %s\t\t%s\n", key, torsion));
}
return map;
}
use of ffx.potential.bonded.Bond in project ffx by mjschnie.
the class XRayEnergy method getBFactorRestraints.
/**
* determine similarity and non-zero B factor restraints (done independently
* of getBFactorGradients), affects atomic gradients
*
* @return energy of the restraint
*/
public double getBFactorRestraints() {
Atom a1, a2;
double b1, b2, bdiff;
double[] anisou1 = null;
double[] anisou2 = null;
double gradb;
double det1, det2;
double[] gradu = new double[6];
double e = 0.0;
for (Atom a : activeAtomArray) {
double biso = a.getTempFactor();
// ignore hydrogens!!!
if (a.getAtomicNumber() == 1) {
continue;
}
if (a.getAnisou(null) == null) {
// isotropic B restraint
// non-zero restraint: -kTln[Z], Z is ADP partition function
e += -3.0 * kTbNonzero * Math.log(biso / (4.0 * Math.PI));
gradb = -3.0 * kTbNonzero / biso;
a.addToTempFactorGradient(gradb);
// similarity harmonic restraint
ArrayList<Bond> bonds = a.getBonds();
for (Bond b : bonds) {
if (a.compareTo(b.getAtom(0)) == 0) {
a1 = b.getAtom(0);
a2 = b.getAtom(1);
} else {
a1 = b.getAtom(1);
a2 = b.getAtom(0);
}
if (a2.getAtomicNumber() == 1) {
continue;
}
if (a2.getIndex() < a1.getIndex()) {
continue;
}
if (!a1.getAltLoc().equals(' ') && !a1.getAltLoc().equals('A') && a2.getAltLoc().equals(' ')) {
continue;
}
b1 = a1.getTempFactor();
b2 = a2.getTempFactor();
// transform B similarity restraints to U scale
bdiff = b2u(b1 - b2);
e += kTbSimWeight * Math.pow(bdiff, 2.0);
gradb = 2.0 * kTbSimWeight * bdiff;
a1.addToTempFactorGradient(gradb);
a2.addToTempFactorGradient(-gradb);
}
} else {
// anisotropic B restraint
anisou1 = a.getAnisou(anisou1);
det1 = determinant3(anisou1);
// non-zero restraint: -kTln[Z], Z is ADP partition function
e += u2b(-kTbNonzero * Math.log(det1 * eightPI2 * Math.PI));
gradu[0] = u2b(-kTbNonzero * ((anisou1[1] * anisou1[2] - anisou1[5] * anisou1[5]) / det1));
gradu[1] = u2b(-kTbNonzero * ((anisou1[0] * anisou1[2] - anisou1[4] * anisou1[4]) / det1));
gradu[2] = u2b(-kTbNonzero * ((anisou1[0] * anisou1[1] - anisou1[3] * anisou1[3]) / det1));
gradu[3] = u2b(-kTbNonzero * ((2.0 * (anisou1[4] * anisou1[5] - anisou1[3] * anisou1[2])) / det1));
gradu[4] = u2b(-kTbNonzero * ((2.0 * (anisou1[3] * anisou1[5] - anisou1[4] * anisou1[1])) / det1));
gradu[5] = u2b(-kTbNonzero * ((2.0 * (anisou1[3] * anisou1[4] - anisou1[5] * anisou1[0])) / det1));
a.addToAnisouGradient(gradu);
// similarity harmonic restraint based on determinants
ArrayList<Bond> bonds = a.getBonds();
for (Bond b : bonds) {
if (a.compareTo(b.getAtom(0)) == 0) {
a1 = b.getAtom(0);
a2 = b.getAtom(1);
} else {
a1 = b.getAtom(1);
a2 = b.getAtom(0);
}
if (a2.getAtomicNumber() == 1) {
continue;
}
if (a2.getIndex() < a1.getIndex()) {
continue;
}
if (a2.getAnisou(null) == null) {
continue;
}
if (!a1.getAltLoc().equals(' ') && !a1.getAltLoc().equals('A') && a2.getAltLoc().equals(' ')) {
continue;
}
anisou2 = a2.getAnisou(anisou2);
det2 = determinant3(anisou2);
bdiff = det1 - det2;
e += eightPI23 * kTbSimWeight * Math.pow(bdiff, 2.0);
gradb = eightPI23 * 2.0 * kTbSimWeight * bdiff;
// parent atom
gradu[0] = gradb * (anisou1[1] * anisou1[2] - anisou1[5] * anisou1[5]);
gradu[1] = gradb * (anisou1[0] * anisou1[2] - anisou1[4] * anisou1[4]);
gradu[2] = gradb * (anisou1[0] * anisou1[1] - anisou1[3] * anisou1[3]);
gradu[3] = gradb * (2.0 * (anisou1[4] * anisou1[5] - anisou1[3] * anisou1[2]));
gradu[4] = gradb * (2.0 * (anisou1[3] * anisou1[5] - anisou1[4] * anisou1[1]));
gradu[5] = gradb * (2.0 * (anisou1[3] * anisou1[4] - anisou1[5] * anisou1[0]));
a1.addToAnisouGradient(gradu);
// bonded atom
gradu[0] = gradb * (anisou2[5] * anisou2[5] - anisou2[1] * anisou2[2]);
gradu[1] = gradb * (anisou2[4] * anisou2[4] - anisou2[0] * anisou2[2]);
gradu[2] = gradb * (anisou2[3] * anisou2[3] - anisou2[0] * anisou2[1]);
gradu[3] = gradb * (2.0 * (anisou2[3] * anisou2[2] - anisou2[4] * anisou2[5]));
gradu[4] = gradb * (2.0 * (anisou2[4] * anisou2[1] - anisou2[3] * anisou2[5]));
gradu[5] = gradb * (2.0 * (anisou2[5] * anisou2[0] - anisou2[3] * anisou2[4]));
a2.addToAnisouGradient(gradu);
}
}
}
return e;
}
Aggregations