Merge pull request #42 from pyscal/fix_wrong_repeats

Fix wrong repeats

.bumpversion.cfg

@@ -1,5 +1,5 @@
 [bumpversion]
-current_version = 3.2.2
+current_version = 3.2.3
 commit = True
 tag = True
 

setup.py

@@ -7,7 +7,7 @@
 
 setup(
     name='pyscal3',
-    version='3.2.2',
+    version='3.2.3',
     author='Sarath Menon',
     author_email='[email protected]',
     description='Python library written in C++ for calculation of local atomic structural environment',

src/pyscal3/centrosymmetry.cpp

@@ -28,6 +28,7 @@ void calculate_centrosymmetry(py::dict& atoms,
     vector<double> centrosymmetry(nop);
 
     for (int ti=0; ti<nop; ti++){
+		temp.clear();
     	int count = 0;
     	for (int i=0; i<neighbors[ti].size(); i++){
     		for (int j=i+1; j<neighbors[ti].size(); j++){

src/pyscal3/core.py

@@ -392,41 +392,30 @@ def atoms(self, atoms):
         """
         Set atoms
         """
-        #mod = False
+        #first task here is to find the size of the box needed
+        if np.sum(self.box) == 0:
+            raise ValueError("Simulation box should be initialized before filling with atoms")
+
+        box = copy.deepcopy(self.box)
+
         if(len(atoms['positions']) < 200):
-            #we need to estimate a rough idea
             needed_atoms = 200 - len(atoms['positions'])
-            #get a rough cell
-            #print(needed_atoms)
             needed_cells = np.ceil(needed_atoms/len(atoms['positions']))
-            nx = int(needed_cells**(1/3))
-            if nx < 2:
-                nx = 2
-            if np.sum(self.box) == 0:
-                raise ValueError("Simulation box should be initialized before atoms")
-            atoms, box = operations.repeat(self, (nx, nx, nx), atoms=atoms, ghost=True, return_atoms=True)
-            self.actual_box = self.box.copy()
-            self.internal_box = box
-            #mod = True
-
-        #we also check if each dimension is small
-        box = self.internal_box.copy()
-        repeats = [1,1,1]
-        for count, side in enumerate(box):
-            val = np.linalg.norm(side)
-            if val < 10:
-                repeats[count] = int(np.ceil(10/val))
-
-        #if any side is greater run
-        prod = np.prod(np.array(repeats))
-        if prod > 1:
-            #we need to scale sides
-            atoms, box = operations.repeat(self, repeats, atoms=atoms, 
-                ghost=True, return_atoms=True)
+            nx = max(int(needed_cells**(1/3)), 2)
+            scaled_box = operations.get_scaled_box(box, (nx, nx, nx))
+            #now check if this box is enough
+            repeats = [1,1,1]
+            for count, side in enumerate(scaled_box):
+                val = np.linalg.norm(side)
+                if val < 10:
+                    repeats[count] = int(np.ceil(10/val))
+            #now finalize repeats
+            repeats = np.array(repeats)*nx
+            #and finally scale
+            atoms, box = operations.repeat(self, repeats, atoms=atoms, ghost=True, return_atoms=True)            
             self.actual_box = self.box.copy()
             self.internal_box = box
 
-
         self._atoms = atoms
 
 

src/pyscal3/operations/operations.py

@@ -181,6 +181,12 @@ def repeat(system, repetitions,
         system.atoms = atoms
         return system
 
+def get_scaled_box(box, repetitions):
+    box = copy.deepcopy(box)
+    box[0] = repetitions[0]*np.array(box[0])
+    box[1] = repetitions[1]*np.array(box[1])
+    box[2] = repetitions[2]*np.array(box[2])
+    return box
 
 def embed_in_cubic_box(system, input_box=None, 
     return_box=False):

tests/test_centro.py

@@ -2,6 +2,7 @@
 import os,sys,inspect
 import numpy as np
 import pyscal3.core as pc
+from ase.build import bulk
 
 def test_cs_ges():
     sys = pc.System.create.lattice.bcc(repetitions = [7, 7, 7], lattice_constant=4.00)
@@ -10,3 +11,11 @@ def test_cs_ges():
 
     q = sys.calculate.centrosymmetry(nmax=4)
     assert np.round(np.mean(np.array(q)), decimals=2) > 0.00
+
+def test_hcp():
+    ti_hcp = bulk("Ti", orthorhombic=True)
+    sys = pc.System(ti_hcp, format='ase')
+    q = sys.calculate.centrosymmetry(nmax=12)
+    assert np.round(np.mean(np.array(q)), decimals=2) == 8.7
+
+