

In [1]:

    
from pylab import*
import glob, os
from os.path import join
from os.path import expanduser, join, split
from scipy import integrate

lammpsHeaderType = dtype([("step", 'int32'),("nAtoms", 'int32'),
                          ("bounds", float, (9,)),
                          ("nColumns", 'int32'),("nChunks", 'int32'),
                          ("chunkLength", 'int32')])

dataType = dtype([("id", float),("position", float, (3,)),
                  ("energy", float),("force", float),])

def loadAtoms(fileName):
    fileName = expanduser(fileName)
    lammpsFile = open(fileName, "rb")
    lammpsHeader = fromfile(lammpsFile, dtype=lammpsHeaderType, count=1)
    atoms = fromfile(lammpsFile, dtype=dataType)
    lammpsFile.close()
    return lammpsHeader, atoms



In [86]:

    
rawDataPath = "/home/milad/kurs/qmd"
stateFiles = glob.glob1(rawDataPath,'*.lmp')
nStateFiles = len(stateFiles)

data = []
for i in range(0, nStateFiles):
    stateFile = join(rawDataPath,"state" + "%04d" % i  + ".lmp")
    if os.path.exists(stateFile):
        header, atoms =loadAtoms(stateFile)
        data.append(atoms)

data = array(data)

nStates = data.shape[0]
nAtoms  = data.shape[1]

HHx = []
HOx = []
for s in range(0, nStates):
    ER = array([norm(data[s][0][1] - data[s][1][1]), data[s][0][2], data[s][0][3]])
    HHx.append(ER)
    
    if(data.shape[1] > 2):
        ER = array([norm(data[s][2][1] - data[s][0][1]), data[s][0][2], data[s][0][3]]) 
        HOx.append(ER)
        ER = array([norm(data[s][2][1] - data[s][1][1]), data[s][0][2], data[s][0][3]])
        HOx.append(ER)
    
HHx = array(HHx) 
HOx = array(HOx)

print "number of states: ", nStateFiles
print "number of atoms: ", nAtoms









    



number of states:  300
number of atoms:  2

Potential Energy Surface



In [8]:

    
#%matplotlib gtk
rcParams["figure.figsize"] = (12,8)
rcParams["font.size"] = 15

y_int = integrate.cumtrapz(HHx[:,2], HHx[:,0], initial=0)

plot(HHx[:,0],HHx[:,1], '--g')
plot(HHx[:,0],y_int+HHx[0,1], '--r')
title("Energy vs HH-bond length")
xlabel("$X_{HH}$ [$a_0$]")
ylabel("E[a.u.]")

#if(data.shape[1] > 2):
 #   figure()
  #  plot(HOx[:,0],HOx[:,1], 'o')
   # title("Energy vs HO-bond length")
    #xlabel("$X_{OH}$ [$a_0$]")
    #ylabel("E[a.u.]")









    



---------------------------------------------------------------------------
NameError                                 Traceback (most recent call last)
<ipython-input-8-6c7e5473d4b7> in <module>()
      3 rcParams["font.size"] = 15
      4 
----> 5 y_int = integrate.cumtrapz(HHx[:,2], HHx[:,0], initial=0)
      6 
      7 plot(HHx[:,0],HHx[:,1], '--g')

NameError: name 'integrate' is not defined

Force



In [88]:

    
plot(HHx[:,0],HHx[:,2], '--g')
title("Froce vs HH-bond length")
xlabel("$X_{HH}$ [$a_0$]")
ylabel("F[a.u.]")

#if(data.shape[1] > 2):
 #   figure()
 #   plot(HOx[:,0],HOx[:,2], 'o')
 #   title("Force vs HO-bond length")
 #   xlabel("$X_{OH}$ [$a_0$]")
 #   ylabel("F[a.u.]")









    Out[88]:





<matplotlib.text.Text at 0xc28f150>



In [89]:

    
data1 = loadtxt("/home/milad/kurs/energyForce.dat")
X = data1[:,0]
E = data1[:,1]
F = data1[:,2]



In [90]:

    
plot(X,E, 'g--')
#F_int = integrate.cumtrapz(F,X, initial=0.0)
#plot(X,F_int, 'b--')
title("Energy vs HH-bond length")
xlabel("$X_{HH}$ [$a_0$]")
ylabel("E[a.u.]")









    Out[90]:





<matplotlib.text.Text at 0xc222e10>



In [91]:

    
Fn = array(diff(E[:-1])/diff(X[:-1]))
axhline(y=0.0,color='r')
axvline(x=1.385,color='r')
print max(Fn - F[:-2])
plot(X,F, 'g--')
plot(X[:-2],Fn, 'b--')
title("Force vs HH-bond length")
xlabel("$X_{HH}$ [$a_0$]")
ylabel("F[a.u.]")









    



0.0205788254096






    Out[91]:





<matplotlib.text.Text at 0xcba3150>



In [92]:

    
figure()
step = 1.0
t = arange(0,len(HHx[:,0]),1)* 2.42e-2 * step 
plot(t, HHx[:,0])
title("HH-bond length vs time")
xlabel("Time [fs]")
ylabel("$X_{HH}$[$a_0$]")

if(data.shape[1] > 2):
    t = arange(0,len(HOx[:,0]),1)* 2.42e-2 * step 
    figure()
    plot(t, HOx[:,0])
    title("HO-bond length vs time")
    xlabel("Time [fs]")
    ylabel("$X_{OH}$ [$a_0$]")



In [16]:

    
a = array([])



In [18]:

    
plot(a[:,0],a[:,1],'--b',label ="num")
plot(a[:,0],a[:,2],'--g',label ="analytic")
legend()









    Out[18]:





<matplotlib.legend.Legend at 0x55107d0>