From 4191627f83469d253e6ee0bb6058ea4d5cff1e35 Mon Sep 17 00:00:00 2001 From: kamischi <karl-michael.schindler@web.de> Date: Wed, 12 Jun 2024 14:49:40 +0200 Subject: [PATCH] first try with WFW version --- .gitignore | 1 + source/WfW seteps/Notes.txt | 41 ++++++++++++++ source/WfW seteps/myEels2-seteps.f90 | 51 +++++++++++++++++ source/WfW seteps/myEels20-seteps.f90 | 58 +++++++++++++++++++ source/WfW seteps/myEels3-seteps.f90 | 53 ++++++++++++++++++ source/WfW seteps/myEels4-seteps.f90 | 53 ++++++++++++++++++ tests/seteps/Makefile | 42 +++++++++----- tests/seteps/setepsdriver-WFW.f90 | 80 +++++++++++++++++++++++++++ 8 files changed, 366 insertions(+), 13 deletions(-) create mode 100644 source/WfW seteps/Notes.txt create mode 100755 source/WfW seteps/myEels2-seteps.f90 create mode 100755 source/WfW seteps/myEels20-seteps.f90 create mode 100755 source/WfW seteps/myEels3-seteps.f90 create mode 100755 source/WfW seteps/myEels4-seteps.f90 create mode 100644 tests/seteps/setepsdriver-WFW.f90 diff --git a/.gitignore b/.gitignore index 82806a6..ccf9b96 100644 --- a/.gitignore +++ b/.gitignore @@ -22,3 +22,4 @@ tests/oxlog/eelsou tests/setepslog/eelsf90 tests/seteps/build-f90 tests/seteps/build-f77 +tests/seteps/build-WFW diff --git a/source/WfW seteps/Notes.txt b/source/WfW seteps/Notes.txt new file mode 100644 index 0000000..b6e68f1 --- /dev/null +++ b/source/WfW seteps/Notes.txt @@ -0,0 +1,41 @@ +Changes to seteps by WFW + +myEels2 → myEels3 +----------------- + +Change for Drude term + +38,43c39,44 +< else ! Negative TO mode means: treat as term added to epsilon +< if (osc(3,j) > 0) then ! it is a Lorentz oscillator +< addeps = addeps + osc(1,j)**2 * osc(2,j) /dcmplx(osc(1,j)**2 + wn2, wn*osc(3,j)) +< else ! it is a Drude term +< addeps = addeps - osc(1,j)**2/dcmplx(wn2, -1*wn*osc(3,j)) +< end if +--- +> +> else if (osc(1,j) < 0.) then! Negative TO mode means: _Additive_ Lorentz oscillator with Q +> addeps = addeps + osc(1,j)**2 * osc(2,j) /dcmplx(osc(1,j)**2 + wn2, wn*osc(3,j)) +> +> else ! osc(1,j) = 0 -> it is a Drude term +> addeps = addeps - dcmplx(osc(1,j+m)**2, wn*(osc(3,j)-osc(3,j+m))) /dcmplx(wn2, wn*osc(3,j)) + +myEels3 → myEels4 +----------------- + +Fix of final formula + +49c49 +< eps(l) = epsinf(l) * nomi / deno + addeps +--- +> eps(l) = epsinf(l) * (nomi / deno + addeps) ! brackets changed by HHe 230915 + +myEels4 → myEels20 +------------------ + +Fix of sign for additive Lorentz oscillator + +41c43 +< addeps = addeps + osc(1,j)**2 * osc(2,j) /dcmplx(osc(1,j)**2 + wn2, wn*osc(3,j)) +--- +> addeps = addeps + osc(1,j)**2 * osc(2,j) /dcmplx(osc(1,j)**2 - wn2, -1*wn*osc(3,j)) ! Sign of imaginary part changed (WFW) diff --git a/source/WfW seteps/myEels2-seteps.f90 b/source/WfW seteps/myEels2-seteps.f90 new file mode 100755 index 0000000..de49c75 --- /dev/null +++ b/source/WfW seteps/myEels2-seteps.f90 @@ -0,0 +1,51 @@ +subroutine seteps(nos, osc_size, osc, epsinf, wn, nLayer, eps) + + ! ****************************************************************** + ! * set up long-wavelength dielectric functions of the layers for * + ! * the present frequency wn (in cm**-1) * + ! ****************************************************************** + + ! implicit none + integer, intent(in) :: nLayer + integer, dimension(nLayer),intent(in) :: nos + integer, intent(in) :: osc_size + double precision, dimension(3,osc_size),intent(in) :: osc + !f2py depend(osc_size) osc + double precision, dimension(nLayer),intent(in) :: epsinf + double precision, intent(in) :: wn + double complex, dimension(nLayer), intent(out) :: eps + !f2py depend(nLayer) nos, epsinf, eps + + double complex :: nomi, deno, addeps + double precision :: wn2, b + integer :: flag + logical debugFirstRun + common /control/ debugFirstRun + + j = 0 + do l = 1, nLayer ! loop over different thin film layers + m = nos(l)/2 ! m number of TO modes = offset to reach the LO mode in the joint TO-LO list + nomi = dcmplx(1.0d0, 0.0d0) + deno = dcmplx(1.0d0, 0.0d0) + addeps = dcmplx(0.0d0, 0.0d0) + wn2 = wn**2 + do k = 1, m ! loop over all TO modes + j = j + 1 + if (osc(1,j) > 0.) then ! positive TO mode: 'Kurosa' form + b = wn/osc(1, j+m) + nomi = nomi * osc(1, j+m)**2 * (1.0 - b * dcmplx(b, osc(3,j+m)/osc(1, j+m)) ) + deno =deno * (osc(1,j)**2 - wn * dcmplx( wn, osc(3,j) ) ) + else ! Negative TO mode means: treat as term added to epsilon + if (osc(3,j) > 0) then ! it is a Lorentz oscillator + addeps = addeps + osc(1,j)**2 * osc(2,j) /dcmplx(osc(1,j)**2 + wn2, wn*osc(3,j)) + else ! it is a Drude term + addeps = addeps - osc(1,j)**2/dcmplx(wn2, -1*wn*osc(3,j)) + end if + end if + enddo + j = j+m ! we have already looped over the LO modes, therefore increase the index + eps(l) = epsinf(l) * nomi / deno + addeps + enddo + debugFirstRun = .false. + return +end subroutine seteps diff --git a/source/WfW seteps/myEels20-seteps.f90 b/source/WfW seteps/myEels20-seteps.f90 new file mode 100755 index 0000000..96ebd18 --- /dev/null +++ b/source/WfW seteps/myEels20-seteps.f90 @@ -0,0 +1,58 @@ +module seteps_mod +contains +subroutine seteps(nos, osc_size, osc, epsinf, wn, nLayer, eps) + + ! ****************************************************************** + ! * set up long-wavelength dielectric functions of the layers for * + ! * the present frequency wn (in cm**-1) * + ! ****************************************************************** + + implicit none + + integer, intent(in) :: nLayer + integer, dimension(nLayer),intent(in) :: nos + integer, intent(in) :: osc_size + double precision, dimension(3,osc_size),intent(in) :: osc + !f2py depend(osc_size) osc + double precision, dimension(nLayer),intent(in) :: epsinf + double precision, intent(in) :: wn + double complex, dimension(nLayer), intent(out) :: eps + !f2py depend(nLayer) nos, epsinf, eps + + double complex :: nomi, deno, addeps + double precision :: wn2, b + integer j, k, l, m + logical debugFirstRun + + common /control/ debugFirstRun + + j = 0 + do l = 1, nLayer ! loop over different thin film layers + m = nos(l)/2 ! m number of TO modes = offset to reach the LO mode in the joint TO-LO list + nomi = dcmplx(1.0d0, 0.0d0) + deno = dcmplx(1.0d0, 0.0d0) + addeps = dcmplx(0.0d0, 0.0d0) + wn2 = wn**2 + do k = 1, m ! loop over all TO modes + j = j + 1 + + if (osc(1,j) > 0.) then ! positive TO mode: 'Kurosawa' form: _Multiplicative_ phonon mode + b = wn/osc(1, j+m) + nomi = nomi * osc(1, j+m)**2 * (1.0 - b * dcmplx(b, osc(3,j+m)/osc(1, j+m)) ) + deno =deno * (osc(1,j)**2 - wn * dcmplx( wn, osc(3,j) ) ) + + else if (osc(1,j) < 0.) then! Negative TO mode means: _Additive_ Lorentz oscillator with Q + addeps = addeps + osc(1,j)**2 * osc(2,j) /dcmplx(osc(1,j)**2 - wn2, -1*wn*osc(3,j)) ! Sign of imaginary part changed (WFW) + + else ! osc(1,j) = 0 -> it is a Drude term + addeps = addeps - dcmplx(osc(1,j+m)**2, wn*(osc(3,j)-osc(3,j+m))) /dcmplx(wn2, wn*osc(3,j)) + end if + + enddo + j = j+m ! we have already looped over the LO modes, therefore increase the index + eps(l) = epsinf(l) * (nomi / deno + addeps) ! brackets changed by HHe 230915 + enddo + debugFirstRun = .false. + return +end subroutine seteps +end module seteps_mod diff --git a/source/WfW seteps/myEels3-seteps.f90 b/source/WfW seteps/myEels3-seteps.f90 new file mode 100755 index 0000000..4c41acc --- /dev/null +++ b/source/WfW seteps/myEels3-seteps.f90 @@ -0,0 +1,53 @@ +subroutine seteps(nos, osc_size, osc, epsinf, wn, nLayer, eps) + + ! ****************************************************************** + ! * set up long-wavelength dielectric functions of the layers for * + ! * the present frequency wn (in cm**-1) * + ! ****************************************************************** + + ! implicit none + integer, intent(in) :: nLayer + integer, dimension(nLayer),intent(in) :: nos + integer, intent(in) :: osc_size + double precision, dimension(3,osc_size),intent(in) :: osc + !f2py depend(osc_size) osc + double precision, dimension(nLayer),intent(in) :: epsinf + double precision, intent(in) :: wn + double complex, dimension(nLayer), intent(out) :: eps + !f2py depend(nLayer) nos, epsinf, eps + + double complex :: nomi, deno, addeps + double precision :: wn2, b + integer :: flag + logical debugFirstRun + common /control/ debugFirstRun + + j = 0 + do l = 1, nLayer ! loop over different thin film layers + m = nos(l)/2 ! m number of TO modes = offset to reach the LO mode in the joint TO-LO list + nomi = dcmplx(1.0d0, 0.0d0) + deno = dcmplx(1.0d0, 0.0d0) + addeps = dcmplx(0.0d0, 0.0d0) + wn2 = wn**2 + do k = 1, m ! loop over all TO modes + j = j + 1 + + if (osc(1,j) > 0.) then ! positive TO mode: 'Kurosawa' form: _Multiplicative_ phonon mode + b = wn/osc(1, j+m) + nomi = nomi * osc(1, j+m)**2 * (1.0 - b * dcmplx(b, osc(3,j+m)/osc(1, j+m)) ) + deno =deno * (osc(1,j)**2 - wn * dcmplx( wn, osc(3,j) ) ) + + else if (osc(1,j) < 0.) then! Negative TO mode means: _Additive_ Lorentz oscillator with Q + addeps = addeps + osc(1,j)**2 * osc(2,j) /dcmplx(osc(1,j)**2 + wn2, wn*osc(3,j)) + + else ! osc(1,j) = 0 -> it is a Drude term + addeps = addeps - dcmplx(osc(1,j+m)**2, wn*(osc(3,j)-osc(3,j+m))) /dcmplx(wn2, wn*osc(3,j)) + end if + + enddo + j = j+m ! we have already looped over the LO modes, therefore increase the index + eps(l) = epsinf(l) * nomi / deno + addeps + enddo + debugFirstRun = .false. + return +end subroutine seteps diff --git a/source/WfW seteps/myEels4-seteps.f90 b/source/WfW seteps/myEels4-seteps.f90 new file mode 100755 index 0000000..92b13a8 --- /dev/null +++ b/source/WfW seteps/myEels4-seteps.f90 @@ -0,0 +1,53 @@ +subroutine seteps(nos, osc_size, osc, epsinf, wn, nLayer, eps) + + ! ****************************************************************** + ! * set up long-wavelength dielectric functions of the layers for * + ! * the present frequency wn (in cm**-1) * + ! ****************************************************************** + + ! implicit none + integer, intent(in) :: nLayer + integer, dimension(nLayer),intent(in) :: nos + integer, intent(in) :: osc_size + double precision, dimension(3,osc_size),intent(in) :: osc + !f2py depend(osc_size) osc + double precision, dimension(nLayer),intent(in) :: epsinf + double precision, intent(in) :: wn + double complex, dimension(nLayer), intent(out) :: eps + !f2py depend(nLayer) nos, epsinf, eps + + double complex :: nomi, deno, addeps + double precision :: wn2, b + integer :: flag + logical debugFirstRun + common /control/ debugFirstRun + + j = 0 + do l = 1, nLayer ! loop over different thin film layers + m = nos(l)/2 ! m number of TO modes = offset to reach the LO mode in the joint TO-LO list + nomi = dcmplx(1.0d0, 0.0d0) + deno = dcmplx(1.0d0, 0.0d0) + addeps = dcmplx(0.0d0, 0.0d0) + wn2 = wn**2 + do k = 1, m ! loop over all TO modes + j = j + 1 + + if (osc(1,j) > 0.) then ! positive TO mode: 'Kurosawa' form: _Multiplicative_ phonon mode + b = wn/osc(1, j+m) + nomi = nomi * osc(1, j+m)**2 * (1.0 - b * dcmplx(b, osc(3,j+m)/osc(1, j+m)) ) + deno =deno * (osc(1,j)**2 - wn * dcmplx( wn, osc(3,j) ) ) + + else if (osc(1,j) < 0.) then! Negative TO mode means: _Additive_ Lorentz oscillator with Q + addeps = addeps + osc(1,j)**2 * osc(2,j) /dcmplx(osc(1,j)**2 + wn2, wn*osc(3,j)) + + else ! osc(1,j) = 0 -> it is a Drude term + addeps = addeps - dcmplx(osc(1,j+m)**2, wn*(osc(3,j)-osc(3,j+m))) /dcmplx(wn2, wn*osc(3,j)) + end if + + enddo + j = j+m ! we have already looped over the LO modes, therefore increase the index + eps(l) = epsinf(l) * (nomi / deno + addeps) ! brackets changed by HHe 230915 + enddo + debugFirstRun = .false. + return +end subroutine seteps diff --git a/tests/seteps/Makefile b/tests/seteps/Makefile index e2dbee5..a104149 100644 --- a/tests/seteps/Makefile +++ b/tests/seteps/Makefile @@ -13,43 +13,59 @@ ifeq "$(GFORTAN_VERSION_GE_5)" "1" endif # implicit rules for compiling fortran files +build-f77/%.o: ../../source/split-f77/%.f + $(FC) $(FFLAGS) -c -o $@ $< + build-f90/%.o: ../../source/f90/%.f90 $(FC) $(FFLAGS) -c -Jbuild-f90 -o $@ $< build-f90/%_mod.mod: ../../source/f90/%.f90 - $(FC) $(FFLAGS) -c -Jbuild-f90 -o build-f90/seteps.o $< + $(FC) $(FFLAGS) -c -Jbuild-f90 -o build-f90/seteps.o $< -build-f77/%.o: ../../source/split-f77/%.f - $(FC) $(FFLAGS) -c -o $@ $< +MYEELS20_SETEPS = ../../source/WFW seteps/myEels20-seteps.f90 + +build-WFW/seteps.o: $(wildcard ../../source/WFW\\ seteps/myEels20-seteps.f90) + $(FC) $(FFLAGS) -c -Jbuild-WFW -o $@ '$(MYEELS20_SETEPS)' +build-WFW/seteps_mod.mod: $(wildcard ../../source/WFW\\ seteps/myEels20-seteps.f90) + $(FC) $(FFLAGS) -c -Jbuild-WFW -o build-WFW/seteps.o '$(MYEELS20_SETEPS)' -all: build-f90/setepsdriver build-f77/setepsdriver +all: build-f77/setepsdriver build-f90/setepsdriver build-WFW/setepsdriver + +build-f77/setepsdriver: setepsdriver-f77.f90 build-f77/seteps.o + $(FC) $(FFLAGS) -o $@ $^ build-f90/setepsdriver: setepsdriver-f90.f90 build-f90/seteps.o build-f90/seteps_mod.mod $(FC) $(FFLAGS) -Jbuild-f90 -o $@ $< build-f90/seteps.o -build-f77/setepsdriver: setepsdriver-f77.f90 build-f77/seteps.o - $(FC) $(FFLAGS) -o $@ $^ +build-WFW/setepsdriver: setepsdriver-WFW.f90 build-WFW/seteps.o build-WFW/seteps_mod.mod + $(FC) $(FFLAGS) -Jbuild-WFW -o $@ $< build-WFW/seteps.o -test: build-f90/setepsdriver build-f77/setepsdriver +test: build-f90/setepsdriver build-f77/setepsdriver build-WFW/setepsdriver for NUMBER in 001 002 003 004 005 006 007 008 009 010 011 012 013 014 015 016 017 018 019 020 021 022 023 024 025 ; do \ echo " eels input file number: "$$NUMBER ; \ - cp ../../tests/inputFiles/eelsin$$NUMBER build-f90/eelsin ; \ - build-f90/setepsdriver ; \ - mv build-f90/epsout results-f90/epsout$$NUMBER ; \ - rm build-f90/eelsin ; \ cp ../../tests/inputFiles/eelsin$$NUMBER build-f77/eelsin ; \ build-f77/setepsdriver ; \ mv build-f77/epsout results-f77/epsout$$NUMBER ; \ rm build-f77/eelsin ; \ + cp ../../tests/inputFiles/eelsin$$NUMBER build-f90/eelsin ; \ + build-f90/setepsdriver ; \ + mv build-f90/epsout results-f90/epsout$$NUMBER ; \ + rm build-f90/eelsin ; \ + cp ../../tests/inputFiles/eelsin$$NUMBER build-WFW/eelsin ; \ + build-WFW/setepsdriver ; \ + mv build-WFW/epsout results-WFW/epsout$$NUMBER ; \ + rm build-WFW/eelsin ; \ done clean: cleanbuild cleanresults cleanbuild: - rm -rf build-f90/* rm -rf build-f77/* + rm -rf build-f90/* + rm -rf build-WFW/* cleanresults: - rm -rf results-f90/* rm -rf results-f77/* + rm -rf results-f90/* + rm -rf results-WFW/* .PHONY: all clean test diff --git a/tests/seteps/setepsdriver-WFW.f90 b/tests/seteps/setepsdriver-WFW.f90 new file mode 100644 index 0000000..2584bf7 --- /dev/null +++ b/tests/seteps/setepsdriver-WFW.f90 @@ -0,0 +1,80 @@ +program setepsdriver + + use seteps_mod + + implicit none + + double precision :: wmin, wmax, dw, wn + + double precision, allocatable :: osc(:,:), epsinf(:) + double complex, allocatable :: eps(:) + integer, allocatable :: nos(:) + integer :: jos, nLayer, osc_size, neps, nper + character (len = 10), allocatable :: name(:) + character (len = 10) mode + integer :: old_size_1, old_size_2 + double precision, allocatable :: tmp_osc(:,:) + + integer :: i, j, k, l, noPoints + character (len = 1) tab + + write(*,*) 'setepsdriver-WFW start' + tab = char(9) + noPoints = 200 + + open(unit = 11, file = 'build-WFW/eelsin') + read(11, *) + read(11, *) + read(11, *) + read(11, *) + read(11, *) wmin + read(11, *) wmax + read(11, *) dw + read(11, *) + read(11, *) + read(11, *) nLayer, nper, mode + nLayer = nLayer + if (nper == -1) then + nLayer = nLayer + 1 + nper = 1 + endif + jos = 0 + allocate (name(nLayer), epsinf(nLayer), nos(nLayer)) + do l = 1, nLayer + if (l <= nLayer) then + read(11, *) name(l) + endif + read(11, *) epsinf(l), nos(l) + if (nos(l) <= 0) then + else + do j = 1, nos(l) + if (.not. allocated(osc)) then + allocate(osc(3, nos(l))) + else if (j == 1) then + old_size_1 = size(osc, 1) + old_size_2 = size(osc, 2) + allocate(tmp_osc(old_size_1, old_size_2 + nos(l))) + tmp_osc(1:old_size_1, 1:old_size_2) = osc + deallocate(osc) + call move_alloc(tmp_osc, osc) + endif + jos = jos + 1 + read(11, *) (osc(k, jos), k = 1, 3) + enddo + endif + enddo + close (unit = 11) + + noPoints = 1 + int((wmax - wmin) / dw) + wn = wmin + allocate (eps(nLayer)) + + open(unit = 12, file = 'build-WFW/epsout') + do i = 0, noPoints + call seteps(nos, osc_size, osc, epsinf, wn, nLayer, eps) + write (12, '(30g15.7)') wn, (eps(j), j = 1, neps) + wn = wn + dw + end do + close (unit = 12) + +end program setepsdriver -- GitLab