prepareComplexMolecules - set up proteins

various helper functions for setting up systems containing complex molecules such as proteins

espressopp.tools.findConstrainedBonds(atomPids, bondtypes, bondtypeparams, masses, massCutoff = 1.1)

Finds all heavyatom-hydrogen bonds in a given list of particle IDs, and outputs a list describing the bonds, in a format suitable for use with the RATTLE algorithm for constrained bonds

Parameters:
  • atomPids (list of int) – list of pids of atoms between which to search for bonds
  • bondtypes (dict, key: (int,int), value: int) – dictionary mapping from tuple of pids to bondtypeid, e.g. as returned by tools.gromacs.read()
  • bondtypeparams (dict, key: int, value: espressopp bond type) – dictionary mapping from bondtypeid to class storing parameters of that bond type, e.g. as returned by tools.gromacs.read()
  • masses (list of float) – list of masses, e.g. as returned by tools.gromacs.read()
  • massCutoff (float) – for identifying light atoms (hydrogens), default 1.1 mass units, can also be increased e.g. for use with deuterated systems
Returns:

hydrogenIDs - list of pids (integer) of light atoms (hydrogens)
constrainedBondsDict - dict, keys: pid (integer) of heavy atom, values: list of pids of light atoms that are bonded to it
constrainedBondsList - list of lists, one entry for each constrained bond with format: [pid of heavy atom, pid of light atom, bond distance, mass of heavy atom ,mass of light atom]

Can then be used with RATTLE, e.g.

>>> rattle = espressopp.integrator.Rattle(system, maxit = 1000, tol = 1e-6, rptol = 1e-6)
>>> rattle.addConstrainedBonds(constrainedBondsList)
>>> integrator.addExtension(rattle)
espressopp.tools.getInternalNonbondedInteractions(atExclusions, pidlist)

Gets the non-bonded pairs within a list of particle indices, excluding those which are in a supplied list of exclusions. Useful for example for getting the internal atomistic non-bonded interactions in a coarse-grained particle and adding them as a fixedpairlist

Parameters:
  • atExclusions (list of 2-tuples of int) – list of excluded pairs
  • pidlist (list of int) – list of pids among which to create pairs
Returns:

list of pairs which are not in atExclusions
Return type:

list of 2-tuples of int
espressopp.tools.readSimpleSystem(filename, nparticles, header)

Read in a column-formatted file containing information about the particles in a simple system, for example a coarsegrained protein.

This function expects the input file to have between 2 and 5 columns. The number of columns in the file is automatically detected. The function reads each column into a list and returns the lists. Column types are interpreted as follows: 2 columns: float, float 3 columns: float, float, int 4 columns: float, float, int, str 5 columns: float, float, int, str, str

For example in the case of a coarsegrained protein model, these could be: mass, charge, corresponding atomistic index, beadname, beadtype

Parameters:
  • filename (string) – name of file to open and read
  • nparticles (int) – number of particles in file
  • header (int) – number of lines to skip at start of file (default 0)

Returns: between 2 and 5 lists

espressopp.tools.applyBoreschRestraints(system, restraintAtoms, restraintK, restraintR0)

Applies restraints between ligand and protein as defined in Boresch et al, JPCB 2003, 107, 9535-9551 The restraints (one bond, two angles and three dihedrals) are applied between three ligand atoms A,B,C and three protein atoms a,b,c.

In espressopp, the potential for harmonic bond and angles is k(x-x0)^2, but the harmonic dihedral potential is 0.5*k(x-x0)^2. This is taken care of in this function, i.e. the user should supply the force constants exactly as given in Boresch et al.

Parameters:
  • system (System object) – espressopp system
  • restraintAtoms (python dictionary) – dictionary identifying the six atoms involved in the restraints, key = atom label (one of ‘A’,’B’,’C’,’a’,’b’,’c’), value = atom index
  • restraintK (python dictionary) – dictionary of force constants, key = restraint label (‘aA’ for the bond, ‘baA’ and ‘aAB’ for the angles, ‘aABC’, ‘cbaA’ and ‘baAB’ for the dihedrals), value = force constant (units as for the simulation forcefield, angle and dihedral constants should be in rad^-2)
  • restraintR0 (python dictionary) – dictionary of equilibrium values (distance or angle), key = restraint label, value = distance (in distance units) or angle (in degrees)

Examples of the three dictionaries:

>>> restraintAtoms = {'A':1981,'B':1966,'C':1993,'a':1588,'b':1581,'c':1567}
>>> restraintK = {'aA':4184,'baA':41.84,'aAB':41.84,'aABC':41.84,'cbaA':41.84,'baAB':41.84} #kJ mol-1 nm-2,kJ mol-1 rad-2, all 10 kcal as in JPCB 2003
>>> restraintR0 = {'aA':0.31,'baA':120.0,'aAB':90.0,'aABC':100.0,'cbaA':-170.0,'baAB':-105.0} #nm, degrees