1! Copyright (C) 2001-2012 Quantum ESPRESSO group 2! This file is distributed under the terms of the 3! GNU General Public License. See the file `License' 4! in the root directory of the present distribution, 5! or http://www.gnu.org/copyleft/gpl.txt . 6! 7! 8Module ifconstants 9 ! 10 ! All variables read from file that need dynamical allocation 11 ! 12 USE kinds, ONLY: DP 13 REAL(DP), ALLOCATABLE :: frc(:,:,:,:,:,:,:), tau_blk(:,:), zeu(:,:,:), & 14 m_loc(:,:) 15 ! frc : interatomic force constants in real space 16 ! tau_blk : atomic positions for the original cell 17 ! zeu : effective charges for the original cell 18 ! m_loc: the magnetic moments of each atom 19 INTEGER, ALLOCATABLE :: ityp_blk(:) 20 ! ityp_blk : atomic types for each atom of the original cell 21 ! 22 CHARACTER(LEN=3), ALLOCATABLE :: atm(:) 23end Module ifconstants 24! 25!--------------------------------------------------------------------- 26PROGRAM matdyn 27 !----------------------------------------------------------------------- 28 ! this program calculates the phonon frequencies for a list of generic 29 ! q vectors starting from the interatomic force constants generated 30 ! from the dynamical matrices as written by DFPT phonon code through 31 ! the companion program q2r 32 ! 33 ! matdyn can generate a supercell of the original cell for mass 34 ! approximation calculation. If supercell data are not specified 35 ! in input, the unit cell, lattice vectors, atom types and positions 36 ! are read from the force constant file 37 ! 38 ! Input cards: namelist &input 39 ! flfrc file produced by q2r containing force constants (needed) 40 ! It is the same as in the input of q2r.x (+ the .xml extension 41 ! if the dynamical matrices produced by ph.x were in xml 42 ! format). No default value: must be specified. 43 ! asr (character) indicates the type of Acoustic Sum Rule imposed 44 ! - 'no': no Acoustic Sum Rules imposed (default) 45 ! - 'simple': previous implementation of the asr used 46 ! (3 translational asr imposed by correction of 47 ! the diagonal elements of the force constants matrix) 48 ! - 'crystal': 3 translational asr imposed by optimized 49 ! correction of the force constants (projection). 50 ! - 'one-dim': 3 translational asr + 1 rotational asr 51 ! imposed by optimized correction of the force constants 52 ! (the rotation axis is the direction of periodicity; 53 ! it will work only if this axis considered is one of 54 ! the cartesian axis). 55 ! - 'zero-dim': 3 translational asr + 3 rotational asr 56 ! imposed by optimized correction of the force constants 57 ! Note that in certain cases, not all the rotational asr 58 ! can be applied (e.g. if there are only 2 atoms in a 59 ! molecule or if all the atoms are aligned, etc.). 60 ! In these cases the supplementary asr are cancelled 61 ! during the orthonormalization procedure (see below). 62 ! dos if .true. calculate phonon Density of States (DOS) 63 ! using tetrahedra and a uniform q-point grid (see below) 64 ! NB: may not work properly in noncubic materials 65 ! if .false. calculate phonon bands from the list of q-points 66 ! supplied in input (default) 67 ! nk1,nk2,nk3 uniform q-point grid for DOS calculation (includes q=0) 68 ! (must be specified if dos=.true., ignored otherwise) 69 ! deltaE energy step, in cm^(-1), for DOS calculation: from min 70 ! to max phonon energy (default: 1 cm^(-1) if ndos, see 71 ! below, is not specified) 72 ! ndos number of energy steps for DOS calculations 73 ! (default: calculated from deltaE if not specified) 74 ! fldos output file for dos (default: 'matdyn.dos') 75 ! the dos is in states/cm(-1) plotted vs omega in cm(-1) 76 ! and is normalised to 3*nat, i.e. the number of phonons 77 ! flfrq output file for frequencies (default: 'matdyn.freq') 78 ! flvec output file for normalized phonon displacements 79 ! (default: 'matdyn.modes'). The normalized phonon displacements 80 ! are the eigenvectors divided by the square root of the mass, 81 ! then normalized. As such they are not orthogonal. 82 ! fleig output file for phonon eigenvectors (default: 'matdyn.eig') 83 ! The phonon eigenvectors are the eigenvectors of the dynamical 84 ! matrix. They are orthogonal. 85 ! fldyn output file for dynamical matrix (default: ' ' i.e. not written) 86 ! at supercell lattice vectors - must form a superlattice of the 87 ! original lattice (default: use original cell) 88 ! l1,l2,l3 supercell lattice vectors are original cell vectors times 89 ! l1, l2, l3 respectively (default: 1, ignored if at specified) 90 ! ntyp number of atom types in the supercell (default: ntyp of the 91 ! original cell) 92 ! amass masses of atoms in the supercell (a.m.u.), one per atom type 93 ! (default: use masses read from file flfrc) 94 ! readtau read atomic positions of the supercell from input 95 ! (used to specify different masses) (default: .false.) 96 ! fltau write atomic positions of the supercell to file "fltau" 97 ! (default: fltau=' ', do not write) 98 ! la2F if .true. interpolates also the el-ph coefficients. 99 ! q_in_band_form if .true. the q points are given in band form: 100 ! Only the first and last point of one or more lines 101 ! are given. See below. (default: .false.). 102 ! q_in_cryst_coord if .true. input q points are in crystalline 103 ! coordinates (default: .false.) 104 ! eigen_similarity: use similarity of the displacements to order 105 ! frequencies (default: .false.) 106 ! NB: You cannot use this option with the symmetry 107 ! analysis of the modes. 108 ! fd (logical) if .t. the ifc come from the finite displacement calculation 109 ! na_ifc (logical) add non analitic contributions to the interatomic force 110 ! constants if finite displacement method is used (as in Wang et al. 111 ! Phys. Rev. B 85, 224303 (2012)) [to be used in conjunction with fd.x] 112 ! nosym if .true., no symmetry and no time reversal are imposed 113 ! loto_2d set to .true. to activate two-dimensional treatment of LO-TO splitting. 114 ! loto_disable (logical) if .true. do not apply LO-TO splitting for q=0 115 ! (default: .false.) 116 ! 117 ! if (readtau) atom types and positions in the supercell follow: 118 ! (tau(i,na),i=1,3), ityp(na) 119 ! IF (q_in_band_form.and..not.dos) THEN 120 ! nq ! number of q points 121 ! (q(i,n),i=1,3), nptq nptq is the number of points between this point 122 ! and the next. These points are automatically 123 ! generated. the q points are given in Cartesian 124 ! coordinates, 2pi/a units (a=lattice parameters) 125 ! ELSE, if (.not.dos) : 126 ! nq number of q-points 127 ! (q(i,n), i=1,3) nq q-points in cartesian coordinates, 2pi/a units 128 ! If q = 0, the direction qhat (q=>0) for the non-analytic part 129 ! is extracted from the sequence of q-points as follows: 130 ! qhat = q(n) - q(n-1) or qhat = q(n) - q(n+1) 131 ! depending on which one is available and nonzero. 132 ! For low-symmetry crystals, specify twice q = 0 in the list 133 ! if you want to have q = 0 results for two different directions 134 ! 135 USE kinds, ONLY : DP 136 USE mp, ONLY : mp_bcast 137 USE mp_world, ONLY : world_comm 138 USE mp_global, ONLY : mp_startup, mp_global_end 139 USE environment, ONLY : environment_start, environment_end 140 USE io_global, ONLY : ionode, ionode_id, stdout 141 USE io_dyn_mat, ONLY : read_dyn_mat_param, read_dyn_mat_header, & 142 read_ifc_param, read_ifc 143 USE cell_base, ONLY : at, bg, celldm 144 USE constants, ONLY : RY_TO_THZ, RY_TO_CMM1, amu_ry 145 USE symm_base, ONLY : set_sym 146 USE rap_point_group, ONLY : code_group 147 USE bz_form, ONLY : transform_label_coord 148 USE parser, ONLY : read_line 149 USE rigid, ONLY : dyndiag, nonanal, nonanal_ifc 150 USE el_phon, ONLY : el_ph_nsigma 151 152 USE ifconstants, ONLY : frc, atm, zeu, tau_blk, ityp_blk, m_loc 153 ! 154 IMPLICIT NONE 155 ! 156 INTEGER :: gid 157 ! 158 ! variables *_blk refer to the original cell, other variables 159 ! to the (super)cell (which may coincide with the original cell) 160 ! 161 INTEGER:: nax, nax_blk 162 INTEGER, PARAMETER:: ntypx=10, nrwsx=200 163 REAL(DP), PARAMETER :: eps=1.0d-6 164 INTEGER :: nr1, nr2, nr3, nsc, nk1, nk2, nk3, ibrav 165 CHARACTER(LEN=256) :: flfrc, flfrq, flvec, fltau, fldos, filename, fldyn, & 166 fleig, fildyn, fildyn_prefix 167 CHARACTER(LEN=10) :: asr 168 LOGICAL :: dos, has_zstar, q_in_cryst_coord, eigen_similarity, loto_disable 169 COMPLEX(DP), ALLOCATABLE :: dyn(:,:,:,:), dyn_blk(:,:,:,:), frc_ifc(:,:,:,:) 170 COMPLEX(DP), ALLOCATABLE :: z(:,:) 171 REAL(DP), ALLOCATABLE:: tau(:,:), q(:,:), w2(:,:), freq(:,:), wq(:), & 172 dynq(:,:,:), DOSofE(:) 173 INTEGER, ALLOCATABLE:: ityp(:), itau_blk(:) 174 REAL(DP) :: omega,alat, &! cell parameters and volume 175 at_blk(3,3), bg_blk(3,3), &! original cell 176 omega_blk, &! original cell volume 177 epsil(3,3), &! dielectric tensor 178 amass(ntypx), &! atomic masses 179 amass_blk(ntypx), &! original atomic masses 180 atws(3,3), &! lattice vector for WS initialization 181 rws(0:3,nrwsx) ! nearest neighbor list, rws(0,*) = norm^2 182 ! 183 INTEGER :: nat, nat_blk, ntyp, ntyp_blk, & 184 l1, l2, l3, &! supercell dimensions 185 nrws, &! number of nearest neighbor 186 code_group_old 187 188 INTEGER :: nspin_mag, nqs, ios 189 ! 190 LOGICAL :: readtau, la2F, xmlifc, lo_to_split, na_ifc, fd, nosym, loto_2d 191 ! 192 REAL(DP) :: qhat(3), qh, DeltaE, Emin=0._dp, Emax, E, qq 193 REAL(DP) :: delta, pathL 194 REAL(DP), ALLOCATABLE :: xqaux(:,:) 195 INTEGER, ALLOCATABLE :: nqb(:) 196 INTEGER :: n, i, j, it, nq, nqx, na, nb, ndos, iout, nqtot, iout_dyn, iout_eig 197 LOGICAL, EXTERNAL :: has_xml 198 INTEGER, ALLOCATABLE :: num_rap_mode(:,:) 199 LOGICAL, ALLOCATABLE :: high_sym(:) 200 LOGICAL :: q_in_band_form 201 ! .... variables for band plotting based on similarity of eigenvalues 202 COMPLEX(DP), ALLOCATABLE :: tmp_z(:,:) 203 REAL(DP), ALLOCATABLE :: abs_similarity(:,:), tmp_w2(:) 204 COMPLEX(DP), ALLOCATABLE :: f_of_q(:,:,:,:) 205 INTEGER :: location(1), isig 206 CHARACTER(LEN=6) :: int_to_char 207 LOGICAL, ALLOCATABLE :: mask(:) 208 INTEGER :: npk_label, nch 209 CHARACTER(LEN=3), ALLOCATABLE :: letter(:) 210 INTEGER, ALLOCATABLE :: label_list(:) 211 LOGICAL :: tend, terr 212 CHARACTER(LEN=256) :: input_line, buffer 213 CHARACTER(LEN=10) :: point_label_type 214 CHARACTER(len=80) :: k_points = 'tpiba' 215 ! 216 REAL(DP), external :: dos_gam 217 ! 218 NAMELIST /input/ flfrc, amass, asr, flfrq, flvec, fleig, at, dos, & 219 & fldos, nk1, nk2, nk3, l1, l2, l3, ntyp, readtau, fltau, & 220 & la2F, ndos, DeltaE, q_in_band_form, q_in_cryst_coord, & 221 & eigen_similarity, fldyn, na_ifc, fd, point_label_type, & 222 & nosym, loto_2d, fildyn, fildyn_prefix, el_ph_nsigma, & 223 & loto_disable 224 ! 225 CALL mp_startup() 226 CALL environment_start('MATDYN') 227 ! 228 IF (ionode) CALL input_from_file ( ) 229 ! 230 ! ... all calculations are done by the first cpu 231 ! 232 ! set namelist default 233 ! 234 dos = .FALSE. 235 deltaE = 1.0d0 236 ndos = 1 237 nk1 = 0 238 nk2 = 0 239 nk3 = 0 240 asr ='no' 241 readtau=.FALSE. 242 flfrc=' ' 243 fldos='matdyn.dos' 244 flfrq='matdyn.freq' 245 flvec='matdyn.modes' 246 fleig=' ' 247 fldyn=' ' 248 fltau=' ' 249 fildyn = ' ' 250 fildyn_prefix = ' ' 251 amass(:) =0.d0 252 amass_blk(:) =0.d0 253 at(:,:) = 0.d0 254 ntyp = 0 255 l1=1 256 l2=1 257 l3=1 258 la2F=.false. 259 q_in_band_form=.FALSE. 260 eigen_similarity=.FALSE. 261 q_in_cryst_coord = .FALSE. 262 na_ifc=.FALSE. 263 fd=.FALSE. 264 point_label_type='SC' 265 nosym = .false. 266 loto_2d=.false. 267 el_ph_nsigma=10 268 loto_disable = .false. 269 ! 270 ! 271 IF (ionode) READ (5,input,IOSTAT=ios) 272 CALL mp_bcast(ios, ionode_id, world_comm) 273 CALL errore('matdyn', 'reading input namelist', ABS(ios)) 274 CALL mp_bcast(dos,ionode_id, world_comm) 275 CALL mp_bcast(deltae,ionode_id, world_comm) 276 CALL mp_bcast(ndos,ionode_id, world_comm) 277 CALL mp_bcast(nk1,ionode_id, world_comm) 278 CALL mp_bcast(nk2,ionode_id, world_comm) 279 CALL mp_bcast(nk3,ionode_id, world_comm) 280 CALL mp_bcast(asr,ionode_id, world_comm) 281 CALL mp_bcast(readtau,ionode_id, world_comm) 282 CALL mp_bcast(flfrc,ionode_id, world_comm) 283 CALL mp_bcast(fldos,ionode_id, world_comm) 284 CALL mp_bcast(flfrq,ionode_id, world_comm) 285 CALL mp_bcast(flvec,ionode_id, world_comm) 286 CALL mp_bcast(fleig,ionode_id, world_comm) 287 CALL mp_bcast(fldyn,ionode_id, world_comm) 288 CALL mp_bcast(fltau,ionode_id, world_comm) 289 CALL mp_bcast(fildyn,ionode_id, world_comm) 290 CALL mp_bcast(fildyn_prefix,ionode_id, world_comm) 291 CALL mp_bcast(amass,ionode_id, world_comm) 292 CALL mp_bcast(amass_blk,ionode_id, world_comm) 293 CALL mp_bcast(at,ionode_id, world_comm) 294 CALL mp_bcast(ntyp,ionode_id, world_comm) 295 CALL mp_bcast(l1,ionode_id, world_comm) 296 CALL mp_bcast(l2,ionode_id, world_comm) 297 CALL mp_bcast(l3,ionode_id, world_comm) 298 CALL mp_bcast(na_ifc,ionode_id, world_comm) 299 CALL mp_bcast(fd,ionode_id, world_comm) 300 CALL mp_bcast(la2F,ionode_id, world_comm) 301 CALL mp_bcast(q_in_band_form,ionode_id, world_comm) 302 CALL mp_bcast(eigen_similarity,ionode_id, world_comm) 303 CALL mp_bcast(q_in_cryst_coord,ionode_id, world_comm) 304 CALL mp_bcast(point_label_type,ionode_id, world_comm) 305 CALL mp_bcast(loto_2d,ionode_id, world_comm) 306 CALL mp_bcast(loto_disable,ionode_id, world_comm) 307 CALL mp_bcast(el_ph_nsigma,ionode_id, world_comm) 308 ! 309 IF (loto_2d .AND. loto_disable) CALL errore('matdyn', & 310 'loto_2d and loto_disable cannot be both true', 1) 311 ! 312 ! read force constants 313 ! 314 IF ( trim( fildyn ) /= ' ' ) THEN 315 IF (ionode) THEN 316 WRITE(stdout, *) 317 WRITE(stdout, '(4x,a)') ' fildyn has been provided, running q2r...' 318 END IF 319 CALL do_q2r(fildyn, flfrc, fildyn_prefix, asr, la2F, loto_2d) 320 END IF 321 ! 322 ntyp_blk = ntypx ! avoids fake out-of-bound error 323 xmlifc=has_xml(flfrc) 324 IF (xmlifc) THEN 325 CALL read_dyn_mat_param(flfrc,ntyp_blk,nat_blk) 326 ALLOCATE (m_loc(3,nat_blk)) 327 ALLOCATE (tau_blk(3,nat_blk)) 328 ALLOCATE (ityp_blk(nat_blk)) 329 ALLOCATE (atm(ntyp_blk)) 330 ALLOCATE (zeu(3,3,nat_blk)) 331 CALL read_dyn_mat_header(ntyp_blk, nat_blk, ibrav, nspin_mag, & 332 celldm, at_blk, bg_blk, omega_blk, atm, amass_blk, & 333 tau_blk, ityp_blk, m_loc, nqs, has_zstar, epsil, zeu ) 334 alat=celldm(1) 335 call volume(alat,at_blk(1,1),at_blk(1,2),at_blk(1,3),omega_blk) 336 CALL read_ifc_param(nr1,nr2,nr3) 337 ALLOCATE(frc(nr1,nr2,nr3,3,3,nat_blk,nat_blk)) 338 CALL read_ifc(nr1,nr2,nr3,nat_blk,frc) 339 ELSE 340 CALL readfc ( flfrc, nr1, nr2, nr3, epsil, nat_blk, & 341 ibrav, alat, at_blk, ntyp_blk, & 342 amass_blk, omega_blk, has_zstar) 343 ENDIF 344 ! 345 CALL recips ( at_blk(1,1),at_blk(1,2),at_blk(1,3), & 346 bg_blk(1,1),bg_blk(1,2),bg_blk(1,3) ) 347 ! 348 ! set up (super)cell 349 ! 350 if (ntyp < 0) then 351 call errore ('matdyn','wrong ntyp ', abs(ntyp)) 352 else if (ntyp == 0) then 353 ntyp=ntyp_blk 354 end if 355 ! 356 ! masses (for mass approximation) 357 ! 358 DO it=1,ntyp 359 IF (amass(it) < 0.d0) THEN 360 CALL errore ('matdyn','wrong mass in the namelist',it) 361 ELSE IF (amass(it) == 0.d0) THEN 362 IF (it.LE.ntyp_blk) THEN 363 WRITE (stdout,'(a,i3,a,a)') ' mass for atomic type ',it, & 364 & ' not given; uses mass from file ',flfrc 365 amass(it) = amass_blk(it) 366 ELSE 367 CALL errore ('matdyn','missing mass in the namelist',it) 368 END IF 369 END IF 370 END DO 371 ! 372 ! lattice vectors 373 ! 374 IF (SUM(ABS(at(:,:))) == 0.d0) THEN 375 IF (l1.LE.0 .OR. l2.LE.0 .OR. l3.LE.0) CALL & 376 & errore ('matdyn',' wrong l1,l2 or l3',1) 377 at(:,1) = at_blk(:,1)*DBLE(l1) 378 at(:,2) = at_blk(:,2)*DBLE(l2) 379 at(:,3) = at_blk(:,3)*DBLE(l3) 380 END IF 381 ! 382 CALL check_at(at,bg_blk,alat,omega) 383 ! 384 ! the supercell contains "nsc" times the original unit cell 385 ! 386 nsc = NINT(omega/omega_blk) 387 IF (ABS(omega/omega_blk-nsc) > eps) & 388 CALL errore ('matdyn', 'volume ratio not integer', 1) 389 ! 390 ! read/generate atomic positions of the (super)cell 391 ! 392 nat = nat_blk * nsc 393 nax = nat 394 nax_blk = nat_blk 395 ! 396 ALLOCATE ( tau (3, nat), ityp(nat), itau_blk(nat) ) 397 ! 398 IF (readtau) THEN 399 CALL read_tau & 400 (nat, nat_blk, ntyp, bg_blk, tau, tau_blk, ityp, itau_blk) 401 ELSE 402 CALL set_tau & 403 (nat, nat_blk, at, at_blk, tau, tau_blk, ityp, ityp_blk, itau_blk) 404 ENDIF 405 ! 406 IF (fltau.NE.' ') CALL write_tau (fltau, nat, tau, ityp) 407 ! 408 ! reciprocal lattice vectors 409 ! 410 CALL recips (at(1,1),at(1,2),at(1,3),bg(1,1),bg(1,2),bg(1,3)) 411 ! 412 ! build the WS cell corresponding to the force constant grid 413 ! 414 atws(:,1) = at_blk(:,1)*DBLE(nr1) 415 atws(:,2) = at_blk(:,2)*DBLE(nr2) 416 atws(:,3) = at_blk(:,3)*DBLE(nr3) 417 ! initialize WS r-vectors 418 CALL wsinit(rws,nrwsx,nrws,atws) 419 ! 420 ! end of (super)cell setup 421 ! 422 IF (dos) THEN 423 IF (nk1 < 1 .OR. nk2 < 1 .OR. nk3 < 1) & 424 CALL errore ('matdyn','specify correct q-point grid!',1) 425 nqx = nk1*nk2*nk3 426 ALLOCATE ( q(3,nqx), wq(nqx) ) 427 CALL gen_qpoints (ibrav, at, bg, nat, tau, ityp, nk1, nk2, nk3, & 428 nqx, nq, q, nosym, wq) 429 ELSE 430 ! 431 ! read q-point list 432 ! 433 IF (ionode) READ (5,*) nq 434 CALL mp_bcast(nq, ionode_id, world_comm) 435 ALLOCATE ( q(3,nq) ) 436 IF (.NOT.q_in_band_form) THEN 437 ALLOCATE(wq(nq)) 438 DO n = 1,nq 439 IF (ionode) READ (5,*) (q(i,n),i=1,3) 440 END DO 441 CALL mp_bcast(q, ionode_id, world_comm) 442 ! 443 IF (q_in_cryst_coord) CALL cryst_to_cart(nq,q,bg,+1) 444 ELSE 445 ALLOCATE( nqb(nq) ) 446 ALLOCATE( xqaux(3,nq) ) 447 ALLOCATE( letter(nq) ) 448 ALLOCATE( label_list(nq) ) 449 npk_label=0 450 DO n = 1, nq 451 CALL read_line( input_line, end_of_file = tend, error = terr ) 452 IF (tend) CALL errore('matdyn','Missing lines',1) 453 IF (terr) CALL errore('matdyn','Error reading q points',1) 454 DO j=1,256 ! loop over all characters of input_line 455 IF ( (ICHAR(input_line(j:j)) < 58 .AND. & ! a digit 456 ICHAR(input_line(j:j)) > 47) & 457 .OR.ICHAR(input_line(j:j)) == 43 .OR. & ! the + sign 458 ICHAR(input_line(j:j)) == 45 .OR. & ! the - sign 459 ICHAR(input_line(j:j)) == 46 ) THEN ! a dot . 460! 461! This is a digit, therefore this line contains the coordinates of the 462! k point. We read it and exit from the loop on characters 463! 464 READ(input_line,*) xqaux(1,n), xqaux(2,n), xqaux(3,n), & 465 nqb(n) 466 EXIT 467 ELSEIF ((ICHAR(input_line(j:j)) < 123 .AND. & 468 ICHAR(input_line(j:j)) > 64)) THEN 469! 470! This is a letter, not a space character. We read the next three 471! characters and save them in the letter array, save also which k point 472! it is 473! 474 npk_label=npk_label+1 475 READ(input_line(j:),'(a3)') letter(npk_label) 476 label_list(npk_label)=n 477! 478! now we remove the letters from input_line and read the number of points 479! of the line. The next two line should account for the case in which 480! there is only one space between the letter and the number of points. 481! 482 nch=3 483 IF ( ICHAR(input_line(j+1:j+1))==32 .OR. & 484 ICHAR(input_line(j+2:j+2))==32 ) nch=2 485 buffer=input_line(j+nch:) 486 READ(buffer,*,err=20,iostat=ios) nqb(n) 48720 IF (ios /=0) CALL errore('matdyn',& 488 'problem reading number of points',1) 489 EXIT 490 ENDIF 491 ENDDO 492 ENDDO 493 IF (q_in_cryst_coord) k_points='crystal' 494 IF ( npk_label > 0 ) & 495 CALL transform_label_coord(ibrav, celldm, xqaux, letter, & 496 label_list, npk_label, nq, k_points, point_label_type ) 497 498 DEALLOCATE(letter) 499 DEALLOCATE(label_list) 500 501 CALL mp_bcast(xqaux, ionode_id, world_comm) 502 CALL mp_bcast(nqb, ionode_id, world_comm) 503 IF (q_in_cryst_coord) CALL cryst_to_cart(nq,xqaux,bg,+1) 504 nqtot=SUM(nqb(1:nq-1))+1 505 DO i=1,nq-1 506 IF (nqb(i)==0) nqtot=nqtot+1 507 ENDDO 508 DEALLOCATE(q) 509 ALLOCATE(q(3,nqtot)) 510 ALLOCATE(wq(nqtot)) 511 CALL generate_k_along_lines(nq, xqaux, nqb, q, wq, nqtot) 512 nq=nqtot 513 DEALLOCATE(xqaux) 514 DEALLOCATE(nqb) 515 END IF 516 ! 517 END IF 518 ! 519 IF (asr /= 'no') THEN 520 CALL set_asr (asr, nr1, nr2, nr3, frc, zeu, & 521 nat_blk, ibrav, tau_blk) 522 END IF 523 ! 524 IF (flvec.EQ.' ') THEN 525 iout=0 526 ELSE 527 iout=4 528 IF (ionode) OPEN (unit=iout,file=flvec,status='unknown',form='formatted') 529 END IF 530 531 IF (fldyn.EQ.' ') THEN 532 iout_dyn=0 533 ELSE 534 iout_dyn=44 535 OPEN (unit=iout_dyn,file=fldyn,status='unknown',form='formatted') 536 END IF 537 538 IF (fleig.EQ.' ') THEN 539 iout_eig=0 540 ELSE 541 iout_eig=313 542 IF (ionode) OPEN (unit=iout_eig,file=fleig,status='unknown',form='formatted') 543 END IF 544 545 ALLOCATE ( dyn(3,3,nat,nat), dyn_blk(3,3,nat_blk,nat_blk) ) 546 ALLOCATE ( z(3*nat,3*nat), w2(3*nat,nq), f_of_q(3,3,nat,nat), & 547 dynq(3*nat,nq,nat), DOSofE(nat) ) 548 ALLOCATE ( tmp_w2(3*nat), abs_similarity(3*nat,3*nat), mask(3*nat) ) 549 550 if(la2F.and.ionode) open(unit=300,file='dyna2F',status='unknown') 551 IF (xmlifc) CALL set_sym(nat, tau, ityp, nspin_mag, m_loc ) 552 553 ALLOCATE(num_rap_mode(3*nat,nq)) 554 ALLOCATE(high_sym(nq)) 555 num_rap_mode=-1 556 high_sym=.TRUE. 557 558 DO n=1, nq 559 dyn(:,:,:,:) = (0.d0, 0.d0) 560 561 lo_to_split=.FALSE. 562 f_of_q(:,:,:,:) = (0.d0,0.d0) 563 564 IF(na_ifc) THEN 565 566 qq=sqrt(q(1,n)**2+q(2,n)**2+q(3,n)**2) 567 if(abs(qq) < 1d-8) qq=1.0 568 qhat(1)=q(1,n)/qq 569 qhat(2)=q(2,n)/qq 570 qhat(3)=q(3,n)/qq 571 572 CALL nonanal_ifc (nat,nat_blk,itau_blk,epsil,qhat,zeu,omega,dyn, & 573 nr1, nr2, nr3,f_of_q) 574 END IF 575 576 CALL setupmat (q(1,n), dyn, nat, at, bg, tau, itau_blk, nsc, alat, & 577 dyn_blk, nat_blk, at_blk, bg_blk, tau_blk, omega_blk, & 578 loto_2d, & 579 epsil, zeu, frc, nr1,nr2,nr3, has_zstar, rws, nrws, na_ifc,f_of_q,fd) 580 IF (.not.loto_2d) THEN 581 qhat(1) = q(1,n)*at(1,1)+q(2,n)*at(2,1)+q(3,n)*at(3,1) 582 qhat(2) = q(1,n)*at(1,2)+q(2,n)*at(2,2)+q(3,n)*at(3,2) 583 qhat(3) = q(1,n)*at(1,3)+q(2,n)*at(2,3)+q(3,n)*at(3,3) 584 IF ( ABS( qhat(1) - NINT (qhat(1) ) ) <= eps .AND. & 585 ABS( qhat(2) - NINT (qhat(2) ) ) <= eps .AND. & 586 ABS( qhat(3) - NINT (qhat(3) ) ) <= eps ) THEN 587 ! 588 ! q = 0 : we need the direction q => 0 for the non-analytic part 589 ! 590 IF ( n == 1 ) THEN 591 ! if q is the first point in the list 592 IF ( nq > 1 ) THEN 593 ! one more point 594 qhat(:) = q(:,n) - q(:,n+1) 595 ELSE 596 ! no more points 597 qhat(:) = 0.d0 598 END IF 599 ELSE IF ( n > 1 ) THEN 600 ! if q is not the first point in the list 601 IF ( q(1,n-1)==0.d0 .AND. & 602 q(2,n-1)==0.d0 .AND. & 603 q(3,n-1)==0.d0 .AND. n < nq ) THEN 604 ! if the preceding q is also 0 : 605 qhat(:) = q(:,n) - q(:,n+1) 606 ELSE 607 ! if the preceding q is npt 0 : 608 qhat(:) = q(:,n) - q(:,n-1) 609 END IF 610 END IF 611 qh = SQRT(qhat(1)**2+qhat(2)**2+qhat(3)**2) 612 ! write(*,*) ' qh, has_zstar ',qh, has_zstar 613 IF (qh /= 0.d0) qhat(:) = qhat(:) / qh 614 IF (qh /= 0.d0 .AND. .NOT. has_zstar) THEN 615 CALL infomsg & 616 ('matdyn','Z* not found in file '//TRIM(flfrc)// & 617 ', TO-LO splitting at q=0 will be absent!') 618 ELSEIF (loto_disable) THEN 619 CALL infomsg('matdyn', & 620 'loto_disable is true. Disable LO-TO splitting at q=0.') 621 ELSE 622 lo_to_split=.TRUE. 623 ENDIF 624 ! 625 IF (lo_to_split) CALL nonanal (nat, nat_blk, itau_blk, epsil, qhat, zeu, omega, dyn) 626 ! 627 END IF 628 ! 629 END IF 630 631 if(iout_dyn.ne.0) THEN 632 call write_dyn_on_file(q(1,n),dyn,nat, iout_dyn) 633 if(sum(abs(q(:,n)))==0._dp) call write_epsilon_and_zeu (zeu, epsil, nat, iout_dyn) 634 endif 635 636 637 CALL dyndiag(nat,ntyp,amass,ityp,dyn,w2(1,n),z) 638 ! Atom projection of dynamical matrix 639 DO i = 1, 3*nat 640 DO na = 1, nat 641 dynq(i, n, na) = DOT_PRODUCT(z(3*(na-1)+1:3*na, i), & 642 & z(3*(na-1)+1:3*na, i) ) & 643 & * amu_ry * amass(ityp(na)) 644 END DO 645 END DO 646 IF (ionode.and.iout_eig.ne.0) & 647 & CALL write_eigenvectors(nat,ntyp,amass,ityp,q(1,n),w2(1,n),z,iout_eig) 648 ! 649 ! Cannot use the small group of \Gamma to analize the symmetry 650 ! of the mode if there is an electric field. 651 ! 652 IF (xmlifc.AND..NOT.lo_to_split) THEN 653 WRITE(stdout,'(10x,"xq=",3F8.4)') q(:,n) 654 CALL find_representations_mode_q(nat,ntyp,q(:,n), & 655 w2(:,n),z,tau,ityp,amass, num_rap_mode(:,n), nspin_mag) 656 IF (code_group==code_group_old.OR.high_sym(n-1)) high_sym(n)=.FALSE. 657 code_group_old=code_group 658 ENDIF 659 660 IF (eigen_similarity) THEN 661 ! ... order phonon dispersions using similarity of eigenvalues 662 ! ... Courtesy of Takeshi Nishimatsu, IMR, Tohoku University 663 IF (.NOT.ALLOCATED(tmp_z)) THEN 664 ALLOCATE(tmp_z(3*nat,3*nat)) 665 ELSE 666 abs_similarity = ABS ( MATMUL ( CONJG( TRANSPOSE(z)), tmp_z ) ) 667 mask(:) = .true. 668 DO na=1,3*nat 669 location = maxloc( abs_similarity(:,na), mask(:) ) 670 mask(location(1)) = .false. 671 tmp_w2(na) = w2(location(1),n) 672 tmp_z(:,na) = z(:,location(1)) 673 END DO 674 w2(:,n) = tmp_w2(:) 675 z(:,:) = tmp_z(:,:) 676 END IF 677 tmp_z(:,:) = z(:,:) 678 ENDIF 679 ! 680 if(la2F.and.ionode) then 681 write(300,*) n 682 do na=1,3*nat 683 write(300,*) (z(na,nb),nb=1,3*nat) 684 end do ! na 685 endif 686 ! 687 688 IF (ionode.and.iout.ne.0) CALL writemodes(nat,q(1,n),w2(1,n),z,iout) 689 690 ! 691 END DO !nq 692 DEALLOCATE (tmp_w2, abs_similarity, mask) 693 IF (eigen_similarity) DEALLOCATE(tmp_z) 694 if(la2F.and.ionode) close(300) 695 ! 696 IF(iout .NE. 0.and.ionode) CLOSE(unit=iout) 697 IF(iout_dyn .NE. 0) CLOSE(unit=iout_dyn) 698 IF(iout_eig .NE. 0) CLOSE(unit=iout_eig) 699 ! 700 ALLOCATE (freq(3*nat, nq)) 701 DO n=1,nq 702 ! freq(i,n) = frequencies in cm^(-1), with negative sign if omega^2 is negative 703 DO i=1,3*nat 704 freq(i,n)= SQRT(ABS(w2(i,n))) * RY_TO_CMM1 705 IF (w2(i,n) < 0.0d0) freq(i,n) = -freq(i,n) 706 END DO 707 END DO 708 ! 709 IF(flfrq.NE.' '.and.ionode) THEN 710 OPEN (unit=2,file=flfrq ,status='unknown',form='formatted') 711 WRITE(2, '(" &plot nbnd=",i4,", nks=",i4," /")') 3*nat, nq 712 DO n=1, nq 713 WRITE(2, '(10x,3f10.6)') q(1,n), q(2,n), q(3,n) 714 WRITE(2,'(6f10.4)') (freq(i,n), i=1,3*nat) 715 END DO 716 CLOSE(unit=2) 717 718 OPEN (unit=2,file=trim(flfrq)//'.gp' ,status='unknown',form='formatted') 719 pathL = 0._dp 720 WRITE(2, '(f10.6,3x,999f10.4)') pathL, (freq(i,1), i=1,3*nat) 721 DO n=2, nq 722 pathL=pathL+(SQRT(SUM( (q(:,n)-q(:,n-1))**2 ))) 723 WRITE(2, '(f10.6,3x,999f10.4)') pathL, (freq(i,n), i=1,3*nat) 724 END DO 725 CLOSE(unit=2) 726 727 END IF 728 ! 729 ! If the force constants are in the xml format we write also 730 ! the file with the representations of each mode 731 ! 732 IF (flfrq.NE.' '.AND.xmlifc.AND.ionode) THEN 733 filename=TRIM(flfrq)//'.rap' 734 OPEN (unit=2,file=filename ,status='unknown',form='formatted') 735 WRITE(2, '(" &plot_rap nbnd_rap=",i4,", nks_rap=",i4," /")') 3*nat, nq 736 DO n=1, nq 737 WRITE(2,'(10x,3f10.6,l6)') q(1,n), q(2,n), q(3,n), high_sym(n) 738 WRITE(2,'(6i10)') (num_rap_mode(i,n), i=1,3*nat) 739 END DO 740 CLOSE(unit=2) 741 END IF 742 ! 743 IF (dos) THEN 744 Emin = 0.0d0 745 Emax = 0.0d0 746 DO n=1,nq 747 DO i=1, 3*nat 748 Emin = MIN (Emin, freq(i,n)) 749 Emax = MAX (Emax, freq(i,n)) 750 END DO 751 END DO 752 ! 753 if (ndos > 1) then 754 DeltaE = (Emax - Emin)/(ndos-1) 755 else 756 ndos = NINT ( (Emax - Emin) / DeltaE + 1.51d0 ) 757 end if 758 IF (ionode) OPEN (unit=2,file=fldos,status='unknown',form='formatted') 759 IF (ionode) WRITE (2, *) "# Frequency[cm^-1] DOS PDOS" 760 DO na = 1, nat 761 dynq(1:3*nat,1:nq,na) = dynq(1:3*nat,1:nq,na) * freq(1:3*nat,1:nq) 762 END DO 763 DO n= 1, ndos 764 E = Emin + (n - 1) * DeltaE 765 DO na = 1, nat 766 DOSofE(na) = 0d0 767 DO i = 1, 3*nat 768 DOSofE(na) = DOSofE(na) & 769 & + dos_gam(3*nat, nq, i, dynq(1:3*nat,1:nq,na), freq, E) 770 END DO 771 END DO 772 ! 773 IF (ionode) WRITE (2, '(2ES18.10,1000ES12.4)') E, SUM(DOSofE(1:nat)), DOSofE(1:nat) 774 END DO 775 IF (ionode) CLOSE(unit=2) 776 END IF !dos 777 DEALLOCATE (z, w2, dyn, dyn_blk) 778 ! 779 ! for a2F 780 ! 781 IF(la2F) THEN 782 ! 783 IF (.NOT. dos) THEN 784 DO isig=1,el_ph_nsigma 785 OPEN (unit=200+isig,file='elph.gamma.'//& 786 TRIM(int_to_char(isig)), status='unknown',form='formatted') 787 WRITE(200+isig, '(" &plot nbnd=",i4,", nks=",i4," /")') 3*nat, nq 788 END DO 789 END IF 790 ! 791 ! convert frequencies to Ry 792 ! 793 freq(:,:)= freq(:,:) / RY_TO_CMM1 794 Emin = Emin / RY_TO_CMM1 795 DeltaE=DeltaE/ RY_TO_CMM1 796 ! 797 call a2Fdos (nat, nq, nr1, nr2, nr3, ibrav, at, bg, tau, alat, & 798 nsc, nat_blk, at_blk, bg_blk, itau_blk, omega_blk, & 799 rws, nrws, dos, Emin, DeltaE, ndos, & 800 asr, q, freq,fd, wq) 801 ! 802 IF (.NOT.dos) THEN 803 DO isig=1,el_ph_nsigma 804 CLOSE(UNIT=200+isig) 805 ENDDO 806 ENDIF 807 END IF 808 DEALLOCATE ( freq) 809 DEALLOCATE(num_rap_mode) 810 DEALLOCATE(high_sym) 811 ! 812 813 CALL environment_end('MATDYN') 814 ! 815 CALL mp_global_end() 816 ! 817 STOP 818 ! 819END PROGRAM matdyn 820! 821!----------------------------------------------------------------------- 822SUBROUTINE readfc ( flfrc, nr1, nr2, nr3, epsil, nat, & 823 ibrav, alat, at, ntyp, amass, omega, has_zstar ) 824 !----------------------------------------------------------------------- 825 ! 826 USE kinds, ONLY : DP 827 USE ifconstants,ONLY : tau => tau_blk, ityp => ityp_blk, frc, zeu 828 USE cell_base, ONLY : celldm 829 USE io_global, ONLY : ionode, ionode_id, stdout 830 USE mp, ONLY : mp_bcast 831 USE mp_world, ONLY : world_comm 832 USE constants, ONLY : amu_ry 833 ! 834 IMPLICIT NONE 835 ! I/O variable 836 CHARACTER(LEN=256) :: flfrc 837 INTEGER :: ibrav, nr1,nr2,nr3,nat, ntyp 838 REAL(DP) :: alat, at(3,3), epsil(3,3) 839 LOGICAL :: has_zstar 840 ! local variables 841 INTEGER :: i, j, na, nb, m1,m2,m3 842 INTEGER :: ibid, jbid, nabid, nbbid, m1bid,m2bid,m3bid 843 REAL(DP) :: amass(ntyp), amass_from_file, omega 844 INTEGER :: nt 845 CHARACTER(LEN=3) :: atm 846 ! 847 ! 848 IF (ionode) OPEN (unit=1,file=flfrc,status='old',form='formatted') 849 ! 850 ! read cell data 851 ! 852 IF (ionode)THEN 853 READ(1,*) ntyp,nat,ibrav,(celldm(i),i=1,6) 854 if (ibrav==0) then 855 read(1,*) ((at(i,j),i=1,3),j=1,3) 856 end if 857 ENDIF 858 CALL mp_bcast(ntyp, ionode_id, world_comm) 859 CALL mp_bcast(nat, ionode_id, world_comm) 860 CALL mp_bcast(ibrav, ionode_id, world_comm) 861 CALL mp_bcast(celldm, ionode_id, world_comm) 862 IF (ibrav==0) THEN 863 CALL mp_bcast(at, ionode_id, world_comm) 864 ENDIF 865 ! 866 CALL latgen(ibrav,celldm,at(1,1),at(1,2),at(1,3),omega) 867 alat = celldm(1) 868 at = at / alat ! bring at in units of alat 869 CALL volume(alat,at(1,1),at(1,2),at(1,3),omega) 870 ! 871 ! read atomic types, positions and masses 872 ! 873 DO nt = 1,ntyp 874 IF (ionode) READ(1,*) i,atm,amass_from_file 875 CALL mp_bcast(i,ionode_id, world_comm) 876 CALL mp_bcast(atm,ionode_id, world_comm) 877 CALL mp_bcast(amass_from_file,ionode_id, world_comm) 878 IF (i.NE.nt) CALL errore ('readfc','wrong data read',nt) 879 IF (amass(nt).EQ.0.d0) THEN 880 amass(nt) = amass_from_file/amu_ry 881 ELSE 882 WRITE(stdout,*) 'for atomic type',nt,' mass from file not used' 883 END IF 884 END DO 885 ! 886 ALLOCATE (tau(3,nat), ityp(nat), zeu(3,3,nat)) 887 ! 888 DO na=1,nat 889 IF (ionode) READ(1,*) i,ityp(na),(tau(j,na),j=1,3) 890 CALL mp_bcast(i,ionode_id, world_comm) 891 IF (i.NE.na) CALL errore ('readfc','wrong data read',na) 892 END DO 893 CALL mp_bcast(ityp,ionode_id, world_comm) 894 CALL mp_bcast(tau,ionode_id, world_comm) 895 ! 896 ! read macroscopic variable 897 ! 898 IF (ionode) READ (1,*) has_zstar 899 CALL mp_bcast(has_zstar,ionode_id, world_comm) 900 IF (has_zstar) THEN 901 IF (ionode) READ(1,*) ((epsil(i,j),j=1,3),i=1,3) 902 CALL mp_bcast(epsil,ionode_id, world_comm) 903 IF (ionode) THEN 904 DO na=1,nat 905 READ(1,*) 906 READ(1,*) ((zeu(i,j,na),j=1,3),i=1,3) 907 END DO 908 ENDIF 909 CALL mp_bcast(zeu,ionode_id, world_comm) 910 ELSE 911 zeu (:,:,:) = 0.d0 912 epsil(:,:) = 0.d0 913 END IF 914 ! 915 IF (ionode) READ (1,*) nr1,nr2,nr3 916 CALL mp_bcast(nr1,ionode_id, world_comm) 917 CALL mp_bcast(nr2,ionode_id, world_comm) 918 CALL mp_bcast(nr3,ionode_id, world_comm) 919 ! 920 ! read real-space interatomic force constants 921 ! 922 ALLOCATE ( frc(nr1,nr2,nr3,3,3,nat,nat) ) 923 frc(:,:,:,:,:,:,:) = 0.d0 924 DO i=1,3 925 DO j=1,3 926 DO na=1,nat 927 DO nb=1,nat 928 IF (ionode) READ (1,*) ibid, jbid, nabid, nbbid 929 CALL mp_bcast(ibid,ionode_id, world_comm) 930 CALL mp_bcast(jbid,ionode_id, world_comm) 931 CALL mp_bcast(nabid,ionode_id, world_comm) 932 CALL mp_bcast(nbbid,ionode_id, world_comm) 933 IF(i .NE.ibid .OR. j .NE.jbid .OR. & 934 na.NE.nabid .OR. nb.NE.nbbid) & 935 CALL errore ('readfc','error in reading',1) 936 IF (ionode) READ (1,*) (((m1bid, m2bid, m3bid, & 937 frc(m1,m2,m3,i,j,na,nb), & 938 m1=1,nr1),m2=1,nr2),m3=1,nr3) 939 940 CALL mp_bcast(frc(:,:,:,i,j,na,nb),ionode_id, world_comm) 941 END DO 942 END DO 943 END DO 944 END DO 945 ! 946 IF (ionode) CLOSE(unit=1) 947 ! 948 RETURN 949END SUBROUTINE readfc 950! 951!----------------------------------------------------------------------- 952SUBROUTINE frc_blk(dyn,q,tau,nat,nr1,nr2,nr3,frc,at,bg,rws,nrws,f_of_q,fd) 953 !----------------------------------------------------------------------- 954 ! calculates the dynamical matrix at q from the (short-range part of the) 955 ! force constants 956 ! 957 USE kinds, ONLY : DP 958 USE constants, ONLY : tpi 959 USE io_global, ONLY : stdout 960 ! 961 IMPLICIT NONE 962 INTEGER nr1, nr2, nr3, nat, n1, n2, n3, nr1_, nr2_, nr3_, & 963 ipol, jpol, na, nb, m1, m2, m3, nint, i,j, nrws, nax 964 COMPLEX(DP) dyn(3,3,nat,nat), f_of_q(3,3,nat,nat) 965 REAL(DP) frc(nr1,nr2,nr3,3,3,nat,nat), tau(3,nat), q(3), arg, & 966 at(3,3), bg(3,3), r(3), weight, r_ws(3), & 967 total_weight, rws(0:3,nrws), alat 968 REAL(DP), EXTERNAL :: wsweight 969 REAL(DP),SAVE,ALLOCATABLE :: wscache(:,:,:,:,:) 970 REAL(DP), ALLOCATABLE :: ttt(:,:,:,:,:), tttx(:,:) 971 LOGICAL,SAVE :: first=.true. 972 LOGICAL :: fd 973 ! 974 nr1_=2*nr1 975 nr2_=2*nr2 976 nr3_=2*nr3 977 FIRST_TIME : IF (first) THEN 978 first=.false. 979 ALLOCATE( wscache(-nr3_:nr3_, -nr2_:nr2_, -nr1_:nr1_, nat,nat) ) 980 DO na=1, nat 981 DO nb=1, nat 982 total_weight=0.0d0 983 ! 984 ! SUM OVER R VECTORS IN THE SUPERCELL - VERY VERY VERY SAFE RANGE! 985 ! 986 DO n1=-nr1_,nr1_ 987 DO n2=-nr2_,nr2_ 988 DO n3=-nr3_,nr3_ 989 DO i=1, 3 990 r(i) = n1*at(i,1)+n2*at(i,2)+n3*at(i,3) 991 r_ws(i) = r(i) + tau(i,na)-tau(i,nb) 992 if (fd) r_ws(i) = r(i) + tau(i,nb)-tau(i,na) 993 END DO 994 wscache(n3,n2,n1,nb,na) = wsweight(r_ws,rws,nrws) 995 total_weight=total_weight + wscache(n3,n2,n1,nb,na) 996 ENDDO 997 ENDDO 998 ENDDO 999 IF (ABS(total_weight-nr1*nr2*nr3).GT.1.0d-8) THEN 1000 WRITE(stdout,*) na,nb,total_weight 1001 CALL errore ('frc_blk','wrong total_weight',1) 1002 END IF 1003 ENDDO 1004 ENDDO 1005 ENDIF FIRST_TIME 1006 ! 1007 ALLOCATE(ttt(3,nat,nr1,nr2,nr3)) 1008 ALLOCATE(tttx(3,nat*nr1*nr2*nr3)) 1009 ttt(:,:,:,:,:)=0.d0 1010 1011 DO na=1, nat 1012 DO nb=1, nat 1013 DO n1=-nr1_,nr1_ 1014 DO n2=-nr2_,nr2_ 1015 DO n3=-nr3_,nr3_ 1016 ! 1017 ! SUM OVER R VECTORS IN THE SUPERCELL - VERY VERY SAFE RANGE! 1018 ! 1019 DO i=1, 3 1020 r(i) = n1*at(i,1)+n2*at(i,2)+n3*at(i,3) 1021 END DO 1022 1023 weight = wscache(n3,n2,n1,nb,na) 1024 IF (weight .GT. 0.0d0) THEN 1025 ! 1026 ! FIND THE VECTOR CORRESPONDING TO R IN THE ORIGINAL CELL 1027 ! 1028 m1 = MOD(n1+1,nr1) 1029 IF(m1.LE.0) m1=m1+nr1 1030 m2 = MOD(n2+1,nr2) 1031 IF(m2.LE.0) m2=m2+nr2 1032 m3 = MOD(n3+1,nr3) 1033 IF(m3.LE.0) m3=m3+nr3 1034 ! write(*,'(6i4)') n1,n2,n3,m1,m2,m3 1035 ! 1036 ! FOURIER TRANSFORM 1037 ! 1038 do i=1,3 1039 ttt(i,na,m1,m2,m3)=tau(i,na)+m1*at(i,1)+m2*at(i,2)+m3*at(i,3) 1040 ttt(i,nb,m1,m2,m3)=tau(i,nb)+m1*at(i,1)+m2*at(i,2)+m3*at(i,3) 1041 end do 1042 1043 arg = tpi*(q(1)*r(1) + q(2)*r(2) + q(3)*r(3)) 1044 DO ipol=1, 3 1045 DO jpol=1, 3 1046 dyn(ipol,jpol,na,nb) = dyn(ipol,jpol,na,nb) + & 1047 (frc(m1,m2,m3,ipol,jpol,na,nb)+f_of_q(ipol,jpol,na,nb)) & 1048 *CMPLX(COS(arg),-SIN(arg),kind=DP)*weight 1049 END DO 1050 END DO 1051 1052 END IF 1053 END DO 1054 END DO 1055 END DO 1056 END DO 1057 END DO 1058 ! 1059 RETURN 1060END SUBROUTINE frc_blk 1061! 1062!----------------------------------------------------------------------- 1063SUBROUTINE setupmat (q,dyn,nat,at,bg,tau,itau_blk,nsc,alat, & 1064 & dyn_blk,nat_blk,at_blk,bg_blk,tau_blk,omega_blk, & 1065 & loto_2d, & 1066 & epsil,zeu,frc,nr1,nr2,nr3,has_zstar,rws,nrws,na_ifc,f_of_q,fd) 1067 !----------------------------------------------------------------------- 1068 ! compute the dynamical matrix (the analytic part only) 1069 ! 1070 USE kinds, ONLY : DP 1071 USE constants, ONLY : tpi 1072 USE cell_base, ONLY : celldm 1073 USE rigid, ONLY : rgd_blk 1074 ! 1075 IMPLICIT NONE 1076 ! 1077 ! I/O variables 1078 ! 1079 INTEGER:: nr1, nr2, nr3, nat, nat_blk, nsc, nrws, itau_blk(nat) 1080 REAL(DP) :: q(3), tau(3,nat), at(3,3), bg(3,3), alat, & 1081 epsil(3,3), zeu(3,3,nat_blk), rws(0:3,nrws), & 1082 frc(nr1,nr2,nr3,3,3,nat_blk,nat_blk) 1083 REAL(DP) :: tau_blk(3,nat_blk), at_blk(3,3), bg_blk(3,3), omega_blk 1084 COMPLEX(DP) dyn_blk(3,3,nat_blk,nat_blk), f_of_q(3,3,nat,nat) 1085 COMPLEX(DP) :: dyn(3,3,nat,nat) 1086 LOGICAL :: has_zstar, na_ifc, fd, loto_2d 1087 ! 1088 ! local variables 1089 ! 1090 REAL(DP) :: arg 1091 COMPLEX(DP) :: cfac(nat) 1092 INTEGER :: i,j,k, na,nb, na_blk, nb_blk, iq 1093 REAL(DP) :: qp(3), qbid(3,nsc) ! automatic array 1094 ! 1095 ! 1096 CALL q_gen(nsc,qbid,at_blk,bg_blk,at,bg) 1097 ! 1098 DO iq=1,nsc 1099 ! 1100 DO k=1,3 1101 qp(k)= q(k) + qbid(k,iq) 1102 END DO 1103 ! 1104 dyn_blk(:,:,:,:) = (0.d0,0.d0) 1105 CALL frc_blk (dyn_blk,qp,tau_blk,nat_blk, & 1106 & nr1,nr2,nr3,frc,at_blk,bg_blk,rws,nrws,f_of_q,fd) 1107 IF (has_zstar .and. .not.na_ifc) & 1108 CALL rgd_blk(nr1,nr2,nr3,nat_blk,dyn_blk,qp,tau_blk, & 1109 epsil,zeu,bg_blk,omega_blk,celldm(1), loto_2d,+1.d0) 1110 ! LOTO 2D added celldm(1)=alat to passed arguments 1111 ! 1112 DO na=1,nat 1113 na_blk = itau_blk(na) 1114 DO nb=1,nat 1115 nb_blk = itau_blk(nb) 1116 ! 1117 arg=tpi* ( qp(1) * ( (tau(1,na)-tau_blk(1,na_blk)) - & 1118 (tau(1,nb)-tau_blk(1,nb_blk)) ) + & 1119 qp(2) * ( (tau(2,na)-tau_blk(2,na_blk)) - & 1120 (tau(2,nb)-tau_blk(2,nb_blk)) ) + & 1121 qp(3) * ( (tau(3,na)-tau_blk(3,na_blk)) - & 1122 (tau(3,nb)-tau_blk(3,nb_blk)) ) ) 1123 ! 1124 cfac(nb) = CMPLX(COS(arg),SIN(arg),kind=DP)/nsc 1125 ! 1126 END DO ! nb 1127 ! 1128 DO i=1,3 1129 DO j=1,3 1130 ! 1131 DO nb=1,nat 1132 nb_blk = itau_blk(nb) 1133 dyn(i,j,na,nb) = dyn(i,j,na,nb) + cfac(nb) * & 1134 dyn_blk(i,j,na_blk,nb_blk) 1135 END DO ! nb 1136 ! 1137 END DO ! j 1138 END DO ! i 1139 END DO ! na 1140 ! 1141 END DO ! iq 1142 ! 1143 RETURN 1144END SUBROUTINE setupmat 1145! 1146! 1147!---------------------------------------------------------------------- 1148SUBROUTINE set_asr (asr, nr1, nr2, nr3, frc, zeu, nat, ibrav, tau) 1149 !----------------------------------------------------------------------- 1150 ! 1151 USE kinds, ONLY : DP 1152 USE io_global, ONLY : stdout 1153 ! 1154 IMPLICIT NONE 1155 CHARACTER (LEN=10), intent(in) :: asr 1156 INTEGER, intent(in) :: nr1, nr2, nr3, nat, ibrav 1157 REAL(DP), intent(in) :: tau(3,nat) 1158 REAL(DP), intent(inout) :: frc(nr1,nr2,nr3,3,3,nat,nat), zeu(3,3,nat) 1159 ! 1160 INTEGER :: axis, n, i, j, na, nb, n1,n2,n3, m,p,k,l,q,r, i1,j1,na1 1161 REAL(DP) :: zeu_new(3,3,nat) 1162 REAL(DP), ALLOCATABLE :: frc_new(:,:,:,:,:,:,:) 1163 type vector 1164 real(DP),pointer :: vec(:,:,:,:,:,:,:) 1165 end type vector 1166 ! 1167 type (vector) u(6*3*nat) 1168 ! These are the "vectors" associated with the sum rules on force-constants 1169 ! 1170 integer :: u_less(6*3*nat),n_less,i_less 1171 ! indices of the vectors u that are not independent to the preceding ones, 1172 ! n_less = number of such vectors, i_less = temporary parameter 1173 ! 1174 integer, allocatable :: ind_v(:,:,:) 1175 real(DP), allocatable :: v(:,:) 1176 ! These are the "vectors" associated with symmetry conditions, coded by 1177 ! indicating the positions (i.e. the seven indices) of the non-zero elements (there 1178 ! should be only 2 of them) and the value of that element. We do so in order 1179 ! to limit the amount of memory used. 1180 ! 1181 real(DP), allocatable :: w(:,:,:,:,:,:,:), x(:,:,:,:,:,:,:) 1182 ! temporary vectors and parameters 1183 real(DP) :: scal,norm2, sum 1184 ! 1185 real(DP) :: zeu_u(6*3,3,3,nat) 1186 ! These are the "vectors" associated with the sum rules on effective charges 1187 ! 1188 integer :: zeu_less(6*3),nzeu_less,izeu_less 1189 ! indices of the vectors zeu_u that are not independent to the preceding ones, 1190 ! nzeu_less = number of such vectors, izeu_less = temporary parameter 1191 ! 1192 real(DP) :: zeu_w(3,3,nat), zeu_x(3,3,nat) 1193 ! temporary vectors 1194 1195 ! Initialization. n is the number of sum rules to be considered (if asr.ne.'simple') 1196 ! and 'axis' is the rotation axis in the case of a 1D system 1197 ! (i.e. the rotation axis is (Ox) if axis='1', (Oy) if axis='2' and (Oz) if axis='3') 1198 ! 1199 if((asr.ne.'simple').and.(asr.ne.'crystal').and.(asr.ne.'one-dim') & 1200 .and.(asr.ne.'zero-dim')) then 1201 call errore('set_asr','invalid Acoustic Sum Rule:' // asr, 1) 1202 endif 1203 ! 1204 if(asr.eq.'simple') then 1205 ! 1206 ! Simple Acoustic Sum Rule on effective charges 1207 ! 1208 do i=1,3 1209 do j=1,3 1210 sum=0.0d0 1211 do na=1,nat 1212 sum = sum + zeu(i,j,na) 1213 end do 1214 do na=1,nat 1215 zeu(i,j,na) = zeu(i,j,na) - sum/nat 1216 end do 1217 end do 1218 end do 1219 ! 1220 ! Simple Acoustic Sum Rule on force constants in real space 1221 ! 1222 do i=1,3 1223 do j=1,3 1224 do na=1,nat 1225 sum=0.0d0 1226 do nb=1,nat 1227 do n1=1,nr1 1228 do n2=1,nr2 1229 do n3=1,nr3 1230 sum=sum+frc(n1,n2,n3,i,j,na,nb) 1231 end do 1232 end do 1233 end do 1234 end do 1235 frc(1,1,1,i,j,na,na) = frc(1,1,1,i,j,na,na) - sum 1236 ! write(6,*) ' na, i, j, sum = ',na,i,j,sum 1237 end do 1238 end do 1239 end do 1240 ! 1241 return 1242 ! 1243 end if 1244 if(asr.eq.'crystal') n=3 1245 if(asr.eq.'one-dim') then 1246 ! the direction of periodicity is the rotation axis 1247 ! It will work only if the crystal axis considered is one of 1248 ! the cartesian axis (typically, ibrav=1, 6 or 8, or 4 along the 1249 ! z-direction) 1250 if (nr1*nr2*nr3.eq.1) axis=3 1251 if ((nr1.ne.1).and.(nr2*nr3.eq.1)) axis=1 1252 if ((nr2.ne.1).and.(nr1*nr3.eq.1)) axis=2 1253 if ((nr3.ne.1).and.(nr1*nr2.eq.1)) axis=3 1254 if (((nr1.ne.1).and.(nr2.ne.1)).or.((nr2.ne.1).and. & 1255 (nr3.ne.1)).or.((nr1.ne.1).and.(nr3.ne.1))) then 1256 call errore('set_asr','too many directions of & 1257 & periodicity in 1D system',axis) 1258 endif 1259 if ((ibrav.ne.1).and.(ibrav.ne.6).and.(ibrav.ne.8).and. & 1260 ((ibrav.ne.4).or.(axis.ne.3)) ) then 1261 write(stdout,*) 'asr: rotational axis may be wrong' 1262 endif 1263 write(stdout,'("asr rotation axis in 1D system= ",I4)') axis 1264 n=4 1265 endif 1266 if(asr.eq.'zero-dim') n=6 1267 ! 1268 ! Acoustic Sum Rule on effective charges 1269 ! 1270 ! generating the vectors of the orthogonal of the subspace to project 1271 ! the effective charges matrix on 1272 ! 1273 zeu_u(:,:,:,:)=0.0d0 1274 do i=1,3 1275 do j=1,3 1276 do na=1,nat 1277 zeu_new(i,j,na)=zeu(i,j,na) 1278 enddo 1279 enddo 1280 enddo 1281 ! 1282 p=0 1283 do i=1,3 1284 do j=1,3 1285 ! These are the 3*3 vectors associated with the 1286 ! translational acoustic sum rules 1287 p=p+1 1288 zeu_u(p,i,j,:)=1.0d0 1289 ! 1290 enddo 1291 enddo 1292 ! 1293 if (n.eq.4) then 1294 do i=1,3 1295 ! These are the 3 vectors associated with the 1296 ! single rotational sum rule (1D system) 1297 p=p+1 1298 do na=1,nat 1299 zeu_u(p,i,MOD(axis,3)+1,na)=-tau(MOD(axis+1,3)+1,na) 1300 zeu_u(p,i,MOD(axis+1,3)+1,na)=tau(MOD(axis,3)+1,na) 1301 enddo 1302 ! 1303 enddo 1304 endif 1305 ! 1306 if (n.eq.6) then 1307 do i=1,3 1308 do j=1,3 1309 ! These are the 3*3 vectors associated with the 1310 ! three rotational sum rules (0D system - typ. molecule) 1311 p=p+1 1312 do na=1,nat 1313 zeu_u(p,i,MOD(j,3)+1,na)=-tau(MOD(j+1,3)+1,na) 1314 zeu_u(p,i,MOD(j+1,3)+1,na)=tau(MOD(j,3)+1,na) 1315 enddo 1316 ! 1317 enddo 1318 enddo 1319 endif 1320 ! 1321 ! Gram-Schmidt orthonormalization of the set of vectors created. 1322 ! 1323 nzeu_less=0 1324 do k=1,p 1325 zeu_w(:,:,:)=zeu_u(k,:,:,:) 1326 zeu_x(:,:,:)=zeu_u(k,:,:,:) 1327 do q=1,k-1 1328 r=1 1329 do izeu_less=1,nzeu_less 1330 if (zeu_less(izeu_less).eq.q) r=0 1331 enddo 1332 if (r.ne.0) then 1333 call sp_zeu(zeu_x,zeu_u(q,:,:,:),nat,scal) 1334 zeu_w(:,:,:) = zeu_w(:,:,:) - scal* zeu_u(q,:,:,:) 1335 endif 1336 enddo 1337 call sp_zeu(zeu_w,zeu_w,nat,norm2) 1338 if (norm2.gt.1.0d-16) then 1339 zeu_u(k,:,:,:) = zeu_w(:,:,:) / DSQRT(norm2) 1340 else 1341 nzeu_less=nzeu_less+1 1342 zeu_less(nzeu_less)=k 1343 endif 1344 enddo 1345 ! 1346 ! Projection of the effective charge "vector" on the orthogonal of the 1347 ! subspace of the vectors verifying the sum rules 1348 ! 1349 zeu_w(:,:,:)=0.0d0 1350 do k=1,p 1351 r=1 1352 do izeu_less=1,nzeu_less 1353 if (zeu_less(izeu_less).eq.k) r=0 1354 enddo 1355 if (r.ne.0) then 1356 zeu_x(:,:,:)=zeu_u(k,:,:,:) 1357 call sp_zeu(zeu_x,zeu_new,nat,scal) 1358 zeu_w(:,:,:) = zeu_w(:,:,:) + scal*zeu_u(k,:,:,:) 1359 endif 1360 enddo 1361 ! 1362 ! Final substraction of the former projection to the initial zeu, to get 1363 ! the new "projected" zeu 1364 ! 1365 zeu_new(:,:,:)=zeu_new(:,:,:) - zeu_w(:,:,:) 1366 call sp_zeu(zeu_w,zeu_w,nat,norm2) 1367 write(stdout,'("Norm of the difference between old and new effective ", & 1368 & "charges: ",F25.20)') SQRT(norm2) 1369 ! 1370 ! Check projection 1371 ! 1372 !write(6,'("Check projection of zeu")') 1373 !do k=1,p 1374 ! zeu_x(:,:,:)=zeu_u(k,:,:,:) 1375 ! call sp_zeu(zeu_x,zeu_new,nat,scal) 1376 ! if (DABS(scal).gt.1d-10) write(6,'("k= ",I8," zeu_new|zeu_u(k)= ",F15.10)') k,scal 1377 !enddo 1378 ! 1379 do i=1,3 1380 do j=1,3 1381 do na=1,nat 1382 zeu(i,j,na)=zeu_new(i,j,na) 1383 enddo 1384 enddo 1385 enddo 1386 ! 1387 ! Acoustic Sum Rule on force constants 1388 ! 1389 ! 1390 ! generating the vectors of the orthogonal of the subspace to project 1391 ! the force-constants matrix on 1392 ! 1393 do k=1,18*nat 1394 allocate(u(k) % vec(nr1,nr2,nr3,3,3,nat,nat)) 1395 u(k) % vec (:,:,:,:,:,:,:)=0.0d0 1396 enddo 1397 ALLOCATE (frc_new(nr1,nr2,nr3,3,3,nat,nat)) 1398 do i=1,3 1399 do j=1,3 1400 do na=1,nat 1401 do nb=1,nat 1402 do n1=1,nr1 1403 do n2=1,nr2 1404 do n3=1,nr3 1405 frc_new(n1,n2,n3,i,j,na,nb)=frc(n1,n2,n3,i,j,na,nb) 1406 enddo 1407 enddo 1408 enddo 1409 enddo 1410 enddo 1411 enddo 1412 enddo 1413 ! 1414 p=0 1415 do i=1,3 1416 do j=1,3 1417 do na=1,nat 1418 ! These are the 3*3*nat vectors associated with the 1419 ! translational acoustic sum rules 1420 p=p+1 1421 u(p) % vec (:,:,:,i,j,na,:)=1.0d0 1422 ! 1423 enddo 1424 enddo 1425 enddo 1426 ! 1427 if (n.eq.4) then 1428 do i=1,3 1429 do na=1,nat 1430 ! These are the 3*nat vectors associated with the 1431 ! single rotational sum rule (1D system) 1432 p=p+1 1433 do nb=1,nat 1434 u(p) % vec (:,:,:,i,MOD(axis,3)+1,na,nb)=-tau(MOD(axis+1,3)+1,nb) 1435 u(p) % vec (:,:,:,i,MOD(axis+1,3)+1,na,nb)=tau(MOD(axis,3)+1,nb) 1436 enddo 1437 ! 1438 enddo 1439 enddo 1440 endif 1441 ! 1442 if (n.eq.6) then 1443 do i=1,3 1444 do j=1,3 1445 do na=1,nat 1446 ! These are the 3*3*nat vectors associated with the 1447 ! three rotational sum rules (0D system - typ. molecule) 1448 p=p+1 1449 do nb=1,nat 1450 u(p) % vec (:,:,:,i,MOD(j,3)+1,na,nb)=-tau(MOD(j+1,3)+1,nb) 1451 u(p) % vec (:,:,:,i,MOD(j+1,3)+1,na,nb)=tau(MOD(j,3)+1,nb) 1452 enddo 1453 ! 1454 enddo 1455 enddo 1456 enddo 1457 endif 1458 ! 1459 allocate (ind_v(9*nat*nat*nr1*nr2*nr3,2,7), v(9*nat*nat*nr1*nr2*nr3,2) ) 1460 m=0 1461 do i=1,3 1462 do j=1,3 1463 do na=1,nat 1464 do nb=1,nat 1465 do n1=1,nr1 1466 do n2=1,nr2 1467 do n3=1,nr3 1468 ! These are the vectors associated with the symmetry constraints 1469 q=1 1470 l=1 1471 do while((l.le.m).and.(q.ne.0)) 1472 if ((ind_v(l,1,1).eq.n1).and.(ind_v(l,1,2).eq.n2).and. & 1473 (ind_v(l,1,3).eq.n3).and.(ind_v(l,1,4).eq.i).and. & 1474 (ind_v(l,1,5).eq.j).and.(ind_v(l,1,6).eq.na).and. & 1475 (ind_v(l,1,7).eq.nb)) q=0 1476 if ((ind_v(l,2,1).eq.n1).and.(ind_v(l,2,2).eq.n2).and. & 1477 (ind_v(l,2,3).eq.n3).and.(ind_v(l,2,4).eq.i).and. & 1478 (ind_v(l,2,5).eq.j).and.(ind_v(l,2,6).eq.na).and. & 1479 (ind_v(l,2,7).eq.nb)) q=0 1480 l=l+1 1481 enddo 1482 if ((n1.eq.MOD(nr1+1-n1,nr1)+1).and.(n2.eq.MOD(nr2+1-n2,nr2)+1) & 1483 .and.(n3.eq.MOD(nr3+1-n3,nr3)+1).and.(i.eq.j).and.(na.eq.nb)) q=0 1484 if (q.ne.0) then 1485 m=m+1 1486 ind_v(m,1,1)=n1 1487 ind_v(m,1,2)=n2 1488 ind_v(m,1,3)=n3 1489 ind_v(m,1,4)=i 1490 ind_v(m,1,5)=j 1491 ind_v(m,1,6)=na 1492 ind_v(m,1,7)=nb 1493 v(m,1)=1.0d0/DSQRT(2.0d0) 1494 ind_v(m,2,1)=MOD(nr1+1-n1,nr1)+1 1495 ind_v(m,2,2)=MOD(nr2+1-n2,nr2)+1 1496 ind_v(m,2,3)=MOD(nr3+1-n3,nr3)+1 1497 ind_v(m,2,4)=j 1498 ind_v(m,2,5)=i 1499 ind_v(m,2,6)=nb 1500 ind_v(m,2,7)=na 1501 v(m,2)=-1.0d0/DSQRT(2.0d0) 1502 endif 1503 enddo 1504 enddo 1505 enddo 1506 enddo 1507 enddo 1508 enddo 1509 enddo 1510 ! 1511 ! Gram-Schmidt orthonormalization of the set of vectors created. 1512 ! Note that the vectors corresponding to symmetry constraints are already 1513 ! orthonormalized by construction. 1514 ! 1515 n_less=0 1516 allocate (w(nr1,nr2,nr3,3,3,nat,nat), x(nr1,nr2,nr3,3,3,nat,nat)) 1517 do k=1,p 1518 w(:,:,:,:,:,:,:)=u(k) % vec (:,:,:,:,:,:,:) 1519 x(:,:,:,:,:,:,:)=u(k) % vec (:,:,:,:,:,:,:) 1520 do l=1,m 1521 ! 1522 call sp2(x,v(l,:),ind_v(l,:,:),nr1,nr2,nr3,nat,scal) 1523 do r=1,2 1524 n1=ind_v(l,r,1) 1525 n2=ind_v(l,r,2) 1526 n3=ind_v(l,r,3) 1527 i=ind_v(l,r,4) 1528 j=ind_v(l,r,5) 1529 na=ind_v(l,r,6) 1530 nb=ind_v(l,r,7) 1531 w(n1,n2,n3,i,j,na,nb)=w(n1,n2,n3,i,j,na,nb)-scal*v(l,r) 1532 enddo 1533 enddo 1534 if (k.le.(9*nat)) then 1535 na1=MOD(k,nat) 1536 if (na1.eq.0) na1=nat 1537 j1=MOD((k-na1)/nat,3)+1 1538 i1=MOD((((k-na1)/nat)-j1+1)/3,3)+1 1539 else 1540 q=k-9*nat 1541 if (n.eq.4) then 1542 na1=MOD(q,nat) 1543 if (na1.eq.0) na1=nat 1544 i1=MOD((q-na1)/nat,3)+1 1545 else 1546 na1=MOD(q,nat) 1547 if (na1.eq.0) na1=nat 1548 j1=MOD((q-na1)/nat,3)+1 1549 i1=MOD((((q-na1)/nat)-j1+1)/3,3)+1 1550 endif 1551 endif 1552 do q=1,k-1 1553 r=1 1554 do i_less=1,n_less 1555 if (u_less(i_less).eq.q) r=0 1556 enddo 1557 if (r.ne.0) then 1558 call sp3(x,u(q) % vec (:,:,:,:,:,:,:), i1,na1,nr1,nr2,nr3,nat,scal) 1559 w(:,:,:,:,:,:,:) = w(:,:,:,:,:,:,:) - scal* u(q) % vec (:,:,:,:,:,:,:) 1560 endif 1561 enddo 1562 call sp1(w,w,nr1,nr2,nr3,nat,norm2) 1563 if (norm2.gt.1.0d-16) then 1564 u(k) % vec (:,:,:,:,:,:,:) = w(:,:,:,:,:,:,:) / DSQRT(norm2) 1565 else 1566 n_less=n_less+1 1567 u_less(n_less)=k 1568 endif 1569 enddo 1570 ! 1571 ! Projection of the force-constants "vector" on the orthogonal of the 1572 ! subspace of the vectors verifying the sum rules and symmetry contraints 1573 ! 1574 w(:,:,:,:,:,:,:)=0.0d0 1575 do l=1,m 1576 call sp2(frc_new,v(l,:),ind_v(l,:,:),nr1,nr2,nr3,nat,scal) 1577 do r=1,2 1578 n1=ind_v(l,r,1) 1579 n2=ind_v(l,r,2) 1580 n3=ind_v(l,r,3) 1581 i=ind_v(l,r,4) 1582 j=ind_v(l,r,5) 1583 na=ind_v(l,r,6) 1584 nb=ind_v(l,r,7) 1585 w(n1,n2,n3,i,j,na,nb)=w(n1,n2,n3,i,j,na,nb)+scal*v(l,r) 1586 enddo 1587 enddo 1588 do k=1,p 1589 r=1 1590 do i_less=1,n_less 1591 if (u_less(i_less).eq.k) r=0 1592 enddo 1593 if (r.ne.0) then 1594 x(:,:,:,:,:,:,:)=u(k) % vec (:,:,:,:,:,:,:) 1595 call sp1(x,frc_new,nr1,nr2,nr3,nat,scal) 1596 w(:,:,:,:,:,:,:) = w(:,:,:,:,:,:,:) + scal*u(k)%vec(:,:,:,:,:,:,:) 1597 endif 1598 deallocate(u(k) % vec) 1599 enddo 1600 ! 1601 ! Final substraction of the former projection to the initial frc, to get 1602 ! the new "projected" frc 1603 ! 1604 frc_new(:,:,:,:,:,:,:)=frc_new(:,:,:,:,:,:,:) - w(:,:,:,:,:,:,:) 1605 call sp1(w,w,nr1,nr2,nr3,nat,norm2) 1606 write(stdout,'("Norm of the difference between old and new force-constants:",& 1607 & F25.20)') SQRT(norm2) 1608 ! 1609 ! Check projection 1610 ! 1611 !write(6,'("Check projection IFC")') 1612 !do l=1,m 1613 ! call sp2(frc_new,v(l,:),ind_v(l,:,:),nr1,nr2,nr3,nat,scal) 1614 ! if (DABS(scal).gt.1d-10) write(6,'("l= ",I8," frc_new|v(l)= ",F15.10)') l,scal 1615 !enddo 1616 !do k=1,p 1617 ! x(:,:,:,:,:,:,:)=u(k) % vec (:,:,:,:,:,:,:) 1618 ! call sp1(x,frc_new,nr1,nr2,nr3,nat,scal) 1619 ! if (DABS(scal).gt.1d-10) write(6,'("k= ",I8," frc_new|u(k)= ",F15.10)') k,scal 1620 ! deallocate(u(k) % vec) 1621 !enddo 1622 ! 1623 do i=1,3 1624 do j=1,3 1625 do na=1,nat 1626 do nb=1,nat 1627 do n1=1,nr1 1628 do n2=1,nr2 1629 do n3=1,nr3 1630 frc(n1,n2,n3,i,j,na,nb)=frc_new(n1,n2,n3,i,j,na,nb) 1631 enddo 1632 enddo 1633 enddo 1634 enddo 1635 enddo 1636 enddo 1637 enddo 1638 deallocate (x, w) 1639 deallocate (v, ind_v) 1640 deallocate (frc_new) 1641 ! 1642 return 1643end subroutine set_asr 1644! 1645!---------------------------------------------------------------------- 1646subroutine sp_zeu(zeu_u,zeu_v,nat,scal) 1647 !----------------------------------------------------------------------- 1648 ! 1649 ! does the scalar product of two effective charges matrices zeu_u and zeu_v 1650 ! (considered as vectors in the R^(3*3*nat) space, and coded in the usual way) 1651 ! 1652 USE kinds, ONLY: DP 1653 implicit none 1654 integer i,j,na,nat 1655 real(DP) zeu_u(3,3,nat) 1656 real(DP) zeu_v(3,3,nat) 1657 real(DP) scal 1658 ! 1659 ! 1660 scal=0.0d0 1661 do i=1,3 1662 do j=1,3 1663 do na=1,nat 1664 scal=scal+zeu_u(i,j,na)*zeu_v(i,j,na) 1665 enddo 1666 enddo 1667 enddo 1668 ! 1669 return 1670 ! 1671end subroutine sp_zeu 1672! 1673! 1674!---------------------------------------------------------------------- 1675subroutine sp1(u,v,nr1,nr2,nr3,nat,scal) 1676 !----------------------------------------------------------------------- 1677 ! 1678 ! does the scalar product of two force-constants matrices u and v (considered as 1679 ! vectors in the R^(3*3*nat*nat*nr1*nr2*nr3) space, and coded in the usual way) 1680 ! 1681 USE kinds, ONLY: DP 1682 implicit none 1683 integer nr1,nr2,nr3,i,j,na,nb,n1,n2,n3,nat 1684 real(DP) u(nr1,nr2,nr3,3,3,nat,nat) 1685 real(DP) v(nr1,nr2,nr3,3,3,nat,nat) 1686 real(DP) scal 1687 ! 1688 ! 1689 scal=0.0d0 1690 do i=1,3 1691 do j=1,3 1692 do na=1,nat 1693 do nb=1,nat 1694 do n1=1,nr1 1695 do n2=1,nr2 1696 do n3=1,nr3 1697 scal=scal+u(n1,n2,n3,i,j,na,nb)*v(n1,n2,n3,i,j,na,nb) 1698 enddo 1699 enddo 1700 enddo 1701 enddo 1702 enddo 1703 enddo 1704 enddo 1705 ! 1706 return 1707 ! 1708end subroutine sp1 1709! 1710!---------------------------------------------------------------------- 1711subroutine sp2(u,v,ind_v,nr1,nr2,nr3,nat,scal) 1712 !----------------------------------------------------------------------- 1713 ! 1714 ! does the scalar product of two force-constants matrices u and v (considered as 1715 ! vectors in the R^(3*3*nat*nat*nr1*nr2*nr3) space). u is coded in the usual way 1716 ! but v is coded as explained when defining the vectors corresponding to the 1717 ! symmetry constraints 1718 ! 1719 USE kinds, ONLY: DP 1720 implicit none 1721 integer nr1,nr2,nr3,i,nat 1722 real(DP) u(nr1,nr2,nr3,3,3,nat,nat) 1723 integer ind_v(2,7) 1724 real(DP) v(2) 1725 real(DP) scal 1726 ! 1727 ! 1728 scal=0.0d0 1729 do i=1,2 1730 scal=scal+u(ind_v(i,1),ind_v(i,2),ind_v(i,3),ind_v(i,4),ind_v(i,5),ind_v(i,6), & 1731 ind_v(i,7))*v(i) 1732 enddo 1733 ! 1734 return 1735 ! 1736end subroutine sp2 1737! 1738!---------------------------------------------------------------------- 1739subroutine sp3(u,v,i,na,nr1,nr2,nr3,nat,scal) 1740 !----------------------------------------------------------------------- 1741 ! 1742 ! like sp1, but in the particular case when u is one of the u(k)%vec 1743 ! defined in set_asr (before orthonormalization). In this case most of the 1744 ! terms are zero (the ones that are not are characterized by i and na), so 1745 ! that a lot of computer time can be saved (during Gram-Schmidt). 1746 ! 1747 USE kinds, ONLY: DP 1748 implicit none 1749 integer nr1,nr2,nr3,i,j,na,nb,n1,n2,n3,nat 1750 real(DP) u(nr1,nr2,nr3,3,3,nat,nat) 1751 real(DP) v(nr1,nr2,nr3,3,3,nat,nat) 1752 real(DP) scal 1753 ! 1754 ! 1755 scal=0.0d0 1756 do j=1,3 1757 do nb=1,nat 1758 do n1=1,nr1 1759 do n2=1,nr2 1760 do n3=1,nr3 1761 scal=scal+u(n1,n2,n3,i,j,na,nb)*v(n1,n2,n3,i,j,na,nb) 1762 enddo 1763 enddo 1764 enddo 1765 enddo 1766 enddo 1767 ! 1768 return 1769 ! 1770end subroutine sp3 1771! 1772!----------------------------------------------------------------------- 1773SUBROUTINE q_gen(nsc,qbid,at_blk,bg_blk,at,bg) 1774 !----------------------------------------------------------------------- 1775 ! generate list of q (qbid) that are G-vectors of the supercell 1776 ! but not of the bulk 1777 ! 1778 USE kinds, ONLY : DP 1779 ! 1780 IMPLICIT NONE 1781 INTEGER :: nsc 1782 REAL(DP) qbid(3,nsc), at_blk(3,3), bg_blk(3,3), at(3,3), bg(3,3) 1783 ! 1784 INTEGER, PARAMETER:: nr1=4, nr2=4, nr3=4, & 1785 nrm=(2*nr1+1)*(2*nr2+1)*(2*nr3+1) 1786 REAL(DP), PARAMETER:: eps=1.0d-7 1787 INTEGER :: i, j, k,i1, i2, i3, idum(nrm), iq 1788 REAL(DP) :: qnorm(nrm), qbd(3,nrm) ,qwork(3), delta 1789 LOGICAL lbho 1790 ! 1791 i = 0 1792 DO i1=-nr1,nr1 1793 DO i2=-nr2,nr2 1794 DO i3=-nr3,nr3 1795 i = i + 1 1796 DO j=1,3 1797 qwork(j) = i1*bg(j,1) + i2*bg(j,2) + i3*bg(j,3) 1798 END DO ! j 1799 ! 1800 qnorm(i) = qwork(1)**2 + qwork(2)**2 + qwork(3)**2 1801 ! 1802 DO j=1,3 1803 ! 1804 qbd(j,i) = at_blk(1,j)*qwork(1) + & 1805 at_blk(2,j)*qwork(2) + & 1806 at_blk(3,j)*qwork(3) 1807 END DO ! j 1808 ! 1809 idum(i) = 1 1810 ! 1811 END DO ! i3 1812 END DO ! i2 1813 END DO ! i1 1814 ! 1815 DO i=1,nrm-1 1816 IF (idum(i).EQ.1) THEN 1817 DO j=i+1,nrm 1818 IF (idum(j).EQ.1) THEN 1819 lbho=.TRUE. 1820 DO k=1,3 1821 delta = qbd(k,i)-qbd(k,j) 1822 lbho = lbho.AND. (ABS(NINT(delta)-delta).LT.eps) 1823 END DO ! k 1824 IF (lbho) THEN 1825 IF(qnorm(i).GT.qnorm(j)) THEN 1826 qbd(1,i) = qbd(1,j) 1827 qbd(2,i) = qbd(2,j) 1828 qbd(3,i) = qbd(3,j) 1829 qnorm(i) = qnorm(j) 1830 END IF 1831 idum(j) = 0 1832 END IF 1833 END IF 1834 END DO ! j 1835 END IF 1836 END DO ! i 1837 ! 1838 iq = 0 1839 DO i=1,nrm 1840 IF (idum(i).EQ.1) THEN 1841 iq=iq+1 1842 qbid(1,iq)= bg_blk(1,1)*qbd(1,i) + & 1843 bg_blk(1,2)*qbd(2,i) + & 1844 bg_blk(1,3)*qbd(3,i) 1845 qbid(2,iq)= bg_blk(2,1)*qbd(1,i) + & 1846 bg_blk(2,2)*qbd(2,i) + & 1847 bg_blk(2,3)*qbd(3,i) 1848 qbid(3,iq)= bg_blk(3,1)*qbd(1,i) + & 1849 bg_blk(3,2)*qbd(2,i) + & 1850 bg_blk(3,3)*qbd(3,i) 1851 END IF 1852 END DO ! i 1853 ! 1854 IF (iq.NE.nsc) CALL errore('q_gen',' probably nr1,nr2,nr3 too small ', iq) 1855 RETURN 1856END SUBROUTINE q_gen 1857! 1858!----------------------------------------------------------------------- 1859SUBROUTINE check_at(at,bg_blk,alat,omega) 1860 !----------------------------------------------------------------------- 1861 ! 1862 USE kinds, ONLY : DP 1863 USE io_global, ONLY : stdout 1864 ! 1865 IMPLICIT NONE 1866 ! 1867 REAL(DP) :: at(3,3), bg_blk(3,3), alat, omega 1868 REAL(DP) :: work(3,3) 1869 INTEGER :: i,j 1870 REAL(DP), PARAMETER :: small=1.d-6 1871 ! 1872 work(:,:) = at(:,:) 1873 CALL cryst_to_cart(3,work,bg_blk,-1) 1874 ! 1875 DO j=1,3 1876 DO i =1,3 1877 IF ( ABS(work(i,j)-NINT(work(i,j))) > small) THEN 1878 WRITE (stdout,'(3f9.4)') work(:,:) 1879 CALL errore ('check_at','at not multiple of at_blk',1) 1880 END IF 1881 END DO 1882 END DO 1883 ! 1884 omega =alat**3 * ABS(at(1,1)*(at(2,2)*at(3,3)-at(3,2)*at(2,3))- & 1885 at(1,2)*(at(2,1)*at(3,3)-at(2,3)*at(3,1))+ & 1886 at(1,3)*(at(2,1)*at(3,2)-at(2,2)*at(3,1))) 1887 ! 1888 RETURN 1889END SUBROUTINE check_at 1890! 1891!----------------------------------------------------------------------- 1892SUBROUTINE set_tau (nat, nat_blk, at, at_blk, tau, tau_blk, & 1893 ityp, ityp_blk, itau_blk) 1894 !----------------------------------------------------------------------- 1895 ! 1896 USE kinds, ONLY : DP 1897 ! 1898 IMPLICIT NONE 1899 INTEGER nat, nat_blk,ityp(nat),ityp_blk(nat_blk), itau_blk(nat) 1900 REAL(DP) at(3,3),at_blk(3,3),tau(3,nat),tau_blk(3,nat_blk) 1901 ! 1902 REAL(DP) bg(3,3), r(3) ! work vectors 1903 INTEGER i,i1,i2,i3,na,na_blk 1904 REAL(DP) small 1905 INTEGER NN1,NN2,NN3 1906 PARAMETER (NN1=8, NN2=8, NN3=8, small=1.d-8) 1907 ! 1908 CALL recips (at(1,1),at(1,2),at(1,3),bg(1,1),bg(1,2),bg(1,3)) 1909 ! 1910 na = 0 1911 ! 1912 DO i1 = -NN1,NN1 1913 DO i2 = -NN2,NN2 1914 DO i3 = -NN3,NN3 1915 r(1) = i1*at_blk(1,1) + i2*at_blk(1,2) + i3*at_blk(1,3) 1916 r(2) = i1*at_blk(2,1) + i2*at_blk(2,2) + i3*at_blk(2,3) 1917 r(3) = i1*at_blk(3,1) + i2*at_blk(3,2) + i3*at_blk(3,3) 1918 CALL cryst_to_cart(1,r,bg,-1) 1919 ! 1920 IF ( r(1).GT.-small .AND. r(1).LT.1.d0-small .AND. & 1921 r(2).GT.-small .AND. r(2).LT.1.d0-small .AND. & 1922 r(3).GT.-small .AND. r(3).LT.1.d0-small ) THEN 1923 CALL cryst_to_cart(1,r,at,+1) 1924 ! 1925 DO na_blk=1, nat_blk 1926 na = na + 1 1927 IF (na.GT.nat) CALL errore('set_tau','too many atoms',na) 1928 tau(1,na) = tau_blk(1,na_blk) + r(1) 1929 tau(2,na) = tau_blk(2,na_blk) + r(2) 1930 tau(3,na) = tau_blk(3,na_blk) + r(3) 1931 ityp(na) = ityp_blk(na_blk) 1932 itau_blk(na) = na_blk 1933 END DO 1934 ! 1935 END IF 1936 ! 1937 END DO 1938 END DO 1939 END DO 1940 ! 1941 IF (na.NE.nat) CALL errore('set_tau','too few atoms: increase NNs',na) 1942 ! 1943 RETURN 1944END SUBROUTINE set_tau 1945! 1946!----------------------------------------------------------------------- 1947SUBROUTINE read_tau & 1948 (nat, nat_blk, ntyp, bg_blk, tau, tau_blk, ityp, itau_blk) 1949 !--------------------------------------------------------------------- 1950 ! 1951 USE kinds, ONLY : DP 1952 USE io_global, ONLY : ionode_id, ionode 1953 USE mp, ONLY : mp_bcast 1954 USE mp_world, ONLY : world_comm 1955 ! 1956 IMPLICIT NONE 1957 ! 1958 INTEGER nat, nat_blk, ntyp, ityp(nat),itau_blk(nat) 1959 REAL(DP) bg_blk(3,3),tau(3,nat),tau_blk(3,nat_blk) 1960 ! 1961 REAL(DP) r(3) ! work vectors 1962 INTEGER i,na,na_blk 1963 ! 1964 REAL(DP) small 1965 PARAMETER ( small = 1.d-6 ) 1966 ! 1967 DO na=1,nat 1968 IF (ionode) READ(5,*) (tau(i,na),i=1,3), ityp(na) 1969 CALL mp_bcast(tau(:,na),ionode_id, world_comm) 1970 CALL mp_bcast(ityp(na),ionode_id, world_comm) 1971 IF (ityp(na).LE.0 .OR. ityp(na) .GT. ntyp) & 1972 CALL errore('read_tau',' wrong atomic type', na) 1973 DO na_blk=1,nat_blk 1974 r(1) = tau(1,na) - tau_blk(1,na_blk) 1975 r(2) = tau(2,na) - tau_blk(2,na_blk) 1976 r(3) = tau(3,na) - tau_blk(3,na_blk) 1977 CALL cryst_to_cart(1,r,bg_blk,-1) 1978 IF (ABS( r(1)-NINT(r(1)) ) .LT. small .AND. & 1979 ABS( r(2)-NINT(r(2)) ) .LT. small .AND. & 1980 ABS( r(3)-NINT(r(3)) ) .LT. small ) THEN 1981 itau_blk(na) = na_blk 1982 go to 999 1983 END IF 1984 END DO 1985 CALL errore ('read_tau',' wrong atomic position ', na) 1986999 CONTINUE 1987 END DO 1988 ! 1989 RETURN 1990END SUBROUTINE read_tau 1991! 1992!----------------------------------------------------------------------- 1993SUBROUTINE write_tau(fltau,nat,tau,ityp) 1994 !----------------------------------------------------------------------- 1995 ! 1996 USE kinds, ONLY : DP 1997 USE io_global, ONLY : ionode 1998 ! 1999 IMPLICIT NONE 2000 ! 2001 INTEGER nat, ityp(nat) 2002 REAL(DP) tau(3,nat) 2003 CHARACTER(LEN=*) fltau 2004 ! 2005 INTEGER i,na 2006 ! 2007 IF (.NOT.ionode) RETURN 2008 OPEN (unit=4,file=fltau, status='new') 2009 DO na=1,nat 2010 WRITE(4,'(3(f12.6),i3)') (tau(i,na),i=1,3), ityp(na) 2011 END DO 2012 CLOSE (4) 2013 ! 2014 RETURN 2015END SUBROUTINE write_tau 2016! 2017!----------------------------------------------------------------------- 2018SUBROUTINE gen_qpoints (ibrav, at_, bg_, nat, tau, ityp, nk1, nk2, nk3, & 2019 nqx, nq, q, nosym, wk) 2020 !----------------------------------------------------------------------- 2021 ! 2022 USE kinds, ONLY : DP 2023 USE cell_base, ONLY : at, bg 2024 USE symm_base, ONLY : set_sym_bl, find_sym, s, irt, nsym, & 2025 nrot, t_rev, time_reversal, sname 2026 USE ktetra, ONLY : tetra_init 2027 ! 2028 IMPLICIT NONE 2029 ! input 2030 INTEGER :: ibrav, nat, nk1, nk2, nk3, ityp(*) 2031 REAL(DP) :: at_(3,3), bg_(3,3), tau(3,nat) 2032 LOGICAL :: nosym 2033 ! output 2034 INTEGER :: nqx, nq 2035 REAL(DP) :: q(3,nqx), wk(nqx) 2036 ! local 2037 REAL(DP) :: xqq(3), mdum(3,nat) 2038 LOGICAL :: magnetic_sym=.FALSE., skip_equivalence=.FALSE. 2039 ! 2040 time_reversal = .true. 2041 if (nosym) time_reversal = .false. 2042 t_rev(:) = 0 2043 xqq (:) =0.d0 2044 at = at_ 2045 bg = bg_ 2046 CALL set_sym_bl ( ) 2047 ! 2048 if (nosym) then 2049 nrot = 1 2050 nsym = 1 2051 endif 2052 CALL kpoint_grid ( nrot, time_reversal, skip_equivalence, s, t_rev, bg, nqx, & 2053 0,0,0, nk1,nk2,nk3, nq, q, wk) 2054 ! 2055 CALL find_sym ( nat, tau, ityp, .not.time_reversal, mdum ) 2056 ! 2057 CALL irreducible_BZ (nrot, s, nsym, time_reversal, magnetic_sym, & 2058 at, bg, nqx, nq, q, wk, t_rev) 2059 ! 2060 CALL tetra_init (nsym, s, time_reversal, t_rev, at, bg, nqx, 0, 0, 0, & 2061 nk1, nk2, nk3, nq, q) 2062 ! 2063 RETURN 2064END SUBROUTINE gen_qpoints 2065! 2066!--------------------------------------------------------------------- 2067SUBROUTINE a2Fdos & 2068 (nat, nq, nr1, nr2, nr3, ibrav, at, bg, tau, alat, & 2069 nsc, nat_blk, at_blk, bg_blk, itau_blk, omega_blk, rws, nrws, & 2070 dos, Emin, DeltaE, ndos, asr, q, freq,fd, wq ) 2071 !----------------------------------------------------------------------- 2072 ! 2073 USE kinds, ONLY : DP 2074 USE io_global, ONLY : ionode, ionode_id 2075 USE mp, ONLY : mp_bcast 2076 USE mp_world, ONLY : world_comm 2077 USE mp_images, ONLY : intra_image_comm 2078 USE ifconstants, ONLY : zeu, tau_blk 2079 USE constants, ONLY : pi, RY_TO_THZ 2080 USE constants, ONLY : K_BOLTZMANN_RY 2081 USE el_phon, ONLY : el_ph_nsigma 2082 ! 2083 IMPLICIT NONE 2084 ! 2085 INTEGER, INTENT(in) :: nat, nq, nr1, nr2, nr3, ibrav, ndos 2086 LOGICAL, INTENT(in) :: dos,fd 2087 CHARACTER(LEN=*), INTENT(IN) :: asr 2088 REAL(DP), INTENT(in) :: freq(3*nat,nq), q(3,nq), wq(nq), at(3,3), bg(3,3), & 2089 tau(3,nat), alat, Emin, DeltaE 2090 ! 2091 INTEGER, INTENT(in) :: nsc, nat_blk, itau_blk(nat), nrws 2092 REAL(DP), INTENT(in) :: rws(0:3,nrws), at_blk(3,3), bg_blk(3,3), omega_blk 2093 ! 2094 REAL(DP), ALLOCATABLE :: gamma(:,:), frcg(:,:,:,:,:,:,:) 2095 COMPLEX(DP), ALLOCATABLE :: gam(:,:,:,:), gam_blk(:,:,:,:), z(:,:) 2096 2097 real(DP) :: lambda, dos_a2F(3*nat), temp, dos_ee(el_ph_nsigma), dos_tot, & 2098 deg(el_ph_nsigma), fermi(el_ph_nsigma), E 2099 real(DP), parameter :: eps_w2 = 0.0000001d0 2100 integer :: isig, ifn, n, m, na, nb, nc, nu, nmodes, & 2101 i,j,k, ngauss, jsig, p1, p2, p3, filea2F, ios 2102 character(len=256) :: name 2103 real(DP), external :: dos_gam 2104 CHARACTER(LEN=6) :: int_to_char 2105 ! 2106 ! 2107 nmodes = 3*nat 2108 do isig=1,el_ph_nsigma 2109 filea2F = 60 + isig 2110 name= 'elph_dir/a2Fmatdyn.' // TRIM(int_to_char(filea2F)) 2111 IF (ionode) open(unit=filea2F, file=TRIM(name), & 2112 STATUS = 'unknown', IOSTAT=ios) 2113 CALL mp_bcast(ios, ionode_id, intra_image_comm) 2114 IF (ios /= 0) CALL errore('a2Fdos','problem opening file'//TRIM(name),1) 2115 IF (ionode) & 2116 READ(filea2F,*) deg(isig), fermi(isig), dos_ee(isig) 2117 ENDDO 2118 call mp_bcast(deg, ionode_id, intra_image_comm) 2119 call mp_bcast(fermi, ionode_id, intra_image_comm) 2120 call mp_bcast(dos_ee, ionode_id, intra_image_comm) 2121 ! 2122 IF (ionode) THEN 2123 IF(dos) then 2124 open(unit=400,file='lambda',status='unknown') 2125 write(400,*) 2126 write(400,*) ' Electron-phonon coupling constant, lambda ' 2127 write(400,*) 2128 ELSE 2129 open (unit=20,file='gam.lines' ,status='unknown') 2130 write(20,*) 2131 write(20,*) ' Gamma lines for all modes [THz] ' 2132 write(20,*) 2133 write(6,*) 2134 write(6,*) ' Gamma lines for all modes [Rydberg] ' 2135 write(6,*) 2136 ENDIF 2137 ENDIF 2138 ! 2139 ALLOCATE ( frcg(nr1,nr2,nr3,3,3,nat,nat) ) 2140 ALLOCATE ( gamma(3*nat,nq), gam(3,3,nat,nat), gam_blk(3,3,nat_blk,nat_blk) ) 2141 ALLOCATE ( z(3*nat,3*nat) ) 2142 ! 2143 frcg(:,:,:,:,:,:,:) = 0.d0 2144 DO isig = 1, el_ph_nsigma 2145 filea2F = 60 + isig 2146 CALL readfg ( filea2F, nr1, nr2, nr3, nat, frcg ) 2147 ! 2148 if ( asr /= 'no') then 2149 CALL set_asr (asr, nr1, nr2, nr3, frcg, zeu, nat_blk, ibrav, tau_blk) 2150 endif 2151 ! 2152 IF (ionode) open(unit=300,file='dyna2F',status='old') 2153 ! 2154 do n = 1 ,nq 2155 gam(:,:,:,:) = (0.d0, 0.d0) 2156 IF (ionode) THEN 2157 read(300,*) 2158 do na=1,nmodes 2159 read(300,*) (z(na,m),m=1,nmodes) 2160 end do ! na 2161 ENDIF 2162 CALL mp_bcast(z, ionode_id, world_comm) 2163 2164 ! 2165 CALL setgam (q(1,n), gam, nat, at, bg, tau, itau_blk, nsc, alat, & 2166 gam_blk, nat_blk, at_blk,bg_blk,tau_blk, omega_blk, & 2167 frcg, nr1,nr2,nr3, rws, nrws, fd) 2168 ! 2169 ! here multiply dyn*gam*dyn for gamma and divide by w2 for lambda at given q 2170 ! 2171 do nc = 1, nat 2172 do k =1, 3 2173 p1 = (nc-1)*3+k 2174 nu = p1 2175 gamma(nu,n) = 0.0d0 2176 do i=1,3 2177 do na=1,nat 2178 p2 = (na-1)*3+i 2179 do j=1,3 2180 do nb=1,nat 2181 p3 = (nb-1)*3+j 2182 gamma(nu,n) = gamma(nu,n) + DBLE(conjg(z(p2,p1)) * & 2183 gam(i,j,na,nb) * z(p3,p1)) 2184 enddo ! nb 2185 enddo ! j 2186 enddo ! na 2187 enddo !i 2188 gamma(nu,n) = gamma(nu,n) * pi / 2.0d0 2189 enddo ! k 2190 enddo !nc 2191 ! 2192 ! 2193 EndDo !nq all points in BZ 2194 IF (ionode) close(300) ! file with dyn vectors 2195 ! 2196 ! after we know gamma(q) and lambda(q) calculate DOS(omega) for spectrum a2F 2197 ! 2198 if(dos.and.ionode) then 2199 ! 2200 name='a2F.dos'//int_to_char(isig) 2201 ifn = 200 + isig 2202 open (ifn,file=TRIM(name),status='unknown',form='formatted') 2203 write(ifn,*) 2204 write(ifn,*) '# Eliashberg function a2F (per both spin)' 2205 write(ifn,*) '# frequencies in Rydberg ' 2206 write(ifn,*) '# DOS normalized to E in Rydberg: a2F_total, a2F(mode) ' 2207 write(ifn,*) 2208 ! 2209 ! correction for small frequencies 2210 ! 2211 do n = 1, nq 2212 do i = 1, nmodes 2213 if (freq(i,n).LE.eps_w2) then 2214 gamma(i,n) = 0.0d0 2215 endif 2216 enddo 2217 enddo 2218 ! 2219 lambda = 0.0d0 2220 do n= 1, ndos 2221 ! 2222 E = Emin + (n-1)*DeltaE + 0.5d0*DeltaE 2223 dos_tot = 0.0d0 2224 do j=1,nmodes 2225 ! 2226 dos_a2F(j) = dos_gam(nmodes, nq, j, gamma, freq, E) 2227 dos_a2F(j) = dos_a2F(j) / dos_ee(isig) / 2.d0 / pi 2228 dos_tot = dos_tot + dos_a2F(j) 2229 ! 2230 enddo 2231 lambda = lambda + 2.d0 * dos_tot/E * DeltaE 2232 write (ifn, '(3X,1000E16.6)') E, dos_tot, dos_a2F(1:nmodes) 2233 enddo !ndos 2234 write(ifn,*) " lambda =",lambda,' Delta = ',DeltaE 2235 close (ifn) 2236 ! 2237 ! lambda from alternative way, simple sum. 2238 ! Also Omega_ln is computed 2239 ! 2240 lambda = 0.0_dp 2241 E = 0.0_dp 2242 do n = 1, nq 2243 lambda = lambda & 2244 & + sum(gamma(1:nmodes,n)/freq(1:nmodes,n)**2, & 2245 & freq(1:nmodes,n) > 1.0e-5_dp) * wq(n) 2246 E = E & 2247 & + sum(log(freq(1:nmodes,n)) * gamma(1:nmodes,n)/freq(1:nmodes,n)**2, & 2248 & freq(1:nmodes,n) > 1.0e-5_dp) * wq(n) 2249 end do 2250 E = exp(E / lambda) / K_BOLTZMANN_RY 2251 lambda = lambda / (dos_ee(isig) * pi) 2252 write(400,'(" Broadening ",F8.4," lambda ",F12.4," dos(Ef)",F8.4," omega_ln [K]",F12.4)') & 2253 deg(isig),lambda, dos_ee(isig), E 2254 ! 2255 endif !dos 2256 ! 2257 ! OUTPUT 2258 ! 2259 if(.not.dos.and.ionode) then 2260 write(20,'(" Broadening ",F8.4)') deg(isig) 2261 write( 6,'(" Broadening ",F8.4)') deg(isig) 2262 do n=1, nq 2263 write(20,'(3x,i5)') n 2264 write( 6,'(3x,i5)') n 2265 write(20,'(9F8.4)') (gamma(i,n)*RY_TO_THZ,i=1,3*nat) 2266 write( 6,'(6F12.9)') (gamma(i,n),i=1,3*nat) 2267! 2268! write also in a format that can be read by plotband 2269! 2270 WRITE(200+isig, '(10x,3f10.6)') q(1,n), q(2,n), q(3,n) 2271! 2272! output in GHz 2273! 2274 WRITE(200+isig, '(6f10.4)') (gamma(nu,n)*RY_TO_THZ*1000.0_DP, & 2275 nu=1,3*nat) 2276 end do 2277 endif 2278 ! 2279 ENDDO !isig 2280 ! 2281 DEALLOCATE (z, frcg, gamma, gam, gam_blk ) 2282 ! 2283 IF (ionode) THEN 2284 close(400) !lambda 2285 close(20) 2286 ENDIF 2287 ! 2288 ! 2289 RETURN 2290END SUBROUTINE a2Fdos 2291! 2292!----------------------------------------------------------------------- 2293subroutine setgam (q, gam, nat, at,bg,tau,itau_blk,nsc,alat, & 2294 & gam_blk, nat_blk, at_blk,bg_blk,tau_blk,omega_blk, & 2295 & frcg, nr1,nr2,nr3, rws,nrws, fd) 2296 !----------------------------------------------------------------------- 2297 ! compute the dynamical matrix (the analytic part only) 2298 ! 2299 USE kinds, ONLY : DP 2300 USE constants, ONLY : tpi 2301 implicit none 2302 ! 2303 ! I/O variables 2304 ! 2305 integer :: nr1, nr2, nr3, nat, nat_blk, & 2306 nsc, nrws, itau_blk(nat) 2307 real(DP) :: q(3), tau(3,nat), at(3,3), bg(3,3), alat, rws(0:3,nrws) 2308 real(DP) :: tau_blk(3,nat_blk), at_blk(3,3), bg_blk(3,3), omega_blk, & 2309 frcg(nr1,nr2,nr3,3,3,nat_blk,nat_blk) 2310 COMPLEX(DP) :: gam_blk(3,3,nat_blk,nat_blk),f_of_q(3,3,nat,nat) 2311 COMPLEX(DP) :: gam(3,3,nat,nat) 2312 LOGICAL :: fd 2313 ! 2314 ! local variables 2315 ! 2316 real(DP) :: arg 2317 complex(DP) :: cfac(nat) 2318 integer :: i,j,k, na,nb, na_blk, nb_blk, iq 2319 real(DP) :: qp(3), qbid(3,nsc) ! automatic array 2320 ! 2321 ! 2322 call q_gen(nsc,qbid,at_blk,bg_blk,at,bg) 2323 ! 2324 f_of_q=(0.0_DP,0.0_DP) 2325 do iq=1,nsc 2326 ! 2327 do k=1,3 2328 qp(k)= q(k) + qbid(k,iq) 2329 end do 2330 ! 2331 gam_blk(:,:,:,:) = (0.d0,0.d0) 2332 CALL frc_blk (gam_blk,qp,tau_blk,nat_blk, & 2333 nr1,nr2,nr3,frcg,at_blk,bg_blk,rws,nrws,f_of_q,fd) 2334 ! 2335 do na=1,nat 2336 na_blk = itau_blk(na) 2337 do nb=1,nat 2338 nb_blk = itau_blk(nb) 2339 ! 2340 arg = tpi * ( qp(1) * ( (tau(1,na)-tau_blk(1,na_blk)) - & 2341 (tau(1,nb)-tau_blk(1,nb_blk)) ) + & 2342 qp(2) * ( (tau(2,na)-tau_blk(2,na_blk)) - & 2343 (tau(2,nb)-tau_blk(2,nb_blk)) ) + & 2344 qp(3) * ( (tau(3,na)-tau_blk(3,na_blk)) - & 2345 (tau(3,nb)-tau_blk(3,nb_blk)) ) ) 2346 ! 2347 cfac(nb) = CMPLX(cos(arg),sin(arg), kind=dp)/nsc 2348 ! 2349 end do ! nb 2350 do nb=1,nat 2351 do i=1,3 2352 do j=1,3 2353 nb_blk = itau_blk(nb) 2354 gam(i,j,na,nb) = gam(i,j,na,nb) + cfac(nb) * & 2355 gam_blk(i,j,na_blk,nb_blk) 2356 end do ! j 2357 end do ! i 2358 end do ! nb 2359 end do ! na 2360 ! 2361 end do ! iq 2362 ! 2363 return 2364end subroutine setgam 2365! 2366!-------------------------------------------------------------------- 2367function dos_gam (nbndx, nq, jbnd, gamma, et, ef) 2368 !-------------------------------------------------------------------- 2369 ! calculates weights with the tetrahedron method (Bloechl version) 2370 ! this subroutine is based on tweights.f90 belonging to PW 2371 ! it calculates a2F on the surface of given frequency <=> histogram 2372 ! Band index means the frequency mode here 2373 ! and "et" means the frequency(mode,q-point) 2374 ! 2375 USE kinds, ONLY: DP 2376 USE parameters 2377 USE ktetra, ONLY : ntetra, tetra 2378 implicit none 2379 ! 2380 integer :: nq, nbndx, jbnd 2381 real(DP) :: et(nbndx,nq), gamma(nbndx,nq), func 2382 2383 real(DP) :: ef 2384 real(DP) :: e1, e2, e3, e4, c1, c2, c3, c4, etetra(4) 2385 integer :: ik, ibnd, nt, nk, ns, i, ik1, ik2, ik3, ik4, itetra(4) 2386 2387 real(DP) :: f12,f13,f14,f23,f24,f34, f21,f31,f41,f42,f32,f43 2388 real(DP) :: P1,P2,P3,P4, G, o13, Y1,Y2,Y3,Y4, eps,vol, Tint 2389 real(DP) :: dos_gam 2390 2391 Tint = 0.0d0 2392 o13 = 1.0_dp/3.0_dp 2393 eps = 1.0d-14 2394 vol = 1.0d0/ntetra 2395 P1 = 0.0_dp 2396 P2 = 0.0_dp 2397 P3 = 0.0_dp 2398 P4 = 0.0_dp 2399 do nt = 1, ntetra 2400 ibnd = jbnd 2401 ! 2402 ! etetra are the energies at the vertexes of the nt-th tetrahedron 2403 ! 2404 do i = 1, 4 2405 etetra(i) = et(ibnd, tetra(i,nt)) 2406 enddo 2407 itetra(1) = 0 2408 call hpsort (4,etetra,itetra) 2409 ! 2410 ! ...sort in ascending order: e1 < e2 < e3 < e4 2411 ! 2412 e1 = etetra (1) 2413 e2 = etetra (2) 2414 e3 = etetra (3) 2415 e4 = etetra (4) 2416 ! 2417 ! kp1-kp4 are the irreducible k-points corresponding to e1-e4 2418 ! 2419 ik1 = tetra(itetra(1),nt) 2420 ik2 = tetra(itetra(2),nt) 2421 ik3 = tetra(itetra(3),nt) 2422 ik4 = tetra(itetra(4),nt) 2423 Y1 = gamma(ibnd,ik1)/et(ibnd,ik1) 2424 Y2 = gamma(ibnd,ik2)/et(ibnd,ik2) 2425 Y3 = gamma(ibnd,ik3)/et(ibnd,ik3) 2426 Y4 = gamma(ibnd,ik4)/et(ibnd,ik4) 2427 2428 IF ( e3 < ef .and. ef < e4) THEN 2429 2430 f14 = (ef-e4)/(e1-e4) 2431 f24 = (ef-e4)/(e2-e4) 2432 f34 = (ef-e4)/(e3-e4) 2433 2434 G = 3.0_dp * f14 * f24 * f34 / (e4-ef) 2435 P1 = f14 * o13 2436 P2 = f24 * o13 2437 P3 = f34 * o13 2438 P4 = (3.0_dp - f14 - f24 - f34 ) * o13 2439 2440 ELSE IF ( e2 < ef .and. ef < e3 ) THEN 2441 2442 f13 = (ef-e3)/(e1-e3) 2443 f31 = 1.0_dp - f13 2444 f14 = (ef-e4)/(e1-e4) 2445 f41 = 1.0_dp-f14 2446 f23 = (ef-e3)/(e2-e3) 2447 f32 = 1.0_dp - f23 2448 f24 = (ef-e4)/(e2-e4) 2449 f42 = 1.0_dp - f24 2450 2451 G = 3.0_dp * (f23*f31 + f32*f24) 2452 P1 = f14 * o13 + f13*f31*f23 / G 2453 P2 = f23 * o13 + f24*f24*f32 / G 2454 P3 = f32 * o13 + f31*f31*f23 / G 2455 P4 = f41 * o13 + f42*f24*f32 / G 2456 G = G / (e4-e1) 2457 2458 ELSE IF ( e1 < ef .and. ef < e2 ) THEN 2459 2460 f12 = (ef-e2)/(e1-e2) 2461 f21 = 1.0_dp - f12 2462 f13 = (ef-e3)/(e1-e3) 2463 f31 = 1.0_dp - f13 2464 f14 = (ef-e4)/(e1-e4) 2465 f41 = 1.0_dp - f14 2466 2467 G = 3.0_dp * f21 * f31 * f41 / (ef-e1) 2468 P1 = o13 * (f12 + f13 + f14) 2469 P2 = o13 * f21 2470 P3 = o13 * f31 2471 P4 = o13 * f41 2472 2473 ELSE 2474 2475 G = 0.0_dp 2476 2477 END IF 2478 2479 Tint = Tint + G * (Y1*P1 + Y2*P2 + Y3*P3 + Y4*P4) * vol 2480 2481 enddo ! ntetra 2482 2483 2484 dos_gam = Tint !2 because DOS_ee is per 1 spin 2485 2486 return 2487end function dos_gam 2488! 2489! 2490!----------------------------------------------------------------------- 2491subroutine readfg ( ifn, nr1, nr2, nr3, nat, frcg ) 2492 !----------------------------------------------------------------------- 2493 ! 2494 USE kinds, ONLY : DP 2495 USE io_global, ONLY : ionode, ionode_id, stdout 2496 USE mp, ONLY : mp_bcast 2497 USE mp_world, ONLY : world_comm 2498 implicit none 2499 ! I/O variable 2500 integer, intent(in) :: nr1,nr2,nr3, nat 2501 real(DP), intent(out) :: frcg(nr1,nr2,nr3,3,3,nat,nat) 2502 ! local variables 2503 integer i, j, na, nb, m1,m2,m3, ifn 2504 integer ibid, jbid, nabid, nbbid, m1bid,m2bid,m3bid 2505 ! 2506 ! 2507 IF (ionode) READ (ifn,*) m1, m2, m3 2508 CALL mp_bcast(m1, ionode_id, world_comm) 2509 CALL mp_bcast(m2, ionode_id, world_comm) 2510 CALL mp_bcast(m3, ionode_id, world_comm) 2511 if ( m1 /= nr1 .or. m2 /= nr2 .or. m3 /= nr3) & 2512 call errore('readfg','inconsistent nr1, nr2, nr3 read',1) 2513 do i=1,3 2514 do j=1,3 2515 do na=1,nat 2516 do nb=1,nat 2517 IF (ionode) read (ifn,*) ibid, jbid, nabid, nbbid 2518 CALL mp_bcast(ibid, ionode_id, world_comm) 2519 CALL mp_bcast(jbid, ionode_id, world_comm) 2520 CALL mp_bcast(nabid, ionode_id, world_comm) 2521 CALL mp_bcast(nbbid, ionode_id, world_comm) 2522 2523 if(i.ne.ibid.or.j.ne.jbid.or.na.ne.nabid.or.nb.ne.nbbid) then 2524 write(stdout,*) i,j,na,nb,' <> ', ibid, jbid, nabid, nbbid 2525 call errore ('readfG','error in reading',1) 2526 else 2527 IF (ionode) read (ifn,*) (((m1bid, m2bid, m3bid, & 2528 frcg(m1,m2,m3,i,j,na,nb), & 2529 m1=1,nr1),m2=1,nr2),m3=1,nr3) 2530 endif 2531 CALL mp_bcast(frcg(:,:,:,i,j,na,nb), ionode_id, world_comm) 2532 end do 2533 end do 2534 end do 2535 end do 2536 ! 2537 IF (ionode) close(ifn) 2538 ! 2539 return 2540end subroutine readfg 2541! 2542! 2543SUBROUTINE find_representations_mode_q ( nat, ntyp, xq, w2, u, tau, ityp, & 2544 amass, num_rap_mode, nspin_mag ) 2545 2546 USE kinds, ONLY : DP 2547 USE cell_base, ONLY : at, bg 2548 USE symm_base, ONLY : s, sr, ft, irt, nsym, nrot, t_rev, time_reversal,& 2549 sname, copy_sym, s_axis_to_cart 2550 2551 IMPLICIT NONE 2552 INTEGER, INTENT(IN) :: nat, ntyp, nspin_mag 2553 REAL(DP), INTENT(IN) :: xq(3), amass(ntyp), tau(3,nat) 2554 REAL(DP), INTENT(IN) :: w2(3*nat) 2555 INTEGER, INTENT(IN) :: ityp(nat) 2556 COMPLEX(DP), INTENT(IN) :: u(3*nat,3*nat) 2557 INTEGER, INTENT(OUT) :: num_rap_mode(3*nat) 2558 REAL(DP) :: gi (3, 48), gimq (3), sr_is(3,3,48), rtau(3,48,nat) 2559 INTEGER :: irotmq, nsymq, nsym_is, isym, i, ierr 2560 LOGICAL :: minus_q, search_sym, sym(48), magnetic_sym 2561! 2562! find the small group of q 2563! 2564 time_reversal=.TRUE. 2565 IF (.NOT.time_reversal) minus_q=.FALSE. 2566 2567 sym(1:nsym)=.true. 2568 call smallg_q (xq, 0, at, bg, nsym, s, sym, minus_q) 2569 nsymq=copy_sym(nsym,sym ) 2570 call s_axis_to_cart () 2571 CALL set_giq (xq,s,nsymq,nsym,irotmq,minus_q,gi,gimq) 2572! 2573! if the small group of q is non symmorphic, 2574! search the symmetries only if there are no G such that Sq -> q+G 2575! 2576 search_sym=.TRUE. 2577 IF ( ANY ( ABS(ft(:,1:nsymq)) > 1.0d-8 ) ) THEN 2578 DO isym=1,nsymq 2579 search_sym=( search_sym.and.(abs(gi(1,isym))<1.d-8).and. & 2580 (abs(gi(2,isym))<1.d-8).and. & 2581 (abs(gi(3,isym))<1.d-8) ) 2582 END DO 2583 END IF 2584! 2585! Set the representations tables of the small group of q and 2586! find the mode symmetry 2587! 2588 IF (search_sym) THEN 2589 magnetic_sym=(nspin_mag==4) 2590 CALL prepare_sym_analysis(nsymq,sr,t_rev,magnetic_sym) 2591 sym (1:nsym) = .TRUE. 2592 CALL sgam_lr (at, bg, nsym, s, irt, tau, rtau, nat) 2593 CALL find_mode_sym_new (u, w2, tau, nat, nsymq, s, sr, irt, xq, & 2594 rtau, amass, ntyp, ityp, 1, .FALSE., .FALSE., num_rap_mode, ierr) 2595 2596 ENDIF 2597 RETURN 2598 END SUBROUTINE find_representations_mode_q 2599