LCOV - code coverage report
Current view: top level - wannier - wann_dipole.f (source / functions) Hit Total Coverage
Test: combined.info Lines: 0 109 0.0 %
Date: 2019-09-08 04:53:50 Functions: 0 1 0.0 %

          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             :       module m_wann_dipole
       8             :       use m_juDFT
       9             :       contains
      10           0 :       subroutine wann_dipole(
      11           0 :      >               jspins,omtil,natd,pos,
      12             :      >               amat,
      13           0 :      >               ntype,neq,zatom)
      14             : c***************************************
      15             : c     Calculate electronic polarization.
      16             : c     Frank Freimuth
      17             : c***************************************
      18             :       use m_wann_readcenters
      19             :       implicit none
      20             :       integer,intent(in)           :: jspins
      21             :       real,intent(in)              :: omtil
      22             :       integer,intent(in)           :: natd
      23             :       real,intent(in)              :: pos(3,natd)
      24             :       real,intent(in)              :: amat(3,3)
      25             :       integer,intent(in)           :: ntype
      26             :       integer,intent(in)           :: neq(ntype)
      27             :       real,intent(in)              :: zatom(ntype)
      28             : 
      29             :       integer                      :: nwfs,i
      30             :       integer                      :: num_atoms,num_wann,num_wann2
      31             :       integer                      :: j,k,ind,jj,jspin
      32           0 :       real,allocatable             :: ioncharge(:)
      33           0 :       character(len=2),allocatable :: namat(:)
      34             :       character(len=2)             :: symbol
      35           0 :       real,allocatable             :: wann_centers1(:,:)
      36           0 :       real,allocatable             :: wann_centers2(:,:)
      37             :       integer,allocatable          :: wann_of_at(:)
      38             :       real                         :: polarization,charge
      39             :       integer                      :: num_symbols
      40             :       real                         :: ioni_polari(3)
      41             :       real                         :: shifted_polari(3)
      42             :       real                         :: size_polari
      43             :       real                         :: smallest_polari
      44             :       real                         :: final_polari(3)
      45             :       real                         :: electroni_polari(3,2)
      46             :       real                         :: elemchargmu,bohrtocm
      47             :       character*2                  :: namat2(0:103)
      48             :       character(len=2)             :: spin12(0:2)
      49             :       data spin12/'  ', '.1', '.2'/
      50             :       character(len=6)             :: filename
      51             :       logical                      :: l_file
      52             : 
      53             :       DATA namat2/'va',' H','He','Li','Be',
      54             :      +     ' B',' C',' N',' O',' F','Ne',
      55             :      +     'Na','Mg','Al','Si',' P',' S','Cl','Ar',' K','Ca','Sc','Ti',
      56             :      +     ' V','Cr','Mn','Fe','Co','Ni','Cu','Zn','Ga','Ge','As','Se',
      57             :      +     'Br','Kr','Rb','Sr',' Y','Zr','Nb','Mo','Tc','Ru','Rh','Pd',
      58             :      +     'Ag','Cd','In','Sn','Sb','Te',' I','Xe','Cs','Ba','La','Ce',
      59             :      +     'Pr','Nd','Pm','Sm','Eu','Gd','Tb','Dy','Ho','Er','Tm','Yb',
      60             :      +     'Lu','Hf','Ta',' W','Re','Os','Ir','Pt','Au','Hg','Tl','Pb',
      61             :      +     'Bi','Po','At','Rn','Fr','Ra','Ac','Th','Pa',' U','Np','Pu',
      62             :      +     'Am','Cm','Bk','Cf','Es','Fm','Md','No','Lw'/
      63             : 
      64           0 :       elemchargmu=1.60217646e-13
      65           0 :       bohrtocm=0.529177e-8
      66             : 
      67           0 :       open(666,file='polarization_out')
      68             : 
      69           0 :       num_atoms=0
      70           0 :       do j=1,ntype
      71           0 :          do k=1,neq(j)
      72           0 :             num_atoms=num_atoms+1
      73             :          enddo
      74             :       enddo
      75             : 
      76           0 :       allocate( namat(num_atoms) )
      77           0 :       ind=0
      78           0 :       do j=1,ntype
      79           0 :          do k=1,neq(j)
      80           0 :             ind=ind+1
      81           0 :             namat(ind)=namat2(nint(zatom(j)))
      82             :          enddo
      83             :       enddo
      84             : 
      85           0 :       do j=1,num_atoms
      86           0 :          print*,namat(j)," pos3=",pos(3,j)
      87             :       enddo
      88             : 
      89           0 :       allocate( ioncharge(num_atoms) )
      90           0 :       ioncharge=0.0
      91           0 :       open(400,file='IONS')
      92           0 :       read(400,*)num_symbols
      93           0 :       do j=1,num_symbols
      94           0 :          read(400,fmt=333)symbol,charge
      95           0 :          write(*,fmt=333)symbol,charge
      96           0 :          do k=1,num_atoms
      97           0 :             if(namat(k)==symbol)then
      98           0 :                ioncharge(k)=charge
      99             :             endif
     100             :          enddo
     101             :       enddo
     102           0 :       close(400)
     103             :  333  format(a2,1x,f10.6)
     104             : 
     105           0 :       open(300,file='ioncharge')
     106           0 :       do j=1,num_atoms
     107           0 :          write(300,*)ioncharge(j)
     108             :       enddo
     109           0 :       close(300)
     110             : 
     111           0 :       open(300,file='ioncharge')
     112           0 :       do j=1,num_atoms
     113           0 :          read(300,*)ioncharge(j)
     114             :       enddo
     115           0 :       close(300)
     116             : 
     117           0 :       ioni_polari=0.0
     118           0 :       do j=1,num_atoms
     119             :          ioni_polari(:)=
     120           0 :      &        ioni_polari(:)+ioncharge(j)*pos(:,j)
     121             :       enddo
     122             : 
     123           0 :       ioni_polari=ioni_polari/omtil*elemchargmu/((bohrtocm)**2)
     124             : 
     125           0 :       print*,"ioni_polari=",ioni_polari,"uC/cm2"
     126           0 :       write(666,*) "ioni_polari=",ioni_polari,"uC/cm2"
     127             : 
     128             : c*************************************************
     129             : c..reading the proj.1 / proj.2 / proj file
     130             : c*************************************************
     131           0 :       do j=1,0,-1
     132           0 :          inquire(file=trim('proj'//spin12(j)),exist=l_file)
     133           0 :          if(l_file)then
     134           0 :             filename='proj'//spin12(j)
     135           0 :             exit
     136             :          endif
     137             :       enddo
     138           0 :       if(l_file)then
     139           0 :         open (203,file=trim(filename),status='old')
     140           0 :         rewind (203)
     141             :       else
     142             :          CALL juDFT_error("no proj/proj.1/proj.2",
     143           0 :      +        calledby ="wann_dipole")
     144             :       endif  
     145           0 :       read(203,*)num_wann
     146           0 :       close(203)
     147           0 :       print*,"number of wannier functions= ",num_wann
     148           0 :       allocate(wann_centers1(3,num_wann))
     149             : 
     150           0 :       if(jspins.eq.2)then
     151           0 :         do j=2,0,-1
     152           0 :           inquire(file=trim('proj'//spin12(j)),exist=l_file)
     153           0 :           if(l_file)then
     154           0 :             filename='proj'//spin12(j)
     155           0 :             exit
     156             :           endif
     157             :         enddo
     158           0 :         if(l_file)then
     159           0 :           open (203,file=trim(filename),status='old')
     160           0 :           rewind (203)
     161             :         else
     162             :            CALL juDFT_error("no proj/proj.1/proj.2",calledby
     163           0 :      +          ="wann_dipole")
     164             :         endif  
     165           0 :         read(203,*)num_wann2
     166           0 :         close(203)
     167           0 :         print*,"number of wannier functions spin 2= ",num_wann2
     168           0 :         allocate(wann_centers2(3,num_wann2))
     169             :       endif
     170             : 
     171           0 :       electroni_polari=0.0
     172             :       call wann_readcenters(
     173             :      >         num_wann,1,
     174           0 :      <         wann_centers1(:,:))
     175             : 
     176           0 :       do k=1,num_wann
     177             :          electroni_polari(:,1)=
     178           0 :      &     electroni_polari(:,1)-wann_centers1(:,k)
     179             :       enddo
     180             : 
     181           0 :       if(jspins.eq.2)then
     182             :          call wann_readcenters(
     183             :      >         num_wann2,2,
     184           0 :      <         wann_centers2(:,:))
     185           0 :         do k=1,num_wann2
     186             :            electroni_polari(:,2)=
     187           0 :      &        electroni_polari(:,2)-wann_centers2(:,k)
     188             :         enddo
     189             :       endif   
     190             : 
     191             : 
     192           0 :       if(jspins.eq.1) electroni_polari = electroni_polari*2.0
     193             :       electroni_polari=
     194           0 :      &  electroni_polari/omtil*elemchargmu/((bohrtocm)**2)
     195             : 
     196           0 :       do j=1,jspins
     197           0 :          print*,"spin ",j," electronic_polarization=",
     198           0 :      &                electroni_polari(:,j)
     199           0 :          write(666,*)"spin ",j," electronic_polarization=",
     200           0 :      &                electroni_polari(:,j)
     201             :       enddo
     202             :       
     203           0 :       ioni_polari(:)=ioni_polari(:)+electroni_polari(:,1)
     204           0 :       if(jspins.eq.2)ioni_polari(:)=ioni_polari(:)
     205           0 :      &              +electroni_polari(:,2)
     206             : 
     207           0 :       print*,"total polarization=",ioni_polari(:)
     208           0 :       write(666,*)"total polarization=",ioni_polari(:)
     209             : !     Check if the polarization may be reduced by adding 
     210             : !     contributions due to electrons shifted by primitive
     211             : !     lattice translations.
     212           0 :       final_polari(:)=ioni_polari(:)
     213             :       size_polari=sqrt( (final_polari(1))**2 +
     214             :      +                  (final_polari(2))**2 +
     215           0 :      +                  (final_polari(3))**2 )
     216           0 :       smallest_polari=size_polari
     217           0 :       do i=-5,5
     218           0 :          do j=-5,5
     219           0 :             do k=-5,5
     220             :                shifted_polari(:)=ioni_polari(:)+
     221             :      +                  elemchargmu/((bohrtocm)**2)/omtil *
     222           0 :      *                  (amat(:,1)*k+amat(:,2)*j+amat(:,3)*i)
     223             :                size_polari=sqrt( (shifted_polari(1))**2 +
     224             :      +                           (shifted_polari(2))**2 +
     225           0 :      +                           (shifted_polari(3))**2 )
     226           0 :                if(size_polari.lt.smallest_polari)then
     227           0 :                   final_polari=shifted_polari
     228           0 :                   smallest_polari=size_polari
     229             :                endif
     230             :             enddo
     231             :          enddo
     232             :       enddo
     233           0 :       print*,"final polarization after shifting=",final_polari(:)
     234           0 :       write(666,*)"final polarization after shifting=",final_polari(:)
     235           0 :       close(666)
     236             : 
     237             : 
     238             : 
     239           0 :       end subroutine wann_dipole
     240             :       end module m_wann_dipole

Generated by: LCOV version 1.13