NTA has declared the answer keys and responses of UGC CSIR NET examination June 2021 held in February 2022. Check here: http://testservices.nic.in/examSys21/root/downloadadmitcard.aspx?appFormId=101162111
https://csirnet.nta.nic.in/WebInfo/Page/Page?PageId=1&LangId=P
Source code:
import numpy as np
import matplotlib.pyplot as plt
def Vpot(x):
return 0.5*x**2
a = float(input('Provide lower limit of the domain: '))
b = float(input('Provide upper limit of the domain: '))
N = int(input('Provide number of grid points: '))
x = np.linspace(a,b,N)
h = x[1]-x[0]
T = np.zeros((N-2)**2).reshape(N-2,N-2)
for i in range(N-2):
for j in range(N-2):
if i==j:
T[i,j]= -2
elif np.abs(i-j)==1:
T[i,j]=1
else:
T[i,j]=0
V = np.zeros((N-2)**2).reshape(N-2,N-2)
for i in range(N-2):
for j in range(N-2):
if i==j:
V[i,j]= Vpot(x[i+1])
else:
V[i,j]=0
H = -T/(2*h**2) + V
val,vec=np.linalg.eig(H)
z = np.argsort(val)
z = z[0:4]
energies=0.5*(val[z]/val[z][0])
print(energies)
plt.figure(figsize=(12,10))
for i in range(len(z)):
y = []
y = np.append(y,vec[:,z[i]])
y = np.append(y,0)
y = np.insert(y,0,0)
plt.plot(x,y,lw=3, label="{} ".format(i))
plt.xlabel('x', size=14)
plt.ylabel('$\psi$(x)',size=14)
plt.legend()
plt.title('Wavefunctions for a harmonic oscillator using finite difference method',size=14)
plt.show()
Result:
 |
| Wave functions |
 |
| Executing the code |
The following script uses the Finite Difference Method (FDM) to solve time-independent Schrodinger equation (TISE) to find out the eigenstates of H atom.
// Aim: Solve the s-wave Schrodinger equation for the ground state and the first excited
//state of the hydrogen atom.
//Here, m is the reduced mass of the electron. Obtain the energy eigenvalues and plot
// the corresponding wavefunctions. Remember that the ground state energy of the
//hydrogen atom is -13.6 eV. Take e = 3.795 (eVÅ)^1/2, ħc = 1973 (eVÅ) and m = 0.511x10^6 eV/c2.
close;
clear;
clc;
// declaring constant values
hbarc=1973;//Plancks constant h divided by 2*(pi) called as hbar=(h/2*pi). This factor when multiplied by speed of light c then hbar*c in the units of (eVÅ) comes out to be 1973;
mcsq=0.511*10^(6);//This is mass of electron*c ^2 we call it mcsq in units of (eV);
e = 3.795; // (eVÅ)^1/2
// We hereby input the 'r' values so as to obtain the potential V(r) as a function of 'r')
r_min=0 // in units of angstrom
r_max=18 // in units of angstrom
N = input("Input the number of intervals (should be around 500 to 750 for good computation)")
s = (r_max-r_min)/N; //step size
factor1=-(hbarc^2)/(2*mcsq*s^2); // this factor is (hbar^2*c^2/2m*c^2) divided by s^2 //k=(hbar_c*hbar_c)/(2*m)
// Kinetic energy matrix (Using central difference formula)
T=zeros(N-1,N-1)
for i=1:N-1
T(i,i)=-2
if i<N-1 then
T(i,i+1)=1
T(i+1,i)=1
end
end
T_matrix = factor1*T; // Kinetic Energy Matrix final (scaling done)
//Potential energy matrix
V_matrix=zeros(N-1,N-1)
for i=1:N-1
r(i)=r_min+i*s
V_matrix(i,i)=-(e*e)/r(i);
end
// Hamiltonian matrix
H_matrix=T_matrix+V_matrix
// energy eigenvalue and eigenstates
[u,eigen]=spec(H_matrix);
// displaying of the ground and first excited state energies
disp("Ground state energy (1S orbital) for hydrogen atom is (in eV) : ")
disp(eigen(1,1))
disp("First excited state (2S orbital) energy for hydrogen atom is (in eV): ")
disp(eigen(2,2))
disp("Second excited state (3S orbital) energy for hydrogen atom is (in eV): ")
disp(eigen(3,3))
disp("Third excited state (4S orbital) energy for hydrogen atom is (in eV): ")
disp(eigen(4,4))
rmatrix = [0;r;18]; //including the first and last point at which wavefunction is zero.
//Displaying the first four energy eigen values
figure(1);
//scf()
for n = 1:1:4
plot(rmatrix, eigen(n,n)*ones(N+1,1), 'linewidth',n)
hl=legend(['n = 1';'n= 2';'n=3';'n=4']);
title('Energy Eigen values for Hydrogen atom','fontsize',3)
xlabel('r (Angstrom)','fontsize',3)
ylabel('Eigen Value','fontsize',3)
end
//Plotting the Probability |Psi^2| as a function of r
figure(2);
//scf()
R_wave_1s=u(:,1)./r; //Radial wavefuction
R_wave_1s_final=[0;R_wave_1s;0]; //including the first and last point at which wavefunction is zero.
// plot of probability function
P_wave_1s= R_wave_1s_final.*R_wave_1s_final;
plot(rmatrix,P_wave_1s, 'linewidth',3)
title('Plot of Probability function for 1s orbital', 'fontsize',3)
xlabel('r (Angstrom)','fontsize',3)
ylabel('Probability','fontsize',3)
figure(3);
//scf()
R_wave_2s=u(:,2)./r;
R_wave_2s_final=[0;R_wave_2s;0]; //including the first and last point at which wavefunction is zero.
// plot of probability function
P_wave_2s= R_wave_2s_final.*R_wave_2s_final;
plot(rmatrix,P_wave_2s, 'linewidth',3)
title('Plot of Probability function for 2s orbital','fontsize',3)
xlabel('r (Angstrom)','fontsize',3)
ylabel('Probability','fontsize',3)
//Plotting the Probability density |Psi^2|*r^2 as a function of r
figure(4);
//scf()
plot(rmatrix,[P_wave_1s, P_wave_2s], 'linewidth',3)
title('Plot of Probability function for 1s and 2s orbital','fontsize',3)
xlabel('r (Angstrom)','fontsize',3)
ylabel('Probability','fontsize',3)
//Plotting the Probability of electron |Psi^2 r^2| as a function of r
figure(5);
//scf()
fp = 4*3.14
Rr_wave_1s=u(:,1).*fp; //Radial wavefuction
Rr_wave_1s_final=[0;Rr_wave_1s;0]; //including the first and last point at which wavefunction is zero.
// plot of probability function
Pr_wave_1s= Rr_wave_1s_final.*Rr_wave_1s_final;
plot(rmatrix,Pr_wave_1s, 'linewidth',3)
title('Plot of Probability Probability of finding the e for 1s orbital', 'fontsize',3)
xlabel('r (Angstrom)','fontsize',3)
ylabel('Probability','fontsize',3)
figure(6);
//scf()
Rr_wave_2s=u(:,2).*fp;
Rr_wave_2s_final=[0;Rr_wave_2s;0]; //including the first and last point at which wavefunction is zero.
// plot of probability function
Pr_wave_2s= Rr_wave_2s_final.*Rr_wave_2s_final;
plot(rmatrix,Pr_wave_2s, 'linewidth',3)
title('Plot of Probability of finding the e for 2s orbital','fontsize',3)
xlabel('r (Angstrom)','fontsize',3)
ylabel('Probability','fontsize',3)
figure(7);
//scf()
Rr_wave_3s=u(:,3).*fp;
Rr_wave_3s_final=[0;Rr_wave_3s;0]; //including the first and last point at which wavefunction is zero.
// plot of probability function
Pr_wave_3s= Rr_wave_3s_final.*Rr_wave_3s_final;
plot(rmatrix,Pr_wave_3s, 'linewidth',3)
title('Plot of Probability of finding the e for 3s orbital','fontsize',3)
xlabel('r (Angstrom)','fontsize',3)
ylabel('Probability','fontsize',3)
Source code:
c PROGRAM FOR FITTING GIVEN DATA TO A STRAIGHT LINE USING LEAST SQUARE METHOD
DIMENSION X(100),Y(100)
REAL X,Y
WRITE(*,*)'HOW MANY NUMBERS ARE THERE'
READ(*,*)N
WRITE(*,*)'ENTER DATA POINTS NOW'
READ(*,*)(X(I),Y(I),I=1,N)
SUMX=0.0
SUMY=0.0
SUMXX=0.0
SUMXY=0.0
DO 10 I=1,N
SUMX=SUMX+X(I)
SUMY=SUMY+Y(I)
SUMXX=SUMXX+X(I)**2
SUMXY=SUMXY+X(I)*Y(I)
10 CONTINUE
DENOM=FLOAT(N)*SUMXX-SUMX**2
A=(SUMY*SUMXX-SUMX*SUMXY)/DENOM
B=(FLOAT(N)*SUMXY-SUMX*SUMY)/DENOM
WRITE(*,*)'INTERCEPT OF BEST FITTING LINE ON Y-AXIS IS',A
WRITE(*,*)'SLOPE OF BEST FIT LINE IS',B
WRITE(*,*)'BEST FIT DATA POINTS ARE'
DO 20 I=1,N
Y(I)=A+B*X(I)
WRITE(*,*)X(I),Y(I)
20 CONTINUE
STOP
END
Source code:
c Lagrange's interpolation
dimension x(1000), f(1000)
write(*,*)'x,n'
read(*,*)a,n
write(*,*)'xi, fi'
read(*,*)(x(i),i=1,n)
write(*,*)'f'
read(*,*)(f(i),i=1,n)
c prod=(a-x(j))/(x(i)-x(j))
sum=0
prod=1
do i=1, n
do j=1, n
c prod=(a-x(j))/(x(i)-x(j))
if(j.NE.i) then
prod=(a-x(j))/(x(i)-x(j))
sum=sum+f(i)*prod
endif
enddo
enddo
write(*,*)'ans',sum
stop
end
c Bisection method
c Designed by Physics and Coding.
f(x)=j*x**2+k*x+l
22 write(*,*)'Give the quadratic equation coefficients a,b,c'
read(*,*)j,k,l
if(j.EQ.0) then
write(*,*)'This is a linear equation.'
d=-l/k
write(*,*)'The root is=>',d
goto 99
endif
11 write(*,*)'Give initial guesses x0,x1 & e'
read(*,*)x0,x1,e
y0=f(x0)
y1=f(x1)
i=0
if(y0*y1.GT.0) then
write(*,*)'The initial guesses are not suitable'
write(*,19)'Re-initialising the process...','....','...'
19 format(//,3x,a35,/,3x,a4,/,3x,a4)
goto 11
endif
do i=1,n
if((abs(x1-x0)/x1).GT.e) then
x2=(x1+x0)/2
y2=f(x2)
endif
if(y2*y0.GT.0) then
x0=x2
y0=y2
else
x1=x2
y1=y2
endif
end do
write(*,*)'The root, value of f(x) and no. of iterations are',x2,y2,i
99 write(*,55)
55 format(1x,//,'Would you like to re-calculate?',/)
write(*,*)'Input 1 to recalculate'
write(*,*)'Input 2 to try another equation'
write(*,*)'Input anything other to exit'
read(*,*)a
if(a.EQ.1) then
goto 11
endif
if(a.EQ.2) then
goto 22
endif
end
