module fft_mkl
	use constant,only : nx
	use mkl_dfti
	
	implicit none
	type(dfti_descriptor), POINTER :: dfft_x
	integer(4) :: Status
	contains
		subroutine initial_fft_mkl
		Status = DftiCreateDescriptor(dfft_x, DFTI_DOUBLE,DFTI_COMPLEX,1,nx)
		status = DftiSetValue(dfft_x,DFTI_FORWARD_SCALE,1.d0/dble(nx)) 
		Status = DftiCommitDescriptor(dfft_x)

		end subroutine initial_fft_mkl
		subroutine free_fft_mkl

		Status = DftiFreeDescriptor(dfft_x)
		end  subroutine free_fft_mkl

end module fft_mkl

	subroutine gradient_1d(psi,psi_gd)
	use constant, only : nx,Per,eps,order
	implicit none
	complex(8),intent(in)  :: psi(0:nx)
	complex(8),intent(out) :: psi_gd(0:nx)

	complex(8) :: psi1(0:nx)
	complex(8) :: psi2(0:nx)
	complex(8) :: psi3(0:nx)
	complex(8) :: psi4(0:nx)
	integer(4) :: forward,backward

	integer(4) :: i

	forward = 1
	backward = 0
	if(order.eq.2) then 

		psi1 = psi
		psi2 = psi
		call split_4th_1d(psi1,forward)
		call split_4th_1d(psi2,backward)
		psi_gd = (psi1-psi2)/(2.d0*Per*eps)

	elseif(order.eq.4) then

		psi1 = psi
		call split_4th_1d(psi1,forward)
		psi2 = psi1
		call split_4th_1d(psi2,forward)
		psi_gd = -psi2 + 8.d0*psi1

		psi1 = psi
		call split_4th_1d(psi1,backward)
		psi2 = psi1
		call split_4th_1d(psi2,backward)

		psi_gd = (psi_gd-8.d0*psi1+psi2)/(1.2d1*Per*eps)


	elseif(order.eq.6) then

		psi1 = psi
		call split_4th_1d(psi1,forward)
		psi2 = psi1
		call split_4th_1d(psi2,forward)
		psi3 = psi2
		call split_4th_1d(psi3,forward)
		psi_gd = psi3-9.d0*psi2+4.5d1*psi1

		psi1 = psi
		call split_4th_1d(psi1,backward)
		psi2 = psi1
		call split_4th_1d(psi2,backward)
		psi3 = psi2
		call split_4th_1d(psi3,backward)
		psi_gd = (psi_gd-psi3+9.d0*psi2-4.5d1*psi1)/(6.d1*Per*eps)

	elseif(order.eq.8) then

		psi1 = psi
		call split_4th_1d(psi1,forward)
		psi2 = psi1
		call split_4th_1d(psi2,forward)
		psi3 = psi2
		call split_4th_1d(psi3,forward)
		psi4 = psi3
		call split_4th_1d(psi4,forward)
		psi_gd = 3.d0*psi4-3.2d1*psi3+1.68d2*psi2-6.72d2*psi1

		psi1 = psi
		call split_4th_1d(psi1,backward)
		psi2 = psi1
		call split_4th_1d(psi2,backward)
		psi3 = psi2
		call split_4th_1d(psi3,backward)
		psi4 = psi3
		call split_4th_1d(psi4,backward)
		psi_gd = (psi_gd-(3.d0*psi4-3.2d1*psi3+1.68d2*psi2-6.72d2*psi1))/(-8.4d2*Per*eps)

	endif

	end subroutine

	subroutine split_2nd_1d(psi,flag)
	use constant,only: im,nx,umx,Per,dtMicro, Vtemp,beta,eps
	use fft_mkl
	implicit none
	complex(8) :: psi(0:nx)
	integer(4) :: flag

	integer(4) :: i,k,ntMicro

	ntMicro = nint(Per/dtMicro)

	if(flag.eq.1) then
		do k = 1,ntMicro
			do i = 0,nx-1
				psi(i) = exp(-im*beta*eps*dtMicro/2.d0*abs(psi(i))**2.d0)*psi(i)
			enddo

			Status = DftiComputeForward(dfft_x, psi(0:nx-1))
			do i = 0,nx-1
				psi(i) =exp(-im*dtMicro*Vtemp(i))*psi(i)
			enddo
			Status = DftiComputeBackward(dfft_x, psi(0:nx-1))
	
			do i = 0,nx-1
				psi(i) = exp(-im*beta*eps*dtMicro/2.d0*abs(psi(i))**2.d0)*psi(i)
			enddo
		enddo
	else
		do k = 1,ntMicro
			do i = 0,nx-1
				psi(i) = exp(im*beta*eps*dtMicro/2.d0*abs(psi(i))**2.d0)*psi(i)
			enddo

			Status = DftiComputeForward(dfft_x, psi(0:nx-1))
			do i = 0,nx-1
				psi(i) =exp(im*dtMicro*Vtemp(i))*psi(i)
			enddo
			Status = DftiComputeBackward(dfft_x, psi(0:nx-1))
	
			do i = 0,nx-1
				psi(i) = exp(im*beta*eps*dtMicro/2.d0*abs(psi(i))**2.d0)*psi(i)
			enddo
		enddo

	endif

	end subroutine split_2nd_1d

	subroutine split_4th_1d(psi,flag)
	use constant
	use fft_mkl
	implicit none
	complex(8) :: psi(0:nx)
	integer(4) :: flag

	integer(4) :: i,k,ntMicro
	ntMicro = nint(Per/dtMicro)

	if(flag.eq.1) then
		do k = 1,ntMicro
			do i = 0,nx-1
				psi(i)= exp(-im*beta*eps*2.d0*w1*dtMicro*Vext(i)*abs(psi(i))**2.d0)*psi(i)
			enddo

			Status = DftiComputeForward(dfft_x, psi(0:nx-1))
			do i = 0,nx-1
				psi(i) = exp(-im*2.d0*w2*dtMicro*Vtemp(i))*psi(i)
			enddo
			Status = DftiComputeBackward(dfft_x, psi(0:nx-1))

			do i = 0,nx-1
				psi(i) = exp(-im*beta*eps*2.d0*w3*dtMicro*Vext(i)*abs(psi(i))**2.d0)*psi(i)
			enddo

			Status = DftiComputeForward(dfft_x, psi(0:nx-1))
			do i = 0,nx-1
				psi(i) = exp(-im*2.d0*w4*dtMicro*Vtemp(i))*psi(i)
			enddo
			Status = DftiComputeBackward(dfft_x, psi(0:nx-1))

			do i = 0,nx-1
				psi(i) = exp(-im*beta*eps*2.d0*w3*dtMicro*Vext(i)*abs(psi(i))**2.d0)*psi(i)
			enddo

			Status = DftiComputeForward(dfft_x, psi(0:nx-1))
			do i = 0,nx-1
				psi(i) = exp(-im*2.d0*w2*dtMicro*Vtemp(i))*psi(i)
			enddo
			Status = DftiComputeBackward(dfft_x, psi(0:nx-1))

			do i = 0,nx-1
				psi(i) = exp(-im*beta*eps*2.d0*w1*dtMicro*Vext(i)*abs(psi(i))**2.d0)*psi(i)
			enddo

		enddo

	else
		do k = 1,ntMicro
			do i = 0,nx-1
				psi(i)= exp(im*beta*eps*2.d0*w1*dtMicro*Vext(i)*abs(psi(i))**2.d0)*psi(i)
			enddo

			Status = DftiComputeForward(dfft_x, psi(0:nx-1))
			do i = 0,nx-1
				psi(i) = exp(im*2.d0*w2*dtMicro*Vtemp(i))*psi(i)
			enddo
			Status = DftiComputeBackward(dfft_x, psi(0:nx-1))

			do i = 0,nx-1
				psi(i) = exp(im*beta*eps*2.d0*w3*dtMicro*Vext(i)*abs(psi(i))**2.d0)*psi(i)
			enddo

			Status = DftiComputeForward(dfft_x, psi(0:nx-1))
			do i = 0,nx-1
				psi(i) = exp(im*2.d0*w4*dtMicro*Vtemp(i))*psi(i)
			enddo
			Status = DftiComputeBackward(dfft_x, psi(0:nx-1))

			do i = 0,nx-1
				psi(i) = exp(im*beta*eps*2.d0*w3*dtMicro*Vext(i)*abs(psi(i))**2.d0)*psi(i)
			enddo

			Status = DftiComputeForward(dfft_x, psi(0:nx-1))
			do i = 0,nx-1
				psi(i) = exp(im*2.d0*w2*dtMicro*Vtemp(i))*psi(i)
			enddo
			Status = DftiComputeBackward(dfft_x, psi(0:nx-1))

			do i = 0,nx-1
				psi(i) = exp(im*beta*eps*2.d0*w1*dtMicro*Vext(i)*abs(psi(i))**2.d0)*psi(i)
			enddo
		enddo

	endif

	end subroutine split_4th_1d

