ICODE-ADC/lammps
4
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

git-svn-id: svn://svn.icms.temple.edu/lammps-ro/trunk@10174 f3b2605a-c512-4ea7-a41b-209d697bcdaa

Browse Source
This commit is contained in:
sjplimp
2013-06-28 20:32:49 +00:00
parent 916c06cdbc
commit 191b83cd92
26 changed files with 17657 additions and 0 deletions
Download Patch File Download Diff File Expand all files Collapse all files

119
tools/moltemplate/src/extract_lammps_data.py Executable file
View File

@ -0,0 +1,119 @@
#!/usr/bin/env python
lammps_data_sections = set(['Atoms',
'Masses',
'Bonds',
'Bond Coeffs',
'Angles',
'Angle Coeffs',
'Dihedrals',
'Dihedral Coeffs',
'Impropers',
'Improper Coeffs',
'BondBond Coeffs', # class2 angles
'BondAngle Coeffs', # class2 angles
'MiddleBondTorsion Coeffs', # class2 dihedrals
'EndBondTorsion Coeffs', # class2 dihedrals
'AngleTorsion Coeffs', # class2 dihedrals
'AngleAngleTorsion Coeffs', # class2 dihedrals
'BondBond13 Coeffs', # class2 dihedrals
'AngleAngle Coeffs', # class2 impropers
'Angles By Type', # new. not standard LAMMPS
'Dihedrals By Type',# new. not standard LAMMPS
'Angles By Type']) # new. not standard LAMMPS
def DeleteComments(string,
escape='\\',
comment_char='#'):
escaped_state = False
for i in range(0,len(string)):
if string[i] in escape:
if escaped_state:
escaped_state = False
else:
escaped_state = True
elif string[i] == comment_char:
if not escaped_state:
return string[0:i]
return string
def ExtractDataSection(f,
section_name,
comment_char = '#',
include_section_name = False,
return_line_nums = False):
inside_section = False
if section_name in ('header','Header'): #"Header" section includes beginning
inside_section = True
nonblank_encountered = False
nonheader_encountered = False
i = 0
for line_orig in f:
return_this_line = False
line = DeleteComments(line_orig).strip()
if line in lammps_data_sections:
nonheader_encountered = True
if section_name in ('header', 'Header'):
# The "header" section includes all lines at the beginning of the
# before any other section is encountered.
if nonheader_encountered:
return_this_line = False
else:
return_this_line = True
elif line == section_name:
inside_section = True
nonblank_encountered = False
if include_section_name:
return_this_line = True
# A block of blank lines (which dont immediately follow
# the section_name) signal the end of a section:
elif len(line) == 0:
if inside_section and include_section_name:
return_this_line = True
if nonblank_encountered:
inside_section = False
elif line[0] != comment_char:
if inside_section:
nonblank_encountered = True
return_this_line = True
if return_this_line:
if return_line_nums:
yield i
else:
yield line_orig
i += 1
if __name__ == "__main__":
import sys
lines = sys.stdin.readlines()
exclude_sections = False
if sys.argv[1] == '-n':
exclude_sections = True
del sys.argv[1]
if not exclude_sections:
for section_name in sys.argv[1:]:
for line in ExtractDataSection(lines, section_name):
sys.stdout.write(line)
else:
line_nums_exclude = set([])
for section_name in sys.argv[1:]:
for line_num in ExtractDataSection(lines,
section_name,
include_section_name=True,
return_line_nums=True):
line_nums_exclude.add(line_num)
for i in range(0, len(lines)):
if i not in line_nums_exclude:
sys.stdout.write(lines[i])

2036
tools/moltemplate/src/ltemplify.py Executable file
View File

File diff suppressed because it is too large Load Diff

745
tools/moltemplate/src/lttree.py Executable file
View File

@ -0,0 +1,745 @@
#!/usr/bin/env python
# Author: Andrew Jewett (jewett.aij at g mail)
# http://www.chem.ucsb.edu/~sheagroup
# License: 3-clause BSD License (See LICENSE.TXT)
# Copyright (c) 2011, Regents of the University of California
# All rights reserved.
"""
lttree.py
lttree.py is an extension of the generic ttree.py program.
This version can understand and manipulate ttree-style templates which
are specialized for storing molecule-specific data for use in LAMMPS.
The main difference between lttree.py and ttree.py is:
Unlike ttree.py, lttree.py understands rigid-body movement commands like
"rot()" and "move()" which allows it to reorient and move each copy
of a molecule to a new location. (ttree.py just ignores these commands.
Consequently LAMMPS input file (fragments) created with ttree.py have
invalid (overlapping) atomic coordinates and must be modified or aguemted
later (by loading atomic coordinates from a PDB file or an XYZ file).
lttree.py understands the "Data Atoms" section of a LAMMPS
data file (in addition to the various "atom_styles" which effect it).
Additional LAMMPS-specific features may be added in the future.
"""
import sys
from ttree import *
from lttree_styles import *
from ttree_matrix_stack import *
try:
unicode
except NameError:
# Python 3
basestring = unicode = str
class LttreeSettings(BasicUISettings):
def __init__(self,
user_bindings_x=None,
user_bindings=None,
order_method='by_command'):
BasicUISettings.__init__(self,
user_bindings_x,
user_bindings,
order_method)
# The following new member data indicate which columns store
# LAMMPS-specific information.
# The next 6 members store keep track of the different columns
# of the "Data Atoms" section of a LAMMPS data file:
self.column_names = [] #<--A list of column names (optional)
self.ii_coords=[] #<--A list of triplets of column indexes storing coordinate data
self.ii_vects=[] #<--A list of triplets of column indexes storing directional data
# (such as dipole or ellipsoid orientations)
self.i_atomid=None #<--An integer indicating which column has the atomid
self.i_atomtype=None #<--An integer indicating which column has the atomtype
self.i_molid=None #<--An integer indicating which column has the molid, if applicable
self.infile=None # Name of the outermost file. This is the file
# which was read at the moment parsing begins.
def LttreeParseArgs(argv, settings):
BasicUIParseArgs(argv, settings)
# Loop over the remaining arguments not processed yet.
# These arguments are specific to the lttree.py program
# and are not understood by ttree.py:
i = 1
while i < len(argv):
#sys.stderr.write('argv['+str(i)+'] = \"'+argv[i]+'\"\n')
if ((argv[i].lower() == '-atomstyle') or
(argv[i].lower() == '-atom-style') or
(argv[i].lower() == '-atom_style')):
if i+1 >= len(argv):
raise InputError('Error('+g_program_name+'): The '+argv[i]+' flag should be followed by a LAMMPS\n'
' atom_style name (or single quoted string containing a space-separated\n'
' list of column names such as: atom-ID atom-type q x y z molecule-ID.)\n')
settings.column_names = AtomStyle2ColNames(argv[i+1])
sys.stderr.write('\n \"'+data_atoms+'\" column format:\n')
sys.stderr.write(' '+(' '.join(settings.column_names))+'\n\n')
settings.ii_coords = ColNames2Coords(settings.column_names)
settings.ii_vects = ColNames2Vects(settings.column_names)
settings.i_atomid, settings.i_atomtype, settings.i_molid = ColNames2AidAtypeMolid(settings.column_names)
del(argv[i:i+2])
elif (argv[i].lower() == '-icoord'):
if i+1 >= len(argv):
raise InputError('Error: '+argv[i]+' flag should be followed by list of integers\n'
' corresponding to column numbers for coordinates in\n'
' the \"'+data_atoms+'\" section of a LAMMPS data file.\n')
ilist = argv[i+1].split()
if (len(ilist) % 3) != 0:
raise InputError('Error: '+argv[i]+' flag should be followed by list of integers.\n'
' This is usually a list of 3 integers, but it can contain more.\n'
' The number of cooridnate columns must be divisible by 3,\n'
' (even if the simulation is in 2 dimensions)\n')
settings.iaffinevects = []
for i in range(0, len(ilist)/3):
cols = [int(ilist[3*i])+1,
int(ilist[3*i+1])+1,
int(ilist[3*i+2])+1]
settings.iaffinevects.append(cols)
del(argv[i:i+2])
elif (argv[i].lower() == '-ivect'):
if i+1 >= len(argv):
raise InputError('Error: '+argv[i]+' flag should be followed by list of integers\n'
' corresponding to column numbers for direction vectors in\n'
' the \"'+data_atoms+'\" section of a LAMMPS data file.\n')
ilist = argv[i+1].split()
if (len(ilist) % 3) != 0:
raise InputError('Error: '+argv[i]+' flag should be followed by list of integers.\n'
' This is usually a list of 3 integers, but it can contain more.\n'
' The number of cooridnate columns must be divisible by 3,\n'
' (even if the simulation is in 2 dimensions)\n')
settings.ivects = []
for i in range(0, len(ilist)/3):
cols = [int(ilist[3*i])+1,
int(ilist[3*i+1])+1,
int(ilist[3*i+2])+1]
settings.ivects.append(cols)
del(argv[i:i+2])
elif ((argv[i].lower() == '-iatomid') or
(argv[i].lower() == '-iid') or
(argv[i].lower() == '-iatom-id')):
if ((i+1 >= len(argv)) or (not str.isdigit(argv[i+1]))):
raise InputError('Error: '+argv[i]+' flag should be followed by an integer\n'
' (>=1) indicating which column in the \"'+data_atoms+'\" section of a\n'
' LAMMPS data file contains the atom id number (typically 1).\n'
' (This argument is unnecessary if you use the -atomstyle argument.)\n')
i_atomid = int(argv[i+1])-1
del(argv[i:i+2])
elif ((argv[i].lower() == '-iatomtype') or
(argv[i].lower() == '-itype') or
(argv[i].lower() == '-iatom-type')):
if ((i+1 >= len(argv)) or (not str.isdigit(argv[i+1]))):
raise InputError('Error: '+argv[i]+' flag should be followed by an integer\n'
' (>=1) indicating which column in the \"'+data_atoms+'\" section of a\n'
' LAMMPS data file contains the atom type.\n'
' (This argument is unnecessary if you use the -atomstyle argument.)\n')
i_atomtype = int(argv[i+1])-1
del(argv[i:i+2])
elif ((argv[i].lower() == '-imolid') or
(argv[i].lower() == '-imol') or
(argv[i].lower() == '-imol-id') or
(argv[i].lower() == '-imoleculeid') or
(argv[i].lower() == '-imolecule-id')):
if ((i+1 >= len(argv)) or (not str.isdigit(argv[i+1]))):
raise InputError('Error: '+argv[i]+' flag should be followed by an integer\n'
' (>=1) indicating which column in the \"'+data_atoms+'\" section of a\n'
' LAMMPS data file contains the molecule id number.\n'
' (This argument is unnecessary if you use the -atomstyle argument.)\n')
i_molid = int(argv[i+1])-1
del(argv[i:i+2])
elif ((argv[i][0] == '-') and (__name__ == "__main__")):
#elif (__name__ == "__main__"):
raise InputError('Error('+g_program_name+'):\n'
'Unrecogized command line argument \"'+argv[i]+'\"\n')
else:
i += 1
if __name__ == "__main__":
# Instantiate the lexer we will be using.
# (The lexer's __init__() function requires an openned file.
# Assuming __name__ == "__main__", then the name of that file should
# be the last remaining (unprocessed) argument in the argument list.
# Otherwise, then name of that file will be determined later by the
# python script which imports this module, so we let them handle it.)
if len(argv) == 1:
raise InputError('Error: This program requires at least one argument\n'
' the name of a file containing ttree template commands\n')
elif len(argv) == 2:
try:
settings.lex = TemplateLexer(open(argv[1], 'r'), argv[1]) # Parse text from file
except IOError:
sys.stderr.write('Error: unable to open file\n'
' \"'+argv[1]+'\"\n'
' for reading.\n')
sys.exit(1)
del(argv[1:2])
else:
# if there are more than 2 remaining arguments,
problem_args = ['\"'+arg+'\"' for arg in argv[1:]]
raise InputError('Syntax Error('+g_program_name+'):\n\n'
' Problem with argument list.\n'
' The remaining arguments are:\n\n'
' '+(' '.join(problem_args))+'\n\n'
' (The actual problem may be earlier in the argument list.\n'
' If these arguments are source files, then keep in mind\n'
' that this program can not parse multiple source files.)\n'
' Check the syntax of the entire argument list.\n')
if len(settings.ii_coords) == 0:
sys.stderr.write('########################################################\n'
'## WARNING: atom_style unspecified ##\n'
'## --> \"'+data_atoms+'\" column data has an unknown format ##\n'
'## Assuming atom_style = \"full\" ##\n'
# '########################################################\n'
# '## To specify the \"'+data_atoms+'\" column format you can: ##\n'
# '## 1) Use the -atomstyle \"STYLE\" argument ##\n'
# '## where \"STYLE\" is a string indicating a LAMMPS ##\n'
# '## atom_style, including hybrid styles.(Standard ##\n'
# '## atom styles defined in 2011 are supported.) ##\n'
# '## 2) Use the -atomstyle \"COL_LIST\" argument ##\n'
# '## where \"COL_LIST" is a quoted list of strings ##\n'
# '## indicating the name of each column. ##\n'
# '## Names \"x\",\"y\",\"z\" are interpreted as ##\n'
# '## atomic coordinates. \"mux\",\"muy\",\"muz\" ##\n'
# '## are interpreted as direction vectors. ##\n'
# '## 3) Use the -icoord \"cx cy cz...\" argument ##\n'
# '## where \"cx cy cz\" is a list of integers ##\n'
# '## indicating the column numbers for the x,y,z ##\n'
# '## coordinates of each atom. ##\n'
# '## 4) Use the -ivect \"cmux cmuy cmuz...\" argument ##\n'
# '## where \"cmux cmuy cmuz...\" is a list of ##\n'
# '## integers indicating the column numbers for ##\n'
# '## the vector that determines the direction of a ##\n'
# '## dipole or ellipsoid (ie. a rotateable vector).##\n'
# '## (More than one triplet can be specified. The ##\n'
# '## number of entries must be divisible by 3.) ##\n'
'########################################################\n')
# The default atom_style is "full"
settings.column_names = AtomStyle2ColNames('full')
settings.ii_coords = ColNames2Coords(settings.column_names)
settings.ii_vects = ColNames2Vects(settings.column_names)
settings.i_atomid, settings.i_atomtype, settings.i_molid = ColNames2AidAtypeMolid(settings.column_names)
def TransformAtomText(text, matrix):
""" Apply transformations to the coordinates and other vector degrees
of freedom stored in the \"Data Atoms\" section of a LAMMPS data file.
This is the \"text\" argument.
The \"matrix\" stores the aggregate sum of combined transformations
to be applied.
"""
#sys.stderr.write('matrix_stack.M = \n'+ MatToStr(matrix) + '\n')
lines = text.split('\n')
for i in range(0, len(lines)):
line_orig = lines[i]
ic = line_orig.find('#')
if ic != -1:
line = line_orig[:ic]
comment = ' '+line_orig[ic:].rstrip('\n')
else:
line = line_orig.rstrip('\n')
comment = ''
columns = line.split()
if len(columns) > 0:
if len(columns) == len(settings.column_names)+3:
raise InputError('Error: lttree.py does not yet support integer unit-cell counters \n'
' within the \"'+data_atoms+'\" section of a LAMMPS data file.\n'
' Instead please add the appropriate offsets (these offsets\n'
' should be multiples of the cell size) to the atom coordinates\n'
' in the data file, and eliminate the extra columns. Then try again.\n'
' (If you get this message often, email me and I\'ll fix this limitation.)')
if len(columns) < len(settings.column_names):
raise InputError('Error: The number of columns in your data file does not\n'
' match the LAMMPS atom_style you selected.\n'
' Use the -atomstyle <style> command line argument.\n')
x0 = [0.0, 0.0, 0.0]
x = [0.0, 0.0, 0.0]
# Atomic coordinates transform using "affine" transformations
# (translations plus rotations [or other linear transformations])
for cxcycz in settings.ii_coords:
for d in range(0,3):
x0[d] = float(columns[cxcycz[d]])
AffineTransform(x, matrix, x0) # x = matrix * x0 + b
for d in range(0,3): #("b" is part of "matrix")
columns[cxcycz[d]] = str(x[d])
# Dipole moments and other direction-vectors
# are not effected by translational movement
for cxcycz in settings.ii_vects:
for d in range(0,3):
x0[d] = float(columns[cxcycz[d]])
LinearTransform(x, matrix, x0) # x = matrix * x0
for d in range(0,3):
columns[cxcycz[d]] = str(x[d])
lines[i] = ' '.join(columns) + comment
return '\n'.join(lines)
def CalcCM(text_Atoms,
text_Masses=None,
settings=None):
types2masses = None
# Loop through the "Masses" section: what is the mass of each atom type?
if text_Masses != None:
types2masses = {}
lines = text_Masses.split('\n')
for i in range(0, len(lines)):
line = lines[i]
columns = line.split()
if len(columns) == 2:
atomtype = columns[0]
m = float(columns[1])
types2masses[atomtype] = m
lines = text_Atoms.split('\n')
# Pass 1 through the "Data Atoms" section: Determine each atom's mass
if text_Masses != None:
assert(settings != None)
for i in range(0, len(lines)):
line = lines[i]
columns = line.split()
atomid = columns[settings.i_atomid]
atomtype = columns[settings.i_atomtype]
if atomtype not in types2masses[atomtype]:
raise InputError('Error(lttree): You have neglected to define the mass of atom type: \"'+atomtype+'\"\n'
'Did you specify the mass of every atom type using write(\"Masses\"){}?')
atomid2mass[atomid] = atomtype2mass[atomtype]
# Pass 2 through the "Data Atoms" section: Find the center of mass.
for i in range(0, len(lines)):
line = lines[i]
columns = line.split()
if len(columns) > 0:
if len(columns) == len(settings.column_names)+3:
raise InputError('Error: lttree.py does not yet support integer unit-cell counters (ix, iy, iz)\n'
' within the \"'+data_atoms+'\" section of a LAMMPS data file.\n'
' Instead please add the appropriate offsets (these offsets\n'
' should be multiples of the cell size) to the atom coordinates\n'
' in the data file, and eliminate the extra columns. Then try again.\n'
' (If you get this message often, email me and I\'ll fix this limitation.)')
if len(columns) != len(settings.column_names):
raise InputError('Error: The number of columns in your data file does not\n'
' match the LAMMPS atom_style you selected.\n'
' Use the -atomstyle <style> command line argument.\n')
x = [0.0, 0.0, 0.0]
if atomids2masses != None:
m = atomids2masses[atomid]
else:
m = 1.0
tot_m += m
for cxcycz in settings.ii_coords:
for d in range(0,3):
x[d] = float(columns[cxcycz[d]])
tot_x[d] += x[d]
# Note: dipole moments and other direction vectors don't effect
# the center of mass. So I commented out the loop below.
#for cxcycz in settings.ii_vects:
# for d in range(0,3):
# v[d] = float(columns[cxcycz[d]])
lines[i] = ' '.join(columns)
xcm = [0.0, 0.0, 0.0]
for d in range(0,3):
xcm[d] = tot_x[d] / tot_m
return xcm
def _ExecCommands(command_list,
index,
global_files_content,
settings,
matrix_stack,
current_scope_id=None,
substitute_vars=True):
"""
_ExecCommands():
The argument "commands" is a nested list of lists of
"Command" data structures (defined in ttree.py).
Carry out the write() and write_once() commands (which
write out the contents of the templates contain inside them).
Instead of writing the files, save their contents in a string.
The argument "global_files_content" should be of type defaultdict(list)
It is an associative array whose key is a string (a filename)
and whose value is a lists of strings (of rendered templates).
"""
files_content = defaultdict(list)
postprocessing_commands = []
while index < len(command_list):
command = command_list[index]
index += 1
# For debugging only
if ((not isinstance(command, StackableCommand)) and
(not isinstance(command, ScopeCommand)) and
(not isinstance(command, WriteFileCommand))):
sys.stderr.write(str(command)+'\n')
if isinstance(command, PopCommand):
assert(current_scope_id != None)
if command.context_node == None:
command.context_node = current_scope_id
if isinstance(command, PopRightCommand):
matrix_stack.PopRight(which_stack = command.context_node)
elif isinstance(command, PopLeftCommand):
matrix_stack.PopLeft(which_stack = command.context_node)
else:
assert(False)
elif isinstance(command, PushCommand):
assert(current_scope_id != None)
if command.context_node == None:
command.context_node = current_scope_id
# Some commands are post-processing commands, and must be
# carried out AFTER all the text has been rendered. For example
# the "movecm(0,0,0)" waits until all of the coordinates have
# been rendered, calculates the center-of-mass, and then applies
# a translation moving the center of mass to the origin (0,0,0).
# We need to figure out which of these commands need to be
# postponed, and which commands can be carried out now.
# ("now"=pushing transformation matrices onto the matrix stack).
# UNFORTUNATELY POSTPONING SOME COMMANDS MAKES THE CODE UGLY
transform_list = command.contents.split('.')
transform_blocks = []
i_post_process = -1
# Example: Suppose:
#command.contents = '.rot(30,0,0,1).movecm(0,0,0).rot(45,1,0,0).scalecm(2.0).move(-2,1,0)'
# then
#transform_list = ['rot(30,0,0,1)', 'movecm(0,0,0)', 'rot(45,1,0,0)', 'scalecm(2.0)', 'move(-2,1,0)']
# Note: the first command 'rot(30,0,0,1)' is carried out now.
# The remaining commands are carried out during post-processing,
# (when processing the "ScopeEnd" command.
#
# We break up the commands into "blocks" separated by center-
# of-mass transformations ('movecm', 'rotcm', or 'scalecm')
#
# transform_blocks = ['.rot(30,0,0,1)',
# '.movecm(0,0,0).rot(45,1,0,0)',
# '.scalecm(2.0).move(-2,1,0)']
i = 0
while i < len(transform_list):
transform_block = ''
while i < len(transform_list):
transform = transform_list[i]
i += 1
if transform != '':
transform_block += '.' + transform
transform = transform.split('(')[0]
if ((transform == 'movecm') or
(transform == 'rotcm') or
(transform == 'scalecm')):
break
transform_blocks.append(transform_block)
if len(postprocessing_commands) == 0:
# The first block (before movecm, rotcm, or scalecm)
# can be executed now by modifying the matrix stack.
if isinstance(command, PushRightCommand):
matrix_stack.PushCommandsRight(transform_blocks[0].strip('.'),
command.srcloc,
which_stack=command.context_node)
elif isinstance(command, PushLeftCommand):
matrix_stack.PushCommandsLeft(transform_blocks[0].strip('.'),
command.srcloc,
which_stack=command.context_node)
# Everything else must be saved for later.
postprocessing_blocks = transform_blocks[1:]
else:
# If we already encountered a "movecm" "rotcm" or "scalecm"
# then all of the command blocks must be handled during
# postprocessing.
postprocessing_blocks = transform_blocks
for transform_block in postprocessing_blocks:
assert(isinstance(block, basestring))
if isinstance(command, PushRightCommand):
postprocessing_commands.append(PushRightCommand(transform_block,
command.srcloc,
command.context_node))
elif isinstance(command, PushLeftCommand):
postprocessing_commands.append(PushLeftCommand(transform_block,
command.srcloc,
command.context_node))
elif isinstance(command, WriteFileCommand):
# --- Throw away lines containin references to deleted variables:---
# First: To edit the content of a template,
# you need to make a deep local copy of it
tmpl_list = []
for entry in command.tmpl_list:
if isinstance(entry, TextBlock):
tmpl_list.append(TextBlock(entry.text,
entry.srcloc)) #, entry.srcloc_end))
else:
tmpl_list.append(entry)
# Now throw away lines with deleted variables
DeleteLinesWithBadVars(tmpl_list)
# --- Now render the text ---
text = Render(tmpl_list,
substitute_vars)
# ---- Coordinates of the atoms, must be rotated
# and translated after rendering.
# In addition, other vectors (dipoles, ellipsoid orientations)
# must be processed.
# This requires us to re-parse the contents of this text
# (after it has been rendered), and apply these transformations
# before passing them on to the caller.
if command.filename == data_atoms:
text = TransformAtomText(text, matrix_stack.M)
files_content[command.filename].append(text)
elif isinstance(command, ScopeBegin):
if isinstance(command.node, InstanceObj):
if ((command.node.children != None) and
(len(command.node.children) > 0)):
matrix_stack.PushStack(command.node)
# "command_list" is a long list of commands.
# ScopeBegin and ScopeEnd are (usually) used to demarcate/enclose
# the commands which are issued for a single class or
# class instance. _ExecCommands() carries out the commands for
# a single class/instance. If we reach a ScopeBegin(),
# then recursively process the commands belonging to the child.
index = _ExecCommands(command_list,
index,
files_content,
settings,
matrix_stack,
command.node,
substitute_vars)
elif isinstance(command, ScopeEnd):
if data_atoms in files_content:
for ppcommand in postprocessing_commands:
if data_masses in files_content:
xcm = CalcCM(files_content[data_atoms],
files_content[data_masses],
settings)
else:
xcm = CalcCM(files_content[data_atoms])
if isinstance(ppcommand, PushRightCommand):
matrix_stack.PushCommandsRight(ppcommand.contents,
ppcommand.srcloc,
xcm,
which_stack=command.context_node)
elif isinstance(ppcommand, PushLeftCommand):
matrix_stack.PushCommandsLeft(ppcommand.contents,
ppcommand.srcloc,
xcm,
which_stack=command.context_node)
files_content[data_atoms] = \
TransformAtomText(files_content[data_atoms],
matrix_stack.M)
for ppcommand in postprocessing_commands:
matrix_stack.Pop(which_stack = command.context_node)
#(same as PopRight())
if isinstance(command.node, InstanceObj):
if ((command.node.children != None) and
(len(command.node.children) > 0)):
matrix_stack.PopStack()
# "ScopeEnd" means we're done with this class/instance.
break
else:
assert(False)
# no other command types allowed at this point
# After processing the commands in this list,
# merge the templates with the callers template list
for filename, tmpl_list in files_content.items():
global_files_content[filename] += \
files_content[filename]
return index
def ExecCommands(commands,
files_content,
settings,
substitute_vars=True):
matrix_stack = MultiAffineStack()
index = _ExecCommands(commands,
0,
files_content,
settings,
matrix_stack,
None,
substitute_vars)
assert(index == len(commands))
def WriteFiles(files_content, suffix='', write_to_stdout=True):
for filename, str_list in files_content.items():
if filename != None:
out_file = None
if filename == '':
if write_to_stdout:
out_file = sys.stdout
else:
out_file = open(filename+suffix, 'a')
if out_file != None:
out_file.write(''.join(str_list))
if filename != '':
out_file.close()
if __name__ == "__main__":
"""
This is is a "main module" wrapper for invoking lttree.py
as a stand alone program. This program:
1)reads a ttree file,
2)constructs a tree of class definitions (g_objectdefs)
3)constructs a tree of instantiated class objects (g_objects),
4)automatically assigns values to the variables,
5)and carries out the "write" commands to write the templates a file(s).
"""
g_program_name = __file__.split('/')[-1] # ='lttree.py'
g_date_str = '2013-2-15'
g_version_str = '0.71'
####### Main Code Below: #######
sys.stderr.write(g_program_name+' v'+g_version_str+' '+g_date_str+' ')
sys.stderr.write('\n(python version '+str(sys.version)+')\n')
if sys.version < '2.7':
raise InputError('Error: Alas, you must upgrade to a newever version of python.')
try:
#settings = BasicUISettings()
#BasicUIParseArgs(sys.argv, settings)
settings = LttreeSettings()
LttreeParseArgs(sys.argv, settings)
# Data structures to store the class definitionss and instances
g_objectdefs = StaticObj('', None) # The root of the static tree
# has name '' (equivalent to '/')
g_objects = InstanceObj('', None) # The root of the instance tree
# has name '' (equivalent to '/')
# A list of commands to carry out
g_static_commands = []
g_instance_commands = []
BasicUI(settings,
g_objectdefs,
g_objects,
g_static_commands,
g_instance_commands)
# Interpret the the commands. (These are typically write() or
# write_once() commands, rendering templates into text.
# This step also handles coordinate transformations and delete commands.
# Coordinate transformations can be applied to the rendered text
# as a post-processing step.
sys.stderr.write(' done\nbuilding templates...')
files_content = defaultdict(list)
ExecCommands(g_static_commands,
files_content,
settings,
False)
ExecCommands(g_instance_commands,
files_content,
settings,
False)
# Finally: write the rendered text to actual files.
# Erase the files that will be written to:
sys.stderr.write(' done\nwriting templates...')
EraseTemplateFiles(g_static_commands)
EraseTemplateFiles(g_instance_commands)
# Write the files as templates
# (with the original variable names present)
WriteFiles(files_content, suffix=".template", write_to_stdout=False)
# Write the files with the variables substituted by values
sys.stderr.write(' done\nbuilding and rendering templates...')
files_content = defaultdict(list)
ExecCommands(g_static_commands, files_content, settings, True)
ExecCommands(g_instance_commands, files_content, settings, True)
sys.stderr.write(' done\nwriting rendered templates...\n')
WriteFiles(files_content)
sys.stderr.write(' done\n')
# Step 11: Now write the variable bindings/assignments table.
# Now write the variable bindings/assignments table.
sys.stderr.write('writing \"ttree_assignments.txt\" file...')
open('ttree_assignments.txt', 'w').close() # <-- erase previous version.
WriteVarBindingsFile(g_objectdefs)
WriteVarBindingsFile(g_objects)
sys.stderr.write(' done\n')
except (ValueError, InputError) as err:
sys.stderr.write('\n\n'+str(err)+'\n')
sys.exit(-1)

2099
tools/moltemplate/src/lttree_check.py Executable file
View File

File diff suppressed because it is too large Load Diff

335
tools/moltemplate/src/lttree_postprocess.py Executable file
View File

@ -0,0 +1,335 @@
#!/usr/bin/env python
"""
lttree_postprocess.py
This is a stand-alone python script which checks the files created by
lttree.py to insure that the standard instance-variables ($variables)
have all been defined. This script performs a task which is very similar
to the task performed by lttree_check.py. This script attempts to detect
mistakes in the names of $atom, $bond, $angle, $dihedral, $improper, & $mol
variables.
"""
import sys
from lttree_styles import *
from ttree_lex import ExtractCatName
g_program_name = __file__.split('/')[-1] # = 'lttree_postprocess.py'
g_version_str = '0.4'
g_date_str = '2012-12-12'
atom_style = 'full'
ttree_assignments_fname = 'ttree_assignments.txt'
defined_mols = set([])
defined_atoms = set([])
defined_bonds = set([])
defined_angles = set([])
defined_dihedrals = set([])
defined_impropers = set([])
g_no_check_msg = \
'(To override this error, run moltemplate using the \"-nocheck\" argument.)\n'
if len(sys.argv) > 1:
for i in range(0,len(sys.argv)):
if ((sys.argv[i].lower() == '-atomstyle') or
(sys.argv[i].lower() == '-atom-style') or
(sys.argv[i].lower() == '-atom_style')):
if i+1 >= len(sys.argv):
raise InputError('Error('+g_program_name+'): The '+sys.argv[i]+' flag should be followed by a LAMMPS\n'
' atom_style name (or single quoted string containing a space-separated\n'
' list of column names such as: atom-ID atom-type q x y z molecule-ID.)\n')
atom_style = sys.argv[i+1]
elif ((sys.argv[i].lower() == '-ttreeassignments') or
(sys.argv[i].lower() == '-ttree-assignments') or
(sys.argv[i].lower() == '-ttree_assignments')):
if i+1 >= len(sys.argv):
raise InputError('Error('+g_program_name+'): The '+sys.argv[i]+' flag should be followed by \n'
' a file containing the variable bindings created by ttree/moltemplate.\n')
ttree_assignments_fname = sys.argv[i+1]
else:
pass # ignore other arguments (they are intended for lttree.py)
atom_column_names = AtomStyle2ColNames(atom_style)
i_atomid = 0
i_molid = -1
for i in range(0,len(atom_column_names)):
if atom_column_names[i].lower() == 'atom-id':
i_atomid = i
elif atom_column_names[i].lower() == 'molecule-id':
i_molid = i
i_max_column = max(i_atomid, i_molid)
# The following variables are defined in "lttree_styles.py"
#data_atoms="Data Atoms"
#data_masses="Data Masses"
#data_velocities="Data Velocities"
#data_bonds="Data Bonds"
#data_angles="Data Angles"
#data_dihedrals="Data Dihedrals"
#data_impropers="Data Impropers"
sys.stderr.write(g_program_name+' v'+g_version_str+' '+g_date_str+'\n')
try:
# ------------ defined_atoms ------------
try:
f = open(data_atoms+'.template', 'r')
except:
raise InputError('Error('+g_program_name+'): Unable to open file\n'+
'\"'+data_atoms+'.template\"\n'
' for reading. (Do your files lack a \"'+data_atoms+'\" section?)\n'
+g_no_check_msg+'\n')
for line_orig in f:
ic = line_orig.find('#')
if ic != -1:
line = line_orig[:ic]
else:
line = line_orig.rstrip('\n')
tokens = line.strip().split()
if len(tokens) == 0:
pass
elif len(tokens) <= i_max_column:
raise InputError('Error('+g_program_name+'): The following line from\n'
' "\"'+data_atoms+'.template\" has bad format:\n\n'
+line_orig+'\n'
' This my probably an internal error. (Feel free to contact the developer.)\n'
+g_no_check_msg+'\n')
else:
defined_atoms.add(tokens[i_atomid])
if i_molid != -1:
defined_mols.add(tokens[i_molid])
f.close()
# ------------ defined_bonds ------------
try:
f = open(data_bonds+'.template', 'r')
for line_orig in f:
ic = line_orig.find('#')
if ic != -1:
line = line_orig[:ic]
else:
line = line_orig.rstrip('\n')
tokens = line.strip().split()
if len(tokens) == 0:
pass
elif len(tokens) < 4:
raise InputError('Error('+g_program_name+'): The following line from\n'
' "\"'+data_bonds+'.template\" has bad format:\n\n'
+line_orig+'\n'
' This my probably an internal error. (Feel free to contact the developer.)\n'
+g_no_check_msg+'\n')
else:
defined_bonds.add(tokens[0])
f.close()
except:
pass # Defining bonds (stored in the data_bonds file) is optional
# ------------ defined_angles ------------
try:
f = open(data_angles+'.template', 'r')
for line_orig in f:
ic = line_orig.find('#')
if ic != -1:
line = line_orig[:ic]
else:
line = line_orig.rstrip('\n')
tokens = line.strip().split()
if len(tokens) == 0:
pass
elif len(tokens) < 5:
raise InputError('Error('+g_program_name+'): The following line from\n'
' "\"'+data_angles+'.template\" has bad format:\n\n'
+line_orig+'\n'
' This my probably an internal error. (Feel free to contact the developer.)\n'
+g_no_check_msg+'\n')
else:
defined_angles.add(tokens[0])
f.close()
except:
pass # Defining angles (stored in the data_angles file) is optional
# ------------ defined_dihedrals ------------
try:
f = open(data_dihedrals+'.template', 'r')
for line_orig in f:
ic = line_orig.find('#')
if ic != -1:
line = line_orig[:ic]
else:
line = line_orig.rstrip('\n')
tokens = line.strip().split()
if len(tokens) == 0:
pass
elif len(tokens) < 6:
raise InputError('Error('+g_program_name+'): The following line from\n'
' "\"'+data_dihedrals+'.template\" has bad format:\n\n'
+line_orig+'\n'
' This my probably an internal error. (Feel free to contact the developer.)\n'
+g_no_check_msg+'\n')
else:
defined_dihedrals.add(tokens[0])
f.close()
except:
pass #Defining dihedrals (stored in the data_dihedrals file) is optional
# ------------ defined_impropers ------------
try:
f = open(data_impropers+'.template', 'r')
for line_orig in f:
ic = line_orig.find('#')
if ic != -1:
line = line_orig[:ic]
else:
line = line_orig.rstrip('\n')
tokens = line.strip().split()
if len(tokens) == 0:
pass
elif len(tokens) < 6:
raise InputError('Error('+g_program_name+'): The following line from\n'
' "\"'+data_impropers+'.template\" has bad format:\n\n'
+line_orig+'\n'
' This my probably an internal error. (Feel free to contact the developer.)\n'
+g_no_check_msg+'\n')
else:
defined_impropers.add(tokens[0])
f.close()
except:
pass #Defining impropers (stored in the data_impropers file) is optional
# ---- Check ttree_assignments to make sure variables are defined ----
try:
f = open(ttree_assignments_fname, 'r')
except:
raise InputError('Error('+g_program_name+'): Unable to open file\n'+
'\"'+ttree_assignments_fname+'\"\n'
' for reading. (Do your files lack a \"'+data_atoms+'\" section?)\n'
+g_no_check_msg+'\n')
for line_orig in f:
ic = line_orig.find('#')
if ic != -1:
line = line_orig[:ic]
usage_location_str = 'near ' + line_orig[ic+1:]
else:
line = line_orig.rstrip('\n')
usage_location_str = ''
tokens = line.strip().split()
if len(tokens) == 0:
pass
if len(tokens) > 0:
# This file contains a list of variables of the form:
#
# @/atom:MoleculeType1:C 1
# @/atom:MoleculeType1:H 2
# @/atom:MoleculeType2:N 3
# $/atom:molecule1:N1 1
# $/atom:molecule1:C1 2
# :
# $/atom:molecule1141:CH 13578
# $/atom:molecule1142:N3 13579
# :
# We only care about instance variables (which use the '$' prefix)
# Lines corresponding to static variables (which use the '@' prefix)
# are ignored during this pass.
i_prefix = tokens[0].find('$')
if i_prefix != -1:
descr_str = tokens[0][i_prefix+1:]
cat_name = ExtractCatName(descr_str)
if ((cat_name == 'atom') and
(tokens[0] not in defined_atoms)):
raise InputError('Error('+g_program_name+'): '+usage_location_str+'\n'+
' Reference to undefined $atom:\n\n'
' '+tokens[0]+' (<--full name)\n\n'+
' (If that atom belongs to a molecule (or other subunit), make sure that\n'+
' you specified the correct path which leads to it (using / and ..))\n\n'+
g_no_check_msg)
elif ((cat_name == 'bond') and
(tokens[0] not in defined_bonds)):
raise InputError('Error('+g_program_name+'): '+usage_location_str+'\n'+
' Reference to undefined $bond:\n\n'
' '+tokens[0]+' (<--full name)\n\n'+
' (If that bond belongs to a molecule (or other subunit), make sure that\n'+
' you specified the correct path which leads to it (using / and ..))\n\n'+
g_no_check_msg)
elif ((cat_name == 'angle') and
(tokens[0] not in defined_angles)):
raise InputError('Error('+g_program_name+'): '+usage_location_str+'\n'+
' Reference to undefined $angle:\n\n'+
' '+tokens[0]+' (<--full name)\n\n'+
' (If that angle belongs to a molecule (or other subunit), make sure that\n'+
' you specified the correct path which leads to it (using / and ..))\n\n'+
g_no_check_msg)
elif ((cat_name == 'dihedral') and
(tokens[0] not in defined_dihedrals)):
raise InputError('Error('+g_program_name+'): '+usage_location_str+'\n\n'+
' Reference to undefined $dihedral:\n\n'
' '+tokens[0]+' (<--full name)\n\n'+
' (If that dihedral belongs to a molecule (or other subunit), make sure that\n'+
' you specified the correct path which leads to it (using / and ..))\n\n'+
g_no_check_msg)
elif ((cat_name == 'improper') and
(tokens[0] not in defined_impropers)):
raise InputError('Error('+g_program_name+'): '+usage_location_str+'\n'+
' Reference to undefined $improper:\n\n'
' '+tokens[0]+' (<--full name)\n\n'+
' (If that improper belongs to a molecule (or other subunit), make sure that\n'+
' you specified the correct path which leads to it (using / and ..))\n\n'+
g_no_check_msg)
# I used to generate an error when a users defines a $mol
# variable but does not associate any atoms with it (or if the
# user systematically deletes all the atoms in that molecule),
# but I stopped this practice.
# I don't think there is any real need to complain if some
# molecule id numbers are undefined. LAMMPS does not care.
#
#elif ((cat_name == 'mol') and
# (tokens[0] not in defined_mols)):
# raise InputError('Error('+g_program_name+'): '+usage_location_str+'\n'+
# ' Reference to undefined $mol (molecule-ID) variable:\n\n'
# ' '+tokens[0]+' (<--full name)\n\n'+
# ' (If that molecule is part of a larger molecule, then make sure that\n'+
# ' you specified the correct path which leads to it (using / and ..))\n\n'+
# g_no_check_msg)
f.close()
sys.stderr.write(g_program_name+': -- No errors detected. --\n')
exit(0)
except (ValueError, InputError) as err:
sys.stderr.write('\n'+str(err)+'\n')
sys.exit(1)

235
tools/moltemplate/src/lttree_styles.py Normal file
View File

@ -0,0 +1,235 @@
#!/usr/bin/env python
# -*- coding: utf-8 -*-
# Author: Andrew Jewett (jewett.aij at g mail)
# http://www.chem.ucsb.edu/~sheagroup
# License: 3-clause BSD License (See LICENSE.TXT)
# Copyright (c) 2012, Regents of the University of California
# All rights reserved.
from ttree_lex import InputError
# Users of lttree typically generate the following files:
# The variable below refer to file names generated by
# write() and write_once() commands in a lttree-file.
# (I keep changing my mind what I want these names to be.)
data_prefix="Data "
data_prefix_no_space="Data"
data_atoms="Data Atoms"
data_masses="Data Masses"
data_velocities="Data Velocities"
data_bonds="Data Bonds"
data_angles="Data Angles"
data_dihedrals="Data Dihedrals"
data_impropers="Data Impropers"
data_bond_coeffs="Data Bond Coeffs"
data_angle_coeffs="Data Angle Coeffs"
data_dihedral_coeffs="Data Dihedral Coeffs"
data_improper_coeffs="Data Improper Coeffs"
data_pair_coeffs="Data Pair Coeffs"
# interactions-by-type (not id. This is not part of the LAMMPS standard.)
data_angles_by_type="Data Angles By Type"
data_dihedrals_by_type="Data Dihedrals By Type"
data_impropers_by_type="Data Impropers By Type"
# class2 data sections
data_bondbond_coeffs="Data BondBond Coeffs"
data_bondangle_coeffs="Data BondAngle Coeffs"
data_middlebondtorsion_coeffs="Data MiddleBondTorsion Coeffs"
data_endbondtorsion_coeffs="Data EndBondTorsion Coeffs"
data_angletorsion_coeffs="Data AngleTorsion Coeffs"
data_angleangletorsion_coeffs="Data AngleAngleTorsion Coeffs"
data_bondbond13_coeffs="Data BondBond13 Coeffs"
data_angleangle_coeffs="Data AngleAngle Coeffs"
# sections for non-point-like particles:
data_ellipsoids="Data Ellipsoids"
data_lines="Data Lines"
data_triangles="Data Triangles"
# periodic boundary conditions
data_boundary="Data Boundary"
# (for backward compatibility), an older version of this file was named:
data_pbc="Data PBC"
# The files are fragments of a LAMMPS data file (see "read_data").
# In addition, moltemplate may also generate the following files:
in_prefix="In "
in_prefix_no_space="In"
in_init="In Init"
in_settings="In Settings"
in_coords="In Coords"
# These files represent different sections of the LAMMPS input script.
# Atom Styles in LAMMPS as of 2011-7-29
g_style_map = {'angle': ['atom-ID','molecule-ID','atom-type','x','y','z'],
'atomic': ['atom-ID','atom-type','x','y','z'],
'bond': ['atom-ID','molecule-ID','atom-type','x','y','z'],
'charge': ['atom-ID','atom-type','q','x','y','z'],
'dipole': ['atom-ID','atom-type','q','x','y','z','mux','muy','muz'],
'electron': ['atom-ID','atom-type','q','spin','eradius','x','y','z'],
'ellipsoid': ['atom-ID','atom-type','x','y','z','quatw','quati','quatj','quatk'],
'full': ['atom-ID','molecule-ID','atom-type','q','x','y','z'],
'line': ['atom-ID','molecule-ID','atom-type','lineflag','density','x','y','z'],
'meso': ['atom-ID','atom-type','rho','e','cv','x','y','z'],
'molecular': ['atom-ID','molecule-ID','atom-type','x','y','z'],
'peri': ['atom-ID','atom-type','volume','density','x','y','z'],
'sphere': ['atom-ID','atom-type','diameter','density','x','y','z'],
'tri': ['atom-ID','molecule-ID','atom-type','triangleflag','density','x','y','z'],
'wavepacket':['atom-ID','atom-type','charge','spin','eradius','etag','cs_re','cs_im','x','y','z'],
'hybrid': ['atom-ID','atom-type','x','y','z'],
# The following styles were removed from LAMMPS as of 2012-3
'colloid': ['atom-ID','atom-type','x','y','z'],
'granular': ['atom-ID','atom-type','diameter','density','x','y','z']}
def AtomStyle2ColNames(atom_style_string):
atom_style_string = atom_style_string.strip()
if len(atom_style_string) == 0:
raise InputError('Error: Invalid atom_style\n'
' (The atom_style command was followed by an empty string.)\n')
atom_style_args = atom_style_string.split()
atom_style = atom_style_args[0]
hybrid_args = atom_style_args[1:]
if (atom_style not in g_style_map):
if (len(atom_style_args) >= 2):
# If the atom_style_string includes at least 2 words, then we
# interpret this as a list of the individual column names
return atom_style_args
else:
raise InputError('Error: Unrecognized atom_style: \"'+atom_style+'\"\n')
if (atom_style != 'hybrid'):
return g_style_map[atom_style]
else:
column_names = ['atom-ID','atom-type','x','y','z']
if (len(hybrid_args)==0):
raise InputError('Error: atom_style hybrid must be followed by a sub_style.\n')
for sub_style in hybrid_args:
if (sub_style not in g_style_map):
raise InputError('Error: Unrecognized atom_style: \"'+sub_style+'\"\n')
for cname in g_style_map[sub_style]:
if cname not in column_names:
column_names.append(cname)
return column_names
def ColNames2AidAtypeMolid(column_names):
# Because of the diversity of ways that these
# numbers are referred to in the LAMMPS documentation,
# we have to be flexible and allow the user to refer
# to these quantities in a variety of ways.
# Hopefully this covers everything:
if 'atom-ID' in column_names:
i_atomid = column_names.index('atom-ID')
elif 'atom−ID' in column_names: # (− is the character used in the manual)
i_atomid = column_names.index('atom−ID')
elif 'atomID' in column_names:
i_atomid = column_names.index('atomID')
elif 'atomid' in column_names:
i_atomid = column_names.index('atomid')
elif 'id' in column_names:
i_atomid = column_names.index('id')
elif 'atom' in column_names:
i_atomid = column_names.index('atom')
elif '$atom' in column_names:
i_atomid = column_names.index('$atom')
else:
raise InputError('Error: List of column names lacks an \"atom-ID\"\n')
if 'atom-type' in column_names:
i_atomtype = column_names.index('atom-type')
elif 'atom−type' in column_names: # (− hyphen character used in manual)
i_atomtype = column_names.index('atom−type')
elif 'atomtype' in column_names:
i_atomtype = column_names.index('atomtype')
elif 'type' in column_names:
i_atomtype = column_names.index('type')
elif '@atom' in column_names:
i_atomtype = column_names.index('@atom')
else:
raise InputError('Error: List of column names lacks an \"atom-type\"\n')
i_molid = None
if 'molecule-ID' in column_names:
i_molid = column_names.index('molecule-ID')
elif 'molecule−ID' in column_names: # (− hyphen character used in manual)
i_molid = column_names.index('molecule−ID')
elif 'moleculeID' in column_names:
i_molid = column_names.index('moleculeID')
elif 'moleculeid' in column_names:
i_molid = column_names.index('moleculeid')
elif 'molecule' in column_names:
i_molid = column_names.index('molecule')
elif 'molID' in column_names:
i_molid = column_names.index('molID')
elif 'molid' in column_names:
i_molid = column_names.index('molid')
elif 'mol' in column_names:
i_molid = column_names.index('mol')
else:
pass # some atom_types do not have a valid molecule-ID
return i_atomid, i_atomtype, i_molid
def ColNames2Coords(column_names):
""" Which of the columns correspond to coordinates
which must be transformed using rigid-body
(affine: rotation + translation) transformations?
This function outputs a list of lists of triplets of integers.
"""
i_x = None
i_y = None
i_z = None
if 'x' in column_names:
i_x = column_names.index('x')
if 'y' in column_names:
i_y = column_names.index('y')
if 'z' in column_names:
i_z = column_names.index('z')
if (((i_x != None) != (i_y != None)) or
((i_y != None) != (i_z != None)) or
((i_z != None) != (i_x != None))):
raise InputError('Error: custom atom_style list must define x, y, and z.\n')
return [[i_x, i_y, i_z]]
def ColNames2Vects(column_names):
""" Which of the columns correspond to coordinates
which must be transformed using rotations?
Some coordinates like dipole moments and
ellipsoid orientations should only be rotated
(not translated).
This function outputs a list of lists of triplets of integers.
"""
vects = []
i_mux = None
i_muy = None
i_muz = None
if 'mux' in column_names:
i_mux = column_names.index('mux')
if 'muy' in column_names:
i_muy = column_names.index('muy')
if 'muz' in column_names:
i_muz = column_names.index('muz')
if (((i_mux != None) != (i_muy != None)) or
((i_muy != None) != (i_muz != None)) or
((i_muz != None) != (i_mux != None))):
raise InputError('Error: custom atom_style list must define mux, muy, and muz or none.\n')
if i_mux != None:
vects.append([i_mux, i_muy, i_muz])
return vects

1589
tools/moltemplate/src/moltemplate.sh Executable file
View File

File diff suppressed because it is too large Load Diff

41
tools/moltemplate/src/nbody_Angles.py Normal file
View File

@ -0,0 +1,41 @@
from nbody_graph_search import Ugraph
# To find 3-body "angle" interactions, we would use this subgraph:
#
#
# *---*---* => 1st bond connects atoms 0 and 1
# 0 1 2 2nd bond connects atoms 1 and 2
#
bond_pattern = Ugraph([(0,1), (1,2)])
# (Ugraph atom indices begin at 0, not 1)
# The next function eliminates the redundancy between 0-1-2 and 2-1-0:
def canonical_order(match):
"""
When searching for atoms with matching bond patterns GraphMatcher
often returns redundant results. We must define a "canonical_order"
function which sorts the atoms and bonds in a way which is consistent
with the type of N-body interaction being considered.
The atoms (and bonds) in a candidate match are rearranged by the
canonical_order(). Then the re-ordered list of atom and bond ids is
tested against the list of atom/bond ids in the matches-found-so-far,
before it is added to the list of interactions found so far.
(For example, it does not make sense to define a separate 3-body angle
interaction between atoms 1, 2, 3 AND 3, 2, 1. This is the same triplet
of atoms, and the angle between them is the same.)
"""
# match[0][0:2] contains the ID numbers for the 3 atoms in the match
atom0 = match[0][0]
atom1 = match[0][1]
atom2 = match[0][2]
# match[1][0:1] contains the ID numbers for the 2 bonds
bond0 = match[1][0]
bond1 = match[1][1]
if atom0 < atom2:
#return ((atom0, atom1, atom2), (bond0, bond1)) same thing as:
return match
else:
return ((atom2, atom1, atom0), (bond1, bond0))

31
tools/moltemplate/src/nbody_Bonds.py Normal file
View File

@ -0,0 +1,31 @@
from nbody_graph_search import Ugraph
# To find 2-body "bond" interactions, we would use this subgraph:
#
#
# *---* => one bond connects atoms 0 and 1
# 0 1
#
bond_pattern = Ugraph([(0,1)])
# (Ugraph atom indices begin at 0, not 1)
# The next function eliminates the redundancy between 0-1 and 1-0:
def canonical_order(match):
"""
It does not make sense to define a separate bond between atoms 1 and 2,
and between atoms 2 and 1. This function will swap the atoms in the bond
if the first atom > second atom.
"""
# match[0][0:2] contains the ID numbers for the 2 atoms in the match
atom0 = match[0][0]
atom1 = match[0][1]
# match[1][0:1] contains the ID numbers for the 1 bond
bond0 = match[1][0]
if atom0 < atom1:
#return ((atom0, atom1), (bond0)) same thing as:
return match
else:
return ((atom1, atom0), (bond0))

44
tools/moltemplate/src/nbody_Dihedrals.py Normal file
View File

@ -0,0 +1,44 @@
from nbody_graph_search import Ugraph
# To find 4-body "dihedral" interactions, we would use this subgraph:
#
# 1st bond connects atoms 0 and 1
# *---*---*---* => 2nd bond connects atoms 1 and 2
# 0 1 2 3 3rd bond connects atoms 2 and 3
#
bond_pattern = Ugraph([(0,1), (1,2), (2,3)])
# (Ugraph atom indices begin at 0, not 1)
def canonical_order(match):
"""
When searching for atoms with matching bond patterns GraphMatcher
often returns redundant results. We must define a "canonical_order"
function which sorts the atoms and bonds in a way which is consistent
with the type of N-body interaction being considered.
The atoms (and bonds) in a candidate match are rearranged by the
canonical_order(). Then the re-ordered list of atom and bond ids is
tested against the list of atom/bond ids in the matches-found-so-far,
before it is added to the list of interactions found so far.
(For example, it does not make sense to define a separate 4-body dihedral-
angle interaction between atoms 1, 2, 3, 4 AND 4, 3, 2, 1. The dihedral-
angle is not altered when the order of atoms is reversed, so we arbitrarily
choose whichever order causes the first atom to have the lower atomID.)
"""
# match[0][0:3] contains the ID numbers of the 4 atoms in the match
atom0 = match[0][0]
atom1 = match[0][1]
atom2 = match[0][2]
atom3 = match[0][3]
# match[1][0:2] contains the ID numbers of the the 3 bonds
bond0 = match[1][0]
bond1 = match[1][1]
bond2 = match[1][2]
if atom0 < atom3:
#return ((atom0, atom1, atom2, atom3), (bond0, bond1, bond2)) same as:
return match
else:
return ((atom3, atom2, atom1, atom0), (bond2, bond1, bond0))

41
tools/moltemplate/src/nbody_Impropers.py Normal file
View File

@ -0,0 +1,41 @@
from nbody_graph_search import Ugraph
# To find 4-body "improper" interactions, we would use this subgraph:
# 3
# * 1st bond connects atoms 0 and 1
# | => 2nd bond connects atoms 0 and 2
# _.*._ 3rd bond connects atoms 0 and 3
# *' 0 `*
# 1 2
#
bond_pattern = Ugraph([(0,1), (0,2), (0,3)])
# (Ugraph atom indices begin at 0, not 1)
def canonical_order(match):
"""
When searching for atoms with matching bond patterns GraphMatcher
often returns redundant results. We must define a "canonical_order"
function which sorts the atoms and bonds in a way which is consistent
with the type of N-body interaction being considered.
The atoms (and bonds) in a candidate match are rearranged by the
canonical_order(). Then the re-ordered list of atom and bond ids is
tested against the list of atom/bond ids in the matches-found-so-far,
before it is added to the list of interactions found so far.
(For example, it does not make sense to define a separate 4-body improper-
angle interaction between atoms 1, 2, 3, 4 AND 1, 3, 2, 4. The improper-
angle is defined as the angle between planes formed by atoms 1,2,3 & 2,3,4.
This is not effected by swapping the middle pair of atoms so we arbitrarily
sort them so that the second atom has a lower atomID than the third atom.)
"""
atom0 = match[0][0]
atom1 = match[0][1]
atom2 = match[0][2]
atom3 = match[0][3]
if atom1 <= atom2:
#return ((atom0,atom1,atom2,atom3), match[1]) same thing as:
return match
else:
return ((atom0,atom2,atom1,atom3), match[1])

598
tools/moltemplate/src/nbody_by_type.py Executable file
View File

@ -0,0 +1,598 @@
#!/usr/bin/env python
# Author: Andrew Jewett (jewett.aij at g mail)
# http://www.chem.ucsb.edu/~sheagroup
# License: 3-clause BSD License (See LICENSE.TXT)
# Copyright (c) 2011, Regents of the University of California
# All rights reserved.
man_page_text = """
nbody_by_type.py reads a LAMMPS data file (or an excerpt of a LAMMPS)
data file containing bonded many-body interactions by atom type
(and bond type), and generates a list of additional interactions
in LAMMPS format consistent with those type (to the standard out).
Typical Usage:
nbody_by_type.py X < old.data > new.data
--or--
nbody_by_type.py X \\
-atoms atoms.data \\
-bonds bonds.data \\
-nbody X.data \\
-nbodybytype X_by_type.data \\
> new_X.data
In both cases "X" denotes the interaction type, which
is either "Angles", "Dihedrals", or "Impropers".
(Support for other interaction types can be added by the user. See below.)
-------- Example 1 -------
nbody_by_type.py X < old.data > new.data
In this example, nbody_by_type.py reads a LAMMPS data file
"orig.data", and extracts the relevant section ("Angles",
"Dihedrals", or "Impropers"). It also looks a section named "X By Type",
(eg. "Angles By type", "Impropers By type", "Impropers By type")
which contains a list of criteria for automatically defining additional
interactions of that type. For example, this file might contain:
Angle By Type
7 1 2 1 * *
8 2 2 * * *
9 3 4 3 * *
The first column is an interaction type ID.
The next 3 columns are atom type identifiers.
The final 2 columns are bond type identifiers.
The * is a wildcard symbol indicating there is no preference for bond types
in this example. (Optionally, regular expressions can also be used to
define a type match, by enclosing the atom or bond type in / slashes.)
The first line tells us to that there should be a 3-body "Angle"
interaction of type "7" whenever an atom of type 1 is bonded to an atom
of type "2", which is bonded to another atom of type "1" again.
The second line tells us that an angle is defined whenever three atoms
are bonded together and the first two are of type "2".
(Redundant angle interactions are filtered.)
New interactions are created for every group of bonded
atoms which match these criteria if they are bonded together
in the relevant way for that interaction type (as determined by
nbody_X.py), and printed to the standard output. For example,
suppose you are automatically generating 3-body "Angle" interactions using:
nbody_by_type Angles < old.data > new.data
The file "new.data" will be identical to "old.data", however the
"Angles By Type" section will be deleted, and the following lines of
text will be added to the "Angles" section:
394 7 5983 5894 5895
395 7 5984 5895 5896
396 7 5985 5896 5897
: : : : :
847 9 14827 14848 14849
The numbers in the first column are counters which assign a ID to
every interaction of that type, and start where the original "Angles"
data left off (New angle ID numbers do not overlap with old ID numbers).
The text in the second column ("7", "9", ...) matches the text from the
first column of the "Angle By Type" section of the input file.
-------- Example 2 -------
nbody_by_type.py X \\
-atoms atoms.data \\
-bonds bonds.data \\
-nbody X.data \\
-nbodybytype X_by_type.data \\
-prefix "SOMESTRING" -suffix "ANOTHERSTRING" \\
> new_X.data
In particular, for Angle interactions:
nbody_by_type.py Angles \\
-atoms atoms.data \\
-bonds bonds.data \\
-nbody angles.data \\
-nbodybytype angles_by_type.data \\
> new_Angles.data
When run this way, nbody_by_type.py behaves exactly the same way
as in Example 1, however only the lines of text corresponding to
the new generated interactions are printed, (not the entire data file).
Also note, that when run this way, nbody_by_type.py does not read the
LAMMPS data from the standard input. Instead, it reads each section of
the data file from a different file indicated by the arguments following
the "-atoms", "-bonds", "-nbody", and "-nbodybytype" flags.
"Angles" is a 3-body interaction style. So when run this way,
nbody_by_type.py will create a 5 (=3+2) column file (new_Angles.data).
Note: the atom, bond and other IDs/types in need not be integers.
Note: This program must be distributed with several python modules, including:
nbody_Angles.py, nbody_Dihedrals.py, and nbody_Impropers.py. These
contain bond definitions for angular, dihedral, and improper interactions.
(In case any new interaction types are ever added to LAMMPS, with only
a few lines of python it is easy to edit to define new bonded
interaction types by supplying new "nbody_X.py" python module.
Refer to the modules listed above for examples.)
Note: Optional "-prefix" and "-suffix" arguments can be included to decorate
the interaction IDs (first column). For example, -prefix "auto_" and
-suffix "_angle", causes "new_Angles.data" to contain the following text:
auto_394_angle 7 5983 5894 5895
auto_395_angle 7 5984 5895 5896
auto_396_angle 7 5985 5896 5897
: : : : :
auto_847_angle 9 14827 14848 14849
"""
import sys
from extract_lammps_data import *
from nbody_by_type_lib import GenInteractions_str
from ttree_lex import *
from lttree_styles import AtomStyle2ColNames, ColNames2AidAtypeMolid
if sys.version < '2.6':
raise InputError('Error: Using python '+sys.version+'\n'
' Alas, you must upgrade to a newever version of python (2.7 or later).')
elif sys.version < '2.7':
sys.stderr.write('--------------------------------------------------------\n'
'----------------- WARNING: OLD PYTHON VERSION ----------\n'
' This program is untested on your python version ('+sys.version+').\n'
' PLEASE LET ME KNOW IF THIS PROGRAM CRASHES (and upgrade python).\n'
' -Andrew 2012-4-16\n'
'--------------------------------------------------------\n'
'--------------------------------------------------------\n')
from ordereddict import OrderedDict
else:
from collections import OrderedDict
def GenInteractions_lines(lines_atoms,
lines_bonds,
lines_nbody,
lines_nbodybytype,
atom_style,
g_bond_pattern,
canonical_order, #function to sort atoms and bonds
prefix='',
suffix='',
report_progress = False):
column_names = AtomStyle2ColNames(atom_style)
i_atomid, i_atomtype, i_molid = ColNames2AidAtypeMolid(column_names)
atomids_str = []
atomtypes_str = []
for iv in range(0, len(lines_atoms)):
line = lines_atoms[iv].strip()
if '#' in line:
icomment = line.find('#')
line = (line[:icomment]).strip()
if len(line) > 0:
tokens = SplitQuotedString(line)
if ((len(tokens) <= i_atomid) or (len(tokens) <= i_atomtype)):
raise(InputError('Error not enough columns on line '+str(iv+1)+' of \"Atoms\" section.'))
tokens = SplitQuotedString(line)
atomids_str.append(EscCharStrToChar(tokens[i_atomid]))
atomtypes_str.append(EscCharStrToChar(tokens[i_atomtype]))
bondids_str = []
bondtypes_str = []
bond_pairs = []
for ie in range(0, len(lines_bonds)):
line = lines_bonds[ie].strip()
if '#' in line:
icomment = line.find('#')
line = (line[:icomment]).strip()
if len(line) > 0:
tokens = SplitQuotedString(line)
if len(tokens) < 4:
raise(InputError('Error not enough columns on line '+str(ie+1)+' of \"Bonds\" section.'))
bondids_str.append(EscCharStrToChar(tokens[0]))
bondtypes_str.append(EscCharStrToChar(tokens[1]))
bond_pairs.append( (EscCharStrToChar(tokens[2]),
EscCharStrToChar(tokens[3])) )
typepattern_to_coefftypes = []
for i in range(0, len(lines_nbodybytype)):
line = lines_nbodybytype[i].strip()
if '#' in line:
icomment = line.find('#')
line = (line[:icomment]).strip()
if len(line) > 0:
tokens = SplitQuotedString(line)
if ((len(tokens) != 1 + g_bond_pattern.GetNumVerts()) and
(len(tokens) != 1 + g_bond_pattern.GetNumVerts()
+ g_bond_pattern.GetNumEdges())):
raise(InputError('Error: Wrong number of columns in \"By Type\" section of data file.\n'
'Offending line:\n'+
'\"'+line+'\"\n'
'Expected either '+
str(1 + g_bond_pattern.GetNumVerts()) + ' or ' +
str(1 + g_bond_pattern.GetNumVerts() +
g_bond_pattern.GetNumEdges())
+ ' colunms.'))
coefftype = EscCharStrToChar(tokens[0])
typepattern = []
for typestr in tokens[1:]:
if ((len(typestr) >= 2) and
(typestr[0] == '/') and (typestr[-1] == '/')):
regex_str = typestr[1:-1]
typepattern.append( re.compile(regex_str) )
else:
typepattern.append(EscCharStrToChar(typestr))
# If the user neglected to specify the bond types, assume '*'
if len(tokens) == 1 + g_bond_pattern.GetNumVerts():
typepattern += ['*'] * g_bond_pattern.GetNumEdges()
typepattern_to_coefftypes.append([typepattern, coefftype])
coefftype_to_atomids_str = GenInteractions_str(bond_pairs,
g_bond_pattern,
typepattern_to_coefftypes,
canonical_order,
atomids_str,
atomtypes_str,
bondids_str,
bondtypes_str,
report_progress)
lines_nbody_new = []
for coefftype, atomids_list in coefftype_to_atomids_str.items():
for atomids_found in atomids_list:
n = len(lines_nbody) + len(lines_nbody_new) + 1
line = prefix+str(n)+suffix+' '+ \
coefftype+' '+(' '.join(atomids_found))+'\n'
lines_nbody_new.append(line)
return lines_nbody_new
def GenInteractions_files(lines_data,
src_bond_pattern,
fname_atoms,
fname_bonds,
fname_nbody,
fname_nbodybytype,
section_name,
section_name_bytype,
atom_style,
prefix='',
suffix='',
report_progress = False):
if fname_atoms == None:
lines_atoms = [line for line in ExtractDataSection(lines_data, 'Atoms')]
else:
try:
f = open(fname_atoms, 'r')
except:
sys.stderr.write('Error: Unable to open file \"'+fname_atoms+'\" for reading.\n')
sys.exit(-1)
lines_atoms = [line for line in f.readlines()
if ((len(line.strip())>0) and (line.strip()[0] != '#'))]
f.close()
if fname_bonds == None:
lines_bonds = [line for line in ExtractDataSection(lines_data, 'Bonds')]
else:
try:
f = open(fname_bonds, 'r')
except IOError:
sys.stderr.write('Error: Unable to open file \"'+fname_bonds+'\" for reading.\n')
sys.exit(-1)
lines_bonds = [line for line in f.readlines()
if ((len(line.strip())>0) and (line.strip()[0] != '#'))]
f.close()
if fname_nbody == None:
lines_nbody = [line for line in ExtractDataSection(lines_data, section_name)]
else:
try:
f = open(fname_nbody, 'r')
lines_nbody = [line for line in f.readlines()
if ((len(line.strip())>0) and (line.strip()[0] != '#'))]
f.close()
except IOError:
#sys.stderr.write(' (omitting optional file \"'+fname_nbody+'\")\n')
lines_nbody = []
if fname_nbodybytype == None:
lines_nbodybytype=[line for
line in ExtractDataSection(lines_data,
section_name_bytype)]
else:
try:
f = open(fname_nbodybytype, 'r')
except:
sys.stderr.write('Error: Unable to open file \"'+fname_nbodybytype+'\" for reading.\n')
sys.exit(-1)
lines_nbodybytype = [line for line in f.readlines()
if((len(line.strip())>0)and(line.strip()[0]!='#'))]
f.close()
g = __import__(src_bond_pattern) #defines g.bond_pattern, g.canonical_order
return GenInteractions_lines(lines_atoms,
lines_bonds,
lines_nbody,
lines_nbodybytype,
atom_style,
g.bond_pattern,
g.canonical_order,
prefix,
suffix,
report_progress)
if __name__ == "__main__":
g_program_name = __file__.split('/')[-1] # = 'nbody_by_type.py'
g_date_str = '2013-4-16'
g_version_str = '0.11'
####### Main Code Below: #######
sys.stderr.write(g_program_name+' v'+g_version_str+' '+g_date_str+' ')
if sys.version < '3':
sys.stderr.write(' (python version < 3)\n')
else:
sys.stderr.write('\n')
try:
fname_atoms = None
fname_bonds = None
fname_nbody = None
fname_nbodybytype = None
atom_style = 'full'
prefix=''
suffix=''
argv = [arg for arg in sys.argv]
if len(argv) == 1:
raise InputError('Error: Missing argument required.\n'
' The \"'+g_program_name+'\" program requires an argument containing the\n'
' name of a section from a LAMMPS data file storing bonded interactions.\n'
' (For example: "Angles", "Dihedrals", or "Impropers".)\n'
#' Note: The first letter of each section is usually capitalized.)\n'
'\n'
'--------------- general documentation -------------\n'
'\n' + man_page_text + '\n')
section_name = '' # (This will be replaced later.)
# Loop over the remaining arguments not processed yet.
# These arguments are specific to the lttree.py program
# and are not understood by ttree.py:
i = 1
while i < len(argv):
#sys.stderr.write('argv['+str(i)+'] = \"'+argv[i]+'\"\n')
if ((argv[i].lower() == '-?') or
(argv[i].lower() == '--?') or
(argv[i].lower() == '-help') or
(argv[i].lower() == '-help')):
if i+1 >= len(argv):
sys.stdout.write(man_page_text+'\n')
sys.exit(0)
elif argv[i].lower() == '-atoms':
if i+1 >= len(argv):
raise InputError('Error: '+argv[i]+' flag should be followed by a file name containing lines of\n'
' text which might appear in the "Atoms" section of a LAMMPS data file.\n')
fname_atoms = argv[i+1]
del(argv[i:i+2])
elif argv[i].lower() == '-bonds':
if i+1 >= len(argv):
raise InputError('Error: '+argv[i]+' flag should be followed by a file name containing lines of\n'
' text which might appear in the "Bonds" section of a LAMMPS data file.\n')
fname_bonds = argv[i+1]
del(argv[i:i+2])
elif argv[i].lower() == '-nbody':
if i+1 >= len(argv):
raise InputError('Error: '+argv[i]+' flag should be followed by a file name\n')
#raise InputError('Error: '+argv[i]+' flag should be followed by a file name containing lines of\n'
# ' text which might appear in the "'+section_name+' section of a LAMMPS data file.\n')
fname_nbody = argv[i+1]
del(argv[i:i+2])
elif argv[i].lower() == '-nbodybytype':
if i+1 >= len(argv):
raise InputError('Error: '+argv[i]+' flag should be followed by a file name\n')
#raise InputError('Error: '+argv[i]+' flag should be followed by a file name containing\n'
# ' text which might appear in the "'+section_name+' By Type" section\n'
# ' of a LAMMPS data file.\n')
fname_nbodybytype = argv[i+1]
del(argv[i:i+2])
elif ((argv[i].lower() == '-atom-style') or
(argv[i].lower() == '-atom_style')):
if i+1 >= len(argv):
raise InputError('Error: '+argv[i]+' flag should be followed by a an atom_style name.\n'
' (Or single quoted string which includes a space-separated\n'
' list of column names.)\n')
atom_style = argv[i+1]
del(argv[i:i+2])
elif argv[i].lower() == '-prefix':
if i+1 >= len(argv):
raise InputError('Error: '+argv[i]+' flag should be followed by a prefix string\n'
' (a string you want to appear to the left of the integer\n'
' which counts the bonded interactions you have generated.)\n')
prefix = argv[i+1]
del(argv[i:i+2])
elif argv[i].lower() == '-suffix':
if i+1 >= len(argv):
raise InputError('Error: '+argv[i]+' flag should be followed by a suffix string\n'
' (a string you want to appear to the right of the integer\n'
' which counts the bonded interactions you have generated.)\n')
prefix = argv[i+1]
del(argv[i:i+2])
elif argv[i][0] == '-':
raise InputError('Error('+g_program_name+'):\n'
'Unrecogized command line argument \"'+argv[i]+'\"\n')
else:
i += 1
if len(argv) == 1:
raise InputError('Error: Missing argument required.\n'
' The \"'+g_program_name+'\" program requires an argument containing the\n'
' name of a section from a LAMMPS data file storing bonded interactions.\n'
' (For example: "Angles", "Dihedrals", or "Impropers".)\n')
#' Note: The first letter of each section is usually capitalized.)\n'
elif len(argv) == 2:
section_name = argv[1]
bond_pattern_module_name = 'nbody_'+section_name
del(argv[1:2])
else:
# if there are more than 2 remaining arguments,
problem_args = ['\"'+arg+'\"' for arg in argv[1:]]
raise InputError('Syntax Error('+g_program_name+'):\n\n'
' Problem with argument list.\n'
' The remaining arguments are:\n\n'
' '+(' '.join(problem_args))+'\n\n'
' (The actual problem may be earlier in the argument list.)\n')
# ------------ Done parsing argument list ----------
section_name_bytype = section_name + ' By Type'
if (fname_atoms or fname_bonds or fname_nbody or fname_nbodybytype):
output_full_DATA_file = False
lines_data = []
else:
output_full_DATA_file = True
lines_data = sys.stdin.readlines()
# Calculate the interactions and generate a list of lines of text
lines_new_interactions = \
GenInteractions_files(lines_data,
bond_pattern_module_name,
fname_atoms,
fname_bonds,
fname_nbody,
fname_nbodybytype,
section_name,
section_name_bytype,
atom_style,
prefix,
suffix,
report_progress=True)
# Print this text to the standard out.
# Question: Do we write out the entire DATA file,
# or just the portion that was generated by this program?
if not output_full_DATA_file:
# ...then only print out the interactions which were generated
# by this program, omitting any lines from the original data file:
# (This is the way I usually run this program.)
for line in lines_new_interactions:
sys.stdout.write(line)
else:
# ...then print out the entire data file, deleting the "By Type"
# section, and adding the generated lines of text to the corresponding
# If present, update the interaction counter at the beginning
# of the LAMMPS data file. (For example, if if 100 new "Angles"
# interactions were generated, replace "2 Angles" with "102 Angles")
#
for i in range(0, len(lines_data)):
line = lines_data[i].strip()
tokens = SplitQuotedString(line)
# updating the interaction counter
if ((len(tokens) == 2) and (tokens[1] == (section_name).lower())):
tokens[0] = str(int(tokens[0]) + len(lines_new_interactions))
lines_data[i] = ' '.join(tokens) + '\n'
# stop when you come to a section header
elif line in lammps_data_sections:
#"lammps_data_sections" is defined in "extract_lammps_data.py"
break
# locate the appropriate section of the data file
# (storing the type of interactions we just created)
i_nbody_a, i_nbody_b = \
FindDataSection(lines_data, section_name)
if i_nbody_a == -1:
if len(lines_new_interactions) > 0:
# If not found, create a new section at the end of the file,
# containing a section name followed by the list of lines
lines_data += ['\n', section_name+'\n', '\n'] + \
lines_new_interactions + ['\n']
else:
# Insert the new lines into the existing section
lines_data[i_nbody_b:i_nbody_b] = lines_new_interactions
# Figure out where the "By Type" section is located
# (so we skip over it)
i_bytype_a, i_bytype_b = \
FindDataSection(lines_data, section_name_bytype)
in_bytype_section = False
for i in range(0, len(lines_data)):
line = lines_data[i].strip()
# Omit all lines of text in the 'By Type' section (including the
# header and commments or blank lines which immediately follow it.)
if line == section_name_bytype:
in_bytype_section = True
elif i == i_bytype_b:
in_bytype_section = False
if not in_bytype_section:
sys.stdout.write(lines_data[i])
except (ValueError, InputError) as err:
sys.stderr.write('\n'+str(err)+'\n')
sys.exit(-1)

312
tools/moltemplate/src/nbody_by_type_lib.py Normal file
View File

@ -0,0 +1,312 @@
#!/usr/bin/env python
# Author: Andrew Jewett (jewett.aij at g mail)
# http://www.chem.ucsb.edu/~sheagroup
# License: 3-clause BSD License (See LICENSE.TXT)
# Copyright (c) 2012, Regents of the University of California
# All rights reserved.
import sys
from nbody_graph_search import *
#from collections import namedtuple
from collections import defaultdict, OrderedDict
from ttree_lex import MatchesPattern, MatchesAll, InputError
#import gc
def GenInteractions_int(G_system,
g_bond_pattern,
typepattern_to_coefftype,
canonical_order, #function to sort atoms and bonds
atomtypes_int2str,
bondtypes_int2str,
report_progress = False): #print messages to sys.stderr?
"""
GenInteractions() automatically determines a list of interactions
present in a system of bonded atoms (argument "G_system"),
which satisfy the bond topology present in "g_bond_pattern", and
satisfy the atom and bond type requirements in "typepattern_to_coefftype".
Whenever a set of atoms in "G_system" are bonded together in a way which
matches "g_bond_pattern", and when the atom and bond types is consistent
with one of the entries in "typepattern_to_coefftype", the corresponding
list of atoms from G_system is appended to the list of results.
These results (the list of lists of atoms participating in an interaction)
are organized according their corresponding "coefftype", a string
which identifies the type of interaction they obey as explained above.
results are returned as a dictionary using "coefftype" as the lookup key.
Arguments:
-- typepattern_to_coefftype is a list of 2-tuples --
The first element of the 2-tuple is the "typepattern".
It contains a string describing a list of atom types and bond types.
The typepattern is associated with a "coefftype",
which is the second element of the 2-tuple. This is a string
which identifies the type of interaction between the atoms.
Later on, this string can be used to lookup the force field
parameters for this interaction elsewhere.)
-- Arguments: G_system, g_bond_pattern, atomtypes_int2str, bondtypes_int2str --
G_system stores a list of atoms and bonds, and their attributes in
"Ugraph" format. In this format:
Atom ID numbers are represented by indices into the G_system.verts[] list.
Bond ID numbers are represented by indices into the G_system.edges[] list.
Atom types are represented as integers in the G_system.verts[i].attr list.
Bond types are represented as integers in the G_system.edges[i].attr list.
They are converted into strings using
atomtypes_int2str, and bondtypes_int2str.
g_bond_pattern is a graph which specifies the type of bonding between
the atoms required for a match. It is in Ugraph format (however the
atom and bond types are left blank.)
Atom and bond types are supplied by the user in string format. (These
strings typically encode integers, but could be any string in principle.)
The string-version of the ith atom type is stored in
atomtypes_int2str[ G_system.verts[i].attr ]
The string-version of the ith bond type is stored in
bondtypes_int2str[ G_system.edges[i].attr ]
-- The "canonical_order" argument: --
The search for atoms with a given bond pattern often yields
redundant matches. There is no difference for example
between the angle formed between three consecutively
bonded atoms (named, 1, 2, 3, for example), and the
angle between the same atoms in reverse order (3, 2, 1).
However both triplets of atoms will be returned by the subgraph-
matching algorithm when searching for ALL 3-body interactions.)
To eliminate this redundancy, the caller must supply a "canonical_order"
argument. This is a function which sorts the atoms and bonds in a way
which is consistent with the type of N-body interaction being considered.
The atoms (and bonds) in a candidate match are rearranged by the
canonical_order(). Then the re-ordered list of atom and bond ids is
tested against the list of atom/bond ids in the matches-found-so-far,
before it is added.
"""
if report_progress:
startatomid = 0
sys.stderr.write(' Searching for matching bond patterns:\n')
sys.stderr.write(' 0%')
# Figure out which atoms from "G_system" bond together in a way which
# matches the "g_bond_pattern" argument. Organize these matches by
# atom and bond types and store all of the non-redundant ones in
# the "interactions_by_type" variable.
gm = GraphMatcher(G_system, g_bond_pattern)
interactions_by_type = defaultdict(list)
for atombondids in gm.Matches():
# "atombondids" is a tuple.
# atombondids[0] has atomIDs from G_system corresponding to g_bond_pattern
# (These atomID numbers are indices into the G_system.verts[] list.)
# atombondids[1] has bondIDs from G_system corresponding to g_bond_pattern
# (These bondID numbers are indices into the G_system.edges[] list.)
# It's convenient to organize the list of interactions-between-
# atoms in a dictionary indexed by atomtypes and bondtypes.
# (Because many atoms and bonds typically share the same type,
# organizing the results this way makes it faster to check
# whether a given interaction matches a "typepattern" defined
# by the user. We only have to check once for the whole group.)
atombondtypes = \
(tuple([G_system.GetVert(Iv).attr for Iv in atombondids[0]]),
tuple([G_system.GetEdge(Ie).attr for Ie in atombondids[1]]))
interactions_by_type[atombondtypes].append(atombondids)
if report_progress:
# GraphMatcher.Matches() searches for matches in an order
# that selects a different atomid number from G_system,
# starting at 0, and continuing up to the number of atoms (-1)
# in the system (G_system.nv-1), and using this as the first
# atom in the match (ie match[0][0]). This number can be used
# to guess much progress has been made so far.
oldatomid = startatomid
startatomid = atombondids[0][0]
percent_complete = (100 * startatomid) // G_system.GetNumVerts()
# report less often as more progress made
if percent_complete <= 4:
old_pc = (100 * oldatomid) // G_system.GetNumVerts()
if percent_complete > old_pc:
sys.stderr.write(' '+str(percent_complete)+'%')
elif percent_complete <= 8:
pc_d2 = (100 * startatomid) // (2*G_system.GetNumVerts())
oldpc_d2 = (100 * oldatomid) // (2*G_system.GetNumVerts())
if pc_d2 > oldpc_d2:
sys.stderr.write(' '+str(percent_complete)+'%')
elif percent_complete <= 20:
pc_d4 = (100 * startatomid) // (4*G_system.GetNumVerts())
oldpc_d4 = (100 * oldatomid) // (4*G_system.GetNumVerts())
if pc_d4 > oldpc_d4:
sys.stderr.write(' '+str(percent_complete)+'%')
else:
pc_d10 = (100 * startatomid) // (10*G_system.GetNumVerts())
oldpc_d10 = (100 * oldatomid) // (10*G_system.GetNumVerts())
if pc_d10 > oldpc_d10:
sys.stderr.write(' '+str(percent_complete)+'%')
if report_progress:
sys.stderr.write(' 100%\n')
sys.stderr.write(' Looking up atom and bond types...')
# Now test each match to see if the types of atoms and bonds involved match
# any of the type-patterns in the "typepattern_to_coefftype" argument.
# If so, store them in the out_topo list
#coefftype_to_atomids = defaultdict(list)
#abids_to_coefftypes = defaultdict(list)
coefftype_to_atomids = OrderedDict()
abids_to_coefftypes = OrderedDict()
count = 0
for typepattern, coefftype in typepattern_to_coefftype:
if report_progress:
sys.stderr.write(' Checking (atom-types,bond-types) against \n '+str(typepattern)+'-->'+coefftype+'\n')
for atombondtypes, abidslist in interactions_by_type.items():
# express atom & bond types in a tuple of the original string format
types_atoms = [atomtypes_int2str[Iv] for Iv in atombondtypes[0]]
types_bonds = [bondtypes_int2str[Ie] for Ie in atombondtypes[1]]
type_strings = types_atoms + types_bonds
# use string comparisons to check for a match with typepattern
if MatchesAll(type_strings, typepattern): #<-see "ttree_lex.py"
for abids in abidslist:
# Re-order the atoms (and bonds) in a "canonical" way.
# Only add new interactions to the list after re-ordering
# them and checking that they have not been added earlier.
# (...well not when using the same coefftype at least.
# This prevents the same triplet of atoms from
# being used to calculate the bond-angle twice:
# once for 1-2-3 and 3-2-1, for example.)
abids = canonical_order(abids)
redundant = False
if abids in abids_to_coefftypes:
coefftypes = abids_to_coefftypes[abids]
if coefftype in coefftypes:
redundant = True
if not redundant:
# (It's too bad python does not
# have an Ordered defaultdict)
if coefftype in coefftype_to_atomids:
coefftype_to_atomids[coefftype].append(abids[0])
else:
coefftype_to_atomids[coefftype]=[abids[0]]
if abids in abids_to_coefftypes:
abids_to_coefftypes[abids].append(coefftype)
else:
abids_to_coefftypes[abids] = [coefftype]
count += 1
if report_progress:
sys.stderr.write(' done\n (found '+
str(count)+' non-redundant matches)\n')
return coefftype_to_atomids
def GenInteractions_str(bond_pairs,
g_bond_pattern,
typepattern_to_coefftype,
canonical_order, #function to sort atoms and bonds
atomids_str,
atomtypes_str,
bondids_str,
bondtypes_str,
report_progress = False): #print messages to sys.stderr?
assert(len(atomids_str) == len(atomtypes_str))
assert(len(bondids_str) == len(bondtypes_str))
# The atomids and atomtypes and bondtypes are strings.
# First we assign a unique integer id to each string.
atomids_str2int = {}
atomtypes_str2int = {}
atomtypes_int2str = []
atomtype_int = 0
for i in range(0, len(atomids_str)):
if atomids_str[i] in atomids_str2int:
raise InputError('Error: multiple atoms have the same id ('+
str(atomids_str[i])+')')
atomids_str2int[atomids_str[i]] = i
#atomtypes_int = len(atomtypes_int)+1
if (not (atomtypes_str[i] in atomtypes_str2int)):
atomtypes_str2int[atomtypes_str[i]] = atomtype_int
atomtypes_int2str.append(atomtypes_str[i])
atomtype_int += 1
#atomtypes_int.append(atomtype_int)
bondids_str2int = {}
bondtypes_str2int = {}
bondtypes_int2str = []
bondtype_int = 0
for i in range(0, len(bondids_str)):
if bondids_str[i] in bondids_str2int:
raise InputError('Error: multiple bonds have the same id ('+
str(bondids_str[i])+')')
bondids_str2int[bondids_str[i]] = i
#bondtype_int = len(bondtypes_int)+1
if (not (bondtypes_str[i] in bondtypes_str2int)):
bondtypes_str2int[bondtypes_str[i]] = bondtype_int
bondtypes_int2str.append(bondtypes_str[i])
bondtype_int += 1
# Now convert "bond_pairs" into the UGraph format
G_system = Ugraph()
for iv in range(0, len(atomtypes_str)):
G_system.AddVertex(iv, atomtypes_str2int[atomtypes_str[iv]])
for ie in range(0, len(bond_pairs)):
atomid1_str = bond_pairs[ie][0]
atomid2_str = bond_pairs[ie][1]
if (atomid1_str not in atomids_str2int):
raise InputError('Error in Bonds Section:\n'
' '+atomid1_str+' is not defined in Atoms section\n')
if (atomid2_str not in atomids_str2int):
raise InputError('Error in Bonds Section:\n'
' '+atomid2_str+' is not defined in Atoms section\n')
G_system.AddEdge(atomids_str2int[atomid1_str],
atomids_str2int[atomid2_str],
bondtypes_str2int[bondtypes_str[ie]])
coefftype_to_atomids_int = GenInteractions_int(G_system,
g_bond_pattern,
typepattern_to_coefftype,
canonical_order,
atomtypes_int2str,
bondtypes_int2str,
report_progress)
coefftype_to_atomids_str = OrderedDict()
for coefftype, atomidss_int in coefftype_to_atomids_int.items():
sys.stderr.write(' processing coefftype='+str(coefftype)+'\n')
for atomids_int in atomidss_int:
if coefftype in coefftype_to_atomids_str:
coefftype_to_atomids_str[coefftype].append(
[atomids_str[iv] for iv in atomids_int])
else:
coefftype_to_atomids_str[coefftype] = \
[[atomids_str[iv] for iv in atomids_int]]
#gc.collect()
return coefftype_to_atomids_str

149
tools/moltemplate/src/nbody_fix_ttree_assignments.py Executable file
View File

@ -0,0 +1,149 @@
#!/usr/bin/env python
"""
nbody_fix_ttree_assignments.py
This is an ugly little script which was not intended to be run by end users.
Typical usage:
nbody_fix_ttree_assignments.py "angles" new_Angles.template \
< ttree_assignments.txt > ttree_assigmnents_new.txt
What it does:
In this example, this program extracts the first column from
"new_Angles.template", and appends to it the first column from
the lines in ttree_assignments.txt containing "@angle:".
Then it adds a second column which is just a sequence of integers
counting upwards.
Finally it inserts this 2-column text into the appropriate place in a
ttree_assignments.txt file, replacing the original @angle variables
with the new ones (followed by the renumbered original).
AWK/GREP equivalent
This program is roughly equivalent to the following lines of awk/grep:
awk 'BEGIN{i=-1} {if(substr($0,0,8)=="$/angle") i=NR; if (i==-1){print $0}}'\
< ttree_assignments.txt > ttree_assignments_new.txt
awk '{print $1}' < new_Angles.template > Angles_column1.txt
grep '$/angle:' ttree_assignments.txt | awk '{print $1}' >> Angles_column1.txt
awk '{print $1 " " NR}' < Angles_column1.txt >> ttree_assignments_new.txt
awk 'BEGIN{found=0;passed=0} {if(substr($0,0,8)=="$/angle") found=1;
else {if (found) {passed=1}}
if (passed) print $0}' \
< ttree_assignments.txt >> ttree_assignments_new.txt
I wrote this python script (instead of using awk) just to handle quoted stings
(and strings with other fancy characters and escape sequences).
"""
import sys
from ttree_lex import SplitQuotedString, EscCharStrToChar, SafelyEncodeString, InputError
g_program_name = __file__.split('/')[-1]
try:
if (len(sys.argv) != 3):
raise InputError('Error running \"'+g_program_name+'\"\n'
' Wrong number of arguments.\n'
' (This is likely a programmer error.\n'
' This script was not intended to be run by end users.)\n')
cat_name = sys.argv[1]
f = open(sys.argv[2])
lines_generated = f.readlines()
f.close()
# Selections are simply lists of 2-tuples (pairs)
#f = open('ttree_assignments.txt','r')
#lines_bindings = f.readlines()
#f.close()
lines_bindings = sys.stdin.readlines()
# Figure out which lines in the 'ttree_assignments.txt' file
# contain the variables of the type you are looking for.
# Make note of the relevant line numbers
i_preexisting_begin = -1
i_preexisting_end = -1
in_section = False
possible_cat_names = set(['$'+cat_name, '$/'+cat_name, '${'+cat_name, '${/'+cat_name])
preexisting_interaction_list = []
for i in range(0, len(lines_bindings)):
line = lines_bindings[i].strip()
tokens = SplitQuotedString(line) #strip comments, handle quotes
if len(tokens) == 2:
before_colon = tokens[0].split(':')[0]
if before_colon in possible_cat_names:
if i_preexisting_begin == -1:
i_preexisting_begin = i
in_section = True
else:
if in_section:
i_preexisting_end = i
in_section = False
if i_preexisting_end == -1:
i_preexisting_end = len(lines_bindings)
if i_preexisting_begin == -1:
for line in lines_bindings:
sys.stdout.write(line)
else:
# write out all the lines in the original file up until the point where
# the variables in the category we are looking for were encountered
for i in range(0, i_preexisting_begin):
sys.stdout.write(lines_bindings[i])
sys.stderr.write(' (adding new lines)\n')
# Now add some new lines (2-column format).
# As with any ttree_assignment.txt file:
# The first column has our generated variable names
# The second column has the counter assigned to that variable
new_counter = 1
for line_orig in lines_generated:
line = line_orig.strip()
if len(line) > 0:
tokens = SplitQuotedString(line) #strip comments, handle quotes
sys.stdout.write(tokens[0]+' '+str(new_counter)+'\n')
new_counter += 1
sys.stderr.write(' (adding pre-exisiting lines)\n')
if i_preexisting_begin != -1:
# Append the original pre-existing interactions of that type, but assign
# them to higher numbers. (Hopefully this helps to make sure that these
# assignments will override any of the automatic/generated assignments.)
# As with any ttree_assignment.txt file:
# The first column has our generated variable names
# The second column has the counter assigned to that variable
#sys.stderr.write(' i_preexisting_begin='+
# str(i_preexisting_begin)+
# ' i_preexisting_end='+str(i_preexisting_end)+'\n')
for i in range(i_preexisting_begin, i_preexisting_end):
line = lines_bindings[i].strip()
tokens = SplitQuotedString(line) #strip comments, handle quotes
if len(tokens) == 2:
sys.stdout.write(tokens[0]+' '+str(new_counter)+'\n')
new_counter += 1
#sys.stderr.write(' (writing pre-exisiting lines)\n')
# write out all the lines in the original file after this point.
for i in range(i_preexisting_end, len(lines_bindings)):
sys.stdout.write(lines_bindings[i])
sys.exit(0)
except (ValueError, InputError) as err:
sys.stderr.write('\n'+str(err)+'\n')
sys.exit(-1)

974
tools/moltemplate/src/nbody_graph_search.py Normal file
View File

@ -0,0 +1,974 @@
# Author: Andrew Jewett (jewett.aij at g mail)
# http://www.chem.ucsb.edu/~sheagroup
# License: 3-clause BSD License (See LICENSE.TXT)
# Copyright (c) 2012, Regents of the University of California
# All rights reserved.
#__all__ = ['Ugraph', 'GraphMatcher', 'DFS', 'GenError', 'GraphError', 'Disconnected', 'NotUndirected']
import sys
import copy
from operator import itemgetter
class GenError(Exception):
"""
An exception class containing string for error reporting.
"""
def __init__(self, err_msg):
self.err_msg = err_msg
def __str__(self):
return self.err_msg
def __repr__(self):
return str(self)
class GraphError(GenError):
"""
An exception class containing a graph and a string for error reporting.
"""
def __init__(self, g, err_msg):
GenError.__init__(self, err_msg)
self.g = g
def __str__(self):
g_str = str(g)
# If the string representation of the graph is too
# large to fit in one screen, truncate it
g_str_lines = g_str.split('\n')
if (len(g_str_lines) > 12):
g_str_lines = g_str_lines[0:12] + [' ...(additional lines not shown)]']
g_str = '\n'.join(g_str_lines)
return 'Problem with graph:\n'+g_str+'\n'+self.err_msg
def __repr__(self):
return str(self)
class Disconnected(GraphError):
def __init__(self, g, err_msg):
GraphError.__init__(self, g, err_msg)
class NotUndirected(GraphError):
def __init__(self, g, err_msg):
GraphError.__init__(self, g, err_msg)
class Edge(object):
def __init__(self,
iv_start, # edge starts here (index into vertex list)
iv_stop, # edge ends here (index into vertex list)
attr=None): # edges have an optional type attribute
self.start = iv_start
self.stop = iv_stop
self.attr = attr
def __str__(self):
return '('+str(self.start)+','+str(self.stop)+')'
def __repr__(self):
return str(self)
class Vertex(object):
def __init__(self, attr=None):
self.attr = attr
class Dgraph(object):
"""
This class is a minimal implementation of a directed graph.
Vertices and edges are accessed by integer index only (beginning at 0).
Multiple edges connecting the same pair of vertices are allowed.
(One would use the AddEdge() member function to accomplish this.)
Both vertices and edges have an optional "attr" attribute.
"""
NULL = -1 # forbidden vertex id number (used several places)
def __init__(self, edgelist=None):
"""
The constructor accepts an optional neighborlist argument.
This is a simple list of neighbors for every vertex in the graph
and it completely defines the topology of the graph.
(Vertex and edge attributes can be specified later.)
Alternatley, you can leave the neighborlist argument blank,
and build the graph one vertex at a time later
using the "AddVertex()" and "AddEdge()" commands.
(AddEdge() commands must be issued strictly after
all vertices have been defined.)
"""
if edgelist == None:
self.verts = []
self.edges = []
self.nv = 0 #integer keeps track of # of vertices = len(self.verts)
self.ne = 0 #integer keeps track of # of edges = len(self.edges)
self.neighbors = [] # The adjacency list.
else:
# Parse the edge-list format:
iv_max = 0 # <-- what's the vertex with the maximum id number?
for i in range(0, len(edgelist)):
iv = edgelist[i][0]
jv = edgelist[i][1]
if ((iv < 0) or (jv < 0)):
raise(GenError('Error in Dgraph.__init__: Negative vertex number pair encountered: ('+str(iv)+','+str(jv)+')'))
if iv > iv_max:
iv_max = iv
if jv > iv_max:
iv_max = jv
self.nv = iv_max+1
self.verts = [Vertex() for iv in range(0, self.nv)]
self.edges = []
self.ne = 0
self.neighbors = [[] for iv in range(0, self.nv)]
for i in range(0, len(edgelist)):
iv = edgelist[i][0]
jv = edgelist[i][1]
self.neighbors[iv].append(self.ne)
self.edges.append(Edge(iv, jv))
self.ne += 1
assert(self.ne == len(self.edges))
self.SortNeighborLists()
def AddVertex(self, iv=-1, attr=None):
"""
Add a vertex to the graph.
(Edges connected to this vertex must be added later using "AddEdge()"
All vertices should be added before "AddEdge()" is ever invoked.)
Optional "attr" argument allows you to set the attribute of this vertex.
(for example, in a molecule this might correspond to the type of atom
in the molecule).
Optional "iv" argument allows you to specify the index of that vertex.
Vertices can be added in any order, but thei vertex id numbers
should eventually fill the range from 0 to self.nv-1.
"""
if iv == -1: # if iv unspecified, put the vertex at the end of the list
iv = self.nv
if iv < self.nv:
self.verts[iv].attr = attr
else:
# In case there is a gap between iv and nv, fill it with blanks
self.verts += ([Vertex()] * ((1 + iv) - self.nv))
self.neighbors += ([[]] * ((1 + iv) - self.nv))
self.verts[iv].attr = attr
self.nv = iv+1
assert(self.nv == len(self.verts))
assert(self.nv == len(self.neighbors))
def AddEdge(self, iv, jv, attr=None, remove_duplicates=False):
"""
Add an edge to graph connecting vertex iv to jv.
(both are integers from 0 to self.nv-1)
This function must not be called until all vertices have been added.
If the edge is already present (and remove_duplicates==True),
no new edge will be added.
"""
if remove_duplicates:
for je in self.neighbors[iv]:
if jv == self.edges[je].stop:
return # In that case, do nothing, the edge is already present
self.edges.append(Edge(iv, jv, attr))
self.neighbors[iv].append(self.ne)
self.ne += 1
assert(self.ne == len(self.edges))
def ReorderVerts(self, vpermutation, invert=False):
"""
This function allows the user to re-order (relabel) the vertices
in a graph, making the necessary changes to the
self.verts, self.edges, and self.neighbors lists.
By default (invert=False). The vpermutation is a list
from 1 to self.nv which is interpreted this way:
iv = vpermutation[iv_orig]
where "iv" and "iv_orig" are the vertex id numbers before
and after the mapping (which also corresponds to its
position in the self.verts and self.neighbors arrays).
"""
assert(len(self.verts) == self.nv)
assert(len(self.edges) == self.ne)
assert(len(vpermutation) == self.nv)
if (invert):
vperm = [-1 for iv in vpermutation]
for iv in range(0, self.nv):
vperm[ vpermutation[iv] ] = iv
else:
vperm = vpermutation
orig_verts = [vert for vert in self.verts]
for iv_old in range(0, self.nv):
iv = vperm[iv_old]
self.verts[iv] = orig_verts[iv_old]
for ie in range(0, self.ne):
self.edges[ie].start = vperm[self.edges[ie].start]
self.edges[ie].stop = vperm[self.edges[ie].stop]
orig_neighbors = [nlist for nlist in self.neighbors]
# self.neighbors is a 2-d array.
# We need to re-sort "self.neighbors" because the first index is
# a vertex id number, and these id numbers have been permuted.
# However, there's no need to sort the contents of each sub-array
# (self.neighbors[iv]), because these are edge id numbers (indices into
# the self.edges[] array). These edge index numbers are never altered.
# (However the entries stored in self.edges were modified earlier.)
for iv_old in range(0, self.nv):
iv = vperm[iv_old]
self.neighbors[iv] = orig_neighbors[iv_old]
# Optional:
self.SortNeighborLists()
def ReorderEdges(self, epermutation, invert=False):
"""
This function allows the user to re-order (relabel) the
edges in a graph, making the necessary changes to the
self.edges and self.neighbors lists.
By default (invert=False). The epermutation is a list
from 1 to self.ne which is interpreted this way:
ie = epermutation[ie_orig]
where "ie" and "ie_orig" are the edge id numbers before
and after the mapping (which also corresponds to that edge's
position in the self.edges array).
(Minor detail: Recall that in this code, Ugraphs
are implemented by placing two (directed) edges between each pair of
connected, adjacent vertices, which point back-and-forth between them.
Consequently the list of edges in self.edges is often typically
twice as large you might expect.)
"""
assert(len(self.verts) == self.nv)
assert(len(self.edges) == self.ne)
assert(len(epermutation) == self.ne)
if (invert):
eperm = [-1 for ie in epermutation]
for ie in range(0, self.ne):
eperm[ epermutation[ie] ] = ie
else:
eperm = epermutation
orig_edges = [edge for edge in self.edges]
for ie_old in range(0, self.ne):
ie = eperm[ie_old]
self.edges[ie] = orig_edges[ie_old]
for iv in range(0, self.nv):
for j in range(0, len(self.neighbors[iv])):
je_old = self.neighbors[iv][j]
self.neighbors[iv][j] = eperm[je_old]
def SortNeighborLists(self):
assert(self.nv == len(self.neighbors))
for iv in range(0, self.nv):
#Back when self.neighbors was just a 2-dimensional list of
#vertex id numbers, then the following line would have worked:
# self.neighbors[iv].sort()
# ugly python code alert:
#Unfortunately, we had to change the format of self.neighbors. Now
#it is a list of indices into the self.edges array ("ie" numbers).
#We want to sort the "ie" numbers by the vertices they point to.
#self.edge[ie].start should point to the current vertex (hopefully).
#self.edge[ie].stop should point to the vertex it's attached to.
#So we want to sort the ie's in self.neighbors by self.edge[ie].stop
#Create a temporary array of 2-tuples (ie, jv)
nlist = [(ie, self.edges[ie].stop)
for ie in self.neighbors[iv]]
self.neighbors[iv] = [ie for ie,jv in sorted(nlist,
key=itemgetter(1))]
def FindEdge(self, istart, istop):
"""
A simple function looks up the edge id number
corresponding to an edge connecting vertex istart to istop.
If not present returns Dgraph.NULL.
"""
iv = istart
for je in self.neighbors[iv]:
jv = self.edges[je].stop
if jv == istop:
return je
return Dgraph.NULL
def GetVert(self, iv):
return self.verts[iv]
def GetEdge(self, ie):
return self.edges[ie]
def GetNumVerts(self):
return self.nv
def GetNumEdges(self):
return self.ne
# Commenting out. I think it's clearer to use python's deepcopy instead
#def makecopy(self):
# new_copy = Ugraph()
# new_copy.verts = [vertex for vertex in self.verts]
# new_copy.edges = [ edge for edge in self.edges]
# new_copy.neighbors = [nlist for nlist in self.neighbors]
# new_copy.nv = self.nv
# new_copy.ne = self.ne
# return new_copy
def __str__(self):
# Print the graph as a list of neighbor-lists.
# (Note: This is the same format as the first argument to __init__().
# The Vertex.attr and Edge.attr attributes are not printed.)
l = ['([']
for iv in range(0, self.nv):
l.append('[')
for j in range(0, len(self.neighbors[iv])):
je = self.neighbors[iv][j]
jv = self.edges[je].stop
l.append(str(jv))
if j < len(self.neighbors[iv])-1:
l.append(', ')
else:
l.append(']')
if iv < self.nv-1:
l.append(',\n ')
else:
l.append(']')
l.append(',\n [')
for ie in range(0, self.ne):
l.append(str(self.edges[ie]))
if ie < self.ne-1:
l.append(', ')
else:
l.append('])\n')
return ''.join(l)
def __repr__(self):
return str(self)
class Ugraph(Dgraph):
"""
This class is a minimal implementation of an undirected graph.
Vertices and edges are accessed by integer index only (beginning at 0).
Multiple edges connecting the same pair of vertices are allowed.
(One would use the AddEdge() member function to accomplish this.)
Both vertices and edges have an optional "attr" attribute.
Undirected graphs (Ugraphs) are represented internally as
directed graphs. This means that for every edge in the Ugraph,
connecting vertex 2 to 3, for example, two edges are stored
internally, (2 -> 3, and 3 -> 2),
Edges which begin and end at the same vertex are stored only once.)
"""
def __init__(self, edgelist=None):
Dgraph.__init__(self, edgelist)
# Now add the extra edges which point in the reverse direction.
neu = self.ne
ned = self.ne
for ieu in range(0, self.ne):
iv = self.edges[ieu].start
jv = self.edges[ieu].stop
if iv != jv:
ned += 1
self.ieu_to_ied = [Dgraph.NULL for ieu in range(0, neu)]
self.ied_to_ieu = [Dgraph.NULL for ied in range(0, ned)]
ied_redundant = neu
for ie in range(0, neu):
iv = self.edges[ie].start
jv = self.edges[ie].stop
attr = self.edges[ie].attr
self.ieu_to_ied[ie] = ie
self.ied_to_ieu[ie] = ie
if iv != jv:
# Then create another edge which points in the reverse direction
Dgraph.AddEdge(self, jv, iv, attr) # <--this increments self.ne
self.ied_to_ieu[ied_redundant] = ie
ied_redundant += 1
self.neu = neu
assert(self.ne == ned)
def AddEdge(self, iv, jv, attr=None, remove_duplicates=False):
"""
Add an edge to an undirected graph connecting vertices iv and jv.
If the edge is already present (and remove_duplicates==True),
no new edge will be added.
Note: Undirected Ugraphs are implemented by creating two separate
digraph edges that conect iv->jv and jv->iv.
"""
self.ieu_to_ied.append( len(self.edges) )
Dgraph.AddEdge(self, iv, jv, attr, remove_duplicates)
self.ied_to_ieu.append( self.neu )
if jv != iv:
Dgraph.AddEdge(self, jv, iv, attr, remove_duplicates)
self.ied_to_ieu.append( self.neu )
self.neu += 1
assert(len(self.ieu_to_ied) == self.neu)
assert(len(self.ied_to_ieu) == len(self.edges))
def ReorderEdges(self, epermutation, invert=False):
Dgraph.ReorderEdges(self, epermutation, invert)
# Now update the
# self.ieu_to_ied and
# self.ied_to_ieu lookup tables:
if (invert): # (first invert the permutation if necessary)
eperm = [-1 for ie in epermutation]
for ie in range(0, self.ne):
eperm[ epermutation[ie] ] = ie
else:
eperm = epermutation
#epermutation.reverse()
ieu_to_ied_orig = [ied for ied in self.ieu_to_ied]
ied_to_ieu_orig = [ieu for ieu in self.ied_to_ieu]
for ieu in range(0, self.neu):
ied_old = ieu_to_ied_orig[ieu]
ied = eperm[ied_old]
self.ieu_to_ied[ieu] = ied
for ied_old in range(0, self.ne):
ieu = ied_to_ieu_orig[ied_old]
ied = eperm[ied_old]
self.ied_to_ieu[ied] = ieu
eperm = epermutation
def LookupDirectedEdgeIdx(self, ieu):
return self.ieu_to_ied[ieu]
def LookupUndirectedEdgeIdx(self, ied):
return self.ied_to_ieu[ied]
#def GetVert(self, iv): <-- (inherited from parent)
# return self.verts[iv]
def GetEdge(self, ieu):
ied = self.ieu_to_ied[ieu]
return self.edges[ied]
#def GetNumVerts(self): <-- (inherited from parent)
# return self.nv
def GetNumEdges(self):
return self.neu
def FindEdge(self, istart, istop):
"""
A simple function looks up the (undirected) edge id number
corresponding to an edge connecting vertices istart and istop.
If not present returns Dgraph.NULL.
To find the corresponding entry in the self.edges[] list,
you can either:
use the LookupDirectedEdge() lookup function
or
you can use the parent-class' version of this function
Dgraph.FindEdge(self, istart, istop) which returns
this number by default.
"""
ied = Dgraph.FindEdge(self, istart, istop)
ieu = self.LookupUndirectedEdgeIdx(ied)
return ieu
def CalcEdgeLookupTable(self):
"""
COMMENT: THIS NEXT FUNCTION IS PROBABLY NOT NECESSARY AND MIGHT BE
REMOVED AT A LATER TIME WHEN I FIGURE OUT A BETTER WAY.
Because undirected graphs (Ugraphs) are implemented as directed graphs
(Dgraphs) with redundant edges, they may have some extra edges which
the user never explicitly asked for.
There is some confusion about whether the i'th edge refers to
the i'th undirected edge that the user explicitly added, or
the i'th directed edge which is stored internally.
(The number of directed edges is usually twice the number of
edges that the user asked for. But not always, because edges
wich start and end at the same vertex are only represented once.)
This function calculates lookup tables to translate between
the two edge numbering systems:
self.ieu_to_ied[ieu] returns a directed edge id number,
(which is an index into the self.edges list)
corresponding to the ieu'th undirected edge
which was explicitly added by the caller.
self.ied_to_ieu[ied] takes a directed edge id number (ied,
an index into the self.edges list)
and returns the undirected edge number,
which is allways <= ied
"""
self.ieu_to_ied = []
self.ied_to_ieu = [Ugraph.NULL for ied in range(0, self.ne)]
for ied in range(0, self.ne):
iv = self.edges[ied].start
jv = self.edges[ied].stop
ieu = len(self.ieu_to_ied)
self.ied_to_ieu[ied] = ieu
if iv <= jv:
self.ieu_to_ied.append(ied)
def SortVertsByDegree(g):
vert_numneighbors = [(iv, len(g.neighbors[iv])) for iv in range(0, g.nv)]
vert_numneighbors.sort(key=itemgetter(1))
order = [vert_numneighbors[iv][0] for iv in range(0, g.nv)]
g.ReorderVerts(order, invert=True)
class DFS(object):
"""
This class contains a member function (Order()) calculates the order
of vertices visited in a depth-first-search over a connected graph.
"""
def __init__(self, g):
self.g = g
self.sv = 0 #integer sv keeps track of how many vertices visited so far
self.se = 0 #integer se keeps track of how many edges visited so far
self.vvisited=[False for iv in range(0, self.g.nv)] # verts visited
self.vorder =[Dgraph.NULL for iv in range(0, self.g.nv)] #search order
self.evisited=[False for ie in range(0, self.g.ne)] # edges visited
self.eorder =[Dgraph.NULL for ie in range(0, self.g.ne)] #search order
def Reset(self):
self.sv = 0
self.se = 0
for iv in range(0, self.g.nv):
self.vvisited[iv] = False
self.vorder[iv] = Dgraph.NULL
for ie in range(0, self.g.ne):
self.evisited[ie] = False
self.eorder[ie] = Dgraph.NULL
def Order(self, starting_node=0):
"""
VisitOrder(starting_node)
generates a list of integers from 0 to self.g.nv-1 (=#vertices minus 1)
which represents the order in which the vertices would be visited
during a Depth-First-Search.
The first vertex visited is specified by the "starting_node" argument
(an integer (from 0 to g.nv-1)).
"""
self.Reset()
# The first vertex to be visited should be the starting_node
self.vorder[0] = starting_node
self.vvisited[starting_node] = True
self.sv = 1
self._Order(starting_node)
if self.sv != self.g.nv:
raise(Disconnected(self.g, "Error(Order): "+
"The input graph is not connected."))
assert(self.se == self.g.ne)
return ([iv for iv in self.vorder], [ie for ie in self.eorder])
#return self.order
def _Order(self, iv):
"""
_Order() is a recursive function which carries out a
Depth-First-Search over the graph "self.g", starting with vertex iv.
"""
for je in self.g.neighbors[iv]:
jv = self.g.edges[je].stop
if not self.evisited[je]:
self.eorder[self.se] = je
self.se += 1
self.evisited[je] = True
if not self.vvisited[jv]:
self.vorder[self.sv] = jv
self.sv += 1
self.vvisited[jv] = True
self._Order(jv)
def IsConnected(self):
self.Reset()
self._Order(0)
return (self.sv == self.g.nv)
def IsCyclic(self):
"""
IsCyclic() returns True if the graph is cyclic (and connected).
(An exception is raised on disconnected graphs.)
This function quits early as soon as a cycle is found.
"""
self.Reset()
if (type(self.g) is Ugraph):
is_cyclic = self._IsCyclicUgraph(0, Dgraph.NULL)
else:
is_cyclic = self._IsCyclic(0)
if ((self.sv != self.g.nv) and (not is_cyclic)):
raise(Disconnected(self.g, "Error(IsCyclic): "+
"The input graph is not connected."))
return is_cyclic
def _IsCyclicUgraph(self, iv, ivprev):
"""
_IsCyclicUgraph() is a recursive function which carries out a
Depth-First-Search over the graph "self.g" to determine whether the
graph is cyclic. This function works on undirected graphs (Ugraphs).
Indirected graphs (Ugraphs) are a special case.
Ugraphs are implemented by using two (redundant) forward/backward edges
connecting each pair of adjacent vertices. This creates trivial loops.
This version of _IsCyclicUgraph() only counts loops between more
distantly connected vertices.
"""
self.sv += 1
self.vvisited[iv] = True
for je in self.g.neighbors[iv]:
jv = self.g.edges[je].stop
if self.vvisited[jv]:
if jv != ivprev:
return True
elif self._IsCyclicUgraph(jv, iv):
return True
return False
def _IsCyclic(self, iv):
"""
_IsCyclic() is a recursive function which carries out a
Depth-First-Search over the graph "self.g" to determine whether
the graph is cyclic.
This function works on directed graphs.
"""
self.sv += 1
self.vvisited[iv] = True
for je in self.g.neighbors[iv]:
jv = self.g.edges[je].stop
if self.vvisited[jv]:
return True
elif self._IsCyclic(jv):
return True
return False
class GraphMatcher(object):
"""
This class is a variant of the VF2 algorithm for searching
for small connected subgraphs (g) within a larger graph (G).
GraphMatcher works on directed or underected graphs (Dgraph or Ugraph).
This particular version is better optimized for detecting subgraph
isomorphisms between two graphs of highly unequal size. It should be
faster in these situations because, the computation required for
each step is independent of the number of vertices in the larger graph
In the original VF2 algorithm, the computation time for each step
is proportional to the number of vertices in the larger graph.
(The distinction matters when one graph is much smaller than the other.)
Limitations: At the moment, the matching process uses a simple
depth-first-search to search the vertices of the small graph "g".
Hence this approach fails when the smaller graph g is disconnected.
(but it can probably be fixed by picking a different algorithm to search
the small graph).
"""
def __init__(self,
G, # The "big" graph
g): # The little graph (number of vertices in g must be <= G)
self.G = G
self.g = copy.deepcopy(g)
if (type(self.G) is Ugraph):
assert(type(self.g) is Ugraph)
# self.G.CalcEdgeLookupTable() <-- not needed anymore
self.sv = 0
self.se = 0
self.voccupiedG = [False for iv in range(0, G.nv)]
self.eoccupiedG = [False for ie in range(0, G.ne)]
self.G_is_too_small = False
if ((g.nv > G.nv) or
(g.ne > G.ne)):
self.G_is_too_small = True
#raise GenErr('Error: The first argument of GraphMatcher(G,g),\n'+
# ' must be at least as large as the second.')
# The list self.iv_to_Iv is the mapping between the graph vertices.
# Iv is an index into the large graph's list of vertices.
# iv is an index into the small graph's list of vertices.
# The mapping is stored in the iv_to_Iv list.
self.iv_to_Iv = [Dgraph.NULL for Iv in range(0, self.g.nv)]
self.ie_to_Ie = [Dgraph.NULL for Ie in range(0, self.g.ne)]
# (This used to be called "core_2" in the VF2 algorithm)
# Due to the large number of recursion limit
self.old_recursion_limit = sys.getrecursionlimit()
expected_max_recursion = self.g.nv
if self.old_recursion_limit < 1.5 * expected_max_recursion:
# Give some breathing room.
sys.setrecursionlimit(int(1.5 * expected_max_recursion))
subgraph_searcher = DFS(self.g)
# Perform a Depth-First-Search on the small graph.
self.vorder_g, self.eorder_g = subgraph_searcher.Order()
# Then re-order the vertices and edgers to
# match the order they were visited.
# Note on permutation order:
# (The DFS.Order() function returns the permutation in this format
# old_index[ new_index ]
# where new_index is the DFS iteration when the vertex/edge was visited
# and old_index is the original vertex/edge order.
# However the ReorderVerts() and ReorderEdges() functions expect
# the permutation to have the opposite order: new_index[ old_index ]
# Hence we set "invert=True", when we invoke these functions.)
self.g.ReorderVerts(self.vorder_g, invert=True)
self.g.ReorderEdges(self.eorder_g, invert=True)
# Initialize state
self.Reset()
def Reset(self):
"""Reinitializes the state of the match-search algorithm.
"""
for iv in range(0, self.g.nv):
self.iv_to_Iv[iv] = Dgraph.NULL
for ie in range(0, self.g.ne):
self.ie_to_Ie[ie] = Dgraph.NULL
for Iv in range(0, self.G.nv):
self.voccupiedG[Iv] = False
for Ie in range(0, self.G.ne):
self.eoccupiedG[Ie] = False
self.se = 0
self.sv = 0
# OPTIONAL: First, do a partial sort for the vertices in the graphs
# based on number of edges emanating from each vertex.
# (This is probably unnecessary for small subgraphs.)
#SortVertsByDegree(self.g)
def Matches(self):
"""
Iterator over all matches between G and g.
Each "match" corresponds to a subgraph of G which is isomorphic to g.
Matches is formatted as a 2-tuple of lists:
(list of vertex ids from G, list of edge ids from G)
The vertex ids in the list are a subset of the integers from 0 to G.nv.
The edge ids in the list are a subset of the integers from 0 to G.ne.
(The corresponding vertices and edges from g are indicated by the order)
"""
self.Reset()
if self.G_is_too_small:
# Then there are fewer verts and edges in G than in g.
# Thus it is impossible for a subgraph of G to be isomorphic to g.
return # return no matches
for Iv in range(0, self.G.nv):
# match vertex Iv from G with vertex 0 from graph g
self.iv_to_Iv[0] = Iv
self.voccupiedG[Iv] = True
# Implementation:
# In this loop we begin the search process
# starting with a different vertex (Iv) from big graph G,
# and matching it with the first vertex (iv=0) from small graph g.
# In this way the match "begins" from vertex Iv in G.
#
# Any matches found which begin from vertex Iv are distinct
# from matches beginning from any other vertex in G.
# Looping over all Iv in G is necessary and sufficient
# to insure that all possible subgraphs of G
# (which are isomorphic to g) are considered.
self.sv = 1 # we have matched one vertex already
self.se = 0 # we haven't matched any edges yet
for match in self.Match():
yield match
self.voccupiedG[Iv] = False
def Match(self):
# self.se represents how many vertices have been matched so far.
# We are done searching if all of the edges from 0 to self.se-1
# from graph g have been selected (matched with edges from graph G).
if self.se == self.g.ne:
# Note: This also gaurantees that all vertices have been visited.
assert(self.sv == self.g.nv)
yield self.ReformatMatch()
else:
# VF2-style recursive loop:
# We know the next edge to be matched is connected to at least
# one previously visited vertex from g which has already been
# been added to the the current match-in-progress.
iv = self.g.edges[self.se].start
Iv = self.iv_to_Iv[iv]
assert(iv < self.sv) # <-- check to verify this is so
# The other vertex may or may not have been visited (matched) yet.
iv_neighbor = self.g.edges[self.se].stop
# Two cases:
# Case 1: edge self.se points to a previously visited vertex from g
# This means we have a loop.
if iv_neighbor < self.sv:
# In that case, then the corresponding edge in G must
# connect the corresponding pair of vertices from G.
# (Which we know have already been assigned to vertices in g
# because both iv and iv_neighbor are < self.sv)
Iv_neighbor = self.iv_to_Iv[iv_neighbor]
# Loop over all of the edges in G which connect this pair
# of vertices (Iv --> Iv_neighbor)
for Je in self.G.neighbors[Iv]:
Jv = self.G.edges[Je].stop
if ((Jv == Iv_neighbor) and
(not self.eoccupiedG[Je])):
# Match edge Je from big graph G with
# edge self.se from small graph g
self.ie_to_Ie[self.se] = Je
self.se += 1
self.eoccupiedG[Je] = True
for match in self.Match():
yield match
self.eoccupiedG[Je] = False
self.se -= 1
self.ie_to_Ie[self.se] = Dgraph.NULL
# Case 2:
else: # this would mean that iv_neighbor >= self.sv
# If iv_neighbor>=self.sv, then this edge points to to a vertex
# in g which has not yet been paired with a vertex from G.
# Loop over all of the edges in G which connect vertex
# Iv from G to new (unvisited) vertices in G
for Je in self.G.neighbors[Iv]:
Jv = self.G.edges[Je].stop
if (not self.voccupiedG[Jv]):
assert(not self.eoccupiedG[Je])
# Match both edge Je with je
# AND vertex Jv with jv
self.ie_to_Ie[self.se] = Je
self.se += 1
self.eoccupiedG[Je] = True
self.iv_to_Iv[self.sv] = Jv
self.sv += 1
self.voccupiedG[Jv] = True
# Then continue the recursion
for match in self.Match():
yield match
self.voccupiedG[Jv] = False
self.sv -= 1
self.iv_to_Iv[self.sv] = Dgraph.NULL
self.eoccupiedG[Je] = False
self.se -= 1
self.ie_to_Ie[self.se] = Dgraph.NULL
def ReformatMatch(self):
# (This is because we are assuming g is connected.
# IT should not have any orphanned vertices.)
# Now return the match:
#
# There are different ways of doing this
# version 1:
#match = (self.iv_to_Iv, self.ie_to_Ie) <-return a pointer to array
# version 2:
#match = ([Iv for Iv in self.iv_to_Iv], <-return a copy of the array
# [Ie for Ie in self.ie_to_Ie])
# version 3:
# Recall that the vertices and edges and g have been re-ordered,
# so sort the list of Iv indices in the order they would be
# matched with the original vertices from the original graph g:
#match = ([self.iv_to_Iv[self.vorder_g[iv]]
# for iv in range(0,self.g.nv)],
# [self.ie_to_Ie[self.eorder_g[ie]]
# for ie in range(0,self.g.ne)])
# version 4: Similar to version 3 above, but we also translate
# the directed edge id list into a shorter undirected
# edge id list.
match_verts = [self.iv_to_Iv[self.vorder_g[iv]]
for iv in range(0,self.g.nv)]
if type(self.g) is Dgraph:
match_edges = [self.ie_to_Ie[self.eorder_g[ie]]
for ie in range(0, self.g.ne)]
else:
#assert(atype(self.g) is Ugraph)
match_edges = [Dgraph.NULL for ieu in range(0, self.g.neu)]
for ie in range(0, self.g.ne):
iv = self.g.edges[ie].start
jv = self.g.edges[ie].stop
if iv <= jv: # <-- avoid duplicating edges (iv,jv) and (jv,iv)
ieu = self.g.LookupUndirectedEdgeIdx(ie)
Ie = self.ie_to_Ie[ie]
Ieu = self.G.LookupUndirectedEdgeIdx(Ie)
match_edges[ieu] = Ieu
return (tuple(match_verts), tuple(match_edges))

74
tools/moltemplate/src/nbody_reorder_atoms.py Executable file
View File

@ -0,0 +1,74 @@
#!/usr/bin/env python
"""
Reorder the atoms in the Angles section of a data file to make sure that
atoms have a "canonical order" (for example the first atom has a lower
id than the last atom, for angle and dihedral interactions.
(This helps us detect potential problems like dupicate Angle interactions.)
"""
import sys
from operator import itemgetter
g_program_name = __file__.split('/')[-1]
in_stream = sys.stdin
section_name = ''
if len(sys.argv) == 2:
section_name = sys.argv[1]
else:
sys.stderr.write('Usage Example:\n\n'
' '+g_program_name+' Angles < angles.txt > new_angles.txt\n\n'
' In this example \"angles.txt\" contains only the \"Angles\" section of\n'
' a LAMMPS DATA file. (Either a text-editor, or the \n'
' \"extract_lammps_data.py\" script can be used to select a section from\n'
' a LAMMPS DATA file\n\n'
'Error('+g_program_name+'): expected exactly one argument:\n'
' \"Angles\", \"Dihedrals\", or \"Impropers\"\n')
exit(-1)
# Ordering rules are defined in a seperate module named
# nbody_Angles.py, nbody_Dihedrals.py, nbody_Impropers.py
# Load that now.
module_name = 'nbody_'+section_name
g = __import__(module_name) #defines g.bond_pattern, g.canonical_order
# This module defines the graph representing the bond pattern for this type
# of interaction. (The number of vertices and edges for the graph corresponds
# to the number of atoms and bonds in this type of interaction.)
natoms = g.bond_pattern.GetNumVerts()
nbonds = g.bond_pattern.GetNumEdges()
for line_orig in in_stream:
line = line_orig.rstrip('\n')
comment = ''
if '#' in line_orig:
ic = line.find('#')
line = line_orig[:ic]
comment = ' '+line_orig[ic:].rstrip('\n')
tokens = line.strip().split()
swapped = False
if len(tokens) == 2+natoms:
all_integers = True
abids_l = [[0 for i in range(0, natoms)],
[0 for i in range(0, nbonds)]]
for i in range(0, natoms):
if not tokens[2+i].isdigit():
all_integers = False
if all_integers:
for i in range(0, natoms):
abids_l[0][i] = int(tokens[2+i])
else:
for i in range(0, natoms):
abids_l[0][i] = tokens[2+i]
abids = g.canonical_order( (tuple(abids_l[0]), tuple(abids_l[1])) )
for i in range(0, natoms):
tokens[2+i] = str(abids[0][i])
sys.stdout.write(' '.join(tokens)+comment+'\n')

258
tools/moltemplate/src/ordereddict.py Normal file
View File

@ -0,0 +1,258 @@
# Backport of OrderedDict() class that runs on Python 2.4, 2.5, 2.6, 2.7 and pypy.
# Passes Python2.7's test suite and incorporates all the latest updates.
try:
from thread import get_ident as _get_ident
except ImportError:
from dummy_thread import get_ident as _get_ident
try:
from _abcoll import KeysView, ValuesView, ItemsView
except ImportError:
pass
class OrderedDict(dict):
'Dictionary that remembers insertion order'
# An inherited dict maps keys to values.
# The inherited dict provides __getitem__, __len__, __contains__, and get.
# The remaining methods are order-aware.
# Big-O running times for all methods are the same as for regular dictionaries.
# The internal self.__map dictionary maps keys to links in a doubly linked list.
# The circular doubly linked list starts and ends with a sentinel element.
# The sentinel element never gets deleted (this simplifies the algorithm).
# Each link is stored as a list of length three: [PREV, NEXT, KEY].
def __init__(self, *args, **kwds):
'''Initialize an ordered dictionary. Signature is the same as for
regular dictionaries, but keyword arguments are not recommended
because their insertion order is arbitrary.
'''
if len(args) > 1:
raise TypeError('expected at most 1 arguments, got %d' % len(args))
try:
self.__root
except AttributeError:
self.__root = root = [] # sentinel node
root[:] = [root, root, None]
self.__map = {}
self.__update(*args, **kwds)
def __setitem__(self, key, value, dict_setitem=dict.__setitem__):
'od.__setitem__(i, y) <==> od[i]=y'
# Setting a new item creates a new link which goes at the end of the linked
# list, and the inherited dictionary is updated with the new key/value pair.
if key not in self:
root = self.__root
last = root[0]
last[1] = root[0] = self.__map[key] = [last, root, key]
dict_setitem(self, key, value)
def __delitem__(self, key, dict_delitem=dict.__delitem__):
'od.__delitem__(y) <==> del od[y]'
# Deleting an existing item uses self.__map to find the link which is
# then removed by updating the links in the predecessor and successor nodes.
dict_delitem(self, key)
link_prev, link_next, key = self.__map.pop(key)
link_prev[1] = link_next
link_next[0] = link_prev
def __iter__(self):
'od.__iter__() <==> iter(od)'
root = self.__root
curr = root[1]
while curr is not root:
yield curr[2]
curr = curr[1]
def __reversed__(self):
'od.__reversed__() <==> reversed(od)'
root = self.__root
curr = root[0]
while curr is not root:
yield curr[2]
curr = curr[0]
def clear(self):
'od.clear() -> None. Remove all items from od.'
try:
for node in self.__map.itervalues():
del node[:]
root = self.__root
root[:] = [root, root, None]
self.__map.clear()
except AttributeError:
pass
dict.clear(self)
def popitem(self, last=True):
'''od.popitem() -> (k, v), return and remove a (key, value) pair.
Pairs are returned in LIFO order if last is true or FIFO order if false.
'''
if not self:
raise KeyError('dictionary is empty')
root = self.__root
if last:
link = root[0]
link_prev = link[0]
link_prev[1] = root
root[0] = link_prev
else:
link = root[1]
link_next = link[1]
root[1] = link_next
link_next[0] = root
key = link[2]
del self.__map[key]
value = dict.pop(self, key)
return key, value
# -- the following methods do not depend on the internal structure --
def keys(self):
'od.keys() -> list of keys in od'
return list(self)
def values(self):
'od.values() -> list of values in od'
return [self[key] for key in self]
def items(self):
'od.items() -> list of (key, value) pairs in od'
return [(key, self[key]) for key in self]
def iterkeys(self):
'od.iterkeys() -> an iterator over the keys in od'
return iter(self)
def itervalues(self):
'od.itervalues -> an iterator over the values in od'
for k in self:
yield self[k]
def iteritems(self):
'od.iteritems -> an iterator over the (key, value) items in od'
for k in self:
yield (k, self[k])
def update(*args, **kwds):
'''od.update(E, **F) -> None. Update od from dict/iterable E and F.
If E is a dict instance, does: for k in E: od[k] = E[k]
If E has a .keys() method, does: for k in E.keys(): od[k] = E[k]
Or if E is an iterable of items, does: for k, v in E: od[k] = v
In either case, this is followed by: for k, v in F.items(): od[k] = v
'''
if len(args) > 2:
raise TypeError('update() takes at most 2 positional '
'arguments (%d given)' % (len(args),))
elif not args:
raise TypeError('update() takes at least 1 argument (0 given)')
self = args[0]
# Make progressively weaker assumptions about "other"
other = ()
if len(args) == 2:
other = args[1]
if isinstance(other, dict):
for key in other:
self[key] = other[key]
elif hasattr(other, 'keys'):
for key in other.keys():
self[key] = other[key]
else:
for key, value in other:
self[key] = value
for key, value in kwds.items():
self[key] = value
__update = update # let subclasses override update without breaking __init__
__marker = object()
def pop(self, key, default=__marker):
'''od.pop(k[,d]) -> v, remove specified key and return the corresponding value.
If key is not found, d is returned if given, otherwise KeyError is raised.
'''
if key in self:
result = self[key]
del self[key]
return result
if default is self.__marker:
raise KeyError(key)
return default
def setdefault(self, key, default=None):
'od.setdefault(k[,d]) -> od.get(k,d), also set od[k]=d if k not in od'
if key in self:
return self[key]
self[key] = default
return default
def __repr__(self, _repr_running={}):
'od.__repr__() <==> repr(od)'
call_key = id(self), _get_ident()
if call_key in _repr_running:
return '...'
_repr_running[call_key] = 1
try:
if not self:
return '%s()' % (self.__class__.__name__,)
return '%s(%r)' % (self.__class__.__name__, self.items())
finally:
del _repr_running[call_key]
def __reduce__(self):
'Return state information for pickling'
items = [[k, self[k]] for k in self]
inst_dict = vars(self).copy()
for k in vars(OrderedDict()):
inst_dict.pop(k, None)
if inst_dict:
return (self.__class__, (items,), inst_dict)
return self.__class__, (items,)
def copy(self):
'od.copy() -> a shallow copy of od'
return self.__class__(self)
@classmethod
def fromkeys(cls, iterable, value=None):
'''OD.fromkeys(S[, v]) -> New ordered dictionary with keys from S
and values equal to v (which defaults to None).
'''
d = cls()
for key in iterable:
d[key] = value
return d
def __eq__(self, other):
'''od.__eq__(y) <==> od==y. Comparison to another OD is order-sensitive
while comparison to a regular mapping is order-insensitive.
'''
if isinstance(other, OrderedDict):
return len(self)==len(other) and self.items() == other.items()
return dict.__eq__(self, other)
def __ne__(self, other):
return not self == other
# -- the following methods are only used in Python 2.7 --
def viewkeys(self):
"od.viewkeys() -> a set-like object providing a view on od's keys"
return KeysView(self)
def viewvalues(self):
"od.viewvalues() -> an object providing a view on od's values"
return ValuesView(self)
def viewitems(self):
"od.viewitems() -> a set-like object providing a view on od's items"
return ItemsView(self)

138
tools/moltemplate/src/pdbsort.py Executable file
View File

@ -0,0 +1,138 @@
#!/usr/bin/env python
"""
Unfortunately, the lines in a PDB files are not always listed in the
correct order. Software that reads PDB files is expected to re-sort this
data before interpreting it. (One reason I don't like them.)
This script reads a PDB file from the standard input, sorts the lines
according to ChainID, SeqNum (residue-number), Icode (insert-code),
and AtomID (in that order) and prints the result to the standard-out.
Only the ATOM and HETATM records are effected.
All other lines in the PDB file are printed back to the user verbatim.
Note: the "altLoc" column (character #17 on each line) is ignored
and is not used for sorting purposes.
"""
from collections import defaultdict
# In order to specify any amino acid in a PDB file, you must provide 3
# identifiers:
# the ChainID a single letter specifying
# the SeqNum an integer indicating the location within that chain
# the ICode ("insert code" usually 0. I don't know why this number is
# necessary. ..For the record, I never loved the PDB file format.)
class AtomDescr:
def __init__(self, setChainID, setSeqNum, setICode, setAtomID):
self.chainID = setChainID
self.seqNum = setSeqNum
self.iCode = setICode
self.atomID = setAtomID
#I plan to store this information in a python dictionary.
#Unfortunately, in order to use such classes as keys in python dictionaries
#I must define comparison operators, and a hash function.
#In retrospect I figured out it would have been easier just to use tuples
#as dictionary keys. I suppose it was a good excercise to try to do it
#using python classes instead. I'm sorry it made the code so long.
def __le__(self, other):
#return ((self.chainID < other.chainID) or ((self.chainID == other.chainID) and ((self.seqNum < other.seqNum) or ((self.seqNum == other.seqNum) and (self.iCode <= other.iCode)))))))
# instead I'll exploit python's ability to compare tuples
return (self.chainID, self.seqNum, self.iCode) <= (other.chainID, other.seqNum, other.iCode)
def __lt__(self, other):
return (self.chainID, self.seqNum, self.iCode) < (other.chainID, other.seqNum, other.iCode)
def __eq__(self, other):
#return ((self.chainID == x.chainID) and (self.seqNum == x.seqNum) and (self.iCode == x.iCode))
return (self.chainID, self.seqNum, self.iCode) == (other.chainID, other.seqNum, other.iCode)
def __ne__(self, other):
return not __eq__(self, other)
def __gt__(self, other):
return not __le__(self, other)
def __ge__(self, other):
return not __lt__(self, other)
def __cmp__(self, other):
if __lt__(self, other):
return -1
elif __gt__(self, other):
return 1
else:
return 0
def __hash__(self):
numChainIDs = 128
numICodes = 128
i = self.seqNum
i *= numChainIDs
i += ord(self.chainID)
i *= numICodes
i += ord(self.iCode)
i *=10
i += self.atomID
return i
import sys
from operator import attrgetter
g_program_name = __file__.split('/')[-1]
g_version_str = 0.1
g_date_str = 2012-5-06
if len(sys.argv) == 1:
use_all_residues = True
elif len(sys.argv) == 7:
use_all_residues = False
first = AtomDescr(sys.argv[1], int(sys.argv[2]), sys.argv[3], 0)
last = AtomDescr(sys.argv[4], int(sys.argv[5]), sys.argv[6], 2147483647)
else:
sys.stderr.write("Error("+g_program_name+"): This program requires either 0 or 6 arguments.\n"
" By default, the the sequence is extracted from the entire PDB file.\n"
" In that case, no arguments are required.\n"
" Alternately, you can limit the selection to a single interval of\n"
" residues from one of the chains in the PDB file.\n"
" To specify an interval, you must passing 6 arguments to this program.\n"
" This program requires a pair of residues to designate the first and\n"
" last members of the interval. Each residue requires 3 identifiers.\n"
" Consequently the six arguments needed are:\n"
"ChainID_first SeqNum_first ICode_first ChainID_last SeqNum_last ICode_last\n")
exit(-1)
atoms2lines = defaultdict(list)
for line in sys.stdin:
if (line[0:6] == "ATOM ") or (line[0:6] == "HETATM"):
atomID = int(line[6:11])
#atomType = line[12:16]
#altLoc = line[16:17]
iCode = line[26:27]
#resType = line[17:20]
chainID = line[21:22]
seqNumStr = line[22:26]
seqNum = int(seqNumStr)
atomdescr = AtomDescr(chainID, int(seqNumStr), iCode, int(atomID))
atoms2lines[atomdescr].append(line.rstrip('\n'))
else:
sys.stdout.write(line)
# Extract an (unordered) list of the atomdescrs of the atoms in the sequence
atomdescrs = [atomdescr for atomdescr in atoms2lines]
# Residues in PDB files are often not listed in order.
# Consequently, we must sort the list by chainID, seqNum, and finnaly iCode:
sequence_of_atomdescrs = sorted(atomdescrs, key=attrgetter('chainID','seqNum','iCode','atomID'))
for atomdescr in sequence_of_atomdescrs:
for line in atoms2lines[atomdescr]:
sys.stdout.write(line+'\n')

144
tools/moltemplate/src/postprocess_input_script.py Executable file
View File

@ -0,0 +1,144 @@
#!/usr/bin/env python
"""
Reorder the integer arguments to the commands in a LAMMPS input
file if these arguments violate LAMMPS order requirements.
We have to do this because the moltemplate.sh script will automatically
assign these integers in a way which may violate these restrictions
and the user has little control over this.
Right now, the only thing this script does is swap the I and J integers in
"pair_coeff I J ..." commands. Later changes may be added.
"""
import sys
lines_orig = []
in_stream = sys.stdin
f = None
fname = None
num_lines_ignore = 0
# Lines from files passed as arguments are read and processed silently.
# (Why? Sometimes it's necessary to read the contents of previous input scripts
# in order to be able to understand a script command which appears later.
# I'm assuming these files will be processed by lammps in the same order. So I
# must insure that moltemplate.sh passes them to this program in that order.
# I'm too lazy to read the "include" commands in input scripts correctly.)
if len(sys.argv) > 1:
for fname in sys.argv[1:]:
f = open(fname, 'r')
in_stream = f
lines_orig += in_stream.readlines()
num_lines_ignore += len(lines_orig)
f.close()
# Lines read from the standard input are read, processed, and printed to stdout
in_stream = sys.stdin
lines_orig += in_stream.readlines()
pair_style_list=[]
swap_occured = False
warn_wildcard = False
i=0
while i < len(lines_orig):
# Read the next logical line
# Any lines ending in '&' should be merged with the next line before breaking
line_orig = ''
while i < len(lines_orig):
line_counter = 1 + i - num_lines_ignore
line_orig += lines_orig[i]
if ((len(line_orig) < 2) or (line_orig[-2:] != '&\n')):
break
i += 1
line = line_orig.replace('&\n','\n').rstrip('\n')
comment = ''
if '#' in line_orig:
ic = line.find('#')
line = line_orig[:ic]
comment = line_orig[ic:] # keep track of comments (put them back later)
tokens = line.strip().split()
if ((len(tokens) >= 2) and (tokens[0] == 'pair_style')):
pair_style_list = tokens[1:]
if ((len(tokens) >= 3) and (tokens[0] == 'pair_coeff')):
if ((tokens[1].isdigit() and (tokens[2].isdigit())) and
(int(tokens[1]) > int(tokens[2]))):
swap_occured = True
tmp = tokens[2]
tokens[2] = tokens[1]
tokens[1] = tmp
if i >= num_lines_ignore:
# polite warning:
sys.stderr.write('swapped pair_coeff order on line '+str(line_counter))
#if (fname != None):
# sys.stderr.write(' of file \"'+fname+'\"')
sys.stderr.write('\n')
# Deal with the "hbond/" pair coeffs.
#
# The hbond/dreiding pair style designates one of the two atom types
# as a donor, and the other as an acceptor (using the 'i','j' flags)
# If swapped atom types eariler, we also need to swap 'i' with 'j'.
#
# If "hbond/dreiding.." pair style is used with "hybrid" or
# "hybrid/overlay" then tokens[3] is the name of the pair style
# and tokens[5] is either 'i' or 'j'.
if len(pair_style_list) > 0:
if ((pair_style_list[0] == 'hybrid') or
(pair_style_list[0] == 'hybrid/overlay')):
if ((tokens[5] == 'i') and (tokens[3][0:6]=='hbond/')):
tokens[5] = 'j'
sys.stderr.write(' (and replaced \"i\" with \"j\")\n')
elif ((tokens[5] == 'j') and (tokens[3][0:6]=='hbond/')):
tokens[5] = 'i'
sys.stderr.write(' (and replaced \"j\" with \"i\")\n')
elif (pair_style_list[0][0:6] == 'hbond/'):
if (tokens[4] == 'i'):
tokens[4] = 'j'
sys.stderr.write(' (and replaced \"i\" with \"j\")\n')
elif (tokens[4] == 'j'):
tokens[4] = 'i'
sys.stderr.write(' (and replaced \"j\" with \"i\")\n')
sys.stdout.write((' '.join(tokens)+comment).replace('\n','&\n')+'\n')
else:
if ((('*' in tokens[1]) or ('*' in tokens[2]))
and
(not (('*' == tokens[1]) and ('*' == tokens[2])))):
warn_wildcard = True
if i >= num_lines_ignore:
sys.stdout.write(line_orig)
else:
if i >= num_lines_ignore:
sys.stdout.write(line_orig)
i += 1
if swap_occured:
sys.stderr.write('\n'
' WARNING: Atom order in some pair_coeff commands was swapped to pacify LAMMPS.\n'
' For some exotic pair_styles such as hbond/dreiding, this is not enough. If you\n'
' use exotic pair_styles, please verify the \"pair_coeff\" commands are correct.\n')
if warn_wildcard:
sys.stderr.write('\n'
' WARNING: The use of wildcard characters (\"*\") in your \"pair_coeff\"\n'
' commands is not recommended.\n'
' (It is safer to specify each interaction pair manually.\n'
' Check every pair_coeff command. Make sure that every atom type in\n'
' the first group is <= atom types in the second group.\n'
' Moltemplate does NOT do this when wildcards are used.)\n'
' If you are using a many-body pair style then ignore this warning.\n')

53
tools/moltemplate/src/remove_duplicate_atoms.py Executable file
View File

@ -0,0 +1,53 @@
#!/usr/bin/env python
"""
Get rid of lines containing duplicate copies of the same atom in the "Atoms"
section of a LAMMPS data file. Duplicate lines which occur later are
preserved and the earlier lines are erased.
The file is read from sys.stdin. This program does not parse the entire
data file. The text from the "Atoms" section of the LAMMPS file must
be extracted in advance before it is sent to this program.)
"""
import sys
in_stream = sys.stdin
f = None
fname = None
if len(sys.argv) == 2:
fname = sys.argv[1]
f = open(fname, 'r')
in_stream = f
atom_ids_in_use = set([])
lines = in_stream.readlines()
# Start at the end of the file and read backwards.
# If duplicate lines exist, eliminate the ones that occur earlier in the file.
i = len(lines)
while i > 0:
i -= 1
line_orig = lines[i]
line = line_orig.rstrip('\n')
if '#' in line_orig:
ic = line.find('#')
line = line_orig[:ic]
tokens = line.strip().split()
if len(tokens) > 0:
atom_id = tokens[0]
if atom_id in atom_ids_in_use:
del lines[i]
else:
atom_ids_in_use.add(atom_id)
else:
del lines[i]
for line in lines:
sys.stdout.write(line)
if f != None:
f.close()

49
tools/moltemplate/src/remove_duplicates_nbody.py Executable file
View File

@ -0,0 +1,49 @@
#!/usr/bin/env python
"""
Get rid of lines containing duplicate bonded nbody interactions in the
corresponding section of a LAMMPS data file (such as bonds, angles,
dihedrals and impropers). Duplicate lines which occur later are
preserved and the earlier lines are erased.
(This program reads from sys.stdin. This program does not parse the entire
data file. The text from the relevant section of the LAMMPS file should be
extracted in advance before it is sent to this program.)
"""
import sys
in_stream = sys.stdin
if len(sys.argv) == 2:
n = int(sys.argv[1])
if (len(sys.argv) != 2) or (n < 1):
sys.stderr.write('Error (remove_duplicates_nbody.py): expected a positive integer argument.\n')
sys.exit(-1)
atom_ids_in_use = set([])
lines = in_stream.readlines()
# Start at the end of the file and read backwards.
# If duplicate lines exist, eliminate the ones that occur earlier in the file.
i = len(lines)
while i > 0:
i -= 1
line_orig = lines[i]
line = line_orig.rstrip('\n')
if '#' in line_orig:
ic = line.find('#')
line = line_orig[:ic]
tokens = line.strip().split()
if len(tokens) == 2+n:
atom_ids = tuple(tokens[2:2+n])
if atom_ids in atom_ids_in_use:
del lines[i]
else:
atom_ids_in_use.add(atom_ids)
elif len(tokens) == 0:
del lines[i]
for line in lines:
sys.stdout.write(line)

69
tools/moltemplate/src/renumber_DATA_first_column.py Executable file
View File

@ -0,0 +1,69 @@
#!/usr/bin/env python
"""
renumber the integers at the beginning of ever line in the file
to make sure these numbers are contiguous.
The file is read from sys.stdin.
This program does not parse an entire LAMMPS data file.
The text from the "Atoms" section
(or "Bonds", or "Angles", or "Dihedrals", or "Impropers" sections)
of the LAMMPS file must be extracted in advance.
"""
import sys
from operator import itemgetter
in_stream = sys.stdin
f = None
fname = None
if len(sys.argv) == 2:
fname = sys.argv[1]
f = open(fname, 'r')
in_stream = f
lines = in_stream.readlines()
column1_iorig_columnsAfter1 = []
# Start at the end of the file and read backwards.
# If duplicate lines exist, eliminate the ones that occur earlier in the file.
i = 0
while i < len(lines):
line_orig = lines[i]
line = line_orig.rstrip('\n')
comment = ''
if '#' in line_orig:
ic = line.find('#')
line = line_orig[:ic]
comment = ' '+line_orig[ic:].rstrip('\n')
tokens = line.strip().split()
if len(tokens) > 0:
if str.isdigit(tokens[0]):
column1 = int(tokens[0])
else:
column1 = tokens[0]
column1_iorig_columnsAfter1.append([column1, i, ' '.join(tokens[1:])+comment])
i += 1
# Sort the list of lines by the first number on each line
column1_iorig_columnsAfter1.sort(key=itemgetter(0))
# Change all of these numbers so that they are consecutive starting at 1
for i in range(0, len(column1_iorig_columnsAfter1)):
column1_iorig_columnsAfter1[i][0] = i+1
# Sort the list of lines so they are back in the original order
column1_iorig_columnsAfter1.sort(key=itemgetter(1))
for i in range(0, len(column1_iorig_columnsAfter1)):
column1 = column1_iorig_columnsAfter1[i][0]
columnsAfter1 = column1_iorig_columnsAfter1[i][2]
sys.stdout.write(str(column1) +' ' + columnsAfter1 + '\n')
if f != None:
f.close()

4717
tools/moltemplate/src/ttree.py Executable file
View File

File diff suppressed because it is too large Load Diff

2074
tools/moltemplate/src/ttree_lex.py Normal file
View File

File diff suppressed because it is too large Load Diff

618
tools/moltemplate/src/ttree_matrix_stack.py Normal file
View File

@ -0,0 +1,618 @@
# Author: Andrew Jewett (jewett.aij at g mail)
# http://www.chem.ucsb.edu/~sheagroup
# License: 3-clause BSD License (See LICENSE.TXT)
# Copyright (c) 2012, Regents of the University of California
# All rights reserved.
from collections import deque
from array import array
from ttree_lex import *
#import sys
def MultMat(dest, A, B):
""" Multiply two matrices together. Store result in "dest".
"""
I = len(A)
J = len(B[0])
K = len(B) # or len(A[0])
for i in range(0,I):
for j in range(0,J):
dest[i][j] = 0.0
for k in range(0,K):
dest[i][j] += A[i][k] * B[k][j]
def MatToStr(M):
strs = []
for i in range(0, len(M)):
for j in range(0, len(M[i])):
strs.append(str(M[i][j])+' ')
strs.append('\n')
return(''.join(strs))
def LinTransform(dest, M, x):
""" Multiply matrix M by 1-dimensioal array (vector) "x" (from the right).
Store result in 1-dimensional array "dest".
In this function, wetreat "x" and "dest" as a column vectors.
(Not row vectors.)
"""
I = len(A)
J = len(x)
for i in range(0,I):
dest[i] = 0.0
for j in range(0,J):
dest[i] += M[i][j] * x[j]
def AffineTransform(dest, M, x):
""" This function performs an affine transformation on vector "x".
Multiply 3-dimensional vector "x" by first three columns of 3x4
matrix M. Add to this the final column of M. Store result in "dest":
dest[0] = M[0][0]*x[0] + M[0][1]*x[1] + M[0][2]*x[2] + M[0][3]
dest[1] = M[1][0]*x[0] + M[1][1]*x[1] + M[1][2]*x[2] + M[1][3]
dest[2] = M[2][0]*x[0] + M[2][1]*x[1] + M[2][2]*x[2] + M[2][3]
"""
D = len(M)
#assert(len(M[0]) == D+1)
for i in range(0,D):
dest[i] = 0.0
for j in range(0,D):
dest[i] += M[i][j] * x[j]
dest[i] += M[i][D] #(translation offset stored in final column)
def AffineCompose(dest, M2, M1):
"""
Multiplication for pairs of 3x4 matrices is technically undefined.
However what we want to do is compose two affine transformations: M1 and M2
3x4 matrices are used to define rotations/translations
x' = M[0][0]*x + M[0][1]*y + M[0][2]*z + M[0][3]
y' = M[1][0]*x + M[1][1]*y + M[1][2]*z + M[1][3]
z' = M[2][0]*x + M[2][1]*y + M[2][2]*z + M[2][3]
We want to create a new 3x4 matrix representing an affine transformation
(M2 M1), defined so that when (M2 M1) is applied to vector x, the result is
M2 (M1 x). In other words:
first, affine transformation M1 is applied to to x
then, affine transformation M2 is applied to (M1 x)
"""
D = len(M1)
#assert(len(M1[0]) == D+1)
#assert(len(M2[0]) == D+1)
for i in range(0, D):
dest[i][D] = 0.0
for j in range(0, D+1):
dest[i][j] = 0.0
for k in range(0, D):
dest[i][j] += M2[i][k] * M1[k][j]
dest[i][D] += M2[i][D]
def CopyMat(dest, source):
for i in range(0, len(source)):
for j in range(0, len(source[i])):
dest[i][j] = source[i][j]
class AffineStack(object):
"""
This class defines a matrix stack used to define compositions of affine
transformations of 3 dimensional coordinates (rotation and translation).
Affine transformations are represented using 3x4 matrices.
(Coordinates of atoms are thought of as column vectors: [[x],[y],[z]],
although they are represented internally in the more ordinary way [x,y,z].
To aplly an affine transformation to a vector, multiply the vector
by the matrix, from the left-hand side, as explained in the comments for:
AffineTransform(dest, M, x)
Note: The last column of the 3x4 matrix stores a translational offset.
This bears similarity with the original OpenGL matrix stack
http://content.gpwiki.org/index.php/OpenGL:Tutorials:Theory
(OpenGL uses 4x4 matrices. We don't need the final row of these matrices,
because in OpenGL, these rows are used for perspective transformations.)
http://en.wikipedia.org/wiki/Homogeneous_coordinates#Use_in_computer_graphics
"""
def __init__(self):
self.stack = None
self.M = None
self._tmp = None
self.Clear()
def Clear(self):
self.stack = deque([])
self.M = [[1.0, 0.0, 0.0, 0.0],
[0.0, 1.0, 0.0, 0.0],
[0.0, 0.0, 1.0, 0.0]] # (identity, initially)
self._tmp = [[1.0, 0.0, 0.0, 0.0],
[0.0, 1.0, 0.0, 0.0],
[0.0, 0.0, 1.0, 0.0]]
def PushRight(self, M):
#Push a copy of matrix self.M onto the stack
# We make no distinction between "right" and "left" here.
# All transformations are pushed onto the stack in the same way.
# (The "right" and "left" refer to whether the matrix is multiplied
# on the right of left hand side. Because not all matrices need be
# invertible, we require that matrices be popped from the stack
# in the reverse order they were pushed. This prevents the ability
# to push and pop matrices to either end of the stack in an arbitrary
# order (like append(), appendleft(), pop(), popleft()).)
self.stack.append([[self.M[i][j] for j in range(0,len(self.M[i]))]
for i in range(0,len(self.M))])
# The "Right" and "Left" refer to whether the new matrix is multiplied
# on the right or left side of the culmulatie matrix product.
AffineCompose(self._tmp, self.M, M) # Afterwards, self._tmp = self.M * M
#sys.stderr.write('DEBUG: PushLeft()\n' +
# MatToStr(self._tmp) + '\n = \n' +
# MatToStr(M) + '\n * \n' +
# MatToStr(self.M) + '\n')
CopyMat(self.M, self._tmp) # Copy self._tmp into self.M
def PushLeft(self, M):
#Push a copy of matrix self.M onto the stack
# We make no distinction between right and left here.
# All transformations are pushed onto the stack in the same way.
# (The "right" and "left" refer to whether the matrix is multiplied
# on the right of left hand side. Because not all matrices need be
# invertible, we require that matrices be popped from the stack
# in the reverse order they were pushed. This prevents the ability
# to push and pop matrices to either end of the stack in an arbitrary
# order (like append(), appendleft(), pop(), popleft()).)
self.stack.append([[self.M[i][j] for j in range(0,len(self.M[i]))]
for i in range(0,len(self.M))])
# The "Right" and "Left" refer to whether the new matrix is multiplied
# on the right or left side of the culmulatie matrix product.
AffineCompose(self._tmp, M, self.M) # Afterwards, self._tmp = M * self.M
#sys.stderr.write('DEBUG: PushLeft()\n' +
# MatToStr(self._tmp) + '\n = \n' +
# MatToStr(M) + '\n * \n' +
# MatToStr(self.M) + '\n')
CopyMat(self.M, self._tmp) # Copy self.tmp into self.M
def Pop(self):
CopyMat(self.M, self.stack.pop())
# (No need to return a matrix,"self.M",after popping.
# The caller can directly access self.M later.)
#return self.M
def PopRight(self):
self.Pop()
def PopLeft(self):
self.Pop()
def PushCommandsRight(self,
text, # text containing affine transformation commands
# The next two arguments are optional:
src_loc = OSrcLoc(), # for debugging
xcm = None): # position of center of object
"""Generate affine transformation matrices from simple text commands
(such as \"rotcm(90,0,0,1)\" and \"move(0,5.0,0)".
Chains of "rotcm", "movecm", "rot", and "move" commands
can also be strung together:
\"rotcm(90,0,0,1).move(0,5.0,0)\"
Commands ending in \"cm\" are carried out relative to center-of-mass
(average position) of the object, and consequently require
an additional argument (\"xcm\").
"""
self.PushRight(AffineStack.CommandsToMatrix(text, src_loc, xcm))
def PushCommandsLeft(self,
text, # text containing affine transformation commands
# The next two arguments are optional:
src_loc = OSrcLoc(), # for debugging
xcm = None): # position of center of object
self.PushLeft(AffineStack.CommandsToMatrix(text, src_loc, xcm))
def __len__(self):
return 1 + len(self.stack)
@staticmethod
def CommandsToMatrix(text, # text containing affine transformation commands
src_loc = OSrcLoc(), # for debugging
xcm = None): # position of center of object
Mdest=[[1.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0]]
M =[[1.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0]]
Mtmp =[[1.0, 0.0, 0.0, 0.0], [0.0, 1.0, 0.0, 0.0], [0.0, 0.0, 1.0, 0.0]]
transform_commands = text.split(').')
for transform_str in transform_commands:
if transform_str.find('move(') == 0:
i_paren_close = transform_str.find(')')
if i_paren_close == -1:
i_paren_close = len(transform_str)
args = transform_str[5:i_paren_close].split(',')
if (len(args) != 3):
raise InputError('Error near '+ErrorLeader(src_loc.infile, src_loc.lineno)+':\n'
' Invalid command: \"'+transform_str+'\"\n'
' This command requires 3 numerical arguments.')
M = [[1.0, 0.0, 0.0, float(args[0])],
[0.0, 1.0, 0.0, float(args[1])],
[0.0, 0.0, 1.0, float(args[2])]]
AffineCompose(Mtmp, M, Mdest)
CopyMat(Mdest, Mtmp)
#if transform_str.find('movecm(') == 0:
# # assert(xcm != None)
# i_paren_close = transform_str.find(')')
# if i_paren_close == -1:
# i_paren_close = len(transform_str)
# args = transform_str[8:i_paren_close].split(',')
# if (len(args) != 3):
# raise InputError('Error near '+ErrorLeader(src_loc.infile, src_loc.lineno)+':\n'
# ' Invalid command: \"'+transform_str+'\"\n'
# ' This command requires 3 numerical arguments.')
# M = [[1.0, 0.0, 0.0, float(args[0])-(xcm[0])],
# [0.0, 1.0, 0.0, float(args[1])-(xcm[1])],
# [0.0, 0.0, 1.0, float(args[2])-(xcm[2])]]
# AffineCompose(Mtmp, M, Mdest)
# CopyMat(Mdest, Mtmp)
elif transform_str.find('rot(') == 0:
i_paren_close = transform_str.find(')')
if i_paren_close == -1:
i_paren_close = len(transform_str)
args = transform_str[4:i_paren_close].split(',')
center_v = None
if (len(args) == 7):
center_v = [float(args[4]), float(args[5]), float(args[6])]
elif (len(args) != 4):
raise InputError('Error near '+ErrorLeader(src_loc.infile, src_loc.lineno)+':\n'
' Invalid command: \"'+transform_str+'\"\n'
' This command requires either 4 or 7 numerical arguments. Either:\n'
' rot(angle, axisX, axisY, axiZ) or \n'
' rot(angle, axisX, axisY, axiZ, centerX, centerY, centerZ)')
M[0][3] = 0.0 #RotMatAXYZ() only modifies 3x3 submatrix of M
M[1][3] = 0.0 #The remaining final column must be zeroed by hand
M[2][3] = 0.0
RotMatAXYZ(M,
float(args[0])*math.pi/180.0,
float(args[1]),
float(args[2]),
float(args[3]))
if (center_v == None):
AffineCompose(Mtmp, M, Mdest)
CopyMat(Mdest, Mtmp)
else:
# Move "center_v" to the origin
moveCentToOrig = [[1.0, 0.0, 0.0, -center_v[0]],
[0.0, 1.0, 0.0, -center_v[1]],
[0.0, 0.0, 1.0, -center_v[2]]]
AffineCompose(Mtmp, moveCentToOrig, Mdest)
CopyMat(Mdest, Mtmp)
# Rotate the coordinates (relative to the origin)
AffineCompose(Mtmp, M, Mdest) # M is the rotation matrix
CopyMat(Mdest, Mtmp)
# Move the origin back to center_v
moveCentBack = [[1.0, 0.0, 0.0, center_v[0]],
[0.0, 1.0, 0.0, center_v[1]],
[0.0, 0.0, 1.0, center_v[2]]]
AffineCompose(Mtmp, moveCentBack, Mdest)
CopyMat(Mdest, Mtmp)
# # elif transform_str.find('rotcm(') == 0:
# # assert(xcm != None)
# # i_paren_close = transform_str.find(')')
# # if i_paren_close == -1:
# # i_paren_close = len(transform_str)
# # args = transform_str[6:i_paren_close].split(',')
# # if (len(args) != 4):
# # raise InputError('Error near '+ErrorLeader(src_loc.infile, src_loc.lineno)+':\n'
# # ' Invalid command: \"'+transform_str+'\"\n'
# # ' This command requires 4 numerical arguments.')
# #
# # moveCMtoOrig = [[1.0, 0.0, 0.0, -xcm[0]],
# # [0.0, 1.0, 0.0, -xcm[1]],
# # [0.0, 0.0, 1.0, -xcm[2]]]
# # AffineCompose(Mtmp, moveCMtoOrig, Mdest)
# # CopyMat(Mdest, Mtmp)
# # M[0][3] = 0.0#RotMatAXYZ() only modifies 3x3 submatrix of M
# # M[1][3] = 0.0#The remaining final column must be zeroed by hand
# # M[2][3] = 0.0
# # RotMatAXYZ(M,
# # float(args[0])*math.pi/180.0,
# # float(args[1]),
# # float(args[2]),
# # float(args[3]))
# # AffineCompose(Mtmp, M, Mdest)
# # CopyMat(Mdest, Mtmp)
# # moveCmBack = [[1.0, 0.0, 0.0, xcm[0]],
# # [0.0, 1.0, 0.0, xcm[1]],
# # [0.0, 0.0, 1.0, xcm[2]]]
# # AffineCompose(Mtmp, moveCmBack, Mdest)
# # CopyMat(Mdest, Mtmp)
elif transform_str.find('scale(') == 0:
i_paren_close = transform_str.find(')')
if i_paren_close == -1:
i_paren_close = len(transform_str)
args = transform_str[6:i_paren_close].split(',')
if (len(args) == 1):
scale_v = [float(args[0]), float(args[0]), float(args[0])]
center_v = [0.0, 0.0, 0.0]
elif (len(args) == 3):
scale_v = [float(args[0]), float(args[1]), float(args[2])]
center_v = [0.0, 0.0, 0.0]
elif (len(args) == 4):
scale_v = [float(args[0]), float(args[0]), float(args[0])]
center_v = [float(args[1]), float(args[2]), float(args[3])]
elif (len(args) == 6):
scale_v = [float(args[0]), float(args[1]), float(args[2])]
center_v = [float(args[3]), float(args[4]), float(args[5])]
else:
raise InputError('Error near '+ErrorLeader(src_loc.infile, src_loc.lineno)+':\n'
' Invalid command: \"'+transform_str+'\"\n'
' This command requires either 1, 3, 4, or 6 numerical arguments. Either:\n'
' scale(ratio), or \n'
' scale(ratioX, ratioY, ratioZ),\n'
' scale(ratio, centerX, centerY, centerZ), or\n'
' scale(ratioX, ratioY, ratioZ, centerX, centerY, centerZ)')
ScaleMat(M, scale_v)
# Now worry about translation:
for d in range(0, 3):
M[d][3] = center_v[d] * (1.0 - scale_v[d])
AffineCompose(Mtmp, M, Mdest)
CopyMat(Mdest, Mtmp)
# # elif transform_str.find('scalecm(') == 0:
# # assert(xcm != None)
# # i_paren_close = transform_str.find(')')
# # if i_paren_close == -1:
# # i_paren_close = len(transform_str)
# # args = transform_str[8:i_paren_close].split(',')
# #
# # moveCMtoOrig = [[1.0, 0.0, 0.0, -xcm[0]],
# # [0.0, 1.0, 0.0, -xcm[1]],
# # [0.0, 0.0, 1.0, -xcm[2]]]
# # AffineCompose(Mtmp, moveCMtoOrig, Mdest)
# # CopyMat(Mdest, Mtmp)
# #
# # M[0][3] = 0.0 #ScaleMat() only modifies 3x3 submatrix of M
# # M[1][3] = 0.0 #The remaining final column must be zeroed by hand
# # M[2][3] = 0.0
# # if (len(args) == 1):
# # ScaleMat(M, args[0])
# # elif (len(args) == 3):
# # ScaleMat(M, args)
# # else:
# # raise InputError('Error near '+ErrorLeader(src_loc.infile, src_loc.lineno)+':\n'
# # ' Invalid command: \"'+transform_str+'\"\n'
# # ' This command requires either 1 or 3 numerical arguments.')
# #
# # AffineCompose(Mtmp, M, Mdest)
# # CopyMat(Mdest, Mtmp)
# # moveCmBack = [[1.0, 0.0, 0.0, xcm[0]],
# # [0.0, 1.0, 0.0, xcm[1]],
# # [0.0, 0.0, 1.0, xcm[2]]]
# # AffineCompose(Mtmp, moveCmBack, Mdest)
# # CopyMat(Mdest, Mtmp)
else:
raise InputError('Error near '+ErrorLeader(src_loc.infile, src_loc.lineno)+':\n'
' Unknown transformation command: \"'+transform_str+'\"\n')
return Mdest
class MultiAffineStack(object):
def __init__(self, which_stack=None):
self.tot_stack = None
self.stack_lookup = None
self.stack_keys = None
self.stacks = None
self.M = None
self.error_if_substack_empty = False
self.Clear()
def Clear(self):
self.tot_stack = AffineStack()
self.stack_lookup = {}
self.stack_keys = deque([])
self.stacks = deque([])
self.M = self.tot_stack.M
self.error_if_substack_empty = False
def _Update(self):
self.tot_stack.Clear()
self.M = self.tot_stack.M
assert(len(self.stacks) > 0)
for stack in self.stacks:
self.tot_stack.PushRight(stack.M)
def PushStack(self, which_stack):
stack = AffineStack()
self.stack_keys.append(which_stack)
self.stack_lookup[which_stack] = stack
self.stacks.append(stack)
self.tot_stack.PushRight(stack.M)
def PopStack(self):
assert(len(self.stacks) > 0)
self.tot_stack.PopRight()
which_stack = self.stack_keys.pop()
del self.stack_lookup[which_stack]
self.stacks.pop()
def Push(self, M, which_stack=None, right_not_left = True):
if len(self.stacks) == 0:
self.PushStack(which_stack)
if which_stack == None:
stack = self.stacks[-1]
if right_not_left:
stack.PushRight(M) # This should copy the matrix M into stack.M
else:
stack.PushLeft(M)
else:
stack = self.stack_lookup[which_stack]
if right_not_left:
stack.PushRight(M)
else:
stack.PushLeft(M)
if stack == self.stacks[-1]:
self.tot_stack.PopRight() #Replace the last matrix on self.tot_stack
# Note: Always use tot_stack.PopRight (even if right_not_left=False)
self.tot_stack.PushRight(stack.M) # with the the updated version.
# Note: We could call self._Update(M) here, but that is slower.
else:
self._Update()
def PushRight(self, M, which_stack=None):
self.Push(M, which_stack, right_not_left=True)
def PushLeft(self, M, which_stack=None):
self.Push(M, which_stack, right_not_left=False)
def PushCommandsRight(self,
text, # text containing affine transformation commands
# The next two arguments are optional:
src_loc = OSrcLoc(), # for debugging
xcm = None,
which_stack=None): # position of center of object
"""Generate affine transformation matrices from simple text commands
(such as \"rotcm(90,0,0,1)\" and \"move(0,5.0,0)".
Chains of "rotcm", "movecm", "rot", and "move" commands
can also be strung together:
\"rotcm(90,0,0,1).move(0,5.0,0)\"
Commands ending in \"cm\" are carried out relative to center-of-mass
(average position) of the object, and consequently require
an additional argument (\"xcm\").
"""
self.PushRight(AffineStack.CommandsToMatrix(text, src_loc, xcm),
which_stack)
def PushCommandsLeft(self,
text, # text containing affine transformation commands
# The next two arguments are optional:
src_loc = OSrcLoc(), # for debugging
xcm = None, # position of center of object
which_stack = None):
self.PushLeft(AffineStack.CommandsToMatrix(text, src_loc, xcm),
which_stack)
def Pop(self, which_stack=None, right_not_left=True):
#empty_stack_error = False
if which_stack == None:
stack = self.stacks[-1]
if len(stack) >= 1:
if right_not_left:
stack.PopRight()
else:
stack.PopLeft()
# Note: We could call self._Update(M) here, but that is slower
self.tot_stack.PopRight() #Replace the last matrix on self.tot_stack
# Note: Always use tot_stack.PopRight (even if right_not_left=False)
self.tot_stack.PushRight(stack.M) # with the the updated version.
else:
assert(False)
# OPTIONAL CODE BELOW AUTOMATICALLY INVOKES self.PopStack() WHEN
# THE stacks[-1].stack IS EMPTY. PROBABLY DOES NOT WORK. IGNORE
# if (not self.error_if_substack_empty):
# if right_not_left:
# assert(len(self.stacks) > 0)
# self.PopStack()
# else:
# assert(False)
# else:
# empty_stack_error = True
else:
stack = self.stack_lookup[which_stack]
if len(stack) > 1:
if right_not_left:
stack.PopRight()
else:
stack.PopLeft()
self._Update()
else:
assert(False)
#empty_stack_error = True
def PopRight(self, which_stack=None):
self.Pop(which_stack, right_not_left=True)
def PopLeft(self, which_stack=None):
self.Pop(which_stack, right_not_left=True)
import math
def ScaleMat(dest, scale):
for i in range(0, len(dest)):
for j in range(0, len(dest[i])):
dest[i][j] = 0.0
if ((type(scale) is float) or (type(scale) is int)):
for i in range(0, len(dest)):
dest[i][i] = scale
else:
for i in range(0, len(dest)):
dest[i][i] = scale[i]
def RotMatAXYZ(dest, angle, axis_x, axis_y, axis_z):
r = math.sqrt(axis_x*axis_x + axis_y*axis_y + axis_z*axis_z)
X = axis_x / r
Y = axis_y / r
Z = axis_z / r
#angle *= math.pi/180.0 # "angle" is assumed to be in degrees
# on second thought, let the caller worry about angle units.
c = math.cos(angle)
s = math.sin(angle)
dest[0][0] = X*X*(1-c) + c
dest[1][1] = Y*Y*(1-c) + c
dest[2][2] = Z*Z*(1-c) + c
dest[0][1] = X*Y*(1-c) - Z*s
dest[0][2] = X*Z*(1-c) + Y*s
dest[1][0] = Y*X*(1-c) + Z*s
dest[2][0] = Z*X*(1-c) - Y*s
dest[1][2] = Y*Z*(1-c) - X*s
dest[2][1] = Z*Y*(1-c) + X*s
# formula from these sources:
# http://inside.mines.edu/~gmurray/ArbitraryAxisRotation/
# also check
# http://www.manpagez.com/man/3/glRotate/
# some pdb test commands:
# from lttree_matrixstack import *
# r = [[1.0,0.0,0.0], [0.0,1.0,0.0], [0.0,0.0,1.0]]
# RotMatAXYZ(r, 90.0, 0.0, 0.0, 1.0)

115
tools/moltemplate/src/ttree_render.py Executable file
View File

@ -0,0 +1,115 @@
#!/usr/bin/env python
man_page_text = """
Usage (example):
ttree_render.py ttree_assignments.txt < file.template > file.rendered
The argument (ttree_assignments.txt) should be a 2-column file containing
ttree-style variables (1st column), and their values (bindings, 2nd column).
This program reads a text file containing ttree-style variables,
substitutes the corresponding values stored in ttree_assignments.txt,
and prints out the new (rendered) text to the standard-out.
"""
import sys
from ttree_lex import *
from ttree import *
import gc
g_filename = __file__.split('/')[-1]
g_module_name = g_filename
if g_filename.rfind('.py') != -1:
g_module_name = g_filename[:g_filename.rfind('.py')]
g_date_str = '2012-9-06'
g_version_str = '0.1'
g_program_name = g_filename
#sys.stderr.write(g_program_name+' v'+g_version_str+' '+g_date_str+' ')
try:
if (len(sys.argv) != 2):
raise InputError('Error running \"'+g_program_name+'\"\n'
' Typical usage:\n'
' ttree_render.py ttree_assignments.txt < file.template > file.rendered\n'
'\n'
' Missing argument.\n'
' Expected the name of a 2-column file containing\n'
' variable names and their bindings (values).\n'
' (This is likely a programmer error.\n'
' This script was not intended to be run by end users.)\n')
bindings_filename = sys.argv[1]
f = open(bindings_filename)
assignments = {}
#BasicUIReadBindingsStream(assignments, f, bindings_filename)
# The line above is robust but it uses far too much memory.
# This for loop below works for most cases.
for line in f:
#tokens = lines.strip().split()
tokens = SplitQuotedString(line.strip()) # like split but handles quotes
if len(tokens) < 2:
continue
assignments[tokens[0]] = tokens[1]
f.close()
gc.collect()
lex = TemplateLexer(sys.stdin, '__standard_input_for_ttree_render__')
lex.var_delim = '$@'
text_block_list = lex.ReadTemplate(simplify_output=True)
output = []
for entry in text_block_list:
assert(isinstance(entry, str))
if ((len(entry) > 1) and (entry[0] in lex.var_delim)):
if '.' in entry:
ic = entry.find('.')
var_name = entry[:ic]
var_suffix = entry[ic:]
else:
var_name = entry
var_suffix = ''
var_name = entry
if var_name not in assignments:
raise(InputError('Error('+g_program_name+')'
#' at '+ErrorLeader(var_ref.src_loc.infile,
# var_ref.src_loc.lineno)+
' unknown variable:\n'
' \"'+var_name+'\"\n'))
else:
var_value = assignments[var_name]
format_fname, args = ExtractFormattingCommands(var_suffix)
if format_fname == 'ljust':
if len(args) == 1:
var_value = var_value.ljust(int(args[0]))
else:
var_value = var_value.ljust(int(args[0]), args[1])
elif format_fname == 'rjust':
if len(args) == 1:
var_value = var_value.rjust(int(args[0]))
else:
var_value = var_value.rjust(int(args[0]), args[1])
output.append(var_value)
else:
output += entry
sys.stdout.write(''.join(output))
except (ValueError, InputError) as err:
sys.stderr.write('\n'+str(err)+'\n')
sys.exit(-1)