1!
2! Copyright (C) 2001-2012 Quantum ESPRESSO 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!--------------------------------------------------------------------
9program dynmat
10!--------------------------------------------------------------------
11!
12!  This program
13!  - reads a dynamical matrix file produced by the phonon code
14!  - adds the nonanalytical part (if Z* and epsilon are read from file),
15!    applies the chosen Acoustic Sum Rule (if q=0)
16!  - diagonalise the dynamical matrix
17!  - calculates IR and Raman cross sections (if Z* and Raman tensors
18!    are read from file, respectively)
19!  - writes the results to files, both for inspection and for plotting
20!
21!  Input data (namelist "input")
22!
23!  fildyn  character input file containing the dynamical matrix
24!                    (default: fildyn='matdyn')
25!  q(3)      real    calculate LO modes (add nonanalytic terms) along
26!                    the direction q (cartesian axis, default: q=(0,0,0) )
27!  amass(nt) real    mass for atom type nt, amu
28!                    (default: amass is read from file fildyn)
29!  asr   character   indicates the type of Acoustic Sum Rule imposed
30!                     - 'no': no Acoustic Sum Rules imposed (default)
31!                     - 'simple':  previous implementation of the asr used
32!                     (3 translational asr imposed by correction of
33!                     the diagonal elements of the dynamical matrix)
34!                     - 'crystal': 3 translational asr imposed by optimized
35!                     correction of the dyn. matrix (projection).
36!                     - 'one-dim': 3 translational asr + 1 rotational asr
37!                     imposed by optimized correction of the dyn. mat. (the
38!                     rotation axis is the direction of periodicity; it
39!                     will work only if this axis considered is one of
40!                     the cartesian axis).
41!                     - 'zero-dim': 3 translational asr + 3 rotational asr
42!                     imposed by optimized correction of the dyn. mat.
43!                     Note that in certain cases, not all the rotational asr
44!                     can be applied (e.g. if there are only 2 atoms in a
45!                     molecule or if all the atoms are aligned, etc.).
46!                     In these cases the supplementary asr are cancelled
47!                     during the orthonormalization procedure (see below).
48!                     Finally, in all cases except 'no' a simple correction
49!                     on the effective charges is performed (same as in the
50!                     previous implementation).
51!  axis    integer    indicates the rotation axis for a 1D system
52!                     (1=Ox, 2=Oy, 3=Oz ; default =3)
53!  lperm   logical    .true. to calculate Gamma-point mode contributions to
54!                     dielectric permittivity tensor
55!                     (default: lperm=.false.)
56!  lplasma logical    .true. to calculate Gamma-point mode effective plasma
57!                     frequencies, automatically triggers lperm = .true.
58!                     (default: lplasma=.false.)
59!  filout  character output file containing phonon frequencies and normalized
60!                    phonon displacements (i.e. eigenvectors divided by the
61!                    square root of the mass and then normalized; they are
62!                    not orthogonal)
63!                    (default: filout='dynmat.out')
64!  fileig  character output file containing phonon frequencies and eigenvectors
65!                    of the dynamical matrix (they are orthogonal)
66!                    (default: fileig=' ')
67!  filmol  character as above, in a format suitable for 'molden'
68!                    (default: filmol='dynmat.mold')
69!  filxsf  character as above, in axsf format suitable for xcrysden
70!                    (default: filxsf='dynmat.axsf')
71!  loto_2d logical set to .true. to activate two-dimensional treatment of LO-TO splitting.
72!
73  USE kinds, ONLY: DP
74  USE mp,         ONLY : mp_bcast
75  USE mp_global,  ONLY : mp_startup, mp_global_end
76  USE mp_world,   ONLY : world_comm
77  USE io_global,  ONLY : ionode, ionode_id, stdout
78  USE environment, ONLY : environment_start, environment_end
79  USE io_dyn_mat,  ONLY : read_dyn_mat_param, read_dyn_mat_header, &
80                         read_dyn_mat, read_dyn_mat_tail
81  USE constants,   ONLY : amu_ry
82  USE dynamical
83  USE rigid,       ONLY: dyndiag, nonanal
84  !
85  implicit none
86  integer, parameter :: ntypx = 10
87  character(len=256):: fildyn, filout, filmol, filxsf, fileig
88  character(len=3) :: atm(ntypx)
89  character(len=10) :: asr
90  logical :: lread, gamma, loto_2d
91  complex(DP), allocatable :: z(:,:)
92  real(DP) :: amass(ntypx), amass_(ntypx), eps0(3,3), a0, omega, &
93       at(3,3), bg(3,3), q(3), q_(3)
94  real(DP), allocatable :: w2(:)
95  integer :: nat, na, nt, ntyp, iout, axis, nspin_mag, ios
96  real(DP) :: celldm(6)
97  logical :: xmldyn, lrigid, lraman, lperm, lplasma
98  logical, external :: has_xml
99  integer :: ibrav, nqs
100  integer, allocatable :: itau(:)
101  namelist /input/ amass, asr, axis, fildyn, filout, filmol, filxsf, &
102                   fileig, lperm, lplasma, q, loto_2d
103  !
104  ! code is parallel-compatible but not parallel
105  !
106  CALL mp_startup()
107  CALL environment_start('DYNMAT')
108  !
109  IF (ionode) CALL input_from_file ( )
110  !
111  asr  = 'no'
112  axis = 3
113  fildyn='matdyn'
114  filout='dynmat.out'
115  filmol='dynmat.mold'
116  filxsf='dynmat.axsf'
117  fileig=' '
118  amass(:)=0.0d0
119  q(:)=0.0d0
120  lperm=.false.
121  lplasma=.false.
122  loto_2d=.false.
123  !
124  IF (ionode) read (5,input, iostat=ios)
125  CALL mp_bcast(ios, ionode_id, world_comm)
126  CALL errore('dynmat', 'reading input namelist', ABS(ios))
127  !
128  CALL mp_bcast(asr,ionode_id, world_comm)
129  CALL mp_bcast(axis,ionode_id, world_comm)
130  CALL mp_bcast(amass,ionode_id, world_comm)
131  CALL mp_bcast(fildyn,ionode_id, world_comm)
132  CALL mp_bcast(filout,ionode_id, world_comm)
133  CALL mp_bcast(filmol,ionode_id, world_comm)
134  CALL mp_bcast(fileig,ionode_id, world_comm)
135  CALL mp_bcast(filxsf,ionode_id, world_comm)
136  CALL mp_bcast(q,ionode_id, world_comm)
137  !
138  IF (ionode) inquire(file=fildyn,exist=lread)
139  CALL mp_bcast(lread, ionode_id, world_comm)
140  IF (lread) THEN
141     IF (ionode) WRITE(6,'(/5x,a,a)') 'Reading Dynamical Matrix from file '&
142                                     , TRIM(fildyn)
143  ELSE
144     CALL errore('dynmat', 'File '//TRIM(fildyn)//' not found', 1)
145  END IF
146  !
147  ntyp = ntypx ! avoids spurious out-of-bound errors
148  xmldyn=has_xml(fildyn)
149  IF (xmldyn) THEN
150     CALL read_dyn_mat_param(fildyn,ntyp,nat)
151     ALLOCATE (m_loc(3,nat))
152     ALLOCATE (tau(3,nat))
153     ALLOCATE (ityp(nat))
154     ALLOCATE (zstar(3,3,nat))
155     ALLOCATE (dchi_dtau(3,3,3,nat) )
156     CALL read_dyn_mat_header(ntyp, nat, ibrav, nspin_mag, &
157             celldm, at, bg, omega, atm, amass_, tau, ityp, &
158             m_loc, nqs, lrigid, eps0, zstar, lraman, dchi_dtau)
159     IF (nqs /= 1) CALL errore('dynmat','only q=0 matrix allowed',1)
160     a0=celldm(1) ! define alat
161     ALLOCATE (dyn(3,3,nat,nat) )
162     CALL read_dyn_mat(nat,1,q_,dyn(:,:,:,:))
163     CALL read_dyn_mat_tail(nat)
164     IF(asr.ne.'no') THEN
165         CALL set_asr ( asr, axis, nat, tau, dyn, zstar )
166     END IF
167     IF (ionode) THEN
168        DO nt=1, ntyp
169           IF (amass(nt) <= 0.0d0) amass(nt)=amass_(nt)
170        END DO
171     END IF
172  ELSE
173     IF (ionode) THEN
174        CALL readmat2 ( fildyn, asr, axis, nat, ntyp, atm, a0, &
175                        at, omega, amass_, eps0, q_ )
176        DO nt=1, ntyp
177           IF (amass(nt) <= 0.0d0) amass(nt)=amass_(nt)/amu_ry
178        END DO
179     END IF
180  ENDIF
181  !
182  IF (ionode) THEN
183     !
184     ! from now on, execute on a single processor
185     !
186     gamma = ( abs( q_(1)**2+q_(2)**2+q_(3)**2 ) < 1.0d-8 )
187     !
188     IF (gamma .and. .not.loto_2d) THEN
189        ALLOCATE (itau(nat))
190        DO na=1,nat
191           itau(na)=na
192        END DO
193        CALL nonanal ( nat, nat, itau, eps0, q, zstar, omega, dyn )
194        DEALLOCATE (itau)
195     END IF
196     !
197     ALLOCATE ( z(3*nat,3*nat), w2(3*nat) )
198     CALL dyndiag(nat,ntyp,amass,ityp,dyn,w2,z)
199     !
200     IF (filout.eq.' ') then
201        iout=6
202     ELSE
203        iout=4
204        OPEN (unit=iout,file=filout,status='unknown',form='formatted')
205     END IF
206     CALL writemodes(nat,q_,w2,z,iout)
207     IF(iout .ne. 6) close(unit=iout)
208     IF (fileig .ne. ' ') THEN
209       OPEN (unit=15,file=TRIM(fileig),status='unknown',form='formatted')
210       CALL write_eigenvectors (nat,ntyp,amass,ityp,q_,w2,z,15)
211       CLOSE (unit=15)
212     ENDIF
213     CALL writemolden (filmol, gamma, nat, atm, a0, tau, ityp, w2, z)
214     CALL writexsf (filxsf, gamma, nat, atm, a0, at, tau, ityp, z)
215     IF (gamma) THEN
216        CALL RamanIR (nat, omega, w2, z, zstar, eps0, dchi_dtau)
217        IF (lperm .OR. lplasma) THEN
218            CALL polar_mode_permittivity(nat,eps0,z,zstar,w2,omega, &
219                                         lplasma)
220            IF ( ABS( q(1)**2+q(2)**2+q(3)**2 ) > 1.0d-8 ) &
221               WRITE(6,'(5x,a)') 'BEWARE: phonon contribution to &
222               & permittivity computed with TO-LO splitting'
223        ENDIF
224     ENDIF
225  ENDIF
226  !
227  IF (xmldyn) THEN
228     DEALLOCATE (m_loc)
229     DEALLOCATE (tau)
230     DEALLOCATE (ityp)
231     DEALLOCATE (zstar)
232     DEALLOCATE (dchi_dtau)
233     DEALLOCATE (dyn)
234  ENDIF
235  !
236  CALL environment_end('DYNMAT')
237  !
238  CALL mp_global_end()
239  !
240end program dynmat
241