Search in sources :

Example 21 with Bond

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;
        }
    }
}
Also used : PolarizeType(ffx.potential.parameters.PolarizeType) ArrayList(java.util.ArrayList) Bond(ffx.potential.bonded.Bond) Atom(ffx.potential.bonded.Atom)

Example 22 with Bond

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);
        }
    }
}
Also used : PolarizeType(ffx.potential.parameters.PolarizeType) Bond(ffx.potential.bonded.Bond) Atom(ffx.potential.bonded.Atom)

Example 23 with Bond

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()));
}
Also used : ROLS(ffx.potential.bonded.ROLS) Angle(ffx.potential.bonded.Angle) Bond(ffx.potential.bonded.Bond) Atom(ffx.potential.bonded.Atom)

Example 24 with Bond

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;
}
Also used : Angle(ffx.potential.bonded.Angle) HashMap(java.util.HashMap) ArrayList(java.util.ArrayList) Bond(ffx.potential.bonded.Bond) Torsion(ffx.potential.bonded.Torsion) Atom(ffx.potential.bonded.Atom)

Example 25 with Bond

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;
}
Also used : Bond(ffx.potential.bonded.Bond) Atom(ffx.potential.bonded.Atom)

Aggregations

Bond (ffx.potential.bonded.Bond)38 Atom (ffx.potential.bonded.Atom)37 IOException (java.io.IOException)12 ArrayList (java.util.ArrayList)12 Polymer (ffx.potential.bonded.Polymer)10 Residue (ffx.potential.bonded.Residue)10 File (java.io.File)10 Angle (ffx.potential.bonded.Angle)9 Crystal (ffx.crystal.Crystal)8 MissingAtomTypeException (ffx.potential.bonded.BondedUtils.MissingAtomTypeException)8 MissingHeavyAtomException (ffx.potential.bonded.BondedUtils.MissingHeavyAtomException)8 MSNode (ffx.potential.bonded.MSNode)7 Molecule (ffx.potential.bonded.Molecule)7 AtomType (ffx.potential.parameters.AtomType)7 BufferedWriter (java.io.BufferedWriter)6 FileWriter (java.io.FileWriter)6 RestraintBond (ffx.potential.bonded.RestraintBond)5 Torsion (ffx.potential.bonded.Torsion)5 CoordRestraint (ffx.potential.nonbonded.CoordRestraint)5 MolecularAssembly (ffx.potential.MolecularAssembly)4