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

Generated by: LCOV version 2.0-1