OpenMD 3.2
Molecular Dynamics in the Open
Loading...
Searching...
No Matches
solvator
1#!@Python3_EXECUTABLE@
2
3# program that reads in a water box and solute xyz (or pdb), merges
4# the two systems, and deletes overlapping molecules
5
6__author__ = "Chris Fennell and Dan Gezelter"
7__copyright__ = "Copyright (c) 2004-present The University of Notre Dame. All Rights Reserved."
8__license__ = "OpenMD"
9
10import argparse
11import sys
12import math
13
14d2tolerance = 7.5625 # distance to start cutting
15fileName = 'solvatedSystem.omd'
16startSnap = 0
17startFrame = 0
18startStunts = 0
19startMeta = 0
20soluteName = 'SOLUTE'
21nSolvent = 0
22
23# Argument parser setup
24parser = argparse.ArgumentParser(description='Carves a solute void in an OpenMD water box',
25 formatter_class=argparse.RawDescriptionHelpFormatter,
26 epilog=
27'''
28Example:
29 solvator -i myWater.omd -p mySolute.pdb -o mySystem.omd
30''')
31parser.add_argument('-f', '--flexible=', type=bool, dest='f',
32 action=argparse.BooleanOptionalAction, default=False,
33 help='include a flexible solute model description in the output file rather than the rigid body solute model description')
34parser.add_argument('-d', action='store', dest='d',
35 type=float, default=2.75, required=False,
36 help='Overlap removal distance (in angstroms)')
37parser.add_argument('-i', type=str, required=True, action='store', dest='i',
38 help='Solvent input file (OpenMD .omd file)')
39parser.add_argument('-n', type=str, help='Name for the solute',
40 required=False, dest='n', action='store',
41 default='SOLUTE')
42parser.add_argument('-o', '--output=', action='store', dest='o',
43 required=False, default='solvatedSystem.omd',
44 help='Carved solvent output file (OpenMD .omd format)')
45parser.add_argument('-p', type=str, help='Solute input file (pdb)',
46 action='store', dest='p')
47parser.add_argument('-x', type=str, help='Solute input file (xyz)',
48 action='store', dest='x')
49args = parser.parse_args()
50
51# Set variables from arguments
52if args.i:
53 solventName = args.i
54else:
55 sys.exit("Error: No solvent box specified\nPlease select a solvent box using the -i flag")
56
57soluteFileName = args.p if args.p else args.x
58fileName = args.o if args.o else fileName
59soluteName = args.n if args.n else soluteName
60
61if not (args.p or args.x):
62 sys.exit("Error: No solute file specified\nPlease select a solute file with the -p or -x flags (pdb or xyz respectively)")
63
64if args.d:
65 dval = args.d
66 d2tolerance = dval * dval
67else:
68 d2tolerance = 7.5625
69
70# Open files
71try:
72 with open(soluteFileName, 'r') as solute_file:
73 solute_data = solute_file.readlines()
74except FileNotFoundError:
75 sys.exit(f"Error: Can't open file {soluteFileName}")
76
77solute_names, solute_x, solute_y, solute_z = [], [], [], []
78
79if args.x:
80 headerLines = 2
81 for line in solute_data:
82 if headerLines > 0:
83 headerLines -= 1
84 else:
85 parts = line.split()
86 solute_names.append(parts[0])
87 solute_x.append(float(parts[1]))
88 solute_y.append(float(parts[2]))
89 solute_z.append(float(parts[3]))
90
91if args.p:
92 for line in solute_data:
93 parts = line.split()
94 if parts[0] == 'ATOM' or parts[0] == 'HETATM':
95 solute_names.append(parts[2])
96 solute_x.append(float(parts[5]))
97 solute_y.append(float(parts[6]))
98 solute_z.append(float(parts[7]))
99
100# Remap solute to the center of the box
101xSol = sum(solute_x) / len(solute_x)
102ySol = sum(solute_y) / len(solute_y)
103zSol = sum(solute_z) / len(solute_z)
104
105solute_x = [x - xSol for x in solute_x]
106solute_y = [y - ySol for y in solute_y]
107solute_z = [z - zSol for z in solute_z]
108
109soluteCount = len(solute_x) if args.f else 1
110solventCount = soluteCount
111
112# Process solvent file and write output
113try:
114 with open(solventName, 'r') as solvent_file, open(fileName, 'w') as out_file:
115 startSnap = startFrame = startStunts = startMeta = 0
116 frameData, metaLines, goodSolventMolecules = [], [], []
117
118 for line in solvent_file:
119 if "</Snapshot" in line:
120 startSnap = 0
121 if "</FrameData" in line:
122 startFrame = 0
123 if "</StuntDoubles" in line:
124 startStunts = 0
125 if "</MetaData" in line:
126 startMeta = 0
127 if "</OpenMD" in line:
128 break
129
130 if startMeta == 1:
131 metaLines.append(line)
132
133 if startSnap == 1:
134 if startFrame == 1:
135 frameData.append(line)
136 if "Hmat" in line:
137 parts = line.split()
138 hxx, hyy, hzz = float(parts[2][:-1]), float(parts[8][:-1]), float(parts[14])
139 if startStunts == 1:
140 parts = line.split()
141 x_val = float(parts[2]) - (hxx * round(float(parts[2]) / hxx))
142 y_val = float(parts[3]) - (hyy * round(float(parts[3]) / hyy))
143 z_val = float(parts[4]) - (hzz * round(float(parts[4]) / hzz))
144
145 saveFlag = True
146 for i in range(len(solute_x)):
147 diff_x = x_val - solute_x[i]
148 diff_y = y_val - solute_y[i]
149 diff_z = z_val - solute_z[i]
150 dist2 = diff_x * diff_x + diff_y * diff_y + diff_z * diff_z
151 if dist2 < d2tolerance:
152 saveFlag = False
153 break
154
155 if saveFlag:
156 saveLine = f"{solventCount}\t{parts[1]}\t{parts[2]}"
157 saveLine += ''.join([f" {p}" for p in parts[3:]])
158 goodSolventMolecules.append(saveLine)
159 solventCount += 1
160
161 if "<Snapshot" in line:
162 startSnap = 1
163 if "<FrameData" in line:
164 startFrame = 1
165 if "<StuntDoubles" in line:
166 startStunts = 1
167 if "<MetaData" in line:
168 startMeta = 1
169
170 nSolvent = len(goodSolventMolecules)
171
172 # Write output file
173 out_file.write("<OpenMD version=2>\n")
174 out_file.write(" <MetaData>\n")
175 out_file.write("\nmolecule{{\n name = \"{}\";\n\n".format(soluteName))
176
177 if args.f:
178 for i in range(len(solute_x)):
179 out_file.write(f" atom[{i}]{{\n type = \"{solute_names[i]}\";\n }}\n")
180 out_file.write("}\n")
181 else:
182 for i in range(len(solute_x)):
183 out_file.write(f" atom[{i}]{{\n type = \"{solute_names[i]}\";\n")
184 out_file.write(f" position({solute_x[i]}, {solute_y[i]}, {solute_z[i]});\n }}\n")
185 out_file.write("\n rigidBody[0]{\n members(")
186 out_file.write(", ".join(str(i) for i in range(len(solute_x))))
187 out_file.write(");\n }\n}\n")
188
189 out_file.write(f"\ncomponent{{\n type = \"{soluteName}\";\n nMol = 1;\n}}\n")
190 for line in metaLines:
191 if "nMol" in line:
192 out_file.write(f" nMol = {nSolvent};\n")
193 else:
194 out_file.write(line)
195
196 out_file.write(" </MetaData>\n")
197 out_file.write(" <Snapshot>\n")
198 out_file.write(" <FrameData>\n")
199
200 for line in frameData:
201 out_file.write(line)
202
203 out_file.write(" </FrameData>\n")
204 out_file.write(" <StuntDoubles>\n")
205
206 if args.f:
207 for i in range(len(solute_x)):
208 out_file.write(f"{i}\tp\t{solute_x[i]} {solute_y[i]} {solute_z[i]}\n")
209 else:
210 out_file.write("0\tpq\t0.0 0.0 0.0 1.0 0.0 0.0 0.0\n")
211
212 for molecule in goodSolventMolecules:
213 out_file.write(f"{molecule}\n")
214
215 out_file.write(" </StuntDoubles>\n")
216 out_file.write(" </Snapshot>\n")
217 out_file.write("</OpenMD>\n")
218
219 print(f"The solvated system \"{fileName}\" was generated.")
220
221except FileNotFoundError:
222 sys.exit(f"Error: Can't open file {solventName}")