Add files via upload

globus-labs · Sep 11, 2023 · 437f672 · 437f672
1 parent 8530c5f
commit 437f672
Show file tree

Hide file tree

Showing 18 changed files with 1,460 additions and 414 deletions.
diff --git a/mofa/simulation/cif2lammps/Dreiding_constants.py b/mofa/simulation/cif2lammps/Dreiding_constants.py
@@ -21,7 +21,8 @@
     "Si3": (0.937, 109.471, 4.27, 0.31, 0.0, 12.0),
     "P_3": (0.890, 93.3, 4.15, 0.32, 0.0, 12.0),
     "S_3": (1.040, 92.1, 4.03, 0.344, 0.0, 12.0),
-    "S_R": (1.040, 92.1, 4.03, 0.344, 0.0, 12.0),  # same as S_3, but R is more useful
+    # same as S_3, but R is more useful
+    "S_R": (1.040, 92.1, 4.03, 0.344, 0.0, 12.0),
     "Cl_": (0.997, 180.0, 3.9503, 0.2833, 0.0, 13.861),
     "Ga3": (1.210, 109.471, 4.39, 0.4, 0.0, 12.0),
     "Ge3": (1.210, 109.471, 4.27, 0.4, 0.0, 12.0),
@@ -37,9 +38,12 @@
     "Ca": (1.940, 90.0, 3.472, 0.05, 0.0, 12.0),
     "Fe": (1.285, 90.0, 4.54, 0.055, 0.0, 12.0),
     "Zn": (1.330, 109.471, 4.54, 0.055, 0.0, 12.0),
-    "Cu": (1.302, 90.0, 4.54, 0.055, 0.0, 12.0),  # R1 taken from UFF, the rest are just a copy of DREIDING Zn
-    "Zr": (1.564, 109.471, 4.54, 0.055, 0.0, 12.0),  # R1 taken from UFF, the rest are just a copy of DREIDING Zn
-    "Cr": (1.345, 90.0, 4.54, 0.055, 0.0, 12.0),  # R1 taken from UFF, the rest are just a copy of DREIDING Zn
+    # R1 taken from UFF, the rest are just a copy of DREIDING Zn
+    "Cu": (1.302, 90.0, 4.54, 0.055, 0.0, 12.0),
+    # R1 taken from UFF, the rest are just a copy of DREIDING Zn
+    "Zr": (1.564, 109.471, 4.54, 0.055, 0.0, 12.0),
+    # R1 taken from UFF, the rest are just a copy of DREIDING Zn
+    "Cr": (1.345, 90.0, 4.54, 0.055, 0.0, 12.0),
     # added for easier IDs
     "O_2_M": (0.560, 120.0, 3.4046, 0.0957, 0.0, 13.483),
     "O_3_M": (0.660, 104.51, 3.4046, 0.0957, 0.0, 13.483)

diff --git a/mofa/simulation/cif2lammps/Dreiding_construction.py b/mofa/simulation/cif2lammps/Dreiding_construction.py
@@ -28,7 +28,10 @@ def type_atoms(self):
             element_symbol = inf['element_symbol']
             nbors = list(SG.neighbors(name))
             nbor_symbols = [SG.nodes[n]['element_symbol'] for n in nbors]
-            bond_types = [SG.get_edge_data(name, n)['bond_type'] for n in nbors]
+            bond_types = [
+                SG.get_edge_data(
+                    name,
+                    n)['bond_type'] for n in nbors]
             mass = mass_key[element_symbol]
 
             # Atom typing for Dreiding, this can be made much more robust with pattern matching,
@@ -55,7 +58,8 @@ def type_atoms(self):
                     # oxygen case is complex with the UFF4MOF oxygen types
                     if element_symbol == 'O':
                         # -OH, for example
-                        if len(nbors) == 2 and 'A' not in bond_types and 'D' not in bond_types and not any(i in metals for i in nbor_symbols):
+                        if len(nbors) == 2 and 'A' not in bond_types and 'D' not in bond_types and not any(
+                                i in metals for i in nbor_symbols):
                             ty = 'O_3'
                             hyb = 'sp3'
                         # furan oxygen, for example
@@ -79,7 +83,10 @@ def type_atoms(self):
                             ty = 'O_3_M'
                             hyb = 'sp2'
                         else:
-                            raise ValueError('Oxygen with neighbors ' + ' '.join(nbor_symbols) + ' is not parametrized')
+                            raise ValueError(
+                                'Oxygen with neighbors ' +
+                                ' '.join(nbor_symbols) +
+                                ' is not parametrized')
                     # sulfur case is simple
                     elif element_symbol == 'S':
                         ty = 'S_' + str(len(nbors) + 1)
@@ -94,15 +101,20 @@ def type_atoms(self):
                     hyb = 'NA'
                 # if no type can be identified
                 else:
-                    raise ValueError('No Dreiding type identified for ' + element_symbol + 'with neighbors ' + ' '.join(nbor_symbols))
+                    raise ValueError(
+                        'No Dreiding type identified for ' +
+                        element_symbol +
+                        'with neighbors ' +
+                        ' '.join(nbor_symbols))
 
             types.append((ty, element_symbol, mass))
             SG.nodes[name]['force_field_type'] = ty
             SG.nodes[name]['hybridization'] = hyb
 
         types = set(types)
         Ntypes = len(types)
-        atom_types = dict((ty[0], i + 1) for i, ty in zip(range(Ntypes), types))
+        atom_types = dict((ty[0], i + 1)
+                          for i, ty in zip(range(Ntypes), types))
         atom_element_symbols = dict((ty[0], ty[1]) for ty in types)
         atom_masses = dict((ty[0], ty[2]) for ty in types)
 
@@ -196,7 +208,8 @@ def dihedral_parameters(self, bond, hybridization, element_symbols, nodes):
             n = 6.0
             V = 1.0
 
-            if (hyb_j == 'sp3' and els_j == 'O') or (hyb_k == 'sp3' and els_k == 'O'):
+            if (hyb_j == 'sp3' and els_j == 'O') or (
+                    hyb_k == 'sp3' and els_k == 'O'):
                 # case (i)
                 phi0 = 180.0
                 n = 2.0
@@ -276,7 +289,12 @@ def pair_parameters(self, charges=False):
             params[ID] = (style, eps_i, sig_i)
             comments[ID] = [a, a]
 
-        self.pair_data = {'params': params, 'style': style, 'special_bonds': sb, 'comments': comments, 'cutoff': cutoff}
+        self.pair_data = {
+            'params': params,
+            'style': style,
+            'special_bonds': sb,
+            'comments': comments,
+            'cutoff': cutoff}
 
     def enumerate_bonds(self):
 
@@ -291,7 +309,8 @@ def enumerate_bonds(self):
             fft_j = SG.nodes[j]['force_field_type']
             bond_type = data['bond_type']
 
-            # look for the bond order, otherwise use the convention based on the bond type
+            # look for the bond order, otherwise use the convention based on
+            # the bond type
             try:
                 bond_order = bond_order_dict[(fft_i, fft_j)]
             except KeyError:
@@ -327,7 +346,14 @@ def enumerate_bonds(self):
             all_bonds[ID] = bonds[b]
             count += len(bonds[b])
 
-        self.bond_data = {'all_bonds': all_bonds, 'params': bond_params, 'style': 'harmonic', 'count': (count, len(all_bonds)), 'comments': bond_comments}
+        self.bond_data = {
+            'all_bonds': all_bonds,
+            'params': bond_params,
+            'style': 'harmonic',
+            'count': (
+                count,
+                len(all_bonds)),
+            'comments': bond_comments}
 
     def enumerate_angles(self):
 
@@ -390,7 +416,14 @@ def enumerate_angles(self):
         else:
             style = 'hybrid ' + ' '.join(styles)
 
-        self.angle_data = {'all_angles': all_angles, 'params': angle_params, 'style': style, 'count': (count, len(all_angles)), 'comments': angle_comments}
+        self.angle_data = {
+            'all_angles': all_angles,
+            'params': angle_params,
+            'style': style,
+            'count': (
+                count,
+                len(all_angles)),
+            'comments': angle_comments}
 
     def enumerate_dihedrals(self):
 
@@ -422,8 +455,10 @@ def enumerate_dihedrals(self):
             element_symbols = (els_j, els_k)
 
             # here I calculate  parameters for each dihedral (I know) but I prefer identifying
-            # those dihedrals before passing to the final dihedral data construction.
-            params = self.dihedral_parameters(bond, hybridization, element_symbols, nodes)
+            # those dihedrals before passing to the final dihedral data
+            # construction.
+            params = self.dihedral_parameters(
+                bond, hybridization, element_symbols, nodes)
 
             if params != 'NA':
                 try:
@@ -444,11 +479,18 @@ def enumerate_dihedrals(self):
             params = dihedral_params[d]
             all_dihedrals[ID] = dihedrals[d]
             indexed_dihedral_params[ID] = list(dihedral_params[d])
-            dihedral_comments[ID] = list(dihedral) + ['bond order=' + str(d[2])]
+            dihedral_comments[ID] = list(
+                dihedral) + ['bond order=' + str(d[2])]
             count += len(dihedrals[d])
 
-        self.dihedral_data = {'all_dihedrals': all_dihedrals, 'params': indexed_dihedral_params,
-                              'style': 'charmm', 'count': (count, len(all_dihedrals)), 'comments': dihedral_comments}
+        self.dihedral_data = {
+            'all_dihedrals': all_dihedrals,
+            'params': indexed_dihedral_params,
+            'style': 'charmm',
+            'count': (
+                count,
+                len(all_dihedrals)),
+            'comments': dihedral_comments}
 
     def enumerate_impropers(self):
 
@@ -491,8 +533,14 @@ def enumerate_impropers(self):
                 all_impropers[ID] = impropers[i]
                 count += len(impropers[i])
 
-        self.improper_data = {'all_impropers': all_impropers, 'params': improper_params,
-                              'style': 'umbrella', 'count': (count, len(all_impropers)), 'comments': improper_comments}
+        self.improper_data = {
+            'all_impropers': all_impropers,
+            'params': improper_params,
+            'style': 'umbrella',
+            'count': (
+                count,
+                len(all_impropers)),
+            'comments': improper_comments}
 
     def compile_force_field(self, charges):
 

diff --git a/mofa/simulation/cif2lammps/UFF4MOF_constants.py b/mofa/simulation/cif2lammps/UFF4MOF_constants.py
@@ -230,8 +230,10 @@
     "No6+3": (1.679, 90.0000, 3.248, 0.011, 12.0, 3.8882, 0.000, 3.4750),
     "Lr6+3": (1.698, 90.0000, 3.236, 0.011, 12.0, 3.8882, 0.000, 3.5000),
     # # Special O-types for using alternate bond orders
-    "O_2_M": (0.634, 120.00, 3.500, 0.060, 14.085, 2.3000, 2.000, 8.7410),  # carboxyllic oxygen bonded to metal
-    "O_3_M": (0.658, 104.51, 3.500, 0.060, 14.085, 2.3000, 0.018, 8.7410)  # this should be used for O in coordinated solvent
+    # carboxyllic oxygen bonded to metal
+    "O_2_M": (0.634, 120.00, 3.500, 0.060, 14.085, 2.3000, 2.000, 8.7410),
+    # this should be used for O in coordinated solvent
+    "O_3_M": (0.658, 104.51, 3.500, 0.060, 14.085, 2.3000, 0.018, 8.7410)
 }
 
 # all M-M bonds have order 0.25 unless otherwise specified by UFF4MOF