diff --git a/scripts/vasprun b/scripts/vasprun
index bd63d22..1f7035f 100755
--- a/scripts/vasprun
+++ b/scripts/vasprun
@@ -43,6 +43,22 @@ if __name__ == "__main__":
                       help="figure' dpi, default: 300", metavar="dpi")
     parser.add_option("-S", "--save-bands", dest="saveBands", action='store_true',
                       help="saving of bands contributing to the plot, default: false", metavar="saveBands")
+    parser.add_option("-E", "--elastic", dest="runElasticCalculations", action='store_true',
+                      help="run calculations for the elastic properties (requires IBRION >=5), default: false")
+    parser.add_option("-C", "--elasticsymmetry", dest="symmetry", default="Triclinic",
+                      choices=['Triclinic', 'Monoclinic', 'Orthorhombic',
+                               'Trigonal6', 'Trigonal7', 'Hexagonal',
+                               'Tetragonal6', 'Tetragonal7', 'Cubic'],
+                      help="symmetry group of the supercell, default: Triclinic;\
+                            For calculation of elastic properties only!\
+                      choices: \"Triclinic\", \"Monoclinic\", \"Orthorhombic\",\
+                               \"Trigonal6\", \"Trigonal7\", \"Hexagonal\",\
+                               \"Tetragonal6\", \"Tetragonal7\", and \"Cubic\"")
+    parser.add_option("-A", "--approximation", dest="approximation", default="voigt",
+                      choices=['voigt','reuss','hill','all'],
+                      help="approximation used for the calculation of the elastic moduli, default: voigt;\
+                             For calculation of elastic properties only!\
+                             choices: \"voigt\", \"reuss\", \"hill\", and \"all\"")
 
     (options, args) = parser.parse_args()
     if options.vasprun is None:
@@ -145,3 +161,9 @@ if __name__ == "__main__":
         print("{:25s} {:12.3f} {:12.3f} {:12.3f}".format('DFPT', eps[0,0], eps[0,1], eps[0,2]))
         print("{:25s} {:12.3f} {:12.3f} {:12.3f}".format('DFPT', eps[1,0], eps[1,1], eps[1,2]))
         print("{:25s} {:12.3f} {:12.3f} {:12.3f}".format('DFPT', eps[2,0], eps[2,1], eps[2,2]))
+    elif options.runElasticCalculations:
+        from vasprun.elastic import elastic
+        solver = elastic.Calculator(options=options,verbosity=0)
+        solver() # running calculations
+        solver.print()
+
diff --git a/setup.py b/setup.py
index 035e75d..68561c3 100644
--- a/setup.py
+++ b/setup.py
@@ -17,7 +17,7 @@
     long_description = long_description,
     long_description_content_type = 'text/markdown',
     url="https://github.com/qzhu2017/vasprun",
-    packages=['vasprun'],
+    packages=['vasprun','vasprun.elastic'],
     scripts=['scripts/vasprun'],
     classifiers=[
         "Programming Language :: Python :: 3",
diff --git a/vasprun/elastic/CMatrices.py b/vasprun/elastic/CMatrices.py
new file mode 100644
index 0000000..ed5093a
--- /dev/null
+++ b/vasprun/elastic/CMatrices.py
@@ -0,0 +1,299 @@
+"""
+      Functions allowing for calculation a 6x6 elastic-constant matrix
+   for 9 different point-symmetry space groups. Make sure that you use
+   the proper ordering of the Cij constants!
+"""
+import numpy as np
+
+class ElasticTensor:
+    DegreesOfFreedom = {
+                        'Triclinic'    : 21,
+                        'Monoclinic'   : 13,
+                        'Orthorhombic' : 9,
+                        'Trigonal6'    : 6,
+                        'Trigonal7'    : 7,
+                        'Hexagonal'    : 5,
+                        'Tetragonal6'  : 6,
+                        'Tetragonal7'  : 7,
+                        'Cubic'        : 3,
+                        }
+    def __init__(self,symmetry=None):
+        self.CMatrix = {
+            'Triclinic':    ElasticTensor.CMatrixTriclinic,
+            'Monoclinic':   ElasticTensor.CMatrixMonoclinic,
+            'Orthorhombic': ElasticTensor.CMatrixOrthorhombic,
+            'Trigonal6':    ElasticTensor.CMatrixTrigonal6,
+            'Trigonal7':    ElasticTensor.CMatrixTrigonal7,
+            'Hexagonal':    ElasticTensor.CMatrixHexagonal,
+            'Tetragonal6':  ElasticTensor.CMatrixTetragonal6,
+            'Tetragonal7':  ElasticTensor.CMatrixTetragonal7,
+            'Cubic':        ElasticTensor.CMatrixCubic
+            }
+
+        self.symmetry = None
+        if isinstance(symmetry,str):
+            self.symmetry = symmetry
+
+    def __call__(self,*args,symmetry=None):
+        if symmetry is not None:
+            return self.CMatrix[symmetry](*args)
+        elif self.symmetry is not None:
+            return self.CMatrix[self.symmetry](*args)
+        else:
+            raise TypeError("You must define symmetry!")
+
+    @staticmethod
+    def CMatrixTriclinic(C11,C12,C13,C14,C15,C16,
+                             C22,C23,C24,C25,C26,
+                                 C33,C34,C35,C36,
+                                     C44,C45,C46,
+                                         C55,C56,
+                                             C66):
+        """
+        Triclinic cell - 21 constants
+        Order:
+          C11,C12,C13,C14,C15,C16,
+              C22,C23,C24,C25,C26,
+                  C33,C34,C35,C36,
+                      C44,C45,C46,
+                          C55,C56,
+                              C66
+        returns an array: [
+            [ C11, C12, C13, C14, C15, C16 ],
+            [ C12, C22, C23, C24, C25, C26 ],
+            [ C13, C23, C33, C34, C35, C36 ],
+            [ C14, C24, C34, C44, C45, C46 ],
+            [ C15, C25, C35, C45, C55, C56 ],
+            [ C16, C26, C36, C46, C56, C66 ],
+            ]
+        """
+        return np.array([
+            [ C11, C12, C13, C14, C15, C16 ],
+            [ C12, C22, C23, C24, C25, C26 ],
+            [ C13, C23, C33, C34, C35, C36 ],
+            [ C14, C24, C34, C44, C45, C46 ],
+            [ C15, C25, C35, C45, C55, C56 ],
+            [ C16, C26, C36, C46, C56, C66 ],
+            ])
+    
+    @staticmethod
+    def CMatrixMonoclinic(C11,C12,C13,    C15,
+                              C22,C23,    C25,
+                                  C33,    C35,
+                                      C44,    C46,
+                                          C55,
+                                              C66):
+        """
+        Monoclinic cell - 13 constants
+        Order:
+          C11,C12,C13,C15,
+              C22,C23,C25,
+                  C33,C35,
+                      C44,C46,
+                          C55,
+                              C66
+        returns an array: [
+            [ C11, C12, C13, 0.0, C15, 0.0 ],
+            [ C12, C22, C23, 0.0, C25, 0.0 ],
+            [ C13, C23, C33, 0.0, C35, 0.0 ],
+            [ 0.0, 0.0, 0.0, C44, 0.0, C46 ],
+            [ C15, C25, C35, 0.0, C55, 0.0 ],
+            [ 0.0, 0.0, 0.0, C46, 0.0, C66 ],
+            ]
+        """
+        return np.array([
+            [ C11, C12, C13, 0.0, C15, 0.0 ],
+            [ C12, C22, C23, 0.0, C25, 0.0 ],
+            [ C13, C23, C33, 0.0, C35, 0.0 ],
+            [ 0.0, 0.0, 0.0, C44, 0.0, C46 ],
+            [ C15, C25, C35, 0.0, C55, 0.0 ],
+            [ 0.0, 0.0, 0.0, C46, 0.0, C66 ],
+            ])
+    
+    @staticmethod
+    def CMatrixOrthorhombic(C11,C12,C13,
+                                C22,C23,
+                            C33,C44,C55,C66):
+        """
+        Orthorhombic cell - 9 constants
+        Order:
+          C11,C12,C13,
+              C22,C23,
+          C33,C44,C55,C66
+        returns an array: [
+            [ C11, C12, C13, 0.0, 0.0, 0.0 ],
+            [ C12, C22, C23, 0.0, 0.0, 0.0 ],
+            [ C13, C23, C33, 0.0, 0.0, 0.0 ],
+            [ 0.0, 0.0, 0.0, C44, 0.0, 0.0 ],
+            [ 0.0, 0.0, 0.0, 0.0, C55, 0.0 ],
+            [ 0.0, 0.0, 0.0, 0.0, 0.0, C66 ],
+            ]
+        """
+        return np.array([
+            [ C11, C12, C13, 0.0, 0.0, 0.0 ],
+            [ C12, C22, C23, 0.0, 0.0, 0.0 ],
+            [ C13, C23, C33, 0.0, 0.0, 0.0 ],
+            [ 0.0, 0.0, 0.0, C44, 0.0, 0.0 ],
+            [ 0.0, 0.0, 0.0, 0.0, C55, 0.0 ],
+            [ 0.0, 0.0, 0.0, 0.0, 0.0, C66 ],
+            ])
+    
+    @staticmethod
+    def CMatrixTrigonal7(C11,C12,C13,C14,C15,
+                                 C33,C44):
+        """
+        Trigonal cell - 7 constants
+        Order
+          C11,C12,C13,C14,C15,
+                  C33,C44
+        returns an array: [
+            [ C11, C12, C13, C14, C15, 0.0 ],
+            [ C12, C11, C13,-C14,-C15, 0.0 ],
+            [ C13, C13, C33, 0.0, 0.0, 0.0 ],
+            [ C14,-C14, 0.0, C44, 0.0,-C15 ],
+            [ C15,-C15, 0.0, 0.0, C44, C14 ],
+            [ 0.0, 0.0, 0.0,-C15, C14, C66 ],
+            ]
+        where C66 = 0.5*(C11-C12)
+        """
+        C66 = 0.5*(C11-C12)
+        return np.array([
+            [ C11, C12, C13, C14, C15, 0.0 ],
+            [ C12, C11, C13,-C14,-C15, 0.0 ],
+            [ C13, C13, C33, 0.0, 0.0, 0.0 ],
+            [ C14,-C14, 0.0, C44, 0.0,-C15 ],
+            [ C15,-C15, 0.0, 0.0, C44, C14 ],
+            [ 0.0, 0.0, 0.0,-C15, C14, C66 ],
+            ])
+    
+    @staticmethod
+    def CMatrixTrigonal6(C11,C12,C13,C14,
+                                 C33,C44):
+        """
+        Trigonal cell - 6 constants
+        Order
+          C11,C12,C13,C14,
+                  C33,C44
+        returns an array: [
+            [ C11, C12, C13, C14, 0.0, 0.0 ],
+            [ C12, C11, C13,-C14, 0.0, 0.0 ],
+            [ C13, C13, C33, 0.0, 0.0, 0.0 ],
+            [ C14,-C14, 0.0, C44, 0.0, 0.0 ],
+            [ 0.0, 0.0, 0.0, 0.0, C44, C14 ],
+            [ 0.0, 0.0, 0.0, 0.0, C14, C66 ],
+            ]
+        where C66 = 0.5*(C11-C12)
+        """
+        C66 = 0.5*(C11-C12)
+        return np.array([
+            [ C11, C12, C13, C14, 0.0, 0.0 ],
+            [ C12, C11, C13,-C14, 0.0, 0.0 ],
+            [ C13, C13, C33, 0.0, 0.0, 0.0 ],
+            [ C14,-C14, 0.0, C44, 0.0, 0.0 ],
+            [ 0.0, 0.0, 0.0, 0.0, C44, C14 ],
+            [ 0.0, 0.0, 0.0, 0.0, C14, C66 ],
+            ])
+    
+    @staticmethod
+    def CMatrixHexagonal(C11,C12,C13,
+                                 C33,C44):
+        """
+        Hexagonal cell - 5 constants
+        Order:
+          C11,C12,C13,
+                  C33,C44
+        returns an array: [
+            [ C11, C12, C13, 0.0, 0.0, 0.0 ],
+            [ C12, C11, C13, 0.0, 0.0, 0.0 ],
+            [ C13, C13, C33, 0.0, 0.0, 0.0 ],
+            [ 0.0, 0.0, 0.0, C44, 0.0, 0.0 ],
+            [ 0.0, 0.0, 0.0, 0.0, C44, 0.0 ],
+            [ 0.0, 0.0, 0.0, 0.0, 0.0, C66 ],
+            ]
+        where C66 = 0.5*(C11-C12)
+        """
+        C66 = 0.5*(C11-C12)
+        return np.array([
+            [ C11, C12, C13, 0.0, 0.0, 0.0 ],
+            [ C12, C11, C13, 0.0, 0.0, 0.0 ],
+            [ C13, C13, C33, 0.0, 0.0, 0.0 ],
+            [ 0.0, 0.0, 0.0, C44, 0.0, 0.0 ],
+            [ 0.0, 0.0, 0.0, 0.0, C44, 0.0 ],
+            [ 0.0, 0.0, 0.0, 0.0, 0.0, C66 ],
+            ])
+    
+    @staticmethod
+    def CMatrixTetragonal7(C11,C12,C13,    C16,
+                                   C33,C44,C66):
+        """
+        Tetragonal cell - 7 constants
+        Order:
+          C11,C12,C13,    C16,
+                  C33,C44,C66
+        returns an array: [
+            [ C11, C12, C13, 0.0, 0.0, C16 ],
+            [ C12, C11, C13, 0.0, 0.0,-C16 ],
+            [ C13, C13, C33, 0.0, 0.0, 0.0 ],
+            [ 0.0, 0.0, 0.0, C44, 0.0, 0.0 ],
+            [ 0.0, 0.0, 0.0, 0.0, C44, 0.0 ],
+            [ C16,-C16, 0.0, 0.0, 0.0, C66 ],
+            ]
+        """
+        return np.array([
+            [ C11, C12, C13, 0.0, 0.0, C16 ],
+            [ C12, C11, C13, 0.0, 0.0,-C16 ],
+            [ C13, C13, C33, 0.0, 0.0, 0.0 ],
+            [ 0.0, 0.0, 0.0, C44, 0.0, 0.0 ],
+            [ 0.0, 0.0, 0.0, 0.0, C44, 0.0 ],
+            [ C16,-C16, 0.0, 0.0, 0.0, C66 ],
+            ])
+    
+    @staticmethod
+    def CMatrixTetragonal6(C11,C12,C13,
+                           C33,C44,C66):
+        """
+        Tetragonal cell - 6 constants
+        Order:
+          C11,C12,C13,
+          C33,C44,C66
+        returns an array: [
+            [ C11, C12, C13, 0.0, 0.0, 0.0 ],
+            [ C12, C11, C13, 0.0, 0.0, 0.0 ],
+            [ C13, C13, C33, 0.0, 0.0, 0.0 ],
+            [ 0.0, 0.0, 0.0, C44, 0.0, 0.0 ],
+            [ 0.0, 0.0, 0.0, 0.0, C44, 0.0 ],
+            [ 0.0, 0.0, 0.0, 0.0, 0.0, C66 ],
+            ]
+        """
+        return np.array([
+            [ C11, C12, C13, 0.0, 0.0, 0.0 ],
+            [ C12, C11, C13, 0.0, 0.0, 0.0 ],
+            [ C13, C13, C33, 0.0, 0.0, 0.0 ],
+            [ 0.0, 0.0, 0.0, C44, 0.0, 0.0 ],
+            [ 0.0, 0.0, 0.0, 0.0, C44, 0.0 ],
+            [ 0.0, 0.0, 0.0, 0.0, 0.0, C66 ],
+            ])
+    
+    @staticmethod
+    def CMatrixCubic(C11,C12,C44):
+        """
+        Cubic cell - 3 constants
+        Order:
+          C11,C12,C44
+        returns an array: [
+            [ C11, C12, C12, 0.0, 0.0, 0.0 ],
+            [ C12, C11, C12, 0.0, 0.0, 0.0 ],
+            [ C12, C12, C11, 0.0, 0.0, 0.0 ],
+            [ 0.0, 0.0, 0.0, C44, 0.0, 0.0 ],
+            [ 0.0, 0.0, 0.0, 0.0, C44, 0.0 ],
+            [ 0.0, 0.0, 0.0, 0.0, 0.0, C44 ],
+            ]
+        """
+        return np.array([
+            [ C11, C12, C12, 0.0, 0.0, 0.0 ],
+            [ C12, C11, C12, 0.0, 0.0, 0.0 ],
+            [ C12, C12, C11, 0.0, 0.0, 0.0 ],
+            [ 0.0, 0.0, 0.0, C44, 0.0, 0.0 ],
+            [ 0.0, 0.0, 0.0, 0.0, C44, 0.0 ],
+            [ 0.0, 0.0, 0.0, 0.0, 0.0, C44 ],
+            ])
diff --git a/vasprun/elastic/__init__.py b/vasprun/elastic/__init__.py
new file mode 100644
index 0000000..227e4b7
--- /dev/null
+++ b/vasprun/elastic/__init__.py
@@ -0,0 +1,16 @@
+"""
+   Submodule for calculation of the elastic properties
+   from the stress--strain relation based on VASP's ```vasprun.xml``` file.
+
+   Module CMatrix contains a class ```ElasticTensor``` 
+   with all 9 different symmetry-based constrains 6x6 elastic tensors.
+
+   Module ```elastic``` contains two classes:
+     (i) ```Calculation``` containing the scipy-based engine 
+     fitting the elastic tensor C to the stress-strain relation.
+     (ii) ```Properites``` with the methods deriving the elastic properties
+     from the elastic tensor (within different approximations)
+
+   Module ```readdata``` contains a class ```XMLParser``` extracting
+   both stress and strain from the ```vasprun.xml```
+"""
diff --git a/vasprun/elastic/elastic.py b/vasprun/elastic/elastic.py
new file mode 100644
index 0000000..70c7f01
--- /dev/null
+++ b/vasprun/elastic/elastic.py
@@ -0,0 +1,156 @@
+"""
+   Module ```elastic``` contains two classes:
+     (i) ```Calculation``` containing the scipy-based engine 
+     fiting the elastic tensor C to the stress-strain relation.
+     (ii) ```Properites``` with the methods deriving the elastic properties
+     from the elastic tensor (within different approximations)
+"""
+import numpy as np
+from vasprun.elastic.CMatrices import ElasticTensor
+from vasprun.elastic.readdata import XMLParser
+from scipy.optimize import curve_fit
+from scipy.optimize import least_squares
+
+class Calculator:
+    def __init__(self,options=None,verbosity=0):
+        if options is None:
+            self.opt   = Options()
+            self.opt,_ = self.opt.parse()
+        else:
+            self.opt = options
+        
+        data = XMLParser(self.opt.vasprun)
+        data.parse()
+        self.epsilons,self.sigmas = data.get_strain_stress()
+        
+        if verbosity > 0:
+            print("Using "+self.opt.symmetry+" symmetry")
+        self.ET = ElasticTensor(self.opt.symmetry)
+        if verbosity > 0:
+            print(self.ET.CMatrix[self.opt.symmetry].__doc__)
+
+    def penalty(self,trialState):
+        result = 0.0
+        for e,s in zip(self.epsilons,self.sigmas):
+            result += np.sum(np.power(np.transpose([s])+np.dot(self.ET(*trialState),np.transpose([e])),2))
+        return result
+
+    def __call__(self):
+        solution        = least_squares(self.penalty,1e3*np.ones(ElasticTensor.DegreesOfFreedom[self.opt.symmetry]))
+        self.resultantSystem = Properties(self.ET(*(0.1*solution.x))) #1e-1 -> moving from kBar to GPa
+        return self.resultantSystem
+
+    def print(self,approximation=None):
+        print("Elastic Tensor C (GPa):")
+        np.set_printoptions(suppress=True, precision=0)
+        print(self.resultantSystem.Cij)
+        print("Compliance Tensor s (1/GPa):")
+        np.set_printoptions(suppress=True, precision=5)
+        print(self.resultantSystem.sij)
+        
+        print("Elastic Properties :")
+        print("Approximation Bulk_modulus(GPa) Shear_modulus(GPa) Young_modulus(GPa) Poisson_ratio Cauchy_pressure(GPa) Pugh_ratio")
+        if approximation is None:
+            result = self.get_moduli(self.opt.approximation)
+        else:
+            result = self.get_moduli(approximation)
+        for r in result:
+            print("%13s %12.0f %18.0f %18.0f %17.3f %16.0f %14.3f"%(r,*result[r]))
+
+    def get_moduli(self,approximation='all'):
+        return self.resultantSystem(approximation)
+
+
+class Properties:
+    maskC11 = np.array([
+                       [ 1/3, 0.0, 0.0, 0.0, 0.0, 0.0, ],
+                       [ 0.0, 1/3, 0.0, 0.0, 0.0, 0.0, ],
+                       [ 0.0, 0.0, 1/3, 0.0, 0.0, 0.0, ],
+                       [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, ],
+                       [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, ],
+                       [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, ],
+                       ])
+    maskC12 = np.array([
+                       [ 0.0, 1/6, 1/6, 0.0, 0.0, 0.0, ],
+                       [ 1/6, 0.0, 1/6, 0.0, 0.0, 0.0, ],
+                       [ 1/6, 1/6, 0.0, 0.0, 0.0, 0.0, ],
+                       [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, ],
+                       [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, ],
+                       [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, ],
+                       ])
+    maskC44 = np.array([
+                       [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, ],
+                       [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, ],
+                       [ 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, ],
+                       [ 0.0, 0.0, 0.0, 1/3, 0.0, 0.0, ],
+                       [ 0.0, 0.0, 0.0, 0.0, 1/3, 0.0, ],
+                       [ 0.0, 0.0, 0.0, 0.0, 0.0, 1/3, ],
+                       ])
+    def setCij(self,Cij):
+        try:
+            Cij[5,5] += 0.0
+            self.Cij = Cij
+        except (TypeError,IndexError):
+            try:
+                self.Cij = np.reshape(np.array(Cij),(6,6))
+            except:
+                raise TypeError('elastic matrix must be either numpy.array or an array-like with at least 6x6 fields')
+        self.sij = np.linalg.inv(self.Cij)
+    
+    def __init__(self,Cij):
+        self.setCij(Cij)
+
+    def average(self,mask):
+        return np.sum(mask*self.Cij)
+
+    def __call__(self,approximation="Hill"):
+        self.elastic_properties = {}
+        if approximation[0] in "AaVv":
+            K,G     = self.voigt()
+            E,nu    = Properties.derivative_properties(K,G)
+            Cauchy  = self.average(self.maskC12)
+            Cauchy -= self.average(self.maskC44)
+            Pugh    = G/K
+            self.elastic_properties['Voigt'] = (K,G,E,nu,Cauchy,Pugh)
+        if approximation[0] in "AaRr":
+            K,G     = self.reuss()
+            E,nu    = Properties.derivative_properties(K,G)
+            Cauchy  = self.average(self.maskC12)
+            Cauchy -= self.average(self.maskC44)
+            Pugh    = G/K
+            self.elastic_properties['Reuss'] = (K,G,E,nu,Cauchy,Pugh)
+        if approximation[0] in "AaHh":
+            KV,GV     = self.voigt()
+            KR,GR     = self.reuss()
+            K       = 0.5*(KV+KR)
+            G       = 0.5*(GV+GR)
+            E,nu    = Properties.derivative_properties(K,G)
+            Cauchy  = self.average(self.maskC12)
+            Cauchy -= self.average(self.maskC44)
+            Pugh    = G/K
+            self.elastic_properties['Hill'] = (K,G,E,nu,Cauchy,Pugh)
+        return self.elastic_properties
+
+    def voigt(self):
+        K = ((self.Cij[0,0] + self.Cij[1,1] + self.Cij[2,2]) 
+       + 2.0*(self.Cij[0,1] + self.Cij[1,2] + self.Cij[2,0]))/9.0
+        G = ((self.Cij[0,0] + self.Cij[1,1] + self.Cij[2,2]) 
+       -     (self.Cij[0,1] + self.Cij[1,2] + self.Cij[2,0])
+       + 4.0*(self.Cij[3,3] + self.Cij[4,4] + self.Cij[5,5]))/15.0
+        return K,G
+
+    def reuss(self):
+        K =      np.reciprocal(
+                (self.sij[0,0] + self.sij[1,1] + self.sij[2,2])
+          + 2.0*(self.sij[0,1] + self.sij[1,2] + self.sij[2,0]))
+        G = 15.0*np.reciprocal(
+            4.0*(self.sij[0,0] + self.sij[1,1] + self.sij[2,2]) 
+          - 4.0*(self.sij[0,1] + self.sij[1,2] + self.sij[2,0])
+          + 3.0*(self.sij[3,3] + self.sij[4,4] + self.sij[5,5]))
+        return K,G
+
+    @staticmethod
+    def derivative_properties(K,G):
+        E = 9.0*K*G/(3*K+G)
+        nu= (3*K-2*G)/(6*K+2*G)
+        return E,nu
diff --git a/vasprun/elastic/readdata.py b/vasprun/elastic/readdata.py
new file mode 100644
index 0000000..4fc7b71
--- /dev/null
+++ b/vasprun/elastic/readdata.py
@@ -0,0 +1,65 @@
+"""
+   Module ```readdata``` contains a class ```XMLParser``` extracting
+   both stress and strain from the ```vasprun.xml```
+"""
+import numpy as np
+import defusedxml.ElementTree as ET
+
+class VasprunError(Exception):
+    def __init__(self, data):    
+        self.data = data
+    def __str__(self):
+        return repr(self.data)
+
+class XMLParser:
+    def __init__(self,vasprun='vasprun.xml'):
+        self.vasprun  = 'vasprun.xml'
+        if isinstance(vasprun,str):
+            self.vasprun = vasprun
+        self.sigmas   = None
+        self.epsilons = None
+
+    def parse(self):
+        root   = ET.parse(self.vasprun)
+        # Checking if it makes sense to even start calculations
+        IBRION = 0 # default value
+        for line in root.find('incar').iter('i'):
+            try:
+                if line.attrib['name'] == 'IBRION':
+                    IBRION=int(line.text)
+            except KeyError:
+                continue
+        if IBRION < 5:
+            raise VasprunError("File "+self.vasprun+" does not contain stain-stress calculations!")
+        self.bases  = []
+        self.stress = []
+        for i,calc in enumerate(root.iter('calculation')):
+            basis = []
+            strss = []
+            for varray in calc.iter('varray'):
+                if varray.attrib['name'] == 'stress':
+                    for line in varray:
+                        strss.append(np.fromstring(line.text,sep=' '))
+            self.stress.append(np.array(strss))
+            for varray in calc.find('structure').find('crystal').iter('varray'):
+                if varray.attrib['name'] == 'basis':
+                    for line in varray:
+                        basis.append(np.fromstring(line.text,sep=' '))
+            self.bases.append(np.array(basis))
+        del root
+
+    def get_strain_stress(self):
+        if self.sigmas is None or self.epsilons is None:
+            self.generate_strain_stress()
+        return self.epsilons,self.sigmas
+
+
+    def generate_strain_stress(self):
+        self.sigmas   = []
+        self.epsilons = []
+        for stress,basis in zip(self.stress[1:],self.bases[1:]):
+            strain = np.dot(basis,np.linalg.inv(self.bases[0]))-np.identity(3)
+            sigma  = stress - self.stress[0]
+            if (np.abs(strain) > 1e-6).any():
+                self.sigmas  .append(np.array([ sigma[0,0], sigma[1,1], sigma[2,2], sigma[0,1], sigma[1,2], sigma[2,0]]))
+                self.epsilons.append(np.array([strain[0,0],strain[1,1],strain[2,2],strain[0,1],strain[1,2],strain[2,0]]))