#!/usr/bin/env python
#* **************************************************************************
#*
#* CDTK, Chemical Dynamics Toolkit
#* A modular system for chemical dynamics applications and more
#*
#* Copyright (C) 2011, 2012, 2013, 2014, 2015, 2016
#* Oriol Vendrell, DESY, <oriol.vendrell@desy.de>
#*
#* Copyright (C) 2017, 2018, 2019
#* Ralph Welsch, DESY, <ralph.welsch@desy.de>
#*
#* Copyright (C) 2020, 2021, 2022, 2023
#* Ludger Inhester, DESY, ludger.inhester@cfel.de
#*
#* This file is part of CDTK.
#*
#* CDTK is free software: you can redistribute it and/or modify
#* it under the terms of the GNU General Public License as published by
#* the Free Software Foundation, either version 3 of the License, or
#* (at your option) any later version.
#*
#* This program is distributed in the hope that it will be useful,
#* but WITHOUT ANY WARRANTY; without even the implied warranty of
#* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
#* GNU General Public License for more details.
#*
#* You should have received a copy of the GNU General Public License
#* along with this program. If not, see <http://www.gnu.org/licenses/>.
#*
#* **************************************************************************
"""
This module helps to transform gromacs topologies to cdtk forcefield classes
Not all possible topology options of a gromacs topology can be handled, specifically no '#input' options are treated.
"""
from optparse import OptionParser
import numpy as np
import CDTK.Tools.Conversion as cv
nm2bohr = cv.an2au * 10 # 18.8973
nm_ps2bohr_au = cv.an2au * 10 / (cv.fs2au * 1000) # 4.571038e-04
[docs]
def get_arg():
parser = OptionParser()
parser = OptionParser(usage="usage: %prog [options] filename")
parser.add_option(
"-t",
"--top_file",
dest="gmxtopFile",
type="str",
default="topol.top",
help="GMX .top file",
)
return parser
[docs]
def parse_top(filename):
def process_defaults(ffparam, line):
fields = line.split()
if len(fields) != 5:
raise ValueError("[default] must have 5 fields")
nbfunc, combrule, genpairs, fudgeLJ, fudgeQQ = fields
fudgeLJ = float(fudgeLJ)
fudgeQQ = float(fudgeQQ)
nbfunc = int(nbfunc)
combrule = int(combrule)
ffparam["fudgeLJ"] = fudgeLJ
ffparam["fudgeQQ"] = fudgeQQ
ffparam["nbfunc"] = nbfunc
if combrule == 1:
ffparam["combrule"] = "geometric"
elif combrule == 2:
ffparam["combrule"] = "Lorentz-Berthelot"
elif combrule == 3:
ffparam["combrule"] = "OPLS"
else:
raise ValueError("unknown combination rule")
def process_atomtypes(ffparam, line):
fields = line.split()
if len(fields) == 8:
atomtype, bondedtype, atomnum, mass, charge, parttype, v, w = fields
elif len(fields) == 7:
atomtype, bondedtype, mass, charge, parttype, v, w = fields
atomnum = 0
elif len(fields) == 6:
atomnum, mass, charge, parttype, v, w = fields
# atomtype = a for a in cv.periodicTable[a]["atomic_number"] == int(atomnum)
raise ValueError("[atomtypes] with 6 fields not supported")
else:
raise ValueError("[atomtypes] must have 6 or 8 fields")
if not "atomtypes" in ffparam.keys():
ffparam["atomtypes"] = {}
if ffparam["combrule"] == "geometric":
ffparam["atomtypes"][atomtype] = {
"bondedtype": bondedtype,
"atomnum": int(atomnum),
"mass": float(mass),
"charge": float(charge),
"parttype": parttype,
"C6": float(v),
"C12": float(w),
}
else:
ffparam["atomtypes"][atomtype] = {
"bondedtype": bondedtype,
"atomnum": int(atomnum),
"mass": float(mass),
"charge": float(charge),
"parttype": parttype,
"sigma": float(v),
"epsilon": float(w),
}
def process_bondtypes(ffparam, line):
fields = line.split()
if len(fields) != 5:
raise ValueError("[bondtypes] must have 5 fields")
ai, aj, func, b, k = fields
if not "bondtypes" in ffparam.keys():
ffparam["bondtypes"] = []
ffparam["bondtypes"].append({
"ti": ai,
"tj": aj,
"func": int(func),
"b": float(b),
"K": float(k),
})
def process_moleculetype(ffparam, line):
fields = line.split()
if len(fields) != 2:
raise ValueError("[moleculetype] must have 2 fields")
if not "moltype" in ffparam.keys():
ffparam["moltype"] = []
name, nexcl = fields
# ffparam["moltype"].append({"name": name, "nexcl": int(nexcl)})
ffparam["nexcl"] = int(nexcl)
def process_atoms(ffparam, line):
fields = line.split()
if len(fields) != 8:
raise ValueError("[atoms] must have 7 fields")
atomnumber, atomtype, resnum, resname, atomname, chargegroup, charge, mass = fields
atomname = atomname.replace('_', '')
if not "atomnames" in ffparam.keys():
ffparam["atomnames"] = []
if atomname in ffparam["atomnames"]:
print(
f"# Warning: atomnames are not unique; Please use unique atom names\n# relabelling {atomname} to {atomname + atomnumber} ")
atomname = atomname + atomnumber
# raise ValueError(f"atom {atomname} already in atoms; please choose unique atomnames")
ffparam["atomnames"].append(atomname)
if not "atomtype" in ffparam.keys():
ffparam["atomtype"] = []
ffparam["atomtype"].append(atomtype)
if not "charges" in ffparam.keys():
ffparam["charges"] = {}
if not "masses" in ffparam.keys():
ffparam["masses"] = {}
if ffparam["combrule"] == "geometric":
if not "C6" in ffparam.keys():
ffparam["C6"] = {}
if not "C12" in ffparam.keys():
ffparam["C12"] = {}
else:
if not "epsilon" in ffparam.keys():
ffparam["epsilon"] = {}
if not "sigma" in ffparam.keys():
ffparam["sigma"] = {}
ffparam["charges"][atomname] = float(charge)
ffparam["masses"][atomname] = float(mass)
if ffparam["combrule"] == "geometric":
C6 = ffparam["atomtypes"][atomtype]["C6"]
C12 = ffparam["atomtypes"][atomtype]["C12"]
ffparam["C6"][atomname] = C6
ffparam["C12"][atomname] = C12
else:
epsilon = ffparam["atomtypes"][atomtype]["epsilon"]
sigma = ffparam["atomtypes"][atomtype]["sigma"]
ffparam["epsilon"][atomname] = epsilon
ffparam["sigma"][atomname] = sigma
def process_bonds(ffparam, line):
fields = line.split()
if not "bonds" in ffparam.keys():
ffparam["bonds"] = []
if len(fields) == 2:
ai, aj = fields
ati = ffparam["atomtype"][int(ai)-1]
atj = ffparam["atomtype"][int(aj)-1]
bti = ffparam["atomtypes"][ati]["bondedtype"]
btj = ffparam["atomtypes"][atj]["bondedtype"]
bondparameters = [bond for bond in ffparam["bondtypes"]
if bond["ti"] == bti and bond["tj"] == btj
or bond["tj"] == bti and bond["ti"] == btj]
if (len(bondparameters) != 1):
raise ValueError(
f" found {len(bondparameters)} matching for bondtype {bti} {btj}")
bondparameters = bondparameters[0]
ffparam["bonds"].append({
"ai": ffparam["atomnames"][int(ai)-1],
"aj": ffparam["atomnames"][int(aj)-1],
"func": bondparameters["func"],
"b": bondparameters["b"],
"K": bondparameters["K"],
})
elif len(fields) == 5:
ai, aj, func, b, k = fields
ffparam["bonds"].append({
"ai": ffparam["atomnames"][int(ai)-1],
"aj": ffparam["atomnames"][int(aj)-1],
"func": int(func),
"b": float(b),
"K": float(k),
})
elif len(fields) != 5:
raise ValueError("[bonds] must have 5 fields")
def process_pairs(ffparam, line):
return
def process_angles(ffparam, line):
fields = line.split()
if not "angles" in ffparam.keys():
ffparam["angles"] = []
if len(fields) == 3:
ai, aj, ak = fields
ati = ffparam["atomtype"][int(ai)-1]
atj = ffparam["atomtype"][int(aj)-1]
atk = ffparam["atomtype"][int(ak)-1]
bti = ffparam["atomtypes"][ati]["bondedtype"]
btj = ffparam["atomtypes"][atj]["bondedtype"]
btk = ffparam["atomtypes"][atk]["bondedtype"]
angleparameters = [angle for angle in ffparam["angletypes"]
if (angle["ti"] == bti and angle["tj"] == btj and angle["tk"] == btk)
or (angle["tk"] == bti and angle["tj"] == btj and angle["ti"] == btk)]
if (len(angleparameters) != 1):
raise ValueError(
f" found {len(angleparameters)} matching for bondtype {bti} {btj}")
angleparameters = angleparameters[0]
ffparam["angles"].append({
"ai": ffparam["atomnames"][int(ai)-1],
"aj": ffparam["atomnames"][int(aj)-1],
"ak": ffparam["atomnames"][int(ak)-1],
"func": angleparameters["func"],
"theta": angleparameters["theta"],
"K": angleparameters["K"],
})
elif len(fields) == 6:
ai, aj, ak, func, theta, k = fields
ffparam["angles"].append({
"ai": ffparam["atomnames"][int(ai)-1],
"aj": ffparam["atomnames"][int(aj)-1],
"ak": ffparam["atomnames"][int(ak)-1],
"func": int(func),
"theta": float(theta),
"K": float(k),
})
else:
raise ValueError("[angles] must have 6 fields")
def process_angletypes(ffparam, line):
fields = line.split()
if len(fields) != 6:
raise ValueError("[angletypes] must have 6 fields")
ai, aj, ak, func, theta, k = fields
if not "angletypes" in ffparam.keys():
ffparam["angletypes"] = []
ffparam["angletypes"].append({
"ti": ai,
"tj": aj,
"tk": ak,
"func": int(func),
"theta": float(theta),
"K": float(k),
})
def process_dihedrals(ffparam, line):
fields = line.split()
if not "dihedrals" in ffparam.keys():
ffparam["dihedrals"] = []
if len(fields) == 4:
ai, aj, ak, al = fields
ati = ffparam["atomtype"][int(ai)-1]
atj = ffparam["atomtype"][int(aj)-1]
atk = ffparam["atomtype"][int(ak)-1]
atl = ffparam["atomtype"][int(al)-1]
bti = ffparam["atomtypes"][ati]["bondedtype"]
btj = ffparam["atomtypes"][atj]["bondedtype"]
btk = ffparam["atomtypes"][atk]["bondedtype"]
btl = ffparam["atomtypes"][atl]["bondedtype"]
dihedralparameters = [angle for angle in ffparam["dihedraltypes"]
if (angle["ti"] == bti and angle["tj"] == btj and angle["tk"] == btk and angle["tl"] == btl)
or (angle["tk"] == bti and angle["tl"] == btj and angle["ti"] == btk and angle["tj"] == btl)
or (angle["tj"] == bti and angle["ti"] == btj and angle["tl"] == btk and angle["tk"] == btl)
or (angle["tl"] == bti and angle["tk"] == btj and angle["tj"] == btk and angle["ti"] == btl)
]
if (len(dihedralparameters) != 1):
raise ValueError(
f" found {len(dihedralparameters)} matching for dihedraltype {bti} {btj} {btk} {btl}")
dihedralparameters = dihedralparameters[0]
ffparam["dihedrals"].append({
"ai": ffparam["atomnames"][int(ai)-1],
"aj": ffparam["atomnames"][int(aj)-1],
"ak": ffparam["atomnames"][int(ak)-1],
"al": ffparam["atomnames"][int(al)-1],
"func": dihedralparameters["func"],
"phiS": dihedralparameters["phiS"],
"K": dihedralparameters["K"],
"n": dihedralparameters["n"],
})
elif len(fields) == 8:
ai, aj, ak, al, func, phi, k, periodicity = fields
ffparam["dihedrals"].append({
"ai": ffparam["atomnames"][int(ai)-1],
"aj": ffparam["atomnames"][int(aj)-1],
"ak": ffparam["atomnames"][int(ak)-1],
"al": ffparam["atomnames"][int(al)-1],
"func": int(func),
"phiS": float(phi),
"K": float(k),
"n": int(periodicity),
})
elif len(fields) == 11:
ai, aj, ak, al, func, C0, C1, C2, C3, C4, C5 = fields
ffparam["dihedrals"].append({
"ai": ffparam["atomnames"][int(ai)-1],
"aj": ffparam["atomnames"][int(aj)-1],
"ak": ffparam["atomnames"][int(ak)-1],
"al": ffparam["atomnames"][int(al)-1],
"func": int(func),
"C0": float(C0),
"C1": float(C1),
"C2": float(C2),
"C3": float(C3),
"C4": float(C4),
"C5": float(C5),
})
else:
raise ValueError(
f"[dihedrals] must have 8 or 13 fields; this line has {len(fields)}")
def process_dihedraltypes(ffparam, line):
fields = line.split()
if not "dihedraltypes" in ffparam.keys():
ffparam["dihedraltypes"] = []
if len(fields) == 8:
ai, aj, ak, al, func, phi, k, periodicity = fields
ffparam["dihedraltypes"].append({
"ti": ai,
"tj": aj,
"tk": ak,
"tl": al,
"func": int(func),
"phiS": float(phi),
"K": float(k),
"n": int(periodicity),
})
elif len(fields) == 11:
ai, aj, ak, al, func, C0, C1, C2, C3, C4, C5 = fields
ffparam["dihedraltypes"].append({
"ti": ai,
"tj": aj,
"tk": ak,
"tl": al,
"func": int(func),
"C0": float(C0),
"C1": float(C1),
"C2": float(C2),
"C3": float(C3),
"C4": float(C4),
"C5": float(C5),
})
else:
raise ValueError(
f"[dihedrals] must have 8/13 fields; this line has {len(fields)}")
def process_system(ffparam, line):
return
def process_molecules(ffparam, line):
return
directives = {"defaults": process_defaults,
"atomtypes": process_atomtypes,
"moleculetype": process_moleculetype,
"atoms": process_atoms,
"bonds": process_bonds,
"angles": process_angles,
"pairs": process_pairs,
"dihedrals": process_dihedrals,
"system": process_system,
"molecules": process_molecules,
"bondtypes": process_bondtypes,
"angletypes": process_angletypes,
"dihedraltypes": process_dihedraltypes,
}
ffparam = {}
with open(filename, "r") as f:
ffparam["name"] = f"parsed from file {filename}"
for line in f:
line = line.strip()
if ";" in line:
line = line[: line.index(';')].strip()
if len(line) == 0:
continue
if line.startswith("#"):
raise NotImplementedError(
" preprocess directive not implemented")
elif line.startswith("["):
current_directive = line.strip("[] ")
continue
if current_directive in directives:
directives[current_directive](ffparam, line)
else:
raise ValueError(
f" section [{current_directive}] not implemented")
return ffparam
[docs]
def start():
args, _ = get_arg().parse_args()
ffparam = parse_top(args.gmxtopFile)
del ffparam["atomtypes"]
del ffparam["atomnames"]
if "bondtypes" in ffparam:
del ffparam["bondtypes"]
if "angletypes" in ffparam:
del ffparam["angletypes"]
if "dihedraltypes" in ffparam:
del ffparam["dihedraltypes"]
import pprint
pp = pprint.PrettyPrinter(
indent=2, sort_dicts=False, width=200, compact=False)
pp.pprint(ffparam)
if __name__ == "__main__":
start()