in process

molsysmt/lib/pbc/__init__.py

@@ -20,6 +20,7 @@
 from .wrap_to_mic import wrap_to_mic
 
 from .wrap_to_pbc import wrap_to_pbc_vector_single_structure
-from .wrap_to_pbc import wrap_to_pbc_single_structure
-from .wrap_to_pbc import wrap_to_pbc
+#from .wrap_to_pbc import wrap_to_pbc_single_structure
+from .wrap_to_pbc import wrap_to_pbc_center
+from .wrap_to_pbc import wrap_to_pbc_no_center
 

molsysmt/lib/pbc/wrap_to_pbc.py

@@ -24,9 +24,9 @@ def wrap_to_pbc_vector_single_structure(vector, box, inv_box, orthogonal):
 
     if orthogonal:
 
-        output[0]=vector[0]-box[0,0]*round(vector[0]/box[0,0])
-        output[1]=vector[1]-box[1,1]*round(vector[1]/box[1,1])
-        output[2]=vector[2]-box[2,2]*round(vector[2]/box[2,2])
+        output[0]=vector[0]-box[0,0]*(vector[0]//box[0,0])
+        output[1]=vector[1]-box[1,1]*(vector[1]//box[1,1])
+        output[2]=vector[2]-box[2,2]*(vector[2]//box[2,2])
 
     else:
 
@@ -35,49 +35,47 @@ def wrap_to_pbc_vector_single_structure(vector, box, inv_box, orthogonal):
         vaux[0]=inv_box[0,0]*vector[0]+inv_box[1,0]*vector[1]+inv_box[2,0]*vector[2]
         vaux[1]=                       inv_box[1,1]*vector[1]+inv_box[2,1]*vector[2]
         vaux[2]=                                              inv_box[2,2]*vector[2]
-        vaux[0]=vaux[0]-round(vaux[0])
-        vaux[1]=vaux[1]-round(vaux[1])
-        vaux[2]=vaux[2]-round(vaux[2])
+        vaux[0]=vaux[0]-(vaux[0]//1.0)
+        vaux[1]=vaux[1]-(vaux[1]//1.0)
+        vaux[2]=vaux[2]-(vaux[2]//1.0)
         output[0]=box[0,0]*vaux[0]+box[1,0]*vaux[1]+box[2,0]*vaux[2]
         output[1]=                 box[1,1]*vaux[1]+box[2,1]*vaux[2]
         output[2]=                                  box[2,2]*vaux[2]
 
     return output
 
 
-arguments=[nb.float64[:,:], # coordinates
-           nb.float64[:,:], # box
-           nb.float64[:], # center
-           nb.boolean, # center_at_origin
-          ]
-output=None
-@nb.njit(make_numba_signature(arguments, output), cache=True)
-def wrap_to_pbc_single_structure(coordinates, box, center, center_at_origin):
-
-    n_atoms = coordinates.shape[0]
-
-    orthogonal = box_is_orthogonal_single_structure(box[:,:])
-    inv_box = inverse_matrix_3x3(box)
-
-    for ii in range(n_atoms):
-        tmp_vect = coordinates[ii,:]-center
-        tmp_vect = wrap_to_pbc_vector_single_structure(tmp_vect, box, inv_box,
-                orthogonal)
-        if not center_at_origin:
-            tmp_vect=tmp_vect+center
-        coordinates[ii,:]=tmp_vect
-
-    pass
-
+#arguments=[nb.float64[:,:], # coordinates
+#           nb.float64[:,:], # box
+#           nb.float64[:], # center
+#           nb.boolean, # center_at_origin
+#          ]
+#output=None
+#@nb.njit(make_numba_signature(arguments, output), cache=True)
+#def wrap_to_pbc_single_structure(coordinates, box, center, center_at_origin):
+#
+#    n_atoms = coordinates.shape[0]
+#
+#    orthogonal = box_is_orthogonal_single_structure(box[:,:])
+#    inv_box = inverse_matrix_3x3(box)
+#
+#    for ii in range(n_atoms):
+#        tmp_vect = coordinates[ii,:]-center
+#        tmp_vect = wrap_to_pbc_vector_single_structure(tmp_vect, box, inv_box,
+#                orthogonal)
+#        if not center_at_origin:
+#            tmp_vect=tmp_vect+center
+#        coordinates[ii,:]=tmp_vect
+#
+#    pass
 
 arguments=[nb.float64[:,:,:], # coordinates
            nb.float64[:,:,:], # box
-           nb.float64[:,:,:], # center
-           nb.boolean, # center_at_origin
+           nb.float64[:,:,:] # center
           ]
 output=None
 @nb.njit(make_numba_signature(arguments, output), cache=True)
-def wrap_to_pbc(coordinates, box, center, center_at_origin):
+def wrap_to_pbc_center(coordinates, box, center):
 
     n_structures, n_atoms = coordinates.shape[:2]
 
@@ -92,12 +90,29 @@ def wrap_to_pbc(coordinates, box, center, center_at_origin):
         for jj in range(n_atoms):
             tmp_vect = coordinates[ii,jj,:]-tmp_center
             tmp_vect = wrap_to_pbc_vector_single_structure(tmp_vect, tmp_box, inv_box, orthogonal)
-            if not center_at_origin:
-                tmp_vect=tmp_vect+tmp_center
             coordinates[ii,jj,:]=tmp_vect
         if not single_structure_center:
             aa+=1
 
     pass
 
+arguments=[nb.float64[:,:,:], # coordinates
+           nb.float64[:,:,:] # box
+          ]
+output=None
+@nb.njit(make_numba_signature(arguments, output), cache=True)
+def wrap_to_pbc_no_center(coordinates, box):
+
+    n_structures, n_atoms = coordinates.shape[:2]
+
+    for ii in range(n_structures):
+        tmp_box = box[ii,:,:]
+        orthogonal = box_is_orthogonal_single_structure(tmp_box)
+        inv_box = inverse_matrix_3x3(tmp_box)
+        for jj in range(n_atoms):
+            tmp_vect = coordinates[ii,jj,:]
+            tmp_vect = wrap_to_pbc_vector_single_structure(tmp_vect, tmp_box, inv_box, orthogonal)
+            coordinates[ii,jj,:]=tmp_vect
+
+    pass
 

molsysmt/pbc/wrap_to_pbc.py

@@ -7,8 +7,8 @@
 
 @digest()
 def wrap_to_pbc(molecular_system, selection='all', structure_indices='all',
-                center_coordinates='[0,0,0] nanometers', center_of_selection=None, weights=None,
-                center_at_origin=True, keep_covalent_bonds=False,
+                center_at_origin=False, center_of_selection='all', weights=None,
+                center_coordinates='[0,0,0] nanometers', keep_covalent_bonds=False,
                 syntax='MolSysMT', engine='MolSysMT', in_place=False):
     """
     To be written soon...
@@ -24,21 +24,32 @@ def wrap_to_pbc(molecular_system, selection='all', structure_indices='all',
         coordinates= get(molecular_system, element='atom', selection=atom_indices, coordinates=True)
         box = get(molecular_system, element='system', structure_indices=structure_indices, box=True)
 
-        if center_of_selection is not None:
+        if center_at_origin:
 
             center_coordinates = get_center(molecular_system, selection=center_of_selection,
                                 weights=weights, structure_indices=structure_indices,
                                 syntax=syntax, engine='MolSysMT')
 
-        coordinates, length_units = puw.get_value_and_unit(coordinates)
-        box = puw.get_value(box, to_unit=length_units)
-        center_coordinates = puw.get_value(center_coordinates, to_unit=length_units)
+            coordinates, length_units = puw.get_value_and_unit(coordinates)
+            box = puw.get_value(box, to_unit=length_units)
+            center_coordinates = puw.get_value(center_coordinates, to_unit=length_units)
 
-        msmlib.pbc.wrap_to_pbc(coordinates, box, center_coordinates, center_at_origin)
+            msmlib.pbc.wrap_to_pbc_center(coordinates, box, center_coordinates)
 
-        coordinates=puw.quantity(coordinates, length_units)
+            coordinates=puw.quantity(coordinates, length_units)
 
-        del(box, center_coordinates)
+            del(box, center_coordinates)
+
+        else:
+
+            coordinates, length_units = puw.get_value_and_unit(coordinates)
+            box = puw.get_value(box, to_unit=length_units)
+
+            msmlib.pbc.wrap_to_pbc_no_center(coordinates, box)
+
+            coordinates=puw.quantity(coordinates, length_units)
+
+            del(box)
 
     else: