Source code for parmed.tools.coarsegrain

"""
This module handles adding coarse graining data to topology files set up by
Luciana Capece. It has been implemented here in parmed by Jason Swails, but all
ideas were Luciana's. Here is a description of what this section does to
topology files:

ANGLE_FORCE_CONSTANT has been rigged to provide an integer index to the specific
parameter found in the given parameter file. This will be used to set up 4
different (new) prmtop sections: ANGLE_COEF_A, ANGLE_COEF_B, ANGLE_COEF_C, and
ANGLE_COEF_D. ANGLE_FORCE_CONSTANT now becomes useless, so I will add a comment
to the topology file saying such.

DIHEDRAL_FORCE_CONSTANT has been rigged in much the same way, creating 8 new
sections: DIHEDRAL_AMPLITUDE_1,2,3,4 and DIHEDRAL_PHASE_1,2,3,4.
DIHEDRAL_FORCE_CONSTANT, DIHEDRAL_PHASE, and DIHEDRAL_PERIODICITY become useless
and are commented as such (though not removed).

This module sets up ANGLE and DIHEDRAL classes to make printing them easier,
then parses the parameter file, then sets the respective sections of the prmtop
file.
"""
from parmed.utils.six.moves import range
import warnings

#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+

[docs]class Angle(object): """ Angle between 3 bonded atoms """ def __init__(self, atom1, atom2, atom3, acoef, bcoef, ccoef, dcoef): self.atom1 = atom1.strip() self.atom2 = atom2.strip() self.atom3 = atom3.strip() self.acoef = float(acoef) self.bcoef = float(bcoef) self.ccoef = float(ccoef) self.dcoef = float(dcoef) def __eq__(self, other): """ Defines if 2 angles are equal: the middle atom has to be the same, and the 2 outer ones have to be the same in either order, or they can be wildcards (X) """ # First see if we are looking to match the angles in forward or reverse reversed_angles = False if self.atom2 != other.atom2: return False if self.atom1 != other.atom1: if self.atom1 != other.atom3: if self.atom1 != 'X' and other.atom1 != 'X': if other.atom3 == 'X': reversed_angles = True else: return False # end if self.atom1 != 'X' and other.atom1 != 'X' else: reversed_angles = True # Now check the rest (non-reversed case first) if not reversed_angles: if (self.atom3 != other.atom3 and self.atom3 != 'X' and other.atom3 != 'X'): return False else: # reversed angles if (self.atom3 != other.atom1 and self.atom3 != 'X' and other.atom1 != 'X'): return False # If we passed all these tests so far, then we must be the same! return True
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]class Dihedral(object): """ Dihedral between 4 bonded atoms """ def __init__(self, atom1, atom2, atom3, atom4, ampl1, ampl2, ampl3, ampl4, phase1, phase2, phase3, phase4): self.atom1, self.atom2 = atom1, atom2 self.atom3, self.atom4 = atom3, atom4 self.ampl1, self.ampl2 = ampl1, ampl2 self.ampl3, self.ampl4 = ampl3, ampl4 self.phase1, self.phase2 = phase1, phase2 self.phase3, self.phase4 = phase3, phase4 def __eq__(self, other): """ 2 Dihedrals are equal if they go in forward/reverse order and all atom types either match or we have wild cards (which can only be on the ends) """ reversed_angles = False if self.atom2 != other.atom2 or self.atom3 != other.atom3: if self.atom2 != other.atom3 or self.atom3 != other.atom2: return False else: reversed_angles = True if not reversed_angles: if self.atom1 != other.atom1: if self.atom1 != 'X' and other.atom1 != 'X': return False if self.atom4 != other.atom4: if self.atom4 != 'X' and other.atom4 != 'X': return False else: if self.atom1 != other.atom4: if self.atom1 != 'X' and other.atom4 != 'X': return False if self.atom4 != other.atom1: if self.atom4 != 'X' and other.atom4 != 'X': return False # There must be no problems if we made it this far return True
#+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
[docs]def addCoarseGrain(parm, param_file): from parmed.tools.exceptions import CoarseGrainError """ Adds coarse graining sections to the prmtop """ # Angles and Dihedrals dictionaries for storing parameters parsed from # parameter file angle_params = {} dihedral_params = {} # First parse the params file, filling dictionaries of angles and dihedrals # whose keys are their integer index values params = open(param_file, 'r') reading_angles = False reading_dihedrals = False file_line = params.readline() while file_line: # Filter out comments if '#' in file_line: file_line = file_line[:file_line.index('#')] # Filter out blank lines file_line = file_line.strip() if len(file_line) == 0: file_line = params.readline() continue # See if we start reading ANGLEs yet: if file_line[:5].upper() == 'ANGLE': reading_angles = True reading_dihedrals = False file_line = params.readline() continue # See if we start reading DIHEDRALs yet if file_line[:8].upper() == 'DIHEDRAL': reading_dihedrals = True reading_angles = False file_line = params.readline() continue # Pass over any INDEX header lines -- they're just there for readability if (reading_dihedrals or reading_angles) and file_line[:5] == 'INDEX': file_line = params.readline() continue if reading_angles: line_parts = file_line.split() try: idx = int(line_parts[0]) acoef = float(line_parts[1]) bcoef = float(line_parts[2]) ccoef = float(line_parts[3]) dcoef = float(line_parts[4]) atom_types = line_parts[5].split('-') angle_params[idx] = Angle(atom_types[0], atom_types[1], atom_types[2], acoef, bcoef, ccoef, dcoef) except ValueError as err: raise CoarseGrainError( 'Unexpected format in Coarse Grain angles. Expected ' 'different data type: %s. See format specification' % err) except IndexError: raise CoarseGrainError( 'Unexpected format in Coarse Grain parameter file. ' 'Expected more data fields on the line. See format ' 'specification') file_line = params.readline() continue if reading_dihedrals: line_parts = file_line.split() try: idx = int(line_parts[0]) ampl1 = float(line_parts[1]) ampl2 = float(line_parts[2]) ampl3 = float(line_parts[3]) ampl4 = float(line_parts[4]) atom_types = line_parts[5].split('-') file_line = params.readline().strip() line_parts = file_line.split() phase1 = float(line_parts[0]) phase2 = float(line_parts[1]) phase3 = float(line_parts[2]) phase4 = float(line_parts[3]) dihedral_params[idx] = Dihedral(atom_types[0], atom_types[1], atom_types[2], atom_types[3], ampl1, ampl2, ampl3, ampl4, phase1, phase2, phase3, phase4) except ValueError as err: raise CoarseGrainError( 'Unexpected format in Coarse Grain dihedrals. Expected ' 'different data type: %s. See format specification' % err) except IndexError: raise CoarseGrainError( 'Unexpected format in Coarse Grain parameter file. ' 'Expected more data fields on the line. See format ' 'specification') file_line = params.readline() continue warnings.warn('Line (%s) ignored in Coarse Grain parameter file' % file_line) file_line = params.readline() # End reading while # Now let's add comments to our topology file! parm.parm_comments['ANGLE_FORCE_CONSTANT'].append( 'This section is only an index and is not used for Coarse ' 'grained topologies') parm.parm_comments['DIHEDRAL_FORCE_CONSTANT'].append( 'This section is only an index and is not used for Coarse ' 'grained topologies') parm.parm_comments['DIHEDRAL_PERIODICITY'].append( 'This section is not used for Coarse grained topologies') parm.parm_comments['DIHEDRAL_PHASE'].append( 'This section is not used for Coarse grained topologies') # Now let's add our new sections parm.add_flag('ANGLE_COEF_A','5E16.8',parm.ptr('numang'), comments='A Coefficient for Coarse grained force field') parm.add_flag('ANGLE_COEF_B','5E16.8',parm.ptr('numang'), comments='B Coefficient for Coarse grained force field') parm.add_flag('ANGLE_COEF_C','5E16.8',parm.ptr('numang'), comments='C Coefficient for Coarse grained force field') parm.add_flag('ANGLE_COEF_D','5E16.8',parm.ptr('numang'), comments='D Coefficient for Coarse grained force field') parm.add_flag('DIHEDRAL_AMPLITUDE_1','5E16.8',parm.ptr('nptra'), comments='1st Dihedral Amplitude for coarse grained force ' 'field') parm.add_flag('DIHEDRAL_AMPLITUDE_2','5E16.8',parm.ptr('nptra'), comments='2nd Dihedral Amplitude for coarse grained force ' 'field' ) parm.add_flag('DIHEDRAL_AMPLITUDE_3','5E16.8',parm.ptr('nptra'), comments='3rd Dihedral Amplitude for coarse grained force ' 'field' ) parm.add_flag('DIHEDRAL_AMPLITUDE_4','5E16.8',parm.ptr('nptra'), comments='4th Dihedral Amplitude for coarse grained force ' 'field' ) parm.add_flag('DIHEDRAL_PHASE_1','5E16.8',parm.ptr('nptra'), comments='1st Dihedral Phase for coarse grained force field') parm.add_flag('DIHEDRAL_PHASE_2','5E16.8',parm.ptr('nptra'), comments='2nd Dihedral Phase for coarse grained force field') parm.add_flag('DIHEDRAL_PHASE_3','5E16.8',parm.ptr('nptra'), comments='3rd Dihedral Phase for coarse grained force field') parm.add_flag('DIHEDRAL_PHASE_4','5E16.8',parm.ptr('nptra'), comments='4th Dihedral Phase for coarse grained force field') for i in range(len(parm.parm_data['ANGLE_FORCE_CONSTANT'])): try: index = int(parm.parm_data['ANGLE_FORCE_CONSTANT'][i]) angl = angle_params[index] except KeyError: raise CoarseGrainError('Missing angle parameters for index %d' % int(index)) parm.parm_data['ANGLE_COEF_A'][i] = angl.acoef parm.parm_data['ANGLE_COEF_B'][i] = angl.bcoef parm.parm_data['ANGLE_COEF_C'][i] = angl.ccoef parm.parm_data['ANGLE_COEF_D'][i] = angl.dcoef for i in range(len(parm.parm_data['DIHEDRAL_FORCE_CONSTANT'])): try: index = int(parm.parm_data['DIHEDRAL_FORCE_CONSTANT'][i]) dihe = dihedral_params[index] except KeyError: raise CoarseGrainError('Missing dihedral parameters for index %d' % int(index)) parm.parm_data['DIHEDRAL_AMPLITUDE_1'][i] = dihe.ampl1 parm.parm_data['DIHEDRAL_AMPLITUDE_2'][i] = dihe.ampl2 parm.parm_data['DIHEDRAL_AMPLITUDE_3'][i] = dihe.ampl3 parm.parm_data['DIHEDRAL_AMPLITUDE_4'][i] = dihe.ampl4 parm.parm_data['DIHEDRAL_PHASE_1'][i] = dihe.phase1 parm.parm_data['DIHEDRAL_PHASE_2'][i] = dihe.phase2 parm.parm_data['DIHEDRAL_PHASE_3'][i] = dihe.phase3 parm.parm_data['DIHEDRAL_PHASE_4'][i] = dihe.phase4