File: dielectric_tdrt.f90

package info (click to toggle)
elkcode 5.4.24-2
  • links: PTS, VCS
  • area: main
  • in suites: buster, sid
  • size: 12,840 kB
  • sloc: f90: 48,415; fortran: 22,457; perl: 965; makefile: 384; sh: 369; python: 105; ansic: 67
file content (147 lines) | stat: -rw-r--r-- 3,662 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147

! Copyright (C) 2017 J. K. Dewhurst, S. Sharma and E. K. U. Gross.
! This file is distributed under the terms of the GNU General Public License.
! See the file COPYING for license details.

subroutine dielectric_tdrt
use modmain
use modtddft
implicit none
! local variables
integer nts,its0,its
integer iw,i,j
real(8) a0,t0,t1,t2
complex(8) z1,z2
character(256) fname
! allocatable arrays
real(8), allocatable :: f1(:),f2(:),g(:),w(:),jt(:,:)
complex(8), allocatable :: ew(:,:),jw(:,:),eps(:)
! external functions
real(8) splint
external splint
! initialise global variables
call init0
call init1
! set up the frequency grid
allocate(w(nwplot))
t1=wplot(2)/dble(nwplot)
do iw=1,nwplot
  w(iw)=t1*dble(iw)
end do
! find the starting time step for the Fourier transform
its0=1
do its=2,ntimes
  if (times(its).gt.t0tdlr) then
    its0=its-1
    exit
  end if
end do
nts=ntimes-its0+1
! compute the electric field from E = -1/c dA/dt and Fourier transform
allocate(f1(ntimes),f2(ntimes),g(ntimes))
allocate(ew(nwplot,3))
t0=-1.d0/solsc
do i=1,3
  f1(:)=afieldt(i,:)
  call fderiv(1,ntimes,times,f1,g)
! constant term corresponding to instantaneous A-field at t=0
  a0=f1(1)
  if (task.eq.480) then
! Fourier transform E(t) numerically
    do iw=1,nwplot
      do its=its0,ntimes
        t1=g(its)
        t2=w(iw)*times(its)
        f1(its)=t1*cos(t2)
        f2(its)=t1*sin(t2)
      end do
      t1=splint(nts,times(its0),f1(its0))
      t2=splint(nts,times(its0),f2(its0))
      ew(iw,i)=t0*cmplx(t1+a0,t2,8)
    end do
  else
! analytic Fourier transform of E(t) assumed to be a delta function at t=0
    t1=t0*a0
    ew(1:nwplot,i)=t1
  end if
end do
! read in the total current from file
allocate(jt(3,ntimes))
call readcurrent(jt)
! divide by the unit cell volume
jt(:,:)=jt(:,:)/omega
! remove constant term from current
do i=1,3
  do its=its0,ntimes
    f1(its)=jt(i,its)
  end do
  t1=splint(nts,times(its0),f1(its0))
  t1=t1/(ntimes*dtimes)
  do its=its0,ntimes
    jt(i,its)=jt(i,its)-t1
  end do
end do
! filter the high-frequency components from the current using an exponentially
! decaying function
do its=its0,ntimes
  t1=times(its)-t0tdlr
  t1=exp(-swidth*t1)
  jt(:,its)=t1*jt(:,its)
end do
! Fourier transform the current
allocate(jw(nwplot,3))
do i=1,3
  do iw=1,nwplot
    do its=its0,ntimes
      t1=jt(i,its)
      t2=w(iw)*times(its)
      f1(its)=t1*cos(t2)
      f2(its)=t1*sin(t2)
    end do
    t1=splint(nts,times(its0),f1(its0))
    t2=splint(nts,times(its0),f2(its0))
    jw(iw,i)=cmplx(t1,t2,8)
  end do
end do
deallocate(f1,f2,g,jt)
! compute the dielectric function and write to file
allocate(eps(nwplot))
do i=1,3
  do j=1,3
    do iw=1,nwplot
      z1=jw(iw,i)
      z2=ew(iw,j)
      t1=abs(dble(z2))+abs(aimag(z2))
      if (t1.gt.1.d-8) then
        z1=z1/z2
      else
        z1=0.d0
      end if
      z1=fourpi*cmplx(-aimag(z1),dble(z1),8)
      z1=z1/w(iw)
      if (i.eq.j) z1=z1+1.d0
      eps(iw)=z1
    end do
    write(fname,'("EPSILON_TDRT_",2I1,".OUT")') i,j
    open(50,file=trim(fname),form='FORMATTED')
    do iw=1,nwplot
      write(50,'(2G18.10)') w(iw),dble(eps(iw))
    end do
    write(50,'("     ")')
    do iw=1,nwplot
      write(50,'(2G18.10)') w(iw),aimag(eps(iw))
    end do
    close(50)
  end do
end do
write(*,*)
write(*,'("Info(dielectric_tdrt):")')
write(*,'(" dielectric tensor determined from real-time evolution")')
write(*,'(" written to EPSILON_TDRT_ij.OUT for components i,j = 1,2,3")')
write(*,*)
write(*,'("(Note that only those components which are not orthogonal to the")')
write(*,'(" applied A-field will be calculated correctly)")')
deallocate(w,ew,jw,eps)
return
end subroutine