************************************************************************
PROGRAM PATH
************************************************************************
* Calculates the time evolution of wavefunction in real time
* using the Feyman Path Integral Method.
* Integration is approximated by Riemann summation
* Time evolution is converted into repeated matrix multiplications
* of the small time step propagator
************************************************************************
* Variables:
*
* K_dt(N,N) short time propagator matrix with N discretization points
*
*
*
************************************************************************
IMPLICIT NONE
INTEGER I,J,JP,K,L,ND
REAL*8 DX, DT, PI, TWOPI, ALPHA, A_HALF, WF_norm, X0, XD
REAL*8 DX_3, X_start, X_av, PSI_norm, W0, W1, U0, U1
REAL*8 K_dt, K_prop, K_w, PSI, PSI_1, C0, SUM0, SUM1
DIMENSION K_dt(0:1,0:1000,0:1000), K_prop(0:1,0:1000,0:1000)
DIMENSION K_w(0:1,0:1000,0:1000)
DIMENSION PSI(0:1,0:1000), PSI_1(0:1,0:1000)
DATA PI/3.141592653589793D0/
ND = 600
X0 = -4.0D0
XD = 4.0D0
ALPHA = 2.0D0
X_start = 0.75D0
A_HALF = ALPHA/2.0D0
DX = (XD-X0)/DFLOAT(ND)
WF_norm = (ALPHA/PI)**0.25D0
DX_3 = DX*DX*DX
TWOPI = 2.0D0*PI
DT = TWOPI/128.0D0
C0 = 1.0D0/DSQRT( 4.0D0*PI*DSIN(DT) )
DO K = 0,ND
PSI(0,K) = WF_norm*DEXP(-A_HALF*(X0+DX*DFLOAT(K)-X_start)**2)
PSI(1,K) = 0.0D0
ENDDO
DO K = 0,ND
DO L = 0,ND
W0 = ( ( (X0+DFLOAT(K)*DX)*(X0+DFLOAT(K)*DX) +
> (X0+DX*DFLOAT(L))*(X0+DX*DFLOAT(L)) ) -
> 2.0D0*(X0+DX*DFLOAT(K))*(X0+DX*DFLOAT(L)) )/(2.0D0*DT) -
> ( (X0+DFLOAT(K)*DX)*(X0+DFLOAT(K)*DX) +
> (X0+DX*DFLOAT(L))*(X0+DX*DFLOAT(L)) +
> 2.0D0*(X0+DX*DFLOAT(K))*(X0+DX*DFLOAT(L)) )*DT*0.125D0
K_dt(0,K,L) = C0*(DSIN(W0) + DCOS(W0))
K_dt(1,K,L) = C0*(DSIN(W0) - DCOS(W0))
ENDDO
ENDDO
X_av = 0.0D0
PSI_norm = 0.0D0
DO K=0,ND
PSI_norm = PSI_norm + PSI(0,K)*PSI(0,K) + PSI(1,K)*PSI(1,K)
X_av = X_av + (X0+DFLOAT(K)*DX)*( PSI(0,K)*PSI(0,K) +
> PSI(1,K)*PSI(1,K) )
ENDDO
PSI_norm = DX*PSI_norm
X_av = DX*X_av
print *, PSI_norm, X_av
DO K = 0,ND
DO L = 0,ND
K_w(0,K,L) = 0.0D0
K_w(1,K,L) = 0.0D0
DO I = 0,ND
W0 = K_dt(0,K,I)
W1 = K_dt(1,K,I)
U0 = K_dt(0,I,L)
U1 = K_dt(1,I,L)
K_w(0,K,L) = K_w(0,K,L) + (W0*U0 - W1*U1)
K_w(1,K,L) = K_w(1,K,L) + (W0*U1 + W1*U0)
ENDDO
ENDDO
ENDDO
DO K=0,ND
DO L=0,ND
K_prop(0,K,L) = DX*K_w(0,K,L)
K_prop(1,K,L) = DX*K_w(1,K,L)
ENDDO
ENDDO
DO K = 0,ND
DO L = 0,ND
K_w(0,K,L) = 0.0D0
K_w(1,K,L) = 0.0D0
DO I = 0,ND
W0 = K_dt(0,K,I)
W1 = K_dt(1,K,I)
U0 = K_prop(0,I,L)
U1 = K_prop(1,I,L)
K_w(0,K,L) = K_w(0,K,L) + (W0*U0 - W1*U1)
K_w(1,K,L) = K_w(1,K,L) + (W0*U1 + W1*U0)
ENDDO
ENDDO
ENDDO
DO K=0,ND
DO L=0,ND
K_prop(0,K,L) = DX*K_w(0,K,L)
K_prop(1,K,L) = DX*K_w(1,K,L)
ENDDO
ENDDO
DO K = 0,ND
DO L = 0,ND
K_w(0,K,L) = 0.0D0
K_w(1,K,L) = 0.0D0
DO I = 0,ND
W0 = K_dt(0,K,I)
W1 = K_dt(1,K,I)
U0 = K_prop(0,I,L)
U1 = K_prop(1,I,L)
K_w(0,K,L) = K_w(0,K,L) + (W0*U0 - W1*U1)
K_w(1,K,L) = K_w(1,K,L) + (W0*U1 + W1*U0)
ENDDO
ENDDO
ENDDO
DO K=0,ND
DO L=0,ND
K_prop(0,K,L) = DX*K_w(0,K,L)
K_prop(1,K,L) = DX*K_w(1,K,L)
ENDDO
ENDDO
DO K=0,ND
SUM0 = 0.0D0
SUM1 = 0.0D0
DO J=0,ND
DO I=0,ND
W0 = K_prop(0,I,J)
W1 = -K_prop(1,I,J)
U0 = K_prop(0,I,K)
U1 = K_prop(1,I,K)
SUM0 = SUM0 + DX*DX*(W0*U0 - W1*U1)
SUM1 = SUM1 + DX*DX*(W0*U1 + W1*U0)
ENDDO
ENDDO
print *, K, SUM0, SUM1
ENDDO
DO J=1,128
DO K=0,ND
PSI_1(0,K) = 0.0D0
PSI_1(1,K) = 0.0D0
DO L=0,ND
W0 = K_prop(0,K,L)
W1 = K_prop(1,K,L)
U0 = PSI(0,L)
U1 = PSI(1,L)
PSI_1(0,K) = PSI_1(0,K) + DX*(W0*U0 - W1*U1)
PSI_1(1,K) = PSI_1(1,K) + DX*(W0*U1 + W1*U0)
ENDDO
ENDDO
DO K=0,ND
PSI(0,K) = PSI_1(0,K)
PSI(1,K) = PSI_1(1,K)
c print *, K, PSI(0,K), PSI(1,K)
ENDDO
X_av = 0.0D0
PSI_norm = 0.0D0
DO K=0,ND
PSI_norm = PSI_norm + PSI(0,K)*PSI(0,K) + PSI(1,K)*PSI(1,K)
X_av = X_av + (X0+DFLOAT(K)*DX)*( PSI(0,K)*PSI(0,K) +
> PSI(1,K)*PSI(1,K) )
ENDDO
PSI_norm = DX*PSI_norm
X_av = DX*X_av
print *, J, PSI_norm, X_av
ENDDO
STOP
END