From ce7eb1b043954b95d9878fd449d818b1c95f1fea Mon Sep 17 00:00:00 2001 From: Mateus Zanotto <39527468+mateuszanotto@users.noreply.github.com> Date: Tue, 13 Sep 2022 13:05:08 -0300 Subject: [PATCH 01/28] Add files via upload --- forcefield.py | 717 ++++++++++++++++++++++++++++++++++++++++++++++++++ initialize.py | 178 +++++++++++++ main.py | 95 +++++++ 3 files changed, 990 insertions(+) create mode 100644 forcefield.py create mode 100644 initialize.py create mode 100644 main.py diff --git a/forcefield.py b/forcefield.py new file mode 100644 index 0000000..a383ac6 --- /dev/null +++ b/forcefield.py @@ -0,0 +1,717 @@ +import numpy as np +# +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 LOGO_HASH +import string + + +class ForceField(): + def __init__(self, job_name, config, mol, neighbors, exclude_all=[]): + self.polar = config.ff._polar + self.mol_name = job_name + self.n_atoms = mol.n_atoms + self.elements = mol.elements + self.q = self.set_charge(mol.non_bonded) + self.residue = config.ff.res_name[:5] + self.comb_rule = mol.non_bonded.comb_rule + self.fudge_lj = mol.non_bonded.fudge_lj + self.fudge_q = mol.non_bonded.fudge_q + self.urey = config.terms.urey + self.n_excl = config.ff.n_excl + self.atom_names = self.get_atom_names() + self.masses = [round(ATOMMASS[i], 5) for i in self.elements] + self.exclusions = self.make_exclusions(mol.non_bonded, neighbors, exclude_all) + self.pairs = self.make_pairs(neighbors, mol.non_bonded) + + if self.polar: + self.polar_title = '_polar' + else: + self.polar_title = '' + +## Amber + def write_amber(self, directory, mol, coords): + atom_ids, unique_at =self.get_atom_types(mol.topo, mol.non_bonded) + self.write_mol2(directory, mol, coords, atom_ids, unique_at) + self.write_frcmod(directory, mol, coords, atom_ids, unique_at) + + def write_mol2(self, directory, mol, coords, atom_ids, unique_at): + 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, atom_ids, unique_at) + self.write_mol2_bond(mol2, mol.topo, mol.terms, atom_ids, unique_at) + + def write_frcmod(self, directory, mol, coords, atom_ids, unique_at): + with open(f"{directory}/{self.mol_name}_qforce{self.polar_title}.frcmod", "w") as frcmod: + self.write_frcmod_mass(frcmod, mol.non_bonded, atom_ids, unique_at) + self.write_frcmod_bonds(frcmod, mol.terms, mol.non_bonded.alpha_map, atom_ids, unique_at) + self.write_frcmod_angles(frcmod, mol.terms, atom_ids, unique_at) + self.write_frcmod_dihedrals(frcmod, mol.terms, atom_ids, unique_at) + self.write_frcmod_nonbond(frcmod, mol.non_bonded, atom_ids, unique_at) + + +# TLeap is a helper program from AMBER to create the topologies. +# Would be good for the final user to create automatically a +# Tleap script with the new atom types. + + # def write_tleap_script(self, directory): + # with open(f"{directory}/tleap_script_{self.mol_name}.in", "w") as tleap: + # tleap.write(f""" + # # To run the script use tleap on amberTools and the following command: + # # $ tleap -f tleap_script_{self.mol_name}.in > tleap_script_{self.mol_name}.o + # + # source leaprc.protein.ff19SB\nsource leaprc.water.opc + # source leaprc.gaff\n# You can source other force fields here if necessary + # loadmol2 {self.mol_name}_qforce{self.polar_title}.mol2 + # + # loadamberparams {self.mol_name}_qforce{self.polar_title}.frcmod + # + # ## Add any other command, for example "solvateoct MOL TIP3PBOX 10.0" + # + # MOL = loadPDB 1lyd.solv.pdb + # saveamberparm MOL {self.mol_name}_qforce{self.polar_title}.prmtop + # saveamberparm MOL {self.mol_name}_qforce{self.polar_title}.inpcrd + # savepdb MOL {self.mol_name}_qforce{self.polar_title}.pdb """) + +############################################## +# MOL2 WRITING # +############################################## + + def write_mol2_title(self, mol2): + mol2.write(LOGO_HASH) + 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(f"esp")## type of charge + mol2.write("\n\n") + + def write_mol2_atom(self, mol2, topo, coords, non_bonded, atom_ids, unique_at): + mol2.write(f"@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} {coords[i_idx-1][1]:10.4f} {coords[i_idx-1][2]:10.4f}") + mol2.write(f" {unique_at[i_idx][atom_ids[i_idx]]} {1} {self.mol_name} {q:10.6f}\n") + + def write_mol2_bond(self, mol2, topo, terms, atom_ids, unique_at): + n_bonds = 1 + mol2.write(f"@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(f"@SUBSTRUCTURE\n") + mol2.write(f"{1:5} {self.mol_name} {1:8} TEMP {0:>8d} **** **** {0:5} ROOT\n") + +############################################## +# FRCMOD WRITING # +############################################## + + def write_frcmod_mass(self, frcmod, non_bonded, atom_ids, unique_at): + frcmod.write(f'{self.mol_name} - frcmod generated by QForce\n') + frcmod.write(f'MASS\n') + for i, mass in enumerate((self.masses), start=1): + frcmod.write(f"{unique_at[i][atom_ids[i]]} {mass}\n") + + def write_frcmod_bonds(self, frcmod, terms, alpha_map, atom_ids, unique_at): + frcmod.write("\nBOND\n") + for bond in terms['bond']: + ids = bond.atomids + 1 + fconst = bond.fconst * (0.239005)/2 # kJ/mol/A^2 -> kcal/mol/A^2 + equ = bond.equ # Ang + frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") + frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}") + frcmod.write(f"{fconst:>10.2f}{equ:>10.2f}\n") + + def write_frcmod_angles(self, frcmod, terms, atom_ids, unique_at): + frcmod.write("\nANGLE\n") + for angle in terms['angle']: + ids = angle.atomids + 1 + fconst = angle.fconst * (0.239005)/2 # kJ/mol/rad^2 -> kcal/mol/rad^2 + equ = np.degrees(angle.equ) # Degrees -> Degrees + + 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"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") + frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") + frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}") + frcmod.write(f"{fconst:>10.2f}{equ:>10.2f}\n") + else: + urey_equ = urey[0].equ + urey_fconst = urey[0].fconst * (0.239005)/2 # kJ/mol/rad^2 -> kcal/mol/rad^2 + frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") + frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") + frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}") + frcmod.write(f"{fconst:>10.2f}{equ:>10.2f}\n") + + def write_frcmod_dihedrals(self, frcmod, terms, atom_ids, unique_at): + 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 * (0.239005) # kJ/mol -> kcal/mol #scale cos + + fconstAmb = (2*fconst/(2**2)) # Vn = (2*fcte/(n**2)) + frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") #Atom# 1 + frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") #Atom 2 + frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}-") #Atom 3 + frcmod.write(f"{unique_at[ids[3]][atom_ids[ids[3]]]:<2}") #Atom 4 + frcmod.write(f" 1 {fconstAmb:>15.2f} {equ:>15.2f} 2\n") #IDIVF=1, fcte, angle, n=2 (n*angle) + + if len(terms['dihedral/flexible']) > 0: + + for dihed in terms['dihedral/flexible']: + ids = dihed.atomids + 1 + c = dihed.equ + + for n in range(3, 0, -1): + if n%2==1: + equ = 180.0 + else: + equ = 0.00 + frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") #Atom 1 + frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") #Atom 2 + frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}-") #Atom 3 + frcmod.write(f"{unique_at[ids[3]][atom_ids[ids[3]]]:<2}") #Atom 4 + frcmod.write(f" 1 {c[n]:>15.2f} {equ:>15.2f} -{n+1}\n") #IDIVF=1, fcte, angle, n=1 (n*angle) + + # Last term without -n + frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") #Atom 1 + frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") #Atom 2 + frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}-") #Atom 3 + frcmod.write(f"{unique_at[ids[3]][atom_ids[ids[3]]]:<2}") #Atom 4 + frcmod.write(f" 1 {c[0]:>15.2f} {0.00:>15.2f} 1\n") #IDIVF=1, fcte, angle, n=1 (n*angle) + + + # 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) # Degrees -> Degrees + fconst = dihed.fconst * (0.239005)/2 # kJ/mol -> kcal/mol -> f/2 as is Vn/2 in amber + frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") #Atom 1 + frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") #Atom 2 + frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}-") #Atom 3 + frcmod.write(f"{unique_at[ids[3]][atom_ids[ids[3]]]:<2}") #Atom 4 + frcmod.write(f" 1 {fconst:>15.2f} {equ:>15.2f} 2\n") #IDIVF=1, fcte, angle, n=2 (n*angle) + + def write_frcmod_nonbond(self, frcmod, non_bonded, atom_ids, unique_at): + gro_atomtypes, gro_nonbonded, gro_1_4 = self.convert_to_gromacs_nonbonded(non_bonded) + 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"{unique_at[i][atom_ids[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") + + +############################################## +# GET ATOM TYPES AND IDS # +############################################## + def get_atom_types(self, topo, non_bonded): + atom_ids={} #dictonary containing original atom types + unique_at={} #dictonary containing new unique atom types + unique_masses={} + + ascii_lowercase = list(string.ascii_lowercase) + ascii_digits = list(string.digits) + ascii_uppercase = list(string.ascii_uppercase) + + for i, (lj_type, a_name) in enumerate(zip(non_bonded.lj_types, self.atom_names), start=1): + atom_ids[i]=lj_type + if self.n_atoms < 36: + if i <= 10: + unique_at[i]={atom_ids[i] : "Q{}".format(ascii_digits[i-1])} + elif i > 10: + unique_at[i]={atom_ids[i] : "Q{}".format(ascii_lowercase[i-11])} + else: + unique_at[i]={atom_ids[i] : "{}{}".format(ascii_uppercase[i-1], ascii_lowercase[i-1])} + + ## DEBUG : check correspondance between atom_ids and unique_at + #for i, lj_type in enumerate(zip(non_bonded.lj_types), start=1): + # print(unique_at) + # for i in range(1,self.n_atoms+1,1): + # print (i, atom_ids[i], unique_at[i][atom_ids[i]]) + # print ("------------------") + + return atom_ids, unique_at + +############################################## + + def write_gromacs(self, directory, mol, coords): + self.write_itp(mol, directory) + self.write_top(directory) + self.write_gro(directory, coords, mol.non_bonded.alpha_map) + + def write_top(self, directory): + with open(f"{directory}/gas{self.polar_title}.top", "w") as top: + # defaults + top.write("\n[ defaults ]\n") + top.write(";nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ\n") + top.write(f"{1:>7}{self.comb_rule:>12}{'yes':>12}{self.fudge_lj:>10}{self.fudge_q:>10}" + "\n\n\n") + + top.write("; Include the molecule ITP\n") + top.write(f'#include "./{self.mol_name}_qforce{self.polar_title}.itp"\n\n\n') + + size = len(self.mol_name) + top.write("[ system ]\n") + top.write(f"; {' '*(size-6)}name\n") + top.write(f"{' '*(6-size)}{self.mol_name}\n\n\n") + + top.write("[ molecules ]\n") + top.write(f"; {' '*(size-10)}compound n_mol\n") + top.write(f"{' '*(10-size)}{self.mol_name} 1\n") + + def write_gro(self, directory, coords, alpha_map, box=[20., 20., 20.]): + n_atoms = self.n_atoms + if self.polar: + n_atoms += len(alpha_map.keys()) + coords_nm = coords*0.1 + with open(f"{directory}/gas{self.polar_title}.gro", "w") as gro: + gro.write(f"{self.mol_name}\n") + gro.write(f"{n_atoms:>6}\n") + for i, (a_name, coord) in enumerate(zip(self.atom_names, coords_nm), start=1): + gro.write(f"{1:>5}{self.residue:<5}") + gro.write(f"{a_name:>5}{i:>5}") + gro.write(f"{coord[0]:>8.3f}{coord[1]:>8.3f}{coord[2]:>8.3f}\n") + if self.polar: + for i, (atom, drude) in enumerate(alpha_map.items(), start=1): + gro.write(f"{2:>5}{self.residue:<5}{f'D{i}':>5}{drude+1:>5}") + gro.write(f"{coords_nm[atom][0]:>8.3f}{coords_nm[atom][1]:>8.3f}") + gro.write(f"{coords_nm[atom][2]:>8.3f}\n") + gro.write(f'{box[0]:>12.5f}{box[1]:>12.5f}{box[2]:>12.5f}\n') + + 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) + self.write_itp_atoms_and_molecule(itp, mol.non_bonded) + if self.polar: + self.write_itp_polarization(itp, mol.non_bonded) + self.write_itp_bonds(itp, mol.terms, mol.non_bonded.alpha_map) + self.write_itp_angles(itp, mol.terms) + self.write_itp_dihedrals(itp, mol.terms) + self.write_itp_pairs(itp) + self.write_itp_exclusions(itp) + itp.write('\n') + + def convert_to_gromacs_nonbonded(self, non_bonded): + a_types, nb_pairs, nb_1_4 = {}, {}, {} + + for pair, val in non_bonded.lj_pairs.items(): + if non_bonded.comb_rule != 1: + if val[0] == 0: + a, b = 0, 0 + else: + a, b = calc_sigma_epsilon(val[0], val[1]) + a *= 0.1 + else: + a = val[0] * 1e-6 + b = val[1] * 1e-12 + + if pair[0] == pair[1]: + a_types[pair[0]] = [a, b] + + nb_pairs[pair] = [a, b] + + for pair, val in non_bonded.lj_1_4.items(): + if non_bonded.comb_rule != 1: + a, b = calc_sigma_epsilon(val[0], val[1]) + a *= 0.1 + else: + a = val[0] * 1e-6 + b = val[1] * 1e-12 + + nb_1_4[pair] = [a, b] + + return a_types, nb_pairs, nb_1_4 + + def write_itp_atoms_and_molecule(self, itp, non_bonded): + gro_atomtypes, gro_nonbonded, gro_1_4 = self.convert_to_gromacs_nonbonded(non_bonded) + + # atom types + itp.write("\n[ atomtypes ]\n") + if self.comb_rule == 1: + itp.write("; name at_num mass charge type c6 c12\n") + else: + 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_params[0]:>12.5e} {lj_params[1]:>12.5e}\n') + + if self.polar: + itp.write(f'{"DP":>8} {0:>8.4f} {0:>8.4f} {"S":>2} {0:>12.5e} {0:>12.5e}\n') + + # non-bonded pair types + itp.write("\n[ nonbond_params ]\n") + if self.comb_rule == 1: + itp.write("; name1 name2 type c6 c12\n") + else: + itp.write("; name1 name2 type sigma epsilon\n") + + for lj_types, lj_params in gro_nonbonded.items(): + itp.write(f'{lj_types[0]:>8} {lj_types[1]:>8} 1') + itp.write(f'{lj_params[0]:>15.5e} {lj_params[1]:>15.5e}\n') + + # Write 1-4 pair types + if self.n_excl == 2 and gro_1_4 != {}: + itp.write("\n[ pairtypes ]\n") + if self.comb_rule == 1: + itp.write("; name1 name2 type c6 c12\n") + else: + itp.write("; name1 name2 type sigma epsilon\n") + for lj_types, lj_params in gro_1_4.items(): + itp.write(f'{lj_types[0]:>8} {lj_types[1]:>8} 1') + itp.write(f'{lj_params[0]:>15.5e} {lj_params[1]:>15.5e}\n') + + # molecule type + space = " "*(len(self.mol_name)-5) + itp.write("\n[ moleculetype ]\n") + itp.write(f";{space}name nrexcl\n") + itp.write(f"{self.mol_name}{3:>7}\n") + + # atoms + itp.write("\n[ atoms ]\n") + itp.write("; nr type resnr resnm atom cgrp charge mass\n") + for i, (lj_type, a_name, q, mass) in enumerate(zip(non_bonded.lj_types, self.atom_names, + self.q, self.masses), start=1): + itp.write(f'{i:>5}{lj_type:>9}{1:>6}{self.residue:>6}{a_name:>7}{i:>5}{q:>11.5f}') + itp.write(f'{mass:>10.5f}\n') + + if self.polar: + for i, (atom, drude) in enumerate(non_bonded.alpha_map.items(), start=1): + itp.write(f'{drude+1:>5}{"DP":>9}{2:>6}{"DRU":>6}{f"D{i}":>7}{atom+1:>5}') + itp.write(f'{-8.:>11.5f}{0.:>10.5f}\n') + + def write_itp_polarization(self, itp, non_bonded): + # polarization + itp.write("\n[ polarization ]\n") + itp.write("; i j f alpha\n") + for atom, drude in non_bonded.alpha_map.items(): + alpha = non_bonded.alpha[atom]*1e-3 + itp.write(f"{atom+1:>6}{drude+1:>6}{1:>6}{alpha:>14.8f}\n") + + # # thole polarization + # if self.thole != []: + # itp.write("\n[ thole_polarization ]\n") + # itp.write("; ai di aj dj f a alpha(i) " + # "alpha(j)\n") + # for tho in self.thole: + # itp.write("{:>6}{:>6}{:>6}{:>6}{:>4}{:>7.2f}{:>14.8f}{:>14.8f}\n".format(*tho)) + + def write_itp_pairs(self, itp): + if self.pairs != []: + itp.write("\n[ pairs ]\n") + itp.write("; ai aj func \n") + for pair in self.pairs: + itp.write(f"{pair[0]+1:>6}{pair[1]+1:>6}{1:>6}\n") + + def write_itp_bonds(self, itp, terms, alpha_map): + itp.write("\n[ bonds ]\n") + itp.write("; ai aj f r0 kb\n") + for bond in terms['bond']: + ids = bond.atomids + 1 + equ = bond.equ * 0.1 + fconst = bond.fconst * 100 + itp.write(f'{ids[0]:>6}{ids[1]:>6}{1:>6}{equ:>10.5f}{fconst:>10.0f}\n') + + if self.polar: + itp.write('; ai aj f - polar connections\n') + for bond in terms['bond']: + a1, a2 = bond.atomids + + if a2 in alpha_map.keys(): + itp.write(f'{a1+1:>6}{alpha_map[a2]+1:>6}{5:>6}\n') + if a1 in alpha_map.keys(): + itp.write(f'{a2+1:>6}{alpha_map[a1]+1:>6}{5:>6}\n') + + def write_itp_angles(self, itp, terms): + itp.write("\n[ angles ]\n") + itp.write("; ai aj ak f th0 kth\n") + for angle in terms['angle']: + ids = angle.atomids + 1 + equ = np.degrees(angle.equ) + fconst = angle.fconst + + 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: + itp.write(f'{ids[0]:>6}{ids[1]:>6}{ids[2]:>6}{1:>6}{equ:>11.3f}{fconst:>13.3f}\n') + else: + urey_equ = urey[0].equ * 0.1 + urey_fconst = urey[0].fconst * 100 + itp.write(f'{ids[0]:>6}{ids[1]:>6}{ids[2]:>6}{5:>6}{equ:>11.3f}{fconst:>13.3f}' + f'{urey_equ:>10.5f}{urey_fconst:>13.3f}\n') + + def write_itp_dihedrals(self, itp, terms): + if len(terms['dihedral']) > 0: + itp.write("\n[ dihedrals ]\n") + + # rigid dihedrals + if len(terms['dihedral/rigid']) > 0: + itp.write("; rigid dihedrals \n") + itp.write("; ai aj ak al f th0 kth\n") + + for dihed in terms['dihedral/rigid']: + ids = dihed.atomids + 1 + equ = np.degrees(dihed.equ) + fconst = dihed.fconst + + itp.write(f'{ids[0]:>6}{ids[1]:>6}{ids[2]:>6}{ids[3]:>6}{2:>6}{equ:>11.3f}') + itp.write(f'{fconst:>13.3f}\n') + + # improper dihedrals + if len(terms['dihedral/improper']) > 0: + itp.write("; improper dihedrals \n") + itp.write("; ai aj ak al f th0 kth\n") + + for dihed in terms['dihedral/improper']: + ids = dihed.atomids + 1 + equ = np.degrees(dihed.equ) + fconst = dihed.fconst + + itp.write(f'{ids[0]:>6}{ids[1]:>6}{ids[2]:>6}{ids[3]:>6}{2:>6}{equ:>11.3f}') + itp.write(f'{fconst:>13.3f}\n') + + # flexible dihedrals + if len(terms['dihedral/flexible']) > 0: + itp.write("; flexible dihedrals \n") + itp.write('; ai aj ak al f c0 c1 c2 c3') + itp.write(' c4 c5\n') + + for dihed in terms['dihedral/flexible']: + ids = dihed.atomids + 1 + c = dihed.equ + + itp.write(f'{ids[0]:>6}{ids[1]:>6}{ids[2]:>6}{ids[3]:>6}{3:>6}{c[0]:>11.3f}') + itp.write(f'{c[1]:>11.3f}{c[2]:>11.3f}{c[3]:>11.3f}{c[4]:>11.3f}{c[5]:>11.3f}\n') + + # inversion dihedrals + if len(terms['dihedral/inversion']) > 0: + itp.write("; inversion dihedrals \n") + itp.write('; ai aj ak al f c0 c1 c2\n') + + for dihed in terms['dihedral/inversion']: + ids = dihed.atomids + 1 + # itp.write(f'{ids[0]:>6}{ids[1]:>6}{ids[2]:>6}{ids[3]:>6}{1:>6}' + # f'{0.0:>11.3f}{dihed.fconst:>11.3f}{3:>11}\n') + c0, c1, c2 = convert_to_inversion_rb(dihed.fconst, dihed.equ) + itp.write(f'{ids[0]:>6}{ids[1]:>6}{ids[2]:>6}{ids[3]:>6}{3:>6}' + f'{c0:>11.3f}{c1:>11.3f}{c2:>11.3f}{0:>11.1f}{0:>11.1f}{0:>11.1f}\n') + + def write_itp_exclusions(self, itp): + if any(len(exclusion) > 0 for exclusion in self.exclusions): + itp.write("\n[ exclusions ]\n") + for i, exclusion in enumerate(self.exclusions): + if len(exclusion) > 0: + exclusion = sorted(set(exclusion)) + itp.write("{} ".format(i+1)) + itp.write(("{} "*len(exclusion)).format(*exclusion)) + itp.write("\n") + + def make_pairs(self, neighbors, non_bonded): + pairs, polar_pairs = [], [] + + for pair in non_bonded.pairs: + pairs.append(sorted(pair)) + + if self.n_excl == 2: + for i in range(self.n_atoms): + for neigh in neighbors[2][i]: + if (i < neigh and [i, neigh] not in pairs and (i, neigh) + not in non_bonded.exclusions): + pairs.append([i, neigh]) + + if self.polar: + for a1, a2 in pairs: + if a2 in non_bonded.alpha_map.keys(): + polar_pairs.append([a1, non_bonded.alpha_map[a2]]) + if a1 in non_bonded.alpha_map.keys(): + polar_pairs.append([a2, non_bonded.alpha_map[a1]]) + if a1 in non_bonded.alpha_map.keys() and a2 in non_bonded.alpha_map.keys(): + polar_pairs.append([non_bonded.alpha_map[a1], non_bonded.alpha_map[a2]]) + + return pairs+polar_pairs + + def make_exclusions(self, non_bonded, neighbors, exclude_all): + exclusions = [[] for _ in range(self.n_atoms)] + + # input exclusions + for exclusion in non_bonded.exclusions: + exclusions[exclusion[0]].append(exclusion[1]+1) + + # exclusions for input pairs + for pair in non_bonded.pairs: + exclusions[pair[0]].append(pair[1]+1) + + # fragment capping atom exclusions + for i in exclude_all: + exclusions[i].extend(np.arange(1, self.n_atoms+1)) + + if self.polar: + exclusions = self.polarize_exclusions(non_bonded.alpha_map, non_bonded.exclusions, + neighbors, exclude_all, exclusions) + + return exclusions + + def polarize_exclusions(self, alpha_map, input_exclusions, neighbors, exclude_all, exclusions): + n_polar_atoms = len(alpha_map.keys()) + exclusions.extend([[] for _ in range(n_polar_atoms)]) + + # input exclusions + for exclu in input_exclusions: + if exclu[0] in alpha_map.keys(): + exclusions[alpha_map[exclu[0]]].append(exclu[1]+1) + if exclu[1] in alpha_map.keys(): + exclusions[alpha_map[exclu[1]]].append(exclu[0]+1) + if exclu[0] in alpha_map.keys() and exclu[1] in alpha_map.keys(): + exclusions[alpha_map[exclu[0]]].append(alpha_map[exclu[1]]+1) + + # fragment capping atom exclusions + for i in exclude_all: + exclusions[i].extend(np.arange(self.n_atoms+1, self.n_atoms+n_polar_atoms+1)) + if i in alpha_map.keys(): + exclusions[alpha_map[i]].extend(np.arange(1, self.n_atoms+n_polar_atoms+1)) + + return exclusions + + def get_atom_names(self): + atom_names = [] + atom_dict = {} + + for i, elem in enumerate(self.elements): + sym = ATOM_SYM[elem] + if sym not in atom_dict.keys(): + atom_dict[sym] = 1 + else: + atom_dict[sym] += 1 + atom_names.append(f'{sym}{atom_dict[sym]}') + return atom_names + + def add_restraints(self, restraints, directory, fc=1000): + with open(f"{directory}/{self.mol_name}_qforce{self.polar_title}.itp", "a") as itp: + itp.write("[ dihedral_restraints ]\n") + itp.write("; ai aj ak al type phi dp kfac\n") + for restraint in restraints: + a1, a2, a3, a4 = restraint[0]+1 + phi = np.degrees(restraint[1]) + itp.write(f'{a1:>5}{a2:>5}{a3:>5}{a4:>5}{1:>5} {phi:>10.4f} 0.0 {fc}\n') + + def set_charge(self, non_bonded): + q = np.copy(non_bonded.q) + if self.polar: + 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])])) diff --git a/initialize.py b/initialize.py new file mode 100644 index 0000000..5bea8c7 --- /dev/null +++ b/initialize.py @@ -0,0 +1,178 @@ +import os +import shutil +from io import StringIO +from types import SimpleNamespace +import pkg_resources +# +from colt import Colt +# +from .qm.qm import QM, implemented_qm_software +from .molecule.terms import Terms +from .dihedral_scan import DihedralScan +from .misc import LOGO + + +class Initialize(Colt): + _user_input = """ +[ff] +# MD software to use +# Currently amber supports only gaff2 atom types +forcefield = 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 + +# Number of first n neighbors to exclude in the forcefield +n_excl = 2 :: int :: [2, 3] + +# Lennard jones method for the forcefield +lennard_jones = opls_auto :: str :: [gromos_auto, gromos, opls_auto, opls, gaff, gaff2, ext] + +# Use externally provided point charges in the file "ext_q" in the job directyory +ext_charges = no :: bool + +# Scale QM charges to account for condensed phase polarization (should be set to 1 for gas phase) +charge_scaling = 1.2 :: float + +# If user chooses ext_charges=True, by default fragments will still use the chosen QM method for +# determining fragment charges. This is to avoid spuriously high charges on capping hydrogens. +# However, using QM charges for fragments and ext_charges for the molecule can also result in +# inaccuracies if these two charges are very different. +use_ext_charges_for_frags = no :: bool + +# Additional exclusions (GROMACS format) +exclusions = :: literal + +# Switch standard non-bonded interactions between two atoms to pair interactions +# (provide atom pairs in each row) +pairs = :: literal + +# Path for the external FF library for Lennard-Jones parameters (GROMACS format). +ext_lj_lib = :: folder, optional + +# Lennard-Jones fudge parameter for 1-4 interactions for when "lennard_jones" is set to "ext". +ext_lj_fudge = :: float, optional + +# Coulomb fudge parameter for 1-4 interactions for when "lennard_jones" is set to "ext". +ext_q_fudge = :: float, optional + +# Lennard-Jones combinations rules for when "lennard_jones" is set to "ext" (GROMACS numbering). +ext_comb_rule = :: int, optional :: [1, 2, 3] + +# Name of the atom type for capping hydrogens for when "lennard_jones" is set to "ext" +ext_h_cap = :: str, optional + +# Set all dihedrals as rigid (no dihedral scans) +all_rigid = no :: bool + +# Use D4 method +_d4 = no :: bool + +# Residue name printed on the force field file (Max 5 characters) +res_name = MOL :: str + +# Polarize a coordinate file and quit (requires itp_file) +_polarize = no :: bool + +# Name of itp file (only needed for polarize option) +_itp_file = itp_file_missing :: str + +# Make a polarizable FF +_polar = no :: bool + +# Scale the C6 dispersion interactions in the polarizable version of the FF +_polar_c6_scale = 0.8 :: float + +# Specifically not scale some of the atoms +_polar_not_scale_c6 = :: literal + +# Manual polarizabilities in the file ext_alpha +_ext_alpha = no :: bool + +""" + + @staticmethod + def _set_config(config): + config['qm'].update(config['qm']['software']) + config['qm'].update({'software': config['qm']['software'].value}) + config.update({key: SimpleNamespace(**val) for key, val in config.items()}) + return SimpleNamespace(**config) + + @classmethod + def _extend_user_input(cls, questions): + questions.generate_block("qm", QM.colt_user_input) + questions.generate_block("scan", DihedralScan.colt_user_input) + questions.generate_cases("software", {key: software.colt_user_input for key, software in + implemented_qm_software.items()}, block='qm') + questions.generate_block("terms", Terms.get_questions()) + + @classmethod + def from_config(cls, config): + return cls._set_config(config) + + @staticmethod + def set_basis(value): + if value.endswith('**'): + return f'{value[:-2]}(D,P)'.upper() + if value.endswith('*'): + return f'{value[:-1]}(D)'.upper() + return value.upper() + + @staticmethod + def set_dispersion(value): + if value.lower() in ["no", "false", "n", "f"]: + return False + return value.upper() + + +def _get_job_info(filename): + job = {} + filename = filename.rstrip('/') + base = os.path.basename(filename) + path = os.path.dirname(filename) + if path != '': + path = f'{path}/' + + if os.path.isfile(filename): + job['coord_file'] = filename + job['name'] = base.split('.')[0] + else: + job['coord_file'] = False + job['name'] = base.split('_qforce')[0] + + job['dir'] = f'{path}{job["name"]}_qforce' + job['frag_dir'] = f'{job["dir"]}/fragments' + job['md_data'] = pkg_resources.resource_filename('qforce', 'data') + os.makedirs(job['dir'], exist_ok=True) + return SimpleNamespace(**job) + + +def _check_and_copy_settings_file(job_dir, config_file): + """ + If options are provided as a file, copy that to job directory. + If options are provided as StringIO, write that to job directory. + """ + + settings_file = os.path.join(job_dir, 'settings.ini') + + if config_file is not None: + if isinstance(config_file, StringIO): + with open(settings_file, 'w') as fh: + config_file.seek(0) + fh.write(config_file.read()) + else: + shutil.copy2(config_file, settings_file) + + return settings_file + + +def initialize(filename, config_file, presets=None): + #print(LOGO) + + job_info = _get_job_info(filename) + settings_file = _check_and_copy_settings_file(job_info.dir, config_file) + + config = Initialize.from_questions(config=settings_file, presets=presets, check_only=True) + + return config, job_info diff --git a/main.py b/main.py new file mode 100644 index 0000000..b14fc33 --- /dev/null +++ b/main.py @@ -0,0 +1,95 @@ +from .polarize import polarize +from .initialize import initialize +from .qm.qm import QM +from .qm.qm_base import HessianOutput +from .forcefield import ForceField +from .molecule import Molecule +from .fragment import fragment +from .dihedral_scan import DihedralScan +from .frequencies import calc_qm_vs_md_frequencies +from .hessian import fit_hessian + +from .misc import check_if_file_exists, LOGO +from colt import from_commandline +from colt.validator import Validator + + +# define new validator +Validator.overwrite_validator("file", check_if_file_exists) + + +@from_commandline(""" +# Input coordinate file mol.ext (ext: pdb, xyz, gro, ...) +# or directory (mol or mol_qforce) name. +file = :: file + +# File name for the optional options. +options = :: file, optional, alias=o +""", description={ + 'logo': LOGO, + 'alias': 'qforce', + 'arg_format': { + 'name': 12, + 'comment': 60, + }, +}) +def run(file, options): + run_qforce(input_arg=file, config=options) + + +def run_qforce(input_arg, ext_q=None, ext_lj=None, config=None, presets=None): + config, job = initialize(input_arg, config, presets) + + if config.ff._polarize: + polarize(job, config.ff) + + qm = QM(job, config.qm) + qm_hessian_out = qm.read_hessian() + + mol = Molecule(config, job, qm_hessian_out, ext_q, ext_lj) + + md_hessian = fit_hessian(config.terms, mol, qm_hessian_out) + + if len(mol.terms['dihedral/flexible']) > 0 and config.scan.do_scan: + fragments = fragment(mol, qm, job, config) + DihedralScan(fragments, mol, job, config) + + calc_qm_vs_md_frequencies(job, qm_hessian_out, md_hessian) + ff = ForceField(job.name, config, mol, mol.topo.neighbors) + + forcefield = config.ff.forcefield + if forcefield in ["gromacs"]: + ff.write_gromacs(job.dir, mol, qm_hessian_out.coords) + + if forcefield in ["amber"]: + ff.write_amber(job.dir, mol, qm_hessian_out.coords) + + print_outcome(job.dir) + + +def run_hessian_fitting_for_external(job_dir, qm_data, ext_q=None, ext_lj=None, + config=None, presets=None): + config, job = initialize(job_dir, config, presets) + + qm_hessian_out = HessianOutput(config.qm.vib_scaling, **qm_data) + + mol = Molecule(config, job, qm_hessian_out, ext_q, ext_lj) + + md_hessian = fit_hessian(config.terms, mol, qm_hessian_out) + 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) + + print_outcome(job.dir) + + return mol.terms + + +def print_outcome(job_dir): + 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('- 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.') From c4a2d148d4a2cd539d61b5a3dac84fdfa70fc18e Mon Sep 17 00:00:00 2001 From: Mateus Zanotto <39527468+mateuszanotto@users.noreply.github.com> Date: Tue, 13 Sep 2022 13:10:09 -0300 Subject: [PATCH 02/28] Delete forcefield.py --- forcefield.py | 717 -------------------------------------------------- 1 file changed, 717 deletions(-) delete mode 100644 forcefield.py diff --git a/forcefield.py b/forcefield.py deleted file mode 100644 index a383ac6..0000000 --- a/forcefield.py +++ /dev/null @@ -1,717 +0,0 @@ -import numpy as np -# -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 LOGO_HASH -import string - - -class ForceField(): - def __init__(self, job_name, config, mol, neighbors, exclude_all=[]): - self.polar = config.ff._polar - self.mol_name = job_name - self.n_atoms = mol.n_atoms - self.elements = mol.elements - self.q = self.set_charge(mol.non_bonded) - self.residue = config.ff.res_name[:5] - self.comb_rule = mol.non_bonded.comb_rule - self.fudge_lj = mol.non_bonded.fudge_lj - self.fudge_q = mol.non_bonded.fudge_q - self.urey = config.terms.urey - self.n_excl = config.ff.n_excl - self.atom_names = self.get_atom_names() - self.masses = [round(ATOMMASS[i], 5) for i in self.elements] - self.exclusions = self.make_exclusions(mol.non_bonded, neighbors, exclude_all) - self.pairs = self.make_pairs(neighbors, mol.non_bonded) - - if self.polar: - self.polar_title = '_polar' - else: - self.polar_title = '' - -## Amber - def write_amber(self, directory, mol, coords): - atom_ids, unique_at =self.get_atom_types(mol.topo, mol.non_bonded) - self.write_mol2(directory, mol, coords, atom_ids, unique_at) - self.write_frcmod(directory, mol, coords, atom_ids, unique_at) - - def write_mol2(self, directory, mol, coords, atom_ids, unique_at): - 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, atom_ids, unique_at) - self.write_mol2_bond(mol2, mol.topo, mol.terms, atom_ids, unique_at) - - def write_frcmod(self, directory, mol, coords, atom_ids, unique_at): - with open(f"{directory}/{self.mol_name}_qforce{self.polar_title}.frcmod", "w") as frcmod: - self.write_frcmod_mass(frcmod, mol.non_bonded, atom_ids, unique_at) - self.write_frcmod_bonds(frcmod, mol.terms, mol.non_bonded.alpha_map, atom_ids, unique_at) - self.write_frcmod_angles(frcmod, mol.terms, atom_ids, unique_at) - self.write_frcmod_dihedrals(frcmod, mol.terms, atom_ids, unique_at) - self.write_frcmod_nonbond(frcmod, mol.non_bonded, atom_ids, unique_at) - - -# TLeap is a helper program from AMBER to create the topologies. -# Would be good for the final user to create automatically a -# Tleap script with the new atom types. - - # def write_tleap_script(self, directory): - # with open(f"{directory}/tleap_script_{self.mol_name}.in", "w") as tleap: - # tleap.write(f""" - # # To run the script use tleap on amberTools and the following command: - # # $ tleap -f tleap_script_{self.mol_name}.in > tleap_script_{self.mol_name}.o - # - # source leaprc.protein.ff19SB\nsource leaprc.water.opc - # source leaprc.gaff\n# You can source other force fields here if necessary - # loadmol2 {self.mol_name}_qforce{self.polar_title}.mol2 - # - # loadamberparams {self.mol_name}_qforce{self.polar_title}.frcmod - # - # ## Add any other command, for example "solvateoct MOL TIP3PBOX 10.0" - # - # MOL = loadPDB 1lyd.solv.pdb - # saveamberparm MOL {self.mol_name}_qforce{self.polar_title}.prmtop - # saveamberparm MOL {self.mol_name}_qforce{self.polar_title}.inpcrd - # savepdb MOL {self.mol_name}_qforce{self.polar_title}.pdb """) - -############################################## -# MOL2 WRITING # -############################################## - - def write_mol2_title(self, mol2): - mol2.write(LOGO_HASH) - 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(f"esp")## type of charge - mol2.write("\n\n") - - def write_mol2_atom(self, mol2, topo, coords, non_bonded, atom_ids, unique_at): - mol2.write(f"@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} {coords[i_idx-1][1]:10.4f} {coords[i_idx-1][2]:10.4f}") - mol2.write(f" {unique_at[i_idx][atom_ids[i_idx]]} {1} {self.mol_name} {q:10.6f}\n") - - def write_mol2_bond(self, mol2, topo, terms, atom_ids, unique_at): - n_bonds = 1 - mol2.write(f"@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(f"@SUBSTRUCTURE\n") - mol2.write(f"{1:5} {self.mol_name} {1:8} TEMP {0:>8d} **** **** {0:5} ROOT\n") - -############################################## -# FRCMOD WRITING # -############################################## - - def write_frcmod_mass(self, frcmod, non_bonded, atom_ids, unique_at): - frcmod.write(f'{self.mol_name} - frcmod generated by QForce\n') - frcmod.write(f'MASS\n') - for i, mass in enumerate((self.masses), start=1): - frcmod.write(f"{unique_at[i][atom_ids[i]]} {mass}\n") - - def write_frcmod_bonds(self, frcmod, terms, alpha_map, atom_ids, unique_at): - frcmod.write("\nBOND\n") - for bond in terms['bond']: - ids = bond.atomids + 1 - fconst = bond.fconst * (0.239005)/2 # kJ/mol/A^2 -> kcal/mol/A^2 - equ = bond.equ # Ang - frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") - frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}") - frcmod.write(f"{fconst:>10.2f}{equ:>10.2f}\n") - - def write_frcmod_angles(self, frcmod, terms, atom_ids, unique_at): - frcmod.write("\nANGLE\n") - for angle in terms['angle']: - ids = angle.atomids + 1 - fconst = angle.fconst * (0.239005)/2 # kJ/mol/rad^2 -> kcal/mol/rad^2 - equ = np.degrees(angle.equ) # Degrees -> Degrees - - 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"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") - frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") - frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}") - frcmod.write(f"{fconst:>10.2f}{equ:>10.2f}\n") - else: - urey_equ = urey[0].equ - urey_fconst = urey[0].fconst * (0.239005)/2 # kJ/mol/rad^2 -> kcal/mol/rad^2 - frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") - frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") - frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}") - frcmod.write(f"{fconst:>10.2f}{equ:>10.2f}\n") - - def write_frcmod_dihedrals(self, frcmod, terms, atom_ids, unique_at): - 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 * (0.239005) # kJ/mol -> kcal/mol #scale cos - - fconstAmb = (2*fconst/(2**2)) # Vn = (2*fcte/(n**2)) - frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") #Atom# 1 - frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") #Atom 2 - frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}-") #Atom 3 - frcmod.write(f"{unique_at[ids[3]][atom_ids[ids[3]]]:<2}") #Atom 4 - frcmod.write(f" 1 {fconstAmb:>15.2f} {equ:>15.2f} 2\n") #IDIVF=1, fcte, angle, n=2 (n*angle) - - if len(terms['dihedral/flexible']) > 0: - - for dihed in terms['dihedral/flexible']: - ids = dihed.atomids + 1 - c = dihed.equ - - for n in range(3, 0, -1): - if n%2==1: - equ = 180.0 - else: - equ = 0.00 - frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") #Atom 1 - frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") #Atom 2 - frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}-") #Atom 3 - frcmod.write(f"{unique_at[ids[3]][atom_ids[ids[3]]]:<2}") #Atom 4 - frcmod.write(f" 1 {c[n]:>15.2f} {equ:>15.2f} -{n+1}\n") #IDIVF=1, fcte, angle, n=1 (n*angle) - - # Last term without -n - frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") #Atom 1 - frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") #Atom 2 - frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}-") #Atom 3 - frcmod.write(f"{unique_at[ids[3]][atom_ids[ids[3]]]:<2}") #Atom 4 - frcmod.write(f" 1 {c[0]:>15.2f} {0.00:>15.2f} 1\n") #IDIVF=1, fcte, angle, n=1 (n*angle) - - - # 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) # Degrees -> Degrees - fconst = dihed.fconst * (0.239005)/2 # kJ/mol -> kcal/mol -> f/2 as is Vn/2 in amber - frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") #Atom 1 - frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") #Atom 2 - frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}-") #Atom 3 - frcmod.write(f"{unique_at[ids[3]][atom_ids[ids[3]]]:<2}") #Atom 4 - frcmod.write(f" 1 {fconst:>15.2f} {equ:>15.2f} 2\n") #IDIVF=1, fcte, angle, n=2 (n*angle) - - def write_frcmod_nonbond(self, frcmod, non_bonded, atom_ids, unique_at): - gro_atomtypes, gro_nonbonded, gro_1_4 = self.convert_to_gromacs_nonbonded(non_bonded) - 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"{unique_at[i][atom_ids[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") - - -############################################## -# GET ATOM TYPES AND IDS # -############################################## - def get_atom_types(self, topo, non_bonded): - atom_ids={} #dictonary containing original atom types - unique_at={} #dictonary containing new unique atom types - unique_masses={} - - ascii_lowercase = list(string.ascii_lowercase) - ascii_digits = list(string.digits) - ascii_uppercase = list(string.ascii_uppercase) - - for i, (lj_type, a_name) in enumerate(zip(non_bonded.lj_types, self.atom_names), start=1): - atom_ids[i]=lj_type - if self.n_atoms < 36: - if i <= 10: - unique_at[i]={atom_ids[i] : "Q{}".format(ascii_digits[i-1])} - elif i > 10: - unique_at[i]={atom_ids[i] : "Q{}".format(ascii_lowercase[i-11])} - else: - unique_at[i]={atom_ids[i] : "{}{}".format(ascii_uppercase[i-1], ascii_lowercase[i-1])} - - ## DEBUG : check correspondance between atom_ids and unique_at - #for i, lj_type in enumerate(zip(non_bonded.lj_types), start=1): - # print(unique_at) - # for i in range(1,self.n_atoms+1,1): - # print (i, atom_ids[i], unique_at[i][atom_ids[i]]) - # print ("------------------") - - return atom_ids, unique_at - -############################################## - - def write_gromacs(self, directory, mol, coords): - self.write_itp(mol, directory) - self.write_top(directory) - self.write_gro(directory, coords, mol.non_bonded.alpha_map) - - def write_top(self, directory): - with open(f"{directory}/gas{self.polar_title}.top", "w") as top: - # defaults - top.write("\n[ defaults ]\n") - top.write(";nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ\n") - top.write(f"{1:>7}{self.comb_rule:>12}{'yes':>12}{self.fudge_lj:>10}{self.fudge_q:>10}" - "\n\n\n") - - top.write("; Include the molecule ITP\n") - top.write(f'#include "./{self.mol_name}_qforce{self.polar_title}.itp"\n\n\n') - - size = len(self.mol_name) - top.write("[ system ]\n") - top.write(f"; {' '*(size-6)}name\n") - top.write(f"{' '*(6-size)}{self.mol_name}\n\n\n") - - top.write("[ molecules ]\n") - top.write(f"; {' '*(size-10)}compound n_mol\n") - top.write(f"{' '*(10-size)}{self.mol_name} 1\n") - - def write_gro(self, directory, coords, alpha_map, box=[20., 20., 20.]): - n_atoms = self.n_atoms - if self.polar: - n_atoms += len(alpha_map.keys()) - coords_nm = coords*0.1 - with open(f"{directory}/gas{self.polar_title}.gro", "w") as gro: - gro.write(f"{self.mol_name}\n") - gro.write(f"{n_atoms:>6}\n") - for i, (a_name, coord) in enumerate(zip(self.atom_names, coords_nm), start=1): - gro.write(f"{1:>5}{self.residue:<5}") - gro.write(f"{a_name:>5}{i:>5}") - gro.write(f"{coord[0]:>8.3f}{coord[1]:>8.3f}{coord[2]:>8.3f}\n") - if self.polar: - for i, (atom, drude) in enumerate(alpha_map.items(), start=1): - gro.write(f"{2:>5}{self.residue:<5}{f'D{i}':>5}{drude+1:>5}") - gro.write(f"{coords_nm[atom][0]:>8.3f}{coords_nm[atom][1]:>8.3f}") - gro.write(f"{coords_nm[atom][2]:>8.3f}\n") - gro.write(f'{box[0]:>12.5f}{box[1]:>12.5f}{box[2]:>12.5f}\n') - - 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) - self.write_itp_atoms_and_molecule(itp, mol.non_bonded) - if self.polar: - self.write_itp_polarization(itp, mol.non_bonded) - self.write_itp_bonds(itp, mol.terms, mol.non_bonded.alpha_map) - self.write_itp_angles(itp, mol.terms) - self.write_itp_dihedrals(itp, mol.terms) - self.write_itp_pairs(itp) - self.write_itp_exclusions(itp) - itp.write('\n') - - def convert_to_gromacs_nonbonded(self, non_bonded): - a_types, nb_pairs, nb_1_4 = {}, {}, {} - - for pair, val in non_bonded.lj_pairs.items(): - if non_bonded.comb_rule != 1: - if val[0] == 0: - a, b = 0, 0 - else: - a, b = calc_sigma_epsilon(val[0], val[1]) - a *= 0.1 - else: - a = val[0] * 1e-6 - b = val[1] * 1e-12 - - if pair[0] == pair[1]: - a_types[pair[0]] = [a, b] - - nb_pairs[pair] = [a, b] - - for pair, val in non_bonded.lj_1_4.items(): - if non_bonded.comb_rule != 1: - a, b = calc_sigma_epsilon(val[0], val[1]) - a *= 0.1 - else: - a = val[0] * 1e-6 - b = val[1] * 1e-12 - - nb_1_4[pair] = [a, b] - - return a_types, nb_pairs, nb_1_4 - - def write_itp_atoms_and_molecule(self, itp, non_bonded): - gro_atomtypes, gro_nonbonded, gro_1_4 = self.convert_to_gromacs_nonbonded(non_bonded) - - # atom types - itp.write("\n[ atomtypes ]\n") - if self.comb_rule == 1: - itp.write("; name at_num mass charge type c6 c12\n") - else: - 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_params[0]:>12.5e} {lj_params[1]:>12.5e}\n') - - if self.polar: - itp.write(f'{"DP":>8} {0:>8.4f} {0:>8.4f} {"S":>2} {0:>12.5e} {0:>12.5e}\n') - - # non-bonded pair types - itp.write("\n[ nonbond_params ]\n") - if self.comb_rule == 1: - itp.write("; name1 name2 type c6 c12\n") - else: - itp.write("; name1 name2 type sigma epsilon\n") - - for lj_types, lj_params in gro_nonbonded.items(): - itp.write(f'{lj_types[0]:>8} {lj_types[1]:>8} 1') - itp.write(f'{lj_params[0]:>15.5e} {lj_params[1]:>15.5e}\n') - - # Write 1-4 pair types - if self.n_excl == 2 and gro_1_4 != {}: - itp.write("\n[ pairtypes ]\n") - if self.comb_rule == 1: - itp.write("; name1 name2 type c6 c12\n") - else: - itp.write("; name1 name2 type sigma epsilon\n") - for lj_types, lj_params in gro_1_4.items(): - itp.write(f'{lj_types[0]:>8} {lj_types[1]:>8} 1') - itp.write(f'{lj_params[0]:>15.5e} {lj_params[1]:>15.5e}\n') - - # molecule type - space = " "*(len(self.mol_name)-5) - itp.write("\n[ moleculetype ]\n") - itp.write(f";{space}name nrexcl\n") - itp.write(f"{self.mol_name}{3:>7}\n") - - # atoms - itp.write("\n[ atoms ]\n") - itp.write("; nr type resnr resnm atom cgrp charge mass\n") - for i, (lj_type, a_name, q, mass) in enumerate(zip(non_bonded.lj_types, self.atom_names, - self.q, self.masses), start=1): - itp.write(f'{i:>5}{lj_type:>9}{1:>6}{self.residue:>6}{a_name:>7}{i:>5}{q:>11.5f}') - itp.write(f'{mass:>10.5f}\n') - - if self.polar: - for i, (atom, drude) in enumerate(non_bonded.alpha_map.items(), start=1): - itp.write(f'{drude+1:>5}{"DP":>9}{2:>6}{"DRU":>6}{f"D{i}":>7}{atom+1:>5}') - itp.write(f'{-8.:>11.5f}{0.:>10.5f}\n') - - def write_itp_polarization(self, itp, non_bonded): - # polarization - itp.write("\n[ polarization ]\n") - itp.write("; i j f alpha\n") - for atom, drude in non_bonded.alpha_map.items(): - alpha = non_bonded.alpha[atom]*1e-3 - itp.write(f"{atom+1:>6}{drude+1:>6}{1:>6}{alpha:>14.8f}\n") - - # # thole polarization - # if self.thole != []: - # itp.write("\n[ thole_polarization ]\n") - # itp.write("; ai di aj dj f a alpha(i) " - # "alpha(j)\n") - # for tho in self.thole: - # itp.write("{:>6}{:>6}{:>6}{:>6}{:>4}{:>7.2f}{:>14.8f}{:>14.8f}\n".format(*tho)) - - def write_itp_pairs(self, itp): - if self.pairs != []: - itp.write("\n[ pairs ]\n") - itp.write("; ai aj func \n") - for pair in self.pairs: - itp.write(f"{pair[0]+1:>6}{pair[1]+1:>6}{1:>6}\n") - - def write_itp_bonds(self, itp, terms, alpha_map): - itp.write("\n[ bonds ]\n") - itp.write("; ai aj f r0 kb\n") - for bond in terms['bond']: - ids = bond.atomids + 1 - equ = bond.equ * 0.1 - fconst = bond.fconst * 100 - itp.write(f'{ids[0]:>6}{ids[1]:>6}{1:>6}{equ:>10.5f}{fconst:>10.0f}\n') - - if self.polar: - itp.write('; ai aj f - polar connections\n') - for bond in terms['bond']: - a1, a2 = bond.atomids - - if a2 in alpha_map.keys(): - itp.write(f'{a1+1:>6}{alpha_map[a2]+1:>6}{5:>6}\n') - if a1 in alpha_map.keys(): - itp.write(f'{a2+1:>6}{alpha_map[a1]+1:>6}{5:>6}\n') - - def write_itp_angles(self, itp, terms): - itp.write("\n[ angles ]\n") - itp.write("; ai aj ak f th0 kth\n") - for angle in terms['angle']: - ids = angle.atomids + 1 - equ = np.degrees(angle.equ) - fconst = angle.fconst - - 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: - itp.write(f'{ids[0]:>6}{ids[1]:>6}{ids[2]:>6}{1:>6}{equ:>11.3f}{fconst:>13.3f}\n') - else: - urey_equ = urey[0].equ * 0.1 - urey_fconst = urey[0].fconst * 100 - itp.write(f'{ids[0]:>6}{ids[1]:>6}{ids[2]:>6}{5:>6}{equ:>11.3f}{fconst:>13.3f}' - f'{urey_equ:>10.5f}{urey_fconst:>13.3f}\n') - - def write_itp_dihedrals(self, itp, terms): - if len(terms['dihedral']) > 0: - itp.write("\n[ dihedrals ]\n") - - # rigid dihedrals - if len(terms['dihedral/rigid']) > 0: - itp.write("; rigid dihedrals \n") - itp.write("; ai aj ak al f th0 kth\n") - - for dihed in terms['dihedral/rigid']: - ids = dihed.atomids + 1 - equ = np.degrees(dihed.equ) - fconst = dihed.fconst - - itp.write(f'{ids[0]:>6}{ids[1]:>6}{ids[2]:>6}{ids[3]:>6}{2:>6}{equ:>11.3f}') - itp.write(f'{fconst:>13.3f}\n') - - # improper dihedrals - if len(terms['dihedral/improper']) > 0: - itp.write("; improper dihedrals \n") - itp.write("; ai aj ak al f th0 kth\n") - - for dihed in terms['dihedral/improper']: - ids = dihed.atomids + 1 - equ = np.degrees(dihed.equ) - fconst = dihed.fconst - - itp.write(f'{ids[0]:>6}{ids[1]:>6}{ids[2]:>6}{ids[3]:>6}{2:>6}{equ:>11.3f}') - itp.write(f'{fconst:>13.3f}\n') - - # flexible dihedrals - if len(terms['dihedral/flexible']) > 0: - itp.write("; flexible dihedrals \n") - itp.write('; ai aj ak al f c0 c1 c2 c3') - itp.write(' c4 c5\n') - - for dihed in terms['dihedral/flexible']: - ids = dihed.atomids + 1 - c = dihed.equ - - itp.write(f'{ids[0]:>6}{ids[1]:>6}{ids[2]:>6}{ids[3]:>6}{3:>6}{c[0]:>11.3f}') - itp.write(f'{c[1]:>11.3f}{c[2]:>11.3f}{c[3]:>11.3f}{c[4]:>11.3f}{c[5]:>11.3f}\n') - - # inversion dihedrals - if len(terms['dihedral/inversion']) > 0: - itp.write("; inversion dihedrals \n") - itp.write('; ai aj ak al f c0 c1 c2\n') - - for dihed in terms['dihedral/inversion']: - ids = dihed.atomids + 1 - # itp.write(f'{ids[0]:>6}{ids[1]:>6}{ids[2]:>6}{ids[3]:>6}{1:>6}' - # f'{0.0:>11.3f}{dihed.fconst:>11.3f}{3:>11}\n') - c0, c1, c2 = convert_to_inversion_rb(dihed.fconst, dihed.equ) - itp.write(f'{ids[0]:>6}{ids[1]:>6}{ids[2]:>6}{ids[3]:>6}{3:>6}' - f'{c0:>11.3f}{c1:>11.3f}{c2:>11.3f}{0:>11.1f}{0:>11.1f}{0:>11.1f}\n') - - def write_itp_exclusions(self, itp): - if any(len(exclusion) > 0 for exclusion in self.exclusions): - itp.write("\n[ exclusions ]\n") - for i, exclusion in enumerate(self.exclusions): - if len(exclusion) > 0: - exclusion = sorted(set(exclusion)) - itp.write("{} ".format(i+1)) - itp.write(("{} "*len(exclusion)).format(*exclusion)) - itp.write("\n") - - def make_pairs(self, neighbors, non_bonded): - pairs, polar_pairs = [], [] - - for pair in non_bonded.pairs: - pairs.append(sorted(pair)) - - if self.n_excl == 2: - for i in range(self.n_atoms): - for neigh in neighbors[2][i]: - if (i < neigh and [i, neigh] not in pairs and (i, neigh) - not in non_bonded.exclusions): - pairs.append([i, neigh]) - - if self.polar: - for a1, a2 in pairs: - if a2 in non_bonded.alpha_map.keys(): - polar_pairs.append([a1, non_bonded.alpha_map[a2]]) - if a1 in non_bonded.alpha_map.keys(): - polar_pairs.append([a2, non_bonded.alpha_map[a1]]) - if a1 in non_bonded.alpha_map.keys() and a2 in non_bonded.alpha_map.keys(): - polar_pairs.append([non_bonded.alpha_map[a1], non_bonded.alpha_map[a2]]) - - return pairs+polar_pairs - - def make_exclusions(self, non_bonded, neighbors, exclude_all): - exclusions = [[] for _ in range(self.n_atoms)] - - # input exclusions - for exclusion in non_bonded.exclusions: - exclusions[exclusion[0]].append(exclusion[1]+1) - - # exclusions for input pairs - for pair in non_bonded.pairs: - exclusions[pair[0]].append(pair[1]+1) - - # fragment capping atom exclusions - for i in exclude_all: - exclusions[i].extend(np.arange(1, self.n_atoms+1)) - - if self.polar: - exclusions = self.polarize_exclusions(non_bonded.alpha_map, non_bonded.exclusions, - neighbors, exclude_all, exclusions) - - return exclusions - - def polarize_exclusions(self, alpha_map, input_exclusions, neighbors, exclude_all, exclusions): - n_polar_atoms = len(alpha_map.keys()) - exclusions.extend([[] for _ in range(n_polar_atoms)]) - - # input exclusions - for exclu in input_exclusions: - if exclu[0] in alpha_map.keys(): - exclusions[alpha_map[exclu[0]]].append(exclu[1]+1) - if exclu[1] in alpha_map.keys(): - exclusions[alpha_map[exclu[1]]].append(exclu[0]+1) - if exclu[0] in alpha_map.keys() and exclu[1] in alpha_map.keys(): - exclusions[alpha_map[exclu[0]]].append(alpha_map[exclu[1]]+1) - - # fragment capping atom exclusions - for i in exclude_all: - exclusions[i].extend(np.arange(self.n_atoms+1, self.n_atoms+n_polar_atoms+1)) - if i in alpha_map.keys(): - exclusions[alpha_map[i]].extend(np.arange(1, self.n_atoms+n_polar_atoms+1)) - - return exclusions - - def get_atom_names(self): - atom_names = [] - atom_dict = {} - - for i, elem in enumerate(self.elements): - sym = ATOM_SYM[elem] - if sym not in atom_dict.keys(): - atom_dict[sym] = 1 - else: - atom_dict[sym] += 1 - atom_names.append(f'{sym}{atom_dict[sym]}') - return atom_names - - def add_restraints(self, restraints, directory, fc=1000): - with open(f"{directory}/{self.mol_name}_qforce{self.polar_title}.itp", "a") as itp: - itp.write("[ dihedral_restraints ]\n") - itp.write("; ai aj ak al type phi dp kfac\n") - for restraint in restraints: - a1, a2, a3, a4 = restraint[0]+1 - phi = np.degrees(restraint[1]) - itp.write(f'{a1:>5}{a2:>5}{a3:>5}{a4:>5}{1:>5} {phi:>10.4f} 0.0 {fc}\n') - - def set_charge(self, non_bonded): - q = np.copy(non_bonded.q) - if self.polar: - 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])])) From a60171e27a05405c1211b1aa8be915336e9cb1b3 Mon Sep 17 00:00:00 2001 From: Mateus Zanotto <39527468+mateuszanotto@users.noreply.github.com> Date: Tue, 13 Sep 2022 13:10:17 -0300 Subject: [PATCH 03/28] Delete initialize.py --- initialize.py | 178 -------------------------------------------------- 1 file changed, 178 deletions(-) delete mode 100644 initialize.py diff --git a/initialize.py b/initialize.py deleted file mode 100644 index 5bea8c7..0000000 --- a/initialize.py +++ /dev/null @@ -1,178 +0,0 @@ -import os -import shutil -from io import StringIO -from types import SimpleNamespace -import pkg_resources -# -from colt import Colt -# -from .qm.qm import QM, implemented_qm_software -from .molecule.terms import Terms -from .dihedral_scan import DihedralScan -from .misc import LOGO - - -class Initialize(Colt): - _user_input = """ -[ff] -# MD software to use -# Currently amber supports only gaff2 atom types -forcefield = 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 - -# Number of first n neighbors to exclude in the forcefield -n_excl = 2 :: int :: [2, 3] - -# Lennard jones method for the forcefield -lennard_jones = opls_auto :: str :: [gromos_auto, gromos, opls_auto, opls, gaff, gaff2, ext] - -# Use externally provided point charges in the file "ext_q" in the job directyory -ext_charges = no :: bool - -# Scale QM charges to account for condensed phase polarization (should be set to 1 for gas phase) -charge_scaling = 1.2 :: float - -# If user chooses ext_charges=True, by default fragments will still use the chosen QM method for -# determining fragment charges. This is to avoid spuriously high charges on capping hydrogens. -# However, using QM charges for fragments and ext_charges for the molecule can also result in -# inaccuracies if these two charges are very different. -use_ext_charges_for_frags = no :: bool - -# Additional exclusions (GROMACS format) -exclusions = :: literal - -# Switch standard non-bonded interactions between two atoms to pair interactions -# (provide atom pairs in each row) -pairs = :: literal - -# Path for the external FF library for Lennard-Jones parameters (GROMACS format). -ext_lj_lib = :: folder, optional - -# Lennard-Jones fudge parameter for 1-4 interactions for when "lennard_jones" is set to "ext". -ext_lj_fudge = :: float, optional - -# Coulomb fudge parameter for 1-4 interactions for when "lennard_jones" is set to "ext". -ext_q_fudge = :: float, optional - -# Lennard-Jones combinations rules for when "lennard_jones" is set to "ext" (GROMACS numbering). -ext_comb_rule = :: int, optional :: [1, 2, 3] - -# Name of the atom type for capping hydrogens for when "lennard_jones" is set to "ext" -ext_h_cap = :: str, optional - -# Set all dihedrals as rigid (no dihedral scans) -all_rigid = no :: bool - -# Use D4 method -_d4 = no :: bool - -# Residue name printed on the force field file (Max 5 characters) -res_name = MOL :: str - -# Polarize a coordinate file and quit (requires itp_file) -_polarize = no :: bool - -# Name of itp file (only needed for polarize option) -_itp_file = itp_file_missing :: str - -# Make a polarizable FF -_polar = no :: bool - -# Scale the C6 dispersion interactions in the polarizable version of the FF -_polar_c6_scale = 0.8 :: float - -# Specifically not scale some of the atoms -_polar_not_scale_c6 = :: literal - -# Manual polarizabilities in the file ext_alpha -_ext_alpha = no :: bool - -""" - - @staticmethod - def _set_config(config): - config['qm'].update(config['qm']['software']) - config['qm'].update({'software': config['qm']['software'].value}) - config.update({key: SimpleNamespace(**val) for key, val in config.items()}) - return SimpleNamespace(**config) - - @classmethod - def _extend_user_input(cls, questions): - questions.generate_block("qm", QM.colt_user_input) - questions.generate_block("scan", DihedralScan.colt_user_input) - questions.generate_cases("software", {key: software.colt_user_input for key, software in - implemented_qm_software.items()}, block='qm') - questions.generate_block("terms", Terms.get_questions()) - - @classmethod - def from_config(cls, config): - return cls._set_config(config) - - @staticmethod - def set_basis(value): - if value.endswith('**'): - return f'{value[:-2]}(D,P)'.upper() - if value.endswith('*'): - return f'{value[:-1]}(D)'.upper() - return value.upper() - - @staticmethod - def set_dispersion(value): - if value.lower() in ["no", "false", "n", "f"]: - return False - return value.upper() - - -def _get_job_info(filename): - job = {} - filename = filename.rstrip('/') - base = os.path.basename(filename) - path = os.path.dirname(filename) - if path != '': - path = f'{path}/' - - if os.path.isfile(filename): - job['coord_file'] = filename - job['name'] = base.split('.')[0] - else: - job['coord_file'] = False - job['name'] = base.split('_qforce')[0] - - job['dir'] = f'{path}{job["name"]}_qforce' - job['frag_dir'] = f'{job["dir"]}/fragments' - job['md_data'] = pkg_resources.resource_filename('qforce', 'data') - os.makedirs(job['dir'], exist_ok=True) - return SimpleNamespace(**job) - - -def _check_and_copy_settings_file(job_dir, config_file): - """ - If options are provided as a file, copy that to job directory. - If options are provided as StringIO, write that to job directory. - """ - - settings_file = os.path.join(job_dir, 'settings.ini') - - if config_file is not None: - if isinstance(config_file, StringIO): - with open(settings_file, 'w') as fh: - config_file.seek(0) - fh.write(config_file.read()) - else: - shutil.copy2(config_file, settings_file) - - return settings_file - - -def initialize(filename, config_file, presets=None): - #print(LOGO) - - job_info = _get_job_info(filename) - settings_file = _check_and_copy_settings_file(job_info.dir, config_file) - - config = Initialize.from_questions(config=settings_file, presets=presets, check_only=True) - - return config, job_info From b02618e4537d2aa0ed46a6ff6f05ea2e8b707376 Mon Sep 17 00:00:00 2001 From: Mateus Zanotto <39527468+mateuszanotto@users.noreply.github.com> Date: Tue, 13 Sep 2022 13:10:26 -0300 Subject: [PATCH 04/28] Delete main.py --- main.py | 95 --------------------------------------------------------- 1 file changed, 95 deletions(-) delete mode 100644 main.py diff --git a/main.py b/main.py deleted file mode 100644 index b14fc33..0000000 --- a/main.py +++ /dev/null @@ -1,95 +0,0 @@ -from .polarize import polarize -from .initialize import initialize -from .qm.qm import QM -from .qm.qm_base import HessianOutput -from .forcefield import ForceField -from .molecule import Molecule -from .fragment import fragment -from .dihedral_scan import DihedralScan -from .frequencies import calc_qm_vs_md_frequencies -from .hessian import fit_hessian - -from .misc import check_if_file_exists, LOGO -from colt import from_commandline -from colt.validator import Validator - - -# define new validator -Validator.overwrite_validator("file", check_if_file_exists) - - -@from_commandline(""" -# Input coordinate file mol.ext (ext: pdb, xyz, gro, ...) -# or directory (mol or mol_qforce) name. -file = :: file - -# File name for the optional options. -options = :: file, optional, alias=o -""", description={ - 'logo': LOGO, - 'alias': 'qforce', - 'arg_format': { - 'name': 12, - 'comment': 60, - }, -}) -def run(file, options): - run_qforce(input_arg=file, config=options) - - -def run_qforce(input_arg, ext_q=None, ext_lj=None, config=None, presets=None): - config, job = initialize(input_arg, config, presets) - - if config.ff._polarize: - polarize(job, config.ff) - - qm = QM(job, config.qm) - qm_hessian_out = qm.read_hessian() - - mol = Molecule(config, job, qm_hessian_out, ext_q, ext_lj) - - md_hessian = fit_hessian(config.terms, mol, qm_hessian_out) - - if len(mol.terms['dihedral/flexible']) > 0 and config.scan.do_scan: - fragments = fragment(mol, qm, job, config) - DihedralScan(fragments, mol, job, config) - - calc_qm_vs_md_frequencies(job, qm_hessian_out, md_hessian) - ff = ForceField(job.name, config, mol, mol.topo.neighbors) - - forcefield = config.ff.forcefield - if forcefield in ["gromacs"]: - ff.write_gromacs(job.dir, mol, qm_hessian_out.coords) - - if forcefield in ["amber"]: - ff.write_amber(job.dir, mol, qm_hessian_out.coords) - - print_outcome(job.dir) - - -def run_hessian_fitting_for_external(job_dir, qm_data, ext_q=None, ext_lj=None, - config=None, presets=None): - config, job = initialize(job_dir, config, presets) - - qm_hessian_out = HessianOutput(config.qm.vib_scaling, **qm_data) - - mol = Molecule(config, job, qm_hessian_out, ext_q, ext_lj) - - md_hessian = fit_hessian(config.terms, mol, qm_hessian_out) - 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) - - print_outcome(job.dir) - - return mol.terms - - -def print_outcome(job_dir): - 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('- 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.') From 627d65fd47ca842a5661a8e196a1319a6f00174e Mon Sep 17 00:00:00 2001 From: Mateus Zanotto <39527468+mateuszanotto@users.noreply.github.com> Date: Tue, 13 Sep 2022 13:11:41 -0300 Subject: [PATCH 05/28] Add files via upload --- qforce/forcefield.py | 334 +++++++++++++++++++++++++++++++++++++------ qforce/initialize.py | 8 +- qforce/main.py | 8 +- 3 files changed, 301 insertions(+), 49 deletions(-) diff --git a/qforce/forcefield.py b/qforce/forcefield.py index 2c71264..a383ac6 100644 --- a/qforce/forcefield.py +++ b/qforce/forcefield.py @@ -5,6 +5,9 @@ from .forces import convert_to_inversion_rb from .misc import LOGO_SEMICOL +from .misc import LOGO_HASH +import string + class ForceField(): def __init__(self, job_name, config, mol, neighbors, exclude_all=[]): @@ -29,6 +32,236 @@ def __init__(self, job_name, config, mol, neighbors, exclude_all=[]): else: self.polar_title = '' +## Amber + def write_amber(self, directory, mol, coords): + atom_ids, unique_at =self.get_atom_types(mol.topo, mol.non_bonded) + self.write_mol2(directory, mol, coords, atom_ids, unique_at) + self.write_frcmod(directory, mol, coords, atom_ids, unique_at) + + def write_mol2(self, directory, mol, coords, atom_ids, unique_at): + 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, atom_ids, unique_at) + self.write_mol2_bond(mol2, mol.topo, mol.terms, atom_ids, unique_at) + + def write_frcmod(self, directory, mol, coords, atom_ids, unique_at): + with open(f"{directory}/{self.mol_name}_qforce{self.polar_title}.frcmod", "w") as frcmod: + self.write_frcmod_mass(frcmod, mol.non_bonded, atom_ids, unique_at) + self.write_frcmod_bonds(frcmod, mol.terms, mol.non_bonded.alpha_map, atom_ids, unique_at) + self.write_frcmod_angles(frcmod, mol.terms, atom_ids, unique_at) + self.write_frcmod_dihedrals(frcmod, mol.terms, atom_ids, unique_at) + self.write_frcmod_nonbond(frcmod, mol.non_bonded, atom_ids, unique_at) + + +# TLeap is a helper program from AMBER to create the topologies. +# Would be good for the final user to create automatically a +# Tleap script with the new atom types. + + # def write_tleap_script(self, directory): + # with open(f"{directory}/tleap_script_{self.mol_name}.in", "w") as tleap: + # tleap.write(f""" + # # To run the script use tleap on amberTools and the following command: + # # $ tleap -f tleap_script_{self.mol_name}.in > tleap_script_{self.mol_name}.o + # + # source leaprc.protein.ff19SB\nsource leaprc.water.opc + # source leaprc.gaff\n# You can source other force fields here if necessary + # loadmol2 {self.mol_name}_qforce{self.polar_title}.mol2 + # + # loadamberparams {self.mol_name}_qforce{self.polar_title}.frcmod + # + # ## Add any other command, for example "solvateoct MOL TIP3PBOX 10.0" + # + # MOL = loadPDB 1lyd.solv.pdb + # saveamberparm MOL {self.mol_name}_qforce{self.polar_title}.prmtop + # saveamberparm MOL {self.mol_name}_qforce{self.polar_title}.inpcrd + # savepdb MOL {self.mol_name}_qforce{self.polar_title}.pdb """) + +############################################## +# MOL2 WRITING # +############################################## + + def write_mol2_title(self, mol2): + mol2.write(LOGO_HASH) + 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(f"esp")## type of charge + mol2.write("\n\n") + + def write_mol2_atom(self, mol2, topo, coords, non_bonded, atom_ids, unique_at): + mol2.write(f"@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} {coords[i_idx-1][1]:10.4f} {coords[i_idx-1][2]:10.4f}") + mol2.write(f" {unique_at[i_idx][atom_ids[i_idx]]} {1} {self.mol_name} {q:10.6f}\n") + + def write_mol2_bond(self, mol2, topo, terms, atom_ids, unique_at): + n_bonds = 1 + mol2.write(f"@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(f"@SUBSTRUCTURE\n") + mol2.write(f"{1:5} {self.mol_name} {1:8} TEMP {0:>8d} **** **** {0:5} ROOT\n") + +############################################## +# FRCMOD WRITING # +############################################## + + def write_frcmod_mass(self, frcmod, non_bonded, atom_ids, unique_at): + frcmod.write(f'{self.mol_name} - frcmod generated by QForce\n') + frcmod.write(f'MASS\n') + for i, mass in enumerate((self.masses), start=1): + frcmod.write(f"{unique_at[i][atom_ids[i]]} {mass}\n") + + def write_frcmod_bonds(self, frcmod, terms, alpha_map, atom_ids, unique_at): + frcmod.write("\nBOND\n") + for bond in terms['bond']: + ids = bond.atomids + 1 + fconst = bond.fconst * (0.239005)/2 # kJ/mol/A^2 -> kcal/mol/A^2 + equ = bond.equ # Ang + frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") + frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}") + frcmod.write(f"{fconst:>10.2f}{equ:>10.2f}\n") + + def write_frcmod_angles(self, frcmod, terms, atom_ids, unique_at): + frcmod.write("\nANGLE\n") + for angle in terms['angle']: + ids = angle.atomids + 1 + fconst = angle.fconst * (0.239005)/2 # kJ/mol/rad^2 -> kcal/mol/rad^2 + equ = np.degrees(angle.equ) # Degrees -> Degrees + + 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"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") + frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") + frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}") + frcmod.write(f"{fconst:>10.2f}{equ:>10.2f}\n") + else: + urey_equ = urey[0].equ + urey_fconst = urey[0].fconst * (0.239005)/2 # kJ/mol/rad^2 -> kcal/mol/rad^2 + frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") + frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") + frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}") + frcmod.write(f"{fconst:>10.2f}{equ:>10.2f}\n") + + def write_frcmod_dihedrals(self, frcmod, terms, atom_ids, unique_at): + 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 * (0.239005) # kJ/mol -> kcal/mol #scale cos + + fconstAmb = (2*fconst/(2**2)) # Vn = (2*fcte/(n**2)) + frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") #Atom# 1 + frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") #Atom 2 + frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}-") #Atom 3 + frcmod.write(f"{unique_at[ids[3]][atom_ids[ids[3]]]:<2}") #Atom 4 + frcmod.write(f" 1 {fconstAmb:>15.2f} {equ:>15.2f} 2\n") #IDIVF=1, fcte, angle, n=2 (n*angle) + + if len(terms['dihedral/flexible']) > 0: + + for dihed in terms['dihedral/flexible']: + ids = dihed.atomids + 1 + c = dihed.equ + + for n in range(3, 0, -1): + if n%2==1: + equ = 180.0 + else: + equ = 0.00 + frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") #Atom 1 + frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") #Atom 2 + frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}-") #Atom 3 + frcmod.write(f"{unique_at[ids[3]][atom_ids[ids[3]]]:<2}") #Atom 4 + frcmod.write(f" 1 {c[n]:>15.2f} {equ:>15.2f} -{n+1}\n") #IDIVF=1, fcte, angle, n=1 (n*angle) + + # Last term without -n + frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") #Atom 1 + frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") #Atom 2 + frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}-") #Atom 3 + frcmod.write(f"{unique_at[ids[3]][atom_ids[ids[3]]]:<2}") #Atom 4 + frcmod.write(f" 1 {c[0]:>15.2f} {0.00:>15.2f} 1\n") #IDIVF=1, fcte, angle, n=1 (n*angle) + + + # 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) # Degrees -> Degrees + fconst = dihed.fconst * (0.239005)/2 # kJ/mol -> kcal/mol -> f/2 as is Vn/2 in amber + frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") #Atom 1 + frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") #Atom 2 + frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}-") #Atom 3 + frcmod.write(f"{unique_at[ids[3]][atom_ids[ids[3]]]:<2}") #Atom 4 + frcmod.write(f" 1 {fconst:>15.2f} {equ:>15.2f} 2\n") #IDIVF=1, fcte, angle, n=2 (n*angle) + + def write_frcmod_nonbond(self, frcmod, non_bonded, atom_ids, unique_at): + gro_atomtypes, gro_nonbonded, gro_1_4 = self.convert_to_gromacs_nonbonded(non_bonded) + 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"{unique_at[i][atom_ids[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") + + +############################################## +# GET ATOM TYPES AND IDS # +############################################## + def get_atom_types(self, topo, non_bonded): + atom_ids={} #dictonary containing original atom types + unique_at={} #dictonary containing new unique atom types + unique_masses={} + + ascii_lowercase = list(string.ascii_lowercase) + ascii_digits = list(string.digits) + ascii_uppercase = list(string.ascii_uppercase) + + for i, (lj_type, a_name) in enumerate(zip(non_bonded.lj_types, self.atom_names), start=1): + atom_ids[i]=lj_type + if self.n_atoms < 36: + if i <= 10: + unique_at[i]={atom_ids[i] : "Q{}".format(ascii_digits[i-1])} + elif i > 10: + unique_at[i]={atom_ids[i] : "Q{}".format(ascii_lowercase[i-11])} + else: + unique_at[i]={atom_ids[i] : "{}{}".format(ascii_uppercase[i-1], ascii_lowercase[i-1])} + + ## DEBUG : check correspondance between atom_ids and unique_at + #for i, lj_type in enumerate(zip(non_bonded.lj_types), start=1): + # print(unique_at) + # for i in range(1,self.n_atoms+1,1): + # print (i, atom_ids[i], unique_at[i][atom_ids[i]]) + # print ("------------------") + + return atom_ids, unique_at + +############################################## + def write_gromacs(self, directory, mol, coords): self.write_itp(mol, directory) self.write_top(directory) @@ -38,9 +271,9 @@ def write_top(self, directory): with open(f"{directory}/gas{self.polar_title}.top", "w") as top: # defaults top.write("\n[ defaults ]\n") - top.write("; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ\n") - top.write(f"{1:>8} {self.comb_rule:>12} {'yes':>12} {self.fudge_lj:>12} " - f"{self.fudge_q:>12}\n\n\n") + top.write(";nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ\n") + top.write(f"{1:>7}{self.comb_rule:>12}{'yes':>12}{self.fudge_lj:>10}{self.fudge_q:>10}" + "\n\n\n") top.write("; Include the molecule ITP\n") top.write(f'#include "./{self.mol_name}_qforce{self.polar_title}.itp"\n\n\n') @@ -164,24 +397,24 @@ def write_itp_atoms_and_molecule(self, itp, non_bonded): # atoms itp.write("\n[ atoms ]\n") - itp.write("; nr type resnr resnm atom cgrp charge mass\n") + itp.write("; nr type resnr resnm atom cgrp charge mass\n") for i, (lj_type, a_name, q, mass) in enumerate(zip(non_bonded.lj_types, self.atom_names, self.q, self.masses), start=1): - itp.write(f'{i:>5} {lj_type:>9} {1:>6} {self.residue:>6} {a_name:>7} {i:>5} ') - itp.write(f'{q:>11.5f} {mass:>10.5f}\n') + itp.write(f'{i:>5}{lj_type:>9}{1:>6}{self.residue:>6}{a_name:>7}{i:>5}{q:>11.5f}') + itp.write(f'{mass:>10.5f}\n') if self.polar: for i, (atom, drude) in enumerate(non_bonded.alpha_map.items(), start=1): - itp.write(f'{drude+1:>5} {"DP":>9} {2:>6} {"DRU":>6} {f"D{i}":>7} {atom+1:>5}') - itp.write(f'{-8.:>11.5f} {0.:>10.5f}\n') + itp.write(f'{drude+1:>5}{"DP":>9}{2:>6}{"DRU":>6}{f"D{i}":>7}{atom+1:>5}') + itp.write(f'{-8.:>11.5f}{0.:>10.5f}\n') def write_itp_polarization(self, itp, non_bonded): # polarization itp.write("\n[ polarization ]\n") - itp.write("; i j f alpha\n") + itp.write("; i j f alpha\n") for atom, drude in non_bonded.alpha_map.items(): alpha = non_bonded.alpha[atom]*1e-3 - itp.write(f"{atom+1:>6} {drude+1:>6} {1:>6} {alpha:>14.8f}\n") + itp.write(f"{atom+1:>6}{drude+1:>6}{1:>6}{alpha:>14.8f}\n") # # thole polarization # if self.thole != []: @@ -194,36 +427,32 @@ def write_itp_polarization(self, itp, non_bonded): def write_itp_pairs(self, itp): if self.pairs != []: itp.write("\n[ pairs ]\n") - itp.write("; ai aj func\n") + itp.write("; ai aj func \n") for pair in self.pairs: - itp.write(f"{pair[0]+1:>6} {pair[1]+1:>6} {1:>6}\n") + itp.write(f"{pair[0]+1:>6}{pair[1]+1:>6}{1:>6}\n") def write_itp_bonds(self, itp, terms, alpha_map): itp.write("\n[ bonds ]\n") - itp.write("; ai aj f r0 k_bond\n") + itp.write("; ai aj f r0 kb\n") for bond in terms['bond']: ids = bond.atomids + 1 equ = bond.equ * 0.1 fconst = bond.fconst * 100 - itp.write(f'{ids[0]:>6} {ids[1]:>6} {1:>6} {equ:>10.5f} {fconst:>10.0f}\n') + itp.write(f'{ids[0]:>6}{ids[1]:>6}{1:>6}{equ:>10.5f}{fconst:>10.0f}\n') if self.polar: - itp.write('; ai aj f - polar connections\n') + itp.write('; ai aj f - polar connections\n') for bond in terms['bond']: a1, a2 = bond.atomids if a2 in alpha_map.keys(): - itp.write(f'{a1+1:>6} {alpha_map[a2]+1:>6} {5:>6}\n') + itp.write(f'{a1+1:>6}{alpha_map[a2]+1:>6}{5:>6}\n') if a1 in alpha_map.keys(): - itp.write(f'{a2+1:>6} {alpha_map[a1]+1:>6} {5:>6}\n') + itp.write(f'{a2+1:>6}{alpha_map[a1]+1:>6}{5:>6}\n') def write_itp_angles(self, itp, terms): itp.write("\n[ angles ]\n") - itp.write("; ai aj ak f theta0 k_theta") - if self.urey: - itp.write(' r0 k_bond') - itp.write('\n') - + itp.write("; ai aj ak f th0 kth\n") for angle in terms['angle']: ids = angle.atomids + 1 equ = np.degrees(angle.equ) @@ -233,13 +462,12 @@ def write_itp_angles(self, itp, terms): urey = [term for term in terms['urey'] if np.array_equal(term.atomids, angle.atomids)] if not self.urey or len(urey) == 0: - itp.write(f'{ids[0]:>6} {ids[1]:>6} {ids[2]:>6} {1:>6} {equ:>10.3f} ' - f'{fconst:>10.3f}\n') + itp.write(f'{ids[0]:>6}{ids[1]:>6}{ids[2]:>6}{1:>6}{equ:>11.3f}{fconst:>13.3f}\n') else: urey_equ = urey[0].equ * 0.1 urey_fconst = urey[0].fconst * 100 - itp.write(f'{ids[0]:>6} {ids[1]:>6} {ids[2]:>6} {5:>6} {equ:>10.3f} ' - f'{fconst:>10.3f} {urey_equ:>10.5f} {urey_fconst:>10.1f}\n') + itp.write(f'{ids[0]:>6}{ids[1]:>6}{ids[2]:>6}{5:>6}{equ:>11.3f}{fconst:>13.3f}' + f'{urey_equ:>10.5f}{urey_fconst:>13.3f}\n') def write_itp_dihedrals(self, itp, terms): if len(terms['dihedral']) > 0: @@ -248,52 +476,54 @@ def write_itp_dihedrals(self, itp, terms): # rigid dihedrals if len(terms['dihedral/rigid']) > 0: itp.write("; rigid dihedrals \n") - itp.write("; ai aj ak al f theta0 k_theta\n") + itp.write("; ai aj ak al f th0 kth\n") for dihed in terms['dihedral/rigid']: ids = dihed.atomids + 1 equ = np.degrees(dihed.equ) fconst = dihed.fconst - itp.write(f'{ids[0]:>6} {ids[1]:>6} {ids[2]:>6} {ids[3]:>6} {2:>6} {equ:>11.3f} ') + itp.write(f'{ids[0]:>6}{ids[1]:>6}{ids[2]:>6}{ids[3]:>6}{2:>6}{equ:>11.3f}') itp.write(f'{fconst:>13.3f}\n') # improper dihedrals if len(terms['dihedral/improper']) > 0: itp.write("; improper dihedrals \n") - itp.write("; ai aj ak al f theta0 k_theta\n") + itp.write("; ai aj ak al f th0 kth\n") for dihed in terms['dihedral/improper']: ids = dihed.atomids + 1 equ = np.degrees(dihed.equ) fconst = dihed.fconst - itp.write(f'{ids[0]:>6} {ids[1]:>6} {ids[2]:>6} {ids[3]:>6} {2:>6} {equ:>11.3f} ') + itp.write(f'{ids[0]:>6}{ids[1]:>6}{ids[2]:>6}{ids[3]:>6}{2:>6}{equ:>11.3f}') itp.write(f'{fconst:>13.3f}\n') # flexible dihedrals if len(terms['dihedral/flexible']) > 0: itp.write("; flexible dihedrals \n") - itp.write('; ai aj ak al f c0 c1 c2') - itp.write(' c3 c4 c5\n') + itp.write('; ai aj ak al f c0 c1 c2 c3') + itp.write(' c4 c5\n') for dihed in terms['dihedral/flexible']: ids = dihed.atomids + 1 c = dihed.equ - itp.write(f'{ids[0]:>6} {ids[1]:>6} {ids[2]:>6} {ids[3]:>6} {3:>6} {c[0]:>11.3f} ') - itp.write(f'{c[1]:>11.3f} {c[2]:>11.3f} {c[3]:>11.3f} {c[4]:>11.3f} {c[5]:>11.3f}\n') + itp.write(f'{ids[0]:>6}{ids[1]:>6}{ids[2]:>6}{ids[3]:>6}{3:>6}{c[0]:>11.3f}') + itp.write(f'{c[1]:>11.3f}{c[2]:>11.3f}{c[3]:>11.3f}{c[4]:>11.3f}{c[5]:>11.3f}\n') # inversion dihedrals if len(terms['dihedral/inversion']) > 0: itp.write("; inversion dihedrals \n") - itp.write('; ai aj ak al f c0 c1 c2\n') + itp.write('; ai aj ak al f c0 c1 c2\n') for dihed in terms['dihedral/inversion']: ids = dihed.atomids + 1 + # itp.write(f'{ids[0]:>6}{ids[1]:>6}{ids[2]:>6}{ids[3]:>6}{1:>6}' + # f'{0.0:>11.3f}{dihed.fconst:>11.3f}{3:>11}\n') c0, c1, c2 = convert_to_inversion_rb(dihed.fconst, dihed.equ) - itp.write(f'{ids[0]:>6} {ids[1]:>6} {ids[2]:>6} {ids[3]:>6} {3:>6}' - f'{c0:>11.3f} {c1:>11.3f} {c2:>11.3f} {0:>11.1f} {0:>11.1f} {0:>11.1f}\n') + itp.write(f'{ids[0]:>6}{ids[1]:>6}{ids[2]:>6}{ids[3]:>6}{3:>6}' + f'{c0:>11.3f}{c1:>11.3f}{c2:>11.3f}{0:>11.1f}{0:>11.1f}{0:>11.1f}\n') def write_itp_exclusions(self, itp): if any(len(exclusion) > 0 for exclusion in self.exclusions): @@ -306,11 +536,20 @@ def write_itp_exclusions(self, itp): itp.write("\n") def make_pairs(self, neighbors, non_bonded): - polar_pairs = [] + pairs, polar_pairs = [], [] + + for pair in non_bonded.pairs: + pairs.append(sorted(pair)) if self.n_excl == 2: + for i in range(self.n_atoms): + for neigh in neighbors[2][i]: + if (i < neigh and [i, neigh] not in pairs and (i, neigh) + not in non_bonded.exclusions): + pairs.append([i, neigh]) + if self.polar: - for a1, a2 in non_bonded.pairs: + for a1, a2 in pairs: if a2 in non_bonded.alpha_map.keys(): polar_pairs.append([a1, non_bonded.alpha_map[a2]]) if a1 in non_bonded.alpha_map.keys(): @@ -318,15 +557,18 @@ def make_pairs(self, neighbors, non_bonded): if a1 in non_bonded.alpha_map.keys() and a2 in non_bonded.alpha_map.keys(): polar_pairs.append([non_bonded.alpha_map[a1], non_bonded.alpha_map[a2]]) - return non_bonded.pairs+polar_pairs + return pairs+polar_pairs def make_exclusions(self, non_bonded, neighbors, exclude_all): exclusions = [[] for _ in range(self.n_atoms)] - # input exclusions for exclusions if outside of n_excl - for a1, a2 in non_bonded.exclusions+non_bonded.pairs: - if all([a2 not in neighbors[i][a1] for i in range(self.n_excl+1)]): - exclusions[a1].append(a2+1) + # input exclusions + for exclusion in non_bonded.exclusions: + exclusions[exclusion[0]].append(exclusion[1]+1) + + # exclusions for input pairs + for pair in non_bonded.pairs: + exclusions[pair[0]].append(pair[1]+1) # fragment capping atom exclusions for i in exclude_all: @@ -375,11 +617,11 @@ def get_atom_names(self): def add_restraints(self, restraints, directory, fc=1000): with open(f"{directory}/{self.mol_name}_qforce{self.polar_title}.itp", "a") as itp: itp.write("[ dihedral_restraints ]\n") - itp.write("; ai aj ak al type phi dp kfac\n") + itp.write("; ai aj ak al type phi dp kfac\n") for restraint in restraints: a1, a2, a3, a4 = restraint[0]+1 phi = np.degrees(restraint[1]) - itp.write(f'{a1:>5} {a2:>5} {a3:>5} {a4:>5} {1:>5} {phi:>10.4f} 0.0 {fc}\n') + itp.write(f'{a1:>5}{a2:>5}{a3:>5}{a4:>5}{1:>5} {phi:>10.4f} 0.0 {fc}\n') def set_charge(self, non_bonded): q = np.copy(non_bonded.q) diff --git a/qforce/initialize.py b/qforce/initialize.py index 69acc89..5bea8c7 100644 --- a/qforce/initialize.py +++ b/qforce/initialize.py @@ -15,6 +15,10 @@ class Initialize(Colt): _user_input = """ [ff] +# MD software to use +# Currently amber supports only gaff2 atom types +forcefield = 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 @@ -23,7 +27,7 @@ class Initialize(Colt): n_excl = 2 :: int :: [2, 3] # Lennard jones method for the forcefield -lennard_jones = opls_auto :: str :: [gromos_auto, gromos, opls_auto, opls, gaff, gaff2, charmm36, ext] +lennard_jones = opls_auto :: str :: [gromos_auto, gromos, opls_auto, opls, gaff, gaff2, ext] # Use externally provided point charges in the file "ext_q" in the job directyory ext_charges = no :: bool @@ -164,7 +168,7 @@ def _check_and_copy_settings_file(job_dir, config_file): def initialize(filename, config_file, presets=None): - print(LOGO) + #print(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..b14fc33 100755 --- a/qforce/main.py +++ b/qforce/main.py @@ -56,7 +56,13 @@ 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) + + forcefield = config.ff.forcefield + if forcefield in ["gromacs"]: + ff.write_gromacs(job.dir, mol, qm_hessian_out.coords) + + if forcefield in ["amber"]: + ff.write_amber(job.dir, mol, qm_hessian_out.coords) print_outcome(job.dir) From f37c4b526e8f780875fac2d6c976617e37f52799 Mon Sep 17 00:00:00 2001 From: Mateus Zanotto <39527468+mateuszanotto@users.noreply.github.com> Date: Tue, 13 Sep 2022 13:12:04 -0300 Subject: [PATCH 06/28] Add files via upload --- bin/qforce | 51 ++++++++++++++++++++++++++++++++++++++++++++++++++- 1 file changed, 50 insertions(+), 1 deletion(-) diff --git a/bin/qforce b/bin/qforce index 81896cd..b98fc2d 100644 --- a/bin/qforce +++ b/bin/qforce @@ -1,7 +1,56 @@ #!/usr/bin/env python3 -from qforce.main import run +import sys +import threading +import time + +class ProgressBarThread(threading.Thread): + def __init__(self, label='Working', delay=0.1): + super(ProgressBarThread, self).__init__() + self.label = label + self.delay = delay # interval between updates + self.running = False + def start(self): + self.running = True + super(ProgressBarThread, self).start() + def run(self): + label = '\r' + self.label + ' ' + while self.running: + # for c in ('-', '\\', '|', '/'): + # for c in ('▌', '▀', '▐', '▄'): + for c in ('◐', '◓', '◑', '◒'): + # for c in ('▙', '▛', '▜', '▟'): + # for c in ('▤', '▧', '▥', '▨'): + + sys.stdout.write(label + c) + sys.stdout.flush() + time.sleep(self.delay) + def stop(self): + self.running = False + self.join() # wait for run() method to terminate + sys.stdout.write('\r' + len(self.label)*' ' + 2*' ' + '\r') # clean-up + sys.stdout.flush() + +def work(): + time.sleep(5) # +print(""" + ____ ______ + / __ \ | ____| + | | | |______| |__ ___ _ __ ___ ___ + | | | |______| __/ _ \| '__/ __/ _ \\ + | |__| | | | | (_) | | | (_| __/ + \___\_\ |_| \___/|_| \___\___| + + Selim Sami + University of Groningen - 2020 + ============================== +""") + +pb_thread = ProgressBarThread(' Initializing') +pb_thread.start() +from qforce.main import run +pb_thread.stop() if __name__ == '__main__': run() From d7f77e9051763905b8c86540ccbcc7ada0b1e5f0 Mon Sep 17 00:00:00 2001 From: Mateus Zanotto <39527468+mateuszanotto@users.noreply.github.com> Date: Tue, 13 Sep 2022 13:35:19 -0300 Subject: [PATCH 07/28] Added .itp (amber ff format) support --- qforce/forcefield.py | 99 ++++++++++++++++++++------------------------ 1 file changed, 45 insertions(+), 54 deletions(-) diff --git a/qforce/forcefield.py b/qforce/forcefield.py index a383ac6..f7a0185 100644 --- a/qforce/forcefield.py +++ b/qforce/forcefield.py @@ -271,9 +271,9 @@ def write_top(self, directory): with open(f"{directory}/gas{self.polar_title}.top", "w") as top: # defaults top.write("\n[ defaults ]\n") - top.write(";nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ\n") - top.write(f"{1:>7}{self.comb_rule:>12}{'yes':>12}{self.fudge_lj:>10}{self.fudge_q:>10}" - "\n\n\n") + top.write("; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ\n") + top.write(f"{1:>8} {self.comb_rule:>12} {'yes':>12} {self.fudge_lj:>12} " + f"{self.fudge_q:>12}\n\n\n") top.write("; Include the molecule ITP\n") top.write(f'#include "./{self.mol_name}_qforce{self.polar_title}.itp"\n\n\n') @@ -397,24 +397,24 @@ def write_itp_atoms_and_molecule(self, itp, non_bonded): # atoms itp.write("\n[ atoms ]\n") - itp.write("; nr type resnr resnm atom cgrp charge mass\n") + itp.write("; nr type resnr resnm atom cgrp charge mass\n") for i, (lj_type, a_name, q, mass) in enumerate(zip(non_bonded.lj_types, self.atom_names, self.q, self.masses), start=1): - itp.write(f'{i:>5}{lj_type:>9}{1:>6}{self.residue:>6}{a_name:>7}{i:>5}{q:>11.5f}') - itp.write(f'{mass:>10.5f}\n') + itp.write(f'{i:>5} {lj_type:>9} {1:>6} {self.residue:>6} {a_name:>7} {i:>5} ') + itp.write(f'{q:>11.5f} {mass:>10.5f}\n') if self.polar: for i, (atom, drude) in enumerate(non_bonded.alpha_map.items(), start=1): - itp.write(f'{drude+1:>5}{"DP":>9}{2:>6}{"DRU":>6}{f"D{i}":>7}{atom+1:>5}') - itp.write(f'{-8.:>11.5f}{0.:>10.5f}\n') + itp.write(f'{drude+1:>5} {"DP":>9} {2:>6} {"DRU":>6} {f"D{i}":>7} {atom+1:>5}') + itp.write(f'{-8.:>11.5f} {0.:>10.5f}\n') def write_itp_polarization(self, itp, non_bonded): # polarization itp.write("\n[ polarization ]\n") - itp.write("; i j f alpha\n") + itp.write("; i j f alpha\n") for atom, drude in non_bonded.alpha_map.items(): alpha = non_bonded.alpha[atom]*1e-3 - itp.write(f"{atom+1:>6}{drude+1:>6}{1:>6}{alpha:>14.8f}\n") + itp.write(f"{atom+1:>6} {drude+1:>6} {1:>6} {alpha:>14.8f}\n") # # thole polarization # if self.thole != []: @@ -427,32 +427,35 @@ def write_itp_polarization(self, itp, non_bonded): def write_itp_pairs(self, itp): if self.pairs != []: itp.write("\n[ pairs ]\n") - itp.write("; ai aj func \n") + itp.write("; ai aj func\n") for pair in self.pairs: - itp.write(f"{pair[0]+1:>6}{pair[1]+1:>6}{1:>6}\n") + itp.write(f"{pair[0]+1:>6} {pair[1]+1:>6} {1:>6}\n") def write_itp_bonds(self, itp, terms, alpha_map): itp.write("\n[ bonds ]\n") - itp.write("; ai aj f r0 kb\n") + itp.write("; ai aj f r0 k_bond\n") for bond in terms['bond']: ids = bond.atomids + 1 equ = bond.equ * 0.1 fconst = bond.fconst * 100 - itp.write(f'{ids[0]:>6}{ids[1]:>6}{1:>6}{equ:>10.5f}{fconst:>10.0f}\n') + itp.write(f'{ids[0]:>6} {ids[1]:>6} {1:>6} {equ:>10.5f} {fconst:>10.0f}\n') if self.polar: - itp.write('; ai aj f - polar connections\n') + itp.write('; ai aj f - polar connections\n') for bond in terms['bond']: a1, a2 = bond.atomids if a2 in alpha_map.keys(): - itp.write(f'{a1+1:>6}{alpha_map[a2]+1:>6}{5:>6}\n') + itp.write(f'{a1+1:>6} {alpha_map[a2]+1:>6} {5:>6}\n') if a1 in alpha_map.keys(): - itp.write(f'{a2+1:>6}{alpha_map[a1]+1:>6}{5:>6}\n') + itp.write(f'{a2+1:>6} {alpha_map[a1]+1:>6} {5:>6}\n') def write_itp_angles(self, itp, terms): itp.write("\n[ angles ]\n") - itp.write("; ai aj ak f th0 kth\n") + itp.write("; ai aj ak f theta0 k_theta") + if self.urey: + itp.write(' r0 k_bond') + itp.write('\n') for angle in terms['angle']: ids = angle.atomids + 1 equ = np.degrees(angle.equ) @@ -462,12 +465,13 @@ def write_itp_angles(self, itp, terms): urey = [term for term in terms['urey'] if np.array_equal(term.atomids, angle.atomids)] if not self.urey or len(urey) == 0: - itp.write(f'{ids[0]:>6}{ids[1]:>6}{ids[2]:>6}{1:>6}{equ:>11.3f}{fconst:>13.3f}\n') + itp.write(f'{ids[0]:>6} {ids[1]:>6} {ids[2]:>6} {1:>6} {equ:>10.3f} ' + f'{fconst:>10.3f}\n') else: urey_equ = urey[0].equ * 0.1 urey_fconst = urey[0].fconst * 100 - itp.write(f'{ids[0]:>6}{ids[1]:>6}{ids[2]:>6}{5:>6}{equ:>11.3f}{fconst:>13.3f}' - f'{urey_equ:>10.5f}{urey_fconst:>13.3f}\n') + itp.write(f'{ids[0]:>6} {ids[1]:>6} {ids[2]:>6} {5:>6} {equ:>10.3f} ' + f'{fconst:>10.3f} {urey_equ:>10.5f} {urey_fconst:>10.1f}\n') def write_itp_dihedrals(self, itp, terms): if len(terms['dihedral']) > 0: @@ -476,54 +480,52 @@ def write_itp_dihedrals(self, itp, terms): # rigid dihedrals if len(terms['dihedral/rigid']) > 0: itp.write("; rigid dihedrals \n") - itp.write("; ai aj ak al f th0 kth\n") + itp.write("; ai aj ak al f theta0 k_theta\n") for dihed in terms['dihedral/rigid']: ids = dihed.atomids + 1 equ = np.degrees(dihed.equ) fconst = dihed.fconst - itp.write(f'{ids[0]:>6}{ids[1]:>6}{ids[2]:>6}{ids[3]:>6}{2:>6}{equ:>11.3f}') + itp.write(f'{ids[0]:>6} {ids[1]:>6} {ids[2]:>6} {ids[3]:>6} {2:>6} {equ:>11.3f} ') itp.write(f'{fconst:>13.3f}\n') # improper dihedrals if len(terms['dihedral/improper']) > 0: itp.write("; improper dihedrals \n") - itp.write("; ai aj ak al f th0 kth\n") + itp.write("; ai aj ak al f theta0 k_theta\n") for dihed in terms['dihedral/improper']: ids = dihed.atomids + 1 equ = np.degrees(dihed.equ) fconst = dihed.fconst - itp.write(f'{ids[0]:>6}{ids[1]:>6}{ids[2]:>6}{ids[3]:>6}{2:>6}{equ:>11.3f}') + itp.write(f'{ids[0]:>6} {ids[1]:>6} {ids[2]:>6} {ids[3]:>6} {2:>6} {equ:>11.3f} ') itp.write(f'{fconst:>13.3f}\n') # flexible dihedrals if len(terms['dihedral/flexible']) > 0: itp.write("; flexible dihedrals \n") - itp.write('; ai aj ak al f c0 c1 c2 c3') - itp.write(' c4 c5\n') + itp.write('; ai aj ak al f c0 c1 c2') + itp.write(' c3 c4 c5\n') for dihed in terms['dihedral/flexible']: ids = dihed.atomids + 1 c = dihed.equ - itp.write(f'{ids[0]:>6}{ids[1]:>6}{ids[2]:>6}{ids[3]:>6}{3:>6}{c[0]:>11.3f}') - itp.write(f'{c[1]:>11.3f}{c[2]:>11.3f}{c[3]:>11.3f}{c[4]:>11.3f}{c[5]:>11.3f}\n') + itp.write(f'{ids[0]:>6} {ids[1]:>6} {ids[2]:>6} {ids[3]:>6} {3:>6} {c[0]:>11.3f} ') + itp.write(f'{c[1]:>11.3f} {c[2]:>11.3f} {c[3]:>11.3f} {c[4]:>11.3f} {c[5]:>11.3f}\n') # inversion dihedrals if len(terms['dihedral/inversion']) > 0: itp.write("; inversion dihedrals \n") - itp.write('; ai aj ak al f c0 c1 c2\n') + itp.write('; ai aj ak al f c0 c1 c2\n') for dihed in terms['dihedral/inversion']: ids = dihed.atomids + 1 - # itp.write(f'{ids[0]:>6}{ids[1]:>6}{ids[2]:>6}{ids[3]:>6}{1:>6}' - # f'{0.0:>11.3f}{dihed.fconst:>11.3f}{3:>11}\n') c0, c1, c2 = convert_to_inversion_rb(dihed.fconst, dihed.equ) - itp.write(f'{ids[0]:>6}{ids[1]:>6}{ids[2]:>6}{ids[3]:>6}{3:>6}' - f'{c0:>11.3f}{c1:>11.3f}{c2:>11.3f}{0:>11.1f}{0:>11.1f}{0:>11.1f}\n') + itp.write(f'{ids[0]:>6} {ids[1]:>6} {ids[2]:>6} {ids[3]:>6} {3:>6}' + f'{c0:>11.3f} {c1:>11.3f} {c2:>11.3f} {0:>11.1f} {0:>11.1f} {0:>11.1f}\n') def write_itp_exclusions(self, itp): if any(len(exclusion) > 0 for exclusion in self.exclusions): @@ -536,20 +538,12 @@ def write_itp_exclusions(self, itp): itp.write("\n") def make_pairs(self, neighbors, non_bonded): - pairs, polar_pairs = [], [] - - for pair in non_bonded.pairs: - pairs.append(sorted(pair)) + polar_pairs = [] if self.n_excl == 2: - for i in range(self.n_atoms): - for neigh in neighbors[2][i]: - if (i < neigh and [i, neigh] not in pairs and (i, neigh) - not in non_bonded.exclusions): - pairs.append([i, neigh]) if self.polar: - for a1, a2 in pairs: + for a1, a2 in non_bonded.pairs: if a2 in non_bonded.alpha_map.keys(): polar_pairs.append([a1, non_bonded.alpha_map[a2]]) if a1 in non_bonded.alpha_map.keys(): @@ -557,18 +551,15 @@ def make_pairs(self, neighbors, non_bonded): if a1 in non_bonded.alpha_map.keys() and a2 in non_bonded.alpha_map.keys(): polar_pairs.append([non_bonded.alpha_map[a1], non_bonded.alpha_map[a2]]) - return pairs+polar_pairs + return non_bonded.pairs+polar_pairs def make_exclusions(self, non_bonded, neighbors, exclude_all): exclusions = [[] for _ in range(self.n_atoms)] - # input exclusions - for exclusion in non_bonded.exclusions: - exclusions[exclusion[0]].append(exclusion[1]+1) - - # exclusions for input pairs - for pair in non_bonded.pairs: - exclusions[pair[0]].append(pair[1]+1) + # input exclusions for exclusions if outside of n_excl + for a1, a2 in non_bonded.exclusions+non_bonded.pairs: + if all([a2 not in neighbors[i][a1] for i in range(self.n_excl+1)]): + exclusions[a1].append(a2+1) # fragment capping atom exclusions for i in exclude_all: @@ -617,11 +608,11 @@ def get_atom_names(self): def add_restraints(self, restraints, directory, fc=1000): with open(f"{directory}/{self.mol_name}_qforce{self.polar_title}.itp", "a") as itp: itp.write("[ dihedral_restraints ]\n") - itp.write("; ai aj ak al type phi dp kfac\n") + itp.write("; ai aj ak al type phi dp kfac\n") for restraint in restraints: a1, a2, a3, a4 = restraint[0]+1 phi = np.degrees(restraint[1]) - itp.write(f'{a1:>5}{a2:>5}{a3:>5}{a4:>5}{1:>5} {phi:>10.4f} 0.0 {fc}\n') + itp.write(f'{a1:>5} {a2:>5} {a3:>5} {a4:>5} {1:>5} {phi:>10.4f} 0.0 {fc}\n') def set_charge(self, non_bonded): q = np.copy(non_bonded.q) From ba2669d559d5fb91bab02a563049bc5cf2ad2ccb Mon Sep 17 00:00:00 2001 From: Mateus Zanotto <39527468+mateuszanotto@users.noreply.github.com> Date: Tue, 13 Sep 2022 13:48:06 -0300 Subject: [PATCH 08/28] Add a progress bar while import modules Add a progress bar while import modules for aesthetical. In the tests I did, it didn't added any computational cost. From 1acfc9cddca25243b4d9eede1ef5c59615d5f03a Mon Sep 17 00:00:00 2001 From: Mateus Zanotto <39527468+mateuszanotto@users.noreply.github.com> Date: Tue, 13 Sep 2022 14:04:43 -0300 Subject: [PATCH 09/28] Added amber ff option Also removed the logo print as it is printed now in the bin/qforce. --- qforce/initialize.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/qforce/initialize.py b/qforce/initialize.py index 5bea8c7..bfe0443 100644 --- a/qforce/initialize.py +++ b/qforce/initialize.py @@ -27,7 +27,7 @@ class Initialize(Colt): n_excl = 2 :: int :: [2, 3] # Lennard jones method for the forcefield -lennard_jones = opls_auto :: str :: [gromos_auto, gromos, opls_auto, opls, gaff, gaff2, ext] +lennard_jones = opls_auto :: str :: [gromos_auto, gromos, opls_auto, opls, gaff, gaff2, charmm36, ext] # Use externally provided point charges in the file "ext_q" in the job directyory ext_charges = no :: bool From 1ea47a4021221d104a2de57a6263c8ddf1f7a0cb Mon Sep 17 00:00:00 2001 From: Mateus Zanotto <39527468+mateuszanotto@users.noreply.github.com> Date: Tue, 13 Sep 2022 14:05:12 -0300 Subject: [PATCH 10/28] Added .itp (amber format) support --- qforce/forcefield.py | 486 +++++++++++++++++++++---------------------- 1 file changed, 243 insertions(+), 243 deletions(-) diff --git a/qforce/forcefield.py b/qforce/forcefield.py index f7a0185..3949c1e 100644 --- a/qforce/forcefield.py +++ b/qforce/forcefield.py @@ -32,236 +32,6 @@ def __init__(self, job_name, config, mol, neighbors, exclude_all=[]): else: self.polar_title = '' -## Amber - def write_amber(self, directory, mol, coords): - atom_ids, unique_at =self.get_atom_types(mol.topo, mol.non_bonded) - self.write_mol2(directory, mol, coords, atom_ids, unique_at) - self.write_frcmod(directory, mol, coords, atom_ids, unique_at) - - def write_mol2(self, directory, mol, coords, atom_ids, unique_at): - 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, atom_ids, unique_at) - self.write_mol2_bond(mol2, mol.topo, mol.terms, atom_ids, unique_at) - - def write_frcmod(self, directory, mol, coords, atom_ids, unique_at): - with open(f"{directory}/{self.mol_name}_qforce{self.polar_title}.frcmod", "w") as frcmod: - self.write_frcmod_mass(frcmod, mol.non_bonded, atom_ids, unique_at) - self.write_frcmod_bonds(frcmod, mol.terms, mol.non_bonded.alpha_map, atom_ids, unique_at) - self.write_frcmod_angles(frcmod, mol.terms, atom_ids, unique_at) - self.write_frcmod_dihedrals(frcmod, mol.terms, atom_ids, unique_at) - self.write_frcmod_nonbond(frcmod, mol.non_bonded, atom_ids, unique_at) - - -# TLeap is a helper program from AMBER to create the topologies. -# Would be good for the final user to create automatically a -# Tleap script with the new atom types. - - # def write_tleap_script(self, directory): - # with open(f"{directory}/tleap_script_{self.mol_name}.in", "w") as tleap: - # tleap.write(f""" - # # To run the script use tleap on amberTools and the following command: - # # $ tleap -f tleap_script_{self.mol_name}.in > tleap_script_{self.mol_name}.o - # - # source leaprc.protein.ff19SB\nsource leaprc.water.opc - # source leaprc.gaff\n# You can source other force fields here if necessary - # loadmol2 {self.mol_name}_qforce{self.polar_title}.mol2 - # - # loadamberparams {self.mol_name}_qforce{self.polar_title}.frcmod - # - # ## Add any other command, for example "solvateoct MOL TIP3PBOX 10.0" - # - # MOL = loadPDB 1lyd.solv.pdb - # saveamberparm MOL {self.mol_name}_qforce{self.polar_title}.prmtop - # saveamberparm MOL {self.mol_name}_qforce{self.polar_title}.inpcrd - # savepdb MOL {self.mol_name}_qforce{self.polar_title}.pdb """) - -############################################## -# MOL2 WRITING # -############################################## - - def write_mol2_title(self, mol2): - mol2.write(LOGO_HASH) - 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(f"esp")## type of charge - mol2.write("\n\n") - - def write_mol2_atom(self, mol2, topo, coords, non_bonded, atom_ids, unique_at): - mol2.write(f"@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} {coords[i_idx-1][1]:10.4f} {coords[i_idx-1][2]:10.4f}") - mol2.write(f" {unique_at[i_idx][atom_ids[i_idx]]} {1} {self.mol_name} {q:10.6f}\n") - - def write_mol2_bond(self, mol2, topo, terms, atom_ids, unique_at): - n_bonds = 1 - mol2.write(f"@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(f"@SUBSTRUCTURE\n") - mol2.write(f"{1:5} {self.mol_name} {1:8} TEMP {0:>8d} **** **** {0:5} ROOT\n") - -############################################## -# FRCMOD WRITING # -############################################## - - def write_frcmod_mass(self, frcmod, non_bonded, atom_ids, unique_at): - frcmod.write(f'{self.mol_name} - frcmod generated by QForce\n') - frcmod.write(f'MASS\n') - for i, mass in enumerate((self.masses), start=1): - frcmod.write(f"{unique_at[i][atom_ids[i]]} {mass}\n") - - def write_frcmod_bonds(self, frcmod, terms, alpha_map, atom_ids, unique_at): - frcmod.write("\nBOND\n") - for bond in terms['bond']: - ids = bond.atomids + 1 - fconst = bond.fconst * (0.239005)/2 # kJ/mol/A^2 -> kcal/mol/A^2 - equ = bond.equ # Ang - frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") - frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}") - frcmod.write(f"{fconst:>10.2f}{equ:>10.2f}\n") - - def write_frcmod_angles(self, frcmod, terms, atom_ids, unique_at): - frcmod.write("\nANGLE\n") - for angle in terms['angle']: - ids = angle.atomids + 1 - fconst = angle.fconst * (0.239005)/2 # kJ/mol/rad^2 -> kcal/mol/rad^2 - equ = np.degrees(angle.equ) # Degrees -> Degrees - - 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"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") - frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") - frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}") - frcmod.write(f"{fconst:>10.2f}{equ:>10.2f}\n") - else: - urey_equ = urey[0].equ - urey_fconst = urey[0].fconst * (0.239005)/2 # kJ/mol/rad^2 -> kcal/mol/rad^2 - frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") - frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") - frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}") - frcmod.write(f"{fconst:>10.2f}{equ:>10.2f}\n") - - def write_frcmod_dihedrals(self, frcmod, terms, atom_ids, unique_at): - 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 * (0.239005) # kJ/mol -> kcal/mol #scale cos - - fconstAmb = (2*fconst/(2**2)) # Vn = (2*fcte/(n**2)) - frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") #Atom# 1 - frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") #Atom 2 - frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}-") #Atom 3 - frcmod.write(f"{unique_at[ids[3]][atom_ids[ids[3]]]:<2}") #Atom 4 - frcmod.write(f" 1 {fconstAmb:>15.2f} {equ:>15.2f} 2\n") #IDIVF=1, fcte, angle, n=2 (n*angle) - - if len(terms['dihedral/flexible']) > 0: - - for dihed in terms['dihedral/flexible']: - ids = dihed.atomids + 1 - c = dihed.equ - - for n in range(3, 0, -1): - if n%2==1: - equ = 180.0 - else: - equ = 0.00 - frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") #Atom 1 - frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") #Atom 2 - frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}-") #Atom 3 - frcmod.write(f"{unique_at[ids[3]][atom_ids[ids[3]]]:<2}") #Atom 4 - frcmod.write(f" 1 {c[n]:>15.2f} {equ:>15.2f} -{n+1}\n") #IDIVF=1, fcte, angle, n=1 (n*angle) - - # Last term without -n - frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") #Atom 1 - frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") #Atom 2 - frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}-") #Atom 3 - frcmod.write(f"{unique_at[ids[3]][atom_ids[ids[3]]]:<2}") #Atom 4 - frcmod.write(f" 1 {c[0]:>15.2f} {0.00:>15.2f} 1\n") #IDIVF=1, fcte, angle, n=1 (n*angle) - - - # 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) # Degrees -> Degrees - fconst = dihed.fconst * (0.239005)/2 # kJ/mol -> kcal/mol -> f/2 as is Vn/2 in amber - frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") #Atom 1 - frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") #Atom 2 - frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}-") #Atom 3 - frcmod.write(f"{unique_at[ids[3]][atom_ids[ids[3]]]:<2}") #Atom 4 - frcmod.write(f" 1 {fconst:>15.2f} {equ:>15.2f} 2\n") #IDIVF=1, fcte, angle, n=2 (n*angle) - - def write_frcmod_nonbond(self, frcmod, non_bonded, atom_ids, unique_at): - gro_atomtypes, gro_nonbonded, gro_1_4 = self.convert_to_gromacs_nonbonded(non_bonded) - 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"{unique_at[i][atom_ids[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") - - -############################################## -# GET ATOM TYPES AND IDS # -############################################## - def get_atom_types(self, topo, non_bonded): - atom_ids={} #dictonary containing original atom types - unique_at={} #dictonary containing new unique atom types - unique_masses={} - - ascii_lowercase = list(string.ascii_lowercase) - ascii_digits = list(string.digits) - ascii_uppercase = list(string.ascii_uppercase) - - for i, (lj_type, a_name) in enumerate(zip(non_bonded.lj_types, self.atom_names), start=1): - atom_ids[i]=lj_type - if self.n_atoms < 36: - if i <= 10: - unique_at[i]={atom_ids[i] : "Q{}".format(ascii_digits[i-1])} - elif i > 10: - unique_at[i]={atom_ids[i] : "Q{}".format(ascii_lowercase[i-11])} - else: - unique_at[i]={atom_ids[i] : "{}{}".format(ascii_uppercase[i-1], ascii_lowercase[i-1])} - - ## DEBUG : check correspondance between atom_ids and unique_at - #for i, lj_type in enumerate(zip(non_bonded.lj_types), start=1): - # print(unique_at) - # for i in range(1,self.n_atoms+1,1): - # print (i, atom_ids[i], unique_at[i][atom_ids[i]]) - # print ("------------------") - - return atom_ids, unique_at - -############################################## - def write_gromacs(self, directory, mol, coords): self.write_itp(mol, directory) self.write_top(directory) @@ -271,9 +41,9 @@ def write_top(self, directory): with open(f"{directory}/gas{self.polar_title}.top", "w") as top: # defaults top.write("\n[ defaults ]\n") - top.write("; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ\n") + top.write("; nbfunc comb-rule gen-pairs fudgeLJ fudgeQQ\n") top.write(f"{1:>8} {self.comb_rule:>12} {'yes':>12} {self.fudge_lj:>12} " - f"{self.fudge_q:>12}\n\n\n") + f"{self.fudge_q:>12}\n\n\n") top.write("; Include the molecule ITP\n") top.write(f'#include "./{self.mol_name}_qforce{self.polar_title}.itp"\n\n\n') @@ -397,7 +167,7 @@ def write_itp_atoms_and_molecule(self, itp, non_bonded): # atoms itp.write("\n[ atoms ]\n") - itp.write("; nr type resnr resnm atom cgrp charge mass\n") + itp.write("; nr type resnr resnm atom cgrp charge mass\n") for i, (lj_type, a_name, q, mass) in enumerate(zip(non_bonded.lj_types, self.atom_names, self.q, self.masses), start=1): itp.write(f'{i:>5} {lj_type:>9} {1:>6} {self.residue:>6} {a_name:>7} {i:>5} ') @@ -411,7 +181,7 @@ def write_itp_atoms_and_molecule(self, itp, non_bonded): def write_itp_polarization(self, itp, non_bonded): # polarization itp.write("\n[ polarization ]\n") - itp.write("; i j f alpha\n") + itp.write("; i j f alpha\n") for atom, drude in non_bonded.alpha_map.items(): alpha = non_bonded.alpha[atom]*1e-3 itp.write(f"{atom+1:>6} {drude+1:>6} {1:>6} {alpha:>14.8f}\n") @@ -427,13 +197,13 @@ def write_itp_polarization(self, itp, non_bonded): def write_itp_pairs(self, itp): if self.pairs != []: itp.write("\n[ pairs ]\n") - itp.write("; ai aj func\n") + itp.write("; ai aj func\n") for pair in self.pairs: itp.write(f"{pair[0]+1:>6} {pair[1]+1:>6} {1:>6}\n") def write_itp_bonds(self, itp, terms, alpha_map): itp.write("\n[ bonds ]\n") - itp.write("; ai aj f r0 k_bond\n") + itp.write("; ai aj f r0 k_bond\n") for bond in terms['bond']: ids = bond.atomids + 1 equ = bond.equ * 0.1 @@ -441,7 +211,7 @@ def write_itp_bonds(self, itp, terms, alpha_map): itp.write(f'{ids[0]:>6} {ids[1]:>6} {1:>6} {equ:>10.5f} {fconst:>10.0f}\n') if self.polar: - itp.write('; ai aj f - polar connections\n') + itp.write('; ai aj f - polar connections\n') for bond in terms['bond']: a1, a2 = bond.atomids @@ -452,7 +222,7 @@ def write_itp_bonds(self, itp, terms, alpha_map): def write_itp_angles(self, itp, terms): itp.write("\n[ angles ]\n") - itp.write("; ai aj ak f theta0 k_theta") + itp.write("; ai aj ak f theta0 k_theta") if self.urey: itp.write(' r0 k_bond') itp.write('\n') @@ -480,7 +250,7 @@ def write_itp_dihedrals(self, itp, terms): # rigid dihedrals if len(terms['dihedral/rigid']) > 0: itp.write("; rigid dihedrals \n") - itp.write("; ai aj ak al f theta0 k_theta\n") + itp.write("; ai aj ak al f theta0 k_theta\n") for dihed in terms['dihedral/rigid']: ids = dihed.atomids + 1 @@ -493,7 +263,7 @@ def write_itp_dihedrals(self, itp, terms): # improper dihedrals if len(terms['dihedral/improper']) > 0: itp.write("; improper dihedrals \n") - itp.write("; ai aj ak al f theta0 k_theta\n") + itp.write("; ai aj ak al f theta0 k_theta\n") for dihed in terms['dihedral/improper']: ids = dihed.atomids + 1 @@ -519,7 +289,7 @@ def write_itp_dihedrals(self, itp, terms): # inversion dihedrals if len(terms['dihedral/inversion']) > 0: itp.write("; inversion dihedrals \n") - itp.write('; ai aj ak al f c0 c1 c2\n') + itp.write('; ai aj ak al f c0 c1 c2\n') for dihed in terms['dihedral/inversion']: ids = dihed.atomids + 1 @@ -608,11 +378,11 @@ def get_atom_names(self): def add_restraints(self, restraints, directory, fc=1000): with open(f"{directory}/{self.mol_name}_qforce{self.polar_title}.itp", "a") as itp: itp.write("[ dihedral_restraints ]\n") - itp.write("; ai aj ak al type phi dp kfac\n") + itp.write("; ai aj ak al type phi dp kfac\n") for restraint in restraints: a1, a2, a3, a4 = restraint[0]+1 phi = np.degrees(restraint[1]) - itp.write(f'{a1:>5} {a2:>5} {a3:>5} {a4:>5} {1:>5} {phi:>10.4f} 0.0 {fc}\n') + itp.write(f'{a1:>5} {a2:>5} {a3:>5} {a4:>5} {1:>5} {phi:>10.4f} 0.0 {fc}\n') def set_charge(self, non_bonded): q = np.copy(non_bonded.q) @@ -620,6 +390,236 @@ def set_charge(self, non_bonded): q[list(non_bonded.alpha_map.keys())] += 8 return q +## Amber + def write_amber(self, directory, mol, coords): + atom_ids, unique_at =self.get_atom_types(mol.topo, mol.non_bonded) + self.write_mol2(directory, mol, coords, atom_ids, unique_at) + self.write_frcmod(directory, mol, coords, atom_ids, unique_at) + + def write_mol2(self, directory, mol, coords, atom_ids, unique_at): + 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, atom_ids, unique_at) + self.write_mol2_bond(mol2, mol.topo, mol.terms, atom_ids, unique_at) + + def write_frcmod(self, directory, mol, coords, atom_ids, unique_at): + with open(f"{directory}/{self.mol_name}_qforce{self.polar_title}.frcmod", "w") as frcmod: + self.write_frcmod_mass(frcmod, mol.non_bonded, atom_ids, unique_at) + self.write_frcmod_bonds(frcmod, mol.terms, mol.non_bonded.alpha_map, atom_ids, unique_at) + self.write_frcmod_angles(frcmod, mol.terms, atom_ids, unique_at) + self.write_frcmod_dihedrals(frcmod, mol.terms, atom_ids, unique_at) + self.write_frcmod_nonbond(frcmod, mol.non_bonded, atom_ids, unique_at) + + +# TLeap is a helper program from AMBER to create the topologies. +# Would be good for the final user to create automatically a +# Tleap script with the new atom types. + + # def write_tleap_script(self, directory): + # with open(f"{directory}/tleap_script_{self.mol_name}.in", "w") as tleap: + # tleap.write(f""" + # # To run the script use tleap on amberTools and the following command: + # # $ tleap -f tleap_script_{self.mol_name}.in > tleap_script_{self.mol_name}.o + # + # source leaprc.protein.ff19SB\nsource leaprc.water.opc + # source leaprc.gaff\n# You can source other force fields here if necessary + # loadmol2 {self.mol_name}_qforce{self.polar_title}.mol2 + # + # loadamberparams {self.mol_name}_qforce{self.polar_title}.frcmod + # + # ## Add any other command, for example "solvateoct MOL TIP3PBOX 10.0" + # + # MOL = loadPDB 1lyd.solv.pdb + # saveamberparm MOL {self.mol_name}_qforce{self.polar_title}.prmtop + # saveamberparm MOL {self.mol_name}_qforce{self.polar_title}.inpcrd + # savepdb MOL {self.mol_name}_qforce{self.polar_title}.pdb """) + +############################################## +# MOL2 WRITING # +############################################## + + def write_mol2_title(self, mol2): + mol2.write(LOGO_HASH) + 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(f"esp")## type of charge + mol2.write("\n\n") + + def write_mol2_atom(self, mol2, topo, coords, non_bonded, atom_ids, unique_at): + mol2.write(f"@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} {coords[i_idx-1][1]:10.4f} {coords[i_idx-1][2]:10.4f}") + mol2.write(f" {unique_at[i_idx][atom_ids[i_idx]]} {1} {self.mol_name} {q:10.6f}\n") + + def write_mol2_bond(self, mol2, topo, terms, atom_ids, unique_at): + n_bonds = 1 + mol2.write(f"@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(f"@SUBSTRUCTURE\n") + mol2.write(f"{1:5} {self.mol_name} {1:8} TEMP {0:>8d} **** **** {0:5} ROOT\n") + +############################################## +# FRCMOD WRITING # +############################################## + + def write_frcmod_mass(self, frcmod, non_bonded, atom_ids, unique_at): + frcmod.write(f'{self.mol_name} - frcmod generated by QForce\n') + frcmod.write(f'MASS\n') + for i, mass in enumerate((self.masses), start=1): + frcmod.write(f"{unique_at[i][atom_ids[i]]} {mass}\n") + + def write_frcmod_bonds(self, frcmod, terms, alpha_map, atom_ids, unique_at): + frcmod.write("\nBOND\n") + for bond in terms['bond']: + ids = bond.atomids + 1 + fconst = bond.fconst * (0.239005)/2 # kJ/mol/A^2 -> kcal/mol/A^2 + equ = bond.equ # Ang + frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") + frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}") + frcmod.write(f"{fconst:>10.2f}{equ:>10.2f}\n") + + def write_frcmod_angles(self, frcmod, terms, atom_ids, unique_at): + frcmod.write("\nANGLE\n") + for angle in terms['angle']: + ids = angle.atomids + 1 + fconst = angle.fconst * (0.239005)/2 # kJ/mol/rad^2 -> kcal/mol/rad^2 + equ = np.degrees(angle.equ) # Degrees -> Degrees + + 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"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") + frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") + frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}") + frcmod.write(f"{fconst:>10.2f}{equ:>10.2f}\n") + else: + urey_equ = urey[0].equ + urey_fconst = urey[0].fconst * (0.239005)/2 # kJ/mol/rad^2 -> kcal/mol/rad^2 + frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") + frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") + frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}") + frcmod.write(f"{fconst:>10.2f}{equ:>10.2f}\n") + + def write_frcmod_dihedrals(self, frcmod, terms, atom_ids, unique_at): + 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 * (0.239005) # kJ/mol -> kcal/mol #scale cos + + fconstAmb = (2*fconst/(2**2)) # Vn = (2*fcte/(n**2)) + frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") #Atom# 1 + frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") #Atom 2 + frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}-") #Atom 3 + frcmod.write(f"{unique_at[ids[3]][atom_ids[ids[3]]]:<2}") #Atom 4 + frcmod.write(f" 1 {fconstAmb:>15.2f} {equ:>15.2f} 2\n") #IDIVF=1, fcte, angle, n=2 (n*angle) + + if len(terms['dihedral/flexible']) > 0: + + for dihed in terms['dihedral/flexible']: + ids = dihed.atomids + 1 + c = dihed.equ + + for n in range(3, 0, -1): + if n%2==1: + equ = 180.0 + else: + equ = 0.00 + frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") #Atom 1 + frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") #Atom 2 + frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}-") #Atom 3 + frcmod.write(f"{unique_at[ids[3]][atom_ids[ids[3]]]:<2}") #Atom 4 + frcmod.write(f" 1 {c[n]:>15.2f} {equ:>15.2f} -{n+1}\n") #IDIVF=1, fcte, angle, n=1 (n*angle) + + # Last term without -n + frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") #Atom 1 + frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") #Atom 2 + frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}-") #Atom 3 + frcmod.write(f"{unique_at[ids[3]][atom_ids[ids[3]]]:<2}") #Atom 4 + frcmod.write(f" 1 {c[0]:>15.2f} {0.00:>15.2f} 1\n") #IDIVF=1, fcte, angle, n=1 (n*angle) + + + # 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) # Degrees -> Degrees + fconst = dihed.fconst * (0.239005)/2 # kJ/mol -> kcal/mol -> f/2 as is Vn/2 in amber + frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") #Atom 1 + frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") #Atom 2 + frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}-") #Atom 3 + frcmod.write(f"{unique_at[ids[3]][atom_ids[ids[3]]]:<2}") #Atom 4 + frcmod.write(f" 1 {fconst:>15.2f} {equ:>15.2f} 2\n") #IDIVF=1, fcte, angle, n=2 (n*angle) + + def write_frcmod_nonbond(self, frcmod, non_bonded, atom_ids, unique_at): + gro_atomtypes, gro_nonbonded, gro_1_4 = self.convert_to_gromacs_nonbonded(non_bonded) + 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"{unique_at[i][atom_ids[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") + + +############################################## +# GET ATOM TYPES AND IDS # +############################################## + def get_atom_types(self, topo, non_bonded): + atom_ids={} #dictonary containing original atom types + unique_at={} #dictonary containing new unique atom types + unique_masses={} + + ascii_lowercase = list(string.ascii_lowercase) + ascii_digits = list(string.digits) + ascii_uppercase = list(string.ascii_uppercase) + + for i, (lj_type, a_name) in enumerate(zip(non_bonded.lj_types, self.atom_names), start=1): + atom_ids[i]=lj_type + if self.n_atoms < 36: + if i <= 10: + unique_at[i]={atom_ids[i] : "Q{}".format(ascii_digits[i-1])} + elif i > 10: + unique_at[i]={atom_ids[i] : "Q{}".format(ascii_lowercase[i-11])} + else: + unique_at[i]={atom_ids[i] : "{}{}".format(ascii_uppercase[i-1], ascii_lowercase[i-1])} + + ## DEBUG : check correspondance between atom_ids and unique_at + #for i, lj_type in enumerate(zip(non_bonded.lj_types), start=1): + # print(unique_at) + # for i in range(1,self.n_atoms+1,1): + # print (i, atom_ids[i], unique_at[i][atom_ids[i]]) + # print ("------------------") + + return atom_ids, unique_at + +############################################## + # bohr2nm = 0.052917721067 # if polar: # alphas = qm.alpha*bohr2nm**3 From d66e274ffab37a015012a66f6c287bce8093b27e Mon Sep 17 00:00:00 2001 From: Mateus Zanotto <39527468+mateuszanotto@users.noreply.github.com> Date: Tue, 13 Sep 2022 14:07:39 -0300 Subject: [PATCH 11/28] Added .itp (amber format) support --- qforce/forcefield.py | 9 +++++---- 1 file changed, 5 insertions(+), 4 deletions(-) diff --git a/qforce/forcefield.py b/qforce/forcefield.py index 3949c1e..8bd0499 100644 --- a/qforce/forcefield.py +++ b/qforce/forcefield.py @@ -167,7 +167,7 @@ def write_itp_atoms_and_molecule(self, itp, non_bonded): # atoms itp.write("\n[ atoms ]\n") - itp.write("; nr type resnr resnm atom cgrp charge mass\n") + itp.write("; nr type resnr resnm atom cgrp charge mass\n") for i, (lj_type, a_name, q, mass) in enumerate(zip(non_bonded.lj_types, self.atom_names, self.q, self.masses), start=1): itp.write(f'{i:>5} {lj_type:>9} {1:>6} {self.residue:>6} {a_name:>7} {i:>5} ') @@ -181,7 +181,7 @@ def write_itp_atoms_and_molecule(self, itp, non_bonded): def write_itp_polarization(self, itp, non_bonded): # polarization itp.write("\n[ polarization ]\n") - itp.write("; i j f alpha\n") + itp.write("; i j f alpha\n") for atom, drude in non_bonded.alpha_map.items(): alpha = non_bonded.alpha[atom]*1e-3 itp.write(f"{atom+1:>6} {drude+1:>6} {1:>6} {alpha:>14.8f}\n") @@ -226,6 +226,7 @@ def write_itp_angles(self, itp, terms): if self.urey: itp.write(' r0 k_bond') itp.write('\n') + for angle in terms['angle']: ids = angle.atomids + 1 equ = np.degrees(angle.equ) @@ -250,7 +251,7 @@ def write_itp_dihedrals(self, itp, terms): # rigid dihedrals if len(terms['dihedral/rigid']) > 0: itp.write("; rigid dihedrals \n") - itp.write("; ai aj ak al f theta0 k_theta\n") + itp.write("; ai aj ak al f theta0 k_theta\n") for dihed in terms['dihedral/rigid']: ids = dihed.atomids + 1 @@ -263,7 +264,7 @@ def write_itp_dihedrals(self, itp, terms): # improper dihedrals if len(terms['dihedral/improper']) > 0: itp.write("; improper dihedrals \n") - itp.write("; ai aj ak al f theta0 k_theta\n") + itp.write("; ai aj ak al f theta0 k_theta\n") for dihed in terms['dihedral/improper']: ids = dihed.atomids + 1 From 5671c9edfd4b93b4d334c53f4f600812fdadf504 Mon Sep 17 00:00:00 2001 From: Mateus Zanotto <39527468+mateuszanotto@users.noreply.github.com> Date: Thu, 15 Sep 2022 10:28:36 -0300 Subject: [PATCH 12/28] Added logo with hash (#) for .mol2 --- qforce/misc.py | 12 ++++++++++++ 1 file changed, 12 insertions(+) diff --git a/qforce/misc.py b/qforce/misc.py index 80a649f..eb2819d 100755 --- a/qforce/misc.py +++ b/qforce/misc.py @@ -27,6 +27,18 @@ ; ============================== """ +LOGO_HASH = """ +# ____ ______ +# / __ \ | ____| +# | | | |______| |__ ___ _ __ ___ ___ +# | | | |______| __/ _ \| '__/ __/ _ \\ +# | |__| | | | | (_) | | | (_| __/ +# \___\_\ |_| \___/|_| \___\___| +# +# Selim Sami +# University of Groningen - 2020 +# ============================== +""" def check_if_file_exists(filename): if not os.path.exists(filename) and not os.path.exists(f'{filename}_qforce'): From 92ca51e021a7181f9e04122c2900a6bb07588020 Mon Sep 17 00:00:00 2001 From: Mateus Zanotto <39527468+mateuszanotto@users.noreply.github.com> Date: Mon, 3 Oct 2022 10:05:27 -0300 Subject: [PATCH 13/28] Update qforce --- bin/qforce | 51 +-------------------------------------------------- 1 file changed, 1 insertion(+), 50 deletions(-) diff --git a/bin/qforce b/bin/qforce index b98fc2d..81896cd 100644 --- a/bin/qforce +++ b/bin/qforce @@ -1,56 +1,7 @@ #!/usr/bin/env python3 -import sys -import threading -import time - -class ProgressBarThread(threading.Thread): - def __init__(self, label='Working', delay=0.1): - super(ProgressBarThread, self).__init__() - self.label = label - self.delay = delay # interval between updates - self.running = False - def start(self): - self.running = True - super(ProgressBarThread, self).start() - def run(self): - label = '\r' + self.label + ' ' - while self.running: - # for c in ('-', '\\', '|', '/'): - # for c in ('▌', '▀', '▐', '▄'): - for c in ('◐', '◓', '◑', '◒'): - # for c in ('▙', '▛', '▜', '▟'): - # for c in ('▤', '▧', '▥', '▨'): - - sys.stdout.write(label + c) - sys.stdout.flush() - time.sleep(self.delay) - def stop(self): - self.running = False - self.join() # wait for run() method to terminate - sys.stdout.write('\r' + len(self.label)*' ' + 2*' ' + '\r') # clean-up - sys.stdout.flush() - -def work(): - time.sleep(5) # - -print(""" - ____ ______ - / __ \ | ____| - | | | |______| |__ ___ _ __ ___ ___ - | | | |______| __/ _ \| '__/ __/ _ \\ - | |__| | | | | (_) | | | (_| __/ - \___\_\ |_| \___/|_| \___\___| - - Selim Sami - University of Groningen - 2020 - ============================== -""") - -pb_thread = ProgressBarThread(' Initializing') -pb_thread.start() from qforce.main import run -pb_thread.stop() + if __name__ == '__main__': run() From da6ab9ff43ac7f7bc6944f20bcb27ae4ecd723d7 Mon Sep 17 00:00:00 2001 From: Mateus Zanotto <39527468+mateuszanotto@users.noreply.github.com> Date: Mon, 3 Oct 2022 10:17:43 -0300 Subject: [PATCH 14/28] Update forcefield.py --- qforce/forcefield.py | 118 ++++++++++++++----------------------------- 1 file changed, 37 insertions(+), 81 deletions(-) diff --git a/qforce/forcefield.py b/qforce/forcefield.py index 8bd0499..f135853 100644 --- a/qforce/forcefield.py +++ b/qforce/forcefield.py @@ -391,7 +391,6 @@ def set_charge(self, non_bonded): q[list(non_bonded.alpha_map.keys())] += 8 return q -## Amber def write_amber(self, directory, mol, coords): atom_ids, unique_at =self.get_atom_types(mol.topo, mol.non_bonded) self.write_mol2(directory, mol, coords, atom_ids, unique_at) @@ -407,39 +406,11 @@ def write_mol2(self, directory, mol, coords, atom_ids, unique_at): def write_frcmod(self, directory, mol, coords, atom_ids, unique_at): with open(f"{directory}/{self.mol_name}_qforce{self.polar_title}.frcmod", "w") as frcmod: self.write_frcmod_mass(frcmod, mol.non_bonded, atom_ids, unique_at) - self.write_frcmod_bonds(frcmod, mol.terms, mol.non_bonded.alpha_map, atom_ids, unique_at) + self.write_frcmod_bonds(frcmod, mol.terms, mol.non_bonded.alpha_map,atom_ids,unique_at) self.write_frcmod_angles(frcmod, mol.terms, atom_ids, unique_at) self.write_frcmod_dihedrals(frcmod, mol.terms, atom_ids, unique_at) self.write_frcmod_nonbond(frcmod, mol.non_bonded, atom_ids, unique_at) - -# TLeap is a helper program from AMBER to create the topologies. -# Would be good for the final user to create automatically a -# Tleap script with the new atom types. - - # def write_tleap_script(self, directory): - # with open(f"{directory}/tleap_script_{self.mol_name}.in", "w") as tleap: - # tleap.write(f""" - # # To run the script use tleap on amberTools and the following command: - # # $ tleap -f tleap_script_{self.mol_name}.in > tleap_script_{self.mol_name}.o - # - # source leaprc.protein.ff19SB\nsource leaprc.water.opc - # source leaprc.gaff\n# You can source other force fields here if necessary - # loadmol2 {self.mol_name}_qforce{self.polar_title}.mol2 - # - # loadamberparams {self.mol_name}_qforce{self.polar_title}.frcmod - # - # ## Add any other command, for example "solvateoct MOL TIP3PBOX 10.0" - # - # MOL = loadPDB 1lyd.solv.pdb - # saveamberparm MOL {self.mol_name}_qforce{self.polar_title}.prmtop - # saveamberparm MOL {self.mol_name}_qforce{self.polar_title}.inpcrd - # savepdb MOL {self.mol_name}_qforce{self.polar_title}.pdb """) - -############################################## -# MOL2 WRITING # -############################################## - def write_mol2_title(self, mol2): mol2.write(LOGO_HASH) mol2.write(f"# Name: {self.mol_name}\n") @@ -449,8 +420,8 @@ def write_mol2_molecule(self, mol2, topo, terms): 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(f"esp")## type of charge + mol2.write(f"{self.n_atoms:8d} {n_bonds:8d} {1:8d} {0:8d} {0:8d}\n") #n_bonds + mol2.write(f"esp") # type of charge mol2.write("\n\n") def write_mol2_atom(self, mol2, topo, coords, non_bonded, atom_ids, unique_at): @@ -458,7 +429,8 @@ def write_mol2_atom(self, mol2, topo, coords, non_bonded, atom_ids, unique_at): 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} {coords[i_idx-1][1]:10.4f} {coords[i_idx-1][2]:10.4f}") + 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" {unique_at[i_idx][atom_ids[i_idx]]} {1} {self.mol_name} {q:10.6f}\n") def write_mol2_bond(self, mol2, topo, terms, atom_ids, unique_at): @@ -471,10 +443,6 @@ def write_mol2_bond(self, mol2, topo, terms, atom_ids, unique_at): mol2.write(f"@SUBSTRUCTURE\n") mol2.write(f"{1:5} {self.mol_name} {1:8} TEMP {0:>8d} **** **** {0:5} ROOT\n") -############################################## -# FRCMOD WRITING # -############################################## - def write_frcmod_mass(self, frcmod, non_bonded, atom_ids, unique_at): frcmod.write(f'{self.mol_name} - frcmod generated by QForce\n') frcmod.write(f'MASS\n') @@ -485,8 +453,8 @@ def write_frcmod_bonds(self, frcmod, terms, alpha_map, atom_ids, unique_at): frcmod.write("\nBOND\n") for bond in terms['bond']: ids = bond.atomids + 1 - fconst = bond.fconst * (0.239005)/2 # kJ/mol/A^2 -> kcal/mol/A^2 - equ = bond.equ # Ang + fconst = bond.fconst * (0.239005)/2 # kJ -> kcal + equ = bond.equ frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}") frcmod.write(f"{fconst:>10.2f}{equ:>10.2f}\n") @@ -495,8 +463,8 @@ def write_frcmod_angles(self, frcmod, terms, atom_ids, unique_at): frcmod.write("\nANGLE\n") for angle in terms['angle']: ids = angle.atomids + 1 - fconst = angle.fconst * (0.239005)/2 # kJ/mol/rad^2 -> kcal/mol/rad^2 - equ = np.degrees(angle.equ) # Degrees -> Degrees + fconst = angle.fconst * (0.239005)/2 + equ = np.degrees(angle.equ) if self.urey: urey = [term for term in terms['urey'] if np.array_equal(term.atomids, @@ -508,7 +476,7 @@ def write_frcmod_angles(self, frcmod, terms, atom_ids, unique_at): frcmod.write(f"{fconst:>10.2f}{equ:>10.2f}\n") else: urey_equ = urey[0].equ - urey_fconst = urey[0].fconst * (0.239005)/2 # kJ/mol/rad^2 -> kcal/mol/rad^2 + urey_fconst = urey[0].fconst * (0.239005)/2 frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}") @@ -528,15 +496,14 @@ def write_frcmod_dihedrals(self, frcmod, terms, atom_ids, unique_at): if dihed.equ == 0.00 or dihed.equ==3.141592653589793: equ = 180.0 + fconst = dihed.fconst * (0.239005) - fconst = dihed.fconst * (0.239005) # kJ/mol -> kcal/mol #scale cos - - fconstAmb = (2*fconst/(2**2)) # Vn = (2*fcte/(n**2)) - frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") #Atom# 1 - frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") #Atom 2 - frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}-") #Atom 3 - frcmod.write(f"{unique_at[ids[3]][atom_ids[ids[3]]]:<2}") #Atom 4 - frcmod.write(f" 1 {fconstAmb:>15.2f} {equ:>15.2f} 2\n") #IDIVF=1, fcte, angle, n=2 (n*angle) + fconstAmb = (2*fconst/(2**2)) + frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") + frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") + frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}-") + frcmod.write(f"{unique_at[ids[3]][atom_ids[ids[3]]]:<2}") + frcmod.write(f" 1 {fconstAmb:>15.2f} {equ:>15.2f} 2\n") if len(terms['dihedral/flexible']) > 0: @@ -549,18 +516,18 @@ def write_frcmod_dihedrals(self, frcmod, terms, atom_ids, unique_at): equ = 180.0 else: equ = 0.00 - frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") #Atom 1 - frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") #Atom 2 - frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}-") #Atom 3 - frcmod.write(f"{unique_at[ids[3]][atom_ids[ids[3]]]:<2}") #Atom 4 - frcmod.write(f" 1 {c[n]:>15.2f} {equ:>15.2f} -{n+1}\n") #IDIVF=1, fcte, angle, n=1 (n*angle) + frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") + frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") + frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}-") + frcmod.write(f"{unique_at[ids[3]][atom_ids[ids[3]]]:<2}") + frcmod.write(f" 1 {c[n]:>15.2f} {equ:>15.2f} -{n+1}\n") # Last term without -n - frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") #Atom 1 - frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") #Atom 2 - frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}-") #Atom 3 - frcmod.write(f"{unique_at[ids[3]][atom_ids[ids[3]]]:<2}") #Atom 4 - frcmod.write(f" 1 {c[0]:>15.2f} {0.00:>15.2f} 1\n") #IDIVF=1, fcte, angle, n=1 (n*angle) + frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") + frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") + frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}-") + frcmod.write(f"{unique_at[ids[3]][atom_ids[ids[3]]]:<2}") + frcmod.write(f" 1 {c[0]:>15.2f} {0.00:>15.2f} 1\n") # improper dihedrals @@ -568,13 +535,13 @@ def write_frcmod_dihedrals(self, frcmod, terms, atom_ids, unique_at): frcmod.write("\nIMPROPER\n") for dihed in terms['dihedral/improper']: ids = dihed.atomids + 1 - equ = np.degrees(dihed.equ) # Degrees -> Degrees - fconst = dihed.fconst * (0.239005)/2 # kJ/mol -> kcal/mol -> f/2 as is Vn/2 in amber - frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") #Atom 1 - frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") #Atom 2 - frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}-") #Atom 3 - frcmod.write(f"{unique_at[ids[3]][atom_ids[ids[3]]]:<2}") #Atom 4 - frcmod.write(f" 1 {fconst:>15.2f} {equ:>15.2f} 2\n") #IDIVF=1, fcte, angle, n=2 (n*angle) + equ = np.degrees(dihed.equ) + fconst = dihed.fconst * (0.239005)/2 + frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") + frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") + frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}-") + frcmod.write(f"{unique_at[ids[3]][atom_ids[ids[3]]]:<2}") + frcmod.write(f" 1 {fconst:>15.2f} {equ:>15.2f} 2\n") def write_frcmod_nonbond(self, frcmod, non_bonded, atom_ids, unique_at): gro_atomtypes, gro_nonbonded, gro_1_4 = self.convert_to_gromacs_nonbonded(non_bonded) @@ -587,13 +554,9 @@ def write_frcmod_nonbond(self, frcmod, non_bonded, atom_ids, unique_at): frcmod.write(f"{gro_atomtypes[atom_ids[i]][1]/4.184:>12.4f}") frcmod.write(f"\t{atom_ids[i]}\n") - -############################################## -# GET ATOM TYPES AND IDS # -############################################## def get_atom_types(self, topo, non_bonded): atom_ids={} #dictonary containing original atom types - unique_at={} #dictonary containing new unique atom types + unique_at={} #dictonary containing new unique atom types unique_masses={} ascii_lowercase = list(string.ascii_lowercase) @@ -608,18 +571,11 @@ def get_atom_types(self, topo, non_bonded): elif i > 10: unique_at[i]={atom_ids[i] : "Q{}".format(ascii_lowercase[i-11])} else: - unique_at[i]={atom_ids[i] : "{}{}".format(ascii_uppercase[i-1], ascii_lowercase[i-1])} - - ## DEBUG : check correspondance between atom_ids and unique_at - #for i, lj_type in enumerate(zip(non_bonded.lj_types), start=1): - # print(unique_at) - # for i in range(1,self.n_atoms+1,1): - # print (i, atom_ids[i], unique_at[i][atom_ids[i]]) - # print ("------------------") + unique_at[i]={atom_ids[i] : "{}{}".format(ascii_uppercase[i-1], + ascii_lowercase[i-1])} return atom_ids, unique_at -############################################## # bohr2nm = 0.052917721067 # if polar: From 9690ff98f7adb5afd197b1500b96b13103e78360 Mon Sep 17 00:00:00 2001 From: Mateus Zanotto <39527468+mateuszanotto@users.noreply.github.com> Date: Mon, 3 Oct 2022 10:18:29 -0300 Subject: [PATCH 15/28] Update initialize.py --- qforce/initialize.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/qforce/initialize.py b/qforce/initialize.py index bfe0443..abc9c1e 100644 --- a/qforce/initialize.py +++ b/qforce/initialize.py @@ -168,7 +168,7 @@ def _check_and_copy_settings_file(job_dir, config_file): def initialize(filename, config_file, presets=None): - #print(LOGO) + print(LOGO) job_info = _get_job_info(filename) settings_file = _check_and_copy_settings_file(job_info.dir, config_file) From b2ca74d3bec3d0f862693b7ca27ffa9593f8ab9d Mon Sep 17 00:00:00 2001 From: Mateus Zanotto <39527468+mateuszanotto@users.noreply.github.com> Date: Mon, 3 Oct 2022 10:20:09 -0300 Subject: [PATCH 16/28] Update main.py --- qforce/main.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/qforce/main.py b/qforce/main.py index b14fc33..0201160 100755 --- a/qforce/main.py +++ b/qforce/main.py @@ -61,7 +61,7 @@ def run_qforce(input_arg, ext_q=None, ext_lj=None, config=None, presets=None): if forcefield in ["gromacs"]: ff.write_gromacs(job.dir, mol, qm_hessian_out.coords) - if forcefield in ["amber"]: + else: ff.write_amber(job.dir, mol, qm_hessian_out.coords) print_outcome(job.dir) From b6b1f0496862468baaf4cef781439cb05d357b26 Mon Sep 17 00:00:00 2001 From: Mateus Zanotto <39527468+mateuszanotto@users.noreply.github.com> Date: Mon, 3 Oct 2022 10:44:18 -0300 Subject: [PATCH 17/28] Update forcefield.py --- qforce/forcefield.py | 1 - 1 file changed, 1 deletion(-) diff --git a/qforce/forcefield.py b/qforce/forcefield.py index f135853..7eb2769 100644 --- a/qforce/forcefield.py +++ b/qforce/forcefield.py @@ -312,7 +312,6 @@ def make_pairs(self, neighbors, non_bonded): polar_pairs = [] if self.n_excl == 2: - if self.polar: for a1, a2 in non_bonded.pairs: if a2 in non_bonded.alpha_map.keys(): From 7924090131fa0c22a17ad7964d08d7cde799e036 Mon Sep 17 00:00:00 2001 From: Mateus Zanotto <39527468+mateuszanotto@users.noreply.github.com> Date: Mon, 3 Oct 2022 10:44:52 -0300 Subject: [PATCH 18/28] Update main.py --- qforce/main.py | 8 +------- 1 file changed, 1 insertion(+), 7 deletions(-) diff --git a/qforce/main.py b/qforce/main.py index 0201160..061f254 100755 --- a/qforce/main.py +++ b/qforce/main.py @@ -56,13 +56,7 @@ 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) - - forcefield = config.ff.forcefield - if forcefield in ["gromacs"]: - ff.write_gromacs(job.dir, mol, qm_hessian_out.coords) - - else: - ff.write_amber(job.dir, mol, qm_hessian_out.coords) + ff.write_ff(config.ff.forcefield, job.dir, mol, qm_hessian_out.coords) print_outcome(job.dir) From 68126438b83d7e475cf2433e0f2b4d1d36ac14ad Mon Sep 17 00:00:00 2001 From: Mateus Zanotto <39527468+mateuszanotto@users.noreply.github.com> Date: Mon, 3 Oct 2022 10:45:17 -0300 Subject: [PATCH 19/28] Update forcefield.py --- qforce/forcefield.py | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/qforce/forcefield.py b/qforce/forcefield.py index 7eb2769..d0ad6ea 100644 --- a/qforce/forcefield.py +++ b/qforce/forcefield.py @@ -32,6 +32,12 @@ def __init__(self, job_name, config, mol, neighbors, exclude_all=[]): else: self.polar_title = '' + def write_ff(self, forcefield, directory, mol, coords): + if forcefield == "gromacs": + self.write_gromacs(directory, mol, coords) + else: + self.write_amber(directory, mol, coords) + def write_gromacs(self, directory, mol, coords): self.write_itp(mol, directory) self.write_top(directory) From 93ad3de56804dca1b78e5a7f98ac6b07931598db Mon Sep 17 00:00:00 2001 From: Selim Sami Date: Fri, 7 Oct 2022 14:06:22 -0700 Subject: [PATCH 20/28] cosmetic changes --- qforce/forcefield.py | 158 ++++++++++--------------------------------- qforce/initialize.py | 6 +- qforce/main.py | 18 ++--- qforce/misc.py | 31 ++------- 4 files changed, 51 insertions(+), 162 deletions(-) diff --git a/qforce/forcefield.py b/qforce/forcefield.py index d0ad6ea..2910185 100644 --- a/qforce/forcefield.py +++ b/qforce/forcefield.py @@ -1,12 +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 LOGO_HASH -import string +from .misc import get_logo class ForceField(): @@ -27,17 +25,16 @@ 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: self.polar_title = '' - def write_ff(self, forcefield, directory, mol, coords): - if forcefield == "gromacs": - self.write_gromacs(directory, mol, coords) - else: - self.write_amber(directory, mol, coords) - def write_gromacs(self, directory, mol, coords): self.write_itp(mol, directory) self.write_top(directory) @@ -84,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) @@ -137,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} ' + '{0:>8.4f} {"A":>5} ') itp.write(f'{lj_params[0]:>12.5e} {lj_params[1]:>12.5e}\n') if self.polar: @@ -397,7 +395,7 @@ def set_charge(self, non_bonded): return q def write_amber(self, directory, mol, coords): - atom_ids, unique_at =self.get_atom_types(mol.topo, mol.non_bonded) + atom_ids, unique_at = self.get_atom_types(mol.topo, mol.non_bonded) self.write_mol2(directory, mol, coords, atom_ids, unique_at) self.write_frcmod(directory, mol, coords, atom_ids, unique_at) @@ -411,13 +409,14 @@ def write_mol2(self, directory, mol, coords, atom_ids, unique_at): def write_frcmod(self, directory, mol, coords, atom_ids, unique_at): with open(f"{directory}/{self.mol_name}_qforce{self.polar_title}.frcmod", "w") as frcmod: self.write_frcmod_mass(frcmod, mol.non_bonded, atom_ids, unique_at) - self.write_frcmod_bonds(frcmod, mol.terms, mol.non_bonded.alpha_map,atom_ids,unique_at) + self.write_frcmod_bonds(frcmod, mol.terms, mol.non_bonded.alpha_map, atom_ids, + unique_at) self.write_frcmod_angles(frcmod, mol.terms, atom_ids, unique_at) self.write_frcmod_dihedrals(frcmod, mol.terms, atom_ids, unique_at) self.write_frcmod_nonbond(frcmod, mol.non_bonded, atom_ids, unique_at) def write_mol2_title(self, mol2): - mol2.write(LOGO_HASH) + mol2.write(get_logo('#')) mol2.write(f"# Name: {self.mol_name}\n") def write_mol2_molecule(self, mol2, topo, terms): @@ -425,14 +424,15 @@ def write_mol2_molecule(self, mol2, topo, terms): 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(f"esp") # type of charge + 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, atom_ids, unique_at): - mol2.write(f"@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("@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}") @@ -440,17 +440,17 @@ def write_mol2_atom(self, mol2, topo, coords, non_bonded, atom_ids, unique_at): def write_mol2_bond(self, mol2, topo, terms, atom_ids, unique_at): n_bonds = 1 - mol2.write(f"@BOND\n") + 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(f"@SUBSTRUCTURE\n") + mol2.write("@SUBSTRUCTURE\n") mol2.write(f"{1:5} {self.mol_name} {1:8} TEMP {0:>8d} **** **** {0:5} ROOT\n") def write_frcmod_mass(self, frcmod, non_bonded, atom_ids, unique_at): frcmod.write(f'{self.mol_name} - frcmod generated by QForce\n') - frcmod.write(f'MASS\n') + frcmod.write('MASS\n') for i, mass in enumerate((self.masses), start=1): frcmod.write(f"{unique_at[i][atom_ids[i]]} {mass}\n") @@ -458,7 +458,7 @@ def write_frcmod_bonds(self, frcmod, terms, alpha_map, atom_ids, unique_at): frcmod.write("\nBOND\n") for bond in terms['bond']: ids = bond.atomids + 1 - fconst = bond.fconst * (0.239005)/2 # kJ -> kcal + fconst = bond.fconst * (0.239005)/2 # kJ -> kcal equ = bond.equ frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}") @@ -498,7 +498,7 @@ def write_frcmod_dihedrals(self, frcmod, terms, atom_ids, unique_at): ids = dihed.atomids + 1 equ = dihed.equ - if dihed.equ == 0.00 or dihed.equ==3.141592653589793: + if dihed.equ == 0.00 or dihed.equ == 3.141592653589793: equ = 180.0 fconst = dihed.fconst * (0.239005) @@ -517,7 +517,7 @@ def write_frcmod_dihedrals(self, frcmod, terms, atom_ids, unique_at): c = dihed.equ for n in range(3, 0, -1): - if n%2==1: + if n % 2 == 1: equ = 180.0 else: equ = 0.00 @@ -534,7 +534,6 @@ def write_frcmod_dihedrals(self, frcmod, terms, atom_ids, unique_at): frcmod.write(f"{unique_at[ids[3]][atom_ids[ids[3]]]:<2}") frcmod.write(f" 1 {c[0]:>15.2f} {0.00:>15.2f} 1\n") - # improper dihedrals if len(terms['dihedral/improper']) > 0: frcmod.write("\nIMPROPER\n") @@ -553,118 +552,29 @@ def write_frcmod_nonbond(self, frcmod, non_bonded, atom_ids, unique_at): frcmod.write("\nNONB\n") # support to gaff and gaff2 if self.comb_rule == 2: - for i in range(1,self.n_atoms+1,1): + for i in range(1, self.n_atoms+1, 1): frcmod.write(f"{unique_at[i][atom_ids[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") def get_atom_types(self, topo, non_bonded): - atom_ids={} #dictonary containing original atom types - unique_at={} #dictonary containing new unique atom types - unique_masses={} + atom_ids = {} # original atom types + unique_at = {} # new unique atom types ascii_lowercase = list(string.ascii_lowercase) ascii_digits = list(string.digits) ascii_uppercase = list(string.ascii_uppercase) for i, (lj_type, a_name) in enumerate(zip(non_bonded.lj_types, self.atom_names), start=1): - atom_ids[i]=lj_type + atom_ids[i] = lj_type if self.n_atoms < 36: if i <= 10: - unique_at[i]={atom_ids[i] : "Q{}".format(ascii_digits[i-1])} + unique_at[i] = {atom_ids[i]: "Q{}".format(ascii_digits[i-1])} elif i > 10: - unique_at[i]={atom_ids[i] : "Q{}".format(ascii_lowercase[i-11])} + unique_at[i] = {atom_ids[i]: "Q{}".format(ascii_lowercase[i-11])} else: - unique_at[i]={atom_ids[i] : "{}{}".format(ascii_uppercase[i-1], - ascii_lowercase[i-1])} + unique_at[i] = {atom_ids[i]: "{}{}".format(ascii_uppercase[i-1], + ascii_lowercase[i-1])} return atom_ids, unique_at - - - # 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])])) diff --git a/qforce/initialize.py b/qforce/initialize.py index abc9c1e..e093cf7 100644 --- a/qforce/initialize.py +++ b/qforce/initialize.py @@ -9,7 +9,7 @@ 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): @@ -17,7 +17,7 @@ class Initialize(Colt): [ff] # MD software to use # Currently amber supports only gaff2 atom types -forcefield = gromacs :: str :: gromacs, amber +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 @@ -168,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 061f254..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_ff(config.ff.forcefield, 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 eb2819d..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 = """ ; ____ ______ ; / __ \ | ____| ; | | | |______| |__ ___ _ __ ___ ___ @@ -27,18 +14,10 @@ ; ============================== """ -LOGO_HASH = """ -# ____ ______ -# / __ \ | ____| -# | | | |______| |__ ___ _ __ ___ ___ -# | | | |______| __/ _ \| '__/ __/ _ \\ -# | |__| | | | | (_) | | | (_| __/ -# \___\_\ |_| \___/|_| \___\___| -# -# Selim Sami -# University of Groningen - 2020 -# ============================== -""" + +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'): From 0a5c98a34a32f3f2f7b91bf0e0e1d129178208be Mon Sep 17 00:00:00 2001 From: Selim Sami Date: Fri, 7 Oct 2022 14:55:40 -0700 Subject: [PATCH 21/28] missing f-string --- qforce/forcefield.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/qforce/forcefield.py b/qforce/forcefield.py index 2910185..41b6ca7 100644 --- a/qforce/forcefield.py +++ b/qforce/forcefield.py @@ -135,7 +135,7 @@ def write_itp_atoms_and_molecule(self, itp, non_bonded): 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} ') + f'{0:>8.4f} {"A":>5} ') itp.write(f'{lj_params[0]:>12.5e} {lj_params[1]:>12.5e}\n') if self.polar: From ebb54fb0c60ff7d88a92524714d4b88a5eca6b16 Mon Sep 17 00:00:00 2001 From: Selim Sami Date: Fri, 7 Oct 2022 15:03:20 -0700 Subject: [PATCH 22/28] bump version --- setup.cfg | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) 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 From b4d2dc95ad85e16c719f193856f62f8394e79a47 Mon Sep 17 00:00:00 2001 From: Mateus Zanotto <39527468+mateuszanotto@users.noreply.github.com> Date: Thu, 27 Oct 2022 13:54:49 -0300 Subject: [PATCH 23/28] Update forcefield.py Removed urey angles from amber file format --- qforce/forcefield.py | 2 -- 1 file changed, 2 deletions(-) diff --git a/qforce/forcefield.py b/qforce/forcefield.py index 41b6ca7..617aa66 100644 --- a/qforce/forcefield.py +++ b/qforce/forcefield.py @@ -480,8 +480,6 @@ def write_frcmod_angles(self, frcmod, terms, atom_ids, unique_at): frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}") frcmod.write(f"{fconst:>10.2f}{equ:>10.2f}\n") else: - urey_equ = urey[0].equ - urey_fconst = urey[0].fconst * (0.239005)/2 frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}") From a9e87391b06976b17937d2ab1bdeaf3db6a07373 Mon Sep 17 00:00:00 2001 From: Mateus Zanotto <39527468+mateuszanotto@users.noreply.github.com> Date: Wed, 5 Apr 2023 14:29:29 -0300 Subject: [PATCH 24/28] Update forcefield.py Added Urey-Bradley potential Removed Fourier series dihedrals --- qforce/forcefield.py | 49 ++++++++++++++++++++++++-------------------- 1 file changed, 27 insertions(+), 22 deletions(-) diff --git a/qforce/forcefield.py b/qforce/forcefield.py index 617aa66..36ac4c0 100644 --- a/qforce/forcefield.py +++ b/qforce/forcefield.py @@ -436,7 +436,7 @@ def write_mol2_atom(self, mol2, topo, coords, non_bonded, atom_ids, unique_at): 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" {unique_at[i_idx][atom_ids[i_idx]]} {1} {self.mol_name} {q:10.6f}\n") + mol2.write(f" {unique_at[i_idx][atom_ids[i_idx]]} {1} {self.mol_name[0:3]} {q:10.6f}\n") def write_mol2_bond(self, mol2, topo, terms, atom_ids, unique_at): n_bonds = 1 @@ -446,7 +446,7 @@ def write_mol2_bond(self, mol2, topo, terms, atom_ids, unique_at): 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} {1:8} TEMP {0:>8d} **** **** {0:5} ROOT\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, atom_ids, unique_at): frcmod.write(f'{self.mol_name} - frcmod generated by QForce\n') @@ -458,17 +458,32 @@ def write_frcmod_bonds(self, frcmod, terms, alpha_map, atom_ids, unique_at): frcmod.write("\nBOND\n") for bond in terms['bond']: ids = bond.atomids + 1 - fconst = bond.fconst * (0.239005)/2 # kJ -> kcal + fconst = bond.fconst * (0.2390057361)/2 # kJ -> kcal equ = bond.equ frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}") frcmod.write(f"{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 * (0.2390057361)/2 + frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") + frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}") + frcmod.write(f"{urey_fconst:>10.2f}{urey_equ:>10.2f}\n") + + def write_frcmod_angles(self, frcmod, terms, atom_ids, unique_at): frcmod.write("\nANGLE\n") for angle in terms['angle']: ids = angle.atomids + 1 - fconst = angle.fconst * (0.239005)/2 + fconst = angle.fconst * (0.2390057361)/2 equ = np.degrees(angle.equ) if self.urey: @@ -478,12 +493,13 @@ def write_frcmod_angles(self, frcmod, terms, atom_ids, unique_at): frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}") - frcmod.write(f"{fconst:>10.2f}{equ:>10.2f}\n") + frcmod.write(f"{fconst:>10.2f}{equ:>10.2f} \n") else: frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}") - frcmod.write(f"{fconst:>10.2f}{equ:>10.2f}\n") + frcmod.write(f"{fconst:>10.2f}{equ:>10.2f} \n") + def write_frcmod_dihedrals(self, frcmod, terms, atom_ids, unique_at): if len(terms['dihedral']) > 0: @@ -499,7 +515,7 @@ def write_frcmod_dihedrals(self, frcmod, terms, atom_ids, unique_at): if dihed.equ == 0.00 or dihed.equ == 3.141592653589793: equ = 180.0 - fconst = dihed.fconst * (0.239005) + fconst = dihed.fconst * (0.2390057361) fconstAmb = (2*fconst/(2**2)) frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") @@ -510,27 +526,16 @@ def write_frcmod_dihedrals(self, frcmod, terms, atom_ids, unique_at): if len(terms['dihedral/flexible']) > 0: + # flexible dihedrals for dihed in terms['dihedral/flexible']: ids = dihed.atomids + 1 - c = dihed.equ - - for n in range(3, 0, -1): - if n % 2 == 1: - equ = 180.0 - else: - equ = 0.00 - frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") - frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") - frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}-") - frcmod.write(f"{unique_at[ids[3]][atom_ids[ids[3]]]:<2}") - frcmod.write(f" 1 {c[n]:>15.2f} {equ:>15.2f} -{n+1}\n") + c = dihed.equ/2 - # Last term without -n frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}-") frcmod.write(f"{unique_at[ids[3]][atom_ids[ids[3]]]:<2}") - frcmod.write(f" 1 {c[0]:>15.2f} {0.00:>15.2f} 1\n") + frcmod.write(f" 1 {c[0]:>15.2f} {180.00:>15.2f} 1\n") # improper dihedrals if len(terms['dihedral/improper']) > 0: @@ -538,7 +543,7 @@ def write_frcmod_dihedrals(self, frcmod, terms, atom_ids, unique_at): for dihed in terms['dihedral/improper']: ids = dihed.atomids + 1 equ = np.degrees(dihed.equ) - fconst = dihed.fconst * (0.239005)/2 + fconst = dihed.fconst * (0.2390057361)/2 frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}-") From e844a73f777b363553af826ee621cc7660208c12 Mon Sep 17 00:00:00 2001 From: Mateus Zanotto <39527468+mateuszanotto@users.noreply.github.com> Date: Thu, 13 Apr 2023 13:59:24 -0300 Subject: [PATCH 25/28] Update forcefield.py Conversion of RB terms --- qforce/forcefield.py | 64 ++++++++++++++++++++++++++++++++++++++------ 1 file changed, 56 insertions(+), 8 deletions(-) diff --git a/qforce/forcefield.py b/qforce/forcefield.py index 36ac4c0..f2aa5ff 100644 --- a/qforce/forcefield.py +++ b/qforce/forcefield.py @@ -530,12 +530,60 @@ def write_frcmod_dihedrals(self, frcmod, terms, atom_ids, unique_at): for dihed in terms['dihedral/flexible']: ids = dihed.atomids + 1 c = dihed.equ/2 + 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(5, 0, -1): + if n % 2 == 1: + equ = 180.0 + else: + equ = 0.00 + frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") + frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") + frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}-") + frcmod.write(f"{unique_at[ids[3]][atom_ids[ids[3]]]:<2}") + frcmod.write(f" 1 {fc[n]:>15.2f} {equ:>15.2f} -{n+1}\n") + # Last term frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}-") frcmod.write(f"{unique_at[ids[3]][atom_ids[ids[3]]]:<2}") - frcmod.write(f" 1 {c[0]:>15.2f} {180.00:>15.2f} 1\n") + frcmod.write(f" 1 {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*c[0]/2 + 0.5*c[2]/2, + 1.0*c[1]/2, + 0.5*c[2]/2] + + for n in range(3, 0, -1): + if n % 2 == 1: + equ = 180.0 + else: + equ = 0.00 + frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") + frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") + frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}-") + frcmod.write(f"{unique_at[ids[3]][atom_ids[ids[3]]]:<2}") + frcmod.write(f" 1 {fc[n]:>15.2f} {equ:>15.2f} -{n+1}\n") + + # Last term + frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") + frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") + frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}-") + frcmod.write(f"{unique_at[ids[3]][atom_ids[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: @@ -571,13 +619,13 @@ def get_atom_types(self, topo, non_bonded): for i, (lj_type, a_name) in enumerate(zip(non_bonded.lj_types, self.atom_names), start=1): atom_ids[i] = lj_type - if self.n_atoms < 36: - if i <= 10: - unique_at[i] = {atom_ids[i]: "Q{}".format(ascii_digits[i-1])} - elif i > 10: - unique_at[i] = {atom_ids[i]: "Q{}".format(ascii_lowercase[i-11])} + if self.n_atoms < 35: + if i <= 9: + unique_at[i] = {atom_ids[i]: "Q{}".format(ascii_digits[i])} + elif i > 9: + unique_at[i] = {atom_ids[i]: "Q{}".format(ascii_lowercase[i-10])} else: - unique_at[i] = {atom_ids[i]: "{}{}".format(ascii_uppercase[i-1], - ascii_lowercase[i-1])} + unique_at[i] = {atom_ids[i]: "{}{}".format(ascii_uppercase[i], + ascii_lowercase[i])} return atom_ids, unique_at From 4b3c4eb1f4038313eed44015e99974d7db4fe16e Mon Sep 17 00:00:00 2001 From: Mateus Zanotto <39527468+mateuszanotto@users.noreply.github.com> Date: Tue, 18 Apr 2023 11:15:31 -0300 Subject: [PATCH 26/28] Update forcefield.py Replaced unique_at[] with ids[] --- qforce/forcefield.py | 127 ++++++++++++++----------------------------- 1 file changed, 42 insertions(+), 85 deletions(-) diff --git a/qforce/forcefield.py b/qforce/forcefield.py index f2aa5ff..91391e1 100644 --- a/qforce/forcefield.py +++ b/qforce/forcefield.py @@ -395,25 +395,24 @@ def set_charge(self, non_bonded): return q def write_amber(self, directory, mol, coords): - atom_ids, unique_at = self.get_atom_types(mol.topo, mol.non_bonded) - self.write_mol2(directory, mol, coords, atom_ids, unique_at) - self.write_frcmod(directory, mol, coords, atom_ids, unique_at) + # atom_ids, unique_at = self.get_atom_types(mol.topo, mol.non_bonded) + self.write_mol2(directory, mol, coords) + self.write_frcmod(directory, mol, coords) - def write_mol2(self, directory, mol, coords, atom_ids, unique_at): + 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, atom_ids, unique_at) - self.write_mol2_bond(mol2, mol.topo, mol.terms, atom_ids, unique_at) + 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, atom_ids, unique_at): + def write_frcmod(self, directory, mol, coords): with open(f"{directory}/{self.mol_name}_qforce{self.polar_title}.frcmod", "w") as frcmod: - self.write_frcmod_mass(frcmod, mol.non_bonded, atom_ids, unique_at) - self.write_frcmod_bonds(frcmod, mol.terms, mol.non_bonded.alpha_map, atom_ids, - unique_at) - self.write_frcmod_angles(frcmod, mol.terms, atom_ids, unique_at) - self.write_frcmod_dihedrals(frcmod, mol.terms, atom_ids, unique_at) - self.write_frcmod_nonbond(frcmod, mol.non_bonded, atom_ids, unique_at) + self.write_frcmod_mass(frcmod, mol.non_bonded) + self.write_frcmod_bonds(frcmod, mol.terms, mol.non_bonded.alpha_map) + self.write_frcmod_angles(frcmod, mol.terms) + self.write_frcmod_dihedrals(frcmod, mol.terms) + self.write_frcmod_nonbond(frcmod, mol.non_bonded) def write_mol2_title(self, mol2): mol2.write(get_logo('#')) @@ -428,7 +427,7 @@ def write_mol2_molecule(self, mol2, topo, terms): mol2.write("esp") # type of charge mol2.write("\n\n") - def write_mol2_atom(self, mol2, topo, coords, non_bonded, atom_ids, unique_at): + 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, @@ -436,9 +435,9 @@ def write_mol2_atom(self, mol2, topo, coords, non_bonded, atom_ids, unique_at): 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" {unique_at[i_idx][atom_ids[i_idx]]} {1} {self.mol_name[0:3]} {q:10.6f}\n") + mol2.write(f" {i_idx} {1} {self.mol_name[0:3]} {q:10.6f}\n") - def write_mol2_bond(self, mol2, topo, terms, atom_ids, unique_at): + def write_mol2_bond(self, mol2, topo, terms): n_bonds = 1 mol2.write("@BOND\n") for bond in terms['bond']: @@ -448,21 +447,19 @@ def write_mol2_bond(self, mol2, topo, terms, atom_ids, unique_at): 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, atom_ids, unique_at): + 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"{unique_at[i][atom_ids[i]]} {mass}\n") + frcmod.write(f"{i} {mass}\n") - def write_frcmod_bonds(self, frcmod, terms, alpha_map, atom_ids, unique_at): + def write_frcmod_bonds(self, frcmod, terms, alpha_map): frcmod.write("\nBOND\n") for bond in terms['bond']: ids = bond.atomids + 1 fconst = bond.fconst * (0.2390057361)/2 # kJ -> kcal equ = bond.equ - frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") - frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}") - frcmod.write(f"{fconst:>10.2f}{equ:>10.2f}\n") + 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 @@ -474,12 +471,10 @@ def write_frcmod_bonds(self, frcmod, terms, alpha_map, atom_ids, unique_at): else: urey_equ = urey[0].equ urey_fconst = urey[0].fconst * (0.2390057361)/2 - frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") - frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}") - frcmod.write(f"{urey_fconst:>10.2f}{urey_equ:>10.2f}\n") + 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, atom_ids, unique_at): + def write_frcmod_angles(self, frcmod, terms): frcmod.write("\nANGLE\n") for angle in terms['angle']: ids = angle.atomids + 1 @@ -490,18 +485,14 @@ def write_frcmod_angles(self, frcmod, terms, atom_ids, unique_at): 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"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") - frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") - frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}") + 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"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") - frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") - frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}") + 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, atom_ids, unique_at): + def write_frcmod_dihedrals(self, frcmod, terms): if len(terms['dihedral']) > 0: frcmod.write("\nDIHE\n") @@ -518,18 +509,15 @@ def write_frcmod_dihedrals(self, frcmod, terms, atom_ids, unique_at): fconst = dihed.fconst * (0.2390057361) fconstAmb = (2*fconst/(2**2)) - frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") - frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") - frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}-") - frcmod.write(f"{unique_at[ids[3]][atom_ids[ids[3]]]:<2}") - frcmod.write(f" 1 {fconstAmb:>15.2f} {equ:>15.2f} 2\n") + 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/2 + c = dihed.equ 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], @@ -537,23 +525,17 @@ def write_frcmod_dihedrals(self, frcmod, terms, atom_ids, unique_at): 0.125*c[4], 0.0625*c[5]] - for n in range(5, 0, -1): + for n in range(3, 0, -1): if n % 2 == 1: equ = 180.0 else: equ = 0.00 - frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") - frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") - frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}-") - frcmod.write(f"{unique_at[ids[3]][atom_ids[ids[3]]]:<2}") - frcmod.write(f" 1 {fc[n]:>15.2f} {equ:>15.2f} -{n+1}\n") + 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"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") - frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") - frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}-") - frcmod.write(f"{unique_at[ids[3]][atom_ids[ids[3]]]:<2}") - frcmod.write(f" 1 {fc[0]:>15.2f} {0.00:>15.2f} 1\n") + 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: @@ -571,17 +553,11 @@ def write_frcmod_dihedrals(self, frcmod, terms, atom_ids, unique_at): equ = 180.0 else: equ = 0.00 - frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") - frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") - frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}-") - frcmod.write(f"{unique_at[ids[3]][atom_ids[ids[3]]]:<2}") + 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"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") - frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") - frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}-") - frcmod.write(f"{unique_at[ids[3]][atom_ids[ids[3]]]:<2}") + 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") @@ -592,40 +568,21 @@ def write_frcmod_dihedrals(self, frcmod, terms, atom_ids, unique_at): ids = dihed.atomids + 1 equ = np.degrees(dihed.equ) fconst = dihed.fconst * (0.2390057361)/2 - frcmod.write(f"{unique_at[ids[0]][atom_ids[ids[0]]]:<2}-") - frcmod.write(f"{unique_at[ids[1]][atom_ids[ids[1]]]:<2}-") - frcmod.write(f"{unique_at[ids[2]][atom_ids[ids[2]]]:<2}-") - frcmod.write(f"{unique_at[ids[3]][atom_ids[ids[3]]]:<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, atom_ids, unique_at): + 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"{unique_at[i][atom_ids[i]]:<8}") + 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") - - def get_atom_types(self, topo, non_bonded): - atom_ids = {} # original atom types - unique_at = {} # new unique atom types - - ascii_lowercase = list(string.ascii_lowercase) - ascii_digits = list(string.digits) - ascii_uppercase = list(string.ascii_uppercase) - - for i, (lj_type, a_name) in enumerate(zip(non_bonded.lj_types, self.atom_names), start=1): - atom_ids[i] = lj_type - if self.n_atoms < 35: - if i <= 9: - unique_at[i] = {atom_ids[i]: "Q{}".format(ascii_digits[i])} - elif i > 9: - unique_at[i] = {atom_ids[i]: "Q{}".format(ascii_lowercase[i-10])} - else: - unique_at[i] = {atom_ids[i]: "{}{}".format(ascii_uppercase[i], - ascii_lowercase[i])} - - return atom_ids, unique_at From a72c1c7d1cd17f35c0e534508e9490b0f58026aa Mon Sep 17 00:00:00 2001 From: Mateus Zanotto <39527468+mateuszanotto@users.noreply.github.com> Date: Fri, 12 May 2023 16:25:46 -0300 Subject: [PATCH 27/28] Update forcefield.py dihedrals converted from kJ to kcal/mol --- qforce/forcefield.py | 36 ++++++++++++++++++------------------ 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/qforce/forcefield.py b/qforce/forcefield.py index 91391e1..2acda29 100644 --- a/qforce/forcefield.py +++ b/qforce/forcefield.py @@ -395,9 +395,9 @@ def set_charge(self, non_bonded): return q def write_amber(self, directory, mol, coords): - # atom_ids, unique_at = self.get_atom_types(mol.topo, mol.non_bonded) + j2cal = 0.2390057361 self.write_mol2(directory, mol, coords) - self.write_frcmod(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: @@ -406,12 +406,12 @@ def write_mol2(self, directory, mol, coords): 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): + 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) - self.write_frcmod_angles(frcmod, mol.terms) - self.write_frcmod_dihedrals(frcmod, mol.terms) + 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): @@ -453,11 +453,11 @@ def write_frcmod_mass(self, frcmod, non_bonded): for i, mass in enumerate((self.masses), start=1): frcmod.write(f"{i} {mass}\n") - def write_frcmod_bonds(self, frcmod, terms, alpha_map): + 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 * (0.2390057361)/2 # kJ -> kcal + 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") @@ -470,15 +470,15 @@ def write_frcmod_bonds(self, frcmod, terms, alpha_map): None else: urey_equ = urey[0].equ - urey_fconst = urey[0].fconst * (0.2390057361)/2 + 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): + def write_frcmod_angles(self, frcmod, terms, j2cal): frcmod.write("\nANGLE\n") for angle in terms['angle']: ids = angle.atomids + 1 - fconst = angle.fconst * (0.2390057361)/2 + fconst = angle.fconst * j2cal/2 equ = np.degrees(angle.equ) if self.urey: @@ -492,7 +492,7 @@ def write_frcmod_angles(self, frcmod, terms): frcmod.write(f"{fconst:>10.2f}{equ:>10.2f} \n") - def write_frcmod_dihedrals(self, frcmod, terms): + def write_frcmod_dihedrals(self, frcmod, terms, j2cal): if len(terms['dihedral']) > 0: frcmod.write("\nDIHE\n") @@ -506,7 +506,7 @@ def write_frcmod_dihedrals(self, frcmod, terms): if dihed.equ == 0.00 or dihed.equ == 3.141592653589793: equ = 180.0 - fconst = dihed.fconst * (0.2390057361) + fconst = dihed.fconst * j2cal fconstAmb = (2*fconst/(2**2)) frcmod.write(f"{ids[0]:<2}-{ids[1]:<2}-{ids[2]:<2}-{ids[3]:<2}") @@ -517,7 +517,7 @@ def write_frcmod_dihedrals(self, frcmod, terms): # flexible dihedrals for dihed in terms['dihedral/flexible']: ids = dihed.atomids + 1 - c = dihed.equ + 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], @@ -544,9 +544,9 @@ def write_frcmod_dihedrals(self, frcmod, terms): ids = dihed.atomids + 1 c0, c1, c2 = convert_to_inversion_rb(dihed.fconst, dihed.equ) - fc = [1.0*c[0]/2 + 0.5*c[2]/2, - 1.0*c[1]/2, - 0.5*c[2]/2] + 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: @@ -567,7 +567,7 @@ def write_frcmod_dihedrals(self, frcmod, terms): for dihed in terms['dihedral/improper']: ids = dihed.atomids + 1 equ = np.degrees(dihed.equ) - fconst = dihed.fconst * (0.2390057361)/2 + 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") From 0da7dcfe6d37d83416056feffd77a4120ca6d0d5 Mon Sep 17 00:00:00 2001 From: Mateus Zanotto <39527468+mateuszanotto@users.noreply.github.com> Date: Fri, 12 May 2023 16:55:51 -0300 Subject: [PATCH 28/28] Update forcefield.py Set a note to turn off urey if using amber --- qforce/forcefield.py | 3 +++ 1 file changed, 3 insertions(+) diff --git a/qforce/forcefield.py b/qforce/forcefield.py index 2acda29..ba05620 100644 --- a/qforce/forcefield.py +++ b/qforce/forcefield.py @@ -395,6 +395,9 @@ def set_charge(self, non_bonded): return q 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)