Mercurial > repos > chemteam > parmconv

comparison parmconv.py @ 3:2b82fc7bec67 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
"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
comparison
equal deleted inserted replaced
2:b500cc25dd15 3:2b82fc7bec67
1 import argparse
2 import io
3 import sys
4 from contextlib import redirect_stdout
5
6 import parmed
7 from parmed import amber, gromacs
8 from parmed.tools.changeradii import ChRad
9
10
11 def parse_command_line(argv):
12 parser = argparse.ArgumentParser()
13 parser.add_argument('--istr', help='input structure', required=True)
14 parser.add_argument('--itop', help='input topology file', required=True)
15 parser.add_argument('--istripmask', help='stripmask')
16 parser.add_argument('--iradii', required=True, help='parmed radii are \
17 GB_RADII amber6,bondi, mbondi, mbondi2, mbondi3')
18 parser.add_argument('--removedihe', action='store_true',
19 default=False, help='remove dihedrals with zero \
20 periodicity')
21 parser.add_argument('--removebox', action='store_true',
22 default=False, help='remove periodic box info')
23 parser.add_argument('--o_prmtop', help='AMBER output topology',
24 required=True)
25 return parser.parse_args()
26
27
28 def get_ids(dihedrals):
29 """
30 goes through dihedrals and looks for any with per=0.
31 returns a reverse sorted list of ids to be removed.
32 """
33 indices = []
34 for k, v in enumerate(dihedrals):
35 f = io.StringIO()
36 with redirect_stdout(f):
37 print(v)
38 if f.getvalue().find("per=0") != -1:
39 indices.append(k)
40 indices.sort(reverse=True)
41 return indices
42
43
44 args = parse_command_line(sys.argv)
45
46 gmx_top = gromacs.GromacsTopologyFile(args.itop)
47 gmx_gro = gromacs.GromacsGroFile.parse(args.istr)
48
49 if not args.removebox:
50 # keep box info
51 gmx_top.box = gmx_gro.box
52 gmx_top.positions = gmx_gro.positions
53
54
55 if args.removedihe:
56 ids_to_remove = get_ids(gmx_top.dihedrals)
57 print("Original number of dihedrals %i" % len(gmx_top.dihedrals))
58 for i in ids_to_remove:
59 gmx_top.dihedrals.pop(i)
60 print("Update number of dihedrals %i" % len(gmx_top.dihedrals))
61
62 if args.istripmask is not None:
63 if args.istripmask == "":
64 pass
65 else:
66 gmx_top.strip(args.istripmask)
67
68 radii = str(args.iradii)
69 parmed.tools.changeRadii(gmx_top, radii)
70 amb_prm = amber.AmberParm.from_structure(gmx_top)
71 parmed.tools.changeRadii(amb_prm, radii)
72
73 if args.removebox:
74 amb_prm.pointers['IFBOX'] = 0
75
76 ChRad(amb_prm, radii)
77 for i, atom in enumerate(amb_prm.atoms):
78 amb_prm.parm_data['RADII'][i] = atom.solvent_radius
79 amb_prm.parm_data['SCREEN'][i] = atom.screen
80
81
82 amb_prm.write_parm(args.o_prmtop)