Skip to content

Amber #62

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
mateuszanotto wants to merge 28 commits into
base: amber
Choose a base branch
from
Open

Amber #62

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
28 commits
Select commit Hold shift + click to select a range
ce7eb1b
Add files via upload
mateuszanotto Sep 13, 2022
c4a2d14
Delete forcefield.py
mateuszanotto Sep 13, 2022
a60171e
Delete initialize.py
mateuszanotto Sep 13, 2022
b02618e
Delete main.py
mateuszanotto Sep 13, 2022
627d65f
Add files via upload
mateuszanotto Sep 13, 2022
f37c4b5
Add files via upload
mateuszanotto Sep 13, 2022
d7f77e9
Added .itp (amber ff format) support
mateuszanotto Sep 13, 2022
ba2669d
Add a progress bar while import modules
mateuszanotto Sep 13, 2022
1acfc9c
Added amber ff option
mateuszanotto Sep 13, 2022
1ea47a4
Added .itp (amber format) support
mateuszanotto Sep 13, 2022
d66e274
Added .itp (amber format) support
mateuszanotto Sep 13, 2022
5671c9e
Added logo with hash (#) for .mol2
mateuszanotto Sep 15, 2022
92ca51e
Update qforce
mateuszanotto Oct 3, 2022
da6ab9f
Update forcefield.py
mateuszanotto Oct 3, 2022
9690ff9
Update initialize.py
mateuszanotto Oct 3, 2022
b2ca74d
Update main.py
mateuszanotto Oct 3, 2022
b6b1f04
Update forcefield.py
mateuszanotto Oct 3, 2022
7924090
Update main.py
mateuszanotto Oct 3, 2022
6812643
Update forcefield.py
mateuszanotto Oct 3, 2022
93ad3de
cosmetic changes
selimsami Oct 7, 2022
0a5c98a
missing f-string
selimsami Oct 7, 2022
ebb54fb
bump version
selimsami Oct 7, 2022
b4d2dc9
Update forcefield.py
mateuszanotto Oct 27, 2022
a9e8739
Update forcefield.py
mateuszanotto Apr 5, 2023
e844a73
Update forcefield.py
mateuszanotto Apr 13, 2023
4b3c4eb
Update forcefield.py
mateuszanotto Apr 18, 2023
a72c1c7
Update forcefield.py
mateuszanotto May 12, 2023
0da7dcf
Update forcefield.py
mateuszanotto May 12, 2023
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
294 changes: 205 additions & 89 deletions qforce/forcefield.py
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -1,9 +1,10 @@
import numpy as np
import string
#
from .elements import ATOM_SYM, ATOMMASS
from .molecule.non_bonded import calc_sigma_epsilon
from .forces import convert_to_inversion_rb
from .misc import LOGO_SEMICOL
from .misc import get_logo


class ForceField():
Expand All @@ -24,6 +25,11 @@ def __init__(self, job_name, config, mol, neighbors, exclude_all=[]):
self.exclusions = self.make_exclusions(mol.non_bonded, neighbors, exclude_all)
self.pairs = self.make_pairs(neighbors, mol.non_bonded)

if config.ff.output_software == 'gromacs':
self.write_ff = self.write_gromacs
elif config.ff.output_software == 'amber':
self.write_ff = self.write_amber

if self.polar:
self.polar_title = '_polar'
else:
Expand Down Expand Up @@ -75,7 +81,7 @@ def write_gro(self, directory, coords, alpha_map, box=[20., 20., 20.]):

def write_itp(self, mol, directory):
with open(f"{directory}/{self.mol_name}_qforce{self.polar_title}.itp", "w") as itp:
itp.write(LOGO_SEMICOL)
itp.write(get_logo(';'))
self.write_itp_atoms_and_molecule(itp, mol.non_bonded)
if self.polar:
self.write_itp_polarization(itp, mol.non_bonded)
Expand Down Expand Up @@ -128,7 +134,8 @@ def write_itp_atoms_and_molecule(self, itp, non_bonded):
itp.write("; name at_num mass charge type sigma epsilon\n")

for lj_type, lj_params in gro_atomtypes.items():
itp.write(f'{lj_type:>8} {non_bonded.lj_atomic_number[lj_type]:>7} {0:>8.4f} {0:>8.4f} {"A":>5} ')
itp.write(f'{lj_type:>8} {non_bonded.lj_atomic_number[lj_type]:>7} {0:>8.4f} '
f'{0:>8.4f} {"A":>5} ')
itp.write(f'{lj_params[0]:>12.5e} {lj_params[1]:>12.5e}\n')

if self.polar:
Expand Down Expand Up @@ -387,89 +394,198 @@ def set_charge(self, non_bonded):
q[list(non_bonded.alpha_map.keys())] += 8
return q

# bohr2nm = 0.052917721067
# if polar:
# alphas = qm.alpha*bohr2nm**3
# drude = {}
# n_drude = 1
# ff.atom_types.append(["DP", 0, 0, "S", 0, 0])

# for i, alpha in enumerate(alphas):
# if alpha > 0:
# drude[i] = mol.topo.n_atoms+n_drude
# ff.atoms[i][6] += 8
# # drude atoms
# ff.atoms.append([drude[i], 'DP', 2, 'MOL', f'D{atoms[i]}',
# i+1, -8., 0.])
# ff.coords.append(ff.coords[i])
# # polarizability
# ff.polar.append([i+1, drude[i], 1, alpha])
# n_drude += 1
# ff.natom = len(ff.atoms)
# for i, alpha in enumerate(alphas):
# if alpha > 0:
# # exclusions for balancing the drude particles
# for j in (mol.topo.neighbors[self.n_excl-2][i] +
# mol.topo.neighbors[self.n_excl-1][i]):
# if alphas[j] > 0:
# ff.exclu[drude[i]-1].extend([drude[j]])
# for j in mol.topo.neighbors[self.n_excl-1][i]:
# ff.exclu[drude[i]-1].extend([j+1])
# ff.exclu[drude[i]-1].sort()
# # thole polarizability
# for neigh in [mol.topo.neighbors[n][i] for n in range(self.n_excl)]:
# for j in neigh:
# if i < j and alphas[j] > 0:
# ff.thole.append([i+1, drude[i], j+1, drude[j], "2", 2.6, alpha,
# alphas[j]])

# def read_itp(self, itp_file):
# with open(itp_file, "r") as itp:
# in_section = []
# bond_atoms, bond_r0, bond_k = [], [], []

# for line in itp:
# low_line = line.lower().strip().replace(" ", "")
# unsplit = line
# line = line.split()
# if low_line == "" or low_line[0] == ";":
# continue
# elif "[" in low_line and "]" in low_line:
# open_bra = low_line.index("[") + 1
# close_bra = low_line.index("]")
# in_section = low_line[open_bra:close_bra]
# elif in_section == "atomtypes":
# self.atom_types.append([line[0], float(line[1]),
# float(line[2]), line[3],
# float(line[4]), float(line[5])])
# self.atype.append(line[0])
# self.c6.append(line[4])
# self.c12.append(line[5])
# elif in_section == "moleculetype":
# self.mol_type = line[0]
# elif in_section == "atoms":
# self.atoms.append([int(line[0]), line[1], int(line[2]),
# line[3], line[4], line[5], float(line[6]),
# float(line[7])])
# self.natom += 1
# elif in_section == "bonds":
# bond_atoms = (line[0:2])
# bond_r0 = (float(line[3]) * 10)
# bond_k = (float(line[4]) / 100)
# self.bond.append([bond_atoms, bond_r0, bond_k])
# self.bonds.append(unsplit)
# elif in_section == "angles":
# angle_atoms = (line[0:3])
# angle_theta0 = float(line[4])
# angle_k = float(line[5])
# self.angle.append([angle_atoms, angle_theta0, angle_k])
# self.angles.append(unsplit)
# elif in_section == "dihedrals":
# self.dihedrals.append(unsplit)
# if line[4] == "3":
# self.rbdihed.append([line[0:4], line[4:]])
# elif line[4] == "2":
# self.idihed.append([line[0:4], line[4:]])

# elif in_section == "pairs":
# self.pairs.append(sorted([int(line[0]), int(line[1])]))
def write_amber(self, directory, mol, coords):
if self.urey == True:
print(f'NOTE: AMBERs parameters file does not support Urey-Bradley terms.')
print(f' Set urey to False! \n')
j2cal = 0.2390057361
self.write_mol2(directory, mol, coords)
self.write_frcmod(directory, mol, coords, j2cal)

def write_mol2(self, directory, mol, coords):
with open(f"{directory}/{self.mol_name}_qforce{self.polar_title}.mol2", "w") as mol2:
self.write_mol2_title(mol2)
self.write_mol2_molecule(mol2, mol.topo, mol.terms)
self.write_mol2_atom(mol2, mol.topo, coords, mol.non_bonded)
self.write_mol2_bond(mol2, mol.topo, mol.terms)

def write_frcmod(self, directory, mol, coords, j2cal):
with open(f"{directory}/{self.mol_name}_qforce{self.polar_title}.frcmod", "w") as frcmod:
self.write_frcmod_mass(frcmod, mol.non_bonded)
self.write_frcmod_bonds(frcmod, mol.terms, mol.non_bonded.alpha_map, j2cal)
self.write_frcmod_angles(frcmod, mol.terms, j2cal)
self.write_frcmod_dihedrals(frcmod, mol.terms, j2cal)
self.write_frcmod_nonbond(frcmod, mol.non_bonded)

def write_mol2_title(self, mol2):
mol2.write(get_logo('#'))
mol2.write(f"# Name: {self.mol_name}\n")

def write_mol2_molecule(self, mol2, topo, terms):
mol2.write(f"@<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("esp") # type of charge
mol2.write("\n\n")

def write_mol2_atom(self, mol2, topo, coords, non_bonded):
mol2.write("@<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}"
f"{coords[i_idx-1][1]:10.4f} {coords[i_idx-1][2]:10.4f}")
mol2.write(f" {i_idx} {1} {self.mol_name[0:3]} {q:10.6f}\n")

def write_mol2_bond(self, mol2, topo, terms):
n_bonds = 1
mol2.write("@<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("@<TRIPOS>SUBSTRUCTURE\n")
mol2.write(f"{1:5} {self.mol_name[0:3]} {1:8} TEMP {0:>8d} **** **** {0:5} ROOT\n")

def write_frcmod_mass(self, frcmod, non_bonded):
frcmod.write(f'{self.mol_name} - frcmod generated by QForce\n')
frcmod.write('MASS\n')
for i, mass in enumerate((self.masses), start=1):
frcmod.write(f"{i} {mass}\n")

def write_frcmod_bonds(self, frcmod, terms, alpha_map, j2cal):
frcmod.write("\nBOND\n")
for bond in terms['bond']:
ids = bond.atomids + 1
fconst = bond.fconst * j2cal/2 # kJ -> kcal
equ = bond.equ
frcmod.write(f"{ids[0]:<2}-{ids[1]:<2} {fconst:>10.2f} {equ:>10.2f}\n")

for angle in terms['angle']:
ids = angle.atomids + 1
if self.urey:
urey = [term for term in terms['urey'] if np.array_equal(term.atomids,
angle.atomids)]
if not self.urey or len(urey) == 0:
None
else:
urey_equ = urey[0].equ
urey_fconst = urey[0].fconst * j2cal/2
frcmod.write(f"{ids[0]:<2}-{ids[2]:<2}{urey_fconst:>10.2f}{urey_equ:>10.2f}\n")


def write_frcmod_angles(self, frcmod, terms, j2cal):
frcmod.write("\nANGLE\n")
for angle in terms['angle']:
ids = angle.atomids + 1
fconst = angle.fconst * j2cal/2
equ = np.degrees(angle.equ)

if self.urey:
urey = [term for term in terms['urey'] if np.array_equal(term.atomids,
angle.atomids)]
if not self.urey or len(urey) == 0:
frcmod.write(f"{ids[0]:<2}-{ids[1]:<2}-{ids[2]:<2}")
frcmod.write(f"{fconst:>10.2f}{equ:>10.2f} \n")
else:
frcmod.write(f"{ids[0]:<2}-{ids[1]:<2}-{ids[2]:<2}")
frcmod.write(f"{fconst:>10.2f}{equ:>10.2f} \n")


def write_frcmod_dihedrals(self, frcmod, terms, j2cal):
if len(terms['dihedral']) > 0:
frcmod.write("\nDIHE\n")

# rigid dihedrals
if len(terms['dihedral/rigid']) > 0:

for dihed in terms['dihedral/rigid']:
ids = dihed.atomids + 1
equ = dihed.equ

if dihed.equ == 0.00 or dihed.equ == 3.141592653589793:
equ = 180.0

fconst = dihed.fconst * j2cal

fconstAmb = (2*fconst/(2**2))
frcmod.write(f"{ids[0]:<2}-{ids[1]:<2}-{ids[2]:<2}-{ids[3]:<2}")
frcmod.write(f" 2 {fconstAmb:>15.2f} {equ:>15.2f} 2\n")

if len(terms['dihedral/flexible']) > 0:

# flexible dihedrals
for dihed in terms['dihedral/flexible']:
ids = dihed.atomids + 1
c = dihed.equ * j2cal
fc = [1.0*c[0] + 0.5*c[2] + 0.3750*c[4],
1.0*c[1] + 0.75*c[3] + 0.6250*c[5],
0.5*c[2] + 0.5*c[4],
0.25*c[3] + 0.3125*c[5],
0.125*c[4],
0.0625*c[5]]

for n in range(3, 0, -1):
if n % 2 == 1:
equ = 180.0
else:
equ = 0.00
frcmod.write(f"{ids[0]:<2}-{ids[1]:<2}-{ids[2]:<2}-{ids[3]:<2}")
frcmod.write(f" 2 {fc[n]:>15.2f} {equ:>15.2f} -{n+1}\n")

# Last term
frcmod.write(f"{ids[0]:<2}-{ids[1]:<2}-{ids[2]:<2}-{ids[3]:<2}")
frcmod.write(f" 2 {fc[0]:>15.2f} {0.00:>15.2f} 1\n")

# inversion dihedrals
if len(terms['dihedral/inversion']) > 0:

for dihed in terms['dihedral/inversion']:
ids = dihed.atomids + 1
c0, c1, c2 = convert_to_inversion_rb(dihed.fconst, dihed.equ)

fc = [(1.0*c0/2 + 0.5*c2/2)* j2cal,
1.0*c1/2* j2cal,
0.5*c2/2* j2cal]

for n in range(3, 0, -1):
if n % 2 == 1:
equ = 180.0
else:
equ = 0.00
frcmod.write(f"{ids[0]:<2}-{ids[1]:<2}-{ids[2]:<2}-{ids[3]:<2}")
frcmod.write(f" 1 {fc[n]:>15.2f} {equ:>15.2f} -{n+1}\n")

# Last term
frcmod.write(f"{ids[0]:<2}-{ids[1]:<2}-{ids[2]:<2}-{ids[3]:<2}")
frcmod.write(f" 1 {fc[0]:>15.2f} {0.00:>15.2f} 1\n")


# improper dihedrals
Copy link
Copy Markdown
Contributor

@xiki-tempula xiki-tempula Dec 21, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Ok, this part is indented incorrectly

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)
Copy link
Copy Markdown
Contributor

@xiki-tempula xiki-tempula Dec 21, 2023

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Amber only accept 180 as the equilibrium value for improper
Warning: Expected Improper Torsion PHASE=180 (0.000000)

fconst = dihed.fconst * j2cal/2
frcmod.write(f"{ids[0]:<2}-{ids[1]:<2}-{ids[2]:<2}-{ids[3]:<2}")
frcmod.write(f" 1 {fconst:>15.2f} {equ:>15.2f} 2\n")

def write_frcmod_nonbond(self, frcmod, non_bonded):
gro_atomtypes, gro_nonbonded, gro_1_4 = self.convert_to_gromacs_nonbonded(non_bonded)
atom_ids = {} # original atom types

for i, (lj_type, a_name) in enumerate(zip(non_bonded.lj_types, self.atom_names), start=1):
atom_ids[i] = lj_type

frcmod.write("\nNONB\n")
# support to gaff and gaff2
if self.comb_rule == 2:
for i in range(1, self.n_atoms+1, 1):
frcmod.write(f"{i:<8}")
frcmod.write(f"{gro_atomtypes[atom_ids[i]][0]*5.612:>12.4f}")
frcmod.write(f"{gro_atomtypes[atom_ids[i]][1]/4.184:>12.4f}")
frcmod.write(f"\t{atom_ids[i]}\n")
8 changes: 6 additions & 2 deletions qforce/initialize.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -9,12 +9,16 @@
from .qm.qm import QM, implemented_qm_software
from .molecule.terms import Terms
from .dihedral_scan import DihedralScan
from .misc import LOGO
from .misc import get_logo


class Initialize(Colt):
_user_input = """
[ff]
# MD software to use
# Currently amber supports only gaff2 atom types
output_software = gromacs :: str :: gromacs, amber

# Number of n equivalent neighbors needed to consider two atoms equivalent
# Negative values turns off equivalence, 0 makes same elements equivalent
n_equiv = 4 :: int
Expand Down Expand Up @@ -164,7 +168,7 @@ def _check_and_copy_settings_file(job_dir, config_file):


def initialize(filename, config_file, presets=None):
print(LOGO)
print(get_logo())

job_info = _get_job_info(filename)
settings_file = _check_and_copy_settings_file(job_info.dir, config_file)
Expand Down
Loading