MOLECULAR DYNAMICS - TIME DEPENDENT BEHAVIIOUR MODELING, SIMULATION and VISUALIZATION

'FIRST  IDENTIFY POTENTIAL, see the Morse potential work earlier -HCl

Dehclvib=7.10647  (10^(-19))  kg (m^2)  (s^(-2))

Evib=2.8535      (10^(-20))  kg (m^2)  (s^(-2))

Dehcl=Dehclvib+Evib

De=7.39182`*^-19

Dehcl/(1.6  (10^(-19)))

c=2.99799 (10^8)  m/s

μ=1.62747 (10^(-27)) kg

νe =2.899 (10^3) /cm

νe =2.899 (10^5) /m

β=PowerExpand[2 Pi c νe (Sqrt[(μ/(2  Dehcl))])]

ββ=β

'get out the units

ββ=1.8118576344890297`*^10

(((1.8118576344890297`*^10))^(-1))

'this is 0.55  angstrems

re=1.275  (10^(-10))

'review

De=Dehcl=((7.39182`*^-19 kg (m^2))/(s^2))

De=7.39182`*^-19

re

'
'note that these two forms are equal

De  (1-  Exp[-ββ  (r-re)])^2=De(1 -2   Exp[-ββ  (r-re) + Exp[- 2 ββ  (r-re)])

[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[

[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[[

'Morse function gives potential

f[r_]=De  (1-  Exp[-ββ  (r-re)])^2

Plot[De  (1-  Exp[-ββ  (r-re)])^2,{r,0.1 (10^(-10)),4.1 (10^(-10))}]

'how does th e potential behave when changing parameters

Plot[{De  (1-  Exp[-ββ  (r-0.6*  re)])^2,De  (1-  Exp[-ββ  (r-1*  re)])^2,De  (1-  Exp[-ββ  (r-1.3*  re)])^2},{r,0.1 (10^(-10)),4.1 (10^(-10))}]

f[r_,k_]:=De  (1-  Exp[-ββ  (r-k re)])^2

Plot3D[f[r,k],{k,0.1,1.2},{r,0.1 (10^(-10)),4.1 (10^(-10))}]

f[n_,x_]:= Abs[((1/Pi)^(1/4)HermiteH[n,x])/(E^(x^2/2)Sqrt[2^n n!])]^2

Plot[Evaluate@
  Append[Table[f[n,x]+n+1/2,{n,0,7}],x^2/2],{x,-4,4},Filling->Table[n->n-1/2,{n,1,8}]]

Plot[De  (1-  Exp[-ββ  (r-re)])^2,{r,0.1 (10^(-10)),4.1 (10^(-10))}]

'SECONDLY , use POTENTIAL TO MODEL FORCE AND SIMULATE DYNAMICS (position, velocity)
''force = first derivatiove of V =  2 De (E^(-(re-re) ββ)) (1-(E^(-(r-re) ββ))) ββ, =2 D beta at re'
'total potential is V=V12 +V23 etc
'here V(HCl)
' we will work on simulation in time and use Verlet algorithm
' h is time shift here (not Planck), size of time step
'say h= 1 femtosec
' for i-th nucleus i = H , Cl  or A, B

' so for each step, n,   at time ( t+n h) we calculate position xi and force and then what else is needed


' xi(t+h) =xi(t) +h vi(t) + h^2 Fi(xi(t))/ 2 mi
'vi(t+h)= vi(t) +  h ( Fi(xi(t+h) )+ Fi(xi(t))

' we will calculate first xi,  then force, then vi  and iteratively again for next step
'you need some initial conditions (real, physical)

ββ

'force= - 2 De (E^(-(r-re) ββ)) (1-(E^(-(r-re) ββ))) ββ

'let us see  time scales for frequencies from spectroscopic values in cm-1 energy

eV to J/h is frequency, all to -1 = time

frequency=2900

(((frequency/8065)*(1.6 10^-19)/(6.6 10^-34))^-1)

tt= (((frequency/8065.5)*(1.6 10^-19)/(6.6*10^-34))^(-1))

'that is 72 fs

'corresponds  in frequency s^-1, Hz

frq=((2899/8065.5)*(1.6 10^-19)/(6.6*10^-34))

tt*frq

' time scale will be

h=1. 10^-15

'CHECK MASSES AND PUT ACCURATE IN

mA=1*1836*9.1 10^(-31)

mB=35.5*1836*9.1 10^(-31)

'FORM  Arrays

'for Dynamics you need position-x, velocity -v , and force-f

Array[xA,256];

Array[vA,256];
Array[fA,256];

Array[xC,256];

Array[xB,256];

Array[vB,256];

Array[vC,256];

Array[fB,256];

Array[fC,256];

Array[rAB,256];

Array[rBC,256];

h=1. 10^-15

re

' A is hydrogen around x=0, B=Cl aroun\d equillibrium position

'inital positions estimate

xA[1]=-0.001 10^(-10)

xB[1]=1.277 10^(-10)

mA

mB

xA[1]

xB[1]

re

rAB[1]=re

fA[1]=- 2 De (E^(-(rAB[1]-re) ββ)) (1-(E^(-(rAB[1]-re) ββ))) ββ

fB[1]=-(- 2 De (E^(-(rAB[1]-re) ββ)) (1-(E^(-(rAB[1]-re) ββ))) ββ)

vA[1]=0.

vB[1]=0.

xA[1]+h * vA[1]+(h^2)* fA[1]/(2 mA)

xB[1]+h * vB[1]+(h^2)* fB[1]/(2 mB)

rAB[n+1]=1.0000000000000001`*^-11-(-1.0000000000000001`*^-11)

- 2 De (E^(-(rAB[n+1]-re) ββ)) (1-(E^(-(rAB[n+1]-re) ββ))) ββ

re

xA[1]=-0.001 10^(-10)

xB[1]=1.2751 10^(-10)

Do[{
xA[n+1]=xA[n]+h * vA[n]+(h^2)* fA[n]/(2 mA),
xB[n+1]=xB[n]+h * vB[n]+(h^2)* fB[n]/(2 mB),
rAB[n+1]=xB[n+1]-xA[n+1],
fA[n+1]=- 2 De (E^(-(rAB[n+1]-re) ββ)) (1-(E^(-(rAB[n+1]-re) ββ))) ββ,
fB[n+1]=-(- 2 De (E^(-(rAB[n+1]-re) ββ)) (1-(E^(-(rAB[n+1]-re) ββ))) ββ),
vA[n+1]=vA[n]+h * (fA[n+1]+fA[n])/(2   mA),
vB[n+1]=vB[n]+h * (fB[n+1]+fB[n])/(2   mB)},
{n,1,150,1}]//N

vA[2]

vbt=Table[vB[i],{i,1,150}]

vat=Table[vA[i],{i,1,150}]

ListPlot[{vat,vbt}]

remember here force is  - dv/dx , from Morse potential, obtained through spectroscpic investigations

force is-((dV )/dx)

Plot[{-De (-2 (E^(-2 (r-re) ββ)) ββ+2 (E^(-(r-re) ββ)) ββ)},{r,0.1 (10^(-10)),4.1 (10^(-10))}]

pot1=Plot[{De  (- 2  Exp[-ββ  (r-re)]+ Exp[-ββ   2 (r-re)])},{r,0.1 (10^(-10)),4.1 (10^(-10))}]

in jouls

De

Do[{
xA[n+1]=xA[n]+h * vA[n]+(h^2)* fA[n]/(2 mA),
xB[n+1]=xB[n]+h * vB[n]+(h^2)* fB[n]/(2 mB),
rAB[n+1]=xB[n+1]-xA[n+1],
fA[n+1]=De (-2 (E^(-2 (rAB[n+1]-re) ββ)) ββ+2 (E^(-(rAB[n+1]-re) ββ)) ββ) ,
fB[n+1]=-(De (-2 (E^(-2 (rAB[n+1]-re) ββ)) ββ+2 (E^(-(rAB[n+1]-re) ββ)) ββ)),
vA[n+1]=vA[n]+h * (fA[n+1]+fA[n])/(2   mA),
vB[n+1]=vB[n]+h * (fB[n+1]+fB[n])/(2   mB)},
{n,1,150,1}]//N

h

vbt=Table[vB[i],{i,1,150}]

vat=Table[vA[i],{i,1,150}]

vA[1]

vB[1]

g1=ListPlot[{vat},PlotJoined->True]

g2=ListPlot[{vbt},PlotJoined->True]

Show[g1,g2]

eV to J/h is frequency, all to -1 = time

(((2899/8065)*(1.6 10^-19)/(6.6 10^-34))^-1)

vbt=Table[vB[i],{i,1,150}]

freqB=Abs[Fourier[vbt]]

ListPlot[freqB,PlotJoined->True,PlotRange->All]

eV to J/h is frequency, all to -1 = time

(((2899/8065)*(1.6 10^-19)/(6.6 10^-34))^-1)

'SEE IMPROVEMENTS