#!/bin/env python

#=============================================================================================
# GLOBAL IMPORTS
#=============================================================================================

import os
import sys
import os.path
import copy
import re
import math

import simtk.unit as units
import simtk.openmm

#=============================================================================================
# Parser for OpenMM validation parameter files
#=============================================================================================

class ValidationParameterFileParser(object):
    """Parses OpenMM validation parameter file

    ValidationParameterFileParser reads and parse a validation parameter file.

    """
    def __init__(self, inFilename, verbose=False):
        """

        ARGUMENTS

        inFilename (string) -  OpenMM validation parameter file

        """

        self._name                           = None
        self._version                        = None
        self._particles                      = 0
        self._numberOfMasses                 = 0
        self._masses                         = None
        defaultBoxLength                     = 1000.0
        self._boxVectors                     = [ [ defaultBoxLength, 0.0, 0.0 ], [ 0.0, defaultBoxLength, 0.0 ],[ 0.0, 0.0, defaultBoxLength, ] ]
        self._cMMotionRemover                = 0

        self._harmonicBondForce              = None
        self._harmonicAngleForce             = None
        self._periodicTorsionForce           = None
        self._rbTorsionForce                 = None
        self._nonbondedForce                 = None
        self._gbsaObcForce                   = None
        self._gbviForce                      = None
        self._cMAPTorsionForce               = None
        self._vec3Arrays                     = {}
        self._scalars                        = {}

        self._constraints                    = None

        self._openMM                         = simtk.openmm

        # open file if it exists; otherwise exit

        if os.path.exists( inFilename ):
            filePtr=open(inFilename)
            if verbose: print "Opened file %s" % inFilename
        else:
            print "Could not open file %s" % inFilename
            os.exit(1)

        self._name                           = os.path.basename( inFilename ).replace('.txt','')
        self._path                           = os.path.dirname( inFilename )

        # track line count

        lineCount           = 1

        # prettify diagnostic output

        firstTab            = 30
        secondTab           = 35

        # loop over lines in file
        # blocks of parameters (e.g., charges, sigmas, epsilons) are read in method calls
        # message is output, if a field is not recognized
 
        while True:

            line = filePtr.readline()
            if line is None: break
            if len(line) == 0: break

            # if verbose: print "   Line <%s> at line %d" %  (line, lineCount)

            if line.startswith('Version'):
                tag, self._version = line.rstrip().split(None)
                if verbose:
                   print "   " +  "Version".rjust(firstTab) + self._version.rjust(secondTab) + " at line " + repr( lineCount ).rjust(8) 
                lineCount += 1

            elif line.startswith('Particles'): 
                tag, particles = line.rstrip().split(None)
                self._particles = int(particles)
                if verbose:
                   print "   " +  "Particles".rjust(firstTab) + repr(self._particles).rjust(secondTab) + " at line " + repr( lineCount ).rjust(8) 
                lineCount += 1

            elif line.startswith('Masses'): 
                tag, masses          = line.rstrip().split(None)
                self._numberOfMasses = int(masses)
                if verbose:
                   print "   " +  "Masses".rjust(firstTab) + repr(self._numberOfMasses).rjust(secondTab) + " at line " + repr( lineCount ).rjust(8) 
                lineCount += 1
                lineCount = self._readMasses( filePtr, lineCount )

            elif line.startswith('Box'): 
                tag, a1,a2,a3,b1,b2,b3,c1,c2,c3 = line.rstrip().split(None)

                aArray           = []
                aArray.append( a1 )
                aArray.append( a2 )
                aArray.append( a3 )

                bArray           = []
                bArray.append( b1 )
                bArray.append( b2 )
                bArray.append( b3 )

                cArray           = []
                cArray.append( c1 )
                cArray.append( c2 )
                cArray.append( c3 )

                self._boxVectors = [ aArray, bArray, cArray ]
                if verbose:
                   print "   " +  "Box".rjust(firstTab) + " ".rjust(secondTab) + " at line " + repr( lineCount ).rjust(8) 
                lineCount += 1

            elif line.startswith('NumberOfForces'): 
                lineCount += 1

            elif line.startswith('CMMotionRemover'): 
                tag, cMMotionRemover  = line.rstrip().split(None)
                self._cMMotionRemover = int(cMMotionRemover)
                if verbose:
                   print "   " +  "CMMotionRemover".rjust(firstTab) + repr(self._cMMotionRemover).rjust(secondTab) + " at line " + repr( lineCount ).rjust(8) 
                lineCount            += 1

            elif line.startswith('HarmonicBondForce'): 
                tag, numberOfHarmonicBonds  = line.rstrip().split(None)
                numberOfHarmonicBonds       = int(numberOfHarmonicBonds)
                if verbose:
                   print "   " +  "Harmonic bonds".rjust(firstTab) + repr( numberOfHarmonicBonds).rjust(secondTab) + " at line " + repr( lineCount ).rjust(8) 
                lineCount                  += 1
                lineCount                   = self._readHarmonicBondForce( numberOfHarmonicBonds, filePtr, lineCount )

            elif line.startswith('HarmonicAngleForce'): 
                tag, numberOfHarmonicAngles  = line.rstrip().split(None)
                numberOfHarmonicAngles       = int(numberOfHarmonicAngles)
                if verbose: 
                   print "   " +  "Harmonic angles".rjust(firstTab) + repr( numberOfHarmonicAngles).rjust(secondTab) + " at line " + repr( lineCount ).rjust(8) 
                lineCount                   += 1
                lineCount                    = self._readHarmonicAngleForce( numberOfHarmonicAngles, filePtr, lineCount )

            elif line.startswith('PeriodicTorsionForce'): 
                tag, numberOfPeriodicTorsions  = line.rstrip().split(None)
                numberOfPeriodicTorsions       = int(numberOfPeriodicTorsions)
                if verbose:
                   print "   " +  "Periodic torsions".rjust(firstTab) + repr(numberOfPeriodicTorsions).rjust(secondTab) + " at line " + repr( lineCount ).rjust(8) 
                lineCount                     += 1
                lineCount                      = self._readPeriodicTorsionForce( numberOfPeriodicTorsions, filePtr, lineCount )

            elif line.startswith('RBTorsionForce'): 
                tag, numberOfRBTorsions        = line.rstrip().split(None)
                numberOfRBTorsions             = int(numberOfRBTorsions)
                if verbose:
                   print "   " +  "RB torsions".rjust(firstTab) + repr(numberOfRBTorsions).rjust(secondTab) + " at line " + repr( lineCount ).rjust(8) 
                lineCount                     += 1
                lineCount                      = self._readRBTorsionForce( numberOfRBTorsions, filePtr, lineCount )

            elif line.startswith('NonbondedForceMethod'): 
                tag, nonbondedForceMethod      = line.rstrip().split(None)
                if verbose:
                   print "   " +  "NonbondedForceMethod".rjust(firstTab) + nonbondedForceMethod.rjust(secondTab) + " at line " + repr( lineCount ).rjust(8) 
                lineCount                     += 1
                if self._nonbondedForce:
                    if nonbondedForceMethod == "NoCutoff":
                        nonbondedForceMethod = self._openMM.NonbondedForce.NoCutoff
                    elif nonbondedForceMethod == "CutoffNonPeriodic":
                        nonbondedForceMethod = self._openMM.NonbondedForce.CutoffNonPeriodic
                    elif nonbondedForceMethod == "CutoffPeriodic":
                        nonbondedForceMethod = self._openMM.NonbondedForce.CutoffPeriodic
                    elif nonbondedForceMethod == "Ewald":
                        nonbondedForceMethod = self._openMM.NonbondedForce.Ewald
                    elif nonbondedForceMethod == "PME":
                        nonbondedForceMethod = self._openMM.NonbondedForce.PME
                    else:
                        raise Exception("Nonbonded method not recognized.")

                    self._nonbondedForce.setNonbondedMethod( nonbondedForceMethod )
                else:
                    raise Exception("Cutoff distance set before nonbondedForce instantiated.")

            elif line.startswith('NonbondedForceExceptions'): 
                tag, numberOfNonbondedExceptionsStr  = line.rstrip().split(None)
                numberOfNonbondedExceptions          = int(numberOfNonbondedExceptionsStr)
                if verbose:
                   print "   " +  "Nonbonded exceptions".rjust(firstTab) + repr(numberOfNonbondedExceptions).rjust(secondTab) + " at line " + repr( lineCount ).rjust(8) 
                lineCount                           += 1
                if self._nonbondedForce:
                    lineCount                        = self._readNonbondedForceExceptions( numberOfNonbondedExceptions, filePtr, lineCount )
                else:
                    raise Exception("Nonbonded exceptions set before nonbondedForce instantiated.")

            elif line.startswith('NonbondedForce'): 
                tag, numberOfNonbondeds        = line.rstrip().split(None)
                numberOfNonbondeds             = int(numberOfNonbondeds)
                if verbose:
                   print "   " +  "Nonbonded".rjust(firstTab) + repr(numberOfNonbondeds).rjust(secondTab) + " at line " + repr( lineCount ).rjust(8) 
                lineCount                     += 1
                lineCount                      = self._readNonbondedForce( numberOfNonbondeds, filePtr, lineCount )

            elif line.startswith('CutoffDistance'): 
                tag, cutoffDistanceStr = line.rstrip().split(None)
                cutoffDistance         = float(cutoffDistanceStr)
                if verbose:
                   print "   " +  "CutoffDistance".rjust(firstTab) + cutoffDistanceStr.rjust(secondTab) + " at line " + repr( lineCount ).rjust(8) 
                lineCount                     += 1
                if self._nonbondedForce:
                    self._nonbondedForce.setCutoffDistance( float(cutoffDistanceStr) )
                else:
                    raise Exception("Cutoff distance set before nonbondedForce instantiated.")

            elif line.startswith('RFDielectric'): 
                tag, reactionFieldDielectricStr = line.rstrip().split(None)
                reactionFieldDielectric         = float(reactionFieldDielectricStr)
                if verbose:
                   print "   " +  "RF Dielectric".rjust(firstTab) + ("%10.3f" % reactionFieldDielectric).rjust(secondTab) + " at line " + repr( lineCount ).rjust(8) 
                lineCount                      += 1
                if self._nonbondedForce:
                    self._nonbondedForce.setReactionFieldDielectric( reactionFieldDielectric )
                else:
                    raise Exception("RF dielectric set before nonbondedForce instantiated.")

            elif line.startswith('EwaldRTolerance'): 
                tag, ewaldRToleranceStr = line.rstrip().split(None)
                ewaldRTolerance         = float(ewaldRToleranceStr)
                if verbose:
                   print "   " +  "Ewald tolerance".rjust(firstTab) + ("%12.3e" % ewaldRTolerance).rjust(secondTab) + " at line " + repr( lineCount ).rjust(8) 
                lineCount               += 1
                if self._nonbondedForce:
                    self._nonbondedForce.setEwaldErrorTolerance( ewaldRTolerance )
                else:
                    raise Exception("RF dielectric set before nonbondedForce instantiated.")

            elif line.startswith('GBSAOBCForce'): 
                tag, numberOfGBSAObc         = line.rstrip().split(None)
                numberOfGBSAObc              = int(numberOfGBSAObc)
                if verbose:
                   print "   " +  "GBSA Obc".rjust(firstTab) + repr(numberOfGBSAObc).rjust(secondTab) + " at line " + repr( lineCount ).rjust(8) 
                lineCount                     += 1
                lineCount                      = self._readGBSAObcForce( numberOfGBSAObc, filePtr, lineCount )

            elif line.startswith('SoluteDielectric'): 
                tag, soluteDielectricStr = line.rstrip().split(None)
                soluteDielectric         = float(soluteDielectricStr)
                if verbose:
                   print "   " +  "SoluteDielectric".rjust(firstTab) + ("%10.3f" % soluteDielectric).rjust(secondTab) + " at line " + repr( lineCount ).rjust(8) 
                lineCount                     += 1
                if self._gbsaObcForce:
                    self._gbsaObcForce.setSoluteDielectric( soluteDielectric )
                else:
                    if self._gbviForce:
                        self._gbviForce.setSoluteDielectric( soluteDielectric )
                    else:
                        raise Exception("Solute dielectric set before gbsaObcForce or gbviForce instantiated.")

            elif line.startswith('SolventDielectric'): 
                tag, solventDielectricStr = line.rstrip().split(None)
                solventDielectric         = float(solventDielectricStr)
                if verbose:
                   print "   " +  "SolventDielectric".rjust(firstTab) + ("%10.3f" % solventDielectric).rjust(secondTab) + " at line " + repr( lineCount ).rjust(8) 
                lineCount                     += 1
                if self._gbsaObcForce:
                    self._gbsaObcForce.setSolventDielectric( solventDielectric )
                else:
                    if self._gbviForce:
                        self._gbviForce.setSolventDielectric( soluteDielectric )
                    else:
                        raise Exception("Solvent dielectric set before gbsaObcForce or gbviForce instantiated.")

            elif line.startswith('Constraints'): 
                tag, numberOfConstraintsStr    = line.rstrip().split(None)
                numberOfConstraints            = int(numberOfConstraintsStr)
                if verbose:
                   print "   " +  "Constraints".rjust(firstTab) + repr(numberOfConstraints).rjust(secondTab) + " at line " + repr( lineCount ).rjust(8) 
                lineCount                     += 1
                lineCount                      = self._readConstraints( numberOfConstraints, filePtr, lineCount )

            elif line.startswith('TorsionTorsionParameters'): 
                tag, numberOfTorsionTorsions_Str = line.rstrip().split(None)
                numberOfTorsionTorsions          = int(numberOfTorsionTorsions_Str)
                if verbose:
                   print "   " +  "TorsionTorsions".rjust(firstTab) + repr(numberOfTorsionTorsions).rjust(secondTab) + " at line " + repr( lineCount ).rjust(8) 
                lineCount                     += 1
                lineCount                      = self._readTorsionTorsions( numberOfTorsionTorsions, filePtr, lineCount )

            elif line.startswith('TorsionTorsionGrids'): 
                tag, numberOfTorsionTorsionGrids_Str = line.rstrip().split(None)
                numberOfTorsionTorsionGrids          = int(numberOfTorsionTorsionGrids_Str)
                if verbose:
                   print "   " +  "TorsionTorsionsGrids".rjust(firstTab) + repr(numberOfTorsionTorsionGrids).rjust(secondTab) + " at line " + repr( lineCount ).rjust(8) 
                lineCount                     += 1
                lineCount                      = self._readTorsionTorsionGrids( numberOfTorsionTorsionGrids, filePtr, lineCount )

            else:
                vec3ArrayNames = [ 'Positions', 'Velocities', 'GromacsHarmonicBondForce', 'GromacsHarmonicAngleForce', 'Forces', 'GromacsPeriodicTorsionForce', 'GromacsRBTorsionForce', 'GromacsNonbondedForce', 'GromacsNonbondedForceExceptions', 'GromacsGbsaObcForce', 'TorsionTorsionForce', 'GromacsProperDihedralForce', 'GromacsProperIDihedralForce', 'GromacsGBPolForce', 'Gromacs14Force', 'GromacsCOUL_LRForce', 'GromacsCheckForce', 'GromacsTotalForce' ]

                hit            = 0
                for key in vec3ArrayNames:
                    if( hit == 0 and line.startswith(key) ):
                        hit                           = 1
                        sys.stdout.flush()
                        entries                       = line.rstrip().split(None)
                        tag                           = entries[0]
                        numberOfPositions_Str         = entries[1]
                        numberOfPositions             = int(numberOfPositions_Str)
                        if verbose:
                           print "   " +  tag.rjust(firstTab) + repr(numberOfPositions).rjust(secondTab) + " at line " + repr( lineCount ).rjust(8) 
                        lineCount                     += 1
                        lineCount                      = self._readVec3Array( numberOfPositions, tag, filePtr, lineCount )
        
                if( hit == 0 ):
                    scalarNames = [ 'PotentialEnergy', 'TorsionTorsionEnergy', 'GromacsHarmonicBondEnergy',  'GromacsHarmonicAngleEnergy', 'GromacsProperIDihedralEnergy', 'GromacsProperDihedralEnergy', 'GromacsGBPolEnergy', 'Gromacs14Energy', 'GromacsCOUL_LREnergy', 'GromacsCheckEnergy', 'GromacsTotalEnergy' ]
                    for key in scalarNames:
                        if( hit == 0 and line.startswith(key) ):
                            hit                           = 1
                            tag, value_Str                = line.rstrip().split(None)
                            value                         = float(value_Str)
                            if verbose:
                                print "   " +  key.rjust(firstTab) + repr(value_Str).rjust(secondTab) + " at line " + repr( lineCount ).rjust(8) 
                            lineCount                     += 1
                            self._scalars[tag]             = value

                if( hit == 0 ):
                    print "Line=<%s> at %5d not recognized." % (line.rstrip(), lineCount )
                    lineCount                     += 1

        filePtr.close()
        if verbose:
            lineCount -= 1
            print "Read " + repr( lineCount ) + " lines."

    #=============================================================================================
    # Read mass block 
    #=============================================================================================

    def _readMasses(self, filePtr, lineCount):
        self._masses = []
        for ii in range(self._numberOfMasses):
            line               = filePtr.readline()
            index, mass_String = line.rstrip().split(None)
            #print "   %5d Mass %s at line %d" %  (ii, mas_String, lineCount)
            self._masses.append(  float(mass_String) )
            lineCount        += 1
        return lineCount

    #=============================================================================================
    # Read constraint block 
    #=============================================================================================

    def _readConstraints(self, numberOfConstraints, filePtr, lineCount):
        self._constraints = []
        for ii in range(numberOfConstraints):
            line                  = filePtr.readline()
            indexStr, particle1Str, particle2Str, constrainedDistanceStr = line.rstrip().split(None)
            index                 = int(indexStr)
            particle1             = int(particle1Str)
            particle2             = int(particle2Str)
            constrainedDistance   = float(constrainedDistanceStr)
            #print "   %5d Constraint %5d [%5d %5d] %14.7f at line %5d" %  (ii, index, particle1, particle2, constrainedDistance, lineCount)
            self._constraints.append( [ particle1, particle2, constrainedDistance] )
            lineCount            += 1
        return lineCount

    #=============================================================================================
    # Read force block 
    #=============================================================================================

    def _readVec3Array(self, numberOfEntries, forceName, filePtr, lineCount):
        force = []
        for ii in range(numberOfEntries):
            line                  = filePtr.readline()
            indexStr, xForceStr, yForceStr, zForceStr = line.rstrip().split(None)
            index                 = int(indexStr)
            forceI                = []
            forceI.append( float(xForceStr) )
            forceI.append( float(yForceStr) )
            forceI.append( float(zForceStr) )
            xForce                = float(xForceStr)
            yForce                = float(yForceStr)
            zForce                = float(zForceStr)
            #print "   %5d Forces %5d [%14.7f %14.7f %14.7f] at line %5d" %  (ii, index, xForce, yForce, zForce, lineCount)
            force.append( forceI )
            lineCount            += 1

        self._vec3Arrays[forceName] = force
        return lineCount

    #=============================================================================================
    # Read Torsion-torsion block 
    #=============================================================================================

    def _readTorsionTorsions(self, numberOfTorsionTorsions, filePtr, lineCount):
        self._torsionTorsions = []
        for ii in range(numberOfTorsionTorsions):
            line                  = filePtr.readline()
            #print "<<%s>>" % line
            #sys.stdout.flush()
            indexStr, p1Str, p2Str, p3Str, p4Str, p5Str, chiralAtomStr, gridIndexStr = line.rstrip().split(None)
            index                 = int(indexStr)
            params                = []
            particle1             = int(p1Str)
            particle2             = int(p2Str)
            particle3             = int(p3Str)
            particle4             = int(p4Str)
            particle5             = int(p5Str)
            chiralAtom            = int(chiralAtomStr)
            gridIndex             = int(gridIndexStr)
            print "   %5d TorsionTorsion %5d [%6d %6d %6d %6d %6d ] %6d %6d at line %5d" %  (ii, index, particle1, particle2, particle3, particle4, particle5, chiralAtom, gridIndex, lineCount)
            self._torsionTorsions.append( [ particle1, particle2, particle3, particle4, particle5, chiralAtom, gridIndex ])
            lineCount            += 1
        return lineCount

    #=============================================================================================
    # Read Torsion-torsion grids 
    #=============================================================================================

    def _readTorsionTorsionGrids(self, numberOfTorsionTorsionGrids, filePtr, lineCount):
        self._torsionTorsionGrids = []
        for ii in range(numberOfTorsionTorsionGrids):

            line                        = filePtr.readline()
            str, indexStr, numX1GridStr, numX2GridStr = line.rstrip().split(None)
            numX1Grid                   = int(numX1GridStr)
            numX2Grid                   = int(numX2GridStr)
            grid                        = []
            self._torsionTorsionGridNx1 = numX1Grid
            self._torsionTorsionGridNx2 = numX2Grid
            for jj in range(numX1Grid):
                for kk in range(numX2Grid):
                    line                  = filePtr.readline()
                    indexStr, x1Str, x2Str, x1ValueStr, x2ValueStr, fValueStr, fx1ValueStr, fx2ValueStr, fx1x2ValueStr   = line.rstrip().split(None)
                    index                 = int(indexStr)
                    params                = []
                    x1                    = int(x1Str)
                    x2                    = int(x2Str)
                    x1Value               = float(x1ValueStr)
                    x2Value               = float(x2ValueStr)
                    fValue                = float(fValueStr)
                    fx1Value              = float(fx1ValueStr)
                    fx2Value              = float(fx2ValueStr)
                    fx1x2Value            = float(fx1x2ValueStr)
                    #print "   %5d TorsionTorsion %5d [%6d %6d %6d %6d %6d ] %6d %6d at line %5d" %  (ii, index, particle1, particle2, particle3, particle4, particle5, chiralAtom, gridIndex, lineCount)
                    grid.append( [ x1, x2, x1Value, x2Value, fValue, fx1Value, fx2Value, fx1x2Value ] )
                    lineCount            += 1
            self._torsionTorsionGrids.append( grid )
        return lineCount

    #=============================================================================================
    # Read HarmonicBond block
    #=============================================================================================

    def _readHarmonicBondForce(self, numberOfHarmonicBonds, filePtr, lineCount):
        self._harmonicBondForce = self._openMM.HarmonicBondForce()
        for ii in range(numberOfHarmonicBonds):
            line              = filePtr.readline()
            indexStr, particle1Str, particle2Str, bondLengthStr, kStr = line.rstrip().split(None)
            index             = int(indexStr)
            particle1         = int(particle1Str)
            particle2         = int(particle2Str)
            bondLength        = float(bondLengthStr)
            k                 = float(kStr)
            #print "   %5d HarmonicBondForce %5d [%5d %5d] %14.7f %14.7f at line %5d\n" %  (ii, index, particle1, particle2, bondLength, k, lineCount)
            self._harmonicBondForce.addBond(particle1, particle2, bondLength, k)
            lineCount        += 1
        return lineCount

    #=============================================================================================
    # Read HarmonicAngle block
    #=============================================================================================

    def _readHarmonicAngleForce(self, numberOfHarmonicAngles, filePtr, lineCount):
        self._harmonicAngleForce = self._openMM.HarmonicAngleForce()
        for ii in range(numberOfHarmonicAngles):
            line              = filePtr.readline()
            indexStr, particle1Str, particle2Str, particle3Str, angleStr, kStr = line.rstrip().split(None)
            index             = int(indexStr)
            particle1         = int(particle1Str)
            particle2         = int(particle2Str)
            particle3         = int(particle3Str)
            angle             = float(angleStr)
            k                 = float(kStr)
            #print "   %5d HarmonicAngleForce %5d [%5d %5d %5d] %14.7f %14.7f at line %5d" %  (ii, index, particle1, particle2, particle3, angle, k, lineCount)
            self._harmonicAngleForce.addAngle(particle1, particle2, particle3, angle, k)
            lineCount        += 1
        return lineCount

    #=============================================================================================
    # Read Periodic Torsion block
    #=============================================================================================

    def _readPeriodicTorsionForce(self, numberOfPeriodicTorsions, filePtr, lineCount):
        self._periodicTorsionForce = self._openMM.PeriodicTorsionForce()
        for ii in range(numberOfPeriodicTorsions):
            line              = filePtr.readline()
            indexStr, particle1Str, particle2Str, particle3Str, particle4Str,  periodicityStr, phaseStr, kStr = line.rstrip().split(None)
            index             = int(indexStr)
            particle1         = int(particle1Str)
            particle2         = int(particle2Str)
            particle3         = int(particle3Str)
            particle4         = int(particle4Str)
            periodicity       = int(periodicityStr)
            phase             = float(phaseStr)
            k                 = float(kStr)
            #print "   %5d PeriodicTorsionForce %5d [%5d %5d %5d %5d] %5d %14.7f %14.7f at line %5d" %  (ii, index, particle1, particle2, particle3, particle4, periodicity, phase, k, lineCount)
            self._periodicTorsionForce.addTorsion(particle1, particle2, particle3, particle4, periodicity, phase, k)
            lineCount        += 1
        return lineCount

    #=============================================================================================
    # Read RB Torsion block
    #=============================================================================================

    def _readRBTorsionForce(self, numberOfRBTorsions, filePtr, lineCount):
        self._rbTorsionForce = self._openMM.RBTorsionForce()
        for ii in range(numberOfRBTorsions):
            line              = filePtr.readline()
            indexStr, particle1Str, particle2Str, particle3Str, particle4Str,  c0Str, c1Str, c2Str, c3Str, c4Str, c5Str = line.rstrip().split(None)
            index             = int(indexStr)
            particle1         = int(particle1Str)
            particle2         = int(particle2Str)
            particle3         = int(particle3Str)
            particle4         = int(particle4Str)
            c0                = float(c0Str)
            c1                = float(c1Str)
            c2                = float(c2Str)
            c3                = float(c3Str)
            c4                = float(c4Str)
            c5                = float(c5Str)
            #print "   %5d RBTorsionForce %5d [%5d %5d %5d %5d] [%14.7f %14.7f %14.7f %14.7f %14.7f %14.7f] at line %5d" %  (ii, index, particle1, particle2, particle3, particle4, c0, c1, c2, c3, c4, c5, lineCount)
            self._rbTorsionForce.addTorsion(particle1, particle2, particle3, particle4, c0, c1, c2, c3, c4, c5)
            lineCount        += 1
        return lineCount

    #=============================================================================================
    # Read Nonbonded block
    #=============================================================================================

    def _readNonbondedForce(self, numberOfNonbondeds, filePtr, lineCount):
        self._nonbondedForce = self._openMM.NonbondedForce()
        for ii in range(numberOfNonbondeds):
            line              = filePtr.readline()
            indexStr, chargeStr, sigmaStr, epsilonStr = line.rstrip().split(None)
            index             = int(indexStr)
            charge            = float(chargeStr)
            sigma             = float(sigmaStr)
            epsilon           = float(epsilonStr)
            #print "   %5d NonbondedForce %5d [%14.7f %14.7f %14.7f] at line %5d" %  (ii, index, charge, sigma, epsilon, lineCount)
            self._nonbondedForce.addParticle(charge, sigma, epsilon)
            lineCount        += 1
        return lineCount

    #=============================================================================================
    # Read Nonbonded exceptions block
    #=============================================================================================

    def _readNonbondedForceExceptions(self, numberOfNonbondedExceptions, filePtr, lineCount):
        for ii in range(numberOfNonbondedExceptions):
            line              = filePtr.readline()
            indexStr, particle1Str, particle2Str, chargeProdStr, sigmaStr, epsilonStr = line.rstrip().split(None)
            index             = int(indexStr)
            particle1         = int(particle1Str)
            particle2         = int(particle2Str)
            chargeProd        = float(chargeProdStr)
            sigma             = float(sigmaStr)
            epsilon           = float(epsilonStr)
            #print "   %5d NonbondedExceptionForce %5d [%5d %5d] [%14.7f %14.7f %14.7f] at line %5d" %  (ii, index, particle1, particle2, chargeProd, sigma, epsilon, lineCount)
            self._nonbondedForce.addException(particle1, particle2, chargeProd, sigma, epsilon)
            lineCount        += 1
        return lineCount

    #=============================================================================================
    # Read GBSAObc block
    #=============================================================================================

    def _readGBSAObcForce(self, numberOfGBSAObc, filePtr, lineCount):
        self._gbsaObcForce = self._openMM.GBSAOBCForce()
        for ii in range(numberOfGBSAObc):
            line              = filePtr.readline()
            indexStr, chargeStr, sigmaStr, epsilonStr = line.rstrip().split(None)
            index             = int(indexStr)
            charge            = float(chargeStr)
            sigma             = float(sigmaStr)
            epsilon           = float(epsilonStr)
            #print "   %5d GBSAObcForce %5d [%14.7f %14.7f %14.7f] at line %5d" %  (ii, index, charge, sigma, epsilon, lineCount)
            self._gbsaObcForce.addParticle(charge, sigma, epsilon)
            lineCount        += 1
        return lineCount

    #=============================================================================================
    # Read Gbvi block
    #=============================================================================================

    def _readGbviForce(self, numberOfGbvi, filePtr, lineCount):
        self._gbviForce = self._openMM.GBVIForce()
        for ii in range(numberOfGbvi):
            line              = filePtr.readline()
            indexStr, chargeStr, radiusStr, gammaStr = line.rstrip().split(None)
            index             = int(indexStr)
            charge            = float(chargeStr)
            radius            = float(radiusStr)
            gamma             = float(gammaStr)
            #print "   %5d GBSAObcForce %5d [%14.7f %14.7f %14.7f] at line %5d" %  (ii, index, charge, radius, gamma, lineCount)
            self._gbviForce.addParticle(charge, radius, gamma)
            lineCount        += 1
        return lineCount

    #=============================================================================================
    # Get number of CMAPTorsionForce entries
    #=============================================================================================

    def getNumberOfCMAPTorsions(self):

        if( self._cMAPTorsionForce is None ):
            self.getCMAPTorsionForce()

        if( self._cMAPTorsionForce is None ):
            return 0
        else:
            return self._cMAPTorsionForce.getNumTorsions()

    #=============================================================================================
    # Get CMAPTorsionForce
    #=============================================================================================

    def getCMAPTorsionForce(self):

        if( self._cMAPTorsionForce is not None ):
            return self._cMAPTorsionForce

        self._cMAPTorsionForce = self._openMM.CMAPTorsionForce()
        for ii in range(len(self._torsionTorsions) ):
            parameters = self._torsionTorsions[ii]
            gridIndex  = parameters[6]
            a1         = parameters[0]
            a2         = parameters[1]
            a3         = parameters[2]
            a4         = parameters[3]
            b1         = parameters[1]
            b2         = parameters[2]
            b3         = parameters[3]
            b4         = parameters[4]
            self._cMAPTorsionForce.addTorsion( gridIndex, a1, a2, a3, a4, b1, b2, b3, b4 )

        for ii in range(len(self._torsionTorsionGrids) ):
            grid   = self._torsionTorsionGrids[ii]
            energy = []
            for jj in range(self._torsionTorsionGridNx1*self._torsionTorsionGridNx2 ):
                energy.append( grid[jj][4] )

            #print "<<EArray %s>>" % repr( energy )
            #sys.stdout.flush()
            self._cMAPTorsionForce.addMap(self._torsionTorsionGridNx1, energy) 

        return self._cMAPTorsionForce

    #=============================================================================================
    # Return name
    #=============================================================================================

    def getName(self):
        return self._name

    #=============================================================================================
    # Return array of masses
    #=============================================================================================

    def getMasses(self):
        return self._masses

    #=============================================================================================
    # Return array of masses
    #=============================================================================================

    def getDefaultBoxVectors(self):
        return self._boxVectors

    #=============================================================================================
    # Get CMMotionRemover count
    #=============================================================================================

    def getCMMotionRemover(self):
        return self._cMMotionRemover

    #=============================================================================================
    # Return number of HarmonicBonds
    #=============================================================================================

    def getNumberOfHarmonicBonds(self):
        if( self._harmonicBondForce ):
            return self._harmonicBondForce.getNumBonds()
        else:
            return 0

    #=============================================================================================
    # Return OpenMM HarmonicBondForce
    #=============================================================================================

    def getHarmonicBondForce(self):
        return self._harmonicBondForce;

    #=============================================================================================
    # Return number of HarmonicAngles
    #=============================================================================================

    def getNumberOfHarmonicAngles(self):
        if( self._harmonicAngleForce ):
            return self._harmonicAngleForce.getNumAngles()
        else:
            return 0

    #=============================================================================================
    # Return OpenMM HarmonicAngleForce
    #=============================================================================================

    def getHarmonicAngleForce(self):
        return self._harmonicAngleForce;

    #=============================================================================================
    # Return number of PeriodicTorsions
    #=============================================================================================

    def getNumberOfPeriodicTorsions(self):
        if( self._periodicTorsionForce ):
            return self._periodicTorsionForce.getNumTorsions()
        else:
            return 0

    #=============================================================================================
    # Return OpenMM PeriodicTorsionForce 
    #=============================================================================================

    def getPeriodicTorsionForce(self):
        return self._periodicTorsionForce;

    #=============================================================================================
    # Return number of RBTorsions
    #=============================================================================================

    def getNumberOfRBTorsions(self):
        if( self._rbTorsionForce ):
            return self._rbTorsionForce.getNumTorsions()
        else:
            return 0

    #=============================================================================================
    # Return OpenMM RBTorsionForce
    #=============================================================================================

    def getRBTorsionForce(self):
        return self._rbTorsionForce;

    #=============================================================================================
    # Return number of Nonbondeds
    #=============================================================================================

    def getNumberOfNonbondeds(self):
        if( self._nonbondedForce ):
            return self._nonbondedForce.getNumParticles()
        else:
            return 0

    #=============================================================================================
    # Return OpenMM NonbondedForce
    #=============================================================================================

    def getNonbondedForce(self):
        return self._nonbondedForce;

    #=============================================================================================
    # Return number of gbsaObc particles
    #=============================================================================================

    def getNumberOfGbsaObcs(self):
        if( self._gbsaObcForce):
            return self._gbsaObcForce.getNumParticles()
        else:
            return 0

    #=============================================================================================
    # Get OpenMM GbsaObcForce
    #=============================================================================================

    def getGbsaObcForce(self):
        return self._gbsaObcForce;

    #=============================================================================================
    # Return number of gbvi particles
    #=============================================================================================

    def getNumberOfGbvis(self):
        if( self._gbviForce):
            return self._gbviForce.getNumParticles()
        else:
            return 0

    #=============================================================================================
    # Get OpenMM GbviForce
    #=============================================================================================

    def getGbviForce(self):
        return self._gbviForce;

    #=============================================================================================
    # Return number of constraints
    #=============================================================================================

    def getNumberOfConstraints(self):
        if( self._constraints ):
            return len(self._constraints)
        else:
            return 0

    #=============================================================================================
    # Get constraints array; each entry contains (particle1 index, particle2 index, distance)
    #=============================================================================================

    def getConstraints(self):
        return self._constraints;

    #=============================================================================================
    # Get a vec3 array
    #=============================================================================================

    def getVec3Array(self, name):
        if( name in self._vec3Arrays ): 
            return self._vec3Arrays[name];
        else:
            print "%s array is not available in forces hash." % (name)
            return None

    #=============================================================================================
    # Get a scalar 
    #=============================================================================================

    def getScalar(self, name):
        if( name in self._scalars ): 
            return self._scalars[name];
        else:
            print "%s is not available in scalar hash." % (name)
            return None