From 437f6727b6937b9b91fd5cfb5e9c98162f691d92 Mon Sep 17 00:00:00 2001 From: doublylinkedlist Date: Mon, 11 Sep 2023 03:50:07 -0500 Subject: [PATCH] Add files via upload --- .../cif2lammps/Dreiding_constants.py | 12 +- .../cif2lammps/Dreiding_construction.py | 82 +++- .../cif2lammps/UFF4MOF_constants.py | 6 +- .../cif2lammps/UFF4MOF_construction.py | 212 +++++++--- mofa/simulation/cif2lammps/UFF_constants.py | 9 +- .../simulation/cif2lammps/UFF_construction.py | 115 +++-- .../cif2lammps/ZIFFF_construction.py | 201 ++++++--- mofa/simulation/cif2lammps/cif2system.py | 393 ++++++++++++++---- mofa/simulation/cif2lammps/gaff.py | 6 +- mofa/simulation/cif2lammps/gaff2.py | 3 +- mofa/simulation/cif2lammps/main_conversion.py | 223 ++++++++-- .../cif2lammps/pymatgen_cif2system.py | 197 +++++++-- .../cif2lammps/small_molecule_constants.py | 19 +- .../cif2lammps/small_molecule_construction.py | 148 +++++-- .../cif2lammps/write_GULP_inputs.py | 21 +- .../cif2lammps/write_lammps_data.py | 173 +++++--- .../cif2lammps/write_molecule_files.py | 18 +- .../cif2lammps/zeoliteFFs_construction.py | 36 +- 18 files changed, 1460 insertions(+), 414 deletions(-) diff --git a/mofa/simulation/cif2lammps/Dreiding_constants.py b/mofa/simulation/cif2lammps/Dreiding_constants.py index d734236a..b86e35be 100644 --- 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 index 9604eb8e..59f19124 100644 --- 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,7 +101,11 @@ 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 @@ -102,7 +113,8 @@ def type_atoms(self): 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 index 0d0069f7..5932b3fc 100644 --- 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 diff --git a/mofa/simulation/cif2lammps/UFF4MOF_construction.py b/mofa/simulation/cif2lammps/UFF4MOF_construction.py index 45f628f3..eafa0368 100644 --- a/mofa/simulation/cif2lammps/UFF4MOF_construction.py +++ b/mofa/simulation/cif2lammps/UFF4MOF_construction.py @@ -97,7 +97,8 @@ def __init__(self, system, cutoff, args): cx = c * np.cos(beta * pi / 180.0) cy = (c * b * np.cos(alpha * pi / 180.0) - bx * cx) / by cz = (c ** 2.0 - cx ** 2.0 - cy ** 2.0) ** 0.5 - self.unit_cell = np.asarray([[ax, ay, az], [bx, by, bz], [cx, cy, cz]]).T + self.unit_cell = np.asarray( + [[ax, ay, az], [bx, by, bz], [cx, cy, cz]]).T def type_atoms(self): @@ -113,7 +114,10 @@ def type_atoms(self): element_symbol = inf['element_symbol'] nbors = [a for a in 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] if len(nbors) > 1: @@ -126,26 +130,32 @@ def type_atoms(self): comp_coords = [] for n in nbors: - dist_n, sym_n = PBC3DF_sym(SG.nodes[n]['fractional_position'], inf['fractional_position']) + dist_n, sym_n = PBC3DF_sym( + SG.nodes[n]['fractional_position'], inf['fractional_position']) dist_n = np.dot(self.unit_cell, dist_n) fcoord = SG.nodes[n]['fractional_position'] + sym_n comp_coords.append(np.dot(self.unit_cell, fcoord)) comp_coords = np.array(comp_coords) - dist_square = superimpose(comp_coords, comparison_geometries['square_planar'], count) - dist_tetrahedral = superimpose(comp_coords, comparison_geometries['tetrahedral'], count) + dist_square = superimpose( + comp_coords, comparison_geometries['square_planar'], count) + dist_tetrahedral = superimpose( + comp_coords, comparison_geometries['tetrahedral'], count) else: angles = [] for n0, n1 in itertools.combinations(nbors, 2): - dist_j, sym_j = PBC3DF_sym(SG.nodes[n0]['fractional_position'], inf['fractional_position']) - dist_k, sym_k = PBC3DF_sym(SG.nodes[n1]['fractional_position'], inf['fractional_position']) + dist_j, sym_j = PBC3DF_sym( + SG.nodes[n0]['fractional_position'], inf['fractional_position']) + dist_k, sym_k = PBC3DF_sym( + SG.nodes[n1]['fractional_position'], inf['fractional_position']) dist_j = np.dot(self.unit_cell, dist_j) - cosine_angle = np.dot(dist_j, dist_k) / (np.linalg.norm(dist_j) * np.linalg.norm(dist_k)) + cosine_angle = np.dot( + dist_j, dist_k) / (np.linalg.norm(dist_j) * np.linalg.norm(dist_k)) if cosine_angle > 1: cosine_angle = 1 @@ -159,11 +169,13 @@ def type_atoms(self): angles = np.array(angles) dist_linear = min(np.abs(angles - 180.0)) dist_triangle = min(np.abs(angles - 120.0)) - dist_square = min(min(np.abs(angles - 90.0)), min(np.abs(angles - 180.0))) + dist_square = min(min(np.abs(angles - 90.0)), + min(np.abs(angles - 180.0))) dist_corner = min(np.abs(angles - 90.0)) dist_tetrahedral = min(np.abs(angles - 109.47)) - if 'A' in bond_types and element_symbol not in ('O', 'H') and element_symbol not in metals: + if 'A' in bond_types and element_symbol not in ( + 'O', 'H') and element_symbol not in metals: ty = element_symbol + '_' + 'R' hyb = 'resonant' else: @@ -208,7 +220,8 @@ def type_atoms(self): elif 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' - # coordinated solvent, same parameters as O_3, but different name to modulate bond orders + # coordinated solvent, same parameters as O_3, but + # different name to modulate bond orders elif len(nbors) in (2, 3) and len([i for i in nbor_symbols if i in metals]) == 1 and 'H' in nbor_symbols: ty = 'O_3_M' hyb = 'sp3' @@ -220,7 +233,9 @@ def type_atoms(self): elif len(nbors) == 2 and 'D' in bond_types and not any(i in metals for i in nbor_symbols): ty = 'O_2' hyb = 'sp2' - # carboxylate oxygen bound to metal node, same parameters as O_2, but different name to modulate bond orders + # carboxylate oxygen bound to metal node, same + # parameters as O_2, but different name to modulate + # bond orders elif len(nbors) == 2 and any(i in metals for i in nbor_symbols) and 'C' in nbor_symbols: ty = 'O_2_M' hyb = 'sp2' @@ -234,10 +249,13 @@ def type_atoms(self): # 3 connected oxygens bonded to metals elif len(nbors) == 3 and len([i for i in nbor_symbols if i in metals]) > 1: # trigonal geometry - if dist_triangle < dist_tetrahedral and not any(i in ['Zr', 'Eu', 'Tb', 'U'] for i in nbor_symbols): + if dist_triangle < dist_tetrahedral and not any( + i in ['Zr', 'Eu', 'Tb', 'U'] for i in nbor_symbols): ty = 'O_2_z' hyb = 'sp2' - # sometimes oxygens in Zr6 and analagous nodes can have near trigonal geometry, still want O_3_f, however + # sometimes oxygens in Zr6 and analagous nodes can + # have near trigonal geometry, still want O_3_f, + # however elif dist_triangle < dist_tetrahedral and any(i in ['Zr', 'Eu', 'Tb', 'U'] for i in nbor_symbols): ty = 'O_3_f' hyb = 'sp2' @@ -299,60 +317,90 @@ def type_atoms(self): # 2 connected, linear if len(nbors) == 2 and dist_linear < 60.0: - options = ('1f1', '4f2', '4+2', '6f3', '6+3', '6+2', '6+4') - ty = typing_loop(options, add_symbol, UFF4MOF_atom_parameters) + options = ( + '1f1', '4f2', '4+2', '6f3', '6+3', '6+2', '6+4') + ty = typing_loop( + options, add_symbol, UFF4MOF_atom_parameters) # incomplete square planar elif len(nbors) == 3 and dist_square < min(dist_tetrahedral, dist_triangle): options = ('4f2', '4+2') - ty = typing_loop(options, add_symbol, UFF4MOF_atom_parameters) + ty = typing_loop( + options, add_symbol, UFF4MOF_atom_parameters) # incomplete tetrahedron elif len(nbors) == 3 and dist_tetrahedral < min(dist_square, dist_triangle): options = ('3f2', '3+2') - ty = typing_loop(options, add_symbol, UFF4MOF_atom_parameters) + ty = typing_loop( + options, add_symbol, UFF4MOF_atom_parameters) # trigonal, only Cu, Zn, Ag elif len(nbors) == 3 and dist_triangle < min(dist_square, dist_tetrahedral): options = ['2f2'] - ty = typing_loop(options, add_symbol, UFF4MOF_atom_parameters) + ty = typing_loop( + options, add_symbol, UFF4MOF_atom_parameters) # 4 connected, square planar elif len(nbors) == 4 and dist_square < dist_tetrahedral: options = ('4f2', '4+2', '6f3', '6+3', '6+2', '6+4') - ty = typing_loop(options, add_symbol, UFF4MOF_atom_parameters) + ty = typing_loop( + options, add_symbol, UFF4MOF_atom_parameters) - # 4 connected, tetrahedral 3+3 does not apply for As, Sb, Tl, Bi + # 4 connected, tetrahedral 3+3 does not apply for As, Sb, + # Tl, Bi elif len(nbors) == 4 and (dist_tetrahedral < dist_square): - options = ('3f2', '3f4', '3+1', '3+2', '3+3', '3+4', '3+5', '3+6') - ty = typing_loop(options, add_symbol, UFF4MOF_atom_parameters) + options = ('3f2', '3f4', '3+1', '3+2', + '3+3', '3+4', '3+5', '3+6') + ty = typing_loop( + options, add_symbol, UFF4MOF_atom_parameters) - # paddlewheels, 5 neighbors if bare, 6 neighbors if pillared, one should be another metal + # paddlewheels, 5 neighbors if bare, 6 neighbors if + # pillared, one should be another metal elif len(nbors) in (5, 6) and any(i in metals for i in nbor_symbols): if (dist_square < dist_tetrahedral): - options = ('4f2', '4+2', '6f3', '6+3', '6+2', '6+4') - ty = typing_loop(options, add_symbol, UFF4MOF_atom_parameters) + options = ( + '4f2', '4+2', '6f3', '6+3', '6+2', '6+4') + ty = typing_loop( + options, add_symbol, UFF4MOF_atom_parameters) else: - options = ('4f2', '4+2', '6f3', '6+3', '6+2', '6+4') - ty = typing_loop(options, add_symbol, UFF4MOF_atom_parameters) - message = 'There is a ' + element_symbol + ' that has a near tetrahedral angle typed as ' + ty + '\n' - message += 'The neighbors are ' + ' '.join(nbor_symbols) + options = ( + '4f2', '4+2', '6f3', '6+3', '6+2', '6+4') + ty = typing_loop( + options, add_symbol, UFF4MOF_atom_parameters) + message = 'There is a ' + element_symbol + \ + ' that has a near tetrahedral angle typed as ' + ty + '\n' + message += 'The neighbors are ' + \ + ' '.join(nbor_symbols) warnings.warn(message) - # M3O(CO2H)6 metals, e.g. MIL-100, paddlewheel options are secondary, followed by 8f4 (should give nearly correct geometry) + # M3O(CO2H)6 metals, e.g. MIL-100, paddlewheel options are + # secondary, followed by 8f4 (should give nearly correct + # geometry) elif len(nbors) in (5, 6) and not any(i in metals for i in nbor_symbols) and (dist_square < dist_tetrahedral): options = ('6f3', '6+2', '6+3', '6+4', '4f2', '4+2') - ty = typing_loop(options, add_symbol, UFF4MOF_atom_parameters) + ty = typing_loop( + options, add_symbol, UFF4MOF_atom_parameters) # 5,6 connected, with near-tetrahedral angles elif len(nbors) in (5, 6) and not any(i in metals for i in nbor_symbols) and (dist_tetrahedral < dist_square): - options = ('8f4', '3f2', '3f4', '3+1', '3+2', '3+3', '3+4', '3+5', '3+6') - ty = typing_loop(options, add_symbol, UFF4MOF_atom_parameters) + options = ( + '8f4', + '3f2', + '3f4', + '3+1', + '3+2', + '3+3', + '3+4', + '3+5', + '3+6') + ty = typing_loop( + options, add_symbol, UFF4MOF_atom_parameters) # highly connected metals (max of 12 neighbors) elif 7 <= len(nbors) <= 14: options = ['8f4'] - ty = typing_loop(options, add_symbol, UFF4MOF_atom_parameters) + ty = typing_loop( + options, add_symbol, UFF4MOF_atom_parameters) # only one type for Bi elif element_symbol in ('As', 'Bi', 'Tl', 'Sb'): @@ -379,14 +427,23 @@ def type_atoms(self): SG.nodes[name]['hybridization'] = hyb if ty is None: - raise ValueError('No UFF4MOF type identified for atom ' + element_symbol + ' with neighbors ' + ' '.join(nbor_symbols)) + raise ValueError( + 'No UFF4MOF type identified for atom ' + + element_symbol + + ' with neighbors ' + + ' '.join(nbor_symbols)) elif ty == 'C_R' and len(nbor_symbols) > 3: - raise ValueError('Too many neighbors for aromatic carbon ' + element_symbol + ' with neighbors ' + ' '.join(nbor_symbols)) + raise ValueError( + 'Too many neighbors for aromatic carbon ' + + element_symbol + + ' with neighbors ' + + ' '.join(nbor_symbols)) types = set(types) types = sorted(types, key=lambda x: x[0]) 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) @@ -412,7 +469,8 @@ def bond_parameters(self, bond, bond_order): # bond-order correction rbo = -0.1332 * (r0_i + r0_j) * np.log(bond_order) # electronegativity correction - ren = r0_i * r0_j * (((np.sqrt(X_i) - np.sqrt(X_j))**2)) / (X_i * r0_i + X_j * r0_j) + ren = r0_i * r0_j * (((np.sqrt(X_i) - np.sqrt(X_j))**2) + ) / (X_i * r0_i + X_j * r0_j) # equilibrium distance r_ij = r0_i + r0_j + rbo - ren r_ij3 = r_ij * r_ij * r_ij @@ -458,7 +516,8 @@ def angle_parameters(self, angle, r_ij, r_jk): r_ik = np.sqrt(r_ij**2.0 + r_jk**2.0 - 2.0 * r_ij * r_jk * cosT0) # force constant - K = ((664.12 * Z1_i * Z1_k) / (r_ik**5.0)) * (3.0 * r_ij * r_jk * (1.0 - cosT0**2.0) - r_ik**2.0 * cosT0) + K = ((664.12 * Z1_i * Z1_k) / (r_ik**5.0)) * \ + (3.0 * r_ij * r_jk * (1.0 - cosT0**2.0) - r_ik**2.0 * cosT0) # general non-linear if theta0_j not in (90.0, 120.0, 180.0): @@ -473,7 +532,8 @@ def angle_parameters(self, angle, r_ij, r_jk): # the 1/2 scaling is needed to correct the LAMMPS angle energy calculation # angle energy for angle_style cosine/periodic is multiplied by 2 in LAMMPS for some reason # see https://github.com/lammps/lammps/blob/master/src/MOLECULE/angle_cosine_periodic.cpp, line 140 - # the 1/2 factor for the UFF fourier angles should be included here as it is not included in LAMMPS + # the 1/2 factor for the UFF fourier angles should be included here as + # it is not included in LAMMPS K *= 0.5 return (angle_style, K, b, n) @@ -600,7 +660,11 @@ def pair_parameters(self, charges=False): params[ID] = (style, D_i, x_i) comments[ID] = [a, a] - self.pair_data = {'params': params, 'style': style, 'special_bonds': sb, 'comments': comments} + self.pair_data = { + 'params': params, + 'style': style, + 'special_bonds': sb, + 'comments': comments} def enumerate_bonds(self): @@ -617,7 +681,8 @@ def enumerate_bonds(self): esi = SG.nodes[i]['element_symbol'] esj = SG.nodes[j]['element_symbol'] - # 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: @@ -625,7 +690,13 @@ def enumerate_bonds(self): bond_order = bond_order_dict[(fft_j, fft_i)] except KeyError: # half for metal-nonmetal - if any(a in metals for a in (esi, esj)) and not all(a in metals for a in (esi, esj)): + if any( + a in metals for a in ( + esi, + esj)) and not all( + a in metals for a in ( + esi, + esj)): bond_order = 0.5 # quarter for metal-metal elif all(a in metals for a in (esi, esj)): @@ -669,7 +740,14 @@ def enumerate_bonds(self): else: style = 'hybrid ' + ' '.join(styles) - self.bond_data = {'all_bonds': all_bonds, 'params': bond_params, 'style': style, 'count': (count, len(all_bonds)), 'comments': bond_comments} + self.bond_data = { + 'all_bonds': all_bonds, + 'params': bond_params, + 'style': style, + 'count': ( + count, + len(all_bonds)), + 'comments': bond_comments} def enumerate_angles(self): @@ -751,7 +829,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): @@ -783,8 +868,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: @@ -805,11 +892,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': 'harmonic', 'count': (count, len(all_dihedrals)), 'comments': dihedral_comments} + self.dihedral_data = { + 'all_dihedrals': all_dihedrals, + 'params': indexed_dihedral_params, + 'style': 'harmonic', + 'count': ( + count, + len(all_dihedrals)), + 'comments': dihedral_comments} def enumerate_impropers(self): @@ -824,7 +918,8 @@ def enumerate_impropers(self): if len(nbors) == 3: fft_i = data['force_field_type'] - fft_nbors = tuple(sorted([SG.nodes[m]['force_field_type'] for m in nbors])) + fft_nbors = tuple( + sorted([SG.nodes[m]['force_field_type'] for m in nbors])) O_2_flag = False # force constant is much larger if j,k, or l is O_2 if 'O_2' in fft_nbors or 'O_2_M' in fft_nbors: @@ -853,12 +948,19 @@ def enumerate_impropers(self): if params is not None: ID += 1 improper_params[ID] = list(params) - improper_comments[ID] = [i[0], 'X', 'X', 'X', 'O_2 present=' + str(O_2_flag)] + improper_comments[ID] = [ + i[0], 'X', 'X', 'X', 'O_2 present=' + str(O_2_flag)] all_impropers[ID] = impropers[i] count += len(impropers[i]) - self.improper_data = {'all_impropers': all_impropers, 'params': improper_params, - 'style': 'fourier', 'count': (count, len(all_impropers)), 'comments': improper_comments} + self.improper_data = { + 'all_impropers': all_impropers, + 'params': improper_params, + 'style': 'fourier', + 'count': ( + count, + len(all_impropers)), + 'comments': improper_comments} def compile_force_field(self, charges=False): diff --git a/mofa/simulation/cif2lammps/UFF_constants.py b/mofa/simulation/cif2lammps/UFF_constants.py index 78f2e427..fd1ff504 100644 --- a/mofa/simulation/cif2lammps/UFF_constants.py +++ b/mofa/simulation/cif2lammps/UFF_constants.py @@ -56,7 +56,8 @@ "Co6+2": (1.241, 90.0, 2.872, 0.014, 12.0, 2.4301, 0.000, 4.105), "Ni4+2": (1.164, 90.0, 2.834, 0.015, 12.0, 2.4301, 0.000, 4.465), "Cu3+1": (1.302, 109.4712, 3.495, 0.005, 12.0, 1.7565, 0.000, 3.729), - "Cu4+1": (1.302, 90.0, 3.495, 0.005, 12.0, 1.7565, 0.000, 3.729), # justed changed to angle to 90.0 for paddlewheels + # justed changed to angle to 90.0 for paddlewheels + "Cu4+1": (1.302, 90.0, 3.495, 0.005, 12.0, 1.7565, 0.000, 3.729), "Zn3+2": (1.193, 109.4712, 2.763, 0.124, 12.0, 1.3084, 0.000, 5.106), "Ga3+3": (1.260, 109.4712, 4.383, 0.415, 11.0, 1.8206, 0.000, 2.999), "Ge3": (1.197, 109.4712, 4.280, 0.379, 12.0, 2.7888, 0.701, 4.051), @@ -136,8 +137,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.0, 3.500, 0.060, 14.085, 2.2998, 2.000, 8.741), # carboxyllic oxygen bonded to metal - "O_3_M": (0.658, 104.51, 3.500, 0.060, 14.085, 2.2998, 0.018, 8.741) # this should be used for O in coordinated solvent + # carboxyllic oxygen bonded to metal + "O_2_M": (0.634, 120.0, 3.500, 0.060, 14.085, 2.2998, 2.000, 8.741), + # this should be used for O in coordinated solvent + "O_3_M": (0.658, 104.51, 3.500, 0.060, 14.085, 2.2998, 0.018, 8.741) } # without MIL-53 type nodes diff --git a/mofa/simulation/cif2lammps/UFF_construction.py b/mofa/simulation/cif2lammps/UFF_construction.py index deba3b54..0d8a6da7 100644 --- a/mofa/simulation/cif2lammps/UFF_construction.py +++ b/mofa/simulation/cif2lammps/UFF_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 UFF, this can be made much more robust with pattern matching, @@ -89,7 +92,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) @@ -103,8 +109,10 @@ def type_atoms(self): hyb = 'sp1' # Metals elif element_symbol in metals: - # Cu paddlewheel, just changed equilibrium angle of Cu3+1 to 90.0 - if len(nbors) == 5 and element_symbol == 'Cu' and any(i in metals for i in nbor_symbols): + # Cu paddlewheel, just changed equilibrium angle of Cu3+1 + # to 90.0 + if len(nbors) == 5 and element_symbol == 'Cu' and any( + i in metals for i in nbor_symbols): ty = element_symbol + '4+1' hyb = 'NA' # M3O(CO2H)6 metals, e.g. MIL-100 @@ -123,7 +131,11 @@ def type_atoms(self): hyb = 'NA' # if no type can be identified else: - raise ValueError('No UFF type identified for ' + element_symbol + 'with neighbors ' + ' '.join(nbor_symbols)) + raise ValueError( + 'No UFF type identified for ' + + element_symbol + + 'with neighbors ' + + ' '.join(nbor_symbols)) types.append((ty, element_symbol, mass)) SG.nodes[name]['force_field_type'] = ty @@ -131,7 +143,8 @@ def type_atoms(self): 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) @@ -155,7 +168,8 @@ def bond_parameters(self, bond, bond_order): # bond-order correction rbo = -0.1332 * (r0_i + r0_j) * np.log(bond_order) # electronegativity correction - ren = r0_i * r0_j * (((np.sqrt(X_i) - np.sqrt(X_j))**2)) / (X_i * r0_i + X_j * r0_j) + ren = r0_i * r0_j * (((np.sqrt(X_i) - np.sqrt(X_j))**2) + ) / (X_i * r0_i + X_j * r0_j) # equilibrium distance r_ij = r0_i + r0_j + rbo - ren r_ij3 = r_ij * r_ij * r_ij @@ -201,7 +215,8 @@ def angle_parameters(self, angle, r_ij, r_jk): r_ik = np.sqrt(r_ij**2.0 + r_jk**2.0 - 2.0 * r_ij * r_jk * cosT0) # force constant - K = ((664.12 * Z1_i * Z1_k) / (r_ik**5.0)) * (3.0 * r_ij * r_jk * (1.0 - cosT0**2.0) - r_ik**2.0 * cosT0) + K = ((664.12 * Z1_i * Z1_k) / (r_ik**5.0)) * \ + (3.0 * r_ij * r_jk * (1.0 - cosT0**2.0) - r_ik**2.0 * cosT0) # general non-linear if theta0_j not in (90.0, 120.0, 180.0): @@ -342,7 +357,12 @@ def pair_parameters(self, charges=False): params[ID] = (style, D_i, x_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): @@ -357,7 +377,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: @@ -393,7 +414,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): @@ -417,9 +445,21 @@ def enumerate_angles(self): fft_j = SG.nodes[j]['force_field_type'] fft_k = SG.nodes[k]['force_field_type'] - octa_metals = ('Al6+3', 'Sc6+3', 'Ti4+2', 'V_4+2', 'V_6+3', 'Cr4+2', - 'Cr6f3', 'Mn6+3', 'Mn4+2', 'Fe6+3', 'Fe4+2', 'Co4+2', - 'Cu4+2', 'Zn4+2') + octa_metals = ( + 'Al6+3', + 'Sc6+3', + 'Ti4+2', + 'V_4+2', + 'V_6+3', + 'Cr4+2', + 'Cr6f3', + 'Mn6+3', + 'Mn4+2', + 'Fe6+3', + 'Fe4+2', + 'Co4+2', + 'Cu4+2', + 'Zn4+2') if fft_j in octa_metals: i_coord = SG.nodes[i]['cartesian_position'] @@ -427,7 +467,8 @@ def enumerate_angles(self): k_coord = SG.nodes[k]['cartesian_position'] ij = i_coord - j_coord jk = j_coord - k_coord - cosine_angle = np.dot(ij, jk) / (np.linalg.norm(ij) * np.linalg.norm(jk)) + cosine_angle = np.dot( + ij, jk) / (np.linalg.norm(ij) * np.linalg.norm(jk)) angle = (180.0 / np.pi) * np.arccos(cosine_angle) sort_ik = sorted([(fft_i, i), (fft_k, k)], key=lambda x: x[0]) @@ -482,7 +523,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): @@ -514,8 +562,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: @@ -536,11 +586,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': 'harmonic', 'count': (count, len(all_dihedrals)), 'comments': dihedral_comments} + self.dihedral_data = { + 'all_dihedrals': all_dihedrals, + 'params': indexed_dihedral_params, + 'style': 'harmonic', + 'count': ( + count, + len(all_dihedrals)), + 'comments': dihedral_comments} def enumerate_impropers(self): @@ -555,7 +612,8 @@ def enumerate_impropers(self): if len(nbors) == 3: fft_i = data['force_field_type'] - fft_nbors = tuple(sorted([SG.nodes[m]['force_field_type'] for m in nbors])) + fft_nbors = tuple( + sorted([SG.nodes[m]['force_field_type'] for m in nbors])) O_2_flag = False # force constant is much larger if j,k, or l is O_2 if 'O_2' in fft_nbors or 'O_2_M' in fft_nbors: @@ -584,12 +642,19 @@ def enumerate_impropers(self): if params is not None: ID += 1 improper_params[ID] = list(params) - improper_comments[ID] = [i[0], 'X', 'X', 'X', 'O_2 present=' + str(O_2_flag)] + improper_comments[ID] = [ + i[0], 'X', 'X', 'X', 'O_2 present=' + str(O_2_flag)] all_impropers[ID] = impropers[i] count += len(impropers[i]) - self.improper_data = {'all_impropers': all_impropers, 'params': improper_params, - 'style': 'fourier', 'count': (count, len(all_impropers)), 'comments': improper_comments} + self.improper_data = { + 'all_impropers': all_impropers, + 'params': improper_params, + 'style': 'fourier', + 'count': ( + count, + len(all_impropers)), + 'comments': improper_comments} def compile_force_field(self, charges=False): diff --git a/mofa/simulation/cif2lammps/ZIFFF_construction.py b/mofa/simulation/cif2lammps/ZIFFF_construction.py index 57fbffc2..d4c2ffb0 100644 --- a/mofa/simulation/cif2lammps/ZIFFF_construction.py +++ b/mofa/simulation/cif2lammps/ZIFFF_construction.py @@ -75,40 +75,46 @@ def read_gaffdat(mode='gaff'): else: raise ValueError('mode must be gaff or gaff2') - gaff_atom_types = [LL.split() for LL in gaff_atom_types.split('\n') if len(LL.split()) > 0] - gaff_atom_types = dict((LL[0], (float(LL[1]), float(LL[2]))) for LL in gaff_atom_types) - - gaff_LJ_parameters = [LL.split() for LL in gaff_LJ_parameters.split('\n') if len(LL.split()) > 0] - gaff_LJ_parameters = dict((LL[0], (float(LL[2]), float(LL[1]))) for LL in gaff_LJ_parameters) - - gaff_bond_list = [(''.join(LL[0:5].split()), - ''.join(LL[5:16].split()), - ''.join(LL[16:26].split())) for LL in gaff_bonds.split('\n') if len(LL.split()) > 0] + gaff_atom_types = [ + LL.split() for LL in gaff_atom_types.split('\n') if len( + LL.split()) > 0] + gaff_atom_types = dict((LL[0], (float(LL[1]), float(LL[2]))) + for LL in gaff_atom_types) + + gaff_LJ_parameters = [ + LL.split() for LL in gaff_LJ_parameters.split('\n') if len( + LL.split()) > 0] + gaff_LJ_parameters = dict( + (LL[0], (float(LL[2]), float(LL[1]))) for LL in gaff_LJ_parameters) + + gaff_bond_list = [(''.join(LL[0:5].split()), ''.join(LL[5:16].split()), ''.join( + LL[16:26].split())) for LL in gaff_bonds.split('\n') if len(LL.split()) > 0] gaff_bonds = {} for LL in gaff_bond_list: bond, K, r0 = LL bond = tuple(sorted(bond.split('-'))) gaff_bonds[bond] = (float(K), float(r0)) - gaff_angle_list = [(''.join(LL[0:8].split()), - ''.join(LL[8:20].split()), - ''.join(LL[20:30].split())) for LL in gaff_angles.split('\n') if len(LL.split()) > 0] + gaff_angle_list = [(''.join(LL[0:8].split()), ''.join(LL[8:20].split()), ''.join( + LL[20:30].split())) for LL in gaff_angles.split('\n') if len(LL.split()) > 0] gaff_angles = {} for LL in gaff_angle_list: angle, K, theta0 = LL angle = tuple(angle.split('-')) gaff_angles[angle] = (float(K), float(r0)) - gaff_dihedral_list = [(''.join(LL[0:11].split()), ''.join(LL[11:19].split()), ''.join(LL[19:31].split()), - ''.join(LL[31:48].split()), ''.join(LL[48:52].split())) for LL in gaff_dihedrals.split('\n') if len(LL.split()) > 0] + gaff_dihedral_list = [(''.join(LL[0:11].split()), ''.join(LL[11:19].split()), ''.join(LL[19:31].split()), ''.join( + LL[31:48].split()), ''.join(LL[48:52].split())) for LL in gaff_dihedrals.split('\n') if len(LL.split()) > 0] gaff_dihedrals = {} for LL in gaff_dihedral_list: dihedral, m, K, psi0, n = LL dihedral = tuple(dihedral.split('-')) try: - gaff_dihedrals[dihedral].extend([float(K) / int(float(m)), int(float(n)), float(psi0)]) + gaff_dihedrals[dihedral].extend( + [float(K) / int(float(m)), int(float(n)), float(psi0)]) except KeyError: - gaff_dihedrals[dihedral] = [float(K) / int(float(m)), int(float(n)), float(psi0)] + gaff_dihedrals[dihedral] = [ + float(K) / int(float(m)), int(float(n)), float(psi0)] gaff_improper_list = [(''.join(LL[0:11].split()), ''.join(LL[11:33].split()), ''.join(LL[33:47].split()), ''.join(LL[47:50].split())) for LL in gaff_impropers.split('\n') if len(LL.split()) > 0] @@ -118,8 +124,13 @@ def read_gaffdat(mode='gaff'): improper = tuple(improper.split('-')) gaff_impropers[improper] = (float(K), int(float(n)), float(psi0)) - gaff_data = {'types': gaff_atom_types, 'bonds': gaff_bonds, 'angles': gaff_angles, - 'dihedrals': gaff_dihedrals, 'impropers': gaff_impropers, 'LJ_parameters': gaff_LJ_parameters} + gaff_data = { + 'types': gaff_atom_types, + 'bonds': gaff_bonds, + 'angles': gaff_angles, + 'dihedrals': gaff_dihedrals, + 'impropers': gaff_impropers, + 'LJ_parameters': gaff_LJ_parameters} return gaff_data @@ -157,7 +168,8 @@ def GAFF_type(atom, data, SG, pure_aromatic_atoms, aromatic_atoms): elif len(nbors) == 3: # aromatic carbons if 'A' in bond_types: - # pure aromatic C are in benzene/pyridine (https://github.com/rsdefever/GAFF-foyer) + # pure aromatic C are in benzene/pyridine + # (https://github.com/rsdefever/GAFF-foyer) if atom in pure_aromatic_atoms: ty = 'ca' # non-pure aromatic C (anything else with aromatic bond) @@ -182,11 +194,13 @@ def GAFF_type(atom, data, SG, pure_aromatic_atoms, aromatic_atoms): ty = 'c3' hyb = 'sp3' - # note that cp, cq, cd, ce, cf, cg, ch, cx, cy, cu, cv, cz are are probably not needed for ZIF-FF + # note that cp, cq, cd, ce, cf, cg, ch, cx, cy, cu, cv, cz are are + # probably not needed for ZIF-FF elif sym == 'H': - # make sure to apply this function to hydrogens after all other atoms have been typed + # make sure to apply this function to hydrogens after all other atoms + # have been typed if len(nbors) != 1: raise ValueError('H with more than one neighbor found') @@ -204,7 +218,8 @@ def GAFF_type(atom, data, SG, pure_aromatic_atoms, aromatic_atoms): if nbor_hyb == 'sp3': ty = 'h1' hyb = 'sp1' - # bonded to non-aromatic sp2 carbon, need to add electron withdrawing cases + # bonded to non-aromatic sp2 carbon, need to add electron + # withdrawing cases elif nbor_hyb == 'sp2': ty = 'h4' hyb = 'sp1' @@ -294,7 +309,11 @@ def GAFF_type(atom, data, SG, pure_aromatic_atoms, aromatic_atoms): pass if (ty is None or hyb is None) and sym != 'H': - raise ValueError('No GAFF atom type identified for element', sym, 'with neighbors', nbor_symbols) + raise ValueError( + 'No GAFF atom type identified for element', + sym, + 'with neighbors', + nbor_symbols) return ty, hyb @@ -314,7 +333,8 @@ def __init__(self, system, cutoff, args): cx = c * np.cos(beta * pi / 180.0) cy = (c * b * np.cos(alpha * pi / 180.0) - bx * cx) / by cz = (c ** 2.0 - cx ** 2.0 - cy ** 2.0) ** 0.5 - self.unit_cell = np.asarray([[ax, ay, az], [bx, by, bz], [cx, cy, cz]]).T + self.unit_cell = np.asarray( + [[ax, ay, az], [bx, by, bz], [cx, cy, cz]]).T NMG = system['graph'].copy() edge_list = list(NMG.edges()) @@ -335,15 +355,20 @@ def __init__(self, system, cutoff, args): SG = NMG.subgraph(linker) all_cycles = nx.simple_cycles(nx.to_directed(SG)) - all_cycles = set([tuple(sorted(cy)) for cy in all_cycles if len(cy) > 4]) + all_cycles = set([tuple(sorted(cy)) + for cy in all_cycles if len(cy) > 4]) for cycle in all_cycles: - # rotate the ring normal vec onto the z-axis to determine coplanarity - fcoords = np.array([system['graph'].nodes[c]['fractional_position'] for c in cycle]) - element_symbols = [system['graph'].nodes[c]['element_symbol'] for c in cycle] + # rotate the ring normal vec onto the z-axis to determine + # coplanarity + fcoords = np.array( + [system['graph'].nodes[c]['fractional_position'] for c in cycle]) + element_symbols = [system['graph'].nodes[c] + ['element_symbol'] for c in cycle] anchor = fcoords[0] - fcoords = np.array([vec - PBC3DF_sym(anchor, vec)[1] for vec in fcoords]) + fcoords = np.array( + [vec - PBC3DF_sym(anchor, vec)[1] for vec in fcoords]) coords = np.dot(self.unit_cell.T, fcoords.T).T coords -= np.average(coords, axis=0) @@ -425,7 +450,8 @@ def type_atoms(self): SG.node[name]['hybridization'] = hyb if element_symbol not in metals: - gaff_ty, gaff_hyb = GAFF_type(name, inf, SG, self.pure_aromatic_atoms, self.aromatic_atoms) + gaff_ty, gaff_hyb = GAFF_type( + name, inf, SG, self.pure_aromatic_atoms, self.aromatic_atoms) imidazolate_gaff_types[ty] = gaff_ty # type non-imidazolate-ring atoms @@ -443,16 +469,21 @@ def type_atoms(self): # imidazolate ring atom adjacent that are not hydrogens if element_symbol != 'H': - if name not in imidazolate_ring_atoms and any(n in imidazolate_ring_atoms for n in nbors): - # special case for ZIF-8 -CH3 group which has an explicit type in ZIF-FF - if element_symbol == 'C' and sorted(nbor_symbols) == ['C', 'H', 'H', 'H'] and 'C1' in nbor_symbols: + if name not in imidazolate_ring_atoms and any( + n in imidazolate_ring_atoms for n in nbors): + # special case for ZIF-8 -CH3 group which has an + # explicit type in ZIF-FF + if element_symbol == 'C' and sorted(nbor_symbols) == [ + 'C', 'H', 'H', 'H'] and 'C1' in nbor_symbols: ty = 'C1' hyb = 'sp3' # other functionalizations are typed for GAFF else: - ty, hyb = GAFF_type(name, inf, SG, self.pure_aromatic_atoms, self.aromatic_atoms) + ty, hyb = GAFF_type( + name, inf, SG, self.pure_aromatic_atoms, self.aromatic_atoms) elif name not in imidazolate_ring_atoms and not any(n in imidazolate_ring_atoms for n in nbors): - ty, hyb = GAFF_type(name, inf, SG, self.pure_aromatic_atoms, self.aromatic_atoms) + ty, hyb = GAFF_type( + name, inf, SG, self.pure_aromatic_atoms, self.aromatic_atoms) types.append((ty, element_symbol, mass)) SG.node[name]['force_field_type'] = ty @@ -478,7 +509,8 @@ def type_atoms(self): ty = 'H3' hyb = 'sp1' else: - ty, hyb = GAFF_type(name, inf, SG, self.pure_aromatic_atoms, self.aromatic_atoms) + ty, hyb = GAFF_type( + name, inf, SG, self.pure_aromatic_atoms, self.aromatic_atoms) types.append((ty, element_symbol, mass)) SG.node[name]['force_field_type'] = ty @@ -486,7 +518,8 @@ def type_atoms(self): 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) @@ -512,13 +545,16 @@ def bond_parameters(self, bond): if all(t in self.ZIFFF_types for t in (i, j)): - k_ij, r_ij = parameter_loop([bond, bond[::-1]], ZIFFF_constants.ZIFFF_bonds) + k_ij, r_ij = parameter_loop( + [bond, bond[::-1]], ZIFFF_constants.ZIFFF_bonds) params = ('harmonic', k_ij, r_ij) else: - new_bond = tuple(sorted([self.imidazolate_gaff_types[ty] if ty in self.imidazolate_gaff_types else ty for ty in bond])) - k_ij, r_ij = parameter_loop([bond, bond[::-1], new_bond, new_bond[::-1]], gaff_bonds) + new_bond = tuple(sorted( + [self.imidazolate_gaff_types[ty] if ty in self.imidazolate_gaff_types else ty for ty in bond])) + k_ij, r_ij = parameter_loop( + [bond, bond[::-1], new_bond, new_bond[::-1]], gaff_bonds) params = ('harmonic', k_ij, r_ij) return params @@ -530,13 +566,16 @@ def angle_parameters(self, angle): if all(t in self.ZIFFF_types for t in (i, j, k)): - K, theta0 = parameter_loop([angle, angle[::-1]], ZIFFF_constants.ZIFFF_angles) + K, theta0 = parameter_loop( + [angle, angle[::-1]], ZIFFF_constants.ZIFFF_angles) params = ('harmonic', K, theta0) else: - new_angle = tuple([self.imidazolate_gaff_types[ty] if ty in self.imidazolate_gaff_types else ty for ty in angle]) - K, theta0 = parameter_loop([angle, angle[::-1], new_angle, new_angle[::-1]], gaff_angles) + new_angle = tuple([self.imidazolate_gaff_types[ty] + if ty in self.imidazolate_gaff_types else ty for ty in angle]) + K, theta0 = parameter_loop( + [angle, angle[::-1], new_angle, new_angle[::-1]], gaff_angles) if K is None: print(new_angle) params = ('harmonic', K, theta0) @@ -550,7 +589,8 @@ def dihedral_parameters(self, dihedral): if all(t in self.ZIFFF_types for t in (i, j, k, LL)): - params = dihedral_parameter_loop([dihedral, dihedral[::-1]], ZIFFF_constants.ZIFFF_dihedrals) + params = dihedral_parameter_loop( + [dihedral, dihedral[::-1]], ZIFFF_constants.ZIFFF_dihedrals) if params is not None: @@ -560,10 +600,19 @@ def dihedral_parameters(self, dihedral): else: X_dihedral = ('X', dihedral[1], dihedral[2], 'X') - new_X_dihedral = tuple([self.imidazolate_gaff_types[ty] if ty in self.imidazolate_gaff_types else ty for ty in X_dihedral]) - new_dihedral = tuple([self.imidazolate_gaff_types[ty] if ty in self.imidazolate_gaff_types else ty for ty in dihedral]) - - options = [X_dihedral, X_dihedral[::-1], new_X_dihedral, new_X_dihedral[::-1], dihedral, dihedral[::-1], new_dihedral, new_dihedral[::-1]] + new_X_dihedral = tuple( + [self.imidazolate_gaff_types[ty] if ty in self.imidazolate_gaff_types else ty for ty in X_dihedral]) + new_dihedral = tuple([self.imidazolate_gaff_types[ty] + if ty in self.imidazolate_gaff_types else ty for ty in dihedral]) + + options = [X_dihedral, + X_dihedral[::-1], + new_X_dihedral, + new_X_dihedral[::-1], + dihedral, + dihedral[::-1], + new_dihedral, + new_dihedral[::-1]] params = dihedral_parameter_loop(options, gaff_dihedrals) if params is not None: @@ -616,7 +665,8 @@ def improper_parameters(self, improper): if params is None: - i, j, k, LL = tuple([self.imidazolate_gaff_types[ty] if ty in self.imidazolate_gaff_types else ty for ty in improper]) + i, j, k, LL = tuple( + [self.imidazolate_gaff_types[ty] if ty in self.imidazolate_gaff_types else ty for ty in improper]) improper_combs = [(j, k, i, LL), (j, LL, i, k), (k, j, i, LL), @@ -667,7 +717,8 @@ def pair_parameters(self, charges=False): style = 'lj/cut/coul/long' sb = 'lj 0.0 0.0 0.5 coul 0.0 0.0 0.6874' else: - warnings.warn('ZIF-FF or any AMBER based force-field always uses charges, your results will be incorrect without them') + warnings.warn( + 'ZIF-FF or any AMBER based force-field always uses charges, your results will be incorrect without them') style = 'lj/cut' sb = 'lj 0.0 0.0 1.0' @@ -685,7 +736,11 @@ def pair_parameters(self, charges=False): comments[a] = [a, a] - self.pair_data = {'params': all_params, 'style': style, 'special_bonds': sb, 'comments': comments} + self.pair_data = { + 'params': all_params, + 'style': style, + 'special_bonds': sb, + 'comments': comments} def enumerate_bonds(self): @@ -721,7 +776,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): @@ -781,7 +843,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): @@ -800,7 +869,8 @@ def enumerate_dihedrals(self): nbors_k = [n for n in SG.neighbors(k) if n != j] il_pairs = list(itertools.product(nbors_j, nbors_k)) - dihedral_list = [(SG.node[p[0]]['force_field_type'], fft_j, fft_k, SG.node[p[1]]['force_field_type']) for p in il_pairs] + dihedral_list = [(SG.node[p[0]]['force_field_type'], fft_j, + fft_k, SG.node[p[1]]['force_field_type']) for p in il_pairs] for dihedral in dihedral_list: @@ -829,8 +899,14 @@ def enumerate_dihedrals(self): dihedral_comments[ID] = list(d) count += len(dihedrals[d]) - self.dihedral_data = {'all_dihedrals': all_dihedrals, 'params': indexed_dihedral_params, - 'style': 'harmonic', 'count': (count, len(all_dihedrals)), 'comments': dihedral_comments} + self.dihedral_data = { + 'all_dihedrals': all_dihedrals, + 'params': indexed_dihedral_params, + 'style': 'harmonic', + 'count': ( + count, + len(all_dihedrals)), + 'comments': dihedral_comments} def enumerate_impropers(self): @@ -844,7 +920,8 @@ def enumerate_impropers(self): if len(nbors) == 3: - nbors = sorted([(m, SG.node[m]['force_field_type']) for m in nbors], key=lambda x: x[1]) + nbors = sorted([(m, SG.node[m]['force_field_type']) + for m in nbors], key=lambda x: x[1]) fft_nbors = [m[1] for m in nbors] nbors = [m[0] for m in nbors] @@ -875,8 +952,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': 'fourier', 'count': (count, len(all_impropers)), 'comments': improper_comments} + self.improper_data = { + 'all_impropers': all_impropers, + 'params': improper_params, + 'style': 'fourier', + 'count': ( + count, + len(all_impropers)), + 'comments': improper_comments} def compile_force_field(self, charges=False): diff --git a/mofa/simulation/cif2lammps/cif2system.py b/mofa/simulation/cif2lammps/cif2system.py index 3d97d409..bf4b83a2 100644 --- a/mofa/simulation/cif2lammps/cif2system.py +++ b/mofa/simulation/cif2lammps/cif2system.py @@ -15,12 +15,112 @@ metals = atomic_data.metals mass_key = atomic_data.mass_key -PT = ['H', 'He', 'Li', 'Be', 'B', 'C', 'N', 'O', 'F', 'Ne', 'Na', 'Mg', 'Al', 'Si', 'P', 'S', 'Cl', 'Ar', - 'K', 'Ca', 'Sc', 'Ti', 'V', 'Cr', 'Mn', 'Fe', 'Co', 'Ni', 'Cu', 'Zn', 'Ga', 'Ge', 'As', 'Se', 'Br', 'Kr', - 'Rb', 'Sr', 'Y', 'Zr', 'Nb', 'Mo', 'Tc', 'Ru', 'Rh', 'Pd', 'Ag', 'Cd', 'In', 'Sn', 'Sb', 'Te', 'I', 'Xe', - 'Cs', 'Ba', 'Hf', 'Ta', 'W', 'Re', 'Os', 'Ir', 'Pt', 'Au', 'Hg', 'Tl', 'Pb', 'Bi', 'Po', 'At', 'Rn', 'Fr', - 'Ra', 'La', 'Ce', 'Pr', 'Nd', 'Pm', 'Sm', 'Eu', 'Gd', 'Tb', 'Dy', 'Ho', 'Er', 'Tm', 'Yb', 'Lu', 'Ac', 'Th', - 'Pa', 'U', 'Np', 'Pu', 'Am', 'Cm', 'Bk', 'Cf', 'Es', 'Fm', 'Md', 'No', 'Lr', 'FG', 'X'] +PT = [ + 'H', + 'He', + 'Li', + 'Be', + 'B', + 'C', + 'N', + 'O', + 'F', + 'Ne', + 'Na', + 'Mg', + 'Al', + 'Si', + 'P', + 'S', + 'Cl', + 'Ar', + 'K', + 'Ca', + 'Sc', + 'Ti', + 'V', + 'Cr', + 'Mn', + 'Fe', + 'Co', + 'Ni', + 'Cu', + 'Zn', + 'Ga', + 'Ge', + 'As', + 'Se', + 'Br', + 'Kr', + 'Rb', + 'Sr', + 'Y', + 'Zr', + 'Nb', + 'Mo', + 'Tc', + 'Ru', + 'Rh', + 'Pd', + 'Ag', + 'Cd', + 'In', + 'Sn', + 'Sb', + 'Te', + 'I', + 'Xe', + 'Cs', + 'Ba', + 'Hf', + 'Ta', + 'W', + 'Re', + 'Os', + 'Ir', + 'Pt', + 'Au', + 'Hg', + 'Tl', + 'Pb', + 'Bi', + 'Po', + 'At', + 'Rn', + 'Fr', + 'Ra', + 'La', + 'Ce', + 'Pr', + 'Nd', + 'Pm', + 'Sm', + 'Eu', + 'Gd', + 'Tb', + 'Dy', + 'Ho', + 'Er', + 'Tm', + 'Yb', + 'Lu', + 'Ac', + 'Th', + 'Pa', + 'U', + 'Np', + 'Pu', + 'Am', + 'Cm', + 'Bk', + 'Cf', + 'Es', + 'Fm', + 'Md', + 'No', + 'Lr', + 'FG', + 'X'] def GCD(a, b): @@ -58,7 +158,8 @@ def iscoord(line): """ identifies coordinates in CIFs """ - if nn(line[0]) in PT and line[1] in PT and False not in map(isfloat, line[2:5]): + if nn(line[0]) in PT and line[1] in PT and False not in map( + isfloat, line[2:5]): return True else: return False @@ -68,7 +169,8 @@ def isbond(line): """ identifies bonding in cifs """ - if nn(line[0]) in PT and nn(line[1]) in PT and isfloat(line[2]) and line[-1] in ('S', 'D', 'T', 'A'): + if nn(line[0]) in PT and nn(line[1]) in PT and isfloat( + line[2]) and line[-1] in ('S', 'D', 'T', 'A'): return True else: return False @@ -79,7 +181,8 @@ def PBC3DF_sym(vec1, vec2): applies periodic boundary to distance between vec1 and vec2 (fractional coordinates) """ dist = vec1 - vec2 - sym_dist = [(-1.0, dim - 1.0) if dim > 0.5 else (1.0, dim + 1.0) if dim < -0.5 else (0, dim) for dim in dist] + sym_dist = [(-1.0, dim - 1.0) if dim > 0.5 else (1.0, dim + 1.0) + if dim < -0.5 else (0, dim) for dim in dist] sym = np.array([s[0] for s in sym_dist]) ndist = np.array([s[1] for s in sym_dist]) @@ -157,7 +260,10 @@ def cif_read(filename, charges=False, add_Zr_bonds=False): net_charge = np.round(np.sum(charge_list), 3) if net_charge > 0.1 and charges: - warnings.warn('A potentially significant net charge of ' + str(net_charge) + ' is being removed') + warnings.warn( + 'A potentially significant net charge of ' + + str(net_charge) + + ' is being removed') if len(charge_list) > 0: remove_net = choice(range(len(charge_list))) @@ -181,35 +287,56 @@ def cif_read(filename, charges=False, add_Zr_bonds=False): if dist < 4.5: count += 1 - bonds.append([names[i], names[j], '.', 'S', np.round(dist, 3)]) + bonds.append( + [names[i], names[j], '.', 'S', np.round(dist, 3)]) print(count, 'Zr-Zr bonds added...') - return elems, names, ccoords, fcoords, charge_list, bonds, (a, b, c, alpha, beta, gamma), unit_cell + return elems, names, ccoords, fcoords, charge_list, bonds, ( + a, b, c, alpha, beta, gamma), unit_cell -def initialize_system(filename, charges=False, small_molecule_cutoff=5, read_pymatgen=False): +def initialize_system( + filename, + charges=False, + small_molecule_cutoff=5, + read_pymatgen=False): if not read_pymatgen: - elems, names, ccoords, fcoords, charge_list, bonds, uc_params, unit_cell = cif_read(filename, charges=charges) + elems, names, ccoords, fcoords, charge_list, bonds, uc_params, unit_cell = cif_read( + filename, charges=charges) else: from .pymatgen_cif2system import cif_read_pymatgen - elems, names, ccoords, fcoords, charge_list, bonds, uc_params, unit_cell = cif_read_pymatgen(filename, charges=charges) + elems, names, ccoords, fcoords, charge_list, bonds, uc_params, unit_cell = cif_read_pymatgen( + filename, charges=charges) A, B, C, alpha, beta, gamma = uc_params G = nx.Graph() index = 0 index_key = {} - for e, n, cc, fc, charge in zip(elems, names, ccoords, fcoords, charge_list): + for e, n, cc, fc, charge in zip( + elems, names, ccoords, fcoords, charge_list): index += 1 - G.add_node(index, element_symbol=e, mol_flag='1', index=index, force_field_type='', cartesian_position=cc, - fractional_position=fc, charge=charge, replication=np.array([0.0, 0.0, 0.0]), duplicated_version_of=None, cif_label=n) + G.add_node(index, + element_symbol=e, + mol_flag='1', + index=index, + force_field_type='', + cartesian_position=cc, + fractional_position=fc, + charge=charge, + replication=np.array([0.0, + 0.0, + 0.0]), + duplicated_version_of=None, + cif_label=e + ("%d" % index)) index_key[n] = index for b in bonds: - dist, sym = PBC3DF_sym(G.nodes[index_key[b[0]]]['fractional_position'], G.nodes[index_key[b[1]]]['fractional_position']) + dist, sym = PBC3DF_sym(G.nodes[index_key[b[0]]]['fractional_position'], + G.nodes[index_key[b[1]]]['fractional_position']) if np.any(sym): sym_code = '1_' + ''.join(map(str, map(int, sym + 5))) @@ -218,7 +345,11 @@ def initialize_system(filename, charges=False, small_molecule_cutoff=5, read_pym dist = np.linalg.norm(np.dot(unit_cell, dist)) - G.add_edge(index_key[b[0]], index_key[b[1]], sym_code=sym_code, bond_type=b[3], length=float(b[-1])) + G.add_edge(index_key[b[0]], + index_key[b[1]], + sym_code=sym_code, + bond_type=b[3], + length=float(b[-1])) print_flag = False for e in G.edges(data=True): @@ -233,18 +364,22 @@ def initialize_system(filename, charges=False, small_molecule_cutoff=5, read_pym bond_type = data['bond_type'] # carboxylate oxygens should have aromatic bonds with C - if len(nbors0_symbols) == 2 and es0 == 'O' and es1 == 'C' and any(i in metals for i in nbors0_symbols) and bond_type != 'A': + if len(nbors0_symbols) == 2 and es0 == 'O' and es1 == 'C' and any( + i in metals for i in nbors0_symbols) and bond_type != 'A': print_flag = True data['bond_type'] = 'A' - if len(nbors1_symbols) == 2 and es1 == 'O' and es0 == 'C' and any(i in metals for i in nbors1_symbols) and bond_type != 'A': + if len(nbors1_symbols) == 2 and es1 == 'O' and es0 == 'C' and any( + i in metals for i in nbors1_symbols) and bond_type != 'A': print_flag = True data['bond_type'] = 'A' # nitro nitrogens should have aromatic bonds with O - if es0 == 'N' and es1 == 'O' and sorted(nbors0_symbols) == ['C', 'O', 'O']: + if es0 == 'N' and es1 == 'O' and sorted( + nbors0_symbols) == ['C', 'O', 'O']: print_flag = True data['bond_type'] = 'A' - if es1 == 'N' and es0 == 'O' and sorted(nbors1_symbols) == ['C', 'O', 'O']: + if es1 == 'N' and es0 == 'O' and sorted( + nbors1_symbols) == ['C', 'O', 'O']: print_flag = True data['bond_type'] = 'A' @@ -264,11 +399,15 @@ def initialize_system(filename, charges=False, small_molecule_cutoff=5, read_pym for es in comp_dict: comp_dict[es] = int(comp_dict[es] / float(counts)) - comp = tuple(sorted([(key, val) for key, val in comp_dict.items()], key=lambda x: x[0])) + comp = tuple( + sorted([(key, val) for key, val in comp_dict.items()], key=lambda x: x[0])) formula = ''.join([str(x) for es in comp for x in es]) components.append((len(elems), formula, S)) - print('there are', len(components), 'components in the system with ( # atoms, formula unit):') + print( + 'there are', + len(components), + 'components in the system with ( # atoms, formula unit):') SM = nx.Graph() framework = nx.Graph() @@ -282,10 +421,12 @@ def initialize_system(filename, charges=False, small_molecule_cutoff=5, read_pym if len(S.nodes()) < small_molecule_cutoff: - node_elems = [(n, data['element_symbol']) for n, data in S.nodes(data=True)] + node_elems = [(n, data['element_symbol']) + for n, data in S.nodes(data=True)] sort_elems = sorted(node_elems, key=lambda x: x[1], reverse=True) sort_index = sorted(node_elems, key=lambda x: x[0], reverse=False) - sort_key = dict((i[0], j[0]) for i, j in zip(sort_index, sort_elems)) + sort_key = dict((i[0], j[0]) + for i, j in zip(sort_index, sort_elems)) add_graph = nx.Graph() for node, elem in sort_elems: @@ -314,7 +455,17 @@ def initialize_system(filename, charges=False, small_molecule_cutoff=5, read_pym data['index'] = index SM = nx.relabel_nodes(SM, sm_remap) - return {'box': (A, B, C, alpha, beta, gamma), 'graph': framework, 'SM_graph': SM, 'max_ind': index} + return { + 'box': ( + A, + B, + C, + alpha, + beta, + gamma), + 'graph': framework, + 'SM_graph': SM, + 'max_ind': index} def duplicate_system(system, replications, small_molecule_cutoff=10): @@ -328,7 +479,16 @@ def duplicate_system(system, replications, small_molecule_cutoff=10): box = system['box'] replications = list(map(int, replications.split('x'))) - replicated_box = (box[0] * replications[0], box[1] * replications[1], box[2] * replications[2], box[3], box[4], box[5]) + replicated_box = ( + box[0] * + replications[0], + box[1] * + replications[1], + box[2] * + replications[2], + box[3], + box[4], + box[5]) pi = np.pi a, b, c, alpha, beta, gamma = replicated_box @@ -343,16 +503,21 @@ def duplicate_system(system, replications, small_molecule_cutoff=10): cz = (c ** 2.0 - cx ** 2.0 - cy ** 2.0) ** 0.5 unit_cell = np.asarray([[ax, ay, az], [bx, by, bz], [cx, cy, cz]]).T - basis_vecs = [np.array([1, 0, 0]), np.array([0, 1, 0]), np.array([0, 0, 1])] - dim0 = [[np.array([0, 0, 0])]] + [[np.array([0, 0, 0])] + [basis_vecs[0] for i in range(r + 1)] for r in range(replications[0] - 1)] - dim1 = [[np.array([0, 0, 0])]] + [[np.array([0, 0, 0])] + [basis_vecs[1] for i in range(r + 1)] for r in range(replications[1] - 1)] - dim2 = [[np.array([0, 0, 0])]] + [[np.array([0, 0, 0])] + [basis_vecs[2] for i in range(r + 1)] for r in range(replications[2] - 1)] + basis_vecs = [np.array([1, 0, 0]), np.array( + [0, 1, 0]), np.array([0, 0, 1])] + dim0 = [[np.array([0, 0, 0])]] + [[np.array([0, 0, 0])] + [basis_vecs[0] + for i in range(r + 1)] for r in range(replications[0] - 1)] + dim1 = [[np.array([0, 0, 0])]] + [[np.array([0, 0, 0])] + [basis_vecs[1] + for i in range(r + 1)] for r in range(replications[1] - 1)] + dim2 = [[np.array([0, 0, 0])]] + [[np.array([0, 0, 0])] + [basis_vecs[2] + for i in range(r + 1)] for r in range(replications[2] - 1)] dim0 = [np.sum([v for v in comb], axis=0) for comb in dim0] dim1 = [np.sum([v for v in comb], axis=0) for comb in dim1] dim2 = [np.sum([v for v in comb], axis=0) for comb in dim2] - trans_vecs = [np.sum(comb, axis=0) for comb in itertools.product(dim0, dim1, dim2)] + trans_vecs = [np.sum(comb, axis=0) + for comb in itertools.product(dim0, dim1, dim2)] trans_vecs = [v for v in trans_vecs if np.any(v)] print('The transformation vectors for the replication are:') @@ -360,8 +525,10 @@ def duplicate_system(system, replications, small_molecule_cutoff=10): print(vec) print('...') - if len(trans_vecs) != replications[0] * replications[1] * replications[2] - 1: - raise ValueError('The number of transformation vectors in the replication is wrong somehow') + if len(trans_vecs) != replications[0] * \ + replications[1] * replications[2] - 1: + raise ValueError( + 'The number of transformation vectors in the replication is wrong somehow') NG = G.copy() edge_remove_list = [] @@ -387,16 +554,27 @@ def duplicate_system(system, replications, small_molecule_cutoff=10): fvec = node_data['fractional_position'] translated_fvec = fvec + trans_vec fvec = np.array([c / d for c, d in zip(fvec, replications)]) - translated_fvec = np.array([c / d for c, d in zip(translated_fvec, replications)]) + translated_fvec = np.array( + [c / d for c, d in zip(translated_fvec, replications)]) NG.nodes[node]['fractional_position'] = fvec equivalency[original_atom].append(new_index) - NG.add_node(new_index, element_symbol=element_symbol, mol_flag=1, index=new_index, force_field_type='', cartesian_position=np.array( - [0.0, 0.0, 0.0]), fractional_position=translated_fvec, charge=charge, duplicated_version_of=original_atom) + NG.add_node(new_index, + element_symbol=element_symbol, + mol_flag=1, + index=new_index, + force_field_type='', + cartesian_position=np.array([0.0, + 0.0, + 0.0]), + fractional_position=translated_fvec, + charge=charge, + duplicated_version_of=original_atom) for node, data in NG.nodes(data=True): - data['cartesian_position'] = np.dot(unit_cell, data['fractional_position']) + data['cartesian_position'] = np.dot( + unit_cell, data['fractional_position']) for n0, n1, edge_data in G.edges(data=True): @@ -428,29 +606,47 @@ def duplicate_system(system, replications, small_molecule_cutoff=10): if abs(dist_e0e1 - length) < 0.075: if np.any(sym_e0e1): - sym_code = '1_' + ''.join(map(str, map(int, sym_e0e1 + 5))) + sym_code = '1_' + \ + ''.join(map(str, map(int, sym_e0e1 + 5))) else: sym_code = '.' - NG.add_edge(eq0, eq1, sym_code=sym_code, bond_type=bond_type, length=dist_e0e1) + NG.add_edge( + eq0, + eq1, + sym_code=sym_code, + bond_type=bond_type, + length=dist_e0e1) if abs(dist_n0e1 - length) < 0.075: if np.any(sym_n0e1): - sym_code = '1_' + ''.join(map(str, map(int, sym_n0e1 + 5))) + sym_code = '1_' + \ + ''.join(map(str, map(int, sym_n0e1 + 5))) else: sym_code = '.' - NG.add_edge(n0, eq1, sym_code=sym_code, bond_type=bond_type, length=dist_n0e1) + NG.add_edge( + n0, + eq1, + sym_code=sym_code, + bond_type=bond_type, + length=dist_n0e1) if abs(dist_e0n1 - length) < 0.075: if np.any(sym_e0n1): - sym_code = '1_' + ''.join(map(str, map(int, sym_e0n1 + 5))) + sym_code = '1_' + \ + ''.join(map(str, map(int, sym_e0n1 + 5))) else: sym_code = '.' - NG.add_edge(eq0, n1, sym_code=sym_code, bond_type=bond_type, length=dist_n0e1) + NG.add_edge( + eq0, + n1, + sym_code=sym_code, + bond_type=bond_type, + length=dist_n0e1) if abs(dist_n0n1 - length) > 0.075: if (n0, n1) not in edge_remove_list: @@ -472,11 +668,15 @@ def duplicate_system(system, replications, small_molecule_cutoff=10): for es in comp_dict: comp_dict[es] = int(comp_dict[es] / float(counts)) - comp = tuple(sorted([(key, val) for key, val in comp_dict.items()], key=lambda x: x[0])) + comp = tuple( + sorted([(key, val) for key, val in comp_dict.items()], key=lambda x: x[0])) formula = ''.join([str(x) for es in comp for x in es]) components.append((len(elems), formula, S)) - print('there are', len(components), 'components in the system with ( # atoms, formula unit):') + print( + 'there are', + len(components), + 'components in the system with ( # atoms, formula unit):') SM = nx.Graph() framework = nx.Graph() for component in components: @@ -504,7 +704,11 @@ def duplicate_system(system, replications, small_molecule_cutoff=10): MI = max([data['index'] for n, data in NG.nodes(data=True)]) - return {'box': replicated_box, 'graph': framework, 'SM_graph': SM, 'max_ind': MI} + return { + 'box': replicated_box, + 'graph': framework, + 'SM_graph': SM, + 'max_ind': MI} def replication_determination(system, replication, cutoff): @@ -527,13 +731,19 @@ def replication_determination(system, replication, cutoff): bvec = np.array([bx, by, bz]) cvec = np.array([cx, cy, cz]) - thetac = np.arccos(np.dot(np.cross(avec, bvec), cvec) / (np.linalg.norm(np.cross(avec, bvec)) * np.linalg.norm(cvec))) + thetac = np.arccos(np.dot(np.cross(avec, bvec), cvec) / + (np.linalg.norm(np.cross(avec, bvec)) * + np.linalg.norm(cvec))) dist2 = np.absolute(np.linalg.norm(cvec) * np.cos(thetac)) - thetab = np.arccos(np.dot(np.cross(avec, cvec), bvec) / (np.linalg.norm(np.cross(avec, cvec)) * np.linalg.norm(bvec))) + thetab = np.arccos(np.dot(np.cross(avec, cvec), bvec) / + (np.linalg.norm(np.cross(avec, cvec)) * + np.linalg.norm(bvec))) dist1 = np.absolute(np.linalg.norm(bvec) * np.cos(thetab)) - thetaa = np.arccos(np.dot(np.cross(cvec, bvec), avec) / (np.linalg.norm(np.cross(cvec, bvec)) * np.linalg.norm(avec))) + thetaa = np.arccos(np.dot(np.cross(cvec, bvec), avec) / + (np.linalg.norm(np.cross(cvec, bvec)) * + np.linalg.norm(avec))) dist0 = np.absolute(np.linalg.norm(avec) * np.cos(thetaa)) if 'min_atoms' in replication: @@ -558,8 +768,13 @@ def replication_determination(system, replication, cutoff): rvals = range(duplications + 1)[1:] shapes = itertools.product(rvals, rvals, rvals) - shapes = [s for s in shapes if functools.reduce((lambda x, y: x * y), s) == duplications] - useable_shapes = [s for s in shapes if min(dist0 * s[0], dist1 * s[1], dist2 * s[2]) >= 2 * cutoff] + shapes = [s for s in shapes if functools.reduce( + (lambda x, y: x * y), s) == duplications] + useable_shapes = [ + s for s in shapes if min( + dist0 * s[0], + dist1 * s[1], + dist2 * s[2]) >= 2 * cutoff] useable_shapes = [s for s in useable_shapes if max([((a * s[0]) / (b * s[1])), ((a * s[0]) / (c * s[2])), ((b * s[1]) / (c * s[2])), @@ -572,14 +787,23 @@ def replication_determination(system, replication, cutoff): print('final duplications:', duplications) print('final number of atoms:', int(duplications * Natoms)) - shape_deviations = [(i, np.std([useable_shapes[i][0] * a, useable_shapes[i][1] * b, useable_shapes[i][2] * c])) for i in range(len(useable_shapes))] + shape_deviations = [(i, np.std([useable_shapes[i][0] * + a, useable_shapes[i][1] * + b, useable_shapes[i][2] * + c])) for i in range(len(useable_shapes))] shape_deviations.sort(key=lambda x: x[1]) selected_shape = useable_shapes[shape_deviations[0][0]] replication = 'x'.join(map(str, selected_shape)) - print('replicating to a', replication, 'cell (' + str(duplications) + ' duplications)...') + print( + 'replicating to a', + replication, + 'cell (' + + str(duplications) + + ' duplications)...') system = duplicate_system(system, replication) - print('the minimum boundary-boundary distance is', min([d * s for d, s in zip(selected_shape, (dist0, dist1, dist2))])) + print('the minimum boundary-boundary distance is', + min([d * s for d, s in zip(selected_shape, (dist0, dist1, dist2))])) replication = 'ma' + str(min_atoms) a, b, c, alpha, beta, gamma = system['box'] @@ -590,9 +814,12 @@ def replication_determination(system, replication, cutoff): yz = np.round((b * c * np.cos(math.radians(alpha)) - xy * xz) / ly, 8) lz = np.round(np.sqrt(c**2 - xz**2 - yz**2), 8) - print('lx =', np.round(lx, 3), '(dim 0 separation = ' + str(np.round(selected_shape[0] * dist0, 3)) + ')') - print('ly =', np.round(ly, 3), '(dim 1 separation = ' + str(np.round(selected_shape[1] * dist1, 3)) + ')') - print('lz =', np.round(lz, 3), '(dim 2 separation = ' + str(np.round(selected_shape[2] * dist2, 3)) + ')') + print('lx =', np.round(lx, 3), '(dim 0 separation = ' + + str(np.round(selected_shape[0] * dist0, 3)) + ')') + print('ly =', np.round(ly, 3), '(dim 1 separation = ' + + str(np.round(selected_shape[1] * dist1, 3)) + ')') + print('lz =', np.round(lz, 3), '(dim 2 separation = ' + + str(np.round(selected_shape[2] * dist2, 3)) + ')') print('alpha =', np.round(alpha, 3)) print('beta =', np.round(beta, 3)) print('gamma =', np.round(gamma, 3)) @@ -613,21 +840,35 @@ def replication_determination(system, replication, cutoff): rvals = range(duplications + 1)[1:] shapes = itertools.product(rvals, rvals, rvals) - shapes = [s for s in shapes if functools.reduce((lambda x, y: x * y), s) == duplications] - useable_shapes = [s for s in shapes if min(dist0 * s[0], dist1 * s[1], dist2 * s[2]) >= 2 * cutoff] + shapes = [s for s in shapes if functools.reduce( + (lambda x, y: x * y), s) == duplications] + useable_shapes = [ + s for s in shapes if min( + dist0 * s[0], + dist1 * s[1], + dist2 * s[2]) >= 2 * cutoff] duplications += 1 duplications -= 1 print('final duplications:', duplications) - shape_deviations = [(i, np.std([useable_shapes[i][0] * a, useable_shapes[i][1] * b, useable_shapes[i][2] * c])) for i in range(len(useable_shapes))] + shape_deviations = [(i, np.std([useable_shapes[i][0] * + a, useable_shapes[i][1] * + b, useable_shapes[i][2] * + c])) for i in range(len(useable_shapes))] shape_deviations.sort(key=lambda x: x[1]) selected_shape = useable_shapes[shape_deviations[0][0]] replication = 'x'.join(map(str, selected_shape)) - print('replicating to a', replication, 'cell (' + str(duplications) + ' duplications)...') + print( + 'replicating to a', + replication, + 'cell (' + + str(duplications) + + ' duplications)...') system = duplicate_system(system, replication) - print('the minimum boundary-boundary distance is', min([d * s for d, s in zip(selected_shape, (dist0, dist1, dist2))])) + print('the minimum boundary-boundary distance is', + min([d * s for d, s in zip(selected_shape, (dist0, dist1, dist2))])) a, b, c, alpha, beta, gamma = system['box'] lx = np.round(a, 8) @@ -637,9 +878,12 @@ def replication_determination(system, replication, cutoff): yz = np.round((b * c * np.cos(math.radians(alpha)) - xy * xz) / ly, 8) lz = np.round(np.sqrt(c**2 - xz**2 - yz**2), 8) - print('lx =', np.round(lx, 3), '(dim 0 separation = ' + str(np.round(selected_shape[0] * dist0, 3)) + ')') - print('ly =', np.round(ly, 3), '(dim 1 separation = ' + str(np.round(selected_shape[1] * dist1, 3)) + ')') - print('lz =', np.round(lz, 3), '(dim 2 separation = ' + str(np.round(selected_shape[2] * dist2, 3)) + ')') + print('lx =', np.round(lx, 3), '(dim 0 separation = ' + + str(np.round(selected_shape[0] * dist0, 3)) + ')') + print('ly =', np.round(ly, 3), '(dim 1 separation = ' + + str(np.round(selected_shape[1] * dist1, 3)) + ')') + print('lz =', np.round(lz, 3), '(dim 2 separation = ' + + str(np.round(selected_shape[2] * dist2, 3)) + ')') print('alpha =', np.round(alpha, 3)) print('beta =', np.round(beta, 3)) print('gamma =', np.round(gamma, 3)) @@ -671,7 +915,8 @@ def write_cif_from_system(system, filename): with open(filename, 'w') as out: out.write('data_' + filename[0:-4] + '\n') - out.write('_audit_creation_date ' + datetime.datetime.today().strftime('%Y-%m-%d') + '\n') + out.write('_audit_creation_date ' + + datetime.datetime.today().strftime('%Y-%m-%d') + '\n') out.write("_audit_creation_method 'cif2lammps'" + '\n') out.write("_symmetry_space_group_name_H-M 'P1'" + '\n') out.write('_symmetry_Int_Tables_number 1' + '\n') @@ -700,7 +945,9 @@ def write_cif_from_system(system, filename): ind = es + str(data['index']) index_dict[n] = ind chg = data['charge'] - out.write('{:7} {:>4} {:>15.6f} {:>15.6f} {:>15.6f} {:>15.6f}'.format(ind, es, vec[0], vec[1], vec[2], chg)) + out.write( + '{:7} {:>4} {:>15.6f} {:>15.6f} {:>15.6f} {:>15.6f}'.format( + ind, es, vec[0], vec[1], vec[2], chg)) out.write('\n') out.write('loop_' + '\n') @@ -717,5 +964,7 @@ def write_cif_from_system(system, filename): dist = np.round(data['length'], 3) bond_type = data['bond_type'] - out.write('{:7} {:>7} {:>7} {:>3}'.format(ind0, ind1, dist, bond_type)) + out.write( + '{:7} {:>7} {:>7} {:>3}'.format( + ind0, ind1, dist, bond_type)) out.write('\n') diff --git a/mofa/simulation/cif2lammps/gaff.py b/mofa/simulation/cif2lammps/gaff.py index 46eb08da..d4f0258e 100644 --- a/mofa/simulation/cif2lammps/gaff.py +++ b/mofa/simulation/cif2lammps/gaff.py @@ -1,6 +1,7 @@ from textwrap import dedent -# AMBER General Force Field for organic molecules (Version 1.4, March 2010) add. info. at the end +# AMBER General Force Field for organic molecules (Version 1.4, March +# 2010) add. info. at the end gaff_atom_types = dedent(""" c 12.01 0.616 Sp2 C carbonyl group c1 12.01 0.360 Sp C @@ -5701,7 +5702,8 @@ na-n2-ca-n2 1.1 180. 2. dac, 10/94 """) -# hw ow 0000. 0000. 4. flag for fast water +# hw ow 0000. 0000. 4. flag for fast +# water # MOD4 RE gaff_LJ_params = dedent(""" diff --git a/mofa/simulation/cif2lammps/gaff2.py b/mofa/simulation/cif2lammps/gaff2.py index c4ef7baf..161c4ab4 100644 --- a/mofa/simulation/cif2lammps/gaff2.py +++ b/mofa/simulation/cif2lammps/gaff2.py @@ -7526,7 +7526,8 @@ na-n2-ca-n2 1.1 180. 2. dac, 10/94 """) -# hw ow 0000. 0000. 4. flag for fast water +# hw ow 0000. 0000. 4. flag for fast +# water # MOD4 RE gaff_LJ_params = dedent(""" diff --git a/mofa/simulation/cif2lammps/main_conversion.py b/mofa/simulation/cif2lammps/main_conversion.py index 5c965e7b..7561cfaa 100644 --- a/mofa/simulation/cif2lammps/main_conversion.py +++ b/mofa/simulation/cif2lammps/main_conversion.py @@ -17,18 +17,44 @@ # add more force field classes here as they are made -def single_conversion(cif, force_field=UFF4MOF, ff_string='UFF4MOF', small_molecule_force_field=None, - outdir='unopt_lammps_data', charges=False, parallel=False, replication='1x1x1', read_cifs_pymatgen=False, add_molecule=None, - small_molecule_file=None): +def single_conversion( + cif, + force_field=UFF4MOF, + ff_string='UFF4MOF', + small_molecule_force_field=None, + outdir='unopt_lammps_data', + charges=False, + parallel=False, + replication='1x1x1', + read_cifs_pymatgen=False, + add_molecule=None, + small_molecule_file=None): print('converting ', cif, '...') - lammps_inputs([cif, force_field, ff_string, small_molecule_force_field, outdir, charges, - replication, read_cifs_pymatgen, add_molecule, small_molecule_file]) - - -def serial_conversion(directory, force_field=UFF4MOF, ff_string='UFF4MOF', small_molecule_force_field=None, - outdir='unopt_lammps_data', charges=False, parallel=False, replication='1x1x1', read_cifs_pymatgen=False, add_molecule=None, - small_molecule_file=None): + lammps_inputs([cif, + force_field, + ff_string, + small_molecule_force_field, + outdir, + charges, + replication, + read_cifs_pymatgen, + add_molecule, + small_molecule_file]) + + +def serial_conversion( + directory, + force_field=UFF4MOF, + ff_string='UFF4MOF', + small_molecule_force_field=None, + outdir='unopt_lammps_data', + charges=False, + parallel=False, + replication='1x1x1', + read_cifs_pymatgen=False, + add_molecule=None, + small_molecule_file=None): try: os.mkdir(outdir) @@ -40,15 +66,32 @@ def serial_conversion(directory, force_field=UFF4MOF, ff_string='UFF4MOF', small cifs = sorted(glob.glob(directory + os.sep + '*.cif')) for cif in cifs: print('converting ', cif, '...') - lammps_inputs([cif, force_field, ff_string, small_molecule_force_field, outdir, charges, - replication, read_cifs_pymatgen, add_molecule, small_molecule_file]) + lammps_inputs([cif, + force_field, + ff_string, + small_molecule_force_field, + outdir, + charges, + replication, + read_cifs_pymatgen, + add_molecule, + small_molecule_file]) print('--- cifs in', directory, 'converted and placed in', outdir, '---') -def parallel_conversion(directory, force_field=UFF4MOF, ff_string='UFF4MOF', small_molecule_force_field=None, - outdir='unopt_lammps_data', charges=False, parallel=True, replication='1x1x1', read_cifs_pymatgen=False, add_molecule=None, - small_molecule_file=None): +def parallel_conversion( + directory, + force_field=UFF4MOF, + ff_string='UFF4MOF', + small_molecule_force_field=None, + outdir='unopt_lammps_data', + charges=False, + parallel=True, + replication='1x1x1', + read_cifs_pymatgen=False, + add_molecule=None, + small_molecule_file=None): try: os.mkdir(outdir) @@ -58,8 +101,16 @@ def parallel_conversion(directory, force_field=UFF4MOF, ff_string='UFF4MOF', sma print('conversion running on ' + str(multiprocessing.cpu_count()) + ' cores') cifs = sorted(glob.glob(directory + os.sep + '*.cif')) - args = [[cif, force_field, ff_string, small_molecule_force_field, outdir, charges, - replication, read_cifs_pymatgen, add_molecule, small_molecule_file] for cif in cifs] + args = [[cif, + force_field, + ff_string, + small_molecule_force_field, + outdir, + charges, + replication, + read_cifs_pymatgen, + add_molecule, + small_molecule_file] for cif in cifs] pool = Pool(multiprocessing.cpu_count()) pool.map_async(lammps_inputs, args) # results_par = pool.map_async(lammps_inputs, args) @@ -69,8 +120,15 @@ def parallel_conversion(directory, force_field=UFF4MOF, ff_string='UFF4MOF', sma print('--- cifs in', directory, 'converted and placed in', outdir, '---') -def parallel_GULP_conversion(directory, force_field=UFF4MOF, outdir='GULP_inputs', charges=False, - parallel=True, replication='1x1x1', GULP=True, noautobond=True): +def parallel_GULP_conversion( + directory, + force_field=UFF4MOF, + outdir='GULP_inputs', + charges=False, + parallel=True, + replication='1x1x1', + GULP=True, + noautobond=True): try: os.mkdir(outdir) @@ -80,7 +138,8 @@ def parallel_GULP_conversion(directory, force_field=UFF4MOF, outdir='GULP_inputs print('conversion running on ' + str(multiprocessing.cpu_count()) + ' cores') cifs = sorted(glob.glob(directory + os.sep + '*.cif')) - args = [[cif, force_field, outdir, charges, replication, noautobond] for cif in cifs] + args = [[cif, force_field, outdir, charges, replication, noautobond] + for cif in cifs] pool = Pool(multiprocessing.cpu_count()) pool.map_async(GULP_inputs, args) # results_par = pool.map_async(GULP_inputs, args) @@ -92,23 +151,90 @@ def parallel_GULP_conversion(directory, force_field=UFF4MOF, outdir='GULP_inputs def run_conversion(): - parser = argparse.ArgumentParser(description='Optional arguments for running cif2lammps') - parser.add_argument('--cifs', action='store', dest='directory', type=str, required=True, help='the cifs to convert') - parser.add_argument('--force_field', action='store', dest='ff_string', type=str, required=False, default='UFF4MOF', help='the force field to use') - parser.add_argument('--small_molecule_force_field', action='store', dest='sm_ff_string', type=str, - required=False, default='TraPPE', help='the force field to use for small molecules') - parser.add_argument('--outdir', action='store', dest='outdir', type=str, required=False, - default='unopt_lammps_data', help='where to write the lammps inputs') - parser.add_argument('--charges', action='store_true', dest='charges', required=False, default=False, help='switch on charges') - parser.add_argument('--replication', action='store', dest='replication', type=str, required=False, default='1x1x1', help='replications to use') - parser.add_argument('--GULP', action='store_true', dest='GULP', required=False, default=False, help='write GULP inputs instead of LAMMPS') - parser.add_argument('--parallel', action='store_true', dest='parallel', required=False, default=False, help='switch on parallel conversion') - parser.add_argument('--read_cifs_pymatgen', action='store_true', dest='read_cifs_pymatgen', - required=False, default=False, help='use ASE to read CIF inputs') - parser.add_argument('--add_molecule', action='store', dest='add_molecule', required=False, - default=None, help='name of a molecule to write molecule file for') - parser.add_argument('--small_molecule_file', action='store', dest='sm_file', type=str, required=False, default=None, - help='a cif, xyz, or pdb of small molecules to be added, should be in cifs folder') + parser = argparse.ArgumentParser( + description='Optional arguments for running cif2lammps') + parser.add_argument( + '--cifs', + action='store', + dest='directory', + type=str, + required=True, + help='the cifs to convert') + parser.add_argument( + '--force_field', + action='store', + dest='ff_string', + type=str, + required=False, + default='UFF4MOF', + help='the force field to use') + parser.add_argument( + '--small_molecule_force_field', + action='store', + dest='sm_ff_string', + type=str, + required=False, + default='TraPPE', + help='the force field to use for small molecules') + parser.add_argument( + '--outdir', + action='store', + dest='outdir', + type=str, + required=False, + default='unopt_lammps_data', + help='where to write the lammps inputs') + parser.add_argument( + '--charges', + action='store_true', + dest='charges', + required=False, + default=False, + help='switch on charges') + parser.add_argument( + '--replication', + action='store', + dest='replication', + type=str, + required=False, + default='1x1x1', + help='replications to use') + parser.add_argument( + '--GULP', + action='store_true', + dest='GULP', + required=False, + default=False, + help='write GULP inputs instead of LAMMPS') + parser.add_argument( + '--parallel', + action='store_true', + dest='parallel', + required=False, + default=False, + help='switch on parallel conversion') + parser.add_argument( + '--read_cifs_pymatgen', + action='store_true', + dest='read_cifs_pymatgen', + required=False, + default=False, + help='use ASE to read CIF inputs') + parser.add_argument( + '--add_molecule', + action='store', + dest='add_molecule', + required=False, + default=None, + help='name of a molecule to write molecule file for') + parser.add_argument( + '--small_molecule_file', + action='store', + dest='sm_file', + type=str, + required=False, + default=None, + help='a cif, xyz, or pdb of small molecules to be added, should be in cifs folder') args = parser.parse_args() print(args) @@ -120,12 +246,24 @@ def run_conversion(): else: add_molecule = args.add_molecule - ff_dict = {'UFF4MOF': UFF4MOF, 'UFF': UFF, 'Dreiding': Dreiding, 'MZHB': MZHB, 'ZIFFF': ZIFFF} + ff_dict = { + 'UFF4MOF': UFF4MOF, + 'UFF': UFF, + 'Dreiding': Dreiding, + 'MZHB': MZHB, + 'ZIFFF': ZIFFF} force_field = ff_dict[args.ff_string] - optional_arguments = {'force_field': force_field, 'ff_string': args.ff_string, 'small_molecule_force_field': args.sm_ff_string, - 'outdir': args.outdir, 'charges': args.charges, 'replication': args.replication, 'read_cifs_pymatgen': args.read_cifs_pymatgen, - 'add_molecule': add_molecule, 'small_molecule_file': args.sm_file} + optional_arguments = { + 'force_field': force_field, + 'ff_string': args.ff_string, + 'small_molecule_force_field': args.sm_ff_string, + 'outdir': args.outdir, + 'charges': args.charges, + 'replication': args.replication, + 'read_cifs_pymatgen': args.read_cifs_pymatgen, + 'add_molecule': add_molecule, + 'small_molecule_file': args.sm_file} if args.GULP: print('converting to GULP format...') @@ -140,4 +278,5 @@ def run_conversion(): if __name__ == '__main__': start_time = time.time() run_conversion() - print("conversion took %s seconds " % np.round((time.time() - start_time), 3)) + print("conversion took %s seconds " % + np.round((time.time() - start_time), 3)) diff --git a/mofa/simulation/cif2lammps/pymatgen_cif2system.py b/mofa/simulation/cif2lammps/pymatgen_cif2system.py index 58ea9ced..a3c0f58c 100644 --- a/mofa/simulation/cif2lammps/pymatgen_cif2system.py +++ b/mofa/simulation/cif2lammps/pymatgen_cif2system.py @@ -39,12 +39,112 @@ def M(vec1, vec2): return eye(3) -PT = ['H', 'He', 'Li', 'Be', 'B', 'C', 'N', 'O', 'F', 'Ne', 'Na', 'Mg', 'Al', 'Si', 'P', 'S', 'Cl', 'Ar', - 'K', 'Ca', 'Sc', 'Ti', 'V', 'Cr', 'Mn', 'Fe', 'Co', 'Ni', 'Cu', 'Zn', 'Ga', 'Ge', 'As', 'Se', 'Br', 'Kr', - 'Rb', 'Sr', 'Y', 'Zr', 'Nb', 'Mo', 'Tc', 'Ru', 'Rh', 'Pd', 'Ag', 'Cd', 'In', 'Sn', 'Sb', 'Te', 'I', 'Xe', - 'Cs', 'Ba', 'Hf', 'Ta', 'W', 'Re', 'Os', 'Ir', 'Pt', 'Au', 'Hg', 'Tl', 'Pb', 'Bi', 'Po', 'At', 'Rn', 'Fr', - 'Ra', 'La', 'Ce', 'Pr', 'Nd', 'Pm', 'Sm', 'Eu', 'Gd', 'Tb', 'Dy', 'Ho', 'Er', 'Tm', 'Yb', 'Lu', 'Ac', 'Th', - 'Pa', 'U', 'Np', 'Pu', 'Am', 'Cm', 'Bk', 'Cf', 'Es', 'Fm', 'Md', 'No', 'Lr', 'FG', 'X'] +PT = [ + 'H', + 'He', + 'Li', + 'Be', + 'B', + 'C', + 'N', + 'O', + 'F', + 'Ne', + 'Na', + 'Mg', + 'Al', + 'Si', + 'P', + 'S', + 'Cl', + 'Ar', + 'K', + 'Ca', + 'Sc', + 'Ti', + 'V', + 'Cr', + 'Mn', + 'Fe', + 'Co', + 'Ni', + 'Cu', + 'Zn', + 'Ga', + 'Ge', + 'As', + 'Se', + 'Br', + 'Kr', + 'Rb', + 'Sr', + 'Y', + 'Zr', + 'Nb', + 'Mo', + 'Tc', + 'Ru', + 'Rh', + 'Pd', + 'Ag', + 'Cd', + 'In', + 'Sn', + 'Sb', + 'Te', + 'I', + 'Xe', + 'Cs', + 'Ba', + 'Hf', + 'Ta', + 'W', + 'Re', + 'Os', + 'Ir', + 'Pt', + 'Au', + 'Hg', + 'Tl', + 'Pb', + 'Bi', + 'Po', + 'At', + 'Rn', + 'Fr', + 'Ra', + 'La', + 'Ce', + 'Pr', + 'Nd', + 'Pm', + 'Sm', + 'Eu', + 'Gd', + 'Tb', + 'Dy', + 'Ho', + 'Er', + 'Tm', + 'Yb', + 'Lu', + 'Ac', + 'Th', + 'Pa', + 'U', + 'Np', + 'Pu', + 'Am', + 'Cm', + 'Bk', + 'Cf', + 'Es', + 'Fm', + 'Md', + 'No', + 'Lr', + 'FG', + 'X'] def nl(string): @@ -102,8 +202,23 @@ def cif_read_pymatgen(filename, charges=False, coplanarity_tolerance=0.1): unit_cell = atoms.get_cell() inv_uc = inv(unit_cell.T) elements = atoms.get_chemical_symbols() - small_skin_metals = ('Ce', 'Pr', 'Nd', 'Pm', 'Sm', 'Eu', 'Gd', 'Tb', 'Dy', 'Ho', 'Er', 'Tm', 'Yb', 'Lu', - 'Ba', 'La') + small_skin_metals = ( + 'Ce', + 'Pr', + 'Nd', + 'Pm', + 'Sm', + 'Eu', + 'Gd', + 'Tb', + 'Dy', + 'Ho', + 'Er', + 'Tm', + 'Yb', + 'Lu', + 'Ba', + 'La') if any(i in elements for i in ('Zn')): skin = 0.30 @@ -118,7 +233,11 @@ def cif_read_pymatgen(filename, charges=False, coplanarity_tolerance=0.1): charge_list = [0.0 for a in atoms] cutoffs = neighborlist.natural_cutoffs(atoms) - NL = neighborlist.NewPrimitiveNeighborList(cutoffs, use_scaled_positions=False, self_interaction=False, skin=skin) # default atom cutoffs work well + NL = neighborlist.NewPrimitiveNeighborList( + cutoffs, + use_scaled_positions=False, + self_interaction=False, + skin=skin) # default atom cutoffs work well NL.build([True, True, True], unit_cell, atoms.get_positions()) G = nx.Graph() @@ -134,7 +253,8 @@ def cif_read_pymatgen(filename, charges=False, coplanarity_tolerance=0.1): jsym = atoms[j].symbol - if (isym not in metals) and (jsym not in metals) and not any(e == 'X' for e in [isym, jsym]): + if (isym not in metals) and (jsym not in metals) and not any( + e == 'X' for e in [isym, jsym]): try: bond = bonds.CovalentBond(struct[i.index], struct[j]) bond_order = bond.get_bond_order() @@ -145,10 +265,23 @@ def cif_read_pymatgen(filename, charges=False, coplanarity_tolerance=0.1): else: bond_order = 0.5 - bond_length = get_distances(i.position, p2=atoms[j].position, cell=unit_cell, pbc=[True, True, True]) + bond_length = get_distances( + i.position, + p2=atoms[j].position, + cell=unit_cell, + pbc=[ + True, + True, + True]) bond_length = np.round(bond_length[1][0][0], 3) - G.add_edge(i.index, j, bond_length=bond_length, bond_order=bond_order, bond_type='', pymatgen_bond_order=bond_order) + G.add_edge( + i.index, + j, + bond_length=bond_length, + bond_order=bond_order, + bond_type='', + pymatgen_bond_order=bond_order) NMG = G.copy() edge_list = list(NMG.edges()) @@ -168,9 +301,11 @@ def cif_read_pymatgen(filename, charges=False, coplanarity_tolerance=0.1): nbor_symbols = [G.nodes[n]['element_symbol'] for n in nbors] nonmetal_nbor_symbols = [n for n in nbor_symbols if n not in metals] - # remove C-M bonds if C is also bonded to carboxylate atoms, these are almost always wrong + # remove C-M bonds if C is also bonded to carboxylate atoms, these are + # almost always wrong for n, nsym in zip(nbors, nbor_symbols): - if isym == 'C' and sorted(nonmetal_nbor_symbols) == ['C', 'O', 'O'] and nsym in metals: + if isym == 'C' and sorted(nonmetal_nbor_symbols) == [ + 'C', 'O', 'O'] and nsym in metals: G.remove_edge(i, n) # intial bond typing, guessed from rounding pymatgen bond orders @@ -185,7 +320,8 @@ def cif_read_pymatgen(filename, charges=False, coplanarity_tolerance=0.1): isym = data['element_symbol'] nbors = list(G.neighbors(i)) nbor_symbols = [G.nodes[n]['element_symbol'] for n in nbors] - nonmetal_nbor_symbols = [n for n in nbor_symbols if n not in metals] + nonmetal_nbor_symbols = [ + n for n in nbor_symbols if n not in metals] CB = nx.cycle_basis(SG) check_cycles = True @@ -217,12 +353,14 @@ def cif_read_pymatgen(filename, charges=False, coplanarity_tolerance=0.1): elif 2.00 <= bond_order < 3.00: bond_order = round(bond_order) - # bonds between two disparate cycles or cycles and non-cycles should have order 1.0 + # bonds between two disparate cycles or cycles and non-cycles + # should have order 1.0 if check_cycles and cyloc is not None: if n not in CB[cyloc]: bond_order = 1.0 - if any(i in metals for i in nbor_symbols) and isym == 'O' and nsym == 'C': + if any( + i in metals for i in nbor_symbols) and isym == 'O' and nsym == 'C': bond_order = 1.5 if nsym in metals: @@ -244,16 +382,19 @@ def cif_read_pymatgen(filename, charges=False, coplanarity_tolerance=0.1): edge_data['bond_type'] = bond_types[bond_order] all_cycles = nx.simple_cycles(nx.to_directed(SG)) - all_cycles = set([tuple(sorted(cy)) for cy in all_cycles if len(cy) > 4]) + all_cycles = set([tuple(sorted(cy)) + for cy in all_cycles if len(cy) > 4]) # # assign aromatic bond orders as 1.5 (in most cases they will be already) for cycle in all_cycles: - # rotate the ring normal vec onto the z-axis to determine coplanarity + # rotate the ring normal vec onto the z-axis to determine + # coplanarity coords = np.array([G.nodes[c]['position'] for c in cycle]) fcoords = np.dot(inv_uc, coords.T).T anchor = fcoords[0] - fcoords = np.array([vec - PBC3DF_sym(anchor, vec)[1] for vec in fcoords]) + fcoords = np.array([vec - PBC3DF_sym(anchor, vec)[1] + for vec in fcoords]) coords = np.dot(unit_cell.T, fcoords.T).T coords -= np.average(coords, axis=0) @@ -281,10 +422,16 @@ def cif_read_pymatgen(filename, charges=False, coplanarity_tolerance=0.1): bond_orders = [G[i][n]['bond_order'] for n in G.neighbors(i)] total_bond_order = np.sum(bond_orders) bond_orders = [str(o) for o in bond_orders] - nbor_symbols = ' '.join([G.nodes[n]['element_symbol'] for n in G.neighbors(i)]) + nbor_symbols = ' '.join([G.nodes[n]['element_symbol'] + for n in G.neighbors(i)]) if isym not in metals and total_bond_order != valencies[isym]: - message = ' '.join([str(isym), 'has total bond order', str(total_bond_order), 'with neighbors', nbor_symbols, 'and bond orders'] + bond_orders) + message = ' '.join([str(isym), + 'has total bond order', + str(total_bond_order), + 'with neighbors', + nbor_symbols, + 'and bond orders'] + bond_orders) warnings.warn(message) elems = atoms.get_chemical_symbols() @@ -301,9 +448,11 @@ def cif_read_pymatgen(filename, charges=False, coplanarity_tolerance=0.1): name0 = sym0 + str(e0) name1 = sym1 + str(e1) - bond_list.append([name0, name1, '.', data['bond_type'], data['bond_length']]) + bond_list.append( + [name0, name1, '.', data['bond_type'], data['bond_length']]) A, B, C = unit_cell.lengths() alpha, beta, gamma = unit_cell.angles() - return elems, names, ccoords, fcoords, charge_list, bond_list, (A, B, C, alpha, beta, gamma), np.asarray(unit_cell).T + return elems, names, ccoords, fcoords, charge_list, bond_list, ( + A, B, C, alpha, beta, gamma), np.asarray(unit_cell).T diff --git a/mofa/simulation/cif2lammps/small_molecule_constants.py b/mofa/simulation/cif2lammps/small_molecule_constants.py index 48cc8f48..b7bd0edf 100644 --- a/mofa/simulation/cif2lammps/small_molecule_constants.py +++ b/mofa/simulation/cif2lammps/small_molecule_constants.py @@ -1,8 +1,10 @@ TraPPE = { 'O2': { 'pair': {'style': 'lj/cut/coul/long', 'vdW': {'O_O2': (0.0974, 3.02), 'O_com': (0.0, 0.0)}, 'charges': {'O_O2': -0.113, 'O_com': 0.226}}, - 'bonds': {('O_O2', 'O_com'): ('harmonic', 100.0, 0.604)}, # molecule should be kept rigid, force constants don't matter - 'angles': {('O_O2', 'O_com', 'O_O2'): ('harmonic', 100.0, 180.0)}, # molecule should be kept rigid, force constants don't matter + # molecule should be kept rigid, force constants don't matter + 'bonds': {('O_O2', 'O_com'): ('harmonic', 100.0, 0.604)}, + # molecule should be kept rigid, force constants don't matter + 'angles': {('O_O2', 'O_com', 'O_O2'): ('harmonic', 100.0, 180.0)}, 'dihedrals': None, 'impropers': None }, @@ -79,10 +81,15 @@ Ions = { 'Cl1': { - 'pair': {'style': 'lj/cut/coul/long', 'vdW': {'Cl_Cl1': (0.22700, 3.51638)}, 'charges': {'Cl_Cl1': -1.0}}, + 'pair': { + 'style': 'lj/cut/coul/long', + 'vdW': { + 'Cl_Cl1': ( + 0.22700, + 3.51638)}, + 'charges': { + 'Cl_Cl1': -1.0}}, 'bonds': None, 'angles': None, 'dihedrals': None, - 'impropers': None - } -} + 'impropers': None}} diff --git a/mofa/simulation/cif2lammps/small_molecule_construction.py b/mofa/simulation/cif2lammps/small_molecule_construction.py index d1378bae..d596ea6e 100644 --- a/mofa/simulation/cif2lammps/small_molecule_construction.py +++ b/mofa/simulation/cif2lammps/small_molecule_construction.py @@ -32,7 +32,10 @@ def add_small_molecules(FF, ff_string): FF.pair_data['special_bonds'] = 'lj/coul 0.0 0.0 1.0' # insert more force fields here if needed else: - raise ValueError('the small molecule force field', ff_string, 'is not defined') + raise ValueError( + 'the small molecule force field', + ff_string, + 'is not defined') SG = FF.system['graph'] SMG = FF.system['SM_graph'] @@ -46,12 +49,19 @@ def add_small_molecules(FF, ff_string): for node, data in SMG.nodes(data=True): # print(node, data) - atoms.append(Atom(data['element_symbol'], data['cartesian_position'])) + atoms.append( + Atom( + data['element_symbol'], + data['cartesian_position'])) atoms.set_cell(FF.system['box']) unit_cell = atoms.get_cell() cutoffs = neighborlist.natural_cutoffs(atoms) - NL = neighborlist.NewPrimitiveNeighborList(cutoffs, use_scaled_positions=False, self_interaction=False, skin=0.10) # shorten the cutoff a bit + NL = neighborlist.NewPrimitiveNeighborList( + cutoffs, + use_scaled_positions=False, + self_interaction=False, + skin=0.10) # shorten the cutoff a bit NL.build([True, True, True], unit_cell, atoms.get_positions()) for i in atoms: @@ -60,9 +70,21 @@ def add_small_molecules(FF, ff_string): for j in nbors: - bond_length = get_distances(i.position, p2=atoms[j].position, cell=unit_cell, pbc=[True, True, True]) + bond_length = get_distances( + i.position, + p2=atoms[j].position, + cell=unit_cell, + pbc=[ + True, + True, + True]) bond_length = np.round(bond_length[1][0][0], 3) - SMG.add_edge(i.index + offset, j + offset, bond_length=bond_length, bond_order='1.0', bond_type='S') + SMG.add_edge( + i.index + offset, + j + offset, + bond_length=bond_length, + bond_order='1.0', + bond_type='S') NMOL = len(list(nx.connected_components(SMG))) print(NMOL, 'small molecules were recovered after bond calculation') @@ -95,7 +117,8 @@ def add_small_molecules(FF, ff_string): mol_flag += 1 comp = sorted(list(comp)) ID_string = sorted([SMG.nodes[n]['element_symbol'] for n in comp]) - ID_string = [(key, len(list(group))) for key, group in groupby(ID_string)] + ID_string = [(key, len(list(group))) + for key, group in groupby(ID_string)] ID_string = ''.join([str(e) for c in ID_string for e in c]) comps.append(ID_string) @@ -108,7 +131,8 @@ def add_small_molecules(FF, ff_string): if ID_string == 'H2O1': SMG.nodes[n]['force_field_type'] = SMG.nodes[n]['element_symbol'] + '_w' else: - SMG.nodes[n]['force_field_type'] = SMG.nodes[n]['element_symbol'] + '_' + ID_string + SMG.nodes[n]['force_field_type'] = SMG.nodes[n]['element_symbol'] + \ + '_' + ID_string # add COM sites where relevant, extend this list as new types are added if ID_string in ('O2', 'N2'): @@ -135,12 +159,22 @@ def add_small_molecules(FF, ff_string): elif ID_string == 'N2': fft = 'N_com' - ndata = {'element_symbol': 'NA', 'mol_flag': mol_flag, 'index': index, 'force_field_type': fft, 'cartesian_position': ccom, - 'fractional_position': fcom, 'charge': 0.0, 'replication': np.array([0.0, 0.0, 0.0]), 'duplicated_version_of': None} + ndata = {'element_symbol': 'NA', + 'mol_flag': mol_flag, + 'index': index, + 'force_field_type': fft, + 'cartesian_position': ccom, + 'fractional_position': fcom, + 'charge': 0.0, + 'replication': np.array([0.0, + 0.0, + 0.0]), + 'duplicated_version_of': None} edata = {'sym_code': None, 'length': None, 'bond_type': None} add_nodes.append([index, ndata]) - add_edges.extend([(index, comp[0], edata), (index, comp[1], edata)]) + add_edges.extend( + [(index, comp[0], edata), (index, comp[1], edata)]) for n, data in add_nodes: SMG.add_node(n, **data) @@ -171,7 +205,8 @@ def add_small_molecules(FF, ff_string): # new_dihedral_types = {} # new_improper_types = {} - for subG, ID_string in zip([SMG.subgraph(c).copy() for c in nx.connected_components(SMG)], comps): + for subG, ID_string in zip([SMG.subgraph(c).copy() + for c in nx.connected_components(SMG)], comps): constants = SM_constants[ID_string] @@ -193,12 +228,14 @@ def add_small_molecules(FF, ff_string): FF.atom_types[fft] = ntypes style = constants['pair']['style'] vdW = constants['pair']['vdW'][fft] - FF.pair_data['params'][FF.atom_types[fft]] = (style, vdW[0], vdW[1]) + FF.pair_data['params'][FF.atom_types[fft]] = ( + style, vdW[0], vdW[1]) FF.pair_data['comments'][FF.atom_types[fft]] = [fft, fft] FF.atom_masses[fft] = mass_key[data['element_symbol']] if 'hybrid' not in FF.pair_data['style'] and style != FF.pair_data['style']: - FF.pair_data['style'] = ' '.join(['hybrid', FF.pair_data['style'], style]) + FF.pair_data['style'] = ' '.join( + ['hybrid', FF.pair_data['style'], style]) elif 'hybrid' in FF.pair_data['style'] and style in FF.pair_data['style']: pass elif 'hybrid' in FF.pair_data['style'] and style not in FF.pair_data['style']: @@ -212,7 +249,8 @@ def add_small_molecules(FF, ff_string): bonds = constants['bonds'] fft_i = SG.nodes[e0]['force_field_type'] fft_j = SG.nodes[e1]['force_field_type'] - # make sure the order corresponds to that in the molecule dictionary + # make sure the order corresponds to that in the molecule + # dictionary bond = tuple(sorted([fft_i, fft_j])) try: @@ -229,17 +267,22 @@ def add_small_molecules(FF, ff_string): used_bonds.append(bond) if 'hybrid' not in FF.bond_data['style'] and style != FF.bond_data['style']: - FF.bond_data['style'] = ' '.join(['hybrid', FF.bond_data['style'], style]) + FF.bond_data['style'] = ' '.join( + ['hybrid', FF.bond_data['style'], style]) elif 'hybrid' in FF.bond_data['style'] and style in FF.bond_data['style']: pass elif 'hybrid' in FF.bond_data['style'] and style not in FF.bond_data['style']: FF.bond_data['style'] += ' ' + style if ty in FF.bond_data['all_bonds']: - FF.bond_data['count'] = (FF.bond_data['count'][0] + 1, FF.bond_data['count'][1] + 1) + FF.bond_data['count'] = ( + FF.bond_data['count'][0] + 1, + FF.bond_data['count'][1] + 1) FF.bond_data['all_bonds'][ty].append((e0, e1)) else: - FF.bond_data['count'] = (FF.bond_data['count'][0] + 1, FF.bond_data['count'][1] + 1) + FF.bond_data['count'] = ( + FF.bond_data['count'][0] + 1, + FF.bond_data['count'][1] + 1) FF.bond_data['all_bonds'][ty] = [(e0, e1)] except KeyError: @@ -267,20 +310,24 @@ def add_small_molecules(FF, ff_string): try: style = angles[angle][0] - FF.angle_data['count'] = (FF.angle_data['count'][0] + 1, FF.angle_data['count'][1]) + FF.angle_data['count'] = ( + FF.angle_data['count'][0] + 1, + FF.angle_data['count'][1]) if angle not in used_angles: ty = ty + 1 new_angle_types[angle] = ty - FF.angle_data['count'] = (FF.angle_data['count'][0], FF.angle_data['count'][1] + 1) + FF.angle_data['count'] = ( + FF.angle_data['count'][0], FF.angle_data['count'][1] + 1) FF.angle_data['params'][ty] = list(angles[angle]) FF.angle_data['comments'][ty] = list(angle) used_angles.append(angle) if 'hybrid' not in FF.angle_data['style'] and style != FF.angle_data['style']: - FF.angle_data['style'] = ' '.join(['hybrid', FF.angle_data['style'], style]) + FF.angle_data['style'] = ' '.join( + ['hybrid', FF.angle_data['style'], style]) elif 'hybrid' in FF.angle_data['style'] and style in FF.angle_data['style']: pass elif 'hybrid' in FF.angle_data['style'] and style not in FF.angle_data['style']: @@ -296,8 +343,12 @@ def add_small_molecules(FF, ff_string): # add new dihedrals - FF.bond_data['count'] = (FF.bond_data['count'][0], len(FF.bond_data['params'])) - FF.angle_data['count'] = (FF.angle_data['count'][0], len(FF.angle_data['params'])) + FF.bond_data['count'] = ( + FF.bond_data['count'][0], len( + FF.bond_data['params'])) + FF.angle_data['count'] = ( + FF.angle_data['count'][0], len( + FF.angle_data['params'])) if 'tip4p' in FF.pair_data['style']: @@ -317,7 +368,8 @@ def add_small_molecules(FF, ff_string): FF.pair_data['H2O_angle_type'] = ty if 'long' in FF.pair_data['style']: - FF.pair_data['M_site_dist'] = 0.1546 # only TIP4P/2005 is implemented + # only TIP4P/2005 is implemented + FF.pair_data['M_site_dist'] = 0.1546 elif 'cut' in FF.pair_data['style'] and ff_string == 'TIP4P_2005_cutoff': FF.pair_data['M_site_dist'] = 0.1546 elif 'cut' in FF.pair_data['style'] and ff_string == 'TIP4P_cutoff': @@ -327,26 +379,37 @@ def add_small_molecules(FF, ff_string): def update_potential(potential_data, new_potential_params, potential_coeff): write_instyles = False - add_styles = set([new_potential_params[ty]['style'] for ty in new_potential_params]) + add_styles = set([new_potential_params[ty]['style'] + for ty in new_potential_params]) for ABS in add_styles: if ABS not in potential_data['style'] and 'hybrid' in potential_data['style']: potential_data['style'] = potential_data['style'] + ' ' + ABS write_instyles = True if ABS not in potential_data['style'] and 'hybrid' not in potential_data['style']: - potential_data['style'] = 'hybrid ' + potential_data['style'] + ' ' + ABS + potential_data['style'] = 'hybrid ' + \ + potential_data['style'] + ' ' + ABS write_instyles = True else: pass if write_instyles: - instyles = {ty: ' ' + new_potential_params[ty]['style'] for ty in new_potential_params} + instyles = { + ty: ' ' + + new_potential_params[ty]['style'] for ty in new_potential_params} else: instyles = {ty: '' for ty in new_potential_params} potential_data['infile_add_lines'] = [] for ty, data in new_potential_params.items(): strparams = ' '.join([str(p) for p in data['params']]) - potential_data['infile_add_lines'].append(potential_coeff + str(ty) + instyles[ty] + ' ' + strparams + ' ' + data['comments']) + potential_data['infile_add_lines'].append( + potential_coeff + + str(ty) + + instyles[ty] + + ' ' + + strparams + + ' ' + + data['comments']) def include_molecule_file(FF, maxIDs, add_molecule): @@ -383,18 +446,23 @@ def include_molecule_file(FF, maxIDs, add_molecule): if angle_params is not None: update_potential(FF.angle_data, angle_params, 'angle_coeff ') if dihedral_params is not None: - update_potential(FF.dihedral_params, dihedral_params, 'dihedral_coeff ') + update_potential( + FF.dihedral_params, + dihedral_params, + 'dihedral_coeff ') if improper_params is not None: update_potential(FF.improper_data, improper_params, 'improper_coeff ') infile_add_lines = ['molecule ' + ' '.join(molnames)] for atom in mass_dict: - infile_add_lines.append('mass ' + str(atom) + ' ' + str(mass_dict[atom])) + infile_add_lines.append('mass ' + + str(atom) + ' ' + str(mass_dict[atom])) seed0 = randint(1, 10000) seed1 = randint(1, 10000) if N > 0: - create_line = ' '.join([str(N), str(seed0), 'NULL', 'mol', molnames[0], str(seed1), 'units', 'box']) + create_line = ' '.join( + [str(N), str(seed0), 'NULL', 'mol', molnames[0], str(seed1), 'units', 'box']) infile_add_lines.append('create_atoms 0 random ' + create_line) return molfile, infile_add_lines, extra_types @@ -422,7 +490,8 @@ def read_small_molecule_file(sm_file, system): ind = max_ind + 1 if fm not in ('pdb', 'xyz', 'cif'): - raise ValueError('only xyz and RASPA pdb formats are supported for small molecule files') + raise ValueError( + 'only xyz and RASPA pdb formats are supported for small molecule files') if fm == 'pdb': print('assuming small molecule file is in RASPA pdb format, if not, too bad...') @@ -435,9 +504,20 @@ def read_small_molecule_file(sm_file, system): for atom in atoms: - SMG.add_node(ind, element_symbol=atom.symbol, mol_flag='1', index=ind, force_field_type='', cartesian_position=atom.position, - fractional_position=atom.scaled_position, charge=0.0, replication=np.array([0.0, 0.0, 0.0]), duplicated_version_of=None) + SMG.add_node(ind, + element_symbol=atom.symbol, + mol_flag='1', + index=ind, + force_field_type='', + cartesian_position=atom.position, + fractional_position=atom.scaled_position, + charge=0.0, + replication=np.array([0.0, + 0.0, + 0.0]), + duplicated_version_of=None) ind += 1 - system['SM_graph'] = nx.compose(SMG, system['SM_graph']) # don't want to overwrite extra framework species already in the cif + # don't want to overwrite extra framework species already in the cif + system['SM_graph'] = nx.compose(SMG, system['SM_graph']) diff --git a/mofa/simulation/cif2lammps/write_GULP_inputs.py b/mofa/simulation/cif2lammps/write_GULP_inputs.py index 79333eef..1355a1c1 100644 --- a/mofa/simulation/cif2lammps/write_GULP_inputs.py +++ b/mofa/simulation/cif2lammps/write_GULP_inputs.py @@ -1,5 +1,6 @@ from __future__ import print_function -from .cif2system import initialize_system, replication_determination # , duplicate_system, write_cif_from_system +# , duplicate_system, write_cif_from_system +from .cif2system import initialize_system, replication_determination from . import atomic_data import os import numpy as np @@ -49,14 +50,23 @@ def isfloat(value): def GULP_inputs(args): - gulp_bond_types = {0.25: 'quarter', 0.5: 'half', 1.0: '', 1.5: 'resonant', 2.0: '', 3.0: ''} + gulp_bond_types = { + 0.25: 'quarter', + 0.5: 'half', + 1.0: '', + 1.5: 'resonant', + 2.0: '', + 3.0: ''} cifname, force_field, outdir, charges, replication, noautobond = args - FF_args = {'FF_parameters': UFF4MOF_atom_parameters, 'bond_orders': UFF4MOF_bond_orders_0} + FF_args = { + 'FF_parameters': UFF4MOF_atom_parameters, + 'bond_orders': UFF4MOF_bond_orders_0} cutoff = 12.5 system = initialize_system(cifname, charges=charges) - system, replication = replication_determination(system, replication, cutoff) + system, replication = replication_determination( + system, replication, cutoff) FF = force_field(system, cutoff, FF_args) FF.compile_force_field(charges=charges) @@ -96,7 +106,8 @@ def GULP_inputs(args): gin.write('\n') gin.write('\n') - bonds = [(b, ty) for ty in FF.bond_data['all_bonds'] for b in FF.bond_data['all_bonds'][ty]] + bonds = [(b, ty) for ty in FF.bond_data['all_bonds'] + for b in FF.bond_data['all_bonds'][ty]] bonds.sort(key=lambda x: x[0][0]) if noautobond: diff --git a/mofa/simulation/cif2lammps/write_lammps_data.py b/mofa/simulation/cif2lammps/write_lammps_data.py index 592a1511..a39bbecf 100644 --- a/mofa/simulation/cif2lammps/write_lammps_data.py +++ b/mofa/simulation/cif2lammps/write_lammps_data.py @@ -1,5 +1,6 @@ from __future__ import print_function -from .cif2system import initialize_system, replication_determination # , write_cif_from_system +# , write_cif_from_system +from .cif2system import initialize_system, replication_determination from .small_molecule_construction import add_small_molecules, include_molecule_file, read_small_molecule_file from . import atomic_data import os @@ -47,15 +48,21 @@ def lammps_inputs(args): # add more forcefields here as they are created if ff_string == 'UFF4MOF': - FF_args = {'FF_parameters': UFF4MOF_atom_parameters, 'bond_orders': UFF4MOF_bond_orders_0} + FF_args = { + 'FF_parameters': UFF4MOF_atom_parameters, + 'bond_orders': UFF4MOF_bond_orders_0} cutoff = 12.50 mixing_rules = 'shift yes mix geometric' elif ff_string == 'UFF': - FF_args = {'FF_parameters': UFF_atom_parameters, 'bond_orders': UFF_bond_orders_0} + FF_args = { + 'FF_parameters': UFF_atom_parameters, + 'bond_orders': UFF_bond_orders_0} cutoff = 12.50 mixing_rules = 'shift yes mix geometric' elif ff_string == 'Dreiding': - FF_args = {'FF_parameters': Dreiding_atom_parameters, 'bond_orders': Dreiding_bond_orders_0} + FF_args = { + 'FF_parameters': Dreiding_atom_parameters, + 'bond_orders': Dreiding_bond_orders_0} cutoff = 12.50 mixing_rules = 'shift yes mix arithmetic' elif ff_string == 'MZHB': @@ -67,13 +74,17 @@ def lammps_inputs(args): cutoff = 12.50 mixing_rules = 'shift yes mix arithmetic' - system = initialize_system(cifname, charges=charges, read_pymatgen=read_pymatgen) + system = initialize_system( + cifname, + charges=charges, + read_pymatgen=read_pymatgen) if sm_file is not None: direc = cifname.split(os.sep)[0] read_small_molecule_file(direc + os.sep + sm_file, system) - system, replication = replication_determination(system, replication, cutoff) + system, replication = replication_determination( + system, replication, cutoff) print('system initialized...') @@ -84,10 +95,12 @@ def lammps_inputs(args): add_small_molecules(FF, sm_ff_string) else: if len(FF.system['SM_graph'].nodes()) != 0: - warnings.warn('extra-framework molecules detected, but no small molecule force field is specified!') + warnings.warn( + 'extra-framework molecules detected, but no small molecule force field is specified!') # write_cif_from_system(system, 'check.cif') - first_line = "Created by Ryther's code on " + str(datetime.datetime.now()) + first_line = "Created by Ryther's code on " + \ + str(datetime.datetime.now()) + ", modified by xyan11@uic.edu. " SG = FF.system['graph'] N_atoms, ty_atoms = (len(SG.nodes()), len(FF.atom_types)) @@ -107,7 +120,8 @@ def lammps_inputs(args): ty_impropers = None if replication != '': - suffix = ''.join(cifname.split('/')[-1].split('.')[0:-1]) + '_' + replication + suffix = ''.join(cifname.split( + '/')[-1].split('.')[0:-1]) + '_' + replication else: suffix = ''.join(cifname.split('/')[-1].split('.')[0:-1]) @@ -118,11 +132,14 @@ def lammps_inputs(args): maxIDs = (ty_atoms, ty_bonds, ty_angles, ty_dihedrals, ty_impropers) if add_molecule is not None: - molfile, infile_add_lines, extra_types = include_molecule_file(FF, maxIDs, add_molecule) + molfile, infile_add_lines, extra_types = include_molecule_file( + FF, maxIDs, add_molecule) with open(outdir + os.sep + 'mol.' + suffix, 'w') as MF: MF.write(molfile) a, b, c, alpha, beta, gamma = system['box'] + cifbox = "# a, b, c, alpha, beta, gamma: " + "%.10f" % a + ", " + "%.10f" % b + ", " + "%.10f" % c + \ + ", " + "%.10f" % alpha + ", " + "%.10f" % beta + ", " + "%.10f" % gamma + " $$$atoms$$$" lx = np.round(a, 8) xy = np.round(b * np.cos(math.radians(gamma)), 8) xz = np.round(c * np.cos(math.radians(beta)), 8) @@ -178,7 +195,8 @@ def lammps_inputs(args): # type data.write(' {:<3}'.format(aty)) - format_string = ' '.join(['{:{w}.{p}f}' if not np.issubdtype(x, np.integer) else '{:{w}}' for x in params[1:]]) + format_string = ' '.join(['{:{w}.{p}f}' if not np.issubdtype( + x, np.integer) else '{:{w}}' for x in params[1:]]) data.write(format_string.format(*params[1:], w=12, p=5)) data.write(' '.join([' # '] + comment)) data.write('\n') @@ -208,7 +226,8 @@ def lammps_inputs(args): data.write('{:<20}'.format(params[0])) try: - format_string = ' '.join(['{:{w}.{p}f}' if not np.issubdtype(x, np.integer) else '{:{w}}' for x in params[1:]]) + format_string = ' '.join(['{:{w}.{p}f}' if not np.issubdtype( + x, np.integer) else '{:{w}}' for x in params[1:]]) except TypeError: format_string = ' '.join(['{:{w}}' for x in params[1:]]) @@ -220,7 +239,8 @@ def lammps_inputs(args): data.write(' {:<3}'.format(bty)) try: - format_string = ' '.join(['{:{w}.{p}f}' if not np.issubdtype(x, np.integer) else '{:{w}}' for x in params[1:]]) + format_string = ' '.join(['{:{w}.{p}f}' if not np.issubdtype( + x, np.integer) else '{:{w}}' for x in params[1:]]) except TypeError: format_string = ' '.join(['{:{w}}' for x in params[1:]]) @@ -246,7 +266,8 @@ def lammps_inputs(args): data.write('{:<20}'.format(params[0])) try: - format_string = ' '.join(['{:{w}.{p}f}' if not np.issubdtype(x, np.integer) else '{:{w}}' for x in params[1:]]) + format_string = ' '.join(['{:{w}.{p}f}' if not np.issubdtype( + x, np.integer) else '{:{w}}' for x in params[1:]]) except TypeError: format_string = ' '.join(['{:{w}}' for x in params[1:]]) @@ -258,7 +279,8 @@ def lammps_inputs(args): data.write(' {:<3}'.format(aty)) try: - format_string = ' '.join(['{:{w}.{p}f}' if not np.issubdtype(x, np.integer) else '{:{w}}' for x in params[1:]]) + format_string = ' '.join(['{:{w}.{p}f}' if not np.issubdtype( + x, np.integer) else '{:{w}}' for x in params[1:]]) except TypeError: format_string = ' '.join(['{:{w}}' for x in params[1:]]) @@ -283,14 +305,16 @@ def lammps_inputs(args): data.write(' {:<3}'.format(dty)) # style needs to be written for hybrid data.write('{:<20}'.format(params[0])) - format_string = ' '.join(['{:{w}.{p}f}' if not np.issubdtype(x, np.integer) else '{:{w}}' for x in params[1:]]) + format_string = ' '.join(['{:{w}.{p}f}' if not np.issubdtype( + x, np.integer) else '{:{w}}' for x in params[1:]]) data.write(format_string.format(*params[1:], w=12, p=5)) data.write(' '.join([' # '] + comment)) data.write('\n') else: # type data.write(' {:<3}'.format(dty)) - format_string = ' '.join(['{:{w}.{p}f}' if not np.issubdtype(x, np.integer) else '{:{w}}' for x in params[1:]]) + format_string = ' '.join(['{:{w}.{p}f}' if not np.issubdtype( + x, np.integer) else '{:{w}}' for x in params[1:]]) data.write(format_string.format(*params[1:], w=12, p=5)) data.write(' '.join([' # '] + comment)) data.write('\n') @@ -312,20 +336,22 @@ def lammps_inputs(args): data.write(' {:<3}'.format(ity)) # style needs to be written for hybrid data.write('{:<20}'.format(params[0])) - format_string = ' '.join(['{:{w}.{p}f}' if not np.issubdtype(x, np.integer) else '{:{w}}' for x in params[1:]]) + format_string = ' '.join(['{:{w}.{p}f}' if not np.issubdtype( + x, np.integer) else '{:{w}}' for x in params[1:]]) data.write(format_string.format(*params[1:], w=12, p=5)) data.write(' '.join([' # '] + comment)) data.write('\n') else: # type data.write(' {:<3}'.format(ity)) - format_string = ' '.join(['{:{w}.{p}f}' if not np.issubdtype(x, np.integer) else '{:{w}}' for x in params[1:]]) + format_string = ' '.join(['{:{w}.{p}f}' if not np.issubdtype( + x, np.integer) else '{:{w}}' for x in params[1:]]) data.write(format_string.format(*params[1:], w=12, p=5)) data.write(' '.join([' # '] + comment)) data.write('\n') data.write('\n') - data.write('Atoms\n') + data.write('Atoms ' + cifbox + '\n') data.write('\n') total_charge = 0.0 @@ -343,20 +369,30 @@ def lammps_inputs(args): atom_cif_labels.append(cif_label) total_charge += charge pos = [np.round(v, 8) for v in atom_data['cartesian_position']] - - data.write('{:>5} {:<5} {:<5} {:8.5f} {:12.5f} {:12.5f} {:12.5f} # {:>5}'.format(index, - atom_data['mol_flag'], - lammps_type, - charge, - pos[0], - pos[1], - pos[2], - cif_label)) + frac_pos = [np.round(v, 8) + for v in atom_data['fractional_position']] + + data.write( + '{:>5} {:<5} {:<5} {:8.5f} {:12.5f} {:12.5f} {:12.5f} # {:>5} {:12.8f} {:12.8f} {:12.8f}'.format( + index, + atom_data['mol_flag'], + lammps_type, + charge, + pos[0], + pos[1], + pos[2], + cif_label, + frac_pos[0], + frac_pos[1], + frac_pos[2], + )) data.write('\n') id2label = dict(zip(atom_ids, atom_cif_labels)) if charges and abs(total_charge) > 0.001: - warnings.warn('There is a potentially significant net charge of ' + str(total_charge)) + warnings.warn( + 'There is a potentially significant net charge of ' + + str(total_charge)) data.write('\n') data.write('Bonds\n') @@ -461,7 +497,12 @@ def lammps_inputs(args): read_data_append_string = '' if add_molecule is not None: read_data_append_string = ' ' - et_strings = ['extra/atom/types', 'extra/bond/types', 'extra/angle/types', 'extra/dihedral/types', 'extra/improper/types'] + et_strings = [ + 'extra/atom/types', + 'extra/bond/types', + 'extra/angle/types', + 'extra/dihedral/types', + 'extra/improper/types'] for st, et in zip(et_strings, extra_types): if et is not None: read_data_append_string += st + ' ' + str(et) + ' ' @@ -501,15 +542,21 @@ def lammps_inputs(args): if 'hybrid' not in pair_style: if 'tip4p' in pair_style: - add_arg = ' ' + ' '.join([ - str(FF.pair_data['O_type']), str(FF.pair_data['H_type']), - str(FF.pair_data['H2O_bond_type']), str(FF.pair_data['H2O_angle_type']), - str(FF.pair_data['M_site_dist']), '12.5', '12.5' - ]) + add_arg = ' ' + ' '.join([str(FF.pair_data['O_type']), + str(FF.pair_data['H_type']), + str(FF.pair_data['H2O_bond_type']), + str(FF.pair_data['H2O_angle_type']), + str(FF.pair_data['M_site_dist']), + '12.5', + '12.5']) FF.pair_data['style'] = FF.pair_data['style'] + add_arg infile.write('pair_style ' + FF.pair_data['style'] + '\n') else: - infile.write('pair_style ' + FF.pair_data['style'] + ' ' + str(FF.cutoff) + '\n') + infile.write('pair_style ' + + FF.pair_data['style'] + + ' ' + + str(FF.cutoff) + + '\n') else: style_list = FF.pair_data['style'].split() @@ -549,7 +596,8 @@ def lammps_inputs(args): infile.write('pair_modify ' + mixing_rules + '\n') infile.write('special_bonds ' + sb + '\n') - # use ewald summation for long range solver unless using pair_style lj/cut/tip4p/long + # use ewald summation for long range solver unless using pair_style + # lj/cut/tip4p/long if 'long' in pair_style and 'tip4p' not in pair_style: infile.write('kspace_style ewald 1.0e-5\n') elif 'tip4p/long' in pair_style: @@ -564,15 +612,20 @@ def lammps_inputs(args): params0 = FF.pair_data['params'][aty0] if len(params0) > 3: - raise ValueError('pair_style hybrid is only supported for lj type vdw interactions (two numeric parameters).') + raise ValueError( + 'pair_style hybrid is only supported for lj type vdw interactions (two numeric parameters).') - params0 = [np.round(x, 6) if isfloat(x) else x for x in params0] + params0 = [np.round(x, 6) if isfloat(x) + else x for x in params0] style0, eps0, sig0 = params0 - comment = ' # ' + ' ' + ' '.join(FF.pair_data['comments'][aty0]) + comment = ' # ' + ' ' + \ + ' '.join(FF.pair_data['comments'][aty0]) line = ['pair_coeff', aty0, aty0, style0, eps0, sig0, comment] - infile.write('{:12} {:<3} {:<3} {:20} {:10.6f} {:10.6f} {:<20}'.format(*line)) + infile.write( + '{:12} {:<3} {:<3} {:20} {:10.6f} {:10.6f} {:<20}'.format( + *line)) infile.write('\n') infile.write('\n') @@ -580,12 +633,14 @@ def lammps_inputs(args): for (aty0, aty1) in combinations(FF.pair_data['params'], 2): params0 = FF.pair_data['params'][aty0] - params0 = [np.round(x, 6) if isfloat(x) else x for x in params0] + params0 = [np.round(x, 6) if isfloat(x) + else x for x in params0] style0, eps0, sig0 = params0 comment0 = FF.pair_data['comments'][aty0][0] params1 = FF.pair_data['params'][aty1] - params1 = [np.round(x, 6) if isfloat(x) else x for x in params1] + params1 = [np.round(x, 6) if isfloat(x) + else x for x in params1] style1, eps1, sig1 = params1 comment1 = FF.pair_data['comments'][aty1][0] @@ -593,7 +648,8 @@ def lammps_inputs(args): # current logic is to use the longer style, this actually works well when using lj/cut for # framework atoms and lj/cut + charge interactions for framework/molecule and molecule/molecule - # interactions. This is mostly what I use pair_style hybrid for. + # interactions. This is mostly what I use pair_style hybrid + # for. style = style0 if len(style0) > len(style1) else style1 if 'geometric' in mixing_rules: @@ -606,19 +662,33 @@ def lammps_inputs(args): pair_eps = np.round(np.sqrt(eps0 * eps1), 6) pair_sig = np.round(np.sqrt(sig0 * sig1), 6) - line = ['pair_coeff', aty0, aty1, style, pair_eps, pair_sig, comments] - - infile.write('{:12} {:<3} {:<3} {:20} {:10.6f} {:10.6f} {:<20}'.format(*line)) + line = [ + 'pair_coeff', + aty0, + aty1, + style, + pair_eps, + pair_sig, + comments] + + infile.write( + '{:12} {:<3} {:<3} {:20} {:10.6f} {:10.6f} {:<20}'.format( + *line)) infile.write('\n') infile.write('\n') infile.write('box tilt large\n') - infile.write('read_data ' + data_name + read_data_append_string + '\n\n') + infile.write( + 'read_data ' + + data_name + + read_data_append_string + + '\n\n') if add_molecule is not None: if 'TIP4P' in add_molecule[1]: - group_line = 'group H2O type ' + str(FF.pair_data['O_type']) + ' ' + str(FF.pair_data['H_type']) + '\n' + group_line = 'group H2O type ' + \ + str(FF.pair_data['O_type']) + ' ' + str(FF.pair_data['H_type']) + '\n' shake_line = 'fix H2O_shake H2O shake 0.0001 50 0 b ' + \ str(FF.pair_data['H2O_bond_type']) + ' a ' + str(FF.pair_data['H2O_angle_type']) + ' mol H2O_mol\n' infile.write(group_line) @@ -627,7 +697,8 @@ def lammps_inputs(args): if sm_ff_string is not None: if 'TIP4P' in sm_ff_string: - group_line = 'group H2O type ' + str(FF.pair_data['O_type']) + ' ' + str(FF.pair_data['H_type']) + '\n' + group_line = 'group H2O type ' + \ + str(FF.pair_data['O_type']) + ' ' + str(FF.pair_data['H_type']) + '\n' shake_line = 'fix H2O_shake H2O shake 0.0001 50 0 b ' + \ str(FF.pair_data['H2O_bond_type']) + ' a ' + \ str(FF.pair_data['H2O_angle_type']) + ' t ' + \ diff --git a/mofa/simulation/cif2lammps/write_molecule_files.py b/mofa/simulation/cif2lammps/write_molecule_files.py index 8896858d..d4e8e122 100644 --- a/mofa/simulation/cif2lammps/write_molecule_files.py +++ b/mofa/simulation/cif2lammps/write_molecule_files.py @@ -24,7 +24,8 @@ def water(last_atom_ID, last_bond_ID, last_angle_ID, model='TIP4P_cutoff'): } LJ_dict = { - # LAMMPS has a special TIP4P pair_style that automatically adds the M site + # LAMMPS has a special TIP4P pair_style that automatically adds the M + # site 'TIP4P_cutoff': {ID_O: ('lj/cut/tip4p/cut', 0.15500, 3.15360), ID_H: ('lj/cut/tip4p/cut', 0.0, 0.0), 'style': 'lj/cut/tip4p/cut', @@ -72,7 +73,8 @@ def water(last_atom_ID, last_bond_ID, last_angle_ID, model='TIP4P_cutoff'): if 'TIP4P' in model: - molfile = dedent(""" # Water molecule. useable for TIP3P or TIP4P in LAMMPS. + molfile = dedent( + """ # Water molecule. useable for TIP3P or TIP4P in LAMMPS. 3 atoms 2 bonds @@ -121,11 +123,13 @@ def water(last_atom_ID, last_bond_ID, last_angle_ID, model='TIP4P_cutoff'): 1 {BT} {BT} {AT} 2 {BT} {BT} {AT} -3 {BT} {BT} {AT}""".format(**locals())).strip() +3 {BT} {BT} {AT}""".format( + **locals())).strip() if 'TIP3P' in model: - molfile = dedent(""" # Water molecule. useable for TIP3P or TIP4P in LAMMPS. + molfile = dedent( + """ # Water molecule. useable for TIP3P or TIP4P in LAMMPS. 3 atoms 2 bonds @@ -156,11 +160,13 @@ def water(last_atom_ID, last_bond_ID, last_angle_ID, model='TIP4P_cutoff'): Angles -1 {AT} 2 1 3""".format(**locals())).strip() +1 {AT} 2 1 3""".format( + **locals())).strip() mass_dict = {ID_O: 15.9994, ID_H: 1.00794} molnames = ('H2O_mol', 'H2O.txt') extra_types = (2, 1, 1, None, None) - return molfile, LJ_params, bond_params, angle_params, molnames, mass_dict, M_site_dist_dict[model], extra_types + return molfile, LJ_params, bond_params, angle_params, molnames, mass_dict, M_site_dist_dict[ + model], extra_types diff --git a/mofa/simulation/cif2lammps/zeoliteFFs_construction.py b/mofa/simulation/cif2lammps/zeoliteFFs_construction.py index 9a7b5fea..e9b80da8 100644 --- a/mofa/simulation/cif2lammps/zeoliteFFs_construction.py +++ b/mofa/simulation/cif2lammps/zeoliteFFs_construction.py @@ -34,14 +34,19 @@ def type_atoms(self): elif element_symbol == 'O': ty = 'O' else: - raise ValueError('No Nicholas type identified for ' + element_symbol + 'with neighbors ' + ' '.join(nbor_symbols)) + raise ValueError( + 'No Nicholas type identified for ' + + element_symbol + + 'with neighbors ' + + ' '.join(nbor_symbols)) types.append((ty, element_symbol, mass)) SG.nodes[name]['force_field_type'] = ty 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) @@ -80,7 +85,11 @@ def pair_parameters(self, charges=True): else: pass - self.pair_data = {'params': params, 'style': style, 'special_bonds': sb, 'comments': comments} + self.pair_data = { + 'params': params, + 'style': style, + 'special_bonds': sb, + 'comments': comments} def bond_parameters(self, bond): @@ -94,7 +103,8 @@ def bond_parameters(self, bond): k_ij = 537.31 / 2.0 r_ij = 1.620 else: - raise ValueError('There is a non Si-O bond, which is not yet parametrized for Nicholas') + raise ValueError( + 'There is a non Si-O bond, which is not yet parametrized for Nicholas') return ('harmonic', k_ij, r_ij) @@ -165,7 +175,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): @@ -221,7 +238,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):