Add get_reaxff_content

Author: PytLab

vaspy/atomco.py

@@ -17,9 +17,10 @@
 
 import numpy as np
 
-from vaspy import VasPy
-from vaspy.errors import CarfileValueError
-from vaspy.functions import *
+from . import VasPy
+from .errors import CarfileValueError
+from .functions import *
+from .elements import chem_elements
 
 
 class AtomCo(VasPy):
@@ -228,6 +229,38 @@ def get_cif_content(self):
 
         return content
 
+    def get_reaxff_content(self):
+        """
+        Get ReaxFF data file content
+        """
+        content = '# Created by VASPy\n\n'
+
+        # Info
+        content += '{} atoms\n{} atom types\n\n'.format(len(self.data),
+                                                        len(self.atom_types))
+        content += '0 25.000 xlo xhi\n0 25.000 ylo yhi\n0 25.000 zlo zhi\n\n'
+
+        # Masses
+        content += 'Masses\n\n'
+        for i, element in enumerate(self.atom_types):
+            if element not in chem_elements:
+                raise ValueError('element {} not in elements.py'.format(element))
+            mass = chem_elements[element]['mass']
+            content += '{} {:.4f}\n'.format(i+1, mass)
+
+        # Coordinate
+        content += '\nAtoms\n\n'
+        cart_coords = self.dir2cart(self.bases, self.data).tolist()
+
+        for i, (component, coord) in enumerate(zip(self.atom_components, cart_coords)):
+            template = '{:>4d}{:>2d}{:>4.1f}{:>9.5f}{:>11.5f}{:>11.5f}\n'
+            idx = i+1
+            type_idx = self.atom_types.index(component) + 1
+            x, y, z = coord
+            content += template.format(idx, type_idx, 0.0, x, y, z)
+
+        return content
+
     def get_volume(self):
         """
         Get volume of slab(Angstrom^3)

vaspy/elements.py

@@ -0,0 +1,16 @@
+''' Module to store chemical elements related data.
+'''
+C12 = 1.99264648e-26  # mass of C12 (kg/atom)
+amu = 1.66053904e-27  # atomic mass unit (kg)
+
+chem_elements = {}
+chem_elements['H'] = dict(index=1, mass=1.00794)
+chem_elements['C'] = dict(index=6, mass=12.0107)
+chem_elements['N'] = dict(index=7, mass=14.0067)
+chem_elements['O'] = dict(index=8, mass=15.9994)
+chem_elements['F'] = dict(index=9, mass=18.9984032)
+chem_elements['Ne'] = dict(index=10, mass=20.1797)
+chem_elements['S'] = dict(index=16, mass=32.065)
+chem_elements['Cl'] = dict(index=17, mass=35.453)
+chem_elements['Ni'] = dict(index=28, mass=58.693)
+

vaspy/matstudio.py

@@ -120,6 +120,7 @@ def get_atom_info(self):
         tf = []
         tf_dict = {}
         atom_names = []
+        atom_components = []
         atom_name_dict = {}
 
         identity_mappings = self.__get_identity_mappings()
@@ -131,6 +132,7 @@ def get_atom_info(self):
             for elem in identity_mapping.iter('Atom3d'):
                 # Atom name and number
                 atom = elem.attrib['Components']
+                atom_components.append(atom)
                 if atom not in natoms_dict:
                     natoms_dict.setdefault(atom, 1)
                     atom_types.append(atom)
@@ -193,6 +195,7 @@ def get_atom_info(self):
         self.atom_types = atom_types
         self.tf = np.array(tf)
         self.atom_names = atom_names
+        self.atom_components = atom_components
         self.atom_names_dict = atom_name_dict
         self.data = np.array(coordinates)