Commit 344ab746 authored by Remi Perenon's avatar Remi Perenon Committed by Remi Perenon

Adding of new job for clustering universe according to van der waals radius

parent 2b5f9cae
Pipeline #4809 passed with stages
in 19 minutes and 6 seconds
......@@ -13,6 +13,7 @@
# **************************************************************************
import cython
import numpy as np
cimport numpy as np
from numpy cimport ndarray
......@@ -23,24 +24,77 @@ cdef extern from "math.h":
@cython.boundscheck(False)
@cython.wraparound(False)
def cpt_cluster_connectivity(ndarray[np.float64_t, ndim=2] coords,
ndarray[np.float64_t, ndim=1] covRadii,
np.float32_t tolerance,
bonds):
ndarray[np.float64_t, ndim=2] cell,
ndarray[np.float64_t, ndim=2] rcell,
ndarray[np.float64_t, ndim=1] covRadii,
np.float32_t scale,
np.float32_t offset,
bonds):
'''
Compute the connectivity of an atom cluster.
'''
cdef int i, j, nbat
cdef float radius, distance
cdef float radiusi, radiusj, radius, r2,
cdef double x, y, z, sdx, sdy, sdz, rx, ry, rz
nbat = len(covRadii)
cdef ndarray[np.float64_t, ndim=2] scaleconfig = np.zeros((nbat,3), dtype = np.float64)
if cell is None:
for i in range(nbat-1):
radiusi = covRadii[i]
for j in range(i+1,nbat):
radiusj = covRadii[j]
radius = scale*(radiusi + radiusj) + offset
r2 = (coords[i,0] - coords[j,0])**2 + (coords[i,1] - coords[j,1])**2 + (coords[i,2] - coords[j,2])**2
if r2 <= (radius*radius):
bonds.append([i,j])
else:
for 0 <= i < nbat:
x = coords[i,0]
y = coords[i,1]
z = coords[i,2]
scaleconfig[i,0] = x*rcell[0,0] + y*rcell[0,1] + z*rcell[0,2]
scaleconfig[i,1] = x*rcell[1,0] + y*rcell[1,1] + z*rcell[1,2]
scaleconfig[i,2] = x*rcell[2,0] + y*rcell[2,1] + z*rcell[2,2]
for 0 <= i < nbat:
sx = scaleconfig[i,0]
sy = scaleconfig[i,1]
sz = scaleconfig[i,2]
radiusi = covRadii[i]
for i + 1 <= j < nbat:
sdx = scaleconfig[j,0] - sx
sdy = scaleconfig[j,1] - sy
sdz = scaleconfig[j,2] - sz
sdx -= round(sdx)
sdy -= round(sdy)
sdz -= round(sdz)
rx = sdx*cell[0,0] + sdy*cell[0,1] + sdz*cell[0,2]
ry = sdx*cell[1,0] + sdy*cell[1,1] + sdz*cell[1,2]
rz = sdx*cell[2,0] + sdy*cell[2,1] + sdz*cell[2,2]
for i in range(nbat-1):
radiusi = covRadii[i] + tolerance
for j in range(i+1,nbat):
radius = radiusi + covRadii[j]
distance = (coords[i,0] - coords[j,0])**2 + (coords[i,1] - coords[j,1])**2 + (coords[i,2] - coords[j,2])**2
if distance <= (radius*radius):
bonds.append([i,j])
r2 = rx*rx + ry*ry + rz*rz
radiusj = covRadii[j]
radius = scale*(radiusi + radiusj) + offset
if r2 <= (radius*radius):
bonds.append([i,j])
return bonds
......@@ -41,21 +41,29 @@ def atomindex_to_moleculeindex(universe):
return d
def brute_formula(molecule, sep='_'):
def brute_formula(atom_list, sep='_'):
contents = {}
for at in molecule.atomList():
for at in atom_list:
contents[at.symbol] = str(int(contents.get(at.symbol,0)) + 1)
return sep.join([''.join(v) for v in sorted(contents.items())])
def build_connectivity(universe ,tolerance=0.05):
def build_connectivity(universe ,scale=1.0, tolerance=0.05):
bonds = []
conf = universe.contiguousObjectConfiguration()
try:
directCell = numpy.array(universe.basisVectors()).astype(numpy.float64)
reverseCell = numpy.array(universe.reciprocalBasisVectors()).astype(numpy.float64)
except AttributeError:
directCell = None
reverseCell = None
scannedObjects = [obj for obj in universe.objectList() if isinstance(obj,AtomCluster)]
singleAtomsObjects = [obj for obj in universe.objectList() if isinstance(obj,Atom) or obj.numberOfAtoms()==1]
if singleAtomsObjects:
......@@ -72,8 +80,8 @@ def build_connectivity(universe ,tolerance=0.05):
covRadii[i] = ELEMENTS[at.symbol.capitalize(),'covalent_radius']
bonds = []
fast_calculation.cpt_cluster_connectivity(coords,covRadii,tolerance,bonds)
fast_calculation.cpt_cluster_connectivity(coords,directCell,reverseCell,covRadii,scale,tolerance,bonds)
for idx1,idx2 in bonds:
if hasattr(atoms[idx1],"bonded_to__"):
atoms[idx1].bonded_to__.append(atoms[idx2])
......@@ -83,7 +91,7 @@ def build_connectivity(universe ,tolerance=0.05):
if hasattr(atoms[idx2],"bonded_to__"):
atoms[idx2].bonded_to__.append(atoms[idx1])
else:
atoms[idx2].bonded_to__ = [atoms[idx1]]
atoms[idx2].bonded_to__ = [atoms[idx1]]
def find_atoms_in_molecule(universe, moleculeName, atomNames, indexes=False):
......@@ -174,7 +182,7 @@ def resolve_undefined_molecules_name(universe):
for obj in universe.objectList():
if isChemicalObject(obj):
if not obj.name.strip():
obj.name = brute_formula(obj,sep="")
obj.name = brute_formula(obj.atomList(),sep="")
def sorted_atoms(universe,attribute=None):
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment