diff --git a/intermol/desmond/desmond_parser.py b/intermol/desmond/desmond_parser.py index 90f1fa4a..8865cbd4 100644 --- a/intermol/desmond/desmond_parser.py +++ b/intermol/desmond/desmond_parser.py @@ -2,6 +2,7 @@ from warnings import warn import math import numpy as np +import shlex import parmed.unit as units from intermol.atom import Atom @@ -15,11 +16,11 @@ from intermol.system import System from intermol.desmond import cmap_parameters + #MRS for old desmond functionality import re import copy - logger = logging.getLogger('InterMolLog') ENGINE = 'desmond' @@ -51,7 +52,9 @@ def save(cms_file, system): # parser helper functions + def end_header_section(blank_section, header, header_lines): + if blank_section: header_lines = list() header_lines.append(header) @@ -62,29 +65,21 @@ def end_header_section(blank_section, header, header_lines): def split_with_quotes(line): - line = list(line) - in_quotes = False - for i, char in enumerate(line): - if char == '"': - in_quotes = not in_quotes - if char == ' ' and in_quotes: - line[i] = '_' - - space_split = "".join(line).split() - for i, sub in enumerate(space_split): - sub = sub.replace('"', '') - space_split[i] = sub.replace('_', ' ') - return space_split + if '"' in line: + elements = shlex.split(line) + for e in elements: + e.replace(' ','_') + else: + elements = line.split() + return elements def create_lookup(forward_dict): return dict((v, k) for k, v in forward_dict.items()) - def create_type(forward_dict): return dict((k, eval(v.__name__ + 'Type')) for k, v in forward_dict.items()) - class DesmondParser(object): """ A class containing methods required to read in a Desmond CMS File @@ -108,7 +103,7 @@ class DesmondParser(object): 'HARM': HarmonicBond } - lookup_desmond_bonds = create_lookup(desmond_bonds) # not unique + lookup_desmond_bonds = create_lookup(desmond_bonds) # not unique - revisit. desmond_bond_types = create_type(desmond_bonds) def canonical_bond(self, bond, params, direction='into', name=None): @@ -137,7 +132,7 @@ def canonical_bond(self, bond, params, direction='into', name=None): elif name == 'HARM': bond.c = False else: - warn("ReadError: Found unsupported bond in Desmond %s" % name) + warn("ReadError: Found unsupported bond in Desmond {:s}".format(name)) return bond else: @@ -392,9 +387,9 @@ def create_kwd_dict(self, forcetype_object, values, optvalues = None): def create_forcetype(self, forcetype_object, paramlist, values, optvalues = None): return forcetype_object(*paramlist, **self.create_kwd_dict(forcetype_object, values, optvalues)) -#LOAD FFIO BLOCKS IN FIRST (CONTAINS TOPOLOGY) - def parse_ffio_block(self,start,end): + #LOAD FFIO BLOCKS IN FIRST (CONTAINS TOPOLOGY) + # read in a ffio_block that isn't ffio_ff and split it into the # commands and the values. # lots of room for additional error checking here, such as whether @@ -403,12 +398,11 @@ def parse_ffio_block(self,start,end): # scroll to the next ffio entry while not 'ffio_' in self.lines[start]: # this is not an ffio block! or, we have reached the end of the file - if ('ffio_' not in self.lines[start]): - start+=1 + if 'ffio_' not in self.lines[start]: + start += 1 if start >= end: - return 'Done with ffio', 0, 0, 0, start + return 'Done with ffio', 0, None, None, None, start - self.lines[start].split()[1] components = re.split('\W', self.lines[start].split()[0]) # get rid of whitespace, split on nonword ff_type = components[0] ff_number = int(components[1]) @@ -427,24 +421,30 @@ def parse_ffio_block(self,start,end): i+=1 i+=1 # step past the end of the block - return ff_type,ff_number,entry_data,entry_values,i + entry_dict = dict() + for j, d in enumerate(entry_data): + entry_dict[d] = j+1 # the first one is the entry number - def store_ffio_data(self, ff_type, ff_number, entry_data, entry_values): + return ff_type, ff_number, entry_data, entry_values, entry_dict, i + + def store_ffio_data(self, ff_type, ff_number, entry_data, entry_values, entry_dict): self.stored_ffio_data[ff_type] = dict() self.stored_ffio_data[ff_type]['ff_type'] = ff_type self.stored_ffio_data[ff_type]['ff_number'] = ff_number self.stored_ffio_data[ff_type]['entry_data'] = entry_data self.stored_ffio_data[ff_type]['entry_values'] = entry_values + self.stored_ffio_data[ff_type]['entry_dict'] = entry_dict - def retrive_ffio_data(self, ff_type): + def retrieve_ffio_data(self, ff_type): return [self.stored_ffio_data[ff_type]['ff_number'], self.stored_ffio_data[ff_type]['entry_data'], - self.stored_ffio_data[ff_type]['entry_values'] + self.stored_ffio_data[ff_type]['entry_values'], + self.stored_ffio_data[ff_type]['entry_dict'] ] def parse_vdwtypes(self, type, current_molecule_type): - ff_number, entry_data, entry_values = self.retrive_ffio_data(type) + ff_number, entry_data, entry_values, entry_dict = self.retrieve_ffio_data(type) # molecule name is at sites, but vdwtypes come # before sites. So we store info in vdwtypes and # edit it later at sites. Eventually, we should @@ -458,7 +458,7 @@ def parse_vdwtypes(self, type, current_molecule_type): self.vdwtypeskeys.append(entry_values[j].split()[1]) def parse_sites(self, type, molname, i, start): - ff_number, entry_data, entry_values = self.retrive_ffio_data(type) + ff_number, entry_data, entry_values, entry_dict = self.retrieve_ffio_data(type) #correlate with atomtypes and atoms in GROMACS logger.debug("Parsing [ sites ] ...") @@ -485,21 +485,21 @@ def parse_sites(self, type, molname, i, start): current_molecule = Molecule(name=molname) # should this be the same molname several as lines up? for j in range(ff_number): - split = entry_values[j].split() - if split[1] == "atom": + values = split_with_quotes(entry_values[j]) + if values[1] == "atom": if ('i_ffio_resnr' in entry_data): - atom = Atom(int(split[0]), split[ivdwtype], - int(split[iresnum]), - split[iresidue]) + atom = Atom(int(values[0]), values[ivdwtype], + int(values[iresnum]), + values[iresidue]) else: # No residuenr, means we will have identical atoms sharing this. - atom = Atom(int(split[0]), split[ivdwtype]) + atom = Atom(int(values[0]), values[ivdwtype]) - atom.atomtype = (0, split[ivdwtype]) - atom.charge = (0, float(split[icharge])*units.elementary_charge) - atom.mass = (0, float(split[imass]) * units.amu) - stemp = float(self.vdwtypes[self.vdwtypeskeys.index(split[ivdwtype])][0]) * units.angstroms #was in angstroms - etemp = float(self.vdwtypes[self.vdwtypeskeys.index(split[ivdwtype])][1]) * units.kilocalorie_per_mole #was in kilocal per mol + atom.atomtype = (0, values[ivdwtype]) + atom.charge = (0, float(values[icharge])*units.elementary_charge) + atom.mass = (0, float(values[imass]) * units.amu) + stemp = float(self.vdwtypes[self.vdwtypeskeys.index(values[ivdwtype])][0]) * units.angstroms #was in angstroms + etemp = float(self.vdwtypes[self.vdwtypeskeys.index(values[ivdwtype])][1]) * units.kilocalorie_per_mole #was in kilocal per mol atom.sigma = (0, stemp) atom.epsilon = (0, etemp) atom.cgnr = cgnr @@ -511,20 +511,20 @@ def parse_sites(self, type, molname, i, start): if self.system.combination_rule == 'Multiply-C6C12': sigma = (etemp/stemp)**(1/6) epsilon = (stemp)/(4*sigma**6) - newAtomType = AtomCType(split[ivdwtypes], #atomtype/name - split[ivdwtype], #bondtype + newAtomType = AtomCType(values[ivdwtypes], #atomtype/name + values[ivdwtype], #bondtype -1, #atomic_number - float(split[imass]) * units.amu, #mass - float(split[icharge]) * units.elementary_charge, #charge--NEED TO CONVERT TO ACTUAL UNIT + float(values[imass]) * units.amu, #mass + float(values[icharge]) * units.elementary_charge, #charge--NEED TO CONVERT TO ACTUAL UNIT 'A', #pcharge...saw this in top--NEED TO CONVERT TO ACTUAL UNITS sigma * units.kilocalorie_per_mole * angstroms**(6), epsilon * units.kilocalorie_per_mole * unit.angstro,s**(12)) elif (self.system.combination_rule == 'Lorentz-Berthelot') or (self.system.combination_rule == 'Multiply-Sigeps'): - newAtomType = AtomSigepsType(split[ivdwtype], #atomtype/name - split[ivdwtype], #bondtype + newAtomType = AtomSigepsType(values[ivdwtype], #atomtype/name + values[ivdwtype], #bondtype -1, #atomic_number - float(split[imass]) * units.amu, #mass--NEED TO CONVERT TO ACTUAL UNITS - float(split[icharge]) * units.elementary_charge, #charge--NEED TO CONVERT TO ACTUAL UNIT + float(values[imass]) * units.amu, #mass--NEED TO CONVERT TO ACTUAL UNITS + float(values[icharge]) * units.elementary_charge, #charge--NEED TO CONVERT TO ACTUAL UNIT 'A', #pcharge...saw this in top--NEED TO CONVERT TO ACTUAL UNITS stemp, etemp) @@ -552,7 +552,8 @@ def parse_sites(self, type, molname, i, start): return self.system._molecule_types[molname] def parse_bonds(self, type, current_molecule_type, i, start): - ff_number, entry_data, entry_values = self.retrive_ffio_data(type) + + ff_number, entry_data, ev, ed = self.retrieve_ffio_data(type) if len(self.bond_blockpos) > 1: #LOADING M_BONDS if self.bond_blockpos[0] < start: @@ -566,13 +567,15 @@ def parse_bonds(self, type, current_molecule_type, i, start): logger.debug("Parsing [ bonds ]...") for j in range(ff_number): - entries = entry_values[j].split() - key = entries[3].upper() - atoms = [int(x) for x in entries[1:3]] + values = split_with_quotes(ev[j]) + key = values[ed['s_ffio_funct']].upper() + atomnames = ['i_ffio_ai','i_ffio_aj'] + atoms = [int(values[ed[a]]) for a in atomnames] bondingtypes = [self.atomlist[atom-1].name for atom in atoms] atoms.extend(bondingtypes) - params = [float(x) for x in entries[4:6]] + cnames = ['r_ffio_c1','r_ffio_c2'] + params = [float(values[ed[x]]) for x in cnames] new_bond = self.create_forcetype(self.desmond_bonds[key], atoms, params) kwds = self.get_parameter_kwds_from_force(new_bond) new_bond = self.canonical_bond(new_bond, kwds, direction = 'into', name = key) @@ -586,31 +589,34 @@ def parse_bonds(self, type, current_molecule_type, i, start): current_molecule_type.bond_forces.add(new_bond) def parse_pairs(self, type, current_molecule_type): - ff_number, entry_data, entry_values = self.retrive_ffio_data(type) + + ff_number, entry_data, ev, ed = self.retrieve_ffio_data(type) logger.debug("Parsing [ pairs ] ...") for j in range(ff_number): ljcorr = False coulcorr = False new_pair = None - split = entry_values[j].split() - atoms = [int(x) for x in split[1:3]] + values = split_with_quotes(ev[j]) + atomnames = ['i_ffio_ai','i_ffio_aj'] + atoms = [int(values[ed[a]]) for a in atomnames] bondingtypes = [self.atomlist[atom-1].name for atom in atoms] params = atoms + bondingtypes - key = split[3].upper() + key = values[ed['s_ffio_funct']].upper() if key == "LJ12_6_SIG_EPSILON": - new_pair = self.create_forcetype(LjSigepsPair, params, [float(x) for x in split[4:6]]) + new_pair = self.create_forcetype(LjSigepsPair, params, + [float(values[ed['r_ffio_c1']]),float(values[ed['r_ffio_c2']])]) elif key == "LJ" or key == "COULOMB": # I think we just need LjSigepsPair, not LjPair? new_pair = self.create_forcetype(LjDefaultPair, params, [0, 0]) if key == "LJ": - ljcorr = float(split[4]) + ljcorr = float(values[ed['r_ffio_c1']]) new_pair.scaleLJ = ljcorr elif key == "COULOMB": - coulcorr = float(split[4]) + coulcorr = float(values[ed['r_ffio_c1']]) new_pair.scaleQQ = coulcorr else: - warn("ReadError: didn't recognize type %s in line %s", split[3], entry_values[j]) + warn("ReadError: didn't recognize type {:s} in line {:s}".format(key, ev[j])) # now, we catch the matches and read them into a single potential pair_match = current_molecule_type.match_pairs(new_pair) @@ -659,46 +665,53 @@ def parse_pairs(self, type, current_molecule_type): # same energy. This could eventually be improved by checking versus the sites. def parse_angles(self, type, current_molecule_type): - ff_number, entry_data, entry_values = self.retrive_ffio_data(type) + ff_number, entry_data, ev, ed = self.retrieve_ffio_data(type) logger.debug("Parsing [ angles ] ...") + for j in range(ff_number): - split = entry_values[j].split() - key = split[4].upper() - atoms = [int(x) for x in split[1:4]] + values = split_with_quotes(ev[j]) + key = values[ed['s_ffio_funct']].upper() + atomnames = ['i_ffio_ai','i_ffio_aj','i_ffio_ak'] + atoms = [int(values[ed[a]]) for a in atomnames] bondingtypes = [self.atomlist[atom-1].name for atom in atoms] atoms.extend(bondingtypes) - kwds = [float(x) for x in split[5:7]] + kwds = [float(values[ed['r_ffio_c1']]), float(values[ed['r_ffio_c2']])] new_angle = self.create_forcetype(self.desmond_angles[key], atoms, kwds) kwds = self.get_parameter_kwds_from_force(new_angle) new_angle = self.canonical_angle(new_angle, kwds, direction = 'into', name = key, molecule_type = current_molecule_type) - if new_angle: current_molecule_type.angle_forces.add(new_angle) def parse_dihedrals(self, type, current_molecule_type): - ff_number, entry_data, entry_values = self.retrive_ffio_data(type) + + ff_number, entry_data, ev, ed = self.retrieve_ffio_data(type) logger.debug("Parsing [ dihedrals ] ...") for j in range(ff_number): - split = entry_values[j].split() + values = split_with_quotes(ev[j]) new_dihedral = None dihedral_type = None - atoms = [int(x) for x in split[1:5]] + key = values[ed['s_ffio_funct']].upper() + atomnames = ['i_ffio_ai','i_ffio_aj','i_ffio_ak','i_ffio_al'] + atoms = [int(values[ed[a]]) for a in atomnames] bondingtypes = [self.atomlist[atom-1].name for atom in atoms] - key = split[5].upper() atoms.extend(bondingtypes) # not sure how to put the following lines in canonical, since it expects keywords, # not strings of variable length. will have to fix later. + cnames = [] + for e in entry_data: + if 'r_ffio_c' in e: + cnames.append(e) # append all of the constants, in order if key == "IMPROPER_HARM": - kwds = [float(split[6]), 2*float(split[7])] + kwds = [float(values[ed[cnames[0]]]), 2*float(values[ed[cnames[1]]])] # harmonic, multiple x2 for desmond convention. elif key == "PROPER_TRIG" or key == "IMPROPER_TRIG": - kwds = [float(x) for x in split[6:14]] + kwds = [float(values[ed[x]]) for x in cnames] elif key == "OPLS_PROPER" or key == "OPLS_IMPROPER": - # next 3 lines definitely not the right way to do it. - opls_kwds = {key: value for key, value in zip("c1 c2 c3 c4".split(), [units.kilocalorie_per_mole * float(s) for s in split[7:11]])} + opls_vals = [float(values[ed[x]]) * units.kilocalorie_per_mole for x in cnames[1:5]] + opls_kwds = dict(zip([x[-2:] for x in cnames[1:5]],opls_vals)) opls_kwds = convert_dihedral_from_fourier_to_trig(opls_kwds) kwds = np.zeros(8) # will fill this in later. new_dihedral = self.create_forcetype(self.desmond_dihedrals[key], atoms, kwds) @@ -715,44 +728,52 @@ def parse_dihedrals(self, type, current_molecule_type): current_molecule_type.dihedral_forces.add(new_dihedral) def parse_torsion_torsion(self, type, current_molecule_type): - ff_number, entry_data, entry_values = self.retrive_ffio_data(type) + + ff_number, entry_data, ev, ed = self.retrieve_ffio_data(type) logger.debug("Parsing [ torsion-torsion ] ...") + for j in range(ff_number): - split = entry_values[j].split() + values = split_with_quotes(ev[j]) new_torsiontorsion = None - key = split[9].upper() + key = values[ed['s_ffio_funct']].upper() if key == "CMAP": # we shouldn't need to try/accept because there are no units. - new_torsiontorsion = TorsionTorsionCMAP(int(split[1]), - int(split[2]), - int(split[3]), - int(split[4]), - int(split[5]), - int(split[6]), - int(split[7]), - int(split[8]), - 'cmap', - int(split[10])) + new_torsiontorsion = TorsionTorsionCMAP(int(values[ed['i_ffio_ai']]), + int(values[ed['i_ffio_aj']]), + int(values[ed['i_ffio_ak']]), + int(values[ed['i_ffio_al']]), + int(values[ed['i_ffio_am']]), + int(values[ed['i_ffio_an']]), + int(values[ed['i_ffio_ao']]), + int(values[ed['i_ffio_ap']]), + 'cmap', + int(values[ed['i_ffio_c1']])) else: - warn("ReadError: found unsupported torsion-torsion type in: %s" % str(line[i])) + warn("ReadError: found unsupported torsion-torsion type in: {:s}".format(str(line[i]))) if new_torsiontorsion: current_molecule_type.torsiontorsion_forces.add(new_torsiontorsion) def parse_exclusions(self, type, current_molecule_type): - ff_number, entry_data, entry_values = self.retrive_ffio_data(type) + + ff_number, entry_data, ev, entry_dict = self.retrieve_ffio_data(type) logger.debug("Parsing [ exclusions ] ...") + + # currently, assumes no comments, could be dangerous? for j in range(ff_number): - temp = entry_values[j].split() + temp = split_with_quotes(ev[j]) temp.remove(temp[0]) current_molecule_type.exclusions.add(tuple([int(x) for x in temp])) def parse_restraints(self, type, current_molecule_type): - ff_number, entry_data, entry_values = self.retrive_ffio_data(type) + + ff_number, entry_data, ev, ed = self.retrieve_ffio_data(type) logger.debug("Warning: Parsing [ restraints] not yet implemented") def parse_constraints(self, type, current_molecule_type): - ff_number, entry_data, entry_values = self.retrive_ffio_data(type) + + ff_number, entry_data, ev, ed = self.retrieve_ffio_data(type) logger.debug("Parsing [ constraints ] ...") + ctype = 1 funct_pos = 0 atompos = [] #position of atoms in constraints; spread all over the place @@ -772,28 +793,28 @@ def parse_constraints(self, type, current_molecule_type): for j in range(ff_number): # water constraints actually get written to rigidwater (i.e. settles) constraints. - if 'HOH' in entry_values[j] or 'AH' in entry_values[j]: - split = entry_values[j].split() + if 'HOH' in ev[j] or 'AH' in ev[j]: + values = split_with_quotes(ev[j]) tempatom = [] templength = [] for a in atompos: - if not '<>' in split[a]: - tempatom.append(int(split[a])) + if not '<>' in values[a]: + tempatom.append(int(values[a])) else: tempatom.append(None) for l in lenpos: - if not '<>' in split[l]: - if 'AH' in entry_values[j]: - templength.append(float(split[l])*units.angstroms) # Check units? + if not '<>' in values[l]: + if 'AH' in ev[j]: + templength.append(float(values[l])*units.angstroms) # Check units? else: templength.append(None*units.angstroms) - constr_type = split[funct_pos] + constr_type = values[funct_pos] if 'HOH' in constr_type: - dOH = float(split[lenpos[1]]) - if dOH != float(split[lenpos[2]]): - logger.debug("Warning: second length in a rigid water specification (%s) is not the same as the first (%s)" % (split[lenpos[1]],split[lenpos[2]])) - angle = float(split[lenpos[0]])/(180/math.pi) + dOH = float(values[lenpos[1]]) + if dOH != float(values[lenpos[2]]): + logger.debug("Warning: second length in a rigid water specification {:s} is not the same as the first {:s}".format(values[lenpos[1]],values[lenpos[2]])) + angle = float(values[lenpos[0]])/(180/math.pi) dHH = 2*dOH*math.sin(angle/2) params = [atompos[0], atompos[1], atompos[2], dOH*units.angstroms, dHH*units.angstroms] new_rigidwater = RigidWater(*params) @@ -809,7 +830,7 @@ def parse_constraints(self, type, current_molecule_type): if new_constraint: current_molecule_type.constraints.add(new_constraint) else: - warn("ReadError: found unsupported constraint type %s" % (entry_values[j])) + warn("ReadError: found unsupported constraint type {:s}".format(ev[j])) def load_ffio_block(self, molname, start, end): @@ -826,7 +847,7 @@ def load_ffio_block(self, molname, start, end): # of the ordering later self.stored_ffio_data = {} # dictionary of stored ffio_entries - split = [] + values = [] constraints = [] temp = [] @@ -845,7 +866,7 @@ def load_ffio_block(self, molname, start, end): self.system.nonbonded_function = 1 self.system.genpairs = 'yes' - logger.debug('Parsing [ molecule %s ]'%(molname)) + logger.debug('Parsing [ molecule {:s} ]'.format(molname)) logger.debug('Parsing [ ffio ]') while i < end: @@ -873,13 +894,14 @@ def load_ffio_block(self, molname, start, end): # MISSING: need to identify vdw rule here -- currently assuming LJ12_6_sig_epsilon! # skip to the next ffio entry - while not ("ffio" in self.lines[i]): + while not ('ffio' in self.lines[i]): i+=1 bPreambleRead = True - ff_type, ff_number, entry_data, entry_values, i = self.parse_ffio_block(i, end) + ff_type, ff_number, entry_data, entry_values, entry_dict, i = self.parse_ffio_block(i, end) + self.stored_ffio_types.append(ff_type) - self.store_ffio_data(ff_type,ff_number, entry_data, entry_values) + self.store_ffio_data(ff_type,ff_number, entry_data, entry_values, entry_dict) # Reorder so 'vdwtypes' is first, then 'sites'. Could eventually get some simplification # by putting sites first, but too much rewriting for now. @@ -916,11 +938,19 @@ def loadMBonds(self, lines, start, end, npermol): #adds new bonds for each molec # end: ending of where m_bondsends for each molecule logger.debug("Parsing [ m_bonds ] ...") + bg = False newbond_force = None - split = [] + values = [] i = start bonds = set() + # right now, hard coded with header: + ## First column is bond index # + # i_m_from + # i_m_to + # i_m_order + # which might need to be read at some point. + while i < end: if ':::' in lines[i]: if bg: @@ -929,15 +959,15 @@ def loadMBonds(self, lines, start, end, npermol): #adds new bonds for each molec bg = True i+=1 if bg: - split = lines[i].split() - atomi = int(split[1]) - atomj = int(split[2]) + values = split_with_quotes(lines[i]) + atomi = int(values[1]) + atomj = int(values[2]) bondingtypei = self.atomlist[atomi-1].name bondingtypej = self.atomlist[atomj-1].name params = [atomi, atomj, bondingtypei, bondingtypej] if atomi > npermol: # we've collected the number of atoms per molecule. Exit. break - order = int(split[3]) + order = int(values[3]) kwd = [0, 0] optkwd = {'order': order, 'c': False} new_bond = self.create_forcetype(HarmonicBond, params, kwd, optkwd) @@ -957,8 +987,8 @@ def loadMAtoms(self, lines, start, end, currentMolecule, slength): #adds positio # slength: number of unique atoms in m_atoms, used to calculate repetitions logger.debug("Parsing [ m_atom ] ...") - i = start + i = start bg = False pdbaname = "" aname = "" @@ -975,18 +1005,17 @@ def loadMAtoms(self, lines, start, end, currentMolecule, slength): #adds positio start += 1 for c in self.atom_col_vars: if c in lines[i]: - logger.debug(" Parsing [ %s ] ..." % c) + logger.debug(" Parsing [ {:s} ] ...".format(c)) cols[c] = i - start break i+=1 - atom = None - - newMoleculeAtoms = [] - j = 0 logger.debug(" Parsing atoms...") + atom = None + newMoleculeAtoms = [] molecules = [] + j = 0 while j < mult: newMolecule = copy.deepcopy(currentMolecule) for atom in newMolecule.atoms: @@ -1151,7 +1180,8 @@ def read(self): def write_vdwtypes_and_sites(self, molecule): - #-ADDING VDWTYPES AND SITES + logger.debug(" -Writing vdwtypes...") + i = 0 sites = [] vdwtypes = [] @@ -1163,13 +1193,13 @@ def write_vdwtypes_and_sites(self, molecule): for atom in molecule.atoms: i+=1 if atom.residue_index: - sites.append(' %3d %5s %9.8f %9.8f %2s %1d %4s\n' % ( + sites.append(' {:3d} {:5s} {:9.8f} {:9.8f} {:2s} {:1d} {:4s}\n'.format( i, 'atom', 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' % ( + sites.append(' {:3d} {:5s} {:9.8f} {:9.8f} {:2s}\n'.format( i, 'atom', atom._charge[0].value_in_unit(units.elementary_charge), atom._mass[0].value_in_unit(units.atomic_mass_unit), @@ -1183,12 +1213,14 @@ def write_vdwtypes_and_sites(self, molecule): elif combrule in ['Lorentz-Berthelot','Multiply-Sigeps']: stemp = sig etemp = ep - if ' %2s %18s %8.8f %8.8f\n' % (atom.atomtype[0], "LJ12_6_sig_epsilon", float(stemp), float(etemp)) not in vdwtypes: - vdwtypes.append(' %2s %18s %8.8f %8.8f\n' % (atom.atomtype[0], "LJ12_6_sig_epsilon", float(stemp), float(etemp))) + vdwstring = ' {:2s} {:18s} {:8.8f} {:8.8f}\n'.format(atom.atomtype[0], + "LJ12_6_sig_epsilon", + float(stemp), float(etemp)) + if vdwstring not in vdwtypes: + vdwtypes.append(vdwstring) lines = [] - logger.debug(" -Writing vdwtypes...") - lines.append(" ffio_vdwtypes[%d] {\n"%(len(vdwtypes))) + lines.append(" ffio_vdwtypes[{:d}] {{\n".format(len(vdwtypes))) lines.append(" s_ffio_name\n") lines.append(" s_ffio_funct\n") lines.append(" r_ffio_c1\n") @@ -1197,12 +1229,13 @@ def write_vdwtypes_and_sites(self, molecule): i = 0 for v in vdwtypes: i+=1 - lines.append(' %d%2s'%(i,v)) + lines.append(' {:d}{:2s}'.format(i,v)) lines.append(" :::\n") lines.append(" }\n") logger.debug(" -Writing sites...") - lines.append(" ffio_sites[%d] {\n"%(len(sites))) + + lines.append(" ffio_sites[{:d}] {{\n".format(len(sites))) lines.append(" s_ffio_type\n") lines.append(" r_ffio_charge\n") lines.append(" r_ffio_mass\n") @@ -1212,7 +1245,7 @@ def write_vdwtypes_and_sites(self, molecule): lines.append(" s_ffio_residue\n") lines.append(" :::\n") for s in sites: - lines.append(' %s'%(s)) + lines.append(' {:s}'.format(s)) lines.append(" :::\n") lines.append(" }\n") @@ -1221,7 +1254,6 @@ def write_vdwtypes_and_sites(self, molecule): def write_bonds(self, moleculetype): - #-ADDING BONDS logger.debug(" -Writing bonds...") dlines = list() @@ -1244,14 +1276,14 @@ def write_bonds(self, moleculetype): for nbond, name in enumerate(names): i += 1 converted_bond = self.desmond_bonds[name](*atoms, **paramlists[nbond]) - line = ' %d %d %d %s' %(i, atoms[0], atoms[1], name) + line = ' {:d} {:d} {:d} {:s}'.format(i, atoms[0], atoms[1], name) bond_params = self.get_parameter_list_from_force(converted_bond) param_units = self.unitvars[converted_bond.__class__.__name__] for param, param_unit in zip(bond_params, param_units): - line += " %15.8f" % (param.value_in_unit(param_unit)) + line += " {:15.8f}".format(param.value_in_unit(param_unit)) line += '\n' dlines.append(line) - header = " ffio_bonds[%d] {\n" % (i) + header = " ffio_bonds[{:d}] {{\n".format(i) hlines = end_header_section(i==0,header,hlines) dlines.append(" :::\n") @@ -1260,7 +1292,7 @@ def write_bonds(self, moleculetype): return hlines def write_angles(self, moleculetype): - #-ADDING ANGLES + logger.debug(" -Writing angles...") dlines = list() @@ -1284,15 +1316,15 @@ def write_angles(self, moleculetype): for nangle, name in enumerate(names): i+=1 converted_angle = self.desmond_angles[name](*atoms, **paramlists[nangle]) - line = ' %d %d %d %d %s' % (i, atoms[0], atoms[1], atoms[2], name) + line = ' {:d} {:d} {:d} {:d} {:s}'.format(i, atoms[0], atoms[1], atoms[2], name) angle_params = self.get_parameter_list_from_force(converted_angle) param_units = self.unitvars[converted_angle.__class__.__name__] for param, param_unit in zip(angle_params, param_units): - line += " %15.8f" % (param.value_in_unit(param_unit)) + line += " {:15.8f}".format(param.value_in_unit(param_unit)) line += '\n' dlines.append(line) - header = " ffio_angles[%d] {\n" % (i) + header = " ffio_angles[{:d}] {{\n".format(i) hlines = end_header_section(i==0,header,hlines) dlines.append(" :::\n") @@ -1301,8 +1333,9 @@ def write_angles(self, moleculetype): return hlines def write_dihedrals(self, moleculetype): - #-ADDING DIHEDRALS + logger.debug(" -Writing dihedrals...") + dlines = list() hlines = list() hlines.append(" ffio_dihedrals_placeholder\n") @@ -1316,31 +1349,31 @@ def write_dihedrals(self, moleculetype): hmax = 8 # assume the maximum number of dihedral terms (8) to simplify things for now for ih in range(hmax): - hlines.append(" r_ffio_c%d\n" %(ih)) + hlines.append(" r_ffio_c{:d}\n".format(ih)) hlines.append(" :::\n") i = 0 - #sorting by first index + + # sorting dihedrals by first index dihedrallist = sorted(list(moleculetype.dihedral_forces), key=lambda x: (x.atom1, x.atom2, x.atom3, x.atom4)) - # first, identify the number of terms we will print for dihedral in dihedrallist: atoms = [dihedral.atom1,dihedral.atom2,dihedral.atom3,dihedral.atom4] kwds = self.get_parameter_kwds_from_force(dihedral) names, paramlists = self.canonical_dihedral(dihedral, kwds, direction = 'from') for ndihedrals, name in enumerate(names): i+=1 - line = ' %d %d %d %d %d %s' %(i, atoms[0], atoms[1], atoms[2], atoms[3], name) + line = ' {:d} {:d} {:d} {:d} {:d} {:s}'.format(i, atoms[0], atoms[1], atoms[2], atoms[3], name) converted_dihedral= self.desmond_dihedrals[name](*atoms,**paramlists[ndihedrals]) dihedral_params = self.get_parameter_list_from_force(converted_dihedral) param_units = self.unitvars[converted_dihedral.__class__.__name__] for param, param_unit in zip(dihedral_params, param_units): - line += " %15.8f" % (param.value_in_unit(param_unit)) + line += " {:15.8}".format(param.value_in_unit(param_unit)) for j in range(8-len(dihedral_params)): - line += " %6.3f" % (0.0) + line += " {:6.3f}".format(0.0) line += '\n' dlines.append(line) - header = " ffio_dihedrals[%d] {\n" % (i) + header = " ffio_dihedrals[{:d}] {{\n".format(i) hlines = end_header_section(i==0,header,hlines) dlines.append(" :::\n") @@ -1350,7 +1383,7 @@ def write_dihedrals(self, moleculetype): def write_torsion_torsion(self, moleculetype): - # adding TORSION-TORSION terms + # currently just cmap logger.debug(" -Writing torsion-torsions...") hlines = list() @@ -1372,14 +1405,14 @@ def write_torsion_torsion(self, moleculetype): for torsiontorsion in moleculetype.torsiontorsion_forces: i+=1 # only type of torsion/torsion is CMAP currently - dlines.append(' %d %d %d %d %d %d %d %d %d %s %d\n' % ( + dlines.append(' {:d} {:d} {:d} {:d} {:d} {:d} {:d} {:d} {:d} {:s} {:d}\n'.format( i, int(torsiontorsion.atom1), int(torsiontorsion.atom2), int(torsiontorsion.atom3), int(torsiontorsion.atom4), int(torsiontorsion.atom5), int(torsiontorsion.atom6), int(torsiontorsion.atom7), int(torsiontorsion.atom8), 'cmap', torsiontorsion.chart)) - header = " ffio_torsion_torsion[%d] {\n"%(i) + header = " ffio_torsion_torsion[{:d}] {{\n".format(i) hlines = end_header_section(i==0,header,hlines) dlines.append(" :::\n") @@ -1402,9 +1435,9 @@ def write_torsion_torsion(self, moleculetype): return hlines def write_exclusions(self, moleculetype): - #ADDING EXCLUSIONS - i = 0 + logger.debug(" -Writing exclusions...") + hlines = list() dlines = list() hlines.append(" ffio_exclusions_placeholder\n") @@ -1412,6 +1445,7 @@ def write_exclusions(self, moleculetype): hlines.append(" i_ffio_aj\n") hlines.append(" :::\n") + i = 0 if moleculetype.nrexcl == 0: # Should probably be determined entirely by the bonds, @@ -1422,7 +1456,7 @@ def write_exclusions(self, moleculetype): exclusionlist = sorted(list(moleculetype.exclusions), key=lambda x: (x[0], x[1])) for exclusion in moleculetype.exclusions: i+=1 - dlines.append(' %d %d %d\n'%(i, int(exclusion[0]), int(exclusion[1]))) + dlines.append(' {:d} {:d} {:d}\n'.format(i, int(exclusion[0]), int(exclusion[1]))) else: if moleculetype.nrexcl > 4: @@ -1472,10 +1506,10 @@ def write_exclusions(self, moleculetype): for a in atomexclude: if (a > atom): i+=1 - dlines.append(' %d %d %d\n' % (i, atom, a)) + dlines.append(' {:d} {:d} {:d}\n'.format(i, atom, a)) - header = " ffio_exclusions[%d] {\n"%(i) + header = " ffio_exclusions[{:d}] {{\n".format(i) hlines = end_header_section(i==0,header,hlines) dlines.append(" :::\n") @@ -1485,7 +1519,6 @@ def write_exclusions(self, moleculetype): def write_pairs(self, moleculetype): - #-ADDING PAIRS logger.debug(" -Writing pairs...") dlines = list() @@ -1497,10 +1530,10 @@ def write_pairs(self, moleculetype): hlines.append(" r_ffio_c1\n") hlines.append(" r_ffio_c2\n") hlines.append(" :::\n") - i = 0 + i = 0 for pair in sorted(list(moleculetype.pair_forces), key=lambda x: (x.atom1, x.atom2)): - atoms = ' %d %d ' % (pair.atom1, pair.atom2) + atoms = ' {:d} {:d} '.format(pair.atom1, pair.atom2) # first, the COUL part. if pair.__class__ in (LjDefaultPair, LjqDefaultPair, LjSigepsPair, LjCPair): # the first two appear to be duplicates: consider merging. @@ -1509,11 +1542,11 @@ def write_pairs(self, moleculetype): else: scaleQQ = self.system.coulomb_correction i += 1 - dlines += ' %d %s Coulomb %10.8f <>\n' % (i, atoms, scaleQQ) + dlines += ' {:d} {:s} Coulomb {:10.8f} <>\n'.format(i, atoms, scaleQQ) elif pair._class in (LjqSigepsPair, LjqCPair): - warn("Desmond does not support pairtype %s!",pair.__class__.__name__ ) # may not be true? + warn("Desmond does not support pairtype {:s}!".format(pair.__class__.__name__ )) # may not be true? else: - warn("Unknown pair type %s!",pair.__class__.__name__ ) + warn("Unknown pair type {:s}!".format(pair.__class__.__name__)) # now the LJ part. if pair.__class__ in (LjDefaultPair,LjqDefaultPair): @@ -1522,7 +1555,7 @@ def write_pairs(self, moleculetype): else: scaleLJ = self.system.lj_correction i += 1 - dlines += ' %d %s LJ %10.8f <>\n' % (i, atoms, scaleLJ) + dlines += ' {:d} {:s} LJ {:10.8f} <>\n'.format(i, atoms, scaleLJ) elif pair.__class__ in (LjSigepsPair, LjqSigepsPair, LjCPair, LjqCPair): # Check logic here -- not clear that we can correctly determine which type it is. # Basically, I think it's whether scaleLJ is defined or not. @@ -1533,13 +1566,13 @@ def write_pairs(self, moleculetype): epsilon = pair.epsilon sigma = pair.sigma i += 1 - dlines += ' %d %s LJ12_6_sig_epsilon %10.8f %10.8f\n' % (i, atoms, + dlines += ' {:d} {:s} LJ12_6_sig_epsilon {:10.8f} {:10.8f}\n'.format(i, atoms, sigma.value_in_unit(units.angstroms), epsilon.value_in_unit(units.kilocalorie_per_mole)) else: - warn("Unknown pair type %s!",pair.__class__.__name__ ) + warn("Unknown pair type {:s}!".format(pair.__class__.__name__)) - header = " ffio_pairs[%d] {\n"%(i) + header = " ffio_pairs[{:d}] {{\n".format(i) hlines = end_header_section(i==0,header,hlines) dlines.append(" :::\n") @@ -1549,7 +1582,6 @@ def write_pairs(self, moleculetype): def write_constraints(self, moleculetype): - #ADDING CONSTRAINTS logger.debug(" -Writing constraints...") isHOH = False @@ -1572,22 +1604,23 @@ def write_constraints(self, moleculetype): clen_max = clen # we now know the maximum length of all constraint types + # not sure we need to sort these, but makes it easier to debug - i = 0 constraintlist = sorted(list(moleculetype.constraints),key=lambda x: x.atom1) dlines = list() hlines = list() + i = 0 for constraint in constraintlist: #calculate the max number of atoms in constraint i+=1 if constraint.type == 'HOH': - cline = ' %d %d %d %d ' % (i,int(constraint.atom1),int(constraint.atom2),int(constraint.atom3)) + cline = ' {:d} {:d} {:d} {:d} '.format(i,int(constraint.atom1),int(constraint.atom2),int(constraint.atom3)) for j in range(alen_max-3): cline += '0 ' cline += constraint.type - cline += ' %10.8f' % (float(constraint.length1.value_in_unit(units.degrees))) - cline += ' %10.8f' % (float(constraint.length2.value_in_unit(units.angstroms))) - cline += ' %10.8f' % (float(constraint.length2.value_in_unit(units.angstroms))) + cline += ' {:10.8f}'.format(float(constraint.length1.value_in_unit(units.degrees))) + cline += ' {:10.8f}'.format(float(constraint.length2.value_in_unit(units.angstroms))) + cline += ' {:10.8f}'.format(float(constraint.length2.value_in_unit(units.angstroms))) for j in range(clen_max-3): cline += ' <>' elif constraint.type[0:2] == 'AH': @@ -1601,14 +1634,14 @@ def write_constraints(self, moleculetype): if hasattr(constraint,atomname): catoms.append(getattr(constraint,atomname)) clengths.append(getattr(constraint,lengthname)) - cline = ' %d ' % i + cline = ' {:d} '.format(i) for j in range(alen): - cline += ' %d ' % int(catoms[j]) + cline += ' {:d} '.format(int(catoms[j])) for j in range(alen,alen_max): cline += ' <> ' cline += constraint.type for j in range(clen): - cline += ' %10.8f' % (float(clengths[j].value_in_unit(units.angstroms))) + cline += ' {:10.8f}'.format(float(clengths[j].value_in_unit(units.angstroms))) for j in range(clen,clen_max): cline += ' <>' cline += '\n' @@ -1619,29 +1652,29 @@ def write_constraints(self, moleculetype): for rigidwater in moleculetype.rigidwaters: i += 1 # Assumes the water arrangement O, H, H, which might not always be the case. Consider adding detection. - cline = ' %d %d %d %d ' % (i, rigidwater.atom1, rigidwater.atom2, rigidwater.atom3) + cline = ' {:d} {:d} {:d} {:d} '.format(i, rigidwater.atom1, rigidwater.atom2, rigidwater.atom3) for j in range(alen_max-3): cline += '0 ' cline += ' HOH ' dOH = rigidwater.dOH.value_in_unit(units.angstroms) dHH = rigidwater.dHH.value_in_unit(units.angstroms) angle = 2.0*math.asin(0.5*dHH/dOH)*(180/math.pi) # could automate conversion. . . - cline += " %.8f %.8f %.8f" % (angle,dOH,dOH) + cline += ' {:.8f} {:.8f} {:.8f}'.format(angle,dOH,dOH) cline += '\n' for j in range(alen,alen_max): cline += ' 0.0' dlines.append(cline) - hlines.append(" ffio_constraints[%d] {\n"%(i)) + hlines.append(" ffio_constraints[{:d}] {{\n".format(i)) if (i==0): hlines.append(" :::\n") else: letters = ['i','j','k','l','m','n','o','p','q'] for j in range(alen_max): - hlines.append(' i_ffio_a%s\n'%letters[j]) + hlines.append(' i_ffio_a{:s}\n'.format(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(' r_ffio_c{:d}\n'.format(j+1)) hlines.append(" :::\n") dlines.append(" :::\n") @@ -1675,7 +1708,7 @@ def write(self): lines.append('f_m_ct {\n') lines.append(' s_m_title\n') for c in self.atom_box_vars: - lines.append(' %s\n' % c) + lines.append(' {:s}\n'.format(c)) lines.append(' s_ffio_ct_type\n') lines.append(' :::\n') @@ -1684,7 +1717,7 @@ def write(self): lines.append(' "Desmond file converted by InterMol"\n') for bi in range(3): for bj in range(3): - lines.append('%22s\n' % float(bv[bi][bj].value_in_unit(units.angstroms))) + lines.append('{:22.11f}\n'.format(float(bv[bi][bj].value_in_unit(units.angstroms)))) lines.append(' full_system\n') #M_ATOM @@ -1693,7 +1726,7 @@ def write(self): lines.append(' # First column is atom index #\n') for vars in self.atom_col_vars: if '_pdb_atom' not in vars: - lines.append(' %s\n' % vars) + lines.append(' {:s}\n'.format(vars)) lines.append(' :::\n') i = 0 @@ -1704,24 +1737,24 @@ def write(self): for molecule in moleculetype.molecules: for atom in molecule.atoms: i += 1 - line = ' %d %d' % (i,1) #HAVE TO PUT THE 1 HERE OR ELSE DESMOND DIES, EVEN THOUGH IT DOESN'T USE IT + line = ' {:d} {:d}'.format(i,1) #HAVE TO PUT THE 1 HERE OR ELSE DESMOND DIES, EVEN THOUGH IT DOESN'T USE IT for j in range(3): - line += " %10.8f" % (float(atom._position[j].value_in_unit(units.angstroms))) - line += " %2d %4s %2d %2s" % ( + line += " {:10.8f}".format(float(atom._position[j].value_in_unit(units.angstroms))) + line += " {:2d} {:4s} {:2d} {:2s}".format( atom.residue_index, - '"%s"'%atom.residue_name, + '"{:s}"'.format(atom.residue_name), atom.atomic_number, - '"%s"'%atom.name) + '"{:s}"'.format(atom.name)) if np.any(atom._velocity): for j in range(3): - line += " %10.8f" % (float(atom._velocity[j].value_in_unit(units.angstroms / units.picoseconds))) + line += " {:10.8f}".format(float(atom._velocity[j].value_in_unit(units.angstroms / units.picoseconds))) else: for j in range(3): - line += " %10.8f" % (0) + line += " {:10.8f}".format(0) lines.append(line + '\n') totalatoms.append(i) - lines[apos] = ' m_atom[%d] {\n'%(i) + lines[apos] = " m_atom[{:d}] {{\n".format(i) lines.append(' :::\n') lines.append(' }\n') @@ -1758,8 +1791,7 @@ def write(self): for bond in bondlist: if bond and bond.order: i += 1 - dlines.append(' %d %d %d %d %d %d\n' - %(i, + dlines.append(' {:d} {:d} {:d} {:d} {:d} {:d}\n'.format(i, bond.atom1 + n*atoms_per_molecule + totalatoms[nmol], bond.atom2 + n*atoms_per_molecule + totalatoms[nmol], int(bond.order), @@ -1768,10 +1800,10 @@ def write(self): elif not bond: nonecnt+=1 if nonecnt > 0: - logger.debug('FOUND %d BONDS THAT DO NOT EXIST' % nonecnt) + logger.debug('FOUND {:d} BONDS THAT DO NOT EXIST'.format(nonecnt)) nmol +=1 - hlines[0] = ' m_bond[%d] {\n' % i + hlines[0] = ' m_bond[{:d}] {{\n'.format(i) if (i > 0): lines.extend(hlines) lines.extend(dlines) @@ -1782,21 +1814,21 @@ def write(self): #WRITE OUT ALL FFIO AND F_M_CT BLOCKS for molecule_name, moleculetype in self.system.molecule_types.items(): - logger.debug('Writing molecule block %s...' % (molecule_name)) + logger.debug('Writing molecule block {:s}...'.format(molecule_name)) #BEGINNING BLOCK logger.debug(" Writing f_m_ct...") lines.append('f_m_ct {\n') lines.append(' s_m_title\n') for c in self.atom_box_vars: - lines.append(' %s\n' % c) + lines.append(' {:s}\n'.format(c)) lines.append(' s_ffio_ct_type\n') lines.append(' :::\n') lines.append(' "' + molecule_name + '"\n') for bi in range(3): for bj in range(3): - lines.append('%22s\n' % float(bv[bi][bj].value_in_unit(units.angstroms))) + lines.append('{:22.11f}\n'.format(float(bv[bi][bj].value_in_unit(units.angstroms)))) lines.append(' solute\n') #M_ATOMS @@ -1806,7 +1838,7 @@ def write(self): lines.append(' # First column is atom index #\n') for vars in self.atom_col_vars: if '_pdb_atom' not in vars: # kludge, have better filter - lines.append(' %s\n' % vars) + lines.append(' {:s}\n'.format(vars)) lines.append(' :::\n') i = 0 @@ -1816,27 +1848,26 @@ def write(self): #NOT SURE WHAT TO PUT FOR MMOD TYPE; 1 is currently used. #This can't be determined currently from the information provided, # unless it is stored previous, nor is it used by desmond - line = ' %d %d' % (i,1) + line = ' {:d} {:d}'.format(i,1) for j in range(3): - line += " %10.8f" % (float(atom._position[j].value_in_unit(units.angstroms))) - line += " %2d %4s %2d %2s" % ( + line += " {:10.8f}".format(float(atom._position[j].value_in_unit(units.angstroms))) + line += " {:2d} {:4s} {:2d} {:2s}".format( atom.residue_index, - '"%s"'%atom.residue_name, + '"{:s}"'.format(atom.residue_name), atom.atomic_number, - '"%s"'%atom.name) + '"{:s}"'.format(atom.name)) if np.any(atom._velocity): for j in range(3): - line += " %10.8f" % (float(atom._velocity[j].value_in_unit(units.angstroms / units.picoseconds))) + line += " {:10.8f}".format(float(atom._velocity[j].value_in_unit(units.angstroms / units.picoseconds))) else: for j in range(3): - line += " %10.8f" % (0) + line += " {:10.8f}".format(0) lines.append(line + '\n') - lines[apos] = ' m_atom[%d] {\n'%(i) + lines[apos] = ' m_atom[{:d}] {{\n'.format(i) lines.append(' :::\n') lines.append(' }\n') - #M_BONDS logger.debug(" Writing m_bonds...") hlines = list() @@ -1863,8 +1894,7 @@ def write(self): for bond in bondlist: if bond and bond.order: i += 1 - dlines.append(' %d %d %d %d %d %d\n' - %(i, + dlines.append(' {:d} {:d} {:d} {:d} {:d} {:d}\n'.format(i, bond.atom1 + n*atoms_per_molecule, bond.atom2 + n*atoms_per_molecule, int(bond.order), @@ -1873,9 +1903,9 @@ def write(self): else: nonecnt+=1 if nonecnt > 0: - logger.debug('FOUND %d BONDS THAT DO NOT EXIST' % nonecnt) + logger.debug('FOUND {:d} BONDS THAT DO NOT EXIST'.format(nonecnt)) - header = ' m_bond[%d] {\n'%i + header = ' m_bond[{:d}] {{\n'.format(i) if (i>0): hlines = end_header_section(False,header,hlines) @@ -1884,7 +1914,6 @@ def write(self): lines.append(' :::\n') lines.append(' }\n') - #FFIO # only need the first molecule molecule = next(iter(moleculetype.molecules)) logger.debug(" Writing ffio...") @@ -1898,7 +1927,7 @@ def write(self): if "Viparr" in molecule_name: lines.append(' Generated by Viparr\n') else: - lines.append(' "%s"\n' % molecule_name) + lines.append(' "{:s}"\n'.format(molecule_name)) #Adding Combination Rule if self.system.combination_rule == 'Multiply-C6C12': diff --git a/intermol/gromacs/__init__.py b/intermol/gromacs/__init__.py index b4d308ea..cd3735a5 100644 --- a/intermol/gromacs/__init__.py +++ b/intermol/gromacs/__init__.py @@ -29,6 +29,7 @@ 'CMAP': ['dihedral','cmap'], 'LJ (SR)': ['vdw total', 'vdw (SR)'], 'LJ-14': ['vdw total','vdw-14'], + 'LJ recip.': ['vdw total', 'vdw (LR)'], 'Disper. corr.': ['vdw total', 'vdw (LR)'], 'Coulomb (SR)': ['coulomb total','coulomb (SR)'], 'Coulomb-14': ['coulomb total','coulomb-14'], diff --git a/intermol/lammps/lammps_parser.py b/intermol/lammps/lammps_parser.py index f8309a53..b1b9a199 100644 --- a/intermol/lammps/lammps_parser.py +++ b/intermol/lammps/lammps_parser.py @@ -412,6 +412,8 @@ def read_data(self, data_file): self.parse_box(line.split(), 1) elif ('zlo' in line) and ('zhi' in line): self.parse_box(line.split(), 2) + elif ('xy xz yz' in line): + self.parse_box(line.split(), 'tilt') # Other headers. else: keyword = line.strip() @@ -608,12 +610,23 @@ def parse_box(self, line, dim): line (str): Current line in input file. dim (int): Dimension specified in line. """ - fields = [float(field) for field in line[:2]] - box_length = fields[1] - fields[0] - if box_length > 0: - self.system.box_vector[dim, dim] = box_length * self.DIST + if dim == 'tilt': + fields = [(float(field) * self.DIST) for field in line[:3]] + dirs = line[3:] + for f,d in zip(fields,dirs): + if d == 'xy': + self.system.box_vector[1,0] = f + elif d == 'xz': + self.system.box_vector[2,0] = f + elif d == 'yz': + self.system.box_vector[2,1] = f else: - raise LammpsError("Negative box length specified in data file.") + fields = [float(field) for field in line[:2]] + box_length = fields[1] - fields[0] + if box_length > 0: + self.system.box_vector[dim, dim] = box_length * self.DIST + else: + raise LammpsError("Negative box length specified in data file.") def parse_masses(self, data_lines): """Read masses from data file.""" @@ -1154,6 +1167,10 @@ def write(self, unit_set='real',nonbonded_style=None): f.write('{0:11.7f} {1:11.7f} zlo zhi\n'.format( z_min, z_min + self.system.box_vector[2][2].value_in_unit( self.DIST))) + f.write('{0:11.7f} {1:11.7f} {2:11.7f} xy xz yz\n'.format( + self.system.box_vector[1][0].value_in_unit(self.DIST), + self.system.box_vector[2][0].value_in_unit(self.DIST), + self.system.box_vector[2][1].value_in_unit(self.DIST))) for mass in mass_list: f.write(mass) diff --git a/intermol/tests/lammps/unit_tests/box_vectors_vacuum/box_vectors_vacuum.input b/intermol/tests/lammps/unit_tests/box_vectors_vacuum/box_vectors_vacuum.input new file mode 100644 index 00000000..cb3352d4 --- /dev/null +++ b/intermol/tests/lammps/unit_tests/box_vectors_vacuum/box_vectors_vacuum.input @@ -0,0 +1,18 @@ +units real +atom_style full + +dimension 3 +boundary p p p + +pair_style lj/cut/coul/cut 30.0 30.0 +pair_modify mix geometric + +bond_style hybrid harmonic morse +angle_style harmonic +special_bonds lj 0.0 0.0 0.5 coul 0.0 0.0 0.5 + +read_data box_vectors_vacuum.lmp + +thermo_style custom ebond eangle edihed eimp epair evdwl ecoul elong etail pe + +run 0 diff --git a/intermol/tests/lammps/unit_tests/box_vectors_vacuum/box_vectors_vacuum.lmp b/intermol/tests/lammps/unit_tests/box_vectors_vacuum/box_vectors_vacuum.lmp new file mode 100644 index 00000000..e2e8fd70 --- /dev/null +++ b/intermol/tests/lammps/unit_tests/box_vectors_vacuum/box_vectors_vacuum.lmp @@ -0,0 +1,97 @@ +topol + +9 atoms +8 bonds +13 angles +0 dihedrals +0 impropers + +4 atom types +4 bond types +6 angle types + + 26.160000 46.266100 xlo xhi + 28.470000 48.576100 ylo yhi + 26.030000 46.136100 zlo zhi + 0.0 -4.50000 0.0 xy xz yz + +Masses + +1 12.0110 +2 1.0080 +3 15.9994 +4 1.0080 + +Pair Coeffs + +1 0.0660 3.5000 +2 0.0300 2.5000 +3 0.1700 3.1200 +4 0.0000 0.0000 + +Bond Coeffs + +1 harmonic 340.00000000 1.09000000 +2 harmonic 553.00000000 0.94500000 +3 morse 320.00000000 1.41000000 1.2 +4 morse 268.00000000 1.52900000 1.3 + +Angle Coeffs + +1 50.00000000 109.50000000 +2 37.50000000 110.70000000 +3 37.50000000 110.70000000 +4 33.00000000 107.80000000 +5 35.00000000 109.50000000 +6 55.00000000 108.50000000 + +Atoms + + 1 1 1 2.0 27.110 29.460 28.030 + 2 1 2 -2.0 27.090 28.470 27.570 + 3 1 2 2.0 26.160 29.540 28.550 + 4 1 3 -2.0 28.040 29.570 29.070 + 5 1 4 2.0 28.830 29.460 28.570 + 6 1 1 -2.0 27.320 30.480 26.870 + 7 1 2 2.0 28.330 30.640 26.490 + 8 1 2 -2.0 27.000 31.450 27.230 + 9 1 2 2.0 26.630 30.350 26.030 + +Velocities + + 1 0.0000 0.0000 0.0000 + 2 0.0000 0.0000 0.0000 + 3 0.0000 0.0000 0.0000 + 4 0.0000 0.0000 0.0000 + 5 0.0000 0.0000 0.0000 + 6 0.0000 0.0000 0.0000 + 7 0.0000 0.0000 0.0000 + 8 0.0000 0.0000 0.0000 + 9 0.0000 0.0000 0.0000 + +Bonds + + 1 1 1 2 + 2 1 6 9 + 3 1 1 3 + 4 1 6 7 + 5 1 6 8 + 6 2 4 5 + 7 3 1 4 + 8 4 1 6 + +Angles + + 1 1 4 1 6 + 2 2 3 1 6 + 3 3 1 6 9 + 4 4 8 6 9 + 5 4 7 6 9 + 6 5 2 1 4 + 7 5 3 1 4 + 8 3 1 6 7 + 9 4 7 6 8 + 10 2 2 1 6 + 11 6 1 4 5 + 12 4 2 1 3 + 13 3 1 6 8 diff --git a/intermol/tests/lammps/unit_tests/box_vectors_vacuum/log.lammps b/intermol/tests/lammps/unit_tests/box_vectors_vacuum/log.lammps new file mode 100644 index 00000000..6d4f3181 --- /dev/null +++ b/intermol/tests/lammps/unit_tests/box_vectors_vacuum/log.lammps @@ -0,0 +1,62 @@ +LAMMPS (1 Feb 2014) +units real +atom_style full + +dimension 3 +boundary p p p + +pair_style lj/cut/coul/cut 30.0 30.0 +pair_modify mix geometric + +bond_style hybrid harmonic morse +angle_style harmonic +special_bonds lj 0.0 0.0 0.5 coul 0.0 0.0 0.5 + +read_data box_vectors.lmp + triclinic box = (26.16 28.47 26.03) to (46.2661 48.5761 46.1361) with tilt (0 -4.5 0) + 1 by 1 by 1 MPI processor grid + reading atoms ... + 9 atoms + reading velocities ... + 9 velocities + scanning bonds ... + 4 = max bonds/atom + scanning angles ... + 6 = max angles/atom + reading bonds ... + 8 bonds + reading angles ... + 13 angles + 4 = max # of 1-2 neighbors + 4 = max # of 1-3 neighbors + 7 = max # of 1-4 neighbors + 8 = max # of special neighbors + +thermo_style custom ebond eangle edihed eimp epair evdwl ecoul elong etail pe + +run 0 +WARNING: No fixes defined, atoms won't move (../verlet.cpp:54) +Memory usage per processor = 4.88161 Mbytes +E_bond E_angle E_dihed E_impro E_pair E_vdwl E_coul E_long E_tail PotEng + 47.912249 4.808182 0 0 354.9457 -0.086795261 355.0325 0 0 407.66613 +Loop time of 0 on 1 procs for 0 steps with 9 atoms + +Pair time (%) = 0 (0) +Bond time (%) = 0 (0) +Neigh time (%) = 0 (0) +Comm time (%) = 0 (0) +Outpt time (%) = 0 (0) +Other time (%) = 0 (0) + +Nlocal: 9 ave 9 max 9 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Nghost: 567 ave 567 max 567 min +Histogram: 1 0 0 0 0 0 0 0 0 0 +Neighs: 767 ave 767 max 767 min +Histogram: 1 0 0 0 0 0 0 0 0 0 + +Total # of neighbors = 767 +Ave neighs/atom = 85.2222 +Ave special neighs/atom = 7.33333 +Neighbor list builds = 0 +Dangerous builds = 0