Commit 90e99a2b authored by Remi Perenon's avatar Remi Perenon
Browse files

Modification of source file to compel with new scientific and MMTK ILL version

parent 9cf2fc3f
...@@ -103,7 +103,7 @@ def mic_generator_2D(double[:,:] pts not None, double border, double[:] box_para ...@@ -103,7 +103,7 @@ def mic_generator_2D(double[:,:] pts not None, double border, double[:] box_para
def mic_generator_3D(double[:,:] pts not None, double border, double[:] box_param): def mic_generator_3D(double[:,:] pts not None, double border, double[:] box_param):
cdef int i, m, new_id, nb_init_points cdef int i, j_int, m, new_id, nb_init_points
cdef double x, y, z, Xmin, Xmax, Ymin, Ymax, Zmin, Zmax, j, k, l, n0, n1, n2 cdef double x, y, z, Xmin, Xmax, Ymin, Ymax, Zmin, Zmax, j, k, l, n0, n1, n2
cdef vector[equiv_t] equivalence cdef vector[equiv_t] equivalence
cdef vector[point] new_points cdef vector[point] new_points
...@@ -168,8 +168,9 @@ def mic_generator_3D(double[:,:] pts not None, double border, double[:] box_para ...@@ -168,8 +168,9 @@ def mic_generator_3D(double[:,:] pts not None, double border, double[:] box_para
for i in range(new_points.size()): for i in range(new_points.size()):
j = i + nb_init_points j = i + nb_init_points
res[j,0] = new_points[i][0] j_int = int(j)
res[j,1] = new_points[i][1] res[j_int,0] = new_points[i][0]
res[j,2] = new_points[i][2] res[j_int,1] = new_points[i][1]
res[j_int,2] = new_points[i][2]
return res, d return res, d
...@@ -13,8 +13,8 @@ def mt(ndarray[np.float64_t, ndim = 2] config not None, ...@@ -13,8 +13,8 @@ def mt(ndarray[np.float64_t, ndim = 2] config not None,
ndarray[np.int32_t, ndim = 3] grid not None, ndarray[np.int32_t, ndim = 3] grid not None,
double resolution,ndarray[np.float64_t, ndim = 1] mini not None): double resolution,ndarray[np.float64_t, ndim = 1] mini not None):
cdef int at, nbatom , atom cdef int at, nbatom, atom, i, j, k
cdef double Xpos, Ypos, Zpos, i, j, k , mx, my, mz cdef double Xpos, Ypos, Zpos, mx, my, mz
mx = mini[0] mx = mini[0]
my = mini[1] my = mini[1]
...@@ -23,7 +23,7 @@ def mt(ndarray[np.float64_t, ndim = 2] config not None, ...@@ -23,7 +23,7 @@ def mt(ndarray[np.float64_t, ndim = 2] config not None,
nbatom = config.shape[0] nbatom = config.shape[0]
for atom in range(nbatom): for atom in range(nbatom):
# The of atoms |i| in the current configuration. # The of atoms |i| in the current configuration.
i = floor((config[atom,0]-mx)/resolution) i = int(floor((config[atom,0]-mx)/resolution))
j = floor((config[atom,1]-my)/resolution) j = int(floor((config[atom,1]-my)/resolution))
k = floor((config[atom,2]-mz)/resolution) k = int(floor((config[atom,2]-mz)/resolution))
grid[i,j,k] += 1 grid[i][j][k] += 1
...@@ -96,8 +96,12 @@ class DipoleAutoCorrelationFunction(IJob): ...@@ -96,8 +96,12 @@ class DipoleAutoCorrelationFunction(IJob):
conf = self.configuration['trajectory']['instance'].universe.contiguousObjectConfiguration() conf = self.configuration['trajectory']['instance'].universe.contiguousObjectConfiguration()
dipoleMoment = numpy.zeros((3,),dtype=numpy.float64) dipoleMoment = numpy.zeros((3,),dtype=numpy.float64)
for idx in self._indexes: for idx in self._indexes:
dipoleMoment += self.configuration["atom_charges"]["charges"][idx]*conf[idx] temp = self.configuration["atom_charges"]["charges"][idx]*conf[idx]
# Loop to sum temp because MMTK redefines __add__ operator, leading to a crash
#dipoleMoment = dipoleMoment + self.configuration["atom_charges"]["charges"][idx]*conf[idx]
for k in range(len(temp)):
dipoleMoment[k] = dipoleMoment[k] + temp[k]
return index, dipoleMoment return index, dipoleMoment
......
...@@ -107,12 +107,12 @@ class MolecularTrace(IJob): ...@@ -107,12 +107,12 @@ class MolecularTrace(IJob):
self.min = numpy.array([minx, miny, minz], dtype = numpy.float64) self.min = numpy.array([minx, miny, minz], dtype = numpy.float64)
self._outputData.add('origin',"line", self.min, units = 'nm') self._outputData.add('origin',"line", self.min, units = 'nm')
self.gdim = numpy.ceil(numpy.array([dimx, dimy, dimz])/self.resolution) self.gdim = numpy.ceil(numpy.array([dimx, dimy, dimz])/self.resolution).astype(numpy.int)
spacing = self.configuration['spatial_resolution']['value'] spacing = self.configuration['spatial_resolution']['value']
self._outputData.add('spacing',"line",numpy.array([spacing, spacing, spacing]), units = 'nm') self._outputData.add('spacing',"line",numpy.array([spacing, spacing, spacing]), units = 'nm')
self.grid = numpy.zeros(self.gdim, dtype = numpy.int32) self.grid = numpy.zeros(self.gdim, dtype = numpy.int32)
self._outputData.add('molecular_trace',"volume", tuple(numpy.ceil(numpy.array([dimx, dimy, dimz])/self.resolution))) self._outputData.add('molecular_trace',"volume", tuple(numpy.ceil(numpy.array([dimx, dimy, dimz])/self.resolution).astype(numpy.int)))
self._indexes = [idx for idxs in self.configuration['atom_selection']['indexes'] for idx in idxs] self._indexes = [idx for idxs in self.configuration['atom_selection']['indexes'] for idx in idxs]
......
#MDANSE : Molecular Dynamics Analysis for Neutron Scattering Experiments #MDANSE : Molecular Dynamics Analysis for Neutron Scattering Experiments
#------------------------------------------------------------------------------------------ #------------------------------------------------------------------------------------------
#Copyright (C) #Copyright (C)
#2015- Eric C. Pellegrini Institut Laue-Langevin #2015- Eric C. Pellegrini Institut Laue-Langevin
#BP 156 #BP 156
#6, rue Jules Horowitz #6, rue Jules Horowitz
#38042 Grenoble Cedex 9 #38042 Grenoble Cedex 9
#France #France
#pellegrini[at]ill.fr #pellegrini[at]ill.fr
#goret[at]ill.fr #goret[at]ill.fr
#aoun[at]ill.fr #aoun[at]ill.fr
# #
#This library is free software; you can redistribute it and/or #This library is free software; you can redistribute it and/or
#modify it under the terms of the GNU Lesser General Public #modify it under the terms of the GNU Lesser General Public
#License as published by the Free Software Foundation; either #License as published by the Free Software Foundation; either
#version 2.1 of the License, or (at your option) any later version. #version 2.1 of the License, or (at your option) any later version.
# #
#This library is distributed in the hope that it will be useful, #This library is distributed in the hope that it will be useful,
#but WITHOUT ANY WARRANTY; without even the implied warranty of #but WITHOUT ANY WARRANTY; without even the implied warranty of
#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU #MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
#Lesser General Public License for more details. #Lesser General Public License for more details.
# #
#You should have received a copy of the GNU Lesser General Public #You should have received a copy of the GNU Lesser General Public
#License along with this library; if not, write to the Free Software #License along with this library; if not, write to the Free Software
#Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA #Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA
''' '''
Created on Mar 30, 2015 Created on Mar 30, 2015
:author: Eric C. Pellegrini :author: Eric C. Pellegrini
''' '''
import collections import collections
import random import random
import numpy import numpy
from MDANSE import REGISTRY from MDANSE import REGISTRY
from MDANSE.Framework.QVectors.LatticeQvectors import LatticeQVectors from MDANSE.Framework.QVectors.LatticeQvectors import LatticeQVectors
class SphericalLatticeQVectors(LatticeQVectors): class SphericalLatticeQVectors(LatticeQVectors):
""" """
""" """
settings = collections.OrderedDict() settings = collections.OrderedDict()
settings['seed'] = ('integer', {"mini":0, "default":0}) settings['seed'] = ('integer', {"mini":0, "default":0})
settings['shells'] = ('range', {"valueType":float, "includeLast":True, "mini":0.0}) settings['shells'] = ('range', {"valueType":float, "includeLast":True, "mini":0.0})
settings['n_vectors'] = ('integer', {"mini":1, "default":50}) settings['n_vectors'] = ('integer', {"mini":1, "default":50})
settings['width'] = ('float', {"mini":1.0e-6, "default":1.0}) settings['width'] = ('float', {"mini":1.0e-6, "default":1.0})
def _generate(self): def _generate(self):
if self._configuration["seed"]["value"] != 0: if self._configuration["seed"]["value"] != 0:
numpy.random.seed(self._configuration["seed"]["value"]) numpy.random.seed(self._configuration["seed"]["value"])
random.seed(self._configuration["seed"]["value"]) random.seed(self._configuration["seed"]["value"])
qMax = self._configuration["shells"]["last"] + 0.5*self._configuration["width"]["value"] qMax = self._configuration["shells"]["last"] + 0.5*self._configuration["width"]["value"]
hklMax = numpy.ceil([qMax/v.length()for v in self._reciprocalBasis]) + 1 hklMax = numpy.ceil([qMax/v.length()for v in self._reciprocalBasis]) + 1
vects = numpy.mgrid[-hklMax[0]:hklMax[0]+1,-hklMax[1]:hklMax[1]+1,-hklMax[2]:hklMax[2]+1] vects = numpy.mgrid[-hklMax[0]:hklMax[0]+1,-hklMax[1]:hklMax[1]+1,-hklMax[2]:hklMax[2]+1]
vects = vects.reshape(3,(2*hklMax[0]+1)*(2*hklMax[1]+1)*(2*hklMax[2]+1)) vects = vects.reshape(3,int(2*hklMax[0]+1)*int(2*hklMax[1]+1)*int(2*hklMax[2]+1))
vects = numpy.dot(self._reciprocalMatrix,vects) vects = numpy.dot(self._reciprocalMatrix,vects)
dists2 = numpy.sum(vects**2,axis=0) dists2 = numpy.sum(vects**2,axis=0)
halfWidth = self._configuration["width"]["value"]/2 halfWidth = self._configuration["width"]["value"]/2
nVectors = self._configuration["n_vectors"]["value"] nVectors = self._configuration["n_vectors"]["value"]
if self._status is not None: if self._status is not None:
self._status.start(self._configuration["shells"]['number']) self._status.start(self._configuration["shells"]['number'])
self._configuration["q_vectors"] = collections.OrderedDict() self._configuration["q_vectors"] = collections.OrderedDict()
for q in self._configuration["shells"]["value"]: for q in self._configuration["shells"]["value"]:
qmin = max(0,q - halfWidth) qmin = max(0,q - halfWidth)
q2low = qmin*qmin q2low = qmin*qmin
q2up = (q + halfWidth)*(q + halfWidth) q2up = (q + halfWidth)*(q + halfWidth)
hits = numpy.where((dists2 >= q2low) & (dists2 <= q2up))[0] hits = numpy.where((dists2 >= q2low) & (dists2 <= q2up))[0]
nHits = len(hits) nHits = len(hits)
if nHits != 0: if nHits != 0:
n = min(nHits,nVectors) n = min(nHits,nVectors)
if nHits > nVectors: if nHits > nVectors:
hits = random.sample(hits,nVectors) hits = random.sample(hits,nVectors)
self._configuration["q_vectors"][q] = {} self._configuration["q_vectors"][q] = {}
self._configuration["q_vectors"][q]['q_vectors'] = vects[:,hits] self._configuration["q_vectors"][q]['q_vectors'] = vects[:,hits]
self._configuration["q_vectors"][q]['n_q_vectors'] = n self._configuration["q_vectors"][q]['n_q_vectors'] = n
self._configuration["q_vectors"][q]['q'] = q self._configuration["q_vectors"][q]['q'] = q
self._configuration["q_vectors"][q]['hkls'] = numpy.rint(numpy.dot(self._invReciprocalMatrix, self._configuration["q_vectors"][q]['hkls'] = numpy.rint(numpy.dot(self._invReciprocalMatrix,
self._configuration["q_vectors"][q]['q_vectors'])) self._configuration["q_vectors"][q]['q_vectors']))
if self._status is not None: if self._status is not None:
if self._status.is_stopped(): if self._status.is_stopped():
return None return None
else: else:
self._status.update() self._status.update()
REGISTRY["spherical_lattice"] = SphericalLatticeQVectors REGISTRY["spherical_lattice"] = SphericalLatticeQVectors
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