diff --git a/intermol/desmond/desmond_parser.py b/intermol/desmond/desmond_parser.py index a4fc55c5..6418359e 100644 --- a/intermol/desmond/desmond_parser.py +++ b/intermol/desmond/desmond_parser.py @@ -1,3 +1,4 @@ +from collections import OrderedDict, defaultdict import logging from warnings import warn import math @@ -111,6 +112,24 @@ class DesmondParser(object): lookup_desmond_bonds = create_lookup(desmond_bonds) # not unique desmond_bond_types = create_type(desmond_bonds) + # desmond virtuals are not quite the same as gromacs virtuals. + # lc2, lc3, out3 are the same. + # fdat3 is generalization of gromacs 3fd and 3fda. Eventually this conversion + # should be moved to a canonicalization file. + # 4fdn is not defined in desmond + + desmond_virtuals = { + 'lc2': TwoVirtual, + 'lc3': ThreeLinearVirtual, + 'fdat3': ThreeFdVirtual, + 'fdat3': ThreeFadVirtual, + 'out3': ThreeOutVirtual, + None: FourFdnVirtual + } + + lookup_desmond_virtuals = create_lookup(desmond_virtuals) # not unique + desmond_virtual_types = create_type(desmond_virtuals) + def canonical_bond(self, bond, params, direction='into', name=None): if direction == 'into': @@ -1162,15 +1181,21 @@ def write_vdwtypes_and_sites(self, molecule): combrule = self.system.combination_rule for atom in molecule.atoms: i+=1 + if atom.ptype == 'A': + atom_type = 'atom' + elif atom.ptype == 'D': + atom_type = 'pseudo' + else: + raise DesmondError('The physical type of an atom is {0}, which is not allowed'.format(atom.ptype)) if atom.residue_index: - sites.append(' %3d %5s %9.8f %9.8f %2s %1d %4s\n' % ( - i, 'atom', + sites.append(' %3d %6s %9.8f %9.8f %2s %1d %4s\n' % ( + i, atom_type, atom._charge[0].value_in_unit(units.elementary_charge), atom._mass[0].value_in_unit(units.atomic_mass_unit), atom.atomtype[0], atom.residue_index, atom.residue_name)) else: sites.append(' %3d %5s %9.8f %9.8f %2s\n' % ( - i, 'atom', + i, atom_type, atom._charge[0].value_in_unit(units.elementary_charge), atom._mass[0].value_in_unit(units.atomic_mass_unit), atom.atomtype[0])) @@ -1649,6 +1674,114 @@ def write_constraints(self, moleculetype): hlines.extend(dlines) return hlines + def write_virtuals(self, moleculetype): + + #ADDING VIRTUALS + logger.debug(" -Writing virtuals...") + + virtuals = defaultdict(list) + virtuallist = sorted(list(moleculetype.virtual_forces),key=lambda x: x.atom1) + + dlines = list() + hlines = list() + + for virtual in virtuallist: + if hasattr(virtual, 'atom5'): + virtuals[4].append(virtual) + elif hasattr(virtual, 'atom4'): + virtuals[3].append(virtual) + else: + virtuals[2].append(virtual) + + virtuals[2] = sorted(virtuals[2], key=lambda x: (x.atom1, x.atom2, x.atom3)) + virtuals[3] = sorted(virtuals[3], key=lambda x: (x.atom1, x.atom2, x.atom3, x.atom4)) + virtuals[4] = sorted(virtuals[4], key=lambda x: (x.atom1, x.atom2, x.atom3, x.atom4, x.atom5)) + + alen_max = 2 # figure out the maximum number of atoms in all of the virtuals. + for j in range(2,5): + if virtuals[j] > 0: + alen_max = j + + i = 0 + + clen_max = 3 # maximum number of parameters involved. We won't try to automatically find. + for n_body_type, vsites in virtuals.items(): + for vsite in vsites: + i += 1 + cline = ' %d ' % i + for n in range(1, n_body_type + 2): + atom = getattr(vsite, 'atom{}'.format(n)) + cline += '{0:7d} '.format(atom) + for j in range(n_body_type + 1, alen_max): + cline += ' <>' + cline += ' {:4s} '.format(self.lookup_desmond_virtuals[vsite.__class__]) + + vsite_params = self.get_parameter_list_from_force(vsite) + param_units = self.unitvars[vsite.__class__.__name__] + for param, unit in zip(vsite_params, param_units): + cline += '{0:18.8e}'.format(param.value_in_unit(unit)) + for j in range(len(vsite_params),clen_max): + cline += ' <>' + cline += '\n' + dlines.append(cline) + + hlines.append(" ffio_virtuals[%d] {\n" % (i)) + if i == 0: + hlines.append(" :::\n") + else: + #decide on alen and clen based on the type of virtual + letters = ['i','j','k','l','m'] + hlines.append(' i_ffio_index\n') + for j in range(alen_max-1): + hlines.append(' i_ffio_a%s\n' % letters[j]) + hlines.append(' s_ffio_funct\n') + for j in range(clen_max): + hlines.append(' r_ffio_c%d\n' %(j+1)) + hlines.append(" :::\n") + + dlines.append(" :::\n") + dlines.append(" }\n") + hlines.extend(dlines) + return hlines + + def write_pseudos(self, moleculetype): + + dlines = list() + hlines = list() + + i = 0 + for molecule in moleculetype.molecules: + for atom in molecule.atoms: + if atom.ptype == 'D': + i += 1 + line = ' %d ' % i + for j in range(3): + line += " %10.8f" % (float(atom._position[j].value_in_unit(units.angstroms))) + line += " %2d %4s %2s" % ( + atom.residue_index, + '"%s"'%atom.residue_name, + '"%s"'%atom.name) + dlines.append(line + '\n') + + hlines.append(" ffio_pseudo[%d] {\n" % (i)) + if i == 0: + hlines.append(" :::\n") + else: + hlines.append(' r_ffio_x_coord\n') + hlines.append(' r_ffio_y_coord\n') + hlines.append(' r_ffio_z_coord\n') + hlines.append(' i_ffio_residue_number\n') + hlines.append(' s_ffio_pdb_residue_name\n') + hlines.append(' s_ffio_atom_name\n') + hlines.append(" :::\n") + + dlines.append(" :::\n") + dlines.append(" }\n") + + hlines.extend(dlines) + + return hlines + def write(self): # Write this topology to file @@ -1703,6 +1836,8 @@ def write(self): for moleculetype in self.system._molecule_types.values(): for molecule in moleculetype.molecules: for atom in molecule.atoms: + if atom.ptype == 'D': # dummy (virtual) atoms do not go here. + continue i += 1 line = ' %d %d' % (i,1) #HAVE TO PUT THE 1 HERE OR ELSE DESMOND DIES, EVEN THOUGH IT DOESN'T USE IT for j in range(3): @@ -1783,7 +1918,8 @@ def write(self): for molecule_name, moleculetype in self.system.molecule_types.items(): logger.debug('Writing molecule block %s...' % (molecule_name)) - #BEGINNING BLOCK + + #BEGINNING BLOCK - should be separate function. logger.debug(" Writing f_m_ct...") lines.append('f_m_ct {\n') @@ -1798,7 +1934,8 @@ def write(self): for bj in range(3): lines.append('%22s\n' % float(bv[bi][bj].value_in_unit(units.angstroms))) lines.append(' solute\n') - #M_ATOMS + + #M_ATOMS - should be separate function logger.debug(" Writing m_atoms...") apos = len(lines) #pos of where m_atom will be; will need to overwite later based on the number of atoms @@ -1812,6 +1949,8 @@ def write(self): i = 0 for molecule in moleculetype.molecules: for atom in molecule.atoms: + if atom.ptype == 'D': # dummy (virtual) atoms do not go here. + continue i += 1 #NOT SURE WHAT TO PUT FOR MMOD TYPE; 1 is currently used. #This can't be determined currently from the information provided, @@ -1836,7 +1975,8 @@ def write(self): lines.append(' :::\n') lines.append(' }\n') - #M_BONDS + + #M_BONDS - should be separate function logger.debug(" Writing m_bonds...") hlines = list() @@ -1918,6 +2058,8 @@ def write(self): lines += self.write_torsion_torsion(moleculetype) lines += self.write_exclusions(moleculetype) lines += self.write_pairs(moleculetype) + lines += self.write_virtuals(moleculetype) + lines += self.write_pseudos(moleculetype) lines += self.write_constraints(moleculetype) #STILL NEED TO ADD RESTRAINTS diff --git a/intermol/gromacs/gromacs_parser.py b/intermol/gromacs/gromacs_parser.py index 3ff749f7..a2ea46a3 100644 --- a/intermol/gromacs/gromacs_parser.py +++ b/intermol/gromacs/gromacs_parser.py @@ -398,6 +398,7 @@ def read(self): self.pairtypes = dict() self.cmaptypes = dict() self.nonbondedtypes = dict() + self.virtualtypes = dict() # Parse the top_filename into a set of plain text, intermediate # TopMoleculeType objects. @@ -797,6 +798,7 @@ def create_atom(self, temp_atom): "{0}. Visually inspect before using.".format(atom)) atom.sigma = (state, intermol_atomtype.sigma) atom.epsilon = (state, intermol_atomtype.epsilon) + atom.ptype = intermol_atomtype.ptype self.current_molecule.add_atom(atom) self.n_atoms_added += 1 @@ -1557,6 +1559,7 @@ def process_nonbond_params(self, line): self.system.nonbonded_types[tuple(nonbonded_vars)] = nonbonded_type def process_virtual_sites(self, line, v_site_type): + """Process a line in a [ virtual_sites? ] category.""" if v_site_type == 'n': raise UnimplementedSetting('Parsing of [ virtual_sitesn ] directives'