00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024 import _ost_mol_mm as mm
00025 import ost
00026 from ost import settings, mol, geom
00027 import os, subprocess, math
00028
00029
00030
00031 def _GetInteraction(functype, atoms, params):
00032 """Get an mm.Interaction with the given func-type and params for the given
00033 atoms (name and types extracted from there)."""
00034 interaction = mm.Interaction(functype)
00035 interaction.SetNames([a.name for a in atoms])
00036 interaction.SetTypes([a.GetStringProp('type') for a in atoms])
00037 interaction.SetParam(params)
00038 return interaction
00039
00040 def _MakeComponentBuildingBlock(eh, ff_dict):
00041 """Take EntityHandle eh (from ParseModifiedPDB) and ff_dict (from
00042 ParseAmberForceField) and return BuildingBlock."""
00043
00044 dist_scale = 1./10.0
00045 angle_scale = math.pi/180.
00046
00047 bond_k_scale = 418.4*2.
00048 angle_k_scale = 4.184
00049
00050
00051 atype_dict = {}
00052 for a in eh.atoms:
00053 atype = a.GetStringProp('type')
00054 if not atype in atype_dict:
00055 atype_dict[atype] = [a.handle]
00056 else:
00057 atype_dict[atype].append(a.handle)
00058
00059
00060 for atype, mass in ff_dict['MASS']:
00061 for a in atype_dict[atype]: a.SetMass(mass)
00062
00063
00064 bb = mm.BuildingBlock()
00065 for a in eh.atoms:
00066 bb.AddAtom(a.name, a.GetStringProp('type'), a.GetCharge(), a.GetMass())
00067
00068
00069 bl = eh.GetBondList()
00070 for b in bl:
00071 a1 = b.GetFirst()
00072 a2 = b.GetSecond()
00073 bond = mm.Interaction(mm.FuncType.HarmonicBond)
00074 bond.SetNames([a1.name, a2.name])
00075 bond.SetTypes([a1.GetStringProp('type'), a2.GetStringProp('type')])
00076 bb.AddBond(bond)
00077 added_bonds = []
00078 for atype1, atype2, d0, k in ff_dict['BOND']:
00079 for a1 in atype_dict[atype1]:
00080 for a2 in atype_dict[atype2]:
00081
00082 if not mol.BondExists(a1, a2): continue
00083 if [a1, a2] in added_bonds or [a2, a1] in added_bonds: continue
00084 added_bonds.append([a1,a2])
00085 params = [d0*dist_scale, k*bond_k_scale]
00086 bond = _GetInteraction(mm.FuncType.HarmonicBond, [a1, a2], params)
00087 bb.AddBond(bond, replace_existing=True)
00088
00089
00090 added_angles = []
00091 for atype1, atype2, atype3, a0, k in ff_dict['ANGL']:
00092
00093 for a2 in atype_dict[atype2]:
00094 for a1 in atype_dict[atype1]:
00095 if not mol.BondExists(a1, a2): continue
00096 for a3 in atype_dict[atype3]:
00097
00098 if not mol.BondExists(a2, a3): continue
00099 if a1 == a3: continue
00100 if [a1, a2, a3] in added_angles or [a3, a2, a1] in added_angles:
00101 continue
00102 added_angles.append([a1, a2, a3])
00103 angle = _GetInteraction(mm.FuncType.HarmonicAngle, [a1, a2, a3],
00104 [a0*angle_scale, k*angle_k_scale*2])
00105 bb.AddAngle(angle)
00106
00107
00108 for atype1, atype2, atype3, atype4, idiv, period, phase, k in ff_dict['DIHE']:
00109
00110 added_dihedrals = []
00111 for a1 in atype_dict[atype1]:
00112 for a2 in atype_dict[atype2]:
00113 if not mol.BondExists(a1, a2): continue
00114 for a3 in atype_dict[atype3]:
00115 if not mol.BondExists(a2, a3): continue
00116 if a1 == a3: continue
00117 for a4 in atype_dict[atype4]:
00118
00119 if not mol.BondExists(a3, a4): continue
00120 if a2 == a4: continue
00121 if [a1, a2, a3, a4] in added_dihedrals or \
00122 [a4, a3, a2, a1] in added_dihedrals: continue
00123 added_dihedrals.append([a1, a2, a3, a4])
00124 dihe = _GetInteraction(mm.FuncType.PeriodicDihedral, [a1, a2, a3, a4],
00125 [period, phase*angle_scale, k*angle_k_scale])
00126 bb.AddDihedral(dihe)
00127
00128
00129 added_impropers = []
00130 for atype1, atype2, atype3, atype4, period, phase, k in ff_dict['IMPR']:
00131
00132 for ac in atype_dict[atype3]:
00133 for a1 in atype_dict[atype1]:
00134 if not mol.BondExists(ac, a1): continue
00135 for a2 in atype_dict[atype2]:
00136 if not mol.BondExists(ac, a2): continue
00137 if a1 == a2: continue
00138 for a4 in atype_dict[atype4]:
00139
00140 if not mol.BondExists(ac, a4): continue
00141 if a2 == a4 or a1 == a4: continue
00142 if [ac, a1, a2, a4] in added_impropers or \
00143 [ac, a1, a4, a2] in added_impropers or \
00144 [ac, a2, a1, a4] in added_impropers or \
00145 [ac, a2, a4, a1] in added_impropers or \
00146 [ac, a4, a1, a2] in added_impropers or \
00147 [ac, a4, a2, a1] in added_impropers: continue
00148 added_impropers.append([ac, a1, a2, a4])
00149 impr = _GetInteraction(mm.FuncType.PeriodicImproper, [a1, a2, ac, a4],
00150 [period, phase*angle_scale, k*angle_k_scale])
00151 bb.AddImproper(impr)
00152 return bb
00153
00154 def _ParseModifiedPDB(filename):
00155 """Read mpdb file produced by antechamber and return tuple of:
00156 - EntityHandle with connectivity, atom types (property 'type') and charges
00157 - Residue name as extracted from the mpdb file
00158 A RuntimeError is raised if the file can contains multiple residues.
00159 """
00160 eh = mol.CreateEntity()
00161 rname = ''
00162 edi = eh.EditXCS(mol.BUFFERED_EDIT)
00163 chain = edi.InsertChain('A')
00164 bond_list = []
00165
00166 with open(filename, 'r') as in_file:
00167 for line in in_file:
00168
00169
00170 if line.startswith('ATOM'):
00171 aname = line[12:17].strip()
00172
00173 if not rname:
00174 rname = line[17:20].strip()
00175 r = edi.AppendResidue(chain, rname)
00176 elif rname != line[17:20].strip():
00177 raise RuntimeError("More than one residue in file " + filename +\
00178 ". Cannot parse!")
00179
00180 charge = float(line[54:66])
00181 atype = line[78:80].strip()
00182 a = edi.InsertAtom(r, aname, geom.Vec3())
00183 a.SetStringProp('type', atype)
00184 a.SetCharge(charge)
00185 elif line.startswith('CONECT'):
00186 ai1 = int(line[6:11])
00187
00188 for i in range(4):
00189 try:
00190 j = 11 + 5*i
00191 ai2 = int(line[j:j+5])
00192
00193 s = set([ai1, ai2])
00194 if not s in bond_list: bond_list.append(s)
00195 except:
00196
00197
00198 continue
00199
00200 for indices in bond_list:
00201 indices = list(indices)
00202 a1 = r.atoms[indices[0]-1]
00203 a2 = r.atoms[indices[1]-1]
00204 edi.Connect(a1, a2)
00205
00206 edi.UpdateICS()
00207 return eh, rname
00208
00209 def _ParseAmberForceField(filename):
00210 """Read frcmod file produced by parmchk2 and return dictionary with all
00211 entries for masses, bonds, angles, dihedrals, impropers and non-bonded (LJ)
00212 interactions. Stored as key/list-of-value pairs:
00213 - 'MASS': [atype, mass]
00214 - 'BOND': [atype1, atype2, d0, k]
00215 - 'ANGL': [atype1, atype2, atype3, a0, k]
00216 - 'DIHE': [atype1, atype2, atype3, atype4, idiv, period, phase, k/idiv]
00217 - 'IMPR': [atype1, atype2, atype3, atype4, period, phase, k]
00218 - 'NONB': [Rvdw, epsilon]
00219 """
00220 keywords = ['MASS', 'BOND', 'ANGL', 'DIHE', 'IMPR', 'NONB']
00221 with open(filename, 'r') as in_file:
00222 ff_dict = {}
00223 for line in in_file:
00224
00225 keyword = line[:4]
00226 if not keyword in keywords: continue
00227
00228 ff_dict[keyword] = []
00229 line = in_file.next()
00230 while len(line.strip()) > 0:
00231
00232 if 'ATTN' in line:
00233 ost.LogWarning('The following line in ' + filename + ' (' + keyword +\
00234 ' section) needs revision:\n' + line.strip())
00235
00236 if keyword == 'MASS':
00237 atype = line[0:2].strip()
00238 s = line[2:].split()
00239 mass = float(s[0])
00240 ff_dict[keyword].append([atype, mass])
00241 elif keyword == 'BOND':
00242 atype1 = line[:2].strip()
00243 atype2 = line[3:5].strip()
00244 s = line[5:].split()
00245 k = float(s[0])
00246 d0 = float(s[1])
00247 ff_dict[keyword].append([atype1, atype2, d0, k])
00248 elif keyword == 'ANGL':
00249 atype1 = line[:2].strip()
00250 atype2 = line[3:5].strip()
00251 atype3 = line[6:8].strip()
00252 s = line[8:].split()
00253 k = float(s[0])
00254 a0 = float(s[1])
00255 ff_dict[keyword].append([atype1, atype2, atype3, a0, k])
00256 elif keyword == 'DIHE':
00257 atype1 = line[:2].strip()
00258 atype2 = line[3:5].strip()
00259 atype3 = line[6:8].strip()
00260 atype4 = line[9:11].strip()
00261 s = line[11:].split()
00262 idiv = float(s[0])
00263 k = float(s[1])
00264 phase = float(s[2])
00265
00266
00267 period = abs(float(s[3]))
00268 ff_dict[keyword].append([atype1, atype2, atype3, atype4, idiv,
00269 period, phase, k/float(idiv)])
00270 elif keyword == 'IMPR':
00271 atype1 = line[:2].strip()
00272 atype2 = line[3:5].strip()
00273 atype3 = line[6:8].strip()
00274 atype4 = line[9:11].strip()
00275 s = line[11:].split()
00276 k = float(s[0])
00277 phase = float(s[1])
00278 period = float(s[2])
00279 ff_dict[keyword].append([atype1, atype2, atype3, atype4, period,
00280 phase, k])
00281 elif keyword == 'NONB':
00282 line = line.strip()
00283 atype = line[0:2].strip()
00284 s = line[2:].split()
00285 Rvdw = float(s[0])
00286 epsilon = float(s[1])
00287 ff_dict[keyword].append([atype, Rvdw, epsilon])
00288
00289 line = in_file.next()
00290 return ff_dict
00291
00292
00293 def RunAntechamber(res_name, filename, format='ccif', amberhome=None,
00294 base_out_dir=None):
00295 """Run Antechamber to guess force field parameters for a given residue name.
00296
00297 This requires an installation of AmberTools (tested with AmberTools15) with
00298 binaries ``antechamber`` and ``parmchk2``.
00299
00300 This has the same restrictions as Antechamber itself and we assume the input
00301 to be uncharged. Note that Antechamber cannot deal with metal ions and other
00302 non-organic elements.
00303
00304 The results are stored in a separate folder named `res_name` within
00305 `base_out_dir` (if given, otherwise the current working directory). The main
00306 output files are ``frcmod`` and ``out.mpdb``. The former contains force field
00307 parameters and masses. The latter maps atom names to atom types and defines
00308 the partial charges. The same output could be obtained as follows:
00309
00310 .. code-block:: console
00311
00312 $ antechamber -i <FILENAME> -fi <FORMAT> -bk '<RES_NAME>' -o out.mol2 -fo mol2 -c bcc -pf yes
00313 $ parmchk2 -i out.mol2 -f mol2 -o frcmod -a Y
00314 $ antechamber -i out.mol2 -fi mol2 -o out.mpdb -fo mpdb -pf yes
00315
00316 The force field parameters can be manually modified if needed. It can for
00317 instance happen that some parameters cannot be identified. Those lines will
00318 be marked with a comment "ATTN, need revision".
00319
00320 :param res_name: Residue name for which we desire force field parameters.
00321 :type res_name: :class:`str`
00322 :param filename: Path to a file which contains the necessary information for
00323 `res_name`. It must include all hydrogens.
00324 :type filename: :class:`str`
00325 :param format: Format of file given with `filename`. Common formats are 'ccif'
00326 for PDB's component dictionary or 'pdb' for a PDB file
00327 containing the desired residue with all hydrogens.
00328 :type format: :class:`str`
00329 :param amberhome: Base path of your AmberTools installation. If not None,
00330 we look for ``antechamber`` and ``parmchk2`` within
00331 ``AMBERHOME/bin`` additionally to the system's ``PATH``.
00332 :type amberhome: :class:`str`
00333 :param base_out_dir: Path to a base path, where the output will be stored.
00334 If None, the current working directory is used.
00335 :type base_out_dir: :class:`str`
00336 """
00337
00338 if amberhome is None:
00339 search_paths = []
00340 else:
00341 search_paths = [os.path.join(amberhome, 'bin')]
00342 try:
00343 antechamber = settings.Locate('antechamber', search_paths=search_paths)
00344 parmchk2 = settings.Locate('parmchk2', search_paths=search_paths)
00345 except settings.FileNotFound as ex:
00346 ost.LogError("Failed to find Antechamber binaries. Make sure you have "
00347 "AmberTools installed!")
00348 raise ex
00349
00350
00351 cwd = os.getcwd()
00352 if base_out_dir is None:
00353 base_out_dir = cwd
00354 out_dir = os.path.abspath(os.path.join(base_out_dir, res_name))
00355 if not os.path.exists(out_dir):
00356
00357 try:
00358 os.makedirs(out_dir)
00359 except Exception as ex:
00360 ost.LogError("Failed to create output directory " + out_dir + "!")
00361 raise ex
00362
00363
00364 os.chdir(out_dir)
00365 try:
00366 cmds = [antechamber + " -i " + filename + " -fi " + format + " -bk " \
00367 + res_name + " -o out.mol2 -fo mol2 -c bcc -pf yes",
00368 parmchk2 + " -i out.mol2 -f mol2 -o frcmod -a Y",
00369 antechamber + " -i out.mol2 -fi mol2 -o out.mpdb -fo mpdb -pf yes"]
00370 all_sout = "Generating force field parameters for " + res_name + "\n"
00371 all_serr = ""
00372 for cmd in cmds:
00373 all_sout += "-"*70 + "\n" + "Stdout of: " + cmd + "\n" + "-"*70 + "\n"
00374 all_serr += "-"*70 + "\n" + "Stderr of: " + cmd + "\n" + "-"*70 + "\n"
00375 job = subprocess.Popen(cmd.split(" "), stdout=subprocess.PIPE,
00376 stderr=subprocess.PIPE)
00377 sout, serr = job.communicate()
00378 all_sout += sout
00379 all_serr += serr
00380 if job.returncode != 0:
00381 ost.LogError("Unsuccessful execution of " + cmd + ". Return code: "\
00382 + str(job.returncode))
00383
00384 with open("console.stdout", "w") as txt_file:
00385 txt_file.write(all_sout)
00386 with open("console.stderr", "w") as txt_file:
00387 txt_file.write(all_serr)
00388 except Exception as ex:
00389 ost.LogError("Failed to excecute antechamber binaries!")
00390 raise ex
00391
00392
00393 os.chdir(cwd)
00394
00395
00396 frcmod_filename = os.path.join(out_dir, 'frcmod')
00397 mpdb_filename = os.path.join(out_dir, 'out.mpdb')
00398 if not os.path.exists(frcmod_filename):
00399 raise RuntimeError("Failed to generate frcmod file with Antechamber!")
00400 if not os.path.exists(mpdb_filename):
00401 raise RuntimeError("Failed to generate out.mpdb file with Antechamber!")
00402
00403 def AddFromFiles(force_field, frcmod_filename, mpdb_filename):
00404 """Add data from a frcmod and an mpdb file to a force field.
00405
00406 This will add a new :class:`~ost.mol.mm.BuildingBlock` to `force_field` for
00407 the residue defined in those files (residue name is extracted from the mpdb
00408 file which can only contain a single residue). Charges for each atom are
00409 extracted from the mpdb file. According to the frcmod file, an
00410 :class:`~ost.mol.mm.Interaction` is added for each bond, angle, dihedral and
00411 improper. Atom types with masses and non-bonded interactions are added to
00412 `force_field` as needed.
00413
00414 :param force_field: A force field object to which the new parameters are
00415 added.
00416 :type force_field: :class:`~ost.mol.mm.Forcefield`
00417 :param frcmod_filename: Path to ``frcmod`` file as generated by ``parmchk2``.
00418 :type frcmod_filename: :class:`str`
00419 :param mpdb_filename: Path to mpdb file as generated by ``antechamber``.
00420 :type mpdb_filename: :class:`str`
00421 :return: The updated force field (same as `force_field`).
00422 :rtype: :class:`~ost.mol.mm.Forcefield`
00423 """
00424
00425 if not os.path.exists(frcmod_filename):
00426 raise RuntimeError("Could not find frcmod file: " + frcmod_filename)
00427 if not os.path.exists(mpdb_filename):
00428 raise RuntimeError("Could not find mpdb file: " + mpdb_filename)
00429
00430 try:
00431 eh, res_name = _ParseModifiedPDB(mpdb_filename)
00432 except Exception as ex:
00433 ost.LogError("Failed to parse mpdb file: " + mpdb_filename)
00434 raise ex
00435 try:
00436 ff_dict = _ParseAmberForceField(frcmod_filename)
00437 except Exception as ex:
00438 ost.LogError("Failed to parse frcmod file: " + frcmod_filename)
00439 raise ex
00440 ost.LogInfo("Adding force field for " + res_name)
00441
00442 for aname, mass in ff_dict['MASS']:
00443 force_field.AddMass(aname, mass)
00444
00445 lj_sigma_scale = 2./10./2**(1./6.)
00446 lj_epsilon_scale = 4.184
00447 for aname, Rvdw, epsilon in ff_dict['NONB']:
00448
00449 if Rvdw == 0 or epsilon == 0:
00450 Rvdw, epsilon = 1, 0
00451 lj = mm.Interaction(mm.FuncType.LJ)
00452 lj.SetTypes([aname])
00453 lj.SetParam([Rvdw*lj_sigma_scale, epsilon*lj_epsilon_scale])
00454 force_field.AddLJ(lj)
00455
00456 bb = _MakeComponentBuildingBlock(eh, ff_dict)
00457 force_field.AddBuildingBlock(res_name, bb)
00458
00459 return force_field
00460
00461 def AddFromPath(force_field, out_dir):
00462 """Add data from a directory created with :meth:`Run` to a force field.
00463 See :meth:`AddFromFiles` for details.
00464
00465 :param force_field: A force field object to which the new parameters are
00466 added.
00467 :type force_field: :class:`~ost.mol.mm.Forcefield`
00468 :param out_dir: Output directory as created with :meth:`Run`. Must contain
00469 files ``frcmod`` and ``out.mpdb``.
00470 :type out_dir: :class:`str`
00471 :return: The updated force field (same as `force_field`).
00472 :rtype: :class:`~ost.mol.mm.Forcefield`
00473 """
00474 frcmod_filename = os.path.join(out_dir, 'frcmod')
00475 mpdb_filename = os.path.join(out_dir, 'out.mpdb')
00476 return AddFromFiles(force_field, frcmod_filename, mpdb_filename)
00477
00478 def AddIon(force_field, res_name, atom_name, atom_mass, atom_charge, lj_sigma,
00479 lj_epsilon):
00480 """Add a single atom as an ion to a force field.
00481
00482 Since Antechamber cannot deal with ions, you can add simple ones easily with
00483 this function. This adds a :class:`~ost.mol.mm.BuildingBlock` to `force_field`
00484 for the given residue name containing a single atom. The atom will have a type
00485 with the same name as the atom name and the given mass, charge and non-bonded
00486 (LJ) interaction parameters.
00487
00488 :param force_field: A force field object to which the ion is added.
00489 :type force_field: :class:`~ost.mol.mm.Forcefield`
00490 :param res_name: Residue name for the ion to be added.
00491 :type res_name: :class:`str`
00492 :param atom_name: Atom name which is also used as atom type name.
00493 :type atom_name: :class:`str`
00494 :param atom_mass: Mass of the atom.
00495 :type atom_mass: :class:`float`
00496 :param atom_charge: Charge of the atom.
00497 :type atom_charge: :class:`float`
00498 :param lj_sigma: The sigma parameter for the non-bonded LJ interaction.
00499 :type lj_sigma: :class:`float` in nm
00500 :param lj_epsilon: The sigma parameter for the non-bonded LJ interaction.
00501 :type lj_epsilon: :class:`float` in kJ/mol
00502 """
00503
00504 force_field.AddMass(atom_name, atom_mass)
00505
00506 bb = mm.BuildingBlock()
00507 bb.AddAtom(atom_name, atom_name, atom_charge, atom_mass)
00508 force_field.AddBuildingBlock(res_name, bb)
00509
00510 lj = mm.Interaction(mm.FuncType.LJ)
00511 lj.SetTypes([atom_name])
00512 lj.SetParam([lj_sigma, lj_epsilon])
00513 force_field.AddLJ(lj)
00514
00515 __all__ = ('RunAntechamber', 'AddFromFiles', 'AddFromPath', 'AddIon',)