0
$\begingroup$

I am hoping that someone could give me some insight about my problem. Currently I want to simulate the evolution of a gaussian wavepacket in a time-dependent potential using the Crank-Nicolson scheme. I have already the potential and wavefunction in dimensionless units. However I have the doubt of which time and frequency units I need to use. I want to run my simulation for a total time of 200 ps, and the frequency that I use is 3 THz. However my question is, can I just put the number 3 (instead of 3^12) as my frequency in THz, and my simulation time just as 200 (which will be in ps)? I am posting my python code with comments only on the relevant part of the question.

import scipy as sp
import numpy as np
from scipy import integrate, sparse, linalg,special
import scipy.sparse.linalg
import pylab as pl
ax=-10. #initial point
bx=10. #final point
tf=200 #final time in picoseconds
nx= 500#numero de puntos en el x-grid
Lx=bx-ax #total length
delx= Lx/nx #x-step size=Lx/nx
dt= 1.e-3#T/mt # time step delta t
omega=3# frequency is 3THz
mt=int((tf/omega)/dt)#
print (' MT = ', mt)
nux= 1j*dt/(4.0*(delx**2)) 

gridx = sp.zeros(nx) #space grid
igridx = sp.array(range(nx)) # time grid
psi = sp.zeros(nx)
pot = sp.zeros(nx)
gridx = delx*(igridx - nx/2)

###############These are just parameters#############
sigma=0.1
alpha= 2.614
D=1
R=2.95
r_e=.9699
R_e=R-2*r_e 
r_0=(1./alpha)*np.arccosh(0.5*np.exp(alpha*R_e/2.0))
x=alpha*gridx
X_e=alpha*R_e
E_barrier=D*( 1-2*np.exp(-alpha*R_e/2.) )**2


psi=(2*sp.pi*sigma**2)**(-0.25) *np.exp( -(x+alpha*r_0)**2 /(4*sigma**2) )


print(nux)
####### A+######################
Adiag=sp.empty(nx,dtype=complex) 
Asup=sp.empty(nx,dtype=complex) 
Asub=sp.empty(nx,dtype=complex) 
Adiag.fill(1+2.0*nux) 
Asup.fill(-nux)
Asub.fill(-nux)
aplus=sp.sparse.spdiags([Adiag,Asup,Asub],[0,1,-1],nx,nx) 
####### matrix A-########
adiag=sp.empty(nx,dtype=complex)
asup=sp.empty(nx,dtype=complex)
asub=sp.empty(nx,dtype=complex)
adiag.fill(1-2.0*nux)
asup.fill(nux)
asub.fill(nux)
aminus=sp.sparse.spdiags([adiag,asup,asub],[0,1,-1],nx,nx) 
#############################################################  

for t in range(mt) :
    R=R_ave+R_ave*0.1*np.sin(omega*t*dt)# omega is 3THz
    R_e=R-2*r_e
    pot=(75*( (1-np.exp(-alpha*(R_e/2 + gridx)))**2 + (1-np.exp(-alpha*(R_e/2 -gridx)))**2   - (1 - 2*np.exp(-alpha*R_e)) ))
    ##############
    psi= sp.exp(-1j*dt*pot)*psi
    psi = sp.sparse.linalg.bicg(aplus, aminus*psi)[0]
$\endgroup$
  • $\begingroup$ Which units are energy and distance in? $\endgroup$ – noah Feb 20 at 21:44
  • $\begingroup$ energy is in kcal/mol, distance is in angstroms, alpha is in inverse angstroms. I established the D=1 because when I divide the potential energy by D then it should be dimensionless. all the pontential units used in the calculation are dimensionless, but I just want to confirm if it is correct to use those values for time $\endgroup$ – Oliver Feb 27 at 16:25
  • $\begingroup$ kcal/mol is not an energy. And what is D? You can't just randomly declare things, they must fit together with other units to make sense. So for any useful help your question needs a lot more details about all the variables involved. $\endgroup$ – noah Feb 27 at 19:46

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.