1! 2! Copyright (C) 2001-2016 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! 8Module dynamical 9 ! 10 ! All variables read from file that need dynamical allocation 11 ! 12 USE kinds, ONLY: DP 13 complex(DP), allocatable :: dyn(:,:,:,:) 14 real(DP), allocatable :: tau(:,:), zstar(:,:,:), dchi_dtau(:,:,:,:), & 15 m_loc(:,:) 16 integer, allocatable :: ityp(:) 17 ! 18end Module dynamical 19! 20!----------------------------------------------------------------------- 21subroutine readmat2 ( fildyn, asr, axis, nat, ntyp, atm, & 22 a0, at, omega, amass, eps0, q ) 23!----------------------------------------------------------------------- 24! 25 USE kinds, ONLY: DP 26 USE dynamical 27 ! 28 implicit none 29 character(len=256), intent(in) :: fildyn 30 character(len=10), intent(in) :: asr 31 integer, intent(in) :: axis 32 integer, intent(inout) :: nat, ntyp 33 character(len=3), intent(out) :: atm(ntyp) 34 real(DP), intent(out) :: amass(ntyp), a0, at(3,3), omega, & 35 eps0(3,3), q(3) 36 ! 37 character(len=80) :: line 38 real(DP) :: celldm(6), dyn0r(3,3,2) 39 integer :: ibrav, nt, na, nb, naa, nbb, i, j, k, ios 40 logical :: qfinito, noraman 41 ! 42 ! 43 noraman=.true. 44 open (unit=1,file=fildyn,status='old',form='formatted') 45 read(1,'(a)') line 46 read(1,'(a)') line 47 read(1,*) ntyp,nat,ibrav,celldm 48 ! 49 if (ibrav==0) then 50 read(1,'(a)') line 51 read(1,*) ((at(i,j),i=1,3),j=1,3) 52 end if 53 ! 54 allocate ( dyn (3,3,nat,nat) ) 55 allocate ( dchi_dtau (3,3,3,nat) ) 56 allocate (zstar(3,3,nat) ) 57 allocate ( tau (3,nat) ) 58 allocate (ityp (nat) ) 59 ! 60 call latgen(ibrav,celldm,at(1,1),at(1,2),at(1,3),omega) 61 a0=celldm(1) ! define alat 62 at = at / a0 ! bring at in units of alat 63 64 do nt=1,ntyp 65 read(1,*) i,atm(nt),amass(nt) 66 end do 67 do na=1,nat 68 read(1,*) i,ityp(na), (tau(j,na),j=1,3) 69 end do 70 read(1,'(a)') line 71 read(1,'(a)') line 72 read(1,'(a)') line 73 read(1,'(a)') line 74 read(line(11:80),*) (q(i),i=1,3) 75 qfinito=q(1).ne.0.0 .or. q(2).ne.0.0 .or. q(3).ne.0.0 76 if (qfinito .and. asr .ne. 'no') & 77 call errore('readmat2','Acoustic Sum Rule for q != 0 ?',1) 78 do na = 1,nat 79 do nb = 1,nat 80 read (1,*) naa, nbb 81 if (na.ne.naa .or. nb.ne.nbb) then 82 call errore ('readmat2','mismatch in reading file',1) 83 end if 84 read (1,*) ((dyn0r(i,j,1), dyn0r(i,j,2), j=1,3), i=1,3) 85 dyn(:,:,na,nb) = CMPLX( dyn0r(:,:,1), dyn0r(:,:,2) ,kind=DP) 86 end do 87 end do 88 write(6,'(5x,a)') '...Force constants read' 89 ! 90 if (.not.qfinito) then 91 ios=0 92 read(1,*,iostat=ios) 93 read(1,'(a)',iostat=ios) line 94 if (ios .ne. 0 .or. line(1:23).ne.' Dielectric Tensor:') then 95 write(6,'(5x,a)') '...epsilon and Z* not read (not found on file)' 96 do na=1,nat 97 do j=1,3 98 do i=1,3 99 zstar(i,j,na)=0.0d0 100 end do 101 end do 102 end do 103 do j=1,3 104 do i=1,3 105 eps0(i,j)=0.0d0 106 end do 107 eps0(j,j)=1.0d0 108 end do 109 else 110 read(1,*) 111 read(1,*) ((eps0(i,j), j=1,3), i=1,3) 112 read(1,*) 113 read(1,*) 114 read(1,*) 115 do na = 1,nat 116 read(1,*) 117 read(1,*) ((zstar(i,j,na), j=1,3),i=1,3) 118 end do 119 write(6,'(5x,a)') '...epsilon and Z* read' 120 20 read(1,'(a)',end=10,err=10) line 121 if (line(1:10) == ' Raman') go to 25 122 go to 20 123 25 read(1,*,end=10,err=10) 124 do na = 1,nat 125 do i = 1, 3 126 read(1,*,end=10,err=10) 127 read(1,*,end=10,err=10) & 128 ((dchi_dtau(k,j,i,na), j=1,3), k=1,3) 129 end do 130 end do 131 write(6,'(5x,a)') '...Raman cross sections read' 132 noraman=.false. 13310 continue 134 end if 135 end if 136 if (noraman) dchi_dtau=0.d0 137 138 if(asr.ne.'no') then 139 call set_asr ( asr, axis, nat, tau, dyn, zstar ) 140 endif 141 142 close(unit=1) 143 144 return 145end subroutine readmat2 146! 147!----------------------------------------------------------------------- 148subroutine RamanIR (nat, omega, w2, z, zstar, eps0, dchi_dtau) 149 !----------------------------------------------------------------------- 150 ! 151 ! write IR and Raman cross sections 152 ! on input: z = eigendisplacements (normalized as <z|M|z>) 153 ! zstar = effective charges (units of e) 154 ! dchi_dtau = derivatives of chi wrt atomic displacement 155 ! (units: A^2) 156 USE kinds, ONLY: DP 157 USE constants, ONLY : fpi, BOHR_RADIUS_ANGS, RY_TO_THZ, RY_TO_CMM1, amu_ry 158 implicit none 159 ! input 160 integer, intent(in) :: nat 161 real(DP) omega, w2(3*nat), zstar(3,3,nat), eps0(3,3), & 162 dchi_dtau(3,3,3,nat), chi(3,3) 163 complex(DP) z(3*nat,3*nat) 164 ! local 165 integer na, nu, ipol, jpol, lpol 166 logical noraman 167 real(DP), allocatable :: infrared(:), raman(:,:,:) 168 real(DP):: polar(3), cm1thz, freq, irfac 169 real(DP):: cmfac, alpha, beta2 170 ! 171 ! 172 cm1thz = RY_TO_THZ/RY_TO_CMM1 173 ! 174 ! conversion factor for IR cross sections from 175 ! (Ry atomic units * e^2) to (Debye/A)^2/amu 176 ! 1 Ry mass unit = 2 * mass of one electron = 2 amu 177 ! 1 e = 4.80324x10^(-10) esu = 4.80324 Debye/A 178 ! (1 Debye = 10^(-18) esu*cm = 0.2081928 e*A) 179 ! 180 irfac = 4.80324d0**2/2.d0*amu_ry 181 ! 182 write (6,'(/5x,"Polarizability (A^3 units)")') 183 ! 184 ! correction to molecular polarizabilities from Clausius-Mossotti formula 185 ! (for anisotropic systems epsilon is replaced by its trace) 186 ! 187 cmfac = 3.d0 / ( 2.d0 + (eps0(1,1) + eps0(2,2) + eps0(3,3))/3.d0 ) 188 ! 189 write (6,'(5x,"multiply by",f9.6," for Clausius-Mossotti correction")') cmfac 190 do jpol=1,3 191 do ipol=1,3 192 if (ipol == jpol) then 193 chi(ipol,jpol) = (eps0(ipol,jpol)-1.d0) 194 else 195 chi(ipol,jpol) = eps0(ipol,jpol) 196 end if 197 end do 198 end do 199 do ipol=1,3 200 write (6,'(5x,3f12.6)') (chi(ipol,jpol)*BOHR_RADIUS_ANGS**3*omega/fpi, & 201 jpol=1,3) 202 end do 203 ! 204 allocate(infrared (3*nat)) 205 allocate(raman(3,3,3*nat)) 206 ! 207 noraman=.true. 208 do nu = 1,3*nat 209 do ipol=1,3 210 polar(ipol)=0.0d0 211 end do 212 do na=1,nat 213 do ipol=1,3 214 do jpol=1,3 215 polar(ipol) = polar(ipol) + & 216 zstar(ipol,jpol,na)*z((na-1)*3+jpol,nu) 217 end do 218 end do 219 end do 220 ! 221 infrared(nu) = 2.d0*(polar(1)**2+polar(2)**2+polar(3)**2)*irfac 222 ! 223 do ipol=1,3 224 do jpol=1,3 225 raman(ipol,jpol,nu)=0.0d0 226 do na=1,nat 227 do lpol=1,3 228 raman(ipol,jpol,nu) = raman(ipol,jpol,nu) + & 229 dchi_dtau(ipol,jpol,lpol,na) * z((na-1)*3+lpol,nu) 230 end do 231 end do 232 noraman=noraman .and. abs(raman(ipol,jpol,nu)).lt.1.d-12 233 end do 234 end do 235 ! Raman cross sections are in units of bohr^4/(Ry mass unit) 236 end do 237 ! 238 write (6,'(/5x,"IR activities are in (D/A)^2/amu units")') 239 if (noraman) then 240 write (6,'(/"# mode [cm-1] [THz] IR")') 241 else 242 write (6,'(5x,"Raman activities are in A^4/amu units")') 243 write (6,'(5x,"multiply Raman by",f9.6," for Clausius-Mossotti", & 244 & " correction")') cmfac**2 245 write (6,'(/"# mode [cm-1] [THz] IR Raman depol.fact")') 246 end if 247 ! 248 do nu = 1,3*nat 249 ! 250 freq = sqrt(abs(w2(nu)))*RY_TO_CMM1 251 if (w2(nu).lt.0.0) freq = -freq 252 ! 253 ! alpha, beta2: see PRB 54, 7830 (1996) and refs quoted therein 254 ! 255 if (noraman) then 256 write (6,'(i5,f10.2,2f10.4)') & 257 nu, freq, freq*cm1thz, infrared(nu) 258 else 259 alpha = (raman(1,1,nu) + raman(2,2,nu) + raman(3,3,nu))/3.d0 260 beta2 = ( (raman(1,1,nu) - raman(2,2,nu))**2 + & 261 (raman(1,1,nu) - raman(3,3,nu))**2 + & 262 (raman(2,2,nu) - raman(3,3,nu))**2 + 6.d0 * & 263 (raman(1,2,nu)**2 + raman(1,3,nu)**2 + raman(2,3,nu)**2) )/2.d0 264 write (6,'(i5,f10.2,2f10.4,f15.4,f10.4)') & 265 nu, freq, freq*cm1thz, infrared(nu), & 266 (45.d0*alpha**2 + 7.0d0*beta2)*amu_ry, & 267 3.d0*beta2/(45.d0*alpha**2 + 4.0d0*beta2) 268 end if 269 end do 270 ! 271 deallocate (raman) 272 deallocate (infrared) 273 return 274 ! 275end subroutine RamanIR 276! 277!---------------------------------------------------------------------- 278subroutine set_asr ( asr, axis, nat, tau, dyn, zeu ) 279 !----------------------------------------------------------------------- 280 ! 281 ! Impose ASR - refined version by Nicolas Mounet 282 ! 283 USE kinds, ONLY: DP 284 implicit none 285 character(len=10), intent(in) :: asr 286 integer, intent(in) :: axis, nat 287 real(DP), intent(in) :: tau(3,nat) 288 real(DP), intent(inout) :: zeu(3,3,nat) 289 complex(DP), intent(inout) :: dyn(3,3,nat,nat) 290 ! 291 integer :: i,j,n,m,p,k,l,q,r,na, nb, na1, i1, j1 292 real(DP), allocatable:: dynr_new(:,:,:,:,:), zeu_new(:,:,:) 293 real(DP), allocatable :: u(:,:,:,:,:) 294 ! These are the "vectors" associated with the sum rules 295 ! 296 integer u_less(6*3*nat),n_less,i_less 297 ! indices of the vectors u that are not independent to the preceding ones, 298 ! n_less = number of such vectors, i_less = temporary parameter 299 ! 300 integer, allocatable :: ind_v(:,:,:) 301 real(DP), allocatable :: v(:,:) 302 ! These are the "vectors" associated with symmetry conditions, coded by 303 ! indicating the positions (i.e. the four indices) of the non-zero elements 304 ! (there should be only 2 of them) and the value of that element. 305 ! We do so in order to use limit the amount of memory used. 306 ! 307 real(DP), allocatable :: w(:,:,:,:), x(:,:,:,:) 308 real(DP) sum, scal, norm2 309 ! temporary vectors and parameters 310 ! 311 real(DP), allocatable :: zeu_u(:,:,:,:) 312 ! These are the "vectors" associated with the sum rules on effective charges 313 ! 314 integer zeu_less(6*3),nzeu_less,izeu_less 315 ! indices of the vectors zeu_u that are not independent to the preceding 316 ! ones, nzeu_less = number of such vectors, izeu_less = temporary parameter 317 ! 318 real(DP), allocatable :: zeu_w(:,:,:), zeu_x(:,:,:) 319 ! temporary vectors 320 ! 321 ! Initialization 322 ! n is the number of sum rules to be considered (if asr.ne.'simple') 323 ! 'axis' is the rotation axis in the case of a 1D system (i.e. the rotation 324 ! axis is (Ox) if axis='1', (Oy) if axis='2' and (Oz) if axis='3') 325 ! 326 if ( (asr.ne.'simple') .and. (asr.ne.'crystal') .and. (asr.ne.'one-dim') & 327 .and.(asr.ne.'zero-dim')) then 328 call errore('set_asr','invalid Acoustic Sum Rule:' // asr, 1) 329 endif 330 if(asr.eq.'crystal') n=3 331 if(asr.eq.'one-dim') then 332 write(6,'("asr rotation axis in 1D system= ",I4)') axis 333 n=4 334 endif 335 if(asr.eq.'zero-dim') n=6 336 ! 337 ! ASR on effective charges 338 ! 339 if(asr.eq.'simple') then 340 do i=1,3 341 do j=1,3 342 sum=0.0d0 343 do na=1,nat 344 sum = sum + zeu(i,j,na) 345 end do 346 do na=1,nat 347 zeu(i,j,na) = zeu(i,j,na) - sum/nat 348 end do 349 end do 350 end do 351 else 352 ! generating the vectors of the orthogonal of the subspace to project 353 ! the effective charges matrix on 354 ! 355 allocate ( zeu_new(3,3,nat) ) 356 allocate (zeu_u(6*3,3,3,nat) ) 357 zeu_u(:,:,:,:)=0.0d0 358 do i=1,3 359 do j=1,3 360 do na=1,nat 361 zeu_new(i,j,na)=zeu(i,j,na) 362 enddo 363 enddo 364 enddo 365 ! 366 p=0 367 do i=1,3 368 do j=1,3 369 ! These are the 3*3 vectors associated with the 370 ! translational acoustic sum rules 371 p=p+1 372 zeu_u(p,i,j,:)=1.0d0 373 ! 374 enddo 375 enddo 376 ! 377 if (n.eq.4) then 378 do i=1,3 379 ! These are the 3 vectors associated with the 380 ! single rotational sum rule (1D system) 381 p=p+1 382 do na=1,nat 383 zeu_u(p,i,MOD(axis,3)+1,na)=-tau(MOD(axis+1,3)+1,na) 384 zeu_u(p,i,MOD(axis+1,3)+1,na)=tau(MOD(axis,3)+1,na) 385 enddo 386 ! 387 enddo 388 endif 389 ! 390 if (n.eq.6) then 391 do i=1,3 392 do j=1,3 393 ! These are the 3*3 vectors associated with the 394 ! three rotational sum rules (0D system - typ. molecule) 395 p=p+1 396 do na=1,nat 397 zeu_u(p,i,MOD(j,3)+1,na)=-tau(MOD(j+1,3)+1,na) 398 zeu_u(p,i,MOD(j+1,3)+1,na)=tau(MOD(j,3)+1,na) 399 enddo 400 ! 401 enddo 402 enddo 403 endif 404 ! 405 ! Gram-Schmidt orthonormalization of the set of vectors created. 406 ! 407 allocate ( zeu_w(3,3,nat), zeu_x(3,3,nat) ) 408 nzeu_less=0 409 do k=1,p 410 zeu_w(:,:,:)=zeu_u(k,:,:,:) 411 zeu_x(:,:,:)=zeu_u(k,:,:,:) 412 do q=1,k-1 413 r=1 414 do izeu_less=1,nzeu_less 415 if (zeu_less(izeu_less).eq.q) r=0 416 enddo 417 if (r.ne.0) then 418 call sp_zeu(zeu_x,zeu_u(q,:,:,:),nat,scal) 419 zeu_w(:,:,:) = zeu_w(:,:,:) - scal* zeu_u(q,:,:,:) 420 endif 421 enddo 422 call sp_zeu(zeu_w,zeu_w,nat,norm2) 423 if (norm2.gt.1.0d-16) then 424 zeu_u(k,:,:,:) = zeu_w(:,:,:) / DSQRT(norm2) 425 else 426 nzeu_less=nzeu_less+1 427 zeu_less(nzeu_less)=k 428 endif 429 enddo 430 ! 431 ! 432 ! Projection of the effective charge "vector" on the orthogonal of the 433 ! subspace of the vectors verifying the sum rules 434 ! 435 zeu_w(:,:,:)=0.0d0 436 do k=1,p 437 r=1 438 do izeu_less=1,nzeu_less 439 if (zeu_less(izeu_less).eq.k) r=0 440 enddo 441 if (r.ne.0) then 442 zeu_x(:,:,:)=zeu_u(k,:,:,:) 443 call sp_zeu(zeu_x,zeu_new,nat,scal) 444 zeu_w(:,:,:) = zeu_w(:,:,:) + scal*zeu_u(k,:,:,:) 445 endif 446 enddo 447 ! 448 ! Final substraction of the former projection to the initial zeu, to get 449 ! the new "projected" zeu 450 ! 451 zeu_new(:,:,:)=zeu_new(:,:,:) - zeu_w(:,:,:) 452 call sp_zeu(zeu_w,zeu_w,nat,norm2) 453 write(6,'(5x,"Acoustic Sum Rule: || Z*(ASR) - Z*(orig)|| = ",ES15.6)') & 454 SQRT(norm2) 455 ! 456 ! Check projection 457 ! 458 !write(6,'("Check projection of zeu")') 459 !do k=1,p 460 ! zeu_x(:,:,:)=zeu_u(k,:,:,:) 461 ! call sp_zeu(zeu_x,zeu_new,nat,scal) 462 ! if (DABS(scal).gt.1d-10) write(6,'("k= ",I8," zeu_new|zeu_u(k)= ",F15.10)') k,scal 463 !enddo 464 ! 465 do i=1,3 466 do j=1,3 467 do na=1,nat 468 zeu(i,j,na)=zeu_new(i,j,na) 469 enddo 470 enddo 471 enddo 472 deallocate (zeu_w, zeu_x) 473 deallocate (zeu_u) 474 deallocate (zeu_new) 475 endif 476 ! 477 ! ASR on dynamical matrix 478 ! 479 if(asr.eq.'simple') then 480 do i=1,3 481 do j=1,3 482 do na=1,nat 483 sum=0.0d0 484 do nb=1,nat 485 if (na.ne.nb) sum=sum + DBLE (dyn(i,j,na,nb)) 486 end do 487 dyn(i,j,na,na) = CMPLX(-sum, 0.d0,kind=DP) 488 end do 489 end do 490 end do 491 ! 492 else 493 ! generating the vectors of the orthogonal of the subspace to project 494 ! the dyn. matrix on 495 ! 496 allocate (u(6*3*nat,3,3,nat,nat)) 497 allocate (dynr_new(2,3,3,nat,nat)) 498 u(:,:,:,:,:)=0.0d0 499 do i=1,3 500 do j=1,3 501 do na=1,nat 502 do nb=1,nat 503 dynr_new(1,i,j,na,nb) = DBLE (dyn(i,j,na,nb) ) 504 dynr_new(2,i,j,na,nb) =AIMAG (dyn(i,j,na,nb) ) 505 enddo 506 enddo 507 enddo 508 enddo 509 ! 510 p=0 511 do i=1,3 512 do j=1,3 513 do na=1,nat 514 ! These are the 3*3*nat vectors associated with the 515 ! translational acoustic sum rules 516 p=p+1 517 do nb=1,nat 518 u(p,i,j,na,nb)=1.0d0 519 enddo 520 ! 521 enddo 522 enddo 523 enddo 524 ! 525 if (n.eq.4) then 526 do i=1,3 527 do na=1,nat 528 ! These are the 3*nat vectors associated with the 529 ! single rotational sum rule (1D system) 530 p=p+1 531 do nb=1,nat 532 u(p,i,axis,na,nb)=0.0d0 533 u(p,i,MOD(axis,3)+1,na,nb)=-tau(MOD(axis+1,3)+1,nb) 534 u(p,i,MOD(axis+1,3)+1,na,nb)=tau(MOD(axis,3)+1,nb) 535 enddo 536 ! 537 enddo 538 enddo 539 endif 540 ! 541 if (n.eq.6) then 542 do i=1,3 543 do j=1,3 544 do na=1,nat 545 ! These are the 3*3*nat vectors associated with the 546 ! three rotational sum rules (0D system - typ. molecule) 547 p=p+1 548 do nb=1,nat 549 u(p,i,j,na,nb)=0.0d0 550 u(p,i,MOD(j,3)+1,na,nb)=-tau(MOD(j+1,3)+1,nb) 551 u(p,i,MOD(j+1,3)+1,na,nb)=tau(MOD(j,3)+1,nb) 552 enddo 553 ! 554 enddo 555 enddo 556 enddo 557 endif 558 ! 559 allocate (ind_v(9*nat*nat,2,4)) 560 allocate (v(9*nat*nat,2)) 561 m=0 562 do i=1,3 563 do j=1,3 564 do na=1,nat 565 do nb=1,nat 566 ! These are the vectors associated with the symmetry constraints 567 q=1 568 l=1 569 do while((l.le.m).and.(q.ne.0)) 570 if ((ind_v(l,1,1).eq.i).and.(ind_v(l,1,2).eq.j).and. & 571 (ind_v(l,1,3).eq.na).and.(ind_v(l,1,4).eq.nb)) q=0 572 if ((ind_v(l,2,1).eq.i).and.(ind_v(l,2,2).eq.j).and. & 573 (ind_v(l,2,3).eq.na).and.(ind_v(l,2,4).eq.nb)) q=0 574 l=l+1 575 enddo 576 if ((i.eq.j).and.(na.eq.nb)) q=0 577 if (q.ne.0) then 578 m=m+1 579 ind_v(m,1,1)=i 580 ind_v(m,1,2)=j 581 ind_v(m,1,3)=na 582 ind_v(m,1,4)=nb 583 v(m,1)=1.0d0/DSQRT(2.0d0) 584 ind_v(m,2,1)=j 585 ind_v(m,2,2)=i 586 ind_v(m,2,3)=nb 587 ind_v(m,2,4)=na 588 v(m,2)=-1.0d0/DSQRT(2.0d0) 589 endif 590 enddo 591 enddo 592 enddo 593 enddo 594 ! 595 ! Gram-Schmidt orthonormalization of the set of vectors created. 596 ! Note that the vectors corresponding to symmetry constraints are already 597 ! orthonormalized by construction. 598 ! 599 allocate ( w(3,3,nat,nat), x(3,3,nat,nat)) 600 n_less=0 601 do k=1,p 602 w(:,:,:,:)=u(k,:,:,:,:) 603 x(:,:,:,:)=u(k,:,:,:,:) 604 do l=1,m 605 ! 606 call sp2(x,v(l,:),ind_v(l,:,:),nat,scal) 607 do r=1,2 608 i=ind_v(l,r,1) 609 j=ind_v(l,r,2) 610 na=ind_v(l,r,3) 611 nb=ind_v(l,r,4) 612 w(i,j,na,nb)=w(i,j,na,nb)-scal*v(l,r) 613 enddo 614 enddo 615 if (k.le.(9*nat)) then 616 na1=MOD(k,nat) 617 if (na1.eq.0) na1=nat 618 j1=MOD((k-na1)/nat,3)+1 619 i1=MOD((((k-na1)/nat)-j1+1)/3,3)+1 620 else 621 q=k-9*nat 622 if (n.eq.4) then 623 na1=MOD(q,nat) 624 if (na1.eq.0) na1=nat 625 i1=MOD((q-na1)/nat,3)+1 626 else 627 na1=MOD(q,nat) 628 if (na1.eq.0) na1=nat 629 j1=MOD((q-na1)/nat,3)+1 630 i1=MOD((((q-na1)/nat)-j1+1)/3,3)+1 631 endif 632 endif 633 do q=1,k-1 634 r=1 635 do i_less=1,n_less 636 if (u_less(i_less).eq.q) r=0 637 enddo 638 if (r.ne.0) then 639 call sp3(x,u(q,:,:,:,:),i1,na1,nat,scal) 640 w(:,:,:,:) = w(:,:,:,:) - scal* u(q,:,:,:,:) 641 endif 642 enddo 643 call sp1(w,w,nat,norm2) 644 if (norm2.gt.1.0d-16) then 645 u(k,:,:,:,:) = w(:,:,:,:) / DSQRT(norm2) 646 else 647 n_less=n_less+1 648 u_less(n_less)=k 649 endif 650 enddo 651 ! 652 ! Projection of the dyn. "vector" on the orthogonal of the 653 ! subspace of the vectors verifying the sum rules and symmetry contraints 654 ! 655 w(:,:,:,:)=0.0d0 656 do l=1,m 657 call sp2(dynr_new(1,:,:,:,:),v(l,:),ind_v(l,:,:),nat,scal) 658 do r=1,2 659 i=ind_v(l,r,1) 660 j=ind_v(l,r,2) 661 na=ind_v(l,r,3) 662 nb=ind_v(l,r,4) 663 w(i,j,na,nb)=w(i,j,na,nb)+scal*v(l,r) 664 enddo 665 enddo 666 do k=1,p 667 r=1 668 do i_less=1,n_less 669 if (u_less(i_less).eq.k) r=0 670 enddo 671 if (r.ne.0) then 672 x(:,:,:,:)=u(k,:,:,:,:) 673 call sp1(x,dynr_new(1,:,:,:,:),nat,scal) 674 w(:,:,:,:) = w(:,:,:,:) + scal* u(k,:,:,:,:) 675 endif 676 enddo 677 ! 678 ! Final substraction of the former projection to the initial dyn, 679 ! to get the new "projected" dyn 680 ! 681 dynr_new(1,:,:,:,:)=dynr_new(1,:,:,:,:) - w(:,:,:,:) 682 call sp1(w,w,nat,norm2) 683 write(6,'(5x,"Acoustic Sum Rule: ||dyn(ASR) - dyn(orig)||= ",ES15.6)') & 684 DSQRT(norm2) 685 ! 686 ! Check projection 687 ! 688 !write(6,'("Check projection")') 689 !do l=1,m 690 ! call sp2(dynr_new(1,:,:,:,:),v(l,:),ind_v(l,:,:),nat,scal) 691 ! if (DABS(scal).gt.1d-10) write(6,'("l= ",I8," dyn|v(l)= ",F15.10)') l,scal 692 !enddo 693 !do k=1,p 694 ! x(:,:,:,:)=u(k,:,:,:,:) 695 ! call sp1(x,dynr_new(1,:,:,:,:),nat,scal) 696 ! if (DABS(scal).gt.1d-10) write(6,'("k= ",I8," dyn|u(k)= ",F15.10)') k,scal 697 !enddo 698 ! 699 deallocate ( w, x ) 700 deallocate ( v ) 701 deallocate ( ind_v ) 702 deallocate ( u ) 703 ! 704 do i=1,3 705 do j=1,3 706 do na=1,nat 707 do nb=1,nat 708 dyn (i,j,na,nb) = & 709 CMPLX(dynr_new(1,i,j,na,nb), dynr_new(2,i,j,na,nb) ,kind=DP) 710 enddo 711 enddo 712 enddo 713 enddo 714 deallocate ( dynr_new ) 715 endif 716 ! 717 return 718end subroutine set_asr 719! 720! 721!---------------------------------------------------------------------- 722subroutine sp_zeu(zeu_u,zeu_v,nat,scal) 723 !----------------------------------------------------------------------- 724 ! 725 ! does the scalar product of two effective charges matrices zeu_u and zeu_v 726 ! (considered as vectors in the R^(3*3*nat) space, and coded in the usual way) 727 ! 728 USE kinds, ONLY: DP 729 implicit none 730 integer i,j,na,nat 731 real(DP) zeu_u(3,3,nat) 732 real(DP) zeu_v(3,3,nat) 733 real(DP) scal 734 ! 735 ! 736 scal=0.0d0 737 do i=1,3 738 do j=1,3 739 do na=1,nat 740 scal=scal+zeu_u(i,j,na)*zeu_v(i,j,na) 741 enddo 742 enddo 743 enddo 744 ! 745 return 746 ! 747end subroutine sp_zeu 748! 749! 750!---------------------------------------------------------------------- 751subroutine sp1(u,v,nat,scal) 752 !----------------------------------------------------------------------- 753 ! 754 ! does the scalar product of two dyn. matrices u and v (considered as 755 ! vectors in the R^(3*3*nat*nat) space, and coded in the usual way) 756 ! 757 USE kinds, ONLY: DP 758 implicit none 759 integer i,j,na,nb,nat 760 real(DP) u(3,3,nat,nat) 761 real(DP) v(3,3,nat,nat) 762 real(DP) scal 763 ! 764 ! 765 scal=0.0d0 766 do i=1,3 767 do j=1,3 768 do na=1,nat 769 do nb=1,nat 770 scal=scal+u(i,j,na,nb)*v(i,j,na,nb) 771 enddo 772 enddo 773 enddo 774 enddo 775 ! 776 return 777 ! 778end subroutine sp1 779! 780!---------------------------------------------------------------------- 781subroutine sp2(u,v,ind_v,nat,scal) 782 !----------------------------------------------------------------------- 783 ! 784 ! does the scalar product of two dyn. matrices u and v (considered as 785 ! vectors in the R^(3*3*nat*nat) space). u is coded in the usual way 786 ! but v is coded as explained when defining the vectors corresponding to the 787 ! symmetry constraints 788 ! 789 USE kinds, ONLY: DP 790 implicit none 791 integer i,nat 792 real(DP) u(3,3,nat,nat) 793 integer ind_v(2,4) 794 real(DP) v(2) 795 real(DP) scal 796 ! 797 ! 798 scal=0.0d0 799 do i=1,2 800 scal=scal+u(ind_v(i,1),ind_v(i,2),ind_v(i,3),ind_v(i,4))*v(i) 801 enddo 802 ! 803 return 804 ! 805end subroutine sp2 806! 807!---------------------------------------------------------------------- 808subroutine sp3(u,v,i,na,nat,scal) 809 !----------------------------------------------------------------------- 810 ! 811 ! like sp1, but in the particular case when u is one of the u(k)%vec 812 ! defined in set_asr (before orthonormalization). In this case most of the 813 ! terms are zero (the ones that are not are characterized by i and na), so 814 ! that a lot of computer time can be saved (during Gram-Schmidt). 815 ! 816 USE kinds, ONLY: DP 817 implicit none 818 integer i,j,na,nb,nat 819 real(DP) u(3,3,nat,nat) 820 real(DP) v(3,3,nat,nat) 821 real(DP) scal 822 ! 823 ! 824 scal=0.0d0 825 do j=1,3 826 do nb=1,nat 827 scal=scal+u(i,j,na,nb)*v(i,j,na,nb) 828 enddo 829 enddo 830 ! 831 return 832 ! 833end subroutine sp3 834! 835!---------------------------------------------------------------------- 836subroutine polar_mode_permittivity( nat, eps0, z, zstar, w2, omega, lplasma) 837 !---------------------------------------------------------------------- 838 839 ! 840 ! Algorithm from Fennie and Rabe, Phys. Rev. B 68, 184111 (2003) 841 ! 842 USE kinds, ONLY: DP 843 USE constants, ONLY : pi, tpi, fpi, eps4, eps8, eps12, & 844 ELECTRON_SI, BOHR_RADIUS_SI, AMU_SI, C_SI, & 845 EPSNOUGHT_SI, AMU_RY, RY_TO_CMM1, RY_TO_THZ 846 implicit none 847 !number of atoms 848 integer, intent(in) :: nat 849 850 !electronic part of the permittivity 851 real(DP), intent(in) :: eps0(3,3) 852 853 !displacement eigenvectors 854 complex(DP), intent(in) :: z(3*nat,3*nat) 855 856 !born effective charges 857 real(DP), intent(in) :: zstar(3,3,nat) 858 859 !square of the phonon frequencies 860 real(DP), intent(in) :: w2(3*nat) 861 862 !cell volume 863 real(DP), intent(in) :: omega 864 865 !calculate effective plasma frequencies 866 logical, intent(in) :: lplasma 867 868 !mode index 869 integer :: imode 870 871 !atom index 872 integer :: iat 873 874 !atom vector component index 875 integer :: iat_component 876 877 !Cartesian direction indices 878 integer :: i, j 879 880 !mode effective charge 881 real(DP) :: meffc(3) 882 883 !total effective plasma frequency 884 real(DP) :: weff_tot, freq 885 886 !polar mode contribution to the permittivity 887 real(DP) :: deps(3,3) 888 889 !combined permittivity 890 real(DP) :: eps_new(3,3) 891 892 !Conversion factor for plasma frequency from Rydberg atomic units to SI 893 real(DP) :: plasma_frequency_si 894 !Conversion factor for permittivity from Rydberg atomic units to SI 895 real(DP) :: permittivity_si 896 897 ! some compiler do not like SQRT in initialization expressions 898 plasma_frequency_si = ELECTRON_SI/sqrt(EPSNOUGHT_SI*BOHR_RADIUS_SI**3*AMU_SI) 899 permittivity_si = plasma_frequency_si**2 / (fpi * pi) 900 901 IF (lplasma) THEN 902 WRITE(6,*) 903 WRITE(6,'("# mode freq Z~*_x Z~*_y Z~*_z & 904 & W_eff deps")') 905 WRITE(6,'("# [cm^-1] [e*Bohr/sqrt(2)] & 906 & [cm^-1] [C^2/J*m^2]")') 907 END IF 908 909 eps_new=eps0 910 911 !Calculate contributions by mode 912 DO imode = 1,3*nat 913 914 ! Calculate the mode effective charge 915 meffc = 0.0d0 916 DO i = 1 , 3 917 DO iat = 1 , nat 918 DO j = 1, 3 919 iat_component = 3*(iat-1) + j 920 921 ! Equation (3) of Finnie and Rabe 922 ! Rydberg units = (e / sqrt(2)) * Bohr 923 meffc(i) = meffc(i) + zstar(i,j,iat)*z(iat_component,imode)* & 924 sqrt(AMU_RY) 925 926 END DO 927 END DO 928 END DO 929 930 ! Calculate the polar mode contribution to the permittivity 931 deps = 0.0d0 932 ! Use only hard modes (frequency^2 > 10^-8 Ry) 933 IF (w2(imode) > eps8) THEN 934 DO i = 1 , 3 935 DO j = 1 , 3 936 937 ! Equation (2) of Finnie and Rabe 938 deps(i,j) = (permittivity_si*eps12**2/omega)*meffc(i)*meffc(j) / & 939 (w2(imode)*RY_TO_THZ**2) 940 941 END DO 942 END DO 943 END IF 944 945 ! Add polar mode contribution to the total permittivity 946 DO i = 1 , 3 947 DO j = 1 , 3 948 eps_new(i,j) = eps_new(i,j) + deps(i,j) 949 END DO 950 END DO 951 952 IF (lplasma) THEN 953 ! Calculate the total effective plasma frequency for the mode 954 weff_tot = 0.0d0 955 DO j = 1, 3 956 weff_tot = weff_tot + meffc(j)*meffc(j) 957 END DO 958 ! Rydberg units = (e / sqrt(2)) / (Bohr * sqrt(2*m_e)) 959 weff_tot = sqrt(weff_tot/omega)/tpi 960 961 !Mode frequency [units of sqrt(Ry)]) 962 freq = sqrt(abs(w2(imode))) 963 IF (w2(imode) < 0.0_DP) freq = -freq 964 965 !write out mode index, mode effective charges, 966 ! mode contribution to permittivity, mode plasma frequency 967 WRITE(6,'(i5,6f14.6)') imode,freq*RY_TO_CMM1,meffc(1),meffc(2),meffc(3), & 968 weff_tot*plasma_frequency_si*eps12*(RY_TO_CMM1 / RY_TO_THZ), & 969 (weff_tot*plasma_frequency_si*eps12)**2/(w2(imode)*RY_TO_THZ**2) 970 END IF 971 END DO 972 973 WRITE(6,*) 974 WRITE(6,'("Electronic dielectric permittivity tensor (F/m units)")') 975 WRITE(6,'(5x,3f12.6)') eps0(1,:) 976 WRITE(6,'(5x,3f12.6)') eps0(2,:) 977 WRITE(6,'(5x,3f12.6)') eps0(3,:) 978 WRITE(6,*) 979 WRITE(6,'(" ... with zone-center polar mode contributions")') 980 WRITE(6,'(5x,3f12.6)') eps_new(1,:) 981 WRITE(6,'(5x,3f12.6)') eps_new(2,:) 982 WRITE(6,'(5x,3f12.6)') eps_new(3,:) 983 WRITE(6,*) 984 985end subroutine polar_mode_permittivity 986