Mercurial > repos > chemteam > parmconv
diff parmconv.py @ 3:2b82fc7bec67 draft
"planemo upload for repository https://github.com/galaxycomputationalchemistry/galaxy-tools-compchem/ commit cc13bd32ef2d80b01dc197a3ca120a4ff9f0dacc"
author | chemteam |
---|---|
date | Fri, 12 Mar 2021 12:32:24 +0000 |
parents | |
children |
line wrap: on
line diff
--- /dev/null Thu Jan 01 00:00:00 1970 +0000 +++ b/parmconv.py Fri Mar 12 12:32:24 2021 +0000 @@ -0,0 +1,82 @@ +import argparse +import io +import sys +from contextlib import redirect_stdout + +import parmed +from parmed import amber, gromacs +from parmed.tools.changeradii import ChRad + + +def parse_command_line(argv): + parser = argparse.ArgumentParser() + parser.add_argument('--istr', help='input structure', required=True) + parser.add_argument('--itop', help='input topology file', required=True) + parser.add_argument('--istripmask', help='stripmask') + parser.add_argument('--iradii', required=True, help='parmed radii are \ + GB_RADII amber6,bondi, mbondi, mbondi2, mbondi3') + parser.add_argument('--removedihe', action='store_true', + default=False, help='remove dihedrals with zero \ + periodicity') + parser.add_argument('--removebox', action='store_true', + default=False, help='remove periodic box info') + parser.add_argument('--o_prmtop', help='AMBER output topology', + required=True) + return parser.parse_args() + + +def get_ids(dihedrals): + """ + goes through dihedrals and looks for any with per=0. + returns a reverse sorted list of ids to be removed. + """ + indices = [] + for k, v in enumerate(dihedrals): + f = io.StringIO() + with redirect_stdout(f): + print(v) + if f.getvalue().find("per=0") != -1: + indices.append(k) + indices.sort(reverse=True) + return indices + + +args = parse_command_line(sys.argv) + +gmx_top = gromacs.GromacsTopologyFile(args.itop) +gmx_gro = gromacs.GromacsGroFile.parse(args.istr) + +if not args.removebox: + # keep box info + gmx_top.box = gmx_gro.box + gmx_top.positions = gmx_gro.positions + + +if args.removedihe: + ids_to_remove = get_ids(gmx_top.dihedrals) + print("Original number of dihedrals %i" % len(gmx_top.dihedrals)) + for i in ids_to_remove: + gmx_top.dihedrals.pop(i) + print("Update number of dihedrals %i" % len(gmx_top.dihedrals)) + +if args.istripmask is not None: + if args.istripmask == "": + pass + else: + gmx_top.strip(args.istripmask) + +radii = str(args.iradii) +parmed.tools.changeRadii(gmx_top, radii) +amb_prm = amber.AmberParm.from_structure(gmx_top) +parmed.tools.changeRadii(amb_prm, radii) + +if args.removebox: + amb_prm.pointers['IFBOX'] = 0 + +ChRad(amb_prm, radii) +for i, atom in enumerate(amb_prm.atoms): + amb_prm.parm_data['RADII'][i] = atom.solvent_radius + amb_prm.parm_data['SCREEN'][i] = atom.screen + + +amb_prm.write_parm(args.o_prmtop)