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