Amber

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()

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=[]):
@@ -309,6 +312,7 @@ 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():
@@ -387,6 +391,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"@<TRIPOS>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"@<TRIPOS>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"@<TRIPOS>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"@<TRIPOS>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

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
@@ -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)

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)
 

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'):