Line data Source code
1 : !--------------------------------------------------------------------------------
2 : ! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
3 : ! This file is part of FLEUR and available as free software under the conditions
4 : ! of the MIT license as expressed in the LICENSE file in more detail.
5 : !--------------------------------------------------------------------------------
6 :
7 : c********************************************************
8 : c Calculate the electric polarization.
9 : c Frank Freimuth, October 2006
10 : c********************************************************
11 : module m_wann_dipole2
12 : use m_juDFT
13 : contains
14 0 : subroutine wann_dipole2(
15 0 : > jspins,pos,omtil,natd,
16 : > l_nocosoc)
17 : use m_wann_readcenters
18 : implicit none
19 : integer,intent(in) :: jspins
20 : real,intent(in) :: omtil
21 : integer,intent(in) :: natd
22 : real,intent(in) :: pos(3,natd)
23 : logical,intent(in) :: l_nocosoc
24 :
25 : integer :: jspin,jspins2
26 : logical :: l_file
27 0 : integer,allocatable:: ind(:)
28 : integer :: nwf,nwfs
29 0 : real,allocatable :: distance(:,:)
30 0 : real :: moment(3,jspins)
31 : real :: moment2(3)
32 : real :: elemchargmu,bohrtocm
33 : integer :: shifted,ishift,ind1,ind2
34 0 : real,allocatable :: centers(:,:)
35 0 : real,allocatable :: spreads(:)
36 :
37 0 : call timestart("wann_dipole2")
38 0 : open(204,file='dipole_out')
39 0 : write(204,*)"*****************************************"
40 0 : write(204,*)"Calculation of the electric polarization."
41 0 : write(204,*)"*****************************************"
42 0 : elemchargmu=1.60217646e-13
43 0 : bohrtocm=0.529177e-8
44 :
45 0 : jspins2=jspins
46 0 : if(l_nocosoc)jspins2=1
47 0 : do jspin=1,jspins2
48 0 : inquire(file='proj',exist=l_file)
49 0 : IF(.NOT.l_file) CALL juDFT_error("proj not found",calledby
50 0 : + ="wann_dipole2")
51 0 : open(203,file='proj',form='formatted')
52 0 : read(203,*)nwfs
53 0 : print*,"nwfs=",nwfs
54 0 : allocate(ind(nwfs))
55 0 : allocate(distance(3,nwfs))
56 0 : allocate(centers(3,nwfs))
57 0 : allocate(spreads(nwfs))
58 0 : do nwf=1,nwfs
59 0 : read(203,*)ind(nwf)
60 0 : read(203,*)
61 : enddo
62 0 : close(203)
63 0 : write(204,*)"spin=",jspin
64 0 : call wann_readcenters(nwfs,jspin,centers,spreads)
65 0 : write(204,*)"*****centers:******"
66 0 : do nwf=1,nwfs
67 0 : write(204,fmt=1000)nwf,centers(:,nwf),spreads(nwf)
68 : enddo
69 0 : moment(:,jspin)=0.0
70 0 : do nwf=1,nwfs
71 0 : distance(:,nwf)=pos(:,ind(nwf))-centers(:,nwf)
72 0 : moment(:,jspin)=moment(:,jspin)+distance(:,nwf)
73 : enddo
74 0 : if((jspins.eq.1).and..not.l_nocosoc)
75 0 : & moment(:,jspin)=moment(:,jspin)*2
76 0 : moment(:,jspin)=moment(:,jspin)/omtil
77 0 : moment(:,jspin)=moment(:,jspin)*elemchargmu/((bohrtocm)**2)
78 0 : write(204,*)"polarization due to shift of wannierorbitals"
79 0 : write(204,*)moment(:,jspin),"uC/cm2"
80 0 : deallocate(ind,distance,centers,spreads)
81 : enddo!jspin
82 : 1000 format(2x,'WF centre and spread',
83 : & i5,2x,'(',f10.6,',',f10.6,',',f10.6,' )',f15.8)
84 0 : inquire(file='dipole',exist=l_file)
85 0 : if(.not.l_file)then
86 0 : write(204,*)"no file dipole found"
87 : else
88 0 : moment2(:)=0.0
89 0 : open(224,file='dipole',form='formatted')
90 0 : read(224,*)shifted
91 0 : do ishift=1,shifted
92 0 : read(224,*)ind1,ind2 !electron taken from 1 and moved to 2
93 0 : moment2(:)=moment2(:)+pos(:,ind1)-pos(:,ind2)
94 : enddo
95 0 : close(224)
96 0 : moment2(:)=moment2(:)/omtil
97 0 : moment2(:)=moment2(:)*elemchargmu/(bohrtocm)**2
98 : write(204,*)"polarization due to electrons
99 0 : & being catched by other atoms"
100 0 : write(204,*)moment2(:),"uC/cm2"
101 :
102 : endif
103 0 : write(204,*)"total moment:"
104 0 : if(jspins2.eq.2)moment(:,1)=moment(:,1)+moment(:,2)
105 0 : if(l_file)then
106 0 : moment(:,1)=moment(:,1)+moment2(:)
107 : endif
108 0 : write(204,*)moment(:,1),"uC/cm2"
109 0 : close(204)
110 :
111 0 : call timestop("wann_dipole2")
112 0 : end subroutine wann_dipole2
113 : end module m_wann_dipole2
|
