Fixed the error on the atoms when we load the same structure multiple…

… times

cellconstructor/Manipulate.py

@@ -188,6 +188,43 @@ def SaveXYZTrajectory(fname, atoms, comments = []):
         atms.save_xyz(fname, comment, overwrite)
 
 
+def load_scf_trajectory(fname):
+    """
+    Load a list of scf structures as a trajectory from the given filename
+
+    Parameters
+    ----------
+        fname : string
+            Path to the file on which the trajectory is saved
+    
+    Results
+    -------
+        trajectory : list of CC.Structure.Structure
+            The list of the structures inside the file name
+    """
+
+    structures = []
+
+    with open(fname, "r") as fp:
+        lines = []
+        for line in fp.readlines():
+            if "CELL_PARAMETERS" in line:
+                if len(lines):
+                    s = Structure()
+                    s.read_scf("".join(lines), read_string = True)
+                    lines = []
+                    structures.append(s)
+            lines.append(line)
+
+        s = Structure()
+        s.read_scf("".join(lines), read_string = True)
+        structures.append(s)
+
+    return structures
+
+
+
+
 def GenerateXYZVideoOfVibrations(dynmat, filename, mode_id, amplitude, dt, N_t, supercell=(1,1,1)):
     """
     XYZ VIDEO

cellconstructor/Structure.py

@@ -299,7 +299,7 @@ def read_xyz(self, filename, alat = False, epsilon = 1e-8, frame_id = 0):
         # Close the xyz file
         xyz.close()
 
-    def read_scf(self, filename, alat=1):
+    def read_scf(self, filename, alat=1, read_string = False):
         """
         Read the given filename in the quantum espresso format.
         Note:
@@ -313,13 +313,21 @@ def read_scf(self, filename, alat=1):
                If present the system will convert both the cell and the atoms position
                by this factor. If it is also specified in the CELL_PARAMETERS line,
                the one specified in the file will be used.
+            - read_string : bool
+                If true the filename is interpreted directly as a scf string, not as a file to be read.
         """
-        # Check if the specified filename exists
-        if not os.path.exists(filename):
-            raise ValueError("File %s does not exist" % filename)
 
-        # Read the input filename
-        fp = open(filename, "r")
+        # Check if the specified filename exists
+        if not read_string:
+            if not os.path.exists(filename):
+                raise ValueError("File %s does not exist" % filename)
+
+            # Read the input filename
+            fp = open(filename, "r")
+            lines = fp.readlines()
+            fp.close()
+        else:
+            lines = filename.split("\n")
 
         n_atoms = 0
         #good_lines = []
@@ -337,7 +345,7 @@ def read_scf(self, filename, alat=1):
         #atom_index = 0
         cell = np.zeros((3,3), dtype = np.float64)
         tmp_coords = []
-        for line in fp.readlines():
+        for line in lines:
             line = line.strip()
 
             #Skipp comments
@@ -362,6 +370,7 @@ def read_scf(self, filename, alat=1):
 
                 continue
             if values[0] == "ATOMIC_POSITIONS":
+                self.atoms = []
                 read_cell = False
                 read_atoms = True
                 if "crystal" in values[1].lower():
@@ -385,7 +394,6 @@ def read_scf(self, filename, alat=1):
                     tmp_coords.append([np.float64(v) for v in values[1:4]])
 
                 n_atoms += 1
-        fp.close()
 
 
 
@@ -1082,7 +1090,7 @@ def set_from_xcoords(self, xcoords):
             self.coords[i,:]  = self.unit_cell.T.dot(xcoords[i,:])
 
 
-    def save_scf(self, filename, alat = 1, avoid_header=False):
+    def save_scf(self, filename, alat = 1, avoid_header=False, crystal = False, get_text = False):
         """
         This methods export the phase in the quantum espresso readable format.
         Of course, only the data reguarding the unit cell and the atomic position will be written.
@@ -1091,13 +1099,17 @@ def save_scf(self, filename, alat = 1, avoid_header=False):
         Parameters
         ----------
             filename : string
-                The name of the file that you want to save.
+                The name of the file that you want to save. If None, no file is generated (in that case, use get_text to get the string of the scf file)
             alat : float, optional
                 If different from 1, both the cell and the coordinates are saved in alat units.
                 It must be in Angstrom.
             avoid_header : bool, optional
                 If true nor the cell neither the ATOMIC_POSITION header is printed.
                 Usefull for the sscha.x code.
+            crystal : bool, optional
+                If true, the atomic coordinates are saved in crystal components
+            get_text : bool
+                If true, the scf file is returned as pure text.
         """
 
         if alat <= 0:
@@ -1121,22 +1133,35 @@ def save_scf(self, filename, alat = 1, avoid_header=False):
 
 
         if not avoid_header:
-            if alat == 1:
-                data.append("ATOMIC_POSITIONS angstrom\n")
+            if crystal:
+                unit_type = "crystal"
             else:
-                data.append("ATOMIC_POSITIONS alat\n")
+                if alat == 1:
+                    unit_type = "angstrom"
+                else:
+                    data.append("ATOMIC_POSITIONS alat\n")
+
+            data.append("ATOMIC_POSITIONS {}\n".format(unit_type))
         for i in range(self.N_atoms):
-            coords = np.copy(self.coords)
-            coords /= alat
+            if not crystal:
+                coords = np.copy(self.coords)
+                coords /= alat
+            else:
+                coords = Methods.covariant_coordinates(self.unit_cell, self.coords)
+
             data.append("%s    %.16f  %.16f  %.16f\n" % (self.atoms[i],
                                                          coords[i, 0],
                                                          coords[i, 1],
                                                          coords[i, 2]))
 
         # Write
-        fdata = open(filename, "w")
-        fdata.writelines(data)
-        fdata.close()
+        if filename is not None:
+            fdata = open(filename, "w")
+            fdata.writelines(data)
+            fdata.close()
+
+        if get_text:
+            return "".join(data)
 
 
     def fix_coords_in_unit_cell(self, delete_copies = True, debug = False):

scripts/view_scf_atoms.py

@@ -0,0 +1,31 @@
+#!python
+
+import cellconstructor as CC, cellconstructor.Structure, cellconstructor.Manipulate
+import ase, ase.visualize
+import sys, os
+
+
+def print_info():
+    TXT = """
+This script is part of CellConstructor.
+Display with ASE the scf file format.
+
+Usage:
+    view_scf_atoms.py filename.scf
+
+filename.scf is the file containing the structure in scf format.
+scf format can be saved from cellconstructor.Structure.Structure class 
+with the method save_scf.
+Or generated by the python-sscha package.
+"""
+    print(TXT)
+
+if __name__ == "__main__":
+    if not len(sys.argv) == 2:
+        print_info()
+        exit()
+
+    all_s = CC.Manipulate.load_scf_trajectory(sys.argv[1])
+    all_atms = [x.get_ase_atoms() for x in all_s]
+
+    ase.visualize.view(all_atms)

setup.py

@@ -71,7 +71,7 @@
        install_requires = ["numpy", "ase", "scipy"],
        license = "MIT",
        include_package_data = True,
-       scripts = ["scripts/symmetrize_dynmat.py", "scripts/cellconstructor_test.py"],
+       scripts = ["scripts/symmetrize_dynmat.py", "scripts/cellconstructor_test.py", "scripts/view_scf_atoms.py"],
        ext_modules = [symph_ext, cc_modules_ext, thirdorder_ext, secondorder_ext]
        )