File: emdplot1d.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 (90 lines) | stat: -rw-r--r-- 2,287 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

! Copyright (C) 2015 D. Ernsting, S. Dugdale and J. K. Dewhurst.
! This file is distributed under the terms of the GNU General Public License.
! See the file COPYING for license details.

subroutine emdplot1d(emds)
use modmain
use modpw
use modomp
implicit none
! arguments
real(4), intent(in) :: emds(nhkmax,nkpt)
! local variables
integer nh(3),ip,n,i,j,nthd
real(8) vl1(3),vl2(3),vl3(3)
real(8) vc1(3),vc2(3),vc3(3),t1
! allocatable arrays
real(8), allocatable :: x(:),f1(:),f2(:)
! external functions
real(8) rfhkintp,fintgt
external rfhkintp,fintgt
! generate the 1D plotting points: use only the first segment
call plotpt1d(bvec,2,npp1d,vvlp1d,vplp1d,dvp1d,dpp1d)
! compute two vectors orthogonal to each other and the plotting vector; these
! are the directions to be used for integration
vl1(:)=vvlp1d(:,2)-vvlp1d(:,1)
call r3mv(bvec,vl1,vc1)
t1=sqrt(vc1(1)**2+vc1(2)**2+vc1(3)**2)
if (t1.lt.epslat) then
  write(*,*)
  write(*,'("Error(emdplot1d): zero length plotting vector")')
  write(*,*)
  stop
end if
vc1(:)=vc1(:)/t1
i=1
do j=2,3
  if (abs(vc1(j)).lt.abs(vc1(i))) i=j
end do
vc2(:)=0.d0
vc2(i)=1.d0
t1=dot_product(vc1,vc2)
vc2(:)=vc2(:)-t1*vc1(:)
t1=sqrt(vc2(1)**2+vc2(2)**2+vc2(3)**2)
vc2(:)=vc2(:)/t1
call r3cross(vc1,vc2,vc3)
! integration directions in lattice coordinates
call r3mv(binv,vc2,vl2)
call r3mv(binv,vc3,vl3)
! determine the number of integration points
nh(:)=int(hkmax*sqrt(avec(1,:)**2+avec(2,:)**2+avec(3,:)**2)/pi)+1
n=2*maxval(nh(:)*ngridk(:))
allocate(x(n))
do i=1,n
  t1=2.d0*dble(i-1)/dble(n-1)-1.d0
  x(i)=t1*hkmax
end do
open(50,file='EMD1D.OUT',form='FORMATTED')
write(*,*)
! loop over plotting points along 1D line
call omp_hold(npp1d,nthd)
!$OMP PARALLEL DEFAULT(SHARED) &
!$OMP PRIVATE(f1,f2,i,j,vl1,t1) &
!$OMP NUM_THREADS(nthd)
!$OMP DO ORDERED
do ip=1,npp1d
  allocate(f1(n),f2(n))
  do i=1,n
    do j=1,n
      vl1(:)=vplp1d(:,ip)+x(i)*vl2(:)+x(j)*vl3(:)
      f1(j)=rfhkintp(vl1,emds)
    end do
    f2(i)=fintgt(-2,n,x,f1)
  end do
  t1=fintgt(-2,n,x,f2)
!$OMP ORDERED
  write(*,'("Info(emdplot1d): done ",I6," of ",I6," points")') ip,npp1d
  write(50,'(2G18.10)') dpp1d(ip),t1
  flush(50)
!$OMP END ORDERED
  deallocate(f1,f2)
end do
!$OMP END DO
!$OMP END PARALLEL
call omp_free(nthd)
close(50)
deallocate(x)
return
end subroutine