From d8d8ab024eaed3e6358a3db9e428ee008fd33c51 Mon Sep 17 00:00:00 2001 From: Diego Prada Date: Mon, 23 Dec 2024 15:45:06 -0600 Subject: [PATCH] in process --- molsysmt/lib/pbc/__init__.py | 5 +- molsysmt/lib/pbc/wrap_to_pbc.py | 85 +++++++++++++++++++-------------- molsysmt/pbc/wrap_to_pbc.py | 29 +++++++---- 3 files changed, 73 insertions(+), 46 deletions(-) diff --git a/molsysmt/lib/pbc/__init__.py b/molsysmt/lib/pbc/__init__.py index 7c47b44bf..ad8f3cae4 100644 --- a/molsysmt/lib/pbc/__init__.py +++ b/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 diff --git a/molsysmt/lib/pbc/wrap_to_pbc.py b/molsysmt/lib/pbc/wrap_to_pbc.py index 689e1c19c..750f8dede 100644 --- a/molsysmt/lib/pbc/wrap_to_pbc.py +++ b/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,9 +35,9 @@ 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] @@ -45,39 +45,37 @@ def wrap_to_pbc_vector_single_structure(vector, box, inv_box, orthogonal): 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 diff --git a/molsysmt/pbc/wrap_to_pbc.py b/molsysmt/pbc/wrap_to_pbc.py index cf8f90115..33aba2d71 100644 --- a/molsysmt/pbc/wrap_to_pbc.py +++ b/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: