1!----------------------------------------------------------------------- 2program fd_raman 3!----------------------------------------------------------------------- 4 use constants 5 use io_files, ONLY : prefix, tmp_dir, psfile, pseudo_dir 6 use io_global, ONLY : stdout, ionode, ionode_id 7 USE mp_global, ONLY : mp_startup 8 USE environment,ONLY : environment_start 9 USE mp, ONLY : mp_bcast 10 USE cell_base, ONLY : tpiba2, alat,omega, at, bg, ibrav, celldm 11 USE ions_base, ONLY : amass, nat, atm, zv, tau, ntyp => nsp, ityp 12 USE kinds, ONLY : dp 13 USE gvecw, ONLY : ecutwfc 14 USE symm_base, ONLY : nsym, nsym_ns, nsym_na, invsym, s, sr, & 15 t_rev, ft, sname 16 USE symme 17 USE fft_base, ONLY : dfftp 18 19 USE parser, ONLY : field_count, read_line, get_field, parse_unit 20 21 implicit none 22 character(len=9) :: code = 'FD_RAMAN' 23 integer :: ios 24 CHARACTER(LEN=256), EXTERNAL :: trimcheck 25 character(len=200) :: pp_file 26 logical :: uspp_spsi, ascii, single_file, raw, disp_only 27 logical :: needwf=.false. 28 INTEGER :: apol, na, nt 29 integer :: nrx1,nrx2,nrx3,nr1,nr2,nr3,nb,nax,natx,inn 30 real(kind=dp) :: r1(3),r2(3),r3(3),rr(3,3) 31 INTEGER :: nclass_ref ! The number of classes of the point group 32 INTEGER :: isym, ipol 33 REAL (dp) :: ft1, ft2, ft3 34 35 integer :: npol, npol_rm, npol1_rm, npol_eps, npol_zeu 36 integer :: ndiag,noffd,nmod,npol1 37 integer :: i,j,p,k,ii,jj,n 38 real*8, allocatable :: F0(:,:),dechi(:,:,:,:),Fd(:,:,:,:),dechi_u(:,:,:,:) 39 real*8, allocatable :: Fij(:,:,:,:),alpha(:,:,:),u(:,:,:),ui(:,:,:) 40 real*8 :: de, de_raman, de_zeu, de_eps,conv,ry,a0 41 CHARACTER(50) :: filemodes 42 logical :: lalpha,lpuma 43 44 real*8, allocatable :: pol0(:),pol(:,:,:),eps(:,:) 45 real*8, allocatable :: zeta(:,:,:) 46 real*8 :: sum 47 48 CHARACTER(len=2) :: prog ! calling program ( PW, CP, WA ) 49 CHARACTER(len=256) :: input_line 50 CHARACTER(len=80) :: card 51 CHARACTER(len=1), EXTERNAL :: capital 52 LOGICAL :: tend, verbose 53 54 NAMELIST /inputfd/ prefix,npol_rm, npol_eps, npol_zeu, ndiag,noffd,nmod,npol1,lpuma, & 55 de_raman, de_eps, de_zeu, filemodes, verbose 56 57 lalpha = .true. 58 lpuma = .false. 59 npol_rm=4 60 npol1 = 2 61 npol_eps=2 62 npol_zeu=2 63 ndiag=3 64 noffd=3 65 de_raman=1.0 66 de_eps=1.0 67 de_zeu=1.0 68 filemodes=' ' 69 verbose=.false. 70 71 ! define conversion constant 72 conv=BOHR_RADIUS_ANGS**2 73 74 CALL mp_startup ( ) 75 CALL environment_start ( code ) 76 77 IF ( ionode ) THEN 78 CALL input_from_file ( ) 79 READ(5,inputfd,IOSTAT=ios) 80 endif 81 if (filemodes .eq. ' ') lalpha=.false. 82 83 !reading the xml file 84 call read_file_new ( needwf ) 85 86 if (ionode) then 87 write(6,*) '**************************************************' 88 write(6,*) '* Info from the preceding pw.x run: *' 89 write(6,*) '**************************************************' 90 write(6,*) '' 91 write(6,*) ' prefix= ',trim(prefix) 92 write(6,*) ' outdir= ',trim(tmp_dir) 93 write(6,*) ' ecutwfc= ',ecutwfc, 'Ry' 94 95 WRITE( stdout, 199) ibrav, alat, omega, nat, ntyp 96 199 FORMAT(5X, & 97 & 'bravais-lattice index = ',I12,/,5X, & 98 & 'lattice parameter (alat) = ',F12.4,' a.u.',/,5X, & 99 & 'unit-cell volume = ',F12.4,' (a.u.)^3',/,5X, & 100 & 'number of atoms/cell = ',I12,/,5X, & 101 & 'number of atomic types = ',I12) 102 ! 103 WRITE( stdout, '(/2(3X,3(2X,"celldm(",I1,")=",F11.6),/))' ) & 104 ( i, celldm(i), i = 1, 6 ) 105 ! 106 107 !lattice vectors in Angs 108 rr = at*alat*bohr_radius_angs 109 r1(:) = rr(:,1) 110 r2(:) = rr(:,2) 111 r3(:) = rr(:,3) 112 WRITE( stdout, '(5X, & 113 & "Lattice vectors: (cart. coord. in Angs)",/, & 114 & 3(15x,"a(",i1,") = (",3f11.6," ) ",/ ) )') (apol, & 115 (rr (ipol, apol) , ipol = 1, 3) , apol = 1, 3) 116 117 !atoms 118 WRITE( stdout, '(/,5x,"Atomic coordiantes")') 119 WRITE( stdout, '(5x,"site n. atom positions (Angs)")') 120 WRITE( stdout, '(6x,i4,8x,a6," tau(",i4,") = (",3f12.7," )")') & 121 (na, atm(ityp(na)), na, (tau(ipol,na)*alat*bohr_radius_angs, ipol=1,3), na=1,nat) 122 123 !symmetries 124 IF (verbose) THEN 125 write(6,*) 126 127 WRITE( stdout, '(36x,"s",24x,"frac. trans.")') 128 DO isym = 1, nsym 129 WRITE( stdout, '(/6x,"isym = ",i2,5x,a45/)') isym, sname(isym) 130 IF ( ft(1,isym)**2 + ft(2,isym)**2 + ft(3,isym)**2 > 1.0d-8 ) THEN 131 ft1 = at(1,1)*ft(1,isym) + at(1,2)*ft(2,isym) + at(1,3)*ft(3,isym) 132 ft2 = at(2,1)*ft(1,isym) + at(2,2)*ft(2,isym) + at(2,3)*ft(3,isym) 133 ft3 = at(3,1)*ft(1,isym) + at(3,2)*ft(2,isym) + at(3,3)*ft(3,isym) 134 WRITE( stdout, '(1x,"cryst.",3x,"s(",i2,") = (",3(i6,5x), & 135 & " ) f =( ",f10.7," )")') & 136 isym, (s(1,ipol,isym),ipol=1,3), ft(1,isym) 137 WRITE( stdout, '(17x," (",3(i6,5x), " ) ( ",f10.7," )")') & 138 (s(2,ipol,isym),ipol=1,3), ft(2,isym) 139 WRITE( stdout, '(17x," (",3(i6,5x), " ) ( ",f10.7," )"/)') & 140 (s(3,ipol,isym),ipol=1,3), ft(3,isym) 141 WRITE( stdout, '(1x,"cart. ",3x,"s(",i2,") = (",3f11.7, & 142 & " ) f =( ",f10.7," )")') & 143 isym, (sr(1,ipol,isym),ipol=1,3), ft1 144 WRITE( stdout, '(17x," (",3f11.7, " ) ( ",f10.7," )")') & 145 (sr(2,ipol,isym),ipol=1,3), ft2 146 WRITE( stdout, '(17x," (",3f11.7, " ) ( ",f10.7," )"/)') & 147 (sr(3,ipol,isym),ipol=1,3), ft3 148 ELSE 149 WRITE( stdout, '(1x,"cryst.",3x,"s(",i2,") = (",3(i6,5x), " )")') & 150 isym, (s (1, ipol, isym) , ipol = 1,3) 151 WRITE( stdout, '(17x," (",3(i6,5x)," )")') (s(2,ipol,isym), ipol=1,3) 152 WRITE( stdout, '(17x," (",3(i6,5x)," )"/)') (s(3,ipol,isym), ipol=1,3) 153 WRITE( stdout, '(1x,"cart. ",3x,"s(",i2,") = (",3f11.7," )")') & 154 isym, (sr (1, ipol,isym) , ipol = 1, 3) 155 WRITE( stdout, '(17x," (",3f11.7," )")') (sr (2, ipol,isym) , ipol = 1, 3) 156 WRITE( stdout, '(17x," (",3f11.7," )"/)') (sr (3, ipol,isym) , ipol = 1, 3) 157 END IF 158 END DO 159 END IF 160 end if 161 162100 CALL read_line( input_line, end_of_file=tend ) 163 ! 164 IF( tend ) GOTO 120 165 IF( input_line == ' ' .OR. input_line(1:1) == '#' .OR. & 166 input_line(1:1) == '!' ) GOTO 100 167 ! 168 READ (input_line, *) card 169 ! 170 DO i = 1, len_trim( input_line ) 171 input_line( i : i ) = capital( input_line( i : i ) ) 172 ENDDO 173 174 IF ( trim(card) == 'RAMAN_TENSOR') THEN 175 176! read forces from input card "RAMAN_TENSOR" 177! Arrigo Calzolari's convention (to be automated) 178! 179! npol (lpuma) = 1,2,3,4, --> -2h, -h, +h, +2h 180! npol (else) = 1,2,3,4, --> -1-1, +1+1, +1,-1, -1+1 181! npol1 (else) = 1,2 --> -h, h 182! ndiag = 1,2,3 --> Ex, Ey, Ez 183! noffd = 1,2,3 --> Exy, Exz, Eyz 184 185 npol=npol_rm 186 allocate(F0(3,nat)) 187 allocate(dechi(3,3,3,nat)) 188 if (lpuma) then 189 allocate(Fd(npol,ndiag,3,nat)) 190 else 191 allocate(Fd(npol1,ndiag,3,nat)) 192 end if 193 allocate(Fij(npol,noffd,3,nat)) 194 195 F0(:,:)=0.0d0 196 dechi(:,:,:,:)=0.0d0 197 Fd(:,:,:,:)=0.0d0 198 Fij(:,:,:,:)=0.0d0 199 200 ! read data from input 201 202 do i=1,nat 203 read(5,*) (F0(k,i), k=1,3) 204 end do 205 206 do ii=1,ndiag 207 do p=1,npol1 208 do i=1,nat 209 read(5,*) (Fd(p,ii,k,i), k=1,3) 210 end do 211 enddo 212 end do 213 214 do ii=1,noffd 215 do p=1,npol 216 do i=1,nat 217 read(5,*) (Fij(p,ii,k,i), k=1,3) 218 end do 219 enddo 220 end do 221 222 dechi(:,:,:,:)=0.0d0 223 de=de_raman 224 do i=1,nat 225 do ii=1,3 226 do k=1,3 227 if (lpuma) then 228 dechi(ii,ii,k,i)=(-1.0*Fd(1,ii,k,i)+16.0*Fd(2,ii,k,i)-30.0*F0(k,i)+16.0*Fd(3,ii,k,i) & 229 -1.0*Fd(4,ii,k,i))/(12.*de**2) 230 else 231 dechi(ii,ii,k,i)=(Fd(1,ii,k,i)-2*F0(k,i)+Fd(2,ii,k,i))/(de**2) 232 end if 233 end do 234 end do 235 end do 236 237! construct d chi/dE1dE2 238 IF (verbose) THEN 239 write(6,*) '**************************************************' 240 WRITE( stdout, '("unsymmetrized Raman tensor")') 241 write(6,*) '**************************************************' 242 END IF 243 244 do i=1,nat 245 do k=1,3 246 if (lpuma) then 247 dechi(1,2,k,i) = (-1.0*Fij(1,1,k,i)+16.0*Fij(2,1,k,i)-30.0*F0(k,i)+16.0*Fij(3,1,k,i) & 248 -1.0*Fij(4,1,k,i))/(12.0*de**2) 249 dechi(1,2,k,i) = 0.5*dechi(1,2,k,i)-0.5*dechi(1,1,k,i)-0.5*dechi(2,2,k,i) 250 dechi(2,1,k,i) = dechi(1,2,k,i) 251 252 dechi(1,3,k,i) = (-1.0*Fij(1,2,k,i)+16.0*Fij(2,2,k,i)-30.0*F0(k,i)+16.0*Fij(3,2,k,i) & 253 -1.0*Fij(4,2,k,i))/(12.0*de**2) 254 dechi(1,3,k,i) = 0.5*dechi(1,3,k,i)-0.5*dechi(1,1,k,i)-0.5*dechi(3,3,k,i) 255 dechi(3,1,k,i) = dechi(1,3,k,i) 256 257 dechi(2,3,k,i) = (-1.0*Fij(1,3,k,i)+16.0*Fij(2,3,k,i)-30.0*F0(k,i)+ 16.0*Fij(3,3,k,i) & 258 -1.0*Fij(4,3,k,i))/(12.0*de**2) 259 dechi(2,3,k,i) = 0.5*dechi(2,3,k,i)-0.5*dechi(2,2,k,i)-0.5*dechi(3,3,k,i) 260 dechi(3,2,k,i) = dechi(2,3,k,i) 261 else 262 dechi(1,2,k,i) = (Fij(1,1,k,i)+Fij(2,1,k,i)-Fij(3,1,k,i)-Fij(4,1,k,i))/(4*de**2) 263 dechi(2,1,k,i) = dechi(1,2,k,i) 264 265 dechi(1,3,k,i) = (Fij(1,2,k,i)+Fij(2,2,k,i)-Fij(3,2,k,i)-Fij(4,2,k,i))/(4*de**2) 266 dechi(3,1,k,i) = dechi(1,3,k,i) 267 268 dechi(2,3,k,i) = (Fij(1,3,k,i)+Fij(2,3,k,i)-Fij(3,3,k,i)-Fij(4,3,k,i))/(4*de**2) 269 dechi(3,2,k,i) = dechi(2,3,k,i) 270 end if 271 end do 272 end do 273 274 do i=1,nat 275 do ii=1,3 276 do jj=1,3 277 do k=1,3 278 dechi(ii,jj,k,i)=dechi(ii,jj,k,i)/omega ! *(-1.0d0) 279 end do 280 end do 281 end do 282 end do 283 284 IF (verbose) THEN 285 do i=1,nat 286 do k = 1, 3 287 write (6,'(5x,"atom # ",i4," pol.",i3)') i,k 288 do ii =1,3 289 write(6,43) (dechi(ii,jj,k,i)*omega, jj=1,3) 290 end do 291 end do 292 write(6,*)' ' 293 end do 294 END IF 295 write(6,*) '**************************************************' 296 WRITE( stdout, '("Raman tensor (A^2)")') 297 write(6,*) '**************************************************' 298 299 ! convert in crystal coordinates 300 301 do na = 1,nat 302 call cart_to_crys_mat3 ( dechi(1,1,1,na) ) 303 end do 304 305 call symtensor3( nat, dechi) 306 307 do i=1,nat 308 do ii=1,3 309 do jj=1,3 310 do k=1,3 311 dechi(ii,jj,k,i)=conv*omega*dechi(ii,jj,k,i) 312 end do 313 end do 314 end do 315 end do 316 317 do i=1,nat 318 do k = 1, 3 319 write (6,'(5x,"atom # ",i4," pol.",i3)') i,k 320 do ii =1,3 321 write(6,34) (dechi(ii,jj,k,i), jj=1,3) 322 end do 323 end do 324 write(6,*)' ' 325 end do 326 327if (lalpha) then 328 write(6,*) '**************************************************' 329 WRITE( stdout, '("Raman alpha tensor")') 330 write(6,*) '**************************************************' 331 332 nmod=3*nat 333 allocate (u(nmod,3,nat)) 334 allocate (ui(nmod,3,nat)) 335 allocate(alpha(3,3,nmod)) 336 allocate(dechi_u(3,3,3,nat)) 337 338 alpha(:,:,:)=0.0d0 339 u(:,:,:)=0.0d0 340 341! read normalized eigenmodes from matdyn.modes 342 343 open (2,file=TRIM(filemodes),form='formatted') 344 345 read(2,*) 346 read(2,*) 347 read(2,*) 348 read(2,*) 349 350 do n=1,nmod 351 read (2,*) 352 do i=1,nat 353 read(2,'(1x,1x,3 (f10.6,1x,f10.6,3x),1x)') (u(n,k,i),ui(n,k,i),k=1,3) 354 do k=1,3 355 u(n,k,i)=u(n,k,i)/Sqrt(amass(ityp(i))) 356 end do 357 end do 358 end do 359 close(2) 360 361 do n=1,nmod 362 do ii=1,3 363 do jj=1,3 364 do i=1,nat 365 do k=1,3 366 dechi_u(ii,jj,k,i)=dechi(ii,jj,k,i)*u(n,k,i) 367 alpha(ii,jj,n)=alpha(ii,jj,n)+dechi_u(ii,jj,k,i) 368 end do 369 end do 370 alpha(ii,jj,n)=Sqrt(omega)*alpha(ii,jj,n) 371 end do 372 end do 373 end do 374 375 write(6,*)'' 376 do n=1,nmod 377 write(6,*) n 378 do ii=1,3 379 write(6,43) (alpha(ii,jj,n), jj=1,3) 380 end do 381 end do 382 deallocate(u,alpha,dechi_u) 383 end if 384 385 38632 format(a,i5,a,i5) 38743 format(3f12.6) 38834 format(3e24.12) 389deallocate(F0,dechi,Fd,Fij) 390 391 ELSEIF ( trim(card) == 'DIELECTRIC_TENSOR') THEN 392 393 write(6,*) '**************************************************' 394 WRITE( stdout, '("Dielectric tensor")') 395 write(6,*) '**************************************************' 396 397 npol=npol_eps 398 399 allocate (pol0(3)) 400 allocate (pol(npol,3,3)) 401 allocate (eps(3,3)) 402! 403 pol0=0.0d0 404 de=de_eps*omega 405 read(5,*) (pol0(i),i=1,3) 406 do i=1,3 407 if (Abs (pol0(i)) .lt. conv) pol0(i)=0.0d0 408 end do 409! 410 pol=0.0d0 411 do p=1,npol 412 do i=1,3 413 read(5,*) (pol(p,i,j), j=1,3) 414 end do 415 end do 416! 417 do i=1,3 418 do j=1,3 419 if (npol==2) then 420 eps(i,j)=(pol(1,i,j)-pol0(j))-(pol(2,i,j)-pol0(j)) 421 eps(i,j)=2*pi*eps(i,j)/de 422 else 423 eps(i,j)=(pol(1,i,j)-pol0(j)) 424 eps(i,j)=4*pi*eps(i,j)/de 425 end if 426 if (i==j) eps(i,j)=eps(i,j)+1.0 427 end do 428 end do 429! 430 call symmatrix ( eps ) 431 432 do i=1,3 433 write(6,33) (eps(i,j),j=1,3) 434 end do 435 33 format(3f10.4) 436 deallocate(pol0,pol,eps) 437 438 ELSEIF ( trim(card) == 'BORN_CHARGES') THEN 439 write(6,*) '**************************************************' 440 WRITE( stdout, '("Born effective charges")') 441 write(6,*) '**************************************************' 442 443! npol=1 --> code reads forces due to Efield along 444! x,y,z directions. In this exact order. 445! npol=2 --> code reads forces due to Efield along 446! x,y,z,-x,-y,-z (in this order) 447 448 npol=npol_zeu 449 de=de_zeu*sqrt(2.0d0) 450 451 allocate (F0(3,nat)) 452 allocate (Fij(npol,nat,3,3)) 453 allocate (zeta(3,3,nat)) 454 455 F0(:,:)=0.0d0 456 Fij(:,:,:,:)=0.0d0 457 zeta=0.0d0 458 459 do k=1,nat 460 read(5,*) (F0(j,k),j=1,3) 461 end do 462 463 do p=1,npol 464 do i=1,3 ! 3 --> x,y,z 465 do k=1,nat 466 read(5,*) (Fij(p,k,i,j), j=1,3) 467 end do 468 end do 469 end do 470 do k=1,nat 471 do i=1,3 472 do j=1,3 473 if (npol==2) then 474 zeta(i,j,k)=(Fij(1,k,i,j)-F0(j,k))-(Fij(2,k,i,j)-F0(j,k)) 475 zeta(i,j,k)=0.5*zeta(i,j,k)/de 476 else 477 zeta(i,j,k)=(Fij(1,k,i,j)-F0(j,k)) 478 zeta(i,j,k)=zeta(i,j,k)/de 479 end if 480 end do 481 end do 482 end do 483 ! impose ASR 484 do i=1,3 485 do j=1,3 486 sum=0.0d0 487 do k=1,nat 488 sum = sum + zeta(i,j,k) 489 end do 490 do k=1,nat 491 zeta(i,j,k) = zeta(i,j,k) - sum/nat 492 end do 493 end do 494 end do 495 ! symmetrize 496 call symtensor (nat, zeta) 497 498 do k=1,nat 499 write(6,54) 'atom ', k 500 do i=1,3 501 write(6,53) (zeta(i,j,k),j=1,3) 502 end do 503 end do 504 53 format(3f14.8) 505 54 format(a,2x,i5) 506 deallocate(F0,Fij,zeta) 507 508 ELSE 509 ! 510 IF ( ionode ) & 511 WRITE( stdout,'(A)') 'Warning: card '//trim(input_line)//' ignored' 512 ! 513 514 END IF 515 GO TO 100 516 517120 CONTINUE 518 519end program fd_raman 520 521 522 SUBROUTINE cart_to_crys_mat3 ( mat3 ) 523 !----------------------------------------------------------------------- 524 ! 525 USE kinds, ONLY : dp 526 USE cell_base, ONLY : at 527 528 IMPLICIT NONE 529 ! 530 REAL(DP), intent(INOUT) :: mat3(3,3,3) 531 ! 532 REAL(DP) :: work(3,3,3) 533 INTEGER :: i,j,k,l,m,n 534 ! 535 work(:,:,:) = 0.0d0 536 DO i = 1, 3 537 DO j = 1, 3 538 DO k = 1, 3 539 DO l = 1, 3 540 DO m = 1, 3 541 DO n = 1, 3 542 work (i, j, k) = work (i, j, k) + & 543 mat3 (l, m, n) * at (l, i) * at (m, j) * at (n, k) 544 END DO 545 END DO 546 END DO 547 END DO 548 END DO 549 END DO 550 mat3(:,:,:) = work (:,:,:) 551 ! 552 END SUBROUTINE cart_to_crys_mat3 553