Search in sources :

Example 11 with Bond

use of ffx.potential.bonded.Bond in project ffx by mjschnie.

the class PDBFilter method writeSIFTFile.

public boolean writeSIFTFile(File saveFile, boolean append, String[] resAndScore) {
    if (saveFile == null) {
        return false;
    }
    if (vdwH) {
        logger.info(" Printing hydrogens to van der Waals centers instead of nuclear locations.");
    }
    if (nSymOp != 0) {
        logger.info(String.format(" Printing atoms with symmetry operator %s", activeMolecularAssembly.getCrystal().spaceGroup.getSymOp(nSymOp).toString()));
    }
    /**
     * Create StringBuilders for ATOM, ANISOU and TER records that can be
     * reused.
     */
    StringBuilder sb = new StringBuilder("ATOM  ");
    StringBuilder anisouSB = new StringBuilder("ANISOU");
    StringBuilder terSB = new StringBuilder("TER   ");
    for (int i = 6; i < 80; i++) {
        sb.append(' ');
        anisouSB.append(' ');
        terSB.append(' ');
    }
    FileWriter fw;
    BufferedWriter bw;
    try {
        File newFile = saveFile;
        if (!append && !noVersioning) {
            newFile = version(saveFile);
        }
        activeMolecularAssembly.setFile(newFile);
        activeMolecularAssembly.setName(newFile.getName());
        if (logWrites) {
            logger.log(Level.INFO, " Saving {0}", newFile.getName());
        }
        fw = new FileWriter(newFile, append);
        bw = new BufferedWriter(fw);
        // =============================================================================
        // The CRYST1 record presents the unit cell parameters, space group, and Z
        // value. If the structure was not determined by crystallographic means, CRYST1
        // simply provides the unitary values, with an appropriate REMARK.
        // 
        // 7 - 15       Real(9.3)     a              a (Angstroms).
        // 16 - 24       Real(9.3)     b              b (Angstroms).
        // 25 - 33       Real(9.3)     c              c (Angstroms).
        // 34 - 40       Real(7.2)     alpha          alpha (degrees).
        // 41 - 47       Real(7.2)     beta           beta (degrees).
        // 48 - 54       Real(7.2)     gamma          gamma (degrees).
        // 56 - 66       LString       sGroup         Space  group.
        // 67 - 70       Integer       z              Z value.
        // =============================================================================
        Crystal crystal = activeMolecularAssembly.getCrystal();
        if (crystal != null && !crystal.aperiodic()) {
            Crystal c = crystal.getUnitCell();
            if (!listMode) {
                bw.write(format("CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f %10s\n", c.a, c.b, c.c, c.alpha, c.beta, c.gamma, padRight(c.spaceGroup.pdbName, 10)));
            } else {
                listOutput.add(format("CRYST1%9.3f%9.3f%9.3f%7.2f%7.2f%7.2f %10s", c.a, c.b, c.c, c.alpha, c.beta, c.gamma, padRight(c.spaceGroup.pdbName, 10)));
            }
        }
        // =============================================================================
        // The SSBOND record identifies each disulfide bond in protein and polypeptide
        // structures by identifying the two residues involved in the bond.
        // The disulfide bond distance is included after the symmetry operations at
        // the end of the SSBOND record.
        // 
        // 8 - 10        Integer         serNum       Serial number.
        // 12 - 14        LString(3)      "CYS"        Residue name.
        // 16             Character       chainID1     Chain identifier.
        // 18 - 21        Integer         seqNum1      Residue sequence number.
        // 22             AChar           icode1       Insertion code.
        // 26 - 28        LString(3)      "CYS"        Residue name.
        // 30             Character       chainID2     Chain identifier.
        // 32 - 35        Integer         seqNum2      Residue sequence number.
        // 36             AChar           icode2       Insertion code.
        // 60 - 65        SymOP           sym1         Symmetry oper for 1st resid
        // 67 - 72        SymOP           sym2         Symmetry oper for 2nd resid
        // 74 – 78        Real(5.2)      Length        Disulfide bond distance
        // 
        // If SG of cysteine is disordered then there are possible alternate linkages.
        // wwPDB practice is to put together all possible SSBOND records. This is
        // problematic because the alternate location identifier is not specified in
        // the SSBOND record.
        // =============================================================================
        int serNum = 1;
        Polymer[] polymers = activeMolecularAssembly.getChains();
        if (polymers != null) {
            for (Polymer polymer : polymers) {
                ArrayList<Residue> residues = polymer.getResidues();
                for (Residue residue : residues) {
                    if (residue.getName().equalsIgnoreCase("CYS")) {
                        List<Atom> cysAtoms = residue.getAtomList();
                        Atom SG1 = null;
                        for (Atom atom : cysAtoms) {
                            if (atom.getName().equalsIgnoreCase("SG")) {
                                SG1 = atom;
                                break;
                            }
                        }
                        List<Bond> bonds = SG1.getBonds();
                        for (Bond bond : bonds) {
                            Atom SG2 = bond.get1_2(SG1);
                            if (SG2.getName().equalsIgnoreCase("SG")) {
                                if (SG1.getIndex() < SG2.getIndex()) {
                                    bond.energy(false);
                                    if (!listMode) {
                                        bw.write(format("SSBOND %3d CYS %1s %4s    CYS %1s %4s %36s %5.2f\n", serNum++, SG1.getChainID().toString(), Hybrid36.encode(4, SG1.getResidueNumber()), SG2.getChainID().toString(), Hybrid36.encode(4, SG2.getResidueNumber()), "", bond.getValue()));
                                    } else {
                                        listOutput.add(format("SSBOND %3d CYS %1s %4s    CYS %1s %4s %36s %5.2f\n", serNum++, SG1.getChainID().toString(), Hybrid36.encode(4, SG1.getResidueNumber()), SG2.getChainID().toString(), Hybrid36.encode(4, SG2.getResidueNumber()), "", bond.getValue()));
                                    }
                                }
                            }
                        }
                    }
                }
            }
        }
        // =============================================================================
        // 
        // 7 - 11        Integer       serial       Atom serial number.
        // 13 - 16        Atom          name         Atom name.
        // 17             Character     altLoc       Alternate location indicator.
        // 18 - 20        Residue name  resName      Residue name.
        // 22             Character     chainID      Chain identifier.
        // 23 - 26        Integer       resSeq       Residue sequence number.
        // 27             AChar         iCode        Code for insertion of residues.
        // 31 - 38        Real(8.3)     x            Orthogonal coordinates for X in Angstroms.
        // 39 - 46        Real(8.3)     y            Orthogonal coordinates for Y in Angstroms.
        // 47 - 54        Real(8.3)     z            Orthogonal coordinates for Z in Angstroms.
        // 55 - 60        Real(6.2)     occupancy    Occupancy.
        // 61 - 66        Real(6.2)     tempFactor   Temperature factor.
        // 77 - 78        LString(2)    element      Element symbol, right-justified.
        // 79 - 80        LString(2)    charge       Charge  on the atom.
        // =============================================================================
        // 1         2         3         4         5         6         7
        // 123456789012345678901234567890123456789012345678901234567890123456789012345678
        // ATOM      1  N   ILE A  16      60.614  71.140 -10.592  1.00  7.38           N
        // ATOM      2  CA  ILE A  16      60.793  72.149  -9.511  1.00  6.91           C
        MolecularAssembly[] molecularAssemblies = this.getMolecularAssemblys();
        int serial = 1;
        // Loop over biomolecular chains
        if (polymers != null) {
            for (Polymer polymer : polymers) {
                currentSegID = polymer.getName();
                currentChainID = polymer.getChainID();
                sb.setCharAt(21, currentChainID);
                // Loop over residues
                ArrayList<Residue> residues = polymer.getResidues();
                for (Residue residue : residues) {
                    String resName = residue.getName();
                    if (resName.length() > 3) {
                        resName = resName.substring(0, 3);
                    }
                    int resID = residue.getResidueNumber();
                    int i = 0;
                    String[] entries = null;
                    for (; i < resAndScore.length; i++) {
                        entries = resAndScore[i].split("\\t");
                        if (!entries[0].equals(entries[0].replaceAll("\\D+", ""))) {
                            String[] subEntries = entries[0].split("[^0-9]");
                            entries[0] = subEntries[0];
                        }
                        if (entries[0].equals(String.valueOf(resID)) && !".".equals(entries[1])) {
                            break;
                        }
                    }
                    sb.replace(17, 20, padLeft(resName.toUpperCase(), 3));
                    sb.replace(22, 26, String.format("%4s", Hybrid36.encode(4, resID)));
                    // Loop over atoms
                    ArrayList<Atom> residueAtoms = residue.getAtomList();
                    boolean altLocFound = false;
                    for (Atom atom : residueAtoms) {
                        if (i != resAndScore.length) {
                            writeSIFTAtom(atom, serial++, sb, anisouSB, bw, entries[1]);
                        } else {
                            writeSIFTAtom(atom, serial++, sb, anisouSB, bw, null);
                        }
                        Character altLoc = atom.getAltLoc();
                        if (altLoc != null && !altLoc.equals(' ')) {
                            altLocFound = true;
                        }
                    }
                    // Write out alternate conformers
                    if (altLocFound) {
                        for (int ma = 1; ma < molecularAssemblies.length; ma++) {
                            MolecularAssembly altMolecularAssembly = molecularAssemblies[ma];
                            Polymer altPolymer = altMolecularAssembly.getPolymer(currentChainID, currentSegID, false);
                            Residue altResidue = altPolymer.getResidue(resName, resID, false);
                            residueAtoms = altResidue.getAtomList();
                            for (Atom atom : residueAtoms) {
                                if (atom.getAltLoc() != null && !atom.getAltLoc().equals(' ') && !atom.getAltLoc().equals('A')) {
                                    if (i != resAndScore.length) {
                                        writeSIFTAtom(atom, serial++, sb, anisouSB, bw, entries[1]);
                                    } else {
                                        writeSIFTAtom(atom, serial++, sb, anisouSB, bw, null);
                                    }
                                }
                            }
                        }
                    }
                }
                terSB.replace(6, 11, String.format("%5s", Hybrid36.encode(5, serial++)));
                terSB.replace(12, 16, "    ");
                terSB.replace(16, 26, sb.substring(16, 26));
                if (!listMode) {
                    bw.write(terSB.toString());
                    bw.newLine();
                } else {
                    listOutput.add(terSB.toString());
                }
            }
        }
        sb.replace(0, 6, "HETATM");
        sb.setCharAt(21, 'A');
        int resID = 1;
        Polymer polymer = activeMolecularAssembly.getPolymer('A', "A", false);
        if (polymer != null) {
            ArrayList<Residue> residues = polymer.getResidues();
            for (Residue residue : residues) {
                int resID2 = residue.getResidueNumber();
                if (resID2 >= resID) {
                    resID = resID2 + 1;
                }
            }
        }
        /**
         * Loop over molecules, ions and then water.
         */
        ArrayList<Molecule> molecules = activeMolecularAssembly.getMolecules();
        for (int i = 0; i < molecules.size(); i++) {
            Molecule molecule = (Molecule) molecules.get(i);
            Character chainID = molecule.getChainID();
            sb.setCharAt(21, chainID);
            String resName = molecule.getResidueName();
            if (resName.length() > 3) {
                resName = resName.substring(0, 3);
            }
            sb.replace(17, 20, padLeft(resName.toUpperCase(), 3));
            sb.replace(22, 26, String.format("%4s", Hybrid36.encode(4, resID)));
            ArrayList<Atom> moleculeAtoms = molecule.getAtomList();
            boolean altLocFound = false;
            for (Atom atom : moleculeAtoms) {
                writeSIFTAtom(atom, serial++, sb, anisouSB, bw, null);
                Character altLoc = atom.getAltLoc();
                if (altLoc != null && !altLoc.equals(' ')) {
                    altLocFound = true;
                }
            }
            // Write out alternate conformers
            if (altLocFound) {
                for (int ma = 1; ma < molecularAssemblies.length; ma++) {
                    MolecularAssembly altMolecularAssembly = molecularAssemblies[ma];
                    MSNode altmolecule = altMolecularAssembly.getMolecules().get(i);
                    moleculeAtoms = altmolecule.getAtomList();
                    for (Atom atom : moleculeAtoms) {
                        if (atom.getAltLoc() != null && !atom.getAltLoc().equals(' ') && !atom.getAltLoc().equals('A')) {
                            writeSIFTAtom(atom, serial++, sb, anisouSB, bw, null);
                        }
                    }
                }
            }
            resID++;
        }
        ArrayList<MSNode> ions = activeMolecularAssembly.getIons();
        for (int i = 0; i < ions.size(); i++) {
            Molecule ion = (Molecule) ions.get(i);
            Character chainID = ion.getChainID();
            sb.setCharAt(21, chainID);
            String resName = ion.getResidueName();
            if (resName.length() > 3) {
                resName = resName.substring(0, 3);
            }
            sb.replace(17, 20, padLeft(resName.toUpperCase(), 3));
            sb.replace(22, 26, String.format("%4s", Hybrid36.encode(4, resID)));
            ArrayList<Atom> ionAtoms = ion.getAtomList();
            boolean altLocFound = false;
            for (Atom atom : ionAtoms) {
                writeSIFTAtom(atom, serial++, sb, anisouSB, bw, null);
                Character altLoc = atom.getAltLoc();
                if (altLoc != null && !altLoc.equals(' ')) {
                    altLocFound = true;
                }
            }
            // Write out alternate conformers
            if (altLocFound) {
                for (int ma = 1; ma < molecularAssemblies.length; ma++) {
                    MolecularAssembly altMolecularAssembly = molecularAssemblies[ma];
                    MSNode altion = altMolecularAssembly.getIons().get(i);
                    ionAtoms = altion.getAtomList();
                    for (Atom atom : ionAtoms) {
                        if (atom.getAltLoc() != null && !atom.getAltLoc().equals(' ') && !atom.getAltLoc().equals('A')) {
                            writeSIFTAtom(atom, serial++, sb, anisouSB, bw, null);
                        }
                    }
                }
            }
            resID++;
        }
        ArrayList<MSNode> waters = activeMolecularAssembly.getWaters();
        for (int i = 0; i < waters.size(); i++) {
            Molecule water = (Molecule) waters.get(i);
            Character chainID = water.getChainID();
            sb.setCharAt(21, chainID);
            String resName = water.getResidueName();
            if (resName.length() > 3) {
                resName = resName.substring(0, 3);
            }
            sb.replace(17, 20, padLeft(resName.toUpperCase(), 3));
            sb.replace(22, 26, String.format("%4s", Hybrid36.encode(4, resID)));
            ArrayList<Atom> waterAtoms = water.getAtomList();
            boolean altLocFound = false;
            for (Atom atom : waterAtoms) {
                writeSIFTAtom(atom, serial++, sb, anisouSB, bw, null);
                Character altLoc = atom.getAltLoc();
                if (altLoc != null && !altLoc.equals(' ')) {
                    altLocFound = true;
                }
            }
            // Write out alternate conformers
            if (altLocFound) {
                for (int ma = 1; ma < molecularAssemblies.length; ma++) {
                    MolecularAssembly altMolecularAssembly = molecularAssemblies[ma];
                    MSNode altwater = altMolecularAssembly.getWaters().get(i);
                    waterAtoms = altwater.getAtomList();
                    for (Atom atom : waterAtoms) {
                        if (atom.getAltLoc() != null && !atom.getAltLoc().equals(' ') && !atom.getAltLoc().equals('A')) {
                            writeSIFTAtom(atom, serial++, sb, anisouSB, bw, null);
                        }
                    }
                }
            }
            resID++;
        }
        if (!listMode) {
            bw.write("END");
            bw.newLine();
        } else {
            listOutput.add("END");
        }
        bw.close();
    } catch (Exception e) {
        String message = "Exception writing to file: " + saveFile.toString();
        logger.log(Level.WARNING, message, e);
        return false;
    }
    return true;
}
Also used : FileWriter(java.io.FileWriter) BufferedWriter(java.io.BufferedWriter) MSNode(ffx.potential.bonded.MSNode) Polymer(ffx.potential.bonded.Polymer) Atom(ffx.potential.bonded.Atom) MissingHeavyAtomException(ffx.potential.bonded.BondedUtils.MissingHeavyAtomException) IOException(java.io.IOException) MissingAtomTypeException(ffx.potential.bonded.BondedUtils.MissingAtomTypeException) Molecule(ffx.potential.bonded.Molecule) MolecularAssembly(ffx.potential.MolecularAssembly) Residue(ffx.potential.bonded.Residue) Bond(ffx.potential.bonded.Bond) File(java.io.File) Crystal(ffx.crystal.Crystal)

Example 12 with Bond

use of ffx.potential.bonded.Bond 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);
}
Also used : NucleicAcid3(ffx.potential.bonded.ResidueEnumerations.NucleicAcid3) HashMap(java.util.HashMap) MissingHeavyAtomException(ffx.potential.bonded.BondedUtils.MissingHeavyAtomException) MissingAtomTypeException(ffx.potential.bonded.BondedUtils.MissingAtomTypeException) MSNode(ffx.potential.bonded.MSNode) AtomType(ffx.potential.parameters.AtomType) ResidueEnumerations.nucleicAcidList(ffx.potential.bonded.ResidueEnumerations.nucleicAcidList) ResidueEnumerations.aminoAcidList(ffx.potential.bonded.ResidueEnumerations.aminoAcidList) List(java.util.List) ArrayList(java.util.ArrayList) AminoAcid3(ffx.potential.bonded.ResidueEnumerations.AminoAcid3) Polymer(ffx.potential.bonded.Polymer) Atom(ffx.potential.bonded.Atom) MissingHeavyAtomException(ffx.potential.bonded.BondedUtils.MissingHeavyAtomException) IOException(java.io.IOException) MissingAtomTypeException(ffx.potential.bonded.BondedUtils.MissingAtomTypeException) Molecule(ffx.potential.bonded.Molecule) Residue(ffx.potential.bonded.Residue) Bond(ffx.potential.bonded.Bond)

Example 13 with Bond

use of ffx.potential.bonded.Bond in project ffx by mjschnie.

the class XYZFilter method writeFile.

/**
 * {@inheritDoc}
 */
@Override
public boolean writeFile(File saveFile, boolean append) {
    if (saveFile == null) {
        return false;
    }
    try {
        File newFile = saveFile;
        if (!append) {
            newFile = version(saveFile);
        }
        activeMolecularAssembly.setFile(newFile);
        activeMolecularAssembly.setName(newFile.getName());
        FileWriter fw;
        if (append && !newFile.exists()) {
            fw = new FileWriter(newFile);
        } else {
            fw = new FileWriter(newFile, append);
        }
        BufferedWriter bw = new BufferedWriter(fw);
        // XYZ File First Line
        int numberOfAtoms = activeMolecularAssembly.getAtomList().size();
        String output = format("%7d  %s\n", numberOfAtoms, activeMolecularAssembly.toString());
        bw.write(output);
        Crystal crystal = activeMolecularAssembly.getCrystal();
        if (!crystal.aperiodic()) {
            Crystal uc = crystal.getUnitCell();
            String params = String.format("%14.8f%14.8f%14.8f%14.8f%14.8f%14.8f\n", uc.a, uc.b, uc.c, uc.alpha, uc.beta, uc.gamma);
            bw.write(params);
        }
        Atom a2;
        StringBuilder line;
        StringBuilder[] lines = new StringBuilder[numberOfAtoms];
        // XYZ File Atom Lines
        ArrayList<Atom> atoms = activeMolecularAssembly.getAtomList();
        Vector3d offset = activeMolecularAssembly.getOffset();
        for (Atom a : atoms) {
            if (vdwH) {
                line = new StringBuilder(String.format("%7d %3s%14.8f%14.8f%14.8f%6d", a.getIndex(), a.getAtomType().name, a.getRedX() - offset.x, a.getRedY() - offset.y, a.getRedZ() - offset.z, a.getType()));
            } else {
                line = new StringBuilder(String.format("%7d %3s%14.8f%14.8f%14.8f%6d", a.getIndex(), a.getAtomType().name, a.getX() - offset.x, a.getY() - offset.y, a.getZ() - offset.z, a.getType()));
            }
            for (Bond b : a.getBonds()) {
                a2 = b.get1_2(a);
                line.append(format("%8d", a2.getIndex()));
            }
            lines[a.getIndex() - 1] = line.append("\n");
        }
        try {
            for (int i = 0; i < numberOfAtoms; i++) {
                bw.write(lines[i].toString());
            }
        } catch (IOException e) {
            String message = format(" Their was an unexpected error writing to %s.", getActiveMolecularSystem().toString());
            logger.log(Level.WARNING, message, e);
            return false;
        }
        bw.close();
        fw.close();
    } catch (IOException e) {
        String message = format(" Their was an unexpected error writing to %s.", getActiveMolecularSystem().toString());
        logger.log(Level.WARNING, message, e);
        return false;
    }
    return true;
}
Also used : FileWriter(java.io.FileWriter) IOException(java.io.IOException) Atom(ffx.potential.bonded.Atom) BufferedWriter(java.io.BufferedWriter) Vector3d(javax.vecmath.Vector3d) Bond(ffx.potential.bonded.Bond) File(java.io.File) Crystal(ffx.crystal.Crystal)

Example 14 with Bond

use of ffx.potential.bonded.Bond in project ffx by mjschnie.

the class BiojavaFilter method assignAtomTypes.

/**
 * Assign force field atoms types to common chemistries using "biotype"
 * records.
 */
public void assignAtomTypes() {
    /**
     * Create a new List to store bonds determined based on 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;
            if (!residues.isEmpty()) {
            // renameNTerminusHydrogens(residues.get(0)); Not safe to use until it distinguishes between true N-termini and N-terminal residues in general.
            }
            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;
                        renameNonstandardHydrogens(residue);
                        break;
                    }
                }
                // Check for a patch.
                if (!aa) {
                    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);
                } catch (MissingHeavyAtomException missingHeavyAtomException) {
                    logger.severe(missingHeavyAtomException.toString());
                } catch (MissingAtomTypeException missingAtomTypeException) {
                    logger.severe(missingAtomTypeException.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 {
                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.");
        }
    }
}
Also used : NucleicAcid3(ffx.potential.bonded.ResidueEnumerations.NucleicAcid3) HashMap(java.util.HashMap) MissingHeavyAtomException(ffx.potential.bonded.BondedUtils.MissingHeavyAtomException) MissingAtomTypeException(ffx.potential.bonded.BondedUtils.MissingAtomTypeException) MSNode(ffx.potential.bonded.MSNode) AtomType(ffx.potential.parameters.AtomType) ResidueEnumerations.nucleicAcidList(ffx.potential.bonded.ResidueEnumerations.nucleicAcidList) ResidueEnumerations.aminoAcidList(ffx.potential.bonded.ResidueEnumerations.aminoAcidList) ArrayList(java.util.ArrayList) List(java.util.List) AminoAcid3(ffx.potential.bonded.ResidueEnumerations.AminoAcid3) Polymer(ffx.potential.bonded.Polymer) HetAtoms(ffx.potential.parsers.PDBFilter.HetAtoms) Atom(ffx.potential.bonded.Atom) IOException(java.io.IOException) MissingHeavyAtomException(ffx.potential.bonded.BondedUtils.MissingHeavyAtomException) MissingAtomTypeException(ffx.potential.bonded.BondedUtils.MissingAtomTypeException) Molecule(ffx.potential.bonded.Molecule) Residue(ffx.potential.bonded.Residue) Bond(ffx.potential.bonded.Bond) SSBond(org.biojava.bio.structure.SSBond)

Example 15 with Bond

use of ffx.potential.bonded.Bond in project ffx by mjschnie.

the class MultipoleType method multipoleTypeFactory.

public static MultipoleType multipoleTypeFactory(Atom atom, ForceField forceField) {
    AtomType atomType = atom.getAtomType();
    if (atomType == null) {
        String message = " Multipoles can only be assigned to atoms that have been typed.";
        logger.severe(message);
        return null;
    }
    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);
        atom.setAxisAtoms(null);
        return multipoleType;
    }
    // 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 = forceField.getMultipoleType(key);
        if (multipoleType != null) {
            atom.setMultipoleType(multipoleType);
            atom.setAxisAtoms(atom2);
            return multipoleType;
        }
    }
    // 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) {
                atom.setMultipoleType(multipoleType);
                atom.setAxisAtoms(atom2, atom3);
                return multipoleType;
            }
        }
    }
    /**
     * 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) {
                    atom.setMultipoleType(multipoleType);
                    atom.setAxisAtoms(atom2, atom3, atom4);
                    return multipoleType;
                }
            }
            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) {
                        atom.setMultipoleType(multipoleType);
                        atom.setAxisAtoms(atom2, atom3, atom4);
                        return multipoleType;
                    }
                }
            }
        }
    }
    /**
     * 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) {
                    atom.setMultipoleType(multipoleType);
                    atom.setAxisAtoms(atom2, atom3);
                    return multipoleType;
                }
                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) {
                            atom.setMultipoleType(multipoleType);
                            atom.setAxisAtoms(atom2, atom3, atom4);
                            return multipoleType;
                        }
                    }
                }
            }
        }
    }
    return null;
}
Also used : Angle(ffx.potential.bonded.Angle) 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