1!
2! Copyright (C) 1996-2016	The SIESTA group
3!  This file is distributed under the terms of the
4!  GNU General Public License: see COPYING in the top directory
5!  or http://www.gnu.org/copyleft/gpl.txt.
6! See Docs/Contributors.txt for a list of contributors.
7!
8      subroutine obc( polR, ucell, dx, nspin)
9C *******************************************************************
10C Writes the born effective charges to file
11C by Tom Archer
12C Modified by Alberto Garcia, April 2007
13C ********* INPUT ***************************************************
14C real*8 dx             : atomic displacements (in Bohr)
15C real*8 polR(3,nspin)  : polarization along lattice vectors(Bohr)
16C real*8 ucell(3,3)     : cell vectors
17C integer nspin         : spin polarized calculation flag
18C ********** BEHAVIOUR **********************************************
19C On the first call (undisplaced coordinates), initial polarization
20C is calculated.
21C The born charges are calculated as the difference between the
22C polarization of the undisplaced and the displaced coordinates
23C divided by the displacement
24C
25C a phase is introduced in the KSV calculation which is removed before
26C the BC matrix is saved.
27C *******************************************************************
28
29      use precision, only: dp
30      use fdf,       only: fdf_string
31      use alloc,     only: re_alloc
32
33      implicit          none
34
35      integer, intent(in)  :: nspin
36      real(dp), intent(in) :: polR(3,nspin),
37     $                        ucell(3,3), dx
38      external io_assign, io_close
39
40c Internal
41
42      integer    igrd, ispin, ix, unit1
43
44      character(len=33), save :: fname
45      logical, save  ::    first_time = .true.
46
47      real(dp) :: dmod(3), phaseR(3), uR(3,3), pR(3), pxyz(3)
48
49      double precision, dimension(:,:), pointer, save :: pres
50
51!---------
52
53      do ispin=1,nspin
54         if (nspin.gt.1) then
55            if (ispin.eq.1) write(6,'(/,a)')
56     .           'obc: Macroscopic polarization for spin Up:'
57            if (ispin.eq.2) write(6,'(/,a)')
58     .           'obc: Macroscopic polarization for spin Down:'
59         endif
60         write(6,'(/,a,a)')
61     .        'obc: Macroscopic polarization per unit cell',
62     .        ' along the lattice vectors (Bohr): '
63         write(6,'(a,3f12.6)') 'obc:',(polR(ix,ispin),ix=1,3)
64      enddo
65
66      if (first_time) then
67
68        fname = trim(fdf_string('SystemLabel','siesta')) // '.BC'
69
70        nullify( pres )
71        call re_alloc( pres, 1, 3, 1, nspin, name='pres',
72     &                 routine='obc' )
73
74        call io_assign(unit1)
75        open( unit1, file=fname, status='unknown' )
76        rewind(unit1)
77        if (nspin.eq.1) then
78           write(unit1,'(a)') 'BC matrix'
79        else
80           write(unit1,'(a)') 'BC matrix for spin-up + spin-down'
81        endif
82        !  Set values of residual polarization
83        do ispin=1,nspin
84           do igrd=1,3
85              pres(igrd,ispin) = polR(igrd,ispin)
86           enddo
87        enddo
88        call io_close(unit1)
89
90        first_time = .false.
91
92      else
93
94          !matrix to convert from lattice vectors to xyz
95         do igrd=1,3
96            dmod(igrd)=dot_product(ucell(:,igrd),ucell(:,igrd))
97            dmod(igrd)=sqrt(dmod(igrd))
98            do ix=1,3
99               uR(ix,igrd)=ucell(ix,igrd)/dmod(igrd)
100            enddo
101         enddo
102
103         !     write BC matrix in lattice vectors
104         pR(:)=0.0_dp
105         pxyz(:)=0.0_dp
106         do ispin=1,nspin
107            do igrd=1,3
108               pR(igrd)=pR(igrd)+(polR(igrd,ispin)-pres(igrd,ispin))
109            enddo
110         enddo
111
112         !remove extra phase
113         do igrd=1,3
114            phaseR(igrd) = nint((pR(igrd)/dmod(igrd))) * dmod(igrd)
115            pR(igrd) = pR(igrd) - phaseR(igrd)
116         enddo
117
118         !get polarization in cartesian coordinates
119         do ix=1,3
120            pxyz(ix) = uR(ix,1)*pR(1)+ uR(ix,2)*pR(2)+uR(ix,3)*pR(3)
121         enddo
122
123         call io_assign(unit1)
124         open( unit1, file=fname, status='old',position="append",
125     $         action="write")
126         write(unit1,'(3f15.7)') (pxyz(ix)/dx,ix=1,3)
127         call io_close(unit1)
128
129      endif
130
131      end
132