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