1!
2! Copyright (C) 2005 PWSCF group
3! This file is distributed under the terms of the
4! GNU General Public License. See the file `License'
5! in the root directory of the present distribution,
6! or http://www.gnu.org/copyleft/gpl.txt .
7!
8!
9!----------------------------------------------------------------------------
10SUBROUTINE report_mag
11      !------------------------------------------------------------------------
12      !! This subroutine prints out information about the local magnetization
13      !! and/or charge, integrated around the atomic positions at points which
14      !! are calculated in make_pointlists.
15      !
16      USE kinds,            ONLY: DP
17      USE ions_base,        ONLY: nat, tau, ityp
18      USE io_global,        ONLY: stdout
19      use constants,        ONLY: pi
20      USE scf,              ONLY: rho
21      USE noncollin_module, ONLY: noncolin, mcons, i_cons
22      USE lsda_mod,         ONLY: nspin
23      !
24      IMPLICIT NONE
25      !
26      REAL(DP) :: theta, phi, norm, norm1
27      INTEGER :: ipol, iat
28      REAL(DP) :: r1_loc(nat), m1_loc(nspin-1,nat)
29      !
30      ! get_local integrates on the previously determined points
31      !
32      CALL get_locals( r1_loc, m1_loc, rho%of_r )
33      !
34      IF (nspin == 2) THEN
35         WRITE( stdout, * )
36         WRITE( stdout, '(5X,"Magnetic moment per site:")' )
37         DO iat = 1, nat
38            WRITE(stdout,1020) iat, r1_loc(iat), m1_loc(1,iat), mcons(1,ityp(iat))
39         END DO
40         !
41      ELSE IF (noncolin) THEN
42         DO iat = 1, nat
43            !
44            ! norm is the length of the magnetic moment vector
45            !
46            norm = DSQRT( m1_loc(1,iat)**2+m1_loc(2,iat)**2+m1_loc(3,iat)**2 )
47            !
48            ! norm1 is the length of the projection of the mm vector into
49            ! the xy plane
50            !
51            norm1 = DSQRT( m1_loc(1,iat)**2+m1_loc(2,iat)**2 )
52            !
53            ! calculate the polar angles of the magnetic moment
54            !
55            IF (norm > 1.d-10) THEN
56               theta = ACOS(m1_loc(3,iat)/norm)
57               IF (norm1 > 1.d-10) THEN
58                  phi = ACOS(m1_loc(1,iat)/norm1)
59                  IF (m1_loc(2,iat) < 0.d0) phi = - phi
60               ELSE
61                  phi = 2.d0*pi
62               ENDIF
63            ELSE
64               theta = 2.d0*pi
65               phi = 2.d0*pi
66            ENDIF
67            !
68            ! go to degrees
69            theta = theta*180.d0/pi
70            phi = phi*180.d0/pi
71            !
72            WRITE( stdout,1010)
73            WRITE( stdout,1011) iat,(tau(ipol,iat),ipol=1,3)
74            WRITE( stdout,1014) r1_loc (iat)
75            WRITE( stdout,1012) (m1_loc(ipol,iat),ipol=1,3)
76            WRITE( stdout,1018) (m1_loc(ipol,iat)/r1_loc(iat),ipol=1,3)
77            WRITE( stdout,1013) norm,theta,phi
78            IF (i_cons==1) THEN
79               WRITE( stdout,1015) (mcons(ipol,ityp(iat)),ipol=1,3)
80            ELSEIF (i_cons==2) THEN
81               WRITE( stdout,1017) 180.d0 * ACOS(mcons(3,ityp(iat)))/pi
82            ENDIF
83            WRITE( stdout,1010)
84         ENDDO
85         !
86      END IF
87      !
88      !
89 1010 FORMAT (/,1x,78('='))
90 1011 FORMAT (5x,'atom number ',i4,' relative position : ',3f9.4)
91 1012 FORMAT (5x,'magnetization :      ',3f12.6)
92 1013 FORMAT (5x,'polar coord.: r, theta, phi [deg] : ',3f12.6)
93 1014 FORMAT (5x,'charge : ',f12.6)
94 1018 FORMAT (5x,'magnetization/charge:',3f12.6)
95 1015 FORMAT (5x,'constrained moment : ',3f12.6)
96 1017 FORMAT (5x,'constrained theta [deg] : ',f12.6)
97 1020 FORMAT (5x,'atom: ',i4,4X,'charge: ',F9.4,4X,'magn: ',F9.4,4X,'constr: ',f9.4)
98      !
99END SUBROUTINE report_mag
100