CoordinationNumber.py 6.24 KB
Newer Older
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
#MDANSE : Molecular Dynamics Analysis for Neutron Scattering Experiments
#------------------------------------------------------------------------------------------
#Copyright (C)
#2015- Eric C. Pellegrini Institut Laue-Langevin
#BP 156
#6, rue Jules Horowitz
#38042 Grenoble Cedex 9
#France
#pellegrini[at]ill.fr
#goret[at]ill.fr
#aoun[at]ill.fr
#
#This library is free software; you can redistribute it and/or
#modify it under the terms of the GNU Lesser General Public
#License as published by the Free Software Foundation; either
#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,
#but WITHOUT ANY WARRANTY; without even the implied warranty of
#MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
#Lesser General Public License for more details.
#
#You should have received a copy of the GNU Lesser General Public
#License along with this library; if not, write to the Free Software
#Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301  USA

''' 
Created on Apr 10, 2015

30
:author: Eric C. Pellegrini
31
32
33
34
35
36
'''

import collections

import numpy

37
from MDANSE import REGISTRY
38
39
40
41
from MDANSE.Framework.Jobs.DistanceHistogram import DistanceHistogram

class CoordinationNumber(DistanceHistogram):
    """
42
43
    The Coordination Number is computed from the pair distribution function for a set of atoms. 
    It describes the total number of neighbours, as a function of distance, from a central atom, or the centre of a group of atoms.
44
45
46
47
    """

    label = "Coordination Number"
    
48
    category = ('Analysis','Structure',)
49
    
50
    ancestor = ["mmtk_trajectory","molecular_viewer"]        
51

52
53
54
55
56
    settings = collections.OrderedDict()
    settings['trajectory'] = ('mmtk_trajectory',{})
    settings['frames'] = ('frames', {'dependencies':{'trajectory':'trajectory'}})
    settings['r_values'] = ('range', {'valueType':float, 'includeLast':True, 'mini':0.0})
    settings['atom_selection'] = ('atom_selection', {'dependencies':{'trajectory':'trajectory'}})
57
    settings['atom_transmutation'] = ('atom_transmutation', {'dependencies':{'trajectory':'trajectory','atom_selection':'atom_selection'}})
58
59
    settings['output_files'] = ('output_files', {'formats':["netcdf","ascii"]})
    settings['running_mode'] = ('running_mode',{})
60
61
62
63
64
65
66
67
                
    def finalize(self):
        """
        Finalizes the calculations (e.g. averaging the total term, output files creations ...).
        """   

        npoints = len(self.configuration['r_values']['mid_points'])

68
        self._outputData.add('r','line', self.configuration['r_values']['mid_points'], units="nm") 
69
70
71
        
        for pair in self._elementsPairs:
            invPair = pair[::-1]
72
73
74
75
76
77
            self._outputData.add("cn_intra_%s%s" % pair,"line", (npoints,), axis='r',units="au")                                                 
            self._outputData.add("cn_inter_%s%s" % pair,"line", (npoints,), axis='r', units="au",)                                                 
            self._outputData.add("cn_total_%s%s" % pair,"line", (npoints,), axis='r', units="au")                                                 
            self._outputData.add("cn_intra_%s%s" % invPair,"line", (npoints,), axis='r', units="au")                                                 
            self._outputData.add("cn_inter_%s%s" % invPair,"line", (npoints,), axis='r', units="au",)                                                 
            self._outputData.add("cn_total_%s%s" % invPair,"line", (npoints,), axis='r', units="au")                                                 
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
        
        nFrames = self.configuration['frames']['number']
        
        densityFactor = 4.0*numpy.pi*self.configuration['r_values']['mid_points']
        
        shellSurfaces = densityFactor*self.configuration['r_values']['mid_points']
        
        shellVolumes  = shellSurfaces*self.configuration['r_values']['step']
  
        self.averageDensity *= 4.0*numpy.pi/nFrames

        r2 = self.configuration['r_values']['mid_points']**2
        dr = self.configuration['r_values']['step']
        
        for k in self._concentrations.keys():
            self._concentrations[k] /= nFrames

95
        nAtomsPerElement = self.configuration['atom_selection'].get_natoms()
96
97
        for pair in self._elementsPairs:
            at1,at2 = pair 
98
99
            ni = nAtomsPerElement[at1]
            nj = nAtomsPerElement[at2]
100
            
101
102
            idi = self.selectedElements.index(at1)
            idj = self.selectedElements.index(at2)
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123

            if idi == idj:
                nij = ni*(ni-1)/2.0    
            else:
                nij = ni*nj
                self.hIntra[idi,idj] += self.hIntra[idj,idi]
                self.hInter[idi,idj] += self.hInter[idj,idi]
            
            fact = nij*nFrames*shellVolumes
             
            self.hIntra[idi,idj,:] /= fact             
            self.hInter[idi,idj,:] /= fact
                         
            cnIntra = numpy.add.accumulate(self.hIntra[idi,idj,:]*r2)*dr
            cnInter = numpy.add.accumulate(self.hInter[idi,idj,:]*r2)*dr
            cnTotal = cnIntra + cnInter
            
            cAlpha = self._concentrations[pair[0]]
            cBeta = self._concentrations[pair[1]]
                        
            invPair = pair[::-1]
eric pellegrini's avatar
eric pellegrini committed
124

125
126
127
128
129
130
131
            self._outputData["cn_intra_%s%s" % pair][:] = self.averageDensity*cBeta*cnIntra
            self._outputData["cn_inter_%s%s" % pair][:] = self.averageDensity*cBeta*cnInter        
            self._outputData["cn_total_%s%s" % pair][:] = self.averageDensity*cBeta*cnTotal 
            self._outputData["cn_intra_%s%s" % invPair][:] = self.averageDensity*cAlpha*cnIntra
            self._outputData["cn_inter_%s%s" % invPair][:] = self.averageDensity*cAlpha*cnInter        
            self._outputData["cn_total_%s%s" % invPair][:] = self.averageDensity*cAlpha*cnTotal 

132
        self._outputData.write(self.configuration['output_files']['root'], self.configuration['output_files']['formats'], self._info)
133
134
135
136
        
        self.configuration['trajectory']['instance'].close()     
  
        DistanceHistogram.finalize(self)
137
138

REGISTRY['cn'] = CoordinationNumber