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      module m_born_charge
9
10      private
11      public :: born_charge
12
13      CONTAINS
14
15      subroutine born_charge()
16      use precision
17      use siesta_options, only:dx
18      use sparse_matrices
19      use parallel,       only: Node, Nodes, IOnode
20      USE m_ksvinit,      only: nkpol, kpol, wgthpol
21      use siesta_geom,    only: na_u, na_s, xa, scell, ucell, isa
22      use atomlist,       only: no_u, indxuo, iphorb, iaorb, lasto
23      use atomlist,       only: rmaxo, no_s, no_l
24      use m_spin,         only: nspin, SPpol, NonCol, SpOrb ! CC RC Added
25      use m_ksv,          only: KSV_pol
26      use siesta_geom,    only: shape
27#ifdef MPI
28      use m_mpi_utils, only: globalize_sum
29#endif
30
31      implicit none
32
33      real(dp) :: polxyz(3,nspin)  ! Automatic, small
34      real(dp) :: polR(3,nspin)    ! Automatic, small
35
36      real(dp) :: qspin(4)         ! Local variable
37
38      integer  :: j, ind, io, ispin
39#ifdef MPI
40      real(dp):: qtmp(4)   ! Temporary for globalized spin charge
41#endif
42
43!     Computes Born-effective charges by finite differences,
44!     using the "force constant calculation" mode of operation.
45!     It is called only for those steps (even) in which the dx
46!     is positive.
47!
48      if (NonCol .or. SpOrb) then
49       if ( IONode ) then
50          write(*,'(/,a)')
51     .         'born_charge: implemented only for spin-polarized '
52          write(*,'(a)')
53     .         'born_charge: or paramagnetic calculations'
54       end if
55
56       return
57
58      endif
59
60      if (nkpol.lt.1) then
61         if (IOnode) write(6,'(/,a,f12.6)')
62     .        'siesta: specify polarization grid for BC calculation'
63         if (IOnode) write(6,'(a,f12.6)')
64     .        'siesta: The Born charge matrix will not be calculated'
65         RETURN
66      endif
67      if (IOnode) write(6,'(/,a,f12.6)')
68     .     'siesta: Calculating polarization. '
69
70!     Find total population of spin up and down
71      if (nspin .ge. 2) then
72         do ispin = 1,nspin
73            qspin(ispin) = 0.0_dp
74            do io = 1,no_l
75               do j = 1,numh(io)
76                  ind = listhptr(io) + j
77                  qspin(ispin) = qspin(ispin)
78     .                 + Dscf(ind,ispin)*S(ind)
79               enddo
80            enddo
81         enddo
82#ifdef MPI
83!     Global reduction of spin components
84         call globalize_sum(qspin(1:nspin),qtmp(1:nspin))
85         qspin(1:nspin) = qtmp(1:nspin)
86#endif
87      endif
88      if (nkpol.gt.0) then
89         ! Note use of xa instead of xa_last, since this routine
90         ! is called at every geometry step, before moving the atoms.
91         !
92         call KSV_pol(na_u, na_s, xa, rmaxo, scell, ucell,
93     .        no_u, no_l, no_s, nspin, qspin,
94     .        maxnh, nkpol, numh, listhptr, listh,
95     .        H, S, xijo, indxuo, isa, iphorb,
96     .        iaorb, lasto, shape,
97     .        nkpol,kpol,wgthpol, polR, polxyz)
98      endif
99      if (nkpol.gt.0.and.IOnode) then
100         call obc( polR, ucell, dx, nspin)
101      endif
102
103      end subroutine born_charge
104
105      end module m_born_charge
106