Entropy

.gitignore

+/test/
\ No newline at end of file

EntropyCalculator.py

+import os
+import ase.io
+import numpy as np
+from scipy import signal, fftpack
+from scipy import optimize 
+from matplotlib import pyplot as plt
+from scipy.integrate import simps
+from scipy.special import wofz
+import copy
+from ase.data import atomic_numbers, chemical_symbols, atomic_masses
+
+AA = 1 * 1e-10
+h = 6.626e-34
+k = 1.3806505 * 1e-23
+eV =  1.602176565*1e-19
+fs = 1e-15
+mm = 1.667 * 1e-27
+
+def dichotomy(f,left,right,eps,maxN=20000):
+    middle = (left+right)/2
+    count=0
+    while right-left>eps or count<maxN:
+        middle = (left+right)/2
+        if f(left)*f(middle)<=0:
+            right=middle
+        else:
+            left=middle
+        count=count+1
+    return middle
+
+class Trajectory:
+    def __init__(self,frames,dt):
+        self.frames=frames
+        self.dt=dt
+        self.Nsteps=len(frames)
+        self.Natoms=len(frames[0])
+        self.pbccoords=self.getPBCpositions()
+        self.velocities=self.getVelocities()
+        self.vacf=self.getVACF()
+        self.meanvolume=np.mean(np.array([a.get_volume() for a in frames]))
+        
+    def getVelocities(self):
+        """
+        velocities=(self.pbccoords[1:]-self.pbccoords[:-1])*AA/self.dt
+        return velocities
+        """
+        if self.frames[0].get_velocities() is None:
+            velocities=(self.pbccoords[1:]-self.pbccoords[:-1])*AA/self.dt
+        else:
+            velocities=np.zeros([self.Nsteps,self.Natoms,3])
+            for i,a in enumerate(frames):
+                velocities[i]=a.get_velocities()*100
+        return velocities
+            
+    def getPBCpositions(self):
+        coords=np.zeros([self.Nsteps,self.Natoms,3])
+        coords[0]=self.frames[0].get_positions()
+        pold = self.frames[0].get_scaled_positions()
+        for i in range(1,self.Nsteps):
+            pnew = self.frames[i].get_scaled_positions() 
+            pnew = pnew - np.round(pnew - pold)
+            coords[i]=np.dot(pnew,self.frames[i].get_cell())
+            pold=pnew
+        return coords
+    
+    def getVACF_(velocities):
+        vacf = np.zeros_like(velocities)
+        for i in range(velocities.shape[1]):
+            vel_tmp = velocities[:, i] 
+            vacf[:, i] = signal.fftconvolve(vel_tmp, vel_tmp[::-1], mode = "full", axes = 0)[len(vel_tmp)-1:]
+            vacf[:, i] /= len(vel_tmp)
+        return vacf
+    
+    def getVACF(self):
+        vacf = np.zeros([self.velocities.shape[0],self.velocities.shape[1]])
+        for i in range(self.velocities.shape[1]):
+            vel_tmp = self.velocities[:, i] 
+            vacf[:, i] = np.sum(signal.fftconvolve(vel_tmp, vel_tmp[::-1], mode = "full", axes = 0)[len(vel_tmp)-1:],axis=1)
+            vacf[:, i] /= len(vel_tmp)*3
+        return vacf
+    
+    def getAverageStructure(self,b,e):
+        average_structure = copy.deepcopy(self.frames[0])
+        average_structure.set_positions(np.sum(self.pbccoords[b:e],axis=0)/(b-e+1))
+        return average_structure
+        
+    def getMSD_t(self,timerange=None,atomsgroup=None):
+        if timerange is None:
+            timerange=np.arange(self.Nsteps)
+        if atomsgroup is None:
+            atomsgroup=np.arange(self.Natoms)
+        msd=(self.pbccoords[timerange]-self.pbccoords[timerange[0]])**2
+        return np.mean(msd[:,atomsgroup],axis=1)
+    
+    def getMSD_atom(self,timerange=None,atomsgroup=None):
+        if timerange is None:
+            timerange=np.arange(self.Nsteps)
+        if atomsgroup is None:
+            atomsgroup=np.arange(self.Natoms)
+        msd=(self.pbccoords[timerange]-self.pbccoords[timerange[0]])**2
+        return np.mean(msd[:,atomsgroup],axis=0)
+    
+    
+class EntrophyCalculator:
+    def __init__(self,vacf,m,N,T,V,dt,phase='liquid',frequency_range = np.linspace(1e-13,40,500)):
+        self.vacf=vacf
+        self.m=m* mm
+        self.N=N
+        self.T=T
+        self.V=V
+        self.dt=dt
+        if phase is None:
+            
+        else:
+            self.phase=phase
+        self.frequency_range = frequency_range
+        self.division=len(frequency_range)
+        self.F=self.GetDoS()
+        self.D=simps(self.vacf,dx=dt)
+        #self.solid_entropy,self.gas_entropy,self.liquid_entropy=self.GetEntrophy()
+        
+    def GetDoS(self):
+        F = np.zeros_like(self.frequency_range)
+        for i in range(self.division):
+            f=np.cos(2*np.pi*self.frequency_range[i]*1e12*np.arange(0,len(self.vacf)*self.dt,self.dt))*self.vacf
+            F[i] = simps(f, dx = self.dt)
+        F*=12*self.m/(k*self.T)*1e12
+        return F
+    
+    def Solve_Gamma(self):
+        def solve_function(gamma):
+            return 2*((1-gamma)**3)/(2-gamma)/(gamma**0.4)-self.delta**0.6
+        return dichotomy(solve_function,1e-5,1,1e-6)
+    
+    def GetfgFg(self):
+        def SolvefgAgBg(fg0,alpha,T,D,M2,M4):
+            def solveAgBg(fg):
+                A=fg*k*T/(self.m)/D
+                ans=(np.sqrt(alpha**2/np.pi)+np.sqrt(A**2/np.pi+alpha**2/4-A**2/4))/(alpha/A-A/alpha)
+                return A*ans*2/np.sqrt(np.pi),ans**2
+            
+            def solve_function(fg):
+                Ag,Bg=solveAgBg(fg)
+                return Bg-(fg*(2*M2*Ag-M4-Ag**2)+M4-M2**2)/(2*fg*(1-fg)*Ag)
+            
+            fg=dichotomy(solve_function,fg0*0.5,fg0*0.99,1e-5)
+            Ag,Bg=solveAgBg(fg)
+            return fg,Ag,Bg
+    
+        def Function_Fg(nu,Ag,Bg):
+            K0 = Ag * np.sqrt(np.pi / Bg) / 2.
+            eta = np.pi * nu / np.sqrt(Bg)
+            Kg = K0 * wofz(-eta)
+            num = Kg + np.conj(Kg)
+            Kg += complex(0.,2. * np.pi * nu)
+            denom = Kg * np.conj(Kg)
+            Fg = 6. * num / denom
+            return Fg
+        
+        M2=simps(1/3*((2*np.pi*self.frequency_range*1e12)**2)*self.F,self.frequency_range)
+        M4=simps(1/3*((2*np.pi*self.frequency_range*1e12)**4)*self.F,self.frequency_range)
+        fg0=(self.gamma ** 0.4) * (self.delta ** 0.6)
+        fg,Ag,Bg = SolvefgAgBg(fg0,self.alpha,self.T,self.D,M2,M4)
+    
+        fgFg = []
+        for i in frequency_range:
+            fgFg.append(fg * Function_Fg(i * 1e12,Ag,Bg).real)
+        fgFg = np.array(fgFg)*1e12
+        self.fg=fg
+        return fgFg
+    
+    def CalculateEntrophy(self):
+        if self.phase=='liquid':
+            self.delta = 8/3*((6/np.pi)**(2/3))*self.D*np.sqrt(np.pi*self.m/(k*self.T))*((self.N/self.V)**(1/3)) 
+            self.gamma = self.Solve_Gamma()
+            self.alpha=(k*self.T)/(self.m*self.D)*(self.gamma**0.4)*(self.delta**0.6)
+            self.fgFg=self.GetfgFg()
+            self.gas_entropy =  self.GetGasEntrophy()
+            self.solid_entropy = self.GetSolidEntrophy(self.F-self.fgFg)
+            self.entropy = self.solid_entropy + self.gas_entropy
+        
+        if self.phase=='solid':
+            self.gas_entropy =  0
+            self.solid_entropy = self.GetSolidEntrophy(self.F)
+            self.entropy = self.solid_entropy
+    
+    def GetGasEntrophy(self):
+        WIG=1/3*(5/2+np.log(((2*np.pi*self.m*k*self.T/h**2)**1.5)*self.V/(self.fg*self.N)))
+        Wx=1/3*(np.log((1+self.gamma+(self.gamma**2)-(self.gamma**3))/((1-self.gamma)**3)) \
+            +self.gamma*(3*self.gamma-4)/((1-self.gamma)** 2))
+        return 3 * self.fg * (WIG + Wx) * k * self.N
+    
+    def GetSolidEntrophy(self,Fs):
+        solid_distribution=((h*self.frequency_range*1e12)/(k*self.T))/\
+        (np.exp((h*self.frequency_range*1e12)/(k*self.T))-1) \
+        -np.log(1-np.exp(-((h*self.frequency_range*1e12)/(k* self.T))))
+        es=self.N*k*Fs*solid_distribution
+        return simps(es,self.frequency_range)
+           
+    def plotFandfgFg(self):
+        plt.figure()
+        plt.plot(self.F)
+        plt.plot(self.fgFg)
+
+if __name__=='__main__':
+    filetype='lammps-dump'
+    filetype='vasp-xdatcar'
+    dt= 0.5 * fs
+    frequency_range=np.linspace(1e-13,30,500)
+  
+    file='Ca'
+    frames = ase.io.read(file, index = ":", format = filetype)
+    traj=Trajectory(frames,dt) 
+    T=600.493
+    V=traj.meanvolume * AA**3
+    calculators=[]
+    for atomic_number in set(frames[0].get_atomic_numbers()):
+        symbol=chemical_symbols[atomic_number]
+        m=ase.data.atomic_masses[atomic_number]
+        vacf=np.mean(np.ndarray.squeeze(traj.vacf[:,np.where(frames[0].get_masses()==m)]),axis=1)
+        N=np.sum(frames[0].get_masses()==m)
+        c=EntrophyCalculator(vacf,m,N,T,V,dt,'liquid',frequency_range)
+        c.CalculateEntrophy()
+        print(symbol,c.entropy)
+        vacfs.append(vacf)
+        calculators.append(c)