diff --git a/qforce/forcefield.py b/qforce/forcefield.py index 2c71264..ba05620 100644 --- a/qforce/forcefield.py +++ b/qforce/forcefield.py @@ -1,9 +1,10 @@ import numpy as np +import string # from .elements import ATOM_SYM, ATOMMASS from .molecule.non_bonded import calc_sigma_epsilon from .forces import convert_to_inversion_rb -from .misc import LOGO_SEMICOL +from .misc import get_logo class ForceField(): @@ -24,6 +25,11 @@ def __init__(self, job_name, config, mol, neighbors, exclude_all=[]): self.exclusions = self.make_exclusions(mol.non_bonded, neighbors, exclude_all) self.pairs = self.make_pairs(neighbors, mol.non_bonded) + if config.ff.output_software == 'gromacs': + self.write_ff = self.write_gromacs + elif config.ff.output_software == 'amber': + self.write_ff = self.write_amber + if self.polar: self.polar_title = '_polar' else: @@ -75,7 +81,7 @@ def write_gro(self, directory, coords, alpha_map, box=[20., 20., 20.]): def write_itp(self, mol, directory): with open(f"{directory}/{self.mol_name}_qforce{self.polar_title}.itp", "w") as itp: - itp.write(LOGO_SEMICOL) + itp.write(get_logo(';')) self.write_itp_atoms_and_molecule(itp, mol.non_bonded) if self.polar: self.write_itp_polarization(itp, mol.non_bonded) @@ -128,7 +134,8 @@ def write_itp_atoms_and_molecule(self, itp, non_bonded): itp.write("; name at_num mass charge type sigma epsilon\n") for lj_type, lj_params in gro_atomtypes.items(): - itp.write(f'{lj_type:>8} {non_bonded.lj_atomic_number[lj_type]:>7} {0:>8.4f} {0:>8.4f} {"A":>5} ') + itp.write(f'{lj_type:>8} {non_bonded.lj_atomic_number[lj_type]:>7} {0:>8.4f} ' + f'{0:>8.4f} {"A":>5} ') itp.write(f'{lj_params[0]:>12.5e} {lj_params[1]:>12.5e}\n') if self.polar: @@ -387,89 +394,198 @@ def set_charge(self, non_bonded): q[list(non_bonded.alpha_map.keys())] += 8 return q - # bohr2nm = 0.052917721067 - # if polar: - # alphas = qm.alpha*bohr2nm**3 - # drude = {} - # n_drude = 1 - # ff.atom_types.append(["DP", 0, 0, "S", 0, 0]) - - # for i, alpha in enumerate(alphas): - # if alpha > 0: - # drude[i] = mol.topo.n_atoms+n_drude - # ff.atoms[i][6] += 8 - # # drude atoms - # ff.atoms.append([drude[i], 'DP', 2, 'MOL', f'D{atoms[i]}', - # i+1, -8., 0.]) - # ff.coords.append(ff.coords[i]) - # # polarizability - # ff.polar.append([i+1, drude[i], 1, alpha]) - # n_drude += 1 - # ff.natom = len(ff.atoms) - # for i, alpha in enumerate(alphas): - # if alpha > 0: - # # exclusions for balancing the drude particles - # for j in (mol.topo.neighbors[self.n_excl-2][i] + - # mol.topo.neighbors[self.n_excl-1][i]): - # if alphas[j] > 0: - # ff.exclu[drude[i]-1].extend([drude[j]]) - # for j in mol.topo.neighbors[self.n_excl-1][i]: - # ff.exclu[drude[i]-1].extend([j+1]) - # ff.exclu[drude[i]-1].sort() - # # thole polarizability - # for neigh in [mol.topo.neighbors[n][i] for n in range(self.n_excl)]: - # for j in neigh: - # if i < j and alphas[j] > 0: - # ff.thole.append([i+1, drude[i], j+1, drude[j], "2", 2.6, alpha, - # alphas[j]]) - - # def read_itp(self, itp_file): - # with open(itp_file, "r") as itp: - # in_section = [] - # bond_atoms, bond_r0, bond_k = [], [], [] - - # for line in itp: - # low_line = line.lower().strip().replace(" ", "") - # unsplit = line - # line = line.split() - # if low_line == "" or low_line[0] == ";": - # continue - # elif "[" in low_line and "]" in low_line: - # open_bra = low_line.index("[") + 1 - # close_bra = low_line.index("]") - # in_section = low_line[open_bra:close_bra] - # elif in_section == "atomtypes": - # self.atom_types.append([line[0], float(line[1]), - # float(line[2]), line[3], - # float(line[4]), float(line[5])]) - # self.atype.append(line[0]) - # self.c6.append(line[4]) - # self.c12.append(line[5]) - # elif in_section == "moleculetype": - # self.mol_type = line[0] - # elif in_section == "atoms": - # self.atoms.append([int(line[0]), line[1], int(line[2]), - # line[3], line[4], line[5], float(line[6]), - # float(line[7])]) - # self.natom += 1 - # elif in_section == "bonds": - # bond_atoms = (line[0:2]) - # bond_r0 = (float(line[3]) * 10) - # bond_k = (float(line[4]) / 100) - # self.bond.append([bond_atoms, bond_r0, bond_k]) - # self.bonds.append(unsplit) - # elif in_section == "angles": - # angle_atoms = (line[0:3]) - # angle_theta0 = float(line[4]) - # angle_k = float(line[5]) - # self.angle.append([angle_atoms, angle_theta0, angle_k]) - # self.angles.append(unsplit) - # elif in_section == "dihedrals": - # self.dihedrals.append(unsplit) - # if line[4] == "3": - # self.rbdihed.append([line[0:4], line[4:]]) - # elif line[4] == "2": - # self.idihed.append([line[0:4], line[4:]]) - - # elif in_section == "pairs": - # self.pairs.append(sorted([int(line[0]), int(line[1])])) + def write_amber(self, directory, mol, coords): + if self.urey == True: + print(f'NOTE: AMBERs parameters file does not support Urey-Bradley terms.') + print(f' Set urey to False! \n') + j2cal = 0.2390057361 + self.write_mol2(directory, mol, coords) + self.write_frcmod(directory, mol, coords, j2cal) + + def write_mol2(self, directory, mol, coords): + with open(f"{directory}/{self.mol_name}_qforce{self.polar_title}.mol2", "w") as mol2: + self.write_mol2_title(mol2) + self.write_mol2_molecule(mol2, mol.topo, mol.terms) + self.write_mol2_atom(mol2, mol.topo, coords, mol.non_bonded) + self.write_mol2_bond(mol2, mol.topo, mol.terms) + + def write_frcmod(self, directory, mol, coords, j2cal): + with open(f"{directory}/{self.mol_name}_qforce{self.polar_title}.frcmod", "w") as frcmod: + self.write_frcmod_mass(frcmod, mol.non_bonded) + self.write_frcmod_bonds(frcmod, mol.terms, mol.non_bonded.alpha_map, j2cal) + self.write_frcmod_angles(frcmod, mol.terms, j2cal) + self.write_frcmod_dihedrals(frcmod, mol.terms, j2cal) + self.write_frcmod_nonbond(frcmod, mol.non_bonded) + + def write_mol2_title(self, mol2): + mol2.write(get_logo('#')) + mol2.write(f"# Name: {self.mol_name}\n") + + def write_mol2_molecule(self, mol2, topo, terms): + mol2.write(f"@MOLECULE\n{self.mol_name}\n") + n_bonds = 0 + for bond in terms['bond']: + n_bonds = n_bonds + 1 + mol2.write(f"{self.n_atoms:8d} {n_bonds:8d} {1:8d} {0:8d} {0:8d}\n") # n_bonds + mol2.write("esp") # type of charge + mol2.write("\n\n") + + def write_mol2_atom(self, mol2, topo, coords, non_bonded): + mol2.write("@ATOM\n") + for i_idx, (lj_type, a_name, q, mass) in enumerate(zip(non_bonded.lj_types, + self.atom_names, self.q, + self.masses), start=1): + + mol2.write(f"{i_idx:8d} {lj_type} {coords[i_idx-1][0]:10.4f}" + f"{coords[i_idx-1][1]:10.4f} {coords[i_idx-1][2]:10.4f}") + mol2.write(f" {i_idx} {1} {self.mol_name[0:3]} {q:10.6f}\n") + + def write_mol2_bond(self, mol2, topo, terms): + n_bonds = 1 + mol2.write("@BOND\n") + for bond in terms['bond']: + ids = bond.atomids + 1 + mol2.write(f'{n_bonds:>6} {ids[0]:>6}{ids[1]:>6} un \n') + n_bonds = n_bonds + 1 + mol2.write("@SUBSTRUCTURE\n") + mol2.write(f"{1:5} {self.mol_name[0:3]} {1:8} TEMP {0:>8d} **** **** {0:5} ROOT\n") + + def write_frcmod_mass(self, frcmod, non_bonded): + frcmod.write(f'{self.mol_name} - frcmod generated by QForce\n') + frcmod.write('MASS\n') + for i, mass in enumerate((self.masses), start=1): + frcmod.write(f"{i} {mass}\n") + + def write_frcmod_bonds(self, frcmod, terms, alpha_map, j2cal): + frcmod.write("\nBOND\n") + for bond in terms['bond']: + ids = bond.atomids + 1 + fconst = bond.fconst * j2cal/2 # kJ -> kcal + equ = bond.equ + frcmod.write(f"{ids[0]:<2}-{ids[1]:<2} {fconst:>10.2f} {equ:>10.2f}\n") + + for angle in terms['angle']: + ids = angle.atomids + 1 + if self.urey: + urey = [term for term in terms['urey'] if np.array_equal(term.atomids, + angle.atomids)] + if not self.urey or len(urey) == 0: + None + else: + urey_equ = urey[0].equ + urey_fconst = urey[0].fconst * j2cal/2 + frcmod.write(f"{ids[0]:<2}-{ids[2]:<2}{urey_fconst:>10.2f}{urey_equ:>10.2f}\n") + + + def write_frcmod_angles(self, frcmod, terms, j2cal): + frcmod.write("\nANGLE\n") + for angle in terms['angle']: + ids = angle.atomids + 1 + fconst = angle.fconst * j2cal/2 + equ = np.degrees(angle.equ) + + if self.urey: + urey = [term for term in terms['urey'] if np.array_equal(term.atomids, + angle.atomids)] + if not self.urey or len(urey) == 0: + frcmod.write(f"{ids[0]:<2}-{ids[1]:<2}-{ids[2]:<2}") + frcmod.write(f"{fconst:>10.2f}{equ:>10.2f} \n") + else: + frcmod.write(f"{ids[0]:<2}-{ids[1]:<2}-{ids[2]:<2}") + frcmod.write(f"{fconst:>10.2f}{equ:>10.2f} \n") + + + def write_frcmod_dihedrals(self, frcmod, terms, j2cal): + if len(terms['dihedral']) > 0: + frcmod.write("\nDIHE\n") + + # rigid dihedrals + if len(terms['dihedral/rigid']) > 0: + + for dihed in terms['dihedral/rigid']: + ids = dihed.atomids + 1 + equ = dihed.equ + + if dihed.equ == 0.00 or dihed.equ == 3.141592653589793: + equ = 180.0 + + fconst = dihed.fconst * j2cal + + fconstAmb = (2*fconst/(2**2)) + frcmod.write(f"{ids[0]:<2}-{ids[1]:<2}-{ids[2]:<2}-{ids[3]:<2}") + frcmod.write(f" 2 {fconstAmb:>15.2f} {equ:>15.2f} 2\n") + + if len(terms['dihedral/flexible']) > 0: + + # flexible dihedrals + for dihed in terms['dihedral/flexible']: + ids = dihed.atomids + 1 + c = dihed.equ * j2cal + fc = [1.0*c[0] + 0.5*c[2] + 0.3750*c[4], + 1.0*c[1] + 0.75*c[3] + 0.6250*c[5], + 0.5*c[2] + 0.5*c[4], + 0.25*c[3] + 0.3125*c[5], + 0.125*c[4], + 0.0625*c[5]] + + for n in range(3, 0, -1): + if n % 2 == 1: + equ = 180.0 + else: + equ = 0.00 + frcmod.write(f"{ids[0]:<2}-{ids[1]:<2}-{ids[2]:<2}-{ids[3]:<2}") + frcmod.write(f" 2 {fc[n]:>15.2f} {equ:>15.2f} -{n+1}\n") + + # Last term + frcmod.write(f"{ids[0]:<2}-{ids[1]:<2}-{ids[2]:<2}-{ids[3]:<2}") + frcmod.write(f" 2 {fc[0]:>15.2f} {0.00:>15.2f} 1\n") + + # inversion dihedrals + if len(terms['dihedral/inversion']) > 0: + + for dihed in terms['dihedral/inversion']: + ids = dihed.atomids + 1 + c0, c1, c2 = convert_to_inversion_rb(dihed.fconst, dihed.equ) + + fc = [(1.0*c0/2 + 0.5*c2/2)* j2cal, + 1.0*c1/2* j2cal, + 0.5*c2/2* j2cal] + + for n in range(3, 0, -1): + if n % 2 == 1: + equ = 180.0 + else: + equ = 0.00 + frcmod.write(f"{ids[0]:<2}-{ids[1]:<2}-{ids[2]:<2}-{ids[3]:<2}") + frcmod.write(f" 1 {fc[n]:>15.2f} {equ:>15.2f} -{n+1}\n") + + # Last term + frcmod.write(f"{ids[0]:<2}-{ids[1]:<2}-{ids[2]:<2}-{ids[3]:<2}") + frcmod.write(f" 1 {fc[0]:>15.2f} {0.00:>15.2f} 1\n") + + + # improper dihedrals + if len(terms['dihedral/improper']) > 0: + frcmod.write("\nIMPROPER\n") + for dihed in terms['dihedral/improper']: + ids = dihed.atomids + 1 + equ = np.degrees(dihed.equ) + fconst = dihed.fconst * j2cal/2 + frcmod.write(f"{ids[0]:<2}-{ids[1]:<2}-{ids[2]:<2}-{ids[3]:<2}") + frcmod.write(f" 1 {fconst:>15.2f} {equ:>15.2f} 2\n") + + def write_frcmod_nonbond(self, frcmod, non_bonded): + gro_atomtypes, gro_nonbonded, gro_1_4 = self.convert_to_gromacs_nonbonded(non_bonded) + atom_ids = {} # original atom types + + for i, (lj_type, a_name) in enumerate(zip(non_bonded.lj_types, self.atom_names), start=1): + atom_ids[i] = lj_type + + frcmod.write("\nNONB\n") + # support to gaff and gaff2 + if self.comb_rule == 2: + for i in range(1, self.n_atoms+1, 1): + frcmod.write(f"{i:<8}") + frcmod.write(f"{gro_atomtypes[atom_ids[i]][0]*5.612:>12.4f}") + frcmod.write(f"{gro_atomtypes[atom_ids[i]][1]/4.184:>12.4f}") + frcmod.write(f"\t{atom_ids[i]}\n") diff --git a/qforce/initialize.py b/qforce/initialize.py index 69acc89..e093cf7 100644 --- a/qforce/initialize.py +++ b/qforce/initialize.py @@ -9,12 +9,16 @@ from .qm.qm import QM, implemented_qm_software from .molecule.terms import Terms from .dihedral_scan import DihedralScan -from .misc import LOGO +from .misc import get_logo class Initialize(Colt): _user_input = """ [ff] +# MD software to use +# Currently amber supports only gaff2 atom types +output_software = gromacs :: str :: gromacs, amber + # Number of n equivalent neighbors needed to consider two atoms equivalent # Negative values turns off equivalence, 0 makes same elements equivalent n_equiv = 4 :: int @@ -164,7 +168,7 @@ def _check_and_copy_settings_file(job_dir, config_file): def initialize(filename, config_file, presets=None): - print(LOGO) + print(get_logo()) job_info = _get_job_info(filename) settings_file = _check_and_copy_settings_file(job_info.dir, config_file) diff --git a/qforce/main.py b/qforce/main.py index 90d163b..c494aac 100755 --- a/qforce/main.py +++ b/qforce/main.py @@ -9,7 +9,7 @@ from .frequencies import calc_qm_vs_md_frequencies from .hessian import fit_hessian -from .misc import check_if_file_exists, LOGO +from .misc import check_if_file_exists, get_logo from colt import from_commandline from colt.validator import Validator @@ -26,7 +26,7 @@ # File name for the optional options. options = :: file, optional, alias=o """, description={ - 'logo': LOGO, + 'logo': get_logo(), 'alias': 'qforce', 'arg_format': { 'name': 12, @@ -56,9 +56,9 @@ def run_qforce(input_arg, ext_q=None, ext_lj=None, config=None, presets=None): calc_qm_vs_md_frequencies(job, qm_hessian_out, md_hessian) ff = ForceField(job.name, config, mol, mol.topo.neighbors) - ff.write_gromacs(job.dir, mol, qm_hessian_out.coords) + ff.write_ff(job.dir, mol, qm_hessian_out.coords) - print_outcome(job.dir) + print_outcome(job.dir, config.ff.output_software) def run_hessian_fitting_for_external(job_dir, qm_data, ext_q=None, ext_lj=None, @@ -73,17 +73,17 @@ def run_hessian_fitting_for_external(job_dir, qm_data, ext_q=None, ext_lj=None, calc_qm_vs_md_frequencies(job, qm_hessian_out, md_hessian) ff = ForceField(job.name, config, mol, mol.topo.neighbors) - ff.write_gromacs(job.dir, mol, qm_hessian_out.coords) + ff.write_ff(job.dir, mol, qm_hessian_out.coords) - print_outcome(job.dir) + print_outcome(job.dir, config.ff.output_software) return mol.terms -def print_outcome(job_dir): +def print_outcome(job_dir, output_software): print(f'Output files can be found in the directory: {job_dir}.') - print('- Q-Force force field parameters in GROMACS format (gas.gro, gas.itp, gas.top).') + print(f'- Q-Force force field parameters in {output_software.upper()} format.') print('- QM vs MM vibrational frequencies, pre-dihedral fitting (frequencies.txt,' ' frequencies.pdf).') print('- Vibrational modes which can be visualized in VMD (frequencies.nmd).') - print('- QM vs MM dihedral profiles (if any) in "fragments" folder as ".pdf" files.') + print('- QM vs MM dihedral profiles (if any) in "fragments" folder as ".pdf" files.\n') diff --git a/qforce/misc.py b/qforce/misc.py index 80a649f..ccf373b 100755 --- a/qforce/misc.py +++ b/qforce/misc.py @@ -1,20 +1,7 @@ import os -import sys -LOGO = """ - ____ ______ - / __ \ | ____| - | | | |______| |__ ___ _ __ ___ ___ - | | | |______| __/ _ \| '__/ __/ _ \\ - | |__| | | | | (_) | | | (_| __/ - \___\_\ |_| \___/|_| \___\___| - - Selim Sami - University of Groningen - 2020 - ============================== -""" -LOGO_SEMICOL = """ +LOGO = """ ; ____ ______ ; / __ \ | ____| ; | | | |______| |__ ___ _ __ ___ ___ @@ -28,6 +15,10 @@ """ +def get_logo(line_start=' '): + return LOGO.replace(';', line_start) + + def check_if_file_exists(filename): if not os.path.exists(filename) and not os.path.exists(f'{filename}_qforce'): raise ValueError(f'"{filename}" does not exist.\n') diff --git a/setup.cfg b/setup.cfg index d3b2a37..8afd6e4 100644 --- a/setup.cfg +++ b/setup.cfg @@ -1,6 +1,6 @@ [metadata] name = qforce -version = 0.6.7 +version = 0.6.8 description = Q-Force: Quantum mechanically augmented molecular force fields long_description = file: README.md long_description_content_type = text/markdown