Automatic Py3k merge for testing

pymatgen/io/lammps/data.py

@@ -241,7 +241,7 @@ class LammpsForceFieldData(LammpsData):
     Args:
         box_size (list): [[x_min,x_max], [y_min,y_max], [z_min,z_max]]
         atomic_masses (list): [ [atom type, atomic mass], ... ]
-        pair_coeffs (list): vdw coefficients,
+        pair_coeffs (list): pair coefficients,
             [[unique id, sigma, epsilon ], ... ]
         bond_coeffs (list): bond coefficients,
             [[unique id, value1, value2 ], ... ]
@@ -419,8 +419,7 @@ def get_atoms_data(mols, mols_number, molecule, atomic_masses_dict,
             charge = 0.0
             if hasattr(topologies[0], "charges"):
                 if topologies[mol_type - 1].charges:
-                    charge = topologies[mol_type - 1].charges[mol_atom_id - 1][
-                        0]
+                    charge = topologies[mol_type - 1].charges[mol_atom_id - 1]
             atoms_data.append([atom_id, mol_type, atom_type, charge,
                                site.x, site.y, site.z])
         return atoms_data, molid_to_atomid
@@ -534,7 +533,7 @@ def from_forcefield_and_topology(mols, mols_number, box_size, molecule,
         angle_coeffs, angle_map = LammpsForceFieldData.get_param_coeff(
             forcefield, "angles")
         pair_coeffs, _ = LammpsForceFieldData.get_param_coeff(forcefield,
-                                                              "vdws")
+                                                              "pairs")
         dihedral_coeffs, dihedral_map = LammpsForceFieldData.get_param_coeff(
             forcefield, "dihedrals")
         improper_coeffs, imdihedral_map = LammpsForceFieldData.get_param_coeff(

pymatgen/io/lammps/force_field.py

@@ -16,50 +16,51 @@
 from monty.json import MSONable
 
 __author__ = 'Kiran Mathew'
+__email__ = 'kmathew@lbl.gov'
+__credits__ = 'Brandon Wood'
 
 
 class ForceField(MSONable):
     """
     Stores force field information.
 
     Args:
-        atoms (OrderedDict): store atomic mass for each atom name.
+        atoms (Dict): store atomic mass for each atom name.
             { "atom name": atom mass, ... }
-        bonds (OrderedDict): store the bond distance (A) and spring constant (
+        bonds (Dict): store the bond distance (A) and spring constant (
             Kcal/molA2) for each bond.
             { ("atom name1", "atom name2"): [spring const, distance], ... }
-        angles (OrderedDict): store the bond angle and spring constant
+        angles (Dict): store the bond angle and spring constant
             (Kcal/mol*radian2).
             { ("atom name1", "atom name2", "atom name3"): [spring const, angle], ... }
-        dihedrals (OrderedDict): store the magnitude of torsion (Kcal/mol).
+        dihedrals (Dict): store the magnitude of torsion (Kcal/mol).
             { ("atom name1", "atom name2", "atom name3", "atom name4"): [
             function type, value, angle], ... }
-        imdihedrals (OrderedDict): store improper dihedral information.
+        imdihedrals (Dict): store improper dihedral information.
             similar to dihedrals but the gaff atom name1 and gaff atom name2
             are marked 'X'
-        vdws (OrderedDict): store the van der waal radius (A) and van der wall
-            depth for a given atom (Kcal/mol). Lennard-Jones parameters.
-            { "atom name": [sigma, epsilon], ... }
+        pairs (Dict): store pair coefficient info.
+            { ("atom name1", "atom name2"): [val1, val2, ..], ... }
     """
 
     def __init__(self, atoms, bonds, angles, dihedrals=None, imdihedrals=None,
-                 vdws=None, masses=None, charges=None):
+                 pairs=None):
         self.atoms = atoms
         self.bonds = bonds
         self.angles = angles
         self.dihedrals = dihedrals
         self.imdihedrals = imdihedrals
-        self.vdws = vdws
-        self.masses = masses
-        self.charges = OrderedDict() if charges is None else charges
+        self.pairs = pairs
 
     @staticmethod
     def from_file(filename):
         """
         Read in forcefield parameters from yaml file and return Forcefield
         object. Basically read in atom name mappings(key='Atoms'), bond
-        coefficients(key='Bond Coeffs'), angle coefficients(key='Angle Coeffs')
-        and the dihedral coefficients(key='Dihedral Coeffs').
+        coefficients(key='Bond Coeffs'), angle coefficients(key='Angle
+        Coeffs'), pair coefficients(key='Pair Coeffs'),
+        dihedral coefficients(key='Dihedral Coeffs') and the
+        improper dihedral coefficients(key='Improper Coeffs').
 
         Args:
             filename (string)
@@ -75,5 +76,9 @@ def from_file(filename):
                 tokens = k.split("-")
                 key = tuple(tokens) if len(tokens) > 1 else k
                 ff_data[coeff_key][key] = v
-        return ForceField(ff_data["Atoms"], ff_data["Bond Coeffs"],
-                          ff_data["Angle Coeffs"], ff_data["Dihedral Coeffs"])
+        return ForceField(ff_data["Atoms"],
+                          ff_data["Bond Coeffs"],
+                          ff_data["Angle Coeffs"],
+                          dihedrals=ff_data.get("Dihedral Coeffs", None),
+                          imdihedrals=ff_data.get("Improper Coeffs", None),
+                          pairs=ff_data.get("Pair Coeffs", None))

pymatgen/io/lammps/tests/test_lammps_forcefield_data.py

@@ -29,7 +29,7 @@ def setUpClass(cls):
                              ((u'hw', u'hw'), [553.0, 1.513])])
         angles = OrderedDict([((u'hw', u'ow', u'hw'), [0.0, 104.52])])
         vdws = OrderedDict([(u'hw', [0.0, 0.0]), (u'ow', [1.7683, 0.1520])])
-        forcefield = ForceField(atoms, bonds, angles, vdws=vdws)
+        forcefield = ForceField(atoms, bonds, angles, pairs=vdws)
         h2o_coords = [[9.626, 6.787, 12.673],
                       [9.626, 8.420, 12.673],
                       [10.203, 7.604, 12.673]]

pymatgen/io/lammps/tests/test_topology_forcefield.py

@@ -12,6 +12,7 @@
 
 __author__ = 'Kiran Mathew'
 __email__ = 'kmathew@lbl.gov'
+__credits__ = 'Brandon Wood'
 
 test_dir = os.path.join(os.path.dirname(__file__), "..", "..", "..", "..",
                         "test_files", "lammps")
@@ -21,6 +22,7 @@ class TestLammpsTopology(unittest.TestCase):
     @classmethod
     def setUpClass(cls):
         polymer = Molecule.from_file(os.path.join(test_dir, "polymer.xyz"))
+        cls.monomer = Molecule.from_file(os.path.join(test_dir, "monomer.xyz"))
         cls.topology = Topology.from_molecule(polymer, tol=0.1)
         cls.forcefield = ForceField.from_file(os.path.join(test_dir,
                                                        "ff_data.yaml"))
@@ -38,24 +40,36 @@ def test_topology(self):
                   [11, 10, 16, ('C', 'C', 'H')], [11, 10, 17, ('C', 'C', 'H')],
                   [10, 11, 12, ('C', 'C', 'O')], [10, 11, 18, ('C', 'C', 'H')],
                   [10, 11, 19, ('C', 'C', 'H')], [16, 10, 17, ('H', 'C', 'H')]]
-        dihedrals = [[6, 14, 8, 15, ('H', 'H', 'C', 'H')],
-                     [8, 9, 10, 11, ('C', 'O', 'C', 'C')],
-                     [8, 9, 10, 16, ('C', 'O', 'C', 'H')],
+        dihedrals = [[8, 9, 10, 16, ('C', 'O', 'C', 'H')],
                      [8, 9, 10, 17, ('C', 'O', 'C', 'H')],
                      [15, 8, 9, 10, ('H', 'C', 'O', 'C')],
                      [9, 10, 11, 12, ('O', 'C', 'C', 'O')],
                      [9, 10, 11, 18, ('O', 'C', 'C', 'H')],
                      [9, 10, 11, 19, ('O', 'C', 'C', 'H')],
                      [10, 11, 12, 13, ('C', 'C', 'O', 'C')],
-                     [10, 11, 19, 22, ('C', 'C', 'H', 'C')]]
+                     [10, 11, 19, 22, ('C', 'C', 'H', 'C')],
+                     [16, 10, 11, 12, ('H', 'C', 'C', 'O')],
+                     [16, 10, 11, 18, ('H', 'C', 'C', 'H')]]
         self.assertEqual(len(self.topology.atoms), 44)
         self.assertEqual(len(self.topology.bonds), 50)
         self.assertEqual(len(self.topology.angles), 123)
-        self.assertEqual(len(self.topology.dihedrals), 177)
+        self.assertEqual(len(self.topology.dihedrals), 175)
         self.assertEqual(self.topology.atoms[17:27], atoms)
         self.assertEqual(self.topology.bonds[17:27], bonds)
         self.assertEqual(self.topology.angles[17:27], angles)
         self.assertEqual(self.topology.dihedrals[17:27], dihedrals)
+        monomer_charges = [-0.118700, -0.279200, -0.032600, -0.118700,
+                          -0.279200, -0.032600, 0.086100, 0.086100,
+                           0.086100, 0.086100,0.086100, 0.086100,
+                           0.086100, 0.086100]
+        monomer_atoms = [['C', 'C'], ['O', 'O'], ['C', 'C'], ['C', 'C'],
+                         ['O', 'O'], ['C', 'C'], ['H', 'H'], ['H', 'H'],
+                         ['H', 'H'], ['H', 'H'], ['H', 'H'], ['H', 'H'],
+                         ['H', 'H'], ['H', 'H']]
+        self.monomer.add_site_property("charge", monomer_charges)
+        monomer_topology = Topology.from_molecule(self.monomer, tol=0.1)
+        self.assertEqual(monomer_topology.atoms, monomer_atoms)
+        self.assertEqual(monomer_topology.charges, monomer_charges)
 
     def test_forcefield(self):
         atoms = {'H': 'H', 'C': 'C', 'O': 'O'}
@@ -71,10 +85,16 @@ def test_forcefield(self):
                      (u'C', u'C', u'O', u'C'): [1.76, 0.67, 0.04, 0.0],
                      (u'H', u'C', u'C', u'H'): [0.0, 0.0, 0.28, 0.0],
                      (u'O', u'C', u'C', u'O'): [0.41, -2.1, -0.6, -0.82]}
+        pairs = {(u'C', u'C'): [14976.0, 0.3236, 637.6],
+                 (u'C', u'H'): [4320.0, 0.2928, 137.6],
+                 (u'C', u'O'): [33702.4, 0.2796, 503.0],
+                 (u'H', u'H'): [2649.6, 0.2674, 27.22],
+                 (u'O', u'O'): [75844.8, 0.2461, 396.9]}
         self.assertEqual(self.forcefield.atoms, atoms)
         self.assertEqual(self.forcefield.bonds, bonds)
         self.assertEqual(self.forcefield.angles, angles)
         self.assertEqual(self.forcefield.dihedrals, dihedrals)
+        self.assertEqual(self.forcefield.pairs, pairs)
 
 
 if __name__ == "__main__":

pymatgen/io/lammps/topology.py

@@ -14,34 +14,35 @@
 from pymatgen.core.bonds import CovalentBond
 
 __author__ = 'Kiran Mathew'
+__email__ = 'kmathew@lbl.gov'
+__credits__ = 'Brandon Wood'
 
 
 class Topology(object):
     """
     Args:
         atoms (list): map atom names to force field(ff) atom name,
             [['c', 'c1'],...]
-        charges (list): List of charges, [0.4, 0.7, ... ]
         bonds (list): List of bonds,
             [[i,j, bond_type], ... ] where i, j are integer(starts from 1)
             atom ids in the molecules and bond_type = (ff atomname_i, ff atomname_j)
         angles (list): List of angles,
             [[i,j,k, angle_type], ... ],
             angle_type = (ff atomname_i, ff atomname_j, ff atomname_k)
+        charges (list): List of charges, [0.4, 0.7, ... ]
         dihedrals (list): List of dihedrals,
             [[i,j,k,l, dihedral_type], ... ]
             dihedral_type = (ff atomname_i, ff atomname_j, ff atomname_k, ff atomname_l)
         imdihedrals (list): List of improper dihedrals,
             [['i,j,k,l, dihedral_type], ... ]
-
     """
 
     def __init__(self, atoms, bonds, angles, charges=None, dihedrals=None,
                  imdihedrals=None):
         self.atoms = atoms
-        self.charges = dict() if charges is None else charges
         self.bonds = bonds
         self.angles = angles
+        self.charges = [] if charges is None else charges
         self.dihedrals = dihedrals
         self.imdihedrals = imdihedrals
 
@@ -82,7 +83,7 @@ def from_molecule(molecule, tol=0.1):
             else:
                 for site1, site2 in itertools.product(bond1_sites,
                                                       bond2_sites):
-                    if CovalentBond.is_bonded(site1, site2):
+                    if CovalentBond.is_bonded(site1, site2, tol=tol):
                         bond1_sites.remove(site1)
                         bond2_sites.remove(site2)
                         dihedral = bond1_sites + [site1,
@@ -99,6 +100,8 @@ def from_molecule(molecule, tol=0.1):
                         break
         atoms = [[str(site.specie), str(site.specie)] for site in molecule]
         bonds = [[molecule.index(b.site1), molecule.index(b.site2),
-                  (str(b.site1.specie), str(b.site2.specie))] for b in
-                 bonds]
-        return Topology(atoms, bonds, angles, dihedrals=dihedrals)
+                  (str(b.site1.specie), str(b.site2.specie))] for b in bonds]
+        charges = None
+        if hasattr(molecule[0], "charge"):
+            charges = [site.charge for site in molecule]
+        return Topology(atoms, bonds, angles, charges=charges, dihedrals=dihedrals)

test_files/lammps/ff_data.yaml

@@ -17,6 +17,12 @@ Bond Coeffs:
   C-O: [1000, 1.4115]
   C-H: [1000, 1.1041]
   C-O: [1000, 1.4115]
+Pair Coeffs:
+  C-C: [14976.0, 0.3236, 637.6]
+  C-H: [4320.0, 0.2928, 137.6]
+  C-O: [33702.4, 0.2796, 503.0]
+  H-H: [2649.6, 0.2674, 27.22]
+  O-O: [75844.8, 0.2461, 396.9]
 Dihedral Coeffs:
   C-C-O-C: [1.76, 0.67, 0.04, 0.0]
   C-C-O-C: [1.76, 0.67, 0.04, 0.0]

test_files/lammps/monomer.xyz

@@ -0,0 +1,16 @@
+14
+	Energy:      13.4653345
+C         -5.22200        1.70985        0.08919
+O         -4.02946        1.97224       -0.59862
+C         -2.99480        1.99980        0.35432
+C         -1.67463        2.28017       -0.35418
+O         -0.63972        2.30623        0.59853
+C          0.55649        2.54191       -0.09259
+H         -5.53111        1.80817       -1.10080
+H         -5.05931        1.64530        1.22476
+H         -3.16613        2.81268        1.09884
+H         -2.90651        1.01331        0.86769
+H         -1.76067        3.26696       -0.86707
+H         -1.50561        1.46708       -1.09919
+H          0.53560        3.62886        0.28856
+H          0.63026        1.40254       -0.34908