scnc
/
quantum-espresso
mirror of https://gitlab.com/QEF/q-e.git
1
0
Fork 0
Code Issues Packages Projects Releases Wiki Activity
quantum-espresso/PHonon/PH/matdyn.f90

3024 lines
105 KiB
Fortran
Raw Permalink Blame History

! Copyright (C) 2001-2012 Quantum ESPRESSO group
! This file is distributed under the terms of the
! GNU General Public License. See the file `License'
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
!
Module ifconstants
!
!! All variables read from file that need dynamical allocation.
!
USE kinds, ONLY: DP
!
REAL(DP), ALLOCATABLE :: frc(:,:,:,:,:,:,:)
!! interatomic force constants in real space
REAL(DP), ALLOCATABLE :: frc_lr(:,:,:,:,:,:,:)
!! long-range part of interatomic force constants in real space
REAL(DP), ALLOCATABLE :: tau_blk(:,:)
!! atomic positions for the original cell
REAL(DP), ALLOCATABLE :: zeu(:,:,:)
!! effective charges for the original cell
REAL(DP), ALLOCATABLE :: m_loc(:,:)
!! the magnetic moments of each atom
INTEGER, ALLOCATABLE :: ityp_blk(:)
!! atomic types for each atom of the original cell
!
CHARACTER(LEN=3), ALLOCATABLE :: atm(:)
!
end Module ifconstants
!
!
!---------------------------------------------------------------------
PROGRAM matdyn
!-----------------------------------------------------------------------
!! This program calculates the phonon frequencies for a list of generic
!! q vectors starting from the interatomic force constants generated
!! from the dynamical matrices as written by DFPT phonon code through
!! the companion program \(\texttt{q2r}\).
!
!! \(\texttt{matdyn}\) can generate a supercell of the original cell for
!! mass approximation calculation. If supercell data are not specified
!! in input, the unit cell, lattice vectors, atom types and positions
!! are read from the force constant file.
!
! Input cards: namelist &input
! flfrc file produced by q2r containing force constants (needed)
! It is the same as in the input of q2r.x (+ the .xml extension
! if the dynamical matrices produced by ph.x were in xml
! format). No default value: must be specified.
! asr (character) indicates the type of Acoustic Sum Rule imposed
! - 'no': no Acoustic Sum Rules imposed (default)
! - 'simple': previous implementation of the asr used
! (3 translational asr imposed by correction of
! the diagonal elements of the force constants matrix)
! - 'crystal': 3 translational asr imposed by optimized
! correction of the force constants (projection).
! - 'all': 3 translational asr + 3 rotational asr + 15 Huang
! conditions for vanishing stress tensor, imposed by
! optimized correction of the force constants (projection).
! Remember to set write_lr to .true. to write long-range
! force constants into file when running q2r and set read_lr
! to .true. when running matdyn for the case of polar system.
! - 'one-dim': 3 translational asr + 1 rotational asr
! imposed by optimized correction of the force constants
! (the rotation axis is the direction of periodicity;
! it will work only if this axis considered is one of
! the cartesian axis).
! - 'zero-dim': 3 translational asr + 3 rotational asr
! imposed by optimized correction of the force constants
! Note that in certain cases, not all the rotational asr
! can be applied (e.g. if there are only 2 atoms in a
! molecule or if all the atoms are aligned, etc.).
! In these cases the supplementary asr are cancelled
! during the orthonormalization procedure (see below).
! huang if .true. (default) Huang conditions for vanishing
! stress tensor are included in asr='all'.
! dos if .true. calculate phonon Density of States (DOS)
! using tetrahedra and a uniform q-point grid (see below)
! NB: may not work properly in noncubic materials
! if .false. calculate phonon bands from the list of q-points
! supplied in input (default)
! nk1,nk2,nk3 uniform q-point grid for DOS calculation (includes q=0)
! (must be specified if dos=.true., ignored otherwise)
! deltaE energy step, in cm^(-1), for DOS calculation: from min
! to max phonon energy (default: 1 cm^(-1) if ndos, see
! below, is not specified)
! ndos number of energy steps for DOS calculations
! (default: calculated from deltaE if not specified)
! degauss DOS broadening (in cm^-1). Default 0 - meaning use tetrahedra
! fldos output file for dos (default: 'matdyn.dos')
! the dos is in states/cm(-1) plotted vs omega in cm(-1)
! and is normalised to 3*nat, i.e. the number of phonons
! flfrq output file for frequencies (default: 'matdyn.freq')
! flvec output file for normalized phonon displacements
! (default: 'matdyn.modes'). The normalized phonon displacements
! are the eigenvectors divided by the square root of the mass,
! then normalized. As such they are not orthogonal.
! fleig output file for phonon eigenvectors (default: 'matdyn.eig')
! The phonon eigenvectors are the eigenvectors of the dynamical
! matrix. They are orthogonal.
! fldyn output file for dynamical matrix (default: ' ' i.e. not written)
! at supercell lattice vectors - must form a superlattice of the
! original lattice (default: use original cell)
! l1,l2,l3 supercell lattice vectors are original cell vectors times
! l1, l2, l3 respectively (default: 1, ignored if at specified)
! ntyp number of atom types in the supercell (default: ntyp of the
! original cell)
! amass masses of atoms in the supercell (a.m.u.), one per atom type
! (default: use masses read from file flfrc)
! readtau read atomic positions of the supercell from input
! (used to specify different masses) (default: .false.)
! fltau write atomic positions of the supercell to file "fltau"
! (default: fltau=' ', do not write)
! la2F if .true. interpolates also the el-ph coefficients.
! q_in_band_form if .true. the q points are given in band form:
! Only the first and last point of one or more lines
! are given. See below. (default: .false.).
! q_in_cryst_coord if .true. input q points are in crystalline
! coordinates (default: .false.)
! eigen_similarity: use similarity of the displacements to order
! frequencies (default: .false.)
! NB: You cannot use this option with the symmetry
! analysis of the modes.
! fd (logical) if .t. the ifc come from the finite displacement calculation
! na_ifc (logical) add non analitic contributions to the interatomic force
! constants if finite displacement method is used (as in Wang et al.
! Phys. Rev. B 85, 224303 (2012)) [to be used in conjunction with fd.x]
! nosym if .true., no symmetry and no time reversal are imposed
! loto_2d set to .true. to activate two-dimensional treatment of LO-TO splitting.
! loto_disable (logical) if .true. do not apply LO-TO splitting for q=0
! (default: .false.)
! read_lr set to .true. to read long-range force constants from file,
! when enforcing asr='all' for polar solids in matdyn.
! write_frc set to .true. to write force constants into file.
! The filename would be flfrc+".matdyn".
!
! if (readtau) atom types and positions in the supercell follow:
! (tau(i,na),i=1,3), ityp(na)
! IF (q_in_band_form.and..not.dos) THEN
! nq ! number of q points
! (q(i,n),i=1,3), nptq nptq is the number of points between this point
! and the next. These points are automatically
! generated. the q points are given in Cartesian
! coordinates, 2pi/a units (a=lattice parameters)
! ELSE, if (.not.dos) :
! nq number of q-points
! (q(i,n), i=1,3) nq q-points in cartesian coordinates, 2pi/a units
! If q = 0, the direction qhat (q=>0) for the non-analytic part
! is extracted from the sequence of q-points as follows:
! qhat = q(n) - q(n-1) or qhat = q(n) - q(n+1)
! depending on which one is available and nonzero.
! For low-symmetry crystals, specify twice q = 0 in the list
! if you want to have q = 0 results for two different directions
!
USE kinds, ONLY : DP
USE mp, ONLY : mp_bcast
USE mp_world, ONLY : world_comm
USE mp_global, ONLY : mp_startup, mp_global_end
USE environment, ONLY : environment_start, environment_end
USE io_global, ONLY : ionode, ionode_id, stdout
USE io_dyn_mat, ONLY : read_dyn_mat_param, read_dyn_mat_header, &
read_ifc_param, read_ifc
USE cell_base, ONLY : at, bg, celldm
USE constants, ONLY : RY_TO_THZ, RY_TO_CMM1, amu_ry
USE symm_base, ONLY : set_sym
USE rap_point_group, ONLY : code_group
USE bz_form, ONLY : transform_label_coord
USE parser, ONLY : read_line
USE rigid, ONLY : dyndiag, nonanal, nonanal_ifc
USE el_phon, ONLY : el_ph_nsigma
USE ifconstants, ONLY : frc, frc_lr, atm, zeu, tau_blk, ityp_blk, m_loc
!
IMPLICIT NONE
!
INTEGER :: gid
!
! variables *_blk refer to the original cell, other variables
! to the (super)cell (which may coincide with the original cell)
!
INTEGER:: nax, nax_blk
INTEGER, PARAMETER:: ntypx=10, nrwsx=200
REAL(DP), PARAMETER :: eps=1.0d-6
INTEGER :: nr1, nr2, nr3, nsc, nk1, nk2, nk3, ibrav
CHARACTER(LEN=256) :: flfrc, flfrq, flvec, fltau, fldos, filename, fldyn, &
fleig, fildyn, fildyn_prefix
CHARACTER(LEN=10) :: asr
LOGICAL :: dos, has_zstar, q_in_cryst_coord, eigen_similarity, loto_disable
COMPLEX(DP), ALLOCATABLE :: dyn(:,:,:,:), dyn_blk(:,:,:,:), frc_ifc(:,:,:,:)
COMPLEX(DP), ALLOCATABLE :: z(:,:)
REAL(DP), ALLOCATABLE:: tau(:,:), q(:,:), w2(:,:), freq(:,:), wq(:), &
dynq(:,:,:), DOSofE(:), zq(:, :, :)
INTEGER, ALLOCATABLE:: ityp(:), itau_blk(:)
REAL(DP) :: omega,alat, &! cell parameters and volume
at_blk(3,3), bg_blk(3,3), &! original cell
omega_blk, &! original cell volume
epsil(3,3), &! dielectric tensor
amass(ntypx), &! atomic masses
amass_blk(ntypx), &! original atomic masses
atws(3,3), &! lattice vector for WS initialization
rws(0:3,nrwsx) ! nearest neighbor list, rws(0,*) = norm^2
REAL(DP) :: alph ! Ewald coefficient for non-analytical term
!
INTEGER :: nat, nat_blk, ntyp, ntyp_blk, &
l1, l2, l3, &! supercell dimensions
nrws, &! number of nearest neighbor
code_group_old
INTEGER :: nspin_mag, nqs, ios, ipol
!
LOGICAL :: readtau, la2F, xmlifc, lo_to_split, na_ifc, fd, nosym, loto_2d, &
huang, read_lr, write_frc
!
REAL(DP) :: qhat(3), qh, DeltaE, Emin=0._dp, Emax, E, qq, degauss
REAL(DP) :: delta, pathL
REAL(DP), ALLOCATABLE :: xqaux(:,:)
INTEGER, ALLOCATABLE :: nqb(:)
INTEGER :: n, i, j, it, nq, nqx, na, nb, ndos, iout, nqtot, iout_dyn, iout_eig
LOGICAL, EXTERNAL :: has_xml
INTEGER, ALLOCATABLE :: num_rap_mode(:,:)
LOGICAL, ALLOCATABLE :: high_sym(:)
LOGICAL :: q_in_band_form
! .... variables for band plotting based on similarity of eigenvalues
COMPLEX(DP), ALLOCATABLE :: tmp_z(:,:)
REAL(DP), ALLOCATABLE :: abs_similarity(:,:), tmp_w2(:)
COMPLEX(DP), ALLOCATABLE :: f_of_q(:,:,:,:)
INTEGER :: location(1), isig
CHARACTER(LEN=6) :: int_to_char
LOGICAL, ALLOCATABLE :: mask(:)
INTEGER :: npk_label, nch
CHARACTER(LEN=3), ALLOCATABLE :: letter(:)
INTEGER, ALLOCATABLE :: label_list(:)
LOGICAL :: tend, terr
CHARACTER(LEN=256) :: input_line, buffer
CHARACTER(LEN=10) :: point_label_type
CHARACTER(len=80) :: k_points = 'tpiba'
!
REAL(DP), external :: dos_broad, dos_gam
!
NAMELIST /input/ flfrc, amass, asr, flfrq, flvec, fleig, at, dos, &
& fldos, nk1, nk2, nk3, l1, l2, l3, ntyp, readtau, fltau, &
& la2F, ndos, DeltaE, degauss, q_in_band_form, q_in_cryst_coord, &
& eigen_similarity, fldyn, na_ifc, fd, point_label_type, &
& nosym, loto_2d, fildyn, fildyn_prefix, el_ph_nsigma, &
& loto_disable, huang, read_lr, write_frc
!
CALL mp_startup()
CALL environment_start('MATDYN')
!
IF (ionode) CALL input_from_file ( )
!
! ... all calculations are done by the first cpu
!
! set namelist default
!
dos = .FALSE.
deltaE = 1.0d0
degauss = 0
ndos = 1
nk1 = 0
nk2 = 0
nk3 = 0
asr ='no'
readtau=.FALSE.
flfrc=' '
fldos='matdyn.dos'
flfrq='matdyn.freq'
flvec='matdyn.modes'
fleig=' '
fldyn=' '
fltau=' '
fildyn = ' '
fildyn_prefix = ' '
amass(:) =0.d0
amass_blk(:) =0.d0
at(:,:) = 0.d0
ntyp = 0
l1=1
l2=1
l3=1
la2F=.false.
q_in_band_form=.FALSE.
eigen_similarity=.FALSE.
q_in_cryst_coord = .FALSE.
na_ifc=.FALSE.
fd=.FALSE.
point_label_type='SC'
nosym = .false.
loto_2d=.false.
el_ph_nsigma=10
loto_disable = .false.
huang = .true.
read_lr = .false.
write_frc = .false.
!
!
IF (ionode) READ (5,input,IOSTAT=ios)
CALL mp_bcast(ios, ionode_id, world_comm)
CALL errore('matdyn', 'reading input namelist', ABS(ios))
CALL mp_bcast(dos,ionode_id, world_comm)
CALL mp_bcast(deltae,ionode_id, world_comm)
CALL mp_bcast(ndos,ionode_id, world_comm)
CALL mp_bcast(degauss,ionode_id, world_comm)
CALL mp_bcast(nk1,ionode_id, world_comm)
CALL mp_bcast(nk2,ionode_id, world_comm)
CALL mp_bcast(nk3,ionode_id, world_comm)
CALL mp_bcast(asr,ionode_id, world_comm)
CALL mp_bcast(readtau,ionode_id, world_comm)
CALL mp_bcast(flfrc,ionode_id, world_comm)
CALL mp_bcast(fldos,ionode_id, world_comm)
CALL mp_bcast(flfrq,ionode_id, world_comm)
CALL mp_bcast(flvec,ionode_id, world_comm)
CALL mp_bcast(fleig,ionode_id, world_comm)
CALL mp_bcast(fldyn,ionode_id, world_comm)
CALL mp_bcast(fltau,ionode_id, world_comm)
CALL mp_bcast(fildyn,ionode_id, world_comm)
CALL mp_bcast(fildyn_prefix,ionode_id, world_comm)
CALL mp_bcast(amass,ionode_id, world_comm)
CALL mp_bcast(amass_blk,ionode_id, world_comm)
CALL mp_bcast(at,ionode_id, world_comm)
CALL mp_bcast(ntyp,ionode_id, world_comm)
CALL mp_bcast(l1,ionode_id, world_comm)
CALL mp_bcast(l2,ionode_id, world_comm)
CALL mp_bcast(l3,ionode_id, world_comm)
CALL mp_bcast(na_ifc,ionode_id, world_comm)
CALL mp_bcast(fd,ionode_id, world_comm)
CALL mp_bcast(la2F,ionode_id, world_comm)
CALL mp_bcast(q_in_band_form,ionode_id, world_comm)
CALL mp_bcast(eigen_similarity,ionode_id, world_comm)
CALL mp_bcast(q_in_cryst_coord,ionode_id, world_comm)
CALL mp_bcast(point_label_type,ionode_id, world_comm)
CALL mp_bcast(loto_2d,ionode_id, world_comm)
CALL mp_bcast(loto_disable,ionode_id, world_comm)
CALL mp_bcast(el_ph_nsigma,ionode_id, world_comm)
CALL mp_bcast(huang,ionode_id, world_comm)
CALL mp_bcast(read_lr,ionode_id, world_comm)
CALL mp_bcast(write_frc,ionode_id, world_comm)
!
IF (loto_2d .AND. loto_disable) CALL errore('matdyn', &
'loto_2d and loto_disable cannot be both true', 1)
!
! read force constants
!
IF ( trim( fildyn ) /= ' ' ) THEN
IF (ionode) THEN
WRITE(stdout, *)
WRITE(stdout, '(4x,a)') ' fildyn has been provided, running q2r...'
END IF
CALL do_q2r(fildyn, flfrc, fildyn_prefix, asr, la2F, loto_2d, read_lr)
END IF
!
ntyp_blk = ntypx ! avoids fake out-of-bound error
xmlifc=has_xml(flfrc)
IF (xmlifc) THEN
CALL read_dyn_mat_param(flfrc,ntyp_blk,nat_blk)
ALLOCATE (m_loc(3,nat_blk))
ALLOCATE (tau_blk(3,nat_blk))
ALLOCATE (ityp_blk(nat_blk))
ALLOCATE (atm(ntyp_blk))
ALLOCATE (zeu(3,3,nat_blk))
CALL read_dyn_mat_header(ntyp_blk, nat_blk, ibrav, nspin_mag, &
celldm, at_blk, bg_blk, omega_blk, atm, amass_blk, &
tau_blk, ityp_blk, m_loc, nqs, has_zstar, epsil, zeu )
alat=celldm(1)
call volume(alat,at_blk(1,1),at_blk(1,2),at_blk(1,3),omega_blk)
CALL read_ifc_param(nr1,nr2,nr3)
ALLOCATE(frc(nr1,nr2,nr3,3,3,nat_blk,nat_blk))
ALLOCATE(frc_lr(nr1,nr2,nr3,3,3,nat_blk,nat_blk))
frc_lr = 0.d0
if (read_lr) THEN
CALL read_ifc(alph,nr1,nr2,nr3,nat_blk,frc,frc_lr)
else
CALL read_ifc(alph,nr1,nr2,nr3,nat_blk,frc)
end if
ELSE
CALL readfc ( flfrc, nr1, nr2, nr3, epsil, nat_blk, &
ibrav, alat, at_blk, ntyp_blk, &
amass_blk, omega_blk, has_zstar, alph, read_lr)
ENDIF
!
CALL recips ( at_blk(1,1),at_blk(1,2),at_blk(1,3), &
bg_blk(1,1),bg_blk(1,2),bg_blk(1,3) )
!
! set up (super)cell
!
if (ntyp < 0) then
call errore ('matdyn','wrong ntyp ', abs(ntyp))
else if (ntyp == 0) then
ntyp=ntyp_blk
end if
!
! masses (for mass approximation)
!
DO it=1,ntyp
IF (amass(it) < 0.d0) THEN
CALL errore ('matdyn','wrong mass in the namelist',it)
ELSE IF (amass(it) == 0.d0) THEN
IF (it.LE.ntyp_blk) THEN
WRITE (stdout,'(a,i3,a,a)') ' mass for atomic type ',it, &
& ' not given; uses mass from file ',flfrc
amass(it) = amass_blk(it)
ELSE
CALL errore ('matdyn','missing mass in the namelist',it)
END IF
END IF
END DO
!
! lattice vectors
!
IF (SUM(ABS(at(:,:))) == 0.d0) THEN
IF (l1.LE.0 .OR. l2.LE.0 .OR. l3.LE.0) CALL &
& errore ('matdyn',' wrong l1,l2 or l3',1)
at(:,1) = at_blk(:,1)*DBLE(l1)
at(:,2) = at_blk(:,2)*DBLE(l2)
at(:,3) = at_blk(:,3)*DBLE(l3)
END IF
!
CALL check_at(at,bg_blk,alat,omega)
!
! the supercell contains "nsc" times the original unit cell
!
nsc = NINT(omega/omega_blk)
IF (ABS(omega/omega_blk-nsc) > eps) &
CALL errore ('matdyn', 'volume ratio not integer', 1)
!
! read/generate atomic positions of the (super)cell
!
nat = nat_blk * nsc
nax = nat
nax_blk = nat_blk
!
ALLOCATE ( tau (3, nat), ityp(nat), itau_blk(nat) )
!
IF (readtau) THEN
CALL read_tau &
(nat, nat_blk, ntyp, bg_blk, tau, tau_blk, ityp, itau_blk)
ELSE
CALL set_tau &
(nat, nat_blk, at, at_blk, tau, tau_blk, ityp, ityp_blk, itau_blk)
ENDIF
!
IF (fltau.NE.' ') CALL write_tau (fltau, nat, tau, ityp)
!
! reciprocal lattice vectors
!
CALL recips (at(1,1),at(1,2),at(1,3),bg(1,1),bg(1,2),bg(1,3))
!
! build the WS cell corresponding to the force constant grid
!
atws(:,1) = at_blk(:,1)*DBLE(nr1)
atws(:,2) = at_blk(:,2)*DBLE(nr2)
atws(:,3) = at_blk(:,3)*DBLE(nr3)
! initialize WS r-vectors
CALL wsinit(rws,nrwsx,nrws,atws)
!
! end of (super)cell setup
!
IF (dos) THEN
IF (nk1 < 1 .OR. nk2 < 1 .OR. nk3 < 1) &
CALL errore ('matdyn','specify correct q-point grid!',1)
nqx = nk1*nk2*nk3
ALLOCATE ( q(3,nqx), wq(nqx) )
CALL gen_qpoints (ibrav, at, bg, nat, tau, ityp, nk1, nk2, nk3, &
nqx, nq, q, nosym, wq)
ELSE
!
! read q-point list
!
IF (ionode) READ (5,*) nq
CALL mp_bcast(nq, ionode_id, world_comm)
ALLOCATE ( q(3,nq) )
IF (.NOT.q_in_band_form) THEN
ALLOCATE(wq(nq))
DO n = 1,nq
IF (ionode) READ (5,*) (q(i,n),i=1,3)
END DO
CALL mp_bcast(q, ionode_id, world_comm)
!
IF (q_in_cryst_coord) CALL cryst_to_cart(nq,q,bg,+1)
ELSE
ALLOCATE( nqb(nq) )
ALLOCATE( xqaux(3,nq) )
ALLOCATE( letter(nq) )
ALLOCATE( label_list(nq) )
npk_label=0
DO n = 1, nq
CALL read_line( input_line, end_of_file = tend, error = terr )
IF (tend) CALL errore('matdyn','Missing lines',1)
IF (terr) CALL errore('matdyn','Error reading q points',1)
DO j=1,256 ! loop over all characters of input_line
IF ( (ICHAR(input_line(j:j)) < 58 .AND. & ! a digit
ICHAR(input_line(j:j)) > 47) &
.OR.ICHAR(input_line(j:j)) == 43 .OR. & ! the + sign
ICHAR(input_line(j:j)) == 45 .OR. & ! the - sign
ICHAR(input_line(j:j)) == 46 ) THEN ! a dot .
!
! This is a digit, therefore this line contains the coordinates of the
! k point. We read it and exit from the loop on characters
!
READ(input_line,*) xqaux(1,n), xqaux(2,n), xqaux(3,n), &
nqb(n)
EXIT
ELSEIF ((ICHAR(input_line(j:j)) < 123 .AND. &
ICHAR(input_line(j:j)) > 64)) THEN
!
! This is a letter, not a space character. We read the next three
! characters and save them in the letter array, save also which k point
! it is
!
npk_label=npk_label+1
READ(input_line(j:),'(a3)') letter(npk_label)
label_list(npk_label)=n
!
! now we remove the letters from input_line and read the number of points
! of the line. The next two line should account for the case in which
! there is only one space between the letter and the number of points.
!
nch=3
IF ( ICHAR(input_line(j+1:j+1))==32 .OR. &
ICHAR(input_line(j+2:j+2))==32 ) nch=2
buffer=input_line(j+nch:)
READ(buffer,*,err=20,iostat=ios) nqb(n)
20 IF (ios /=0) CALL errore('matdyn',&
'problem reading number of points',1)
EXIT
ENDIF
ENDDO
ENDDO
IF (q_in_cryst_coord) k_points='crystal'
IF ( npk_label > 0 ) &
CALL transform_label_coord(ibrav, celldm, xqaux, letter, &
label_list, npk_label, nq, k_points, point_label_type )
DEALLOCATE(letter)
DEALLOCATE(label_list)
CALL mp_bcast(xqaux, ionode_id, world_comm)
CALL mp_bcast(nqb, ionode_id, world_comm)
IF (q_in_cryst_coord) CALL cryst_to_cart(nq,xqaux,bg,+1)
nqtot=SUM(nqb(1:nq-1))+1
DO i=1,nq-1
IF (nqb(i)==0) nqtot=nqtot+1
ENDDO
DEALLOCATE(q)
ALLOCATE(q(3,nqtot))
ALLOCATE(wq(nqtot))
CALL generate_k_along_lines(nq, xqaux, nqb, q, wq, nqtot)
nq=nqtot
DEALLOCATE(xqaux)
DEALLOCATE(nqb)
END IF
!
END IF
!
IF (asr /= 'no') THEN
CALL set_asr (asr, nr1, nr2, nr3, frc, frc_lr, zeu, &
nat_blk, ibrav, tau_blk, at_blk, huang)
END IF
!
IF (write_frc) THEN
CALL writefc(flfrc, xmlifc, has_zstar, nr1, nr2, nr3, ibrav, alat, at_blk, bg_blk, &
ntyp_blk, nat_blk, amass_blk, omega_blk, epsil, alph, nspin_mag, nqs)
END IF
!
IF (flvec.EQ.' ') THEN
iout=0
ELSE
iout=4
IF (ionode) OPEN (unit=iout,file=flvec,status='unknown',form='formatted')
END IF
IF (fldyn.EQ.' ') THEN
iout_dyn=0
ELSE
iout_dyn=44
OPEN (unit=iout_dyn,file=fldyn,status='unknown',form='formatted')
END IF
IF (fleig.EQ.' ') THEN
iout_eig=0
ELSE
iout_eig=313
IF (ionode) OPEN (unit=iout_eig,file=fleig,status='unknown',form='formatted')
END IF
ALLOCATE ( dyn(3,3,nat,nat), dyn_blk(3,3,nat_blk,nat_blk) )
ALLOCATE ( z(3*nat,3*nat), w2(3*nat,nq), f_of_q(3,3,nat,nat), &
dynq(3*nat,nq,nat), DOSofE(nat), zq(nq, 3*nat, 3*nat) )
ALLOCATE ( tmp_w2(3*nat), abs_similarity(3*nat,3*nat), mask(3*nat) )
if(la2F.and.ionode) open(unit=300,file='dyna2F',status='unknown')
IF (xmlifc) CALL set_sym(nat, tau, ityp, nspin_mag, m_loc )
ALLOCATE(num_rap_mode(3*nat,nq))
ALLOCATE(high_sym(nq))
num_rap_mode=-1
high_sym=.TRUE.
zq(:, :, :) = (0.d0, 0.d0)
DO n=1, nq
dyn(:,:,:,:) = (0.d0, 0.d0)
lo_to_split=.FALSE.
f_of_q(:,:,:,:) = (0.d0,0.d0)
IF(na_ifc) THEN
qq=sqrt(q(1,n)**2+q(2,n)**2+q(3,n)**2)
if(abs(qq) < 1d-8) qq=1.0
qhat(1)=q(1,n)/qq
qhat(2)=q(2,n)/qq
qhat(3)=q(3,n)/qq
CALL nonanal_ifc (nat,nat_blk,itau_blk,epsil,qhat,zeu,omega,dyn, &
nr1, nr2, nr3,f_of_q)
END IF
CALL setupmat (q(1,n), dyn, nat, at, bg, tau, itau_blk, nsc, alat, &
dyn_blk, nat_blk, at_blk, bg_blk, tau_blk, omega_blk, &
loto_2d, epsil, zeu, alph, &
frc, nr1,nr2,nr3, has_zstar, rws, nrws, na_ifc,f_of_q,fd)
IF (.not.loto_2d) THEN
qhat(1) = q(1,n)*at(1,1)+q(2,n)*at(2,1)+q(3,n)*at(3,1)
qhat(2) = q(1,n)*at(1,2)+q(2,n)*at(2,2)+q(3,n)*at(3,2)
qhat(3) = q(1,n)*at(1,3)+q(2,n)*at(2,3)+q(3,n)*at(3,3)
IF ( ABS( qhat(1) - NINT (qhat(1) ) ) <= eps .AND. &
ABS( qhat(2) - NINT (qhat(2) ) ) <= eps .AND. &
ABS( qhat(3) - NINT (qhat(3) ) ) <= eps ) THEN
!
! q = 0 : we need the direction q => 0 for the non-analytic part
!
IF ( n == 1 ) THEN
! if q is the first point in the list
IF ( nq > 1 ) THEN
! one more point
qhat(:) = q(:,n) - q(:,n+1)
ELSE
! no more points
qhat(:) = 0.d0
END IF
ELSE IF ( n > 1 ) THEN
! if q is not the first point in the list
IF ( q(1,n-1)==0.d0 .AND. &
q(2,n-1)==0.d0 .AND. &
q(3,n-1)==0.d0 .AND. n < nq ) THEN
! if the preceding q is also 0 :
qhat(:) = q(:,n) - q(:,n+1)
ELSE
! if the preceding q is npt 0 :
qhat(:) = q(:,n) - q(:,n-1)
END IF
END IF
qh = SQRT(qhat(1)**2+qhat(2)**2+qhat(3)**2)
! write(*,*) ' qh, has_zstar ',qh, has_zstar
IF (qh /= 0.d0) qhat(:) = qhat(:) / qh
IF (qh /= 0.d0 .AND. .NOT. has_zstar) THEN
CALL infomsg &
('matdyn','Z* not found in file '//TRIM(flfrc)// &
', TO-LO splitting at q=0 will be absent!')
ELSEIF (loto_disable) THEN
CALL infomsg('matdyn', &
'loto_disable is true. Disable LO-TO splitting at q=0.')
ELSE
lo_to_split=.TRUE.
ENDIF
!
IF (lo_to_split) CALL nonanal (nat, nat_blk, itau_blk, epsil, qhat, zeu, omega, dyn)
!
END IF
!
END IF
if(iout_dyn.ne.0) THEN
call write_dyn_on_file(q(1,n),dyn,nat, iout_dyn)
if(sum(abs(q(:,n)))==0._dp) call write_epsilon_and_zeu (zeu, epsil, nat, iout_dyn)
endif
CALL dyndiag(nat,ntyp,amass,ityp,dyn,w2(1,n),z)
!
! Convert from displacements to eigenvectors (see rigid.f90 :: dyndiag)
do i = 1, 3*nat
do na = 1, nat
do ipol = 1,3
zq(n, (na-1)*3+ipol,i) = z((na-1)*3+ipol,i) * sqrt(amu_ry * amass(ityp(na)))
end do
end do
end do
!
! Atom projection of dynamical matrix
DO i = 1, 3*nat
DO na = 1, nat
dynq(i, n, na) = DOT_PRODUCT(z(3*(na-1)+1:3*na, i), &
& z(3*(na-1)+1:3*na, i) ) &
& * amu_ry * amass(ityp(na))
END DO
END DO
IF (ionode.and.iout_eig.ne.0) &
& CALL write_eigenvectors(nat,ntyp,amass,ityp,q(1,n),w2(1,n),z,iout_eig)
!
! Cannot use the small group of \Gamma to analize the symmetry
! of the mode if there is an electric field.
!
IF (xmlifc.AND..NOT.lo_to_split) THEN
WRITE(stdout,'(10x,"xq=",3F8.4)') q(:,n)
CALL find_representations_mode_q(nat,ntyp,q(:,n), &
w2(:,n),z,tau,ityp,amass, num_rap_mode(:,n), nspin_mag)
IF (code_group==code_group_old.OR.high_sym(n-1)) high_sym(n)=.FALSE.
code_group_old=code_group
ENDIF
IF (eigen_similarity) THEN
! ... order phonon dispersions using similarity of eigenvalues
! ... Courtesy of Takeshi Nishimatsu, IMR, Tohoku University
IF (.NOT.ALLOCATED(tmp_z)) THEN
ALLOCATE(tmp_z(3*nat,3*nat))
ELSE
abs_similarity = ABS ( MATMUL ( CONJG( TRANSPOSE(z)), tmp_z ) )
mask(:) = .true.
DO na=1,3*nat
location = maxloc( abs_similarity(:,na), mask(:) )
mask(location(1)) = .false.
tmp_w2(na) = w2(location(1),n)
tmp_z(:,na) = z(:,location(1))
END DO
w2(:,n) = tmp_w2(:)
z(:,:) = tmp_z(:,:)
END IF
tmp_z(:,:) = z(:,:)
ENDIF
!
if(la2F.and.ionode) then
write(300,*) n
do na=1,3*nat
write(300,*) (z(na,nb),nb=1,3*nat)
end do ! na
endif
!
IF (ionode.and.iout.ne.0) CALL writemodes(nat,q(1,n),w2(1,n),z,iout)
!
END DO !nq
DEALLOCATE (tmp_w2, abs_similarity, mask)
IF (eigen_similarity) DEALLOCATE(tmp_z)
if(la2F.and.ionode) close(300)
!
IF(iout .NE. 0.and.ionode) CLOSE(unit=iout)
IF(iout_dyn .NE. 0) CLOSE(unit=iout_dyn)
IF(iout_eig .NE. 0) CLOSE(unit=iout_eig)
!
ALLOCATE (freq(3*nat, nq))
DO n=1,nq
! freq(i,n) = frequencies in cm^(-1), with negative sign if omega^2 is negative
DO i=1,3*nat
freq(i,n)= SQRT(ABS(w2(i,n))) * RY_TO_CMM1
IF (w2(i,n) < 0.0d0) freq(i,n) = -freq(i,n)
END DO
END DO
!
IF(flfrq.NE.' '.and.ionode) THEN
OPEN (unit=2,file=flfrq ,status='unknown',form='formatted')
WRITE(2, '(" &plot nbnd=",i4,", nks=",i4," /")') 3*nat, nq
DO n=1, nq
WRITE(2, '(10x,3f10.6)') q(1,n), q(2,n), q(3,n)
WRITE(2,'(6f10.4)') (freq(i,n), i=1,3*nat)
END DO
CLOSE(unit=2)
OPEN (unit=2,file=trim(flfrq)//'.gp' ,status='unknown',form='formatted')
pathL = 0._dp
WRITE(2, '(f10.6,3x,999f10.4)') pathL, (freq(i,1), i=1,3*nat)
DO n=2, nq
pathL=pathL+(SQRT(SUM( (q(:,n)-q(:,n-1))**2 )))
WRITE(2, '(f10.6,3x,999f10.4)') pathL, (freq(i,n), i=1,3*nat)
END DO
CLOSE(unit=2)
END IF
!
! If the force constants are in the xml format we write also
! the file with the representations of each mode
!
IF (flfrq.NE.' '.AND.xmlifc.AND.ionode) THEN
filename=TRIM(flfrq)//'.rap'
OPEN (unit=2,file=filename ,status='unknown',form='formatted')
WRITE(2, '(" &plot_rap nbnd_rap=",i4,", nks_rap=",i4," /")') 3*nat, nq
DO n=1, nq
WRITE(2,'(10x,3f10.6,l6)') q(1,n), q(2,n), q(3,n), high_sym(n)
WRITE(2,'(6i10)') (num_rap_mode(i,n), i=1,3*nat)
END DO
CLOSE(unit=2)
END IF
!
IF (dos) THEN
Emin = 0.0d0
Emax = 0.0d0
DO n=1,nq
DO i=1, 3*nat
Emin = MIN (Emin, freq(i,n))
Emax = MAX (Emax, freq(i,n))
END DO
END DO
!
if (ndos > 1) then
DeltaE = (Emax - Emin)/(ndos-1)
else
ndos = NINT ( (Emax - Emin) / DeltaE + 1.51d0 )
end if
IF (ionode) OPEN (unit=2,file=fldos,status='unknown',form='formatted')
IF (ionode) WRITE (2, *) "# Frequency[cm^-1] DOS PDOS"
!
IF (degauss .EQ. 0) THEN
! Use tetrahedra
DO na = 1, nat
dynq(1:3*nat,1:nq,na) = dynq(1:3*nat,1:nq,na) * freq(1:3*nat,1:nq)
END DO
END IF
!
DO n= 1, ndos
E = Emin + (n - 1) * DeltaE
DO na = 1, nat
DOSofE(na) = 0d0
!
IF (degauss .EQ. 0) THEN
! Use tetrahedra
DO i = 1, 3*nat
DOSofE(na) = DOSofE(na) &
& + dos_gam(3*nat, nq, i, dynq(1:3*nat,1:nq,na), freq, E)
END DO
ELSE
! Use broadening
DOSofE(na) = DOSofE(na) &
& + dos_broad(na, 3*nat, nq, freq, zq, wq, E, degauss)
END IF
!
END DO
!
IF (ionode) WRITE (2, '(2ES18.10,1000ES12.4)') E, SUM(DOSofE(1:nat)), DOSofE(1:nat)
END DO
IF (ionode) CLOSE(unit=2)
END IF !dos
DEALLOCATE (z, zq, w2, dyn, dyn_blk)
!
! for a2F
!
IF(la2F) THEN
!
IF (.NOT. dos) THEN
DO isig=1,el_ph_nsigma
OPEN (unit=200+isig,file='elph.gamma.'//&
TRIM(int_to_char(isig)), status='unknown',form='formatted')
WRITE(200+isig, '(" &plot nbnd=",i4,", nks=",i4," /")') 3*nat, nq
END DO
END IF
!
! convert frequencies to Ry
!
freq(:,:)= freq(:,:) / RY_TO_CMM1
Emin = Emin / RY_TO_CMM1
DeltaE=DeltaE/ RY_TO_CMM1
!
call a2Fdos (nat, nq, nr1, nr2, nr3, ibrav, at, bg, tau, alat, &
nsc, nat_blk, at_blk, bg_blk, itau_blk, omega_blk, &
rws, nrws, dos, Emin, DeltaE, ndos, &
asr, q, freq,fd, wq, huang)
!
IF (.NOT.dos) THEN
DO isig=1,el_ph_nsigma
CLOSE(UNIT=200+isig)
ENDDO
ENDIF
END IF
DEALLOCATE ( freq)
DEALLOCATE(num_rap_mode)
DEALLOCATE(high_sym)
DEALLOCATE(frc_lr)
!
CALL environment_end('MATDYN')
!
CALL mp_global_end()
!
IF (asr == 'all') THEN
WRITE(stdout, '(/A117/)') " You are using asr="//'"all"'//", please consider citing "// &
"C. Lin, S. Ponc\'e and N. Marzari, npj Comput Mater 8, 236 (2022)."
ENDIF
!
STOP
!
END PROGRAM matdyn
!
!-----------------------------------------------------------------------
SUBROUTINE readfc ( flfrc, nr1, nr2, nr3, epsil, nat, &
ibrav, alat, at, ntyp, amass, omega, &
has_zstar, alph, read_lr )
!-----------------------------------------------------------------------
!
USE kinds, ONLY : DP
USE ifconstants,ONLY : tau => tau_blk, ityp => ityp_blk, frc, frc_lr, zeu, atm
USE cell_base, ONLY : celldm
USE io_global, ONLY : ionode, ionode_id, stdout
USE mp, ONLY : mp_bcast
USE mp_world, ONLY : world_comm
USE constants, ONLY : amu_ry
!
IMPLICIT NONE
! I/O variable
CHARACTER(LEN=256) :: flfrc
CHARACTER(LEN=80) :: line
INTEGER :: ibrav, nr1,nr2,nr3,nat, ntyp
REAL(DP) :: alat, at(3,3), epsil(3,3), alph
LOGICAL :: has_zstar, read_lr
! local variables
INTEGER :: i, j, na, nb, m1,m2,m3
INTEGER :: ios, ibid, jbid, nabid, nbbid, m1bid,m2bid,m3bid
REAL(DP) :: amass(ntyp), amass_from_file, omega
INTEGER :: nt
!
!
IF (ionode) OPEN (unit=1,file=flfrc,status='old',form='formatted')
!
! read cell data
!
IF (ionode)THEN
READ(1,*) ntyp,nat,ibrav,(celldm(i),i=1,6)
if (ibrav==0) then
read(1,*) ((at(i,j),i=1,3),j=1,3)
end if
ENDIF
CALL mp_bcast(ntyp, ionode_id, world_comm)
CALL mp_bcast(nat, ionode_id, world_comm)
CALL mp_bcast(ibrav, ionode_id, world_comm)
CALL mp_bcast(celldm, ionode_id, world_comm)
IF (ibrav==0) THEN
CALL mp_bcast(at, ionode_id, world_comm)
ENDIF
!
CALL latgen(ibrav,celldm,at(1,1),at(1,2),at(1,3),omega)
alat = celldm(1)
at = at / alat ! bring at in units of alat
CALL volume(alat,at(1,1),at(1,2),at(1,3),omega)
!
! read atomic types, positions and masses
!
ALLOCATE (atm(ntyp))
DO nt = 1,ntyp
IF (ionode) READ(1,*) i,atm(nt),amass_from_file
CALL mp_bcast(i,ionode_id, world_comm)
CALL mp_bcast(amass_from_file,ionode_id, world_comm)
IF (i.NE.nt) CALL errore ('readfc','wrong data read',nt)
IF (amass(nt).EQ.0.d0) THEN
amass(nt) = amass_from_file/amu_ry
ELSE
WRITE(stdout,*) 'for atomic type',nt,' mass from file not used'
END IF
END DO
CALL mp_bcast(atm,ionode_id, world_comm)
CALL mp_bcast(amass,ionode_id, world_comm)
!
ALLOCATE (tau(3,nat), ityp(nat), zeu(3,3,nat))
!
DO na=1,nat
IF (ionode) READ(1,*) i,ityp(na),(tau(j,na),j=1,3)
CALL mp_bcast(i,ionode_id, world_comm)
IF (i.NE.na) CALL errore ('readfc','wrong data read',na)
END DO
CALL mp_bcast(ityp,ionode_id, world_comm)
CALL mp_bcast(tau,ionode_id, world_comm)
!
! read macroscopic variables
!
alph = 1.0_dp
IF (ionode) THEN
READ (1,'(a)') line
READ(line,*,iostat=ios) has_zstar, alph
IF ( ios /= 0 ) READ(line,*) has_zstar
ENDIF
!
CALL mp_bcast(has_zstar,ionode_id, world_comm)
IF (has_zstar) THEN
CALL mp_bcast(alph,ionode_id, world_comm)
IF (ionode) READ(1,*) ((epsil(i,j),j=1,3),i=1,3)
CALL mp_bcast(epsil,ionode_id, world_comm)
IF (ionode) THEN
DO na=1,nat
READ(1,*)
READ(1,*) ((zeu(i,j,na),j=1,3),i=1,3)
END DO
ENDIF
CALL mp_bcast(zeu,ionode_id, world_comm)
ELSE
zeu (:,:,:) = 0.d0
epsil(:,:) = 0.d0
END IF
!
IF (ionode) READ (1,*) nr1,nr2,nr3
CALL mp_bcast(nr1,ionode_id, world_comm)
CALL mp_bcast(nr2,ionode_id, world_comm)
CALL mp_bcast(nr3,ionode_id, world_comm)
!
! read real-space interatomic force constants
!
ALLOCATE ( frc(nr1,nr2,nr3,3,3,nat,nat) )
frc(:,:,:,:,:,:,:) = 0.d0
ALLOCATE ( frc_lr(nr1,nr2,nr3,3,3,nat,nat) )
frc_lr(:,:,:,:,:,:,:) = 0.d0
DO i=1,3
DO j=1,3
DO na=1,nat
DO nb=1,nat
IF (ionode) READ (1,*) ibid, jbid, nabid, nbbid
CALL mp_bcast(ibid,ionode_id, world_comm)
CALL mp_bcast(jbid,ionode_id, world_comm)
CALL mp_bcast(nabid,ionode_id, world_comm)
CALL mp_bcast(nbbid,ionode_id, world_comm)
IF(i .NE.ibid .OR. j .NE.jbid .OR. &
na.NE.nabid .OR. nb.NE.nbbid) &
CALL errore ('readfc','error in reading',1)
IF (read_lr) THEN
IF (ionode) READ (1,*) (((m1bid, m2bid, m3bid, &
frc(m1,m2,m3,i,j,na,nb), &
frc_lr(m1,m2,m3,i,j,na,nb), &
m1=1,nr1),m2=1,nr2),m3=1,nr3)
ELSE
IF (ionode) READ (1,*) (((m1bid, m2bid, m3bid, &
frc(m1,m2,m3,i,j,na,nb), &
m1=1,nr1),m2=1,nr2),m3=1,nr3)
END IF
CALL mp_bcast(frc(:,:,:,i,j,na,nb),ionode_id, world_comm)
CALL mp_bcast(frc_lr(:,:,:,i,j,na,nb),ionode_id, world_comm)
END DO
END DO
END DO
END DO
!
IF (ionode) CLOSE(unit=1)
!
RETURN
END SUBROUTINE readfc
!
!-----------------------------------------------------------------------
SUBROUTINE writefc(flfrc, xmlifc, has_zstar, nr1, nr2, nr3, ibrav, alat, &
at, bg, ntyp, nat, amass, omega, epsil, alph, nspin_mag, nqs)
!---------------------------------------------------------------------
!
USE kinds, ONLY : DP
USE io_global, ONLY : ionode
USE io_dyn_mat, ONLY : write_dyn_mat_header, write_ifc
USE cell_base, ONLY : celldm
USE constants, ONLY : amu_ry
USE ifconstants,ONLY : tau => tau_blk, ityp => ityp_blk, &
frc, frc_lr, zeu, atm, m_loc
!
IMPLICIT NONE
! I/O variable
CHARACTER(LEN=256) :: flfrc
INTEGER :: nr1, nr2, nr3, ibrav, nat, ntyp, nspin_mag, nqs
REAL(DP) :: alat, omega, at(3,3), bg(3,3), amass(ntyp), epsil(3,3)
REAL(DP) :: alph
LOGICAL :: xmlifc, has_zstar
! local variables
INTEGER :: i, j, nt, na, nb, nn, m1, m2, m3, leng
CHARACTER(LEN=256) :: flfrc2
COMPLEX(DP) :: frc2(nr1*nr2*nr3,3,3,nat,nat)
!
leng = LEN_TRIM(flfrc)
flfrc2 = flfrc(1:leng) // ".matdyn"
IF (xmlifc) THEN
frc2 = RESHAPE(frc, (/nr1*nr2*nr3,3,3,nat,nat/))
IF (has_zstar) THEN
CALL write_dyn_mat_header(flfrc2, ntyp, nat, ibrav, nspin_mag, celldm, &
at, bg, omega, atm, amass, tau, ityp, m_loc, nqs, epsil, zeu)
ELSE
CALL write_dyn_mat_header(flfrc2, ntyp, nat, ibrav, nspin_mag, &
celldm, at, bg, omega, atm, amass, tau, ityp, m_loc, nqs)
ENDIF
CALL write_ifc(alph,nr1, nr2, nr3, nat, frc2)
ELSE IF (ionode) THEN
OPEN(unit=1, file=flfrc2, status='unknown', form='formatted')
WRITE(1,'(i3,i5,i4,6f11.7)') ntyp, nat, ibrav, celldm
IF (ibrav == 0) WRITE(1,'(2x,3f15.9)') ((at(i,j),i=1,3),j=1,3)
DO nt = 1, ntyp
WRITE(1,*) nt, " '", atm(nt), "' ", amass(nt)*amu_ry
END DO
DO na = 1, nat
WRITE(1,'(2i5,3f18.10)') na, ityp(na), (tau(j,na),j=1,3)
END DO
WRITE (1,*) has_zstar, alph
IF (has_zstar) THEN
WRITE(1,'(3f24.12)') ((epsil(i,j),j=1,3),i=1,3)
DO na = 1, nat
WRITE(1,'(i5)') na
WRITE(1,'(3f15.7)') ((zeu(i,j,na),j=1,3),i=1,3)
END DO
END IF
WRITE (1,'(4i4)') nr1, nr2, nr3
DO i = 1, 3
DO j = 1,3
DO na = 1, nat
DO nb = 1, nat
WRITE (1,'(4i4)') i, j, na, nb
nn = 0
DO m3 = 1, nr3
DO m2 = 1, nr2
DO m1 = 1, nr1
nn = nn+1
WRITE (1,'(3i4,2x,1pe18.11)') &
m1,m2,m3, DBLE(frc(m1,m2,m3,i,j,na,nb))
END DO
END DO
END DO
END DO
END DO
END DO
END DO
CLOSE(1)
END IF
END SUBROUTINE writefc
!
!-----------------------------------------------------------------------
SUBROUTINE frc_blk(dyn,q,tau,nat,nr1,nr2,nr3,frc,at,bg,rws,nrws,f_of_q,fd)
!-----------------------------------------------------------------------
! calculates the dynamical matrix at q from the (short-range part of the)
! force constants
!
USE kinds, ONLY : DP
USE constants, ONLY : tpi
USE io_global, ONLY : stdout
!
IMPLICIT NONE
INTEGER nr1, nr2, nr3, nat, n1, n2, n3, nr1_, nr2_, nr3_, &
ipol, jpol, na, nb, m1, m2, m3, nint, i,j, nrws, nax
COMPLEX(DP) dyn(3,3,nat,nat), f_of_q(3,3,nat,nat)
REAL(DP) frc(nr1,nr2,nr3,3,3,nat,nat), tau(3,nat), q(3), arg, &
at(3,3), bg(3,3), r(3), weight, r_ws(3), &
total_weight, rws(0:3,nrws), alat
REAL(DP), EXTERNAL :: wsweight
REAL(DP),SAVE,ALLOCATABLE :: wscache(:,:,:,:,:)
LOGICAL,SAVE :: first=.true.
LOGICAL :: fd
!
nr1_=2*nr1
nr2_=2*nr2
nr3_=2*nr3
FIRST_TIME : IF (first) THEN
first=.false.
ALLOCATE( wscache(-nr3_:nr3_, -nr2_:nr2_, -nr1_:nr1_, nat,nat) )
DO na=1, nat
DO nb=1, nat
total_weight=0.0d0
!
! SUM OVER R VECTORS IN THE SUPERCELL - VERY VERY VERY SAFE RANGE!
!
DO n1=-nr1_,nr1_
DO n2=-nr2_,nr2_
DO n3=-nr3_,nr3_
DO i=1, 3
r(i) = n1*at(i,1)+n2*at(i,2)+n3*at(i,3)
r_ws(i) = r(i) + tau(i,na)-tau(i,nb)
if (fd) r_ws(i) = r(i) + tau(i,nb)-tau(i,na)
END DO
wscache(n3,n2,n1,nb,na) = wsweight(r_ws,rws,nrws)
total_weight=total_weight + wscache(n3,n2,n1,nb,na)
ENDDO
ENDDO
ENDDO
IF (ABS(total_weight-nr1*nr2*nr3).GT.1.0d-8) THEN
WRITE(stdout,*) na,nb,total_weight
CALL errore ('frc_blk','wrong total_weight',1)
END IF
ENDDO
ENDDO
ENDIF FIRST_TIME
!
DO na=1, nat
DO nb=1, nat
DO n1=-nr1_,nr1_
DO n2=-nr2_,nr2_
DO n3=-nr3_,nr3_
!
! SUM OVER R VECTORS IN THE SUPERCELL - VERY VERY SAFE RANGE!
!
DO i=1, 3
r(i) = n1*at(i,1)+n2*at(i,2)+n3*at(i,3)
END DO
weight = wscache(n3,n2,n1,nb,na)
IF (weight .GT. 0.0d0) THEN
!
! FIND THE VECTOR CORRESPONDING TO R IN THE ORIGINAL CELL
!
m1 = MOD(n1+1,nr1)
IF(m1.LE.0) m1=m1+nr1
m2 = MOD(n2+1,nr2)
IF(m2.LE.0) m2=m2+nr2
m3 = MOD(n3+1,nr3)
IF(m3.LE.0) m3=m3+nr3
! write(*,'(6i4)') n1,n2,n3,m1,m2,m3
!
! FOURIER TRANSFORM
!
arg = tpi*(q(1)*r(1) + q(2)*r(2) + q(3)*r(3))
DO ipol=1, 3
DO jpol=1, 3
dyn(ipol,jpol,na,nb) = dyn(ipol,jpol,na,nb) + &
(frc(m1,m2,m3,ipol,jpol,na,nb)+f_of_q(ipol,jpol,na,nb)) &
*CMPLX(COS(arg),-SIN(arg),kind=DP)*weight
END DO
END DO
END IF
END DO
END DO
END DO
END DO
END DO
!
RETURN
END SUBROUTINE frc_blk
!
!-----------------------------------------------------------------------
SUBROUTINE setupmat (q,dyn,nat,at,bg,tau,itau_blk,nsc,alat, &
& dyn_blk,nat_blk,at_blk,bg_blk,tau_blk,omega_blk, &
& loto_2d, epsil, zeu, alph, &
& frc,nr1,nr2,nr3,has_zstar,rws,nrws,na_ifc,f_of_q,fd)
!-----------------------------------------------------------------------
! compute the dynamical matrix (the analytic part only)
!
USE kinds, ONLY : DP
USE constants, ONLY : tpi
USE cell_base, ONLY : celldm
USE rigid, ONLY : rgd_blk
!
IMPLICIT NONE
!
! I/O variables
!
INTEGER:: nr1, nr2, nr3, nat, nat_blk, nsc, nrws, itau_blk(nat)
REAL(DP) :: q(3), tau(3,nat), at(3,3), bg(3,3), alat, alph, &
epsil(3,3), zeu(3,3,nat_blk), rws(0:3,nrws), &
frc(nr1,nr2,nr3,3,3,nat_blk,nat_blk)
REAL(DP) :: tau_blk(3,nat_blk), at_blk(3,3), bg_blk(3,3), omega_blk
COMPLEX(DP) dyn_blk(3,3,nat_blk,nat_blk), f_of_q(3,3,nat,nat)
COMPLEX(DP) :: dyn(3,3,nat,nat)
LOGICAL :: has_zstar, na_ifc, fd, loto_2d
!
! local variables
!
REAL(DP) :: arg
COMPLEX(DP) :: cfac(nat)
INTEGER :: i,j,k, na,nb, na_blk, nb_blk, iq
REAL(DP) :: qp(3), qbid(3,nsc) ! automatic array
!
!
CALL q_gen(nsc,qbid,at_blk,bg_blk,at,bg)
!
DO iq=1,nsc
!
DO k=1,3
qp(k)= q(k) + qbid(k,iq)
END DO
!
dyn_blk(:,:,:,:) = (0.d0,0.d0)
CALL frc_blk (dyn_blk,qp,tau_blk,nat_blk, &
& nr1,nr2,nr3,frc,at_blk,bg_blk,rws,nrws,f_of_q,fd)
IF (has_zstar .and. .not.na_ifc) &
CALL rgd_blk(nr1, nr2, nr3, nat_blk, dyn_blk, qp, tau_blk, &
epsil, zeu, alph, bg_blk, omega_blk, celldm(1), loto_2d, +1.d0 )
!
DO na=1,nat
na_blk = itau_blk(na)
DO nb=1,nat
nb_blk = itau_blk(nb)
!
arg=tpi* ( qp(1) * ( (tau(1,na)-tau_blk(1,na_blk)) - &
(tau(1,nb)-tau_blk(1,nb_blk)) ) + &
qp(2) * ( (tau(2,na)-tau_blk(2,na_blk)) - &
(tau(2,nb)-tau_blk(2,nb_blk)) ) + &
qp(3) * ( (tau(3,na)-tau_blk(3,na_blk)) - &
(tau(3,nb)-tau_blk(3,nb_blk)) ) )
!
cfac(nb) = CMPLX(COS(arg),SIN(arg),kind=DP)/nsc
!
END DO ! nb
!
DO i=1,3
DO j=1,3
!
DO nb=1,nat
nb_blk = itau_blk(nb)
dyn(i,j,na,nb) = dyn(i,j,na,nb) + cfac(nb) * &
dyn_blk(i,j,na_blk,nb_blk)
END DO ! nb
!
END DO ! j
END DO ! i
END DO ! na
!
END DO ! iq
!
RETURN
END SUBROUTINE setupmat
!
!
!----------------------------------------------------------------------
SUBROUTINE set_asr (asr, nr1, nr2, nr3, frc, frc_lr, zeu, nat, ibrav, tau_blk, at_blk, huang)
!-----------------------------------------------------------------------
!
USE kinds, ONLY : DP
USE io_global, ONLY : stdout
!
IMPLICIT NONE
CHARACTER (LEN=10), intent(in) :: asr
LOGICAL, intent(in) :: huang
INTEGER, intent(in) :: nr1, nr2, nr3, nat, ibrav
REAL(DP), intent(in) :: tau_blk(3,nat), at_blk(3,3), frc_lr(nr1,nr2,nr3,3,3,nat,nat)
REAL(DP), intent(inout) :: frc(nr1,nr2,nr3,3,3,nat,nat), zeu(3,3,nat)
!
INTEGER :: axis, n, i, j, na, nb, n1,n2,n3, m,p,k,l,q,r, i1,j1,na1, ip,ieq,neq
INTEGER :: huang_set(4,15)
REAL(DP) :: zeu_new(3,3,nat), tau(3,48), rcell(3), r_ws(3), eps
REAL(DP), ALLOCATABLE :: frc_new(:,:,:,:,:,:,:)
parameter (eps=1.0d-8)
type vector
real(DP),pointer :: vec(:,:,:,:,:,:,:)
end type vector
!
type (vector) u(6*3*nat+15)
! These are the "vectors" associated with the sum rules on force-constants
!
integer :: u_less(6*3*nat),n_less,i_less
! indices of the vectors u that are not independent to the preceding ones,
! n_less = number of such vectors, i_less = temporary parameter
!
integer, allocatable :: ind_v(:,:,:)
real(DP), allocatable :: v(:,:)
! These are the "vectors" associated with symmetry conditions, coded by
! indicating the positions (i.e. the seven indices) of the non-zero elements (there
! should be only 2 of them) and the value of that element. We do so in order
! to limit the amount of memory used.
!
real(DP), allocatable :: w(:,:,:,:,:,:,:), x(:,:,:,:,:,:,:)
! temporary vectors and parameters
real(DP) :: scal,norm2, sum
!
real(DP) :: zeu_u(6*3,3,3,nat)
! These are the "vectors" associated with the sum rules on effective charges
!
integer :: zeu_less(6*3),nzeu_less,izeu_less
! indices of the vectors zeu_u that are not independent to the preceding ones,
! nzeu_less = number of such vectors, izeu_less = temporary parameter
!
real(DP) :: zeu_w(3,3,nat), zeu_x(3,3,nat)
! temporary vectors
! Initialization. n is the number of sum rules to be considered (if asr.ne.'simple')
! and 'axis' is the rotation axis in the case of a 1D system
! (i.e. the rotation axis is (Ox) if axis='1', (Oy) if axis='2' and (Oz) if axis='3')
!
if((asr.ne.'simple').and.(asr.ne.'crystal').and.(asr.ne.'one-dim') &
.and.(asr.ne.'zero-dim').and.(asr.ne.'all')) then
call errore('set_asr','invalid Acoustic Sum Rule:' // asr, 1)
endif
!
if(asr.eq.'simple') then
!
! Simple Acoustic Sum Rule on effective charges
!
do i=1,3
do j=1,3
sum=0.0d0
do na=1,nat
sum = sum + zeu(i,j,na)
end do
do na=1,nat
zeu(i,j,na) = zeu(i,j,na) - sum/nat
end do
end do
end do
!
! Simple Acoustic Sum Rule on force constants in real space
!
do i=1,3
do j=1,3
do na=1,nat
sum=0.0d0
do nb=1,nat
do n1=1,nr1
do n2=1,nr2
do n3=1,nr3
sum=sum+frc(n1,n2,n3,i,j,na,nb)
end do
end do
end do
end do
frc(1,1,1,i,j,na,na) = frc(1,1,1,i,j,na,na) - sum
! write(6,*) ' na, i, j, sum = ',na,i,j,sum
end do
end do
end do
!
return
!
end if
if(asr.eq.'crystal') n=3
if(asr.eq.'one-dim') then
! the direction of periodicity is the rotation axis
! It will work only if the crystal axis considered is one of
! the cartesian axis (typically, ibrav=1, 6 or 8, or 4 along the
! z-direction)
if (nr1*nr2*nr3.eq.1) axis=3
if ((nr1.ne.1).and.(nr2*nr3.eq.1)) axis=1
if ((nr2.ne.1).and.(nr1*nr3.eq.1)) axis=2
if ((nr3.ne.1).and.(nr1*nr2.eq.1)) axis=3
if (((nr1.ne.1).and.(nr2.ne.1)).or.((nr2.ne.1).and. &
(nr3.ne.1)).or.((nr1.ne.1).and.(nr3.ne.1))) then
call errore('set_asr','too many directions of &
& periodicity in 1D system',axis)
endif
if ((ibrav.ne.1).and.(ibrav.ne.6).and.(ibrav.ne.8).and. &
((ibrav.ne.4).or.(axis.ne.3)) ) then
write(stdout,*) 'asr: rotational axis may be wrong'
endif
write(stdout,'("asr rotation axis in 1D system= ",I4)') axis
n=4
endif
if(asr.eq.'zero-dim') n=6
if(asr.eq.'all') n=21
!
! Acoustic Sum Rule on effective charges
!
! generating the vectors of the orthogonal of the subspace to project
! the effective charges matrix on
!
zeu_u(:,:,:,:)=0.0d0
do i=1,3
do j=1,3
do na=1,nat
zeu_new(i,j,na)=zeu(i,j,na)
enddo
enddo
enddo
!
p=0
do i=1,3
do j=1,3
! These are the 3*3 vectors associated with the
! translational acoustic sum rules
p=p+1
zeu_u(p,i,j,:)=1.0d0
!
enddo
enddo
!
if (n.eq.4) then
do i=1,3
! These are the 3 vectors associated with the
! single rotational sum rule (1D system)
p=p+1
do na=1,nat
zeu_u(p,i,MOD(axis,3)+1,na)=-tau_blk(MOD(axis+1,3)+1,na)
zeu_u(p,i,MOD(axis+1,3)+1,na)=tau_blk(MOD(axis,3)+1,na)
enddo
!
enddo
endif
!
if (n.eq.6) then
do i=1,3
do j=1,3
! These are the 3*3 vectors associated with the
! three rotational sum rules (0D system - typ. molecule)
p=p+1
do na=1,nat
zeu_u(p,i,MOD(j,3)+1,na)=-tau_blk(MOD(j+1,3)+1,na)
zeu_u(p,i,MOD(j+1,3)+1,na)=tau_blk(MOD(j,3)+1,na)
enddo
!
enddo
enddo
endif
!
! Gram-Schmidt orthonormalization of the set of vectors created.
!
nzeu_less=0
do k=1,p
zeu_w(:,:,:)=zeu_u(k,:,:,:)
zeu_x(:,:,:)=zeu_u(k,:,:,:)
do q=1,k-1
r=1
do izeu_less=1,nzeu_less
if (zeu_less(izeu_less).eq.q) r=0
enddo
if (r.ne.0) then
call sp_zeu(zeu_x,zeu_u(q,:,:,:),nat,scal)
zeu_w(:,:,:) = zeu_w(:,:,:) - scal* zeu_u(q,:,:,:)
endif
enddo
call sp_zeu(zeu_w,zeu_w,nat,norm2)
if (norm2.gt.1.0d-16) then
zeu_u(k,:,:,:) = zeu_w(:,:,:) / DSQRT(norm2)
else
nzeu_less=nzeu_less+1
zeu_less(nzeu_less)=k
endif
enddo
!
! Projection of the effective charge "vector" on the orthogonal of the
! subspace of the vectors verifying the sum rules
!
zeu_w(:,:,:)=0.0d0
do k=1,p
r=1
do izeu_less=1,nzeu_less
if (zeu_less(izeu_less).eq.k) r=0
enddo
if (r.ne.0) then
zeu_x(:,:,:)=zeu_u(k,:,:,:)
call sp_zeu(zeu_x,zeu_new,nat,scal)
zeu_w(:,:,:) = zeu_w(:,:,:) + scal*zeu_u(k,:,:,:)
endif
enddo
!
! Final substraction of the former projection to the initial zeu, to get
! the new "projected" zeu
!
zeu_new(:,:,:)=zeu_new(:,:,:) - zeu_w(:,:,:)
call sp_zeu(zeu_w,zeu_w,nat,norm2)
write(stdout,'(" Norm of the difference between old and new effective ", &
& "charges: ",F25.20/)') SQRT(norm2)
!
! Check projection
!
!write(6,'("Check projection of zeu")')
!do k=1,p
! zeu_x(:,:,:)=zeu_u(k,:,:,:)
! call sp_zeu(zeu_x,zeu_new,nat,scal)
! if (DABS(scal).gt.1d-10) write(6,'("k= ",I8," zeu_new|zeu_u(k)= ",F15.10)') k,scal
!enddo
!
do i=1,3
do j=1,3
do na=1,nat
zeu(i,j,na)=zeu_new(i,j,na)
enddo
enddo
enddo
!
! Acoustic Sum Rule on force constants
!
!
! generating the vectors of the orthogonal of the subspace to project
! the force-constants matrix on
!
do k=1,18*nat+15
allocate(u(k) % vec(nr1,nr2,nr3,3,3,nat,nat))
u(k) % vec (:,:,:,:,:,:,:)=0.0d0
enddo
ALLOCATE (frc_new(nr1,nr2,nr3,3,3,nat,nat))
if (asr.eq.'all') then
frc_new(:,:,:,:,:,:,:)=frc(:,:,:,:,:,:,:)+frc_lr(:,:,:,:,:,:,:)
else
frc_new(:,:,:,:,:,:,:)=frc(:,:,:,:,:,:,:)
endif
!
p=0
do i=1,3
do j=1,3
do na=1,nat
! These are the 3*3*nat vectors associated with the
! translational acoustic sum rules
p=p+1
u(p) % vec (:,:,:,i,j,na,:)=1.0d0
!
enddo
enddo
enddo
!
if (n.eq.4) then
do i=1,3
do na=1,nat
! These are the 3*nat vectors associated with the
! single rotational sum rule (1D system)
p=p+1
do nb=1,nat
u(p) % vec (:,:,:,i,MOD(axis,3)+1,na,nb)=-tau_blk(MOD(axis+1,3)+1,nb)
u(p) % vec (:,:,:,i,MOD(axis+1,3)+1,na,nb)=tau_blk(MOD(axis,3)+1,nb)
enddo
!
enddo
enddo
endif
!
if (n.eq.6) then
do i=1,3
do j=1,3
do na=1,nat
! These are the 3*3*nat vectors associated with the
! three rotational sum rules (0D system - typ. molecule)
p=p+1
do nb=1,nat
u(p) % vec (:,:,:,i,MOD(j,3)+1,na,nb)=-tau_blk(MOD(j+1,3)+1,nb)
u(p) % vec (:,:,:,i,MOD(j+1,3)+1,na,nb)=tau_blk(MOD(j,3)+1,nb)
enddo
!
enddo
enddo
enddo
endif
!
if (n.eq.21) then
!
!! Please consider citing C. Lin, S. Ponc\'e and N. Marzari, npj Comput Mater 8, 236 (2022)
!! if asr='all' is used.
!
! Born-Huang invariance conditions
do i=1,3
do j=1,3
do na=1,nat
! These are 3*3*nat vectors associated with the three
! rotational sum rules (valid for system of any dimension).
! These differ from the case of n=6 (zero-dim), where atom b is
! in supercell and its position is wrapped to its nearest periodic
! image with the correct weight in the case of degeneracy.
p=p+1
do nb=1,nat
do n1=1,nr1
do n2=1,nr2
do n3=1,nr3
rcell=matmul(at_blk,(/n1,n2,n3/)-(/1,1,1/))
r_ws=rcell-tau_blk(:,nb)+tau_blk(:,na)
call ws_all(tau,neq,nr1,nr2,nr3,r_ws,at_blk)
do ieq=1,neq
u(p) % vec (n1,n2,n3,i,MOD(j,3)+1,na,nb)=&
u(p) % vec (n1,n2,n3,i,MOD(j,3)+1,na,nb)-(tau(MOD(j+1,3)+1,ieq)-tau_blk(MOD(j+1,3)+1,na))/DBLE(neq)
u(p) % vec (n1,n2,n3,i,MOD(j+1,3)+1,na,nb)=&
u(p) % vec (n1,n2,n3,i,MOD(j+1,3)+1,na,nb)+(tau(MOD(j,3)+1,ieq)-tau_blk(MOD(j,3)+1,na))/DBLE(neq)
enddo
enddo
enddo
enddo
enddo
!call sp1(frc_new, u(p) % vec (:,:,:,:,:,:,:), nr1,nr2,nr3,nat,scal)
!write(*,*) scal
!
where(abs(u(p) % vec)<eps) u(p) % vec=0.d0
enddo
enddo
enddo
!
! Huang conditions
if (huang) then
do na=1,nat
do nb=1,nat
do n1=1,nr1
do n2=1,nr2
do n3=1,nr3
! These are 15 vectors for vanishing stress tensor
rcell=matmul(at_blk,(/n1,n2,n3/)-(/1,1,1/))
r_ws=rcell-tau_blk(:,nb)+tau_blk(:,na)
call ws_all(tau,neq,nr1,nr2,nr3,r_ws,at_blk)
!
do ieq=1,neq
! 1 1 1 2, yx
huang_set(:,1)=(/1,1,1,2/)
u(p+1) % vec (n1,n2,n3,1,1,na,nb)=&
u(p+1) % vec (n1,n2,n3,1,1,na,nb)-tau(1,ieq)*tau(2,ieq)/DBLE(neq)
u(p+1) % vec (n1,n2,n3,1,2,na,nb)=&
u(p+1) % vec (n1,n2,n3,1,2,na,nb)+tau(1,ieq)*tau(1,ieq)/DBLE(neq)
! 1 1 1 3, zx
huang_set(:,2)=(/1,1,1,3/)
u(p+2) % vec (n1,n2,n3,1,1,na,nb)=&
u(p+2) % vec (n1,n2,n3,1,1,na,nb)-tau(1,ieq)*tau(3,ieq)/DBLE(neq)
u(p+2) % vec (n1,n2,n3,1,3,na,nb)=&
u(p+2) % vec (n1,n2,n3,1,3,na,nb)+tau(1,ieq)*tau(1,ieq)/DBLE(neq)
! 1 1 2 2, xx-yy
huang_set(:,3)=(/1,1,2,2/)
u(p+3) % vec (n1,n2,n3,1,1,na,nb)=&
u(p+3) % vec (n1,n2,n3,1,1,na,nb)-tau(2,ieq)*tau(2,ieq)/DBLE(neq)
u(p+3) % vec (n1,n2,n3,2,2,na,nb)=&
u(p+3) % vec (n1,n2,n3,2,2,na,nb)+tau(1,ieq)*tau(1,ieq)/DBLE(neq)
! 1 1 2 3
huang_set(:,4)=(/1,1,2,3/)
u(p+4) % vec (n1,n2,n3,1,1,na,nb)=&
u(p+4) % vec (n1,n2,n3,1,1,na,nb)-tau(2,ieq)*tau(3,ieq)/DBLE(neq)
u(p+4) % vec (n1,n2,n3,2,3,na,nb)=&
u(p+4) % vec (n1,n2,n3,2,3,na,nb)+tau(1,ieq)*tau(1,ieq)/DBLE(neq)
! 1 1 3 3, xx-zz
huang_set(:,5)=(/1,1,3,3/)
u(p+5) % vec (n1,n2,n3,1,1,na,nb)=&
u(p+5) % vec (n1,n2,n3,1,1,na,nb)-tau(3,ieq)*tau(3,ieq)/DBLE(neq)
u(p+5) % vec (n1,n2,n3,3,3,na,nb)=&
u(p+5) % vec (n1,n2,n3,3,3,na,nb)+tau(1,ieq)*tau(1,ieq)/DBLE(neq)
! 1 2 1 3
huang_set(:,6)=(/1,2,1,3/)
u(p+6) % vec (n1,n2,n3,1,2,na,nb)=&
u(p+6) % vec (n1,n2,n3,1,2,na,nb)-tau(1,ieq)*tau(3,ieq)/DBLE(neq)
u(p+6) % vec (n1,n2,n3,1,3,na,nb)=&
u(p+6) % vec (n1,n2,n3,1,3,na,nb)+tau(1,ieq)*tau(2,ieq)/DBLE(neq)
! 1 2 2 2, xy
huang_set(:,7)=(/1,2,2,2/)
u(p+7) % vec (n1,n2,n3,1,2,na,nb)=&
u(p+7) % vec (n1,n2,n3,1,2,na,nb)-tau(2,ieq)*tau(2,ieq)/DBLE(neq)
u(p+7) % vec (n1,n2,n3,2,2,na,nb)=&
u(p+7) % vec (n1,n2,n3,2,2,na,nb)+tau(1,ieq)*tau(2,ieq)/DBLE(neq)
! 1 2 2 3
huang_set(:,8)=(/1,2,2,3/)
u(p+8) % vec (n1,n2,n3,1,2,na,nb)=&
u(p+8) % vec (n1,n2,n3,1,2,na,nb)-tau(2,ieq)*tau(3,ieq)/DBLE(neq)
u(p+8) % vec (n1,n2,n3,2,3,na,nb)=&
u(p+8) % vec (n1,n2,n3,2,3,na,nb)+tau(1,ieq)*tau(2,ieq)/DBLE(neq)
! 1 2 3 3
huang_set(:,9)=(/1,2,3,3/)
u(p+9) % vec (n1,n2,n3,1,2,na,nb)=&
u(p+9) % vec (n1,n2,n3,1,2,na,nb)-tau(3,ieq)*tau(3,ieq)/DBLE(neq)
u(p+9) % vec (n1,n2,n3,3,3,na,nb)=&
u(p+9) % vec (n1,n2,n3,3,3,na,nb)+tau(1,ieq)*tau(2,ieq)/DBLE(neq)
! 1 3 2 2
huang_set(:,10)=(/1,3,2,2/)
u(p+10) % vec (n1,n2,n3,1,3,na,nb)=&
u(p+10) % vec (n1,n2,n3,1,3,na,nb)-tau(2,ieq)*tau(2,ieq)/DBLE(neq)
u(p+10) % vec (n1,n2,n3,2,2,na,nb)=&
u(p+10) % vec (n1,n2,n3,2,2,na,nb)+tau(1,ieq)*tau(3,ieq)/DBLE(neq)
! 1 3 2 3
huang_set(:,11)=(/1,3,2,3/)
u(p+11) % vec (n1,n2,n3,1,3,na,nb)=&
u(p+11) % vec (n1,n2,n3,1,3,na,nb)-tau(2,ieq)*tau(3,ieq)/DBLE(neq)
u(p+11) % vec (n1,n2,n3,2,3,na,nb)=&
u(p+11) % vec (n1,n2,n3,2,3,na,nb)+tau(1,ieq)*tau(3,ieq)/DBLE(neq)
! 1 3 3 3, xz
huang_set(:,12)=(/1,3,3,3/)
u(p+12) % vec (n1,n2,n3,1,3,na,nb)=&
u(p+12) % vec (n1,n2,n3,1,3,na,nb)-tau(3,ieq)*tau(3,ieq)/DBLE(neq)
u(p+12) % vec (n1,n2,n3,3,3,na,nb)=&
u(p+12) % vec (n1,n2,n3,3,3,na,nb)+tau(1,ieq)*tau(3,ieq)/DBLE(neq)
! 2 2 2 3, zy
huang_set(:,13)=(/2,2,2,3/)
u(p+13) % vec (n1,n2,n3,2,2,na,nb)=&
u(p+13) % vec (n1,n2,n3,2,2,na,nb)-tau(2,ieq)*tau(3,ieq)/DBLE(neq)
u(p+13) % vec (n1,n2,n3,2,3,na,nb)=&
u(p+13) % vec (n1,n2,n3,2,3,na,nb)+tau(2,ieq)*tau(2,ieq)/DBLE(neq)
! 2 2 3 3, yy-zz
huang_set(:,14)=(/2,2,3,3/)
u(p+14) % vec (n1,n2,n3,2,2,na,nb)=&
u(p+14) % vec (n1,n2,n3,2,2,na,nb)-tau(3,ieq)*tau(3,ieq)/DBLE(neq)
u(p+14) % vec (n1,n2,n3,3,3,na,nb)=&
u(p+14) % vec (n1,n2,n3,3,3,na,nb)+tau(2,ieq)*tau(2,ieq)/DBLE(neq)
! 2 3 3 3, yz
huang_set(:,15)=(/2,3,3,3/)
u(p+15) % vec (n1,n2,n3,2,3,na,nb)=&
u(p+15) % vec (n1,n2,n3,2,3,na,nb)-tau(3,ieq)*tau(3,ieq)/DBLE(neq)
u(p+15) % vec (n1,n2,n3,3,3,na,nb)=&
u(p+15) % vec (n1,n2,n3,3,3,na,nb)+tau(2,ieq)*tau(3,ieq)/DBLE(neq)
enddo
!
enddo
enddo
enddo
enddo
enddo
!
do ip=1,15
where(abs(u(p+ip) % vec)<eps) u(p+ip) % vec=0.d0
call sp1(u(p+ip) % vec (:,:,:,:,:,:,:), frc_new, nr1,nr2,nr3,nat,scal)
write(stdout,'(A, 4I4, A, F20.16)') " Huang before: ", huang_set(:,ip), ", residual stress: ", scal
enddo
write(stdout,*)
p=p+15
endif
!
endif
!
allocate (ind_v(9*nat*nat*nr1*nr2*nr3,2,7), v(9*nat*nat*nr1*nr2*nr3,2) )
m=0
do i=1,3
do j=1,3
do na=1,nat
do nb=1,nat
do n1=1,nr1
do n2=1,nr2
do n3=1,nr3
! These are the vectors associated with the symmetry constraints
q=1
l=1
do while((l.le.m).and.(q.ne.0))
if ((ind_v(l,1,1).eq.n1).and.(ind_v(l,1,2).eq.n2).and. &
(ind_v(l,1,3).eq.n3).and.(ind_v(l,1,4).eq.i).and. &
(ind_v(l,1,5).eq.j).and.(ind_v(l,1,6).eq.na).and. &
(ind_v(l,1,7).eq.nb)) q=0
if ((ind_v(l,2,1).eq.n1).and.(ind_v(l,2,2).eq.n2).and. &
(ind_v(l,2,3).eq.n3).and.(ind_v(l,2,4).eq.i).and. &
(ind_v(l,2,5).eq.j).and.(ind_v(l,2,6).eq.na).and. &
(ind_v(l,2,7).eq.nb)) q=0
l=l+1
enddo
if ((n1.eq.MOD(nr1+1-n1,nr1)+1).and.(n2.eq.MOD(nr2+1-n2,nr2)+1) &
.and.(n3.eq.MOD(nr3+1-n3,nr3)+1).and.(i.eq.j).and.(na.eq.nb)) q=0
if (q.ne.0) then
m=m+1
ind_v(m,1,1)=n1
ind_v(m,1,2)=n2
ind_v(m,1,3)=n3
ind_v(m,1,4)=i
ind_v(m,1,5)=j
ind_v(m,1,6)=na
ind_v(m,1,7)=nb
v(m,1)=1.0d0/DSQRT(2.0d0)
ind_v(m,2,1)=MOD(nr1+1-n1,nr1)+1
ind_v(m,2,2)=MOD(nr2+1-n2,nr2)+1
ind_v(m,2,3)=MOD(nr3+1-n3,nr3)+1
ind_v(m,2,4)=j
ind_v(m,2,5)=i
ind_v(m,2,6)=nb
ind_v(m,2,7)=na
v(m,2)=-1.0d0/DSQRT(2.0d0)
endif
enddo
enddo
enddo
enddo
enddo
enddo
enddo
!
! Gram-Schmidt orthonormalization of the set of vectors created.
! Note that the vectors corresponding to symmetry constraints are already
! orthonormalized by construction.
!
n_less=0
allocate (w(nr1,nr2,nr3,3,3,nat,nat), x(nr1,nr2,nr3,3,3,nat,nat))
do k=1,p
w(:,:,:,:,:,:,:)=u(k) % vec (:,:,:,:,:,:,:)
x(:,:,:,:,:,:,:)=u(k) % vec (:,:,:,:,:,:,:)
do l=1,m
!
call sp2(x,v(l,:),ind_v(l,:,:),nr1,nr2,nr3,nat,scal)
do r=1,2
n1=ind_v(l,r,1)
n2=ind_v(l,r,2)
n3=ind_v(l,r,3)
i=ind_v(l,r,4)
j=ind_v(l,r,5)
na=ind_v(l,r,6)
nb=ind_v(l,r,7)
w(n1,n2,n3,i,j,na,nb)=w(n1,n2,n3,i,j,na,nb)-scal*v(l,r)
enddo
enddo
if (k.le.(9*nat)) then
na1=MOD(k,nat)
if (na1.eq.0) na1=nat
j1=MOD((k-na1)/nat,3)+1
i1=MOD((((k-na1)/nat)-j1+1)/3,3)+1
else
q=k-9*nat
if (k.le.(18*nat)) then
if (n.eq.4) then
na1=MOD(q,nat)
if (na1.eq.0) na1=nat
i1=MOD((q-na1)/nat,3)+1
else
na1=MOD(q,nat)
if (na1.eq.0) na1=nat
j1=MOD((q-na1)/nat,3)+1
i1=MOD((((q-na1)/nat)-j1+1)/3,3)+1
endif
endif
endif
do q=1,k-1
r=1
do i_less=1,n_less
if (u_less(i_less).eq.q) r=0
enddo
if (r.ne.0) then
if (k.le.(18*nat)) then
call sp3(x,u(q) % vec (:,:,:,:,:,:,:), i1,na1,nr1,nr2,nr3,nat,scal)
else
call sp1(x,u(q) % vec (:,:,:,:,:,:,:), nr1,nr2,nr3,nat,scal)
endif
w(:,:,:,:,:,:,:) = w(:,:,:,:,:,:,:) - scal * u(q) % vec (:,:,:,:,:,:,:)
endif
enddo
call sp1(w,w,nr1,nr2,nr3,nat,norm2)
if (norm2.gt.eps) then
u(k) % vec (:,:,:,:,:,:,:) = w(:,:,:,:,:,:,:) / DSQRT(norm2)
else
n_less=n_less+1
u_less(n_less)=k
endif
enddo
!
! Projection of the force-constants "vector" on the orthogonal of the
! subspace of the vectors verifying the sum rules and symmetry contraints
!
w(:,:,:,:,:,:,:)=0.0d0
do l=1,m
call sp2(frc_new,v(l,:),ind_v(l,:,:),nr1,nr2,nr3,nat,scal)
do r=1,2
n1=ind_v(l,r,1)
n2=ind_v(l,r,2)
n3=ind_v(l,r,3)
i=ind_v(l,r,4)
j=ind_v(l,r,5)
na=ind_v(l,r,6)
nb=ind_v(l,r,7)
w(n1,n2,n3,i,j,na,nb)=w(n1,n2,n3,i,j,na,nb)+scal*v(l,r)
enddo
enddo
do k=1,p
r=1
do i_less=1,n_less
if (u_less(i_less).eq.k) r=0
enddo
if (r.ne.0) then
x(:,:,:,:,:,:,:)=u(k) % vec (:,:,:,:,:,:,:)
call sp1(x,frc_new,nr1,nr2,nr3,nat,scal)
w(:,:,:,:,:,:,:) = w(:,:,:,:,:,:,:) + scal*u(k)%vec(:,:,:,:,:,:,:)
endif
deallocate(u(k) % vec)
enddo
!
! Final substraction of the former projection to the initial frc, to get
! the new "projected" frc
!
frc_new(:,:,:,:,:,:,:)=frc_new(:,:,:,:,:,:,:) - w(:,:,:,:,:,:,:)
if (asr.eq.'all') frc_new(:,:,:,:,:,:,:)=frc_new(:,:,:,:,:,:,:) - frc_lr(:,:,:,:,:,:,:)
call sp1(w,w,nr1,nr2,nr3,nat,norm2)
write(stdout,'(" Norm of the difference between old and new force-constants:",&
& F25.20/)') SQRT(norm2)
!
! Check projection
!
!write(6,'("Check projection IFC")')
!do l=1,m
! call sp2(frc_new,v(l,:),ind_v(l,:,:),nr1,nr2,nr3,nat,scal)
! if (DABS(scal).gt.1d-10) write(6,'("l= ",I8," frc_new|v(l)= ",F15.10)') l,scal
!enddo
!do k=1,p
! x(:,:,:,:,:,:,:)=u(k) % vec (:,:,:,:,:,:,:)
! call sp1(x,frc_new,nr1,nr2,nr3,nat,scal)
! if (DABS(scal).gt.1d-10) write(6,'("k= ",I8," frc_new|u(k)= ",F15.10)') k,scal
! deallocate(u(k) % vec)
!enddo
!
frc(:,:,:,:,:,:,:)=frc_new(:,:,:,:,:,:,:)
!
deallocate (x, w)
deallocate (v, ind_v)
deallocate (frc_new)
!
return
end subroutine set_asr
!
!----------------------------------------------------------------------
subroutine sp_zeu(zeu_u,zeu_v,nat,scal)
!-----------------------------------------------------------------------
!
! does the scalar product of two effective charges matrices zeu_u and zeu_v
! (considered as vectors in the R^(3*3*nat) space, and coded in the usual way)
!
USE kinds, ONLY: DP
implicit none
integer i,j,na,nat
real(DP) zeu_u(3,3,nat)
real(DP) zeu_v(3,3,nat)
real(DP) scal
!
!
scal=0.0d0
do i=1,3
do j=1,3
do na=1,nat
scal=scal+zeu_u(i,j,na)*zeu_v(i,j,na)
enddo
enddo
enddo
!
return
!
end subroutine sp_zeu
!
!
!----------------------------------------------------------------------
subroutine sp1(u,v,nr1,nr2,nr3,nat,scal)
!-----------------------------------------------------------------------
!
! does the scalar product of two force-constants matrices u and v (considered as
! vectors in the R^(3*3*nat*nat*nr1*nr2*nr3) space, and coded in the usual way)
!
USE kinds, ONLY: DP
implicit none
integer nr1,nr2,nr3,i,j,na,nb,n1,n2,n3,nat
real(DP) u(nr1,nr2,nr3,3,3,nat,nat)
real(DP) v(nr1,nr2,nr3,3,3,nat,nat)
real(DP) scal
!
!
scal=0.0d0
do i=1,3
do j=1,3
do na=1,nat
do nb=1,nat
do n1=1,nr1
do n2=1,nr2
do n3=1,nr3
scal=scal+u(n1,n2,n3,i,j,na,nb)*v(n1,n2,n3,i,j,na,nb)
enddo
enddo
enddo
enddo
enddo
enddo
enddo
!
return
!
end subroutine sp1
!
!----------------------------------------------------------------------
subroutine sp2(u,v,ind_v,nr1,nr2,nr3,nat,scal)
!-----------------------------------------------------------------------
!
! does the scalar product of two force-constants matrices u and v (considered as
! vectors in the R^(3*3*nat*nat*nr1*nr2*nr3) space). u is coded in the usual way
! but v is coded as explained when defining the vectors corresponding to the
! symmetry constraints
!
USE kinds, ONLY: DP
implicit none
integer nr1,nr2,nr3,i,nat
real(DP) u(nr1,nr2,nr3,3,3,nat,nat)
integer ind_v(2,7)
real(DP) v(2)
real(DP) scal
!
!
scal=0.0d0
do i=1,2
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), &
ind_v(i,7))*v(i)
enddo
!
return
!
end subroutine sp2
!
!----------------------------------------------------------------------
subroutine sp3(u,v,i,na,nr1,nr2,nr3,nat,scal)
!-----------------------------------------------------------------------
!
! like sp1, but in the particular case when u is one of the u(k)%vec
! defined in set_asr (before orthonormalization). In this case most of the
! terms are zero (the ones that are not are characterized by i and na), so
! that a lot of computer time can be saved (during Gram-Schmidt).
!
USE kinds, ONLY: DP
implicit none
integer nr1,nr2,nr3,i,j,na,nb,n1,n2,n3,nat
real(DP) u(nr1,nr2,nr3,3,3,nat,nat)
real(DP) v(nr1,nr2,nr3,3,3,nat,nat)
real(DP) scal
!
!
scal=0.0d0
do j=1,3
do nb=1,nat
do n1=1,nr1
do n2=1,nr2
do n3=1,nr3
scal=scal+u(n1,n2,n3,i,j,na,nb)*v(n1,n2,n3,i,j,na,nb)
enddo
enddo
enddo
enddo
enddo
!
return
!
end subroutine sp3
!
!-----------------------------------------------------------------------
SUBROUTINE ws_all(tau,neq,nr1,nr2,nr3,r_ws,at_blk)
!-----------------------------------------------------------------------
!! Used in set_asr
!! Determine the nearest periodic image for an atom pair
!! Calculte the number of degeneracies as weight when
!! on the surface of the optimal Wigner-Seitz cell
!
USE kinds, ONLY: DP
!
IMPLICIT NONE
!
INTEGER, INTENT(in) :: nr1,nr2,nr3
!! Supercell size
INTEGER, INTENT(out) :: neq
!! Weight of the degeneracy
REAL(DP), INTENT(in) :: r_ws(3)
!! Original position in the supercell
REAL(DP), INTENT(in) :: at_blk(3, 3)
!! Lattice vector of the primitive cell
REAL(DP), INTENT(out) :: tau(3, 48)
!! Nearest periodic atomic position
!
! Local variable
INTEGER :: n1, n2, n3
!! Search around the cell to a maximum of 2 around
REAL(DP) :: rnorm
!! Distance betweem origin cell and current cell
REAL(DP) :: eps
!! Tolerence
REAL(DP) :: tau_tmp(3)
!! Position of the current cell
REAL(DP) :: at(3,3)
!! Lattice vector of the supercell
REAL(DP) :: rmin
!! Minimal distance
!
! Large tolerence to find degeneracy position in case of bad relaxation
eps = 1.0d-5
!
at(:, 1) = at_blk(:, 1) * DBLE(nr1)
at(:, 2) = at_blk(:, 2) * DBLE(nr2)
at(:, 3) = at_blk(:, 3) * DBLE(nr3)
!
rmin = HUGE(rmin)
DO n1 = -2,2
DO n2 = -2,2
DO n3 = -2,2
tau_tmp = r_ws + n1 * at(:, 1) + n2 * at(:, 2) + n3 * at(:, 3)
rnorm = NORM2(tau_tmp)
IF (ABS(rnorm - rmin) .gt. eps) THEN
IF (rnorm .lt. rmin) THEN
neq = 1
rmin = rnorm
tau(:, neq) = tau_tmp
ENDIF
ELSE
neq = neq + 1
tau(:, neq) = tau_tmp
ENDIF
ENDDO
ENDDO
ENDDO
!
RETURN
!
!-----------------------------------------------------------------------
END SUBROUTINE ws_all
!-----------------------------------------------------------------------
!
!-----------------------------------------------------------------------
SUBROUTINE q_gen(nsc,qbid,at_blk,bg_blk,at,bg)
!-----------------------------------------------------------------------
! generate list of q (qbid) that are G-vectors of the supercell
! but not of the bulk
!
USE kinds, ONLY : DP
!
IMPLICIT NONE
INTEGER :: nsc
REAL(DP) qbid(3,nsc), at_blk(3,3), bg_blk(3,3), at(3,3), bg(3,3)
!
INTEGER, PARAMETER:: nr1=4, nr2=4, nr3=4, &
nrm=(2*nr1+1)*(2*nr2+1)*(2*nr3+1)
REAL(DP), PARAMETER:: eps=1.0d-7
INTEGER :: i, j, k,i1, i2, i3, idum(nrm), iq
REAL(DP) :: qnorm(nrm), qbd(3,nrm) ,qwork(3), delta
LOGICAL lbho
!
i = 0
DO i1=-nr1,nr1
DO i2=-nr2,nr2
DO i3=-nr3,nr3
i = i + 1
DO j=1,3
qwork(j) = i1*bg(j,1) + i2*bg(j,2) + i3*bg(j,3)
END DO ! j
!
qnorm(i) = qwork(1)**2 + qwork(2)**2 + qwork(3)**2
!
DO j=1,3
!
qbd(j,i) = at_blk(1,j)*qwork(1) + &
at_blk(2,j)*qwork(2) + &
at_blk(3,j)*qwork(3)
END DO ! j
!
idum(i) = 1
!
END DO ! i3
END DO ! i2
END DO ! i1
!
DO i=1,nrm-1
IF (idum(i).EQ.1) THEN
DO j=i+1,nrm
IF (idum(j).EQ.1) THEN
lbho=.TRUE.
DO k=1,3
delta = qbd(k,i)-qbd(k,j)
lbho = lbho.AND. (ABS(NINT(delta)-delta).LT.eps)
END DO ! k
IF (lbho) THEN
IF(qnorm(i).GT.qnorm(j)) THEN
qbd(1,i) = qbd(1,j)
qbd(2,i) = qbd(2,j)
qbd(3,i) = qbd(3,j)
qnorm(i) = qnorm(j)
END IF
idum(j) = 0
END IF
END IF
END DO ! j
END IF
END DO ! i
!
iq = 0
DO i=1,nrm
IF (idum(i).EQ.1) THEN
iq=iq+1
qbid(1,iq)= bg_blk(1,1)*qbd(1,i) + &
bg_blk(1,2)*qbd(2,i) + &
bg_blk(1,3)*qbd(3,i)
qbid(2,iq)= bg_blk(2,1)*qbd(1,i) + &
bg_blk(2,2)*qbd(2,i) + &
bg_blk(2,3)*qbd(3,i)
qbid(3,iq)= bg_blk(3,1)*qbd(1,i) + &
bg_blk(3,2)*qbd(2,i) + &
bg_blk(3,3)*qbd(3,i)
END IF
END DO ! i
!
IF (iq.NE.nsc) CALL errore('q_gen',' probably nr1,nr2,nr3 too small ', iq)
RETURN
END SUBROUTINE q_gen
!
!-----------------------------------------------------------------------
SUBROUTINE check_at(at,bg_blk,alat,omega)
!-----------------------------------------------------------------------
!
USE kinds, ONLY : DP
USE io_global, ONLY : stdout
!
IMPLICIT NONE
!
REAL(DP) :: at(3,3), bg_blk(3,3), alat, omega
REAL(DP) :: work(3,3)
INTEGER :: i,j
REAL(DP), PARAMETER :: small=1.d-6
!
work(:,:) = at(:,:)
CALL cryst_to_cart(3,work,bg_blk,-1)
!
DO j=1,3
DO i =1,3
IF ( ABS(work(i,j)-NINT(work(i,j))) > small) THEN
WRITE (stdout,'(3f9.4)') work(:,:)
CALL errore ('check_at','at not multiple of at_blk',1)
END IF
END DO
END DO
!
omega =alat**3 * ABS(at(1,1)*(at(2,2)*at(3,3)-at(3,2)*at(2,3))- &
at(1,2)*(at(2,1)*at(3,3)-at(2,3)*at(3,1))+ &
at(1,3)*(at(2,1)*at(3,2)-at(2,2)*at(3,1)))
!
RETURN
END SUBROUTINE check_at
!
!-----------------------------------------------------------------------
SUBROUTINE set_tau (nat, nat_blk, at, at_blk, tau, tau_blk, &
ityp, ityp_blk, itau_blk)
!-----------------------------------------------------------------------
!
USE kinds, ONLY : DP
!
IMPLICIT NONE
INTEGER nat, nat_blk,ityp(nat),ityp_blk(nat_blk), itau_blk(nat)
REAL(DP) at(3,3),at_blk(3,3),tau(3,nat),tau_blk(3,nat_blk)
!
REAL(DP) bg(3,3), r(3) ! work vectors
INTEGER i,i1,i2,i3,na,na_blk
REAL(DP) small
INTEGER NN1,NN2,NN3
PARAMETER (NN1=8, NN2=8, NN3=8, small=1.d-8)
!
CALL recips (at(1,1),at(1,2),at(1,3),bg(1,1),bg(1,2),bg(1,3))
!
na = 0
!
DO i1 = -NN1,NN1
DO i2 = -NN2,NN2
DO i3 = -NN3,NN3
r(1) = i1*at_blk(1,1) + i2*at_blk(1,2) + i3*at_blk(1,3)
r(2) = i1*at_blk(2,1) + i2*at_blk(2,2) + i3*at_blk(2,3)
r(3) = i1*at_blk(3,1) + i2*at_blk(3,2) + i3*at_blk(3,3)
CALL cryst_to_cart(1,r,bg,-1)
!
IF ( r(1).GT.-small .AND. r(1).LT.1.d0-small .AND. &
r(2).GT.-small .AND. r(2).LT.1.d0-small .AND. &
r(3).GT.-small .AND. r(3).LT.1.d0-small ) THEN
CALL cryst_to_cart(1,r,at,+1)
!
DO na_blk=1, nat_blk
na = na + 1
IF (na.GT.nat) CALL errore('set_tau','too many atoms',na)
tau(1,na) = tau_blk(1,na_blk) + r(1)
tau(2,na) = tau_blk(2,na_blk) + r(2)
tau(3,na) = tau_blk(3,na_blk) + r(3)
ityp(na) = ityp_blk(na_blk)
itau_blk(na) = na_blk
END DO
!
END IF
!
END DO
END DO
END DO
!
IF (na.NE.nat) CALL errore('set_tau','too few atoms: increase NNs',na)
!
RETURN
END SUBROUTINE set_tau
!
!-----------------------------------------------------------------------
SUBROUTINE read_tau &
(nat, nat_blk, ntyp, bg_blk, tau, tau_blk, ityp, itau_blk)
!---------------------------------------------------------------------
!
USE kinds, ONLY : DP
USE io_global, ONLY : ionode_id, ionode
USE mp, ONLY : mp_bcast
USE mp_world, ONLY : world_comm
!
IMPLICIT NONE
!
INTEGER nat, nat_blk, ntyp, ityp(nat),itau_blk(nat)
REAL(DP) bg_blk(3,3),tau(3,nat),tau_blk(3,nat_blk)
!
REAL(DP) r(3) ! work vectors
INTEGER i,na,na_blk
!
REAL(DP) small
PARAMETER ( small = 1.d-6 )
!
DO na=1,nat
IF (ionode) READ(5,*) (tau(i,na),i=1,3), ityp(na)
CALL mp_bcast(tau(:,na),ionode_id, world_comm)
CALL mp_bcast(ityp(na),ionode_id, world_comm)
IF (ityp(na).LE.0 .OR. ityp(na) .GT. ntyp) &
CALL errore('read_tau',' wrong atomic type', na)
DO na_blk=1,nat_blk
r(1) = tau(1,na) - tau_blk(1,na_blk)
r(2) = tau(2,na) - tau_blk(2,na_blk)
r(3) = tau(3,na) - tau_blk(3,na_blk)
CALL cryst_to_cart(1,r,bg_blk,-1)
IF (ABS( r(1)-NINT(r(1)) ) .LT. small .AND. &
ABS( r(2)-NINT(r(2)) ) .LT. small .AND. &
ABS( r(3)-NINT(r(3)) ) .LT. small ) THEN
itau_blk(na) = na_blk
go to 999
END IF
END DO
CALL errore ('read_tau',' wrong atomic position ', na)
999 CONTINUE
END DO
!
RETURN
END SUBROUTINE read_tau
!
!-----------------------------------------------------------------------
SUBROUTINE write_tau(fltau,nat,tau,ityp)
!-----------------------------------------------------------------------
!
USE kinds, ONLY : DP
USE io_global, ONLY : ionode
!
IMPLICIT NONE
!
INTEGER nat, ityp(nat)
REAL(DP) tau(3,nat)
CHARACTER(LEN=*) fltau
!
INTEGER i,na
!
IF (.NOT.ionode) RETURN
OPEN (unit=4,file=fltau, status='new')
DO na=1,nat
WRITE(4,'(3(f12.6),i3)') (tau(i,na),i=1,3), ityp(na)
END DO
CLOSE (4)
!
RETURN
END SUBROUTINE write_tau
!
!-----------------------------------------------------------------------
SUBROUTINE gen_qpoints (ibrav, at_, bg_, nat, tau, ityp, nk1, nk2, nk3, &
nqx, nq, q, nosym, wk)
!-----------------------------------------------------------------------
!
USE kinds, ONLY : DP
USE cell_base, ONLY : at, bg
USE symm_base, ONLY : set_sym_bl, find_sym, s, irt, nsym, &
nrot, t_rev, time_reversal, sname
USE ktetra, ONLY : tetra_init
!
IMPLICIT NONE
! input
INTEGER :: ibrav, nat, nk1, nk2, nk3, ityp(*)
REAL(DP) :: at_(3,3), bg_(3,3), tau(3,nat)
LOGICAL :: nosym
! output
INTEGER :: nqx, nq
REAL(DP) :: q(3,nqx), wk(nqx)
! local
REAL(DP) :: xqq(3), mdum(3,nat)
LOGICAL :: magnetic_sym=.FALSE., skip_equivalence=.FALSE.
!
time_reversal = .true.
if (nosym) time_reversal = .false.
t_rev(:) = 0
xqq (:) =0.d0
at = at_
bg = bg_
CALL set_sym_bl ( )
!
if (nosym) then
nrot = 1
nsym = 1
endif
CALL kpoint_grid ( nrot, time_reversal, skip_equivalence, s, t_rev, bg, nqx, &
0,0,0, nk1,nk2,nk3, nq, q, wk)
!
CALL find_sym ( nat, tau, ityp, .not.time_reversal, mdum )
!
CALL irreducible_BZ (nrot, s, nsym, time_reversal, magnetic_sym, &
at, bg, nqx, nq, q, wk, t_rev)
!
CALL tetra_init (nsym, s, time_reversal, t_rev, at, bg, nqx, 0, 0, 0, &
nk1, nk2, nk3, nq, q)
!
RETURN
END SUBROUTINE gen_qpoints
!
!---------------------------------------------------------------------
SUBROUTINE a2Fdos &
(nat, nq, nr1, nr2, nr3, ibrav, at, bg, tau, alat, &
nsc, nat_blk, at_blk, bg_blk, itau_blk, omega_blk, rws, nrws, &
dos, Emin, DeltaE, ndos, asr, q, freq, fd, wq, huang)
!-----------------------------------------------------------------------
!
USE kinds, ONLY : DP
USE io_global, ONLY : ionode, ionode_id
USE mp, ONLY : mp_bcast
USE mp_world, ONLY : world_comm
USE mp_images, ONLY : intra_image_comm
USE ifconstants, ONLY : zeu, tau_blk, frc_lr
USE constants, ONLY : pi, RY_TO_THZ
USE constants, ONLY : K_BOLTZMANN_RY
USE el_phon, ONLY : el_ph_nsigma
!
IMPLICIT NONE
!
INTEGER, INTENT(in) :: nat, nq, nr1, nr2, nr3, ibrav, ndos
LOGICAL, INTENT(in) :: dos, fd, huang
CHARACTER(LEN=*), INTENT(IN) :: asr
REAL(DP), INTENT(in) :: freq(3*nat,nq), q(3,nq), wq(nq), at(3,3), bg(3,3), &
tau(3,nat), alat, Emin, DeltaE
!
INTEGER, INTENT(in) :: nsc, nat_blk, itau_blk(nat), nrws
REAL(DP), INTENT(in) :: rws(0:3,nrws), at_blk(3,3), bg_blk(3,3), omega_blk
!
REAL(DP), ALLOCATABLE :: gamma(:,:), frcg(:,:,:,:,:,:,:)
COMPLEX(DP), ALLOCATABLE :: gam(:,:,:,:), gam_blk(:,:,:,:), z(:,:)
real(DP) :: lambda, dos_a2F(3*nat), temp, dos_ee(el_ph_nsigma), dos_tot, &
deg(el_ph_nsigma), fermi(el_ph_nsigma), E
real(DP), parameter :: eps_w2 = 0.0000001d0
integer :: isig, ifn, n, m, na, nb, nc, nu, nmodes, &
i,j,k, ngauss, jsig, p1, p2, p3, filea2F, ios
character(len=256) :: name
real(DP), external :: dos_gam
CHARACTER(LEN=6) :: int_to_char
!
!
nmodes = 3*nat
do isig=1,el_ph_nsigma
filea2F = 60 + isig
name= 'elph_dir/a2Fmatdyn.' // TRIM(int_to_char(filea2F))
IF (ionode) open(unit=filea2F, file=TRIM(name), &
STATUS = 'unknown', IOSTAT=ios)
CALL mp_bcast(ios, ionode_id, intra_image_comm)
IF (ios /= 0) CALL errore('a2Fdos','problem opening file'//TRIM(name),1)
IF (ionode) &
READ(filea2F,*) deg(isig), fermi(isig), dos_ee(isig)
ENDDO
call mp_bcast(deg, ionode_id, intra_image_comm)
call mp_bcast(fermi, ionode_id, intra_image_comm)
call mp_bcast(dos_ee, ionode_id, intra_image_comm)
!
IF (ionode) THEN
IF(dos) then
open(unit=400,file='lambda',status='unknown')
write(400,*)
write(400,*) ' Electron-phonon coupling constant, lambda '
write(400,*)
ELSE
open (unit=20,file='gam.lines' ,status='unknown')
write(20,*)
write(20,*) ' Gamma lines for all modes [THz] '
write(20,*)
write(6,*)
write(6,*) ' Gamma lines for all modes [Rydberg] '
write(6,*)
ENDIF
ENDIF
!
ALLOCATE ( frcg(nr1,nr2,nr3,3,3,nat,nat) )
ALLOCATE ( gamma(3*nat,nq), gam(3,3,nat,nat), gam_blk(3,3,nat_blk,nat_blk) )
ALLOCATE ( z(3*nat,3*nat) )
!
frcg(:,:,:,:,:,:,:) = 0.d0
DO isig = 1, el_ph_nsigma
filea2F = 60 + isig
CALL readfg ( filea2F, nr1, nr2, nr3, nat, frcg )
!
if ( asr /= 'no') then
CALL set_asr (asr, nr1, nr2, nr3, frcg, frc_lr, zeu, nat_blk, ibrav, &
tau_blk, at_blk, huang)
endif
!
IF (ionode) open(unit=300,file='dyna2F',status='old')
!
do n = 1 ,nq
gam(:,:,:,:) = (0.d0, 0.d0)
IF (ionode) THEN
read(300,*)
do na=1,nmodes
read(300,*) (z(na,m),m=1,nmodes)
end do ! na
ENDIF
CALL mp_bcast(z, ionode_id, world_comm)
!
CALL setgam (q(1,n), gam, nat, at, bg, tau, itau_blk, nsc, alat, &
gam_blk, nat_blk, at_blk,bg_blk,tau_blk, omega_blk, &
frcg, nr1,nr2,nr3, rws, nrws, fd)
!
! here multiply dyn*gam*dyn for gamma and divide by w2 for lambda at given q
!
do nc = 1, nat
do k =1, 3
p1 = (nc-1)*3+k
nu = p1
gamma(nu,n) = 0.0d0
do i=1,3
do na=1,nat
p2 = (na-1)*3+i
do j=1,3
do nb=1,nat
p3 = (nb-1)*3+j
gamma(nu,n) = gamma(nu,n) + DBLE(conjg(z(p2,p1)) * &
gam(i,j,na,nb) * z(p3,p1))
enddo ! nb
enddo ! j
enddo ! na
enddo !i
gamma(nu,n) = gamma(nu,n) * pi / 2.0d0
enddo ! k
enddo !nc
!
!
EndDo !nq all points in BZ
IF (ionode) close(300) ! file with dyn vectors
!
! after we know gamma(q) and lambda(q) calculate DOS(omega) for spectrum a2F
!
if(dos.and.ionode) then
!
name='a2F.dos'//int_to_char(isig)
ifn = 200 + isig
open (ifn,file=TRIM(name),status='unknown',form='formatted')
write(ifn,*)
write(ifn,*) '# Eliashberg function a2F (per both spin)'
write(ifn,*) '# frequencies in Rydberg '
write(ifn,*) '# DOS normalized to E in Rydberg: a2F_total, a2F(mode) '
write(ifn,*)
!
! correction for small frequencies
!
do n = 1, nq
do i = 1, nmodes
if (freq(i,n).LE.eps_w2) then
gamma(i,n) = 0.0d0
endif
enddo
enddo
!
lambda = 0.0d0
do n= 1, ndos
!
E = Emin + (n-1)*DeltaE + 0.5d0*DeltaE
dos_tot = 0.0d0
do j=1,nmodes
!
dos_a2F(j) = dos_gam(nmodes, nq, j, gamma, freq, E)
dos_a2F(j) = dos_a2F(j) / dos_ee(isig) / 2.d0 / pi
dos_tot = dos_tot + dos_a2F(j)
!
enddo
lambda = lambda + 2.d0 * dos_tot/E * DeltaE
write (ifn, '(3X,1000E16.6)') E, dos_tot, dos_a2F(1:nmodes)
enddo !ndos
write(ifn,*) " lambda =",lambda,' Delta = ',DeltaE
close (ifn)
!
! lambda from alternative way, simple sum.
! Also Omega_ln is computed
!
lambda = 0.0_dp
E = 0.0_dp
do n = 1, nq
lambda = lambda &
& + sum(gamma(1:nmodes,n)/freq(1:nmodes,n)**2, &
& freq(1:nmodes,n) > 1.0e-5_dp) * wq(n)
E = E &
& + sum(log(freq(1:nmodes,n)) * gamma(1:nmodes,n)/freq(1:nmodes,n)**2, &
& freq(1:nmodes,n) > 1.0e-5_dp) * wq(n)
end do
E = exp(E / lambda) / K_BOLTZMANN_RY
lambda = lambda / (dos_ee(isig) * pi)
write(400,'(" Broadening ",F8.4," lambda ",F12.4," dos(Ef)",F8.4," omega_ln [K]",F12.4)') &
deg(isig),lambda, dos_ee(isig), E
!
endif !dos
!
! OUTPUT
!
if(.not.dos.and.ionode) then
write(20,'(" Broadening ",F8.4)') deg(isig)
write( 6,'(" Broadening ",F8.4)') deg(isig)
do n=1, nq
write(20,'(3x,i5)') n
write( 6,'(3x,i5)') n
write(20,'(9F8.4)') (gamma(i,n)*RY_TO_THZ,i=1,3*nat)
write( 6,'(6F12.9)') (gamma(i,n),i=1,3*nat)
!
! write also in a format that can be read by plotband
!
WRITE(200+isig, '(10x,3f10.6)') q(1,n), q(2,n), q(3,n)
!
! output in GHz
!
WRITE(200+isig, '(6f10.4)') (gamma(nu,n)*RY_TO_THZ*1000.0_DP, &
nu=1,3*nat)
end do
endif
!
ENDDO !isig
!
DEALLOCATE (z, frcg, gamma, gam, gam_blk )
!
IF (ionode) THEN
close(400) !lambda
close(20)
ENDIF
!
!
RETURN
END SUBROUTINE a2Fdos
!
!-----------------------------------------------------------------------
subroutine setgam (q, gam, nat, at,bg,tau,itau_blk,nsc,alat, &
& gam_blk, nat_blk, at_blk,bg_blk,tau_blk,omega_blk, &
& frcg, nr1,nr2,nr3, rws,nrws, fd)
!-----------------------------------------------------------------------
! compute the dynamical matrix (the analytic part only)
!
USE kinds, ONLY : DP
USE constants, ONLY : tpi
implicit none
!
! I/O variables
!
integer :: nr1, nr2, nr3, nat, nat_blk, &
nsc, nrws, itau_blk(nat)
real(DP) :: q(3), tau(3,nat), at(3,3), bg(3,3), alat, rws(0:3,nrws)
real(DP) :: tau_blk(3,nat_blk), at_blk(3,3), bg_blk(3,3), omega_blk, &
frcg(nr1,nr2,nr3,3,3,nat_blk,nat_blk)
COMPLEX(DP) :: gam_blk(3,3,nat_blk,nat_blk),f_of_q(3,3,nat,nat)
COMPLEX(DP) :: gam(3,3,nat,nat)
LOGICAL :: fd
!
! local variables
!
real(DP) :: arg
complex(DP) :: cfac(nat)
integer :: i,j,k, na,nb, na_blk, nb_blk, iq
real(DP) :: qp(3), qbid(3,nsc) ! automatic array
!
!
call q_gen(nsc,qbid,at_blk,bg_blk,at,bg)
!
f_of_q=(0.0_DP,0.0_DP)
do iq=1,nsc
!
do k=1,3
qp(k)= q(k) + qbid(k,iq)
end do
!
gam_blk(:,:,:,:) = (0.d0,0.d0)
CALL frc_blk (gam_blk,qp,tau_blk,nat_blk, &
nr1,nr2,nr3,frcg,at_blk,bg_blk,rws,nrws,f_of_q,fd)
!
do na=1,nat
na_blk = itau_blk(na)
do nb=1,nat
nb_blk = itau_blk(nb)
!
arg = tpi * ( qp(1) * ( (tau(1,na)-tau_blk(1,na_blk)) - &
(tau(1,nb)-tau_blk(1,nb_blk)) ) + &
qp(2) * ( (tau(2,na)-tau_blk(2,na_blk)) - &
(tau(2,nb)-tau_blk(2,nb_blk)) ) + &
qp(3) * ( (tau(3,na)-tau_blk(3,na_blk)) - &
(tau(3,nb)-tau_blk(3,nb_blk)) ) )
!
cfac(nb) = CMPLX(cos(arg),sin(arg), kind=dp)/nsc
!
end do ! nb
do nb=1,nat
do i=1,3
do j=1,3
nb_blk = itau_blk(nb)
gam(i,j,na,nb) = gam(i,j,na,nb) + cfac(nb) * &
gam_blk(i,j,na_blk,nb_blk)
end do ! j
end do ! i
end do ! nb
end do ! na
!
end do ! iq
!
return
end subroutine setgam
!
!--------------------------------------------------------------------
function dos_gam (nbndx, nq, jbnd, gamma, et, ef)
!--------------------------------------------------------------------
! calculates weights with the tetrahedron method (Bloechl version)
! this subroutine is based on tweights.f90 belonging to PW
! it calculates a2F on the surface of given frequency <=> histogram
! Band index means the frequency mode here
! and "et" means the frequency(mode,q-point)
!
USE kinds, ONLY: DP
USE parameters
USE ktetra, ONLY : ntetra, tetra
implicit none
!
integer :: nq, nbndx, jbnd
real(DP) :: et(nbndx,nq), gamma(nbndx,nq), func
real(DP) :: ef
real(DP) :: e1, e2, e3, e4, c1, c2, c3, c4, etetra(4)
integer :: ik, ibnd, nt, nk, ns, i, ik1, ik2, ik3, ik4, itetra(4)
real(DP) :: f12,f13,f14,f23,f24,f34, f21,f31,f41,f42,f32,f43
real(DP) :: P1,P2,P3,P4, G, o13, Y1,Y2,Y3,Y4, eps,vol, Tint
real(DP) :: dos_gam
Tint = 0.0d0
o13 = 1.0_dp/3.0_dp
eps = 1.0d-14
vol = 1.0d0/ntetra
P1 = 0.0_dp
P2 = 0.0_dp
P3 = 0.0_dp
P4 = 0.0_dp
do nt = 1, ntetra
ibnd = jbnd
!
! etetra are the energies at the vertexes of the nt-th tetrahedron
!
do i = 1, 4
etetra(i) = et(ibnd, tetra(i,nt))
enddo
itetra(1) = 0
call hpsort (4,etetra,itetra)
!
! ...sort in ascending order: e1 < e2 < e3 < e4
!
e1 = etetra (1)
e2 = etetra (2)
e3 = etetra (3)
e4 = etetra (4)
!
! kp1-kp4 are the irreducible k-points corresponding to e1-e4
!
ik1 = tetra(itetra(1),nt)
ik2 = tetra(itetra(2),nt)
ik3 = tetra(itetra(3),nt)
ik4 = tetra(itetra(4),nt)
Y1 = gamma(ibnd,ik1)/et(ibnd,ik1)
Y2 = gamma(ibnd,ik2)/et(ibnd,ik2)
Y3 = gamma(ibnd,ik3)/et(ibnd,ik3)
Y4 = gamma(ibnd,ik4)/et(ibnd,ik4)
IF ( e3 < ef .and. ef < e4) THEN
f14 = (ef-e4)/(e1-e4)
f24 = (ef-e4)/(e2-e4)
f34 = (ef-e4)/(e3-e4)
G = 3.0_dp * f14 * f24 * f34 / (e4-ef)
P1 = f14 * o13
P2 = f24 * o13
P3 = f34 * o13
P4 = (3.0_dp - f14 - f24 - f34 ) * o13
ELSE IF ( e2 < ef .and. ef < e3 ) THEN
f13 = (ef-e3)/(e1-e3)
f31 = 1.0_dp - f13
f14 = (ef-e4)/(e1-e4)
f41 = 1.0_dp-f14
f23 = (ef-e3)/(e2-e3)
f32 = 1.0_dp - f23
f24 = (ef-e4)/(e2-e4)
f42 = 1.0_dp - f24
G = 3.0_dp * (f23*f31 + f32*f24)
P1 = f14 * o13 + f13*f31*f23 / G
P2 = f23 * o13 + f24*f24*f32 / G
P3 = f32 * o13 + f31*f31*f23 / G
P4 = f41 * o13 + f42*f24*f32 / G
G = G / (e4-e1)
ELSE IF ( e1 < ef .and. ef < e2 ) THEN
f12 = (ef-e2)/(e1-e2)
f21 = 1.0_dp - f12
f13 = (ef-e3)/(e1-e3)
f31 = 1.0_dp - f13
f14 = (ef-e4)/(e1-e4)
f41 = 1.0_dp - f14
G = 3.0_dp * f21 * f31 * f41 / (ef-e1)
P1 = o13 * (f12 + f13 + f14)
P2 = o13 * f21
P3 = o13 * f31
P4 = o13 * f41
ELSE
G = 0.0_dp
END IF
Tint = Tint + G * (Y1*P1 + Y2*P2 + Y3*P3 + Y4*P4) * vol
enddo ! ntetra
dos_gam = Tint !2 because DOS_ee is per 1 spin
return
end function dos_gam
!
!
!------------------------------------------------------------------------------
function dos_broad(iatom, nbndx, nq, et, zq, wq, Ef, degauss)
!------------------------------------------------------------------------------
!
! Get atom projected ph DOS using Gaussian broadening. Inspired by
! thermo_pw/src/generalized_phdos.f90 (GPL)
!
USE kinds, ONLY : DP
!
IMPLICIT NONE
INTEGER, INTENT(IN) :: iatom, nq, nbndx
REAL(DP), INTENT(IN) :: et(nbndx, nq), zq(nq, nbndx, nbndx), wq(nq), Ef, degauss
REAL(DP) :: ufact, weight
INTEGER :: n, ik, ngauss, ipol, indi, jpol, indj
COMPLEX(DP) :: u1, u2
REAL(DP), EXTERNAL :: w0gauss
REAL(DP) :: dos_broad
!
dos_broad = 0.0_DP
!
ngauss = 0
!
DO ik = 1, nq
DO n = 1, nbndx
weight = w0gauss ( (Ef - et(n,ik) ) / degauss, ngauss)
!
DO ipol=1, 3
indi = 3 * (iatom - 1) + ipol
DO jpol = ipol, 3
indj = 3 * (iatom - 1) + jpol
u1=zq(ik, indi, n)
u2=zq(ik, indj, n)
ufact = DBLE(u1*CONJG(u2))
dos_broad = dos_broad + wq(ik) * ufact * weight
ENDDO
ENDDO
!
ENDDO
ENDDO
!
dos_broad = dos_broad / degauss
!
RETURN
end function dos_broad
!
!
!-----------------------------------------------------------------------
subroutine readfg ( ifn, nr1, nr2, nr3, nat, frcg )
!-----------------------------------------------------------------------
!
USE kinds, ONLY : DP
USE io_global, ONLY : ionode, ionode_id, stdout
USE mp, ONLY : mp_bcast
USE mp_world, ONLY : world_comm
implicit none
! I/O variable
integer, intent(in) :: nr1,nr2,nr3, nat
real(DP), intent(out) :: frcg(nr1,nr2,nr3,3,3,nat,nat)
! local variables
integer i, j, na, nb, m1,m2,m3, ifn
integer ibid, jbid, nabid, nbbid, m1bid,m2bid,m3bid
!
!
IF (ionode) READ (ifn,*) m1, m2, m3
CALL mp_bcast(m1, ionode_id, world_comm)
CALL mp_bcast(m2, ionode_id, world_comm)
CALL mp_bcast(m3, ionode_id, world_comm)
if ( m1 /= nr1 .or. m2 /= nr2 .or. m3 /= nr3) &
call errore('readfg','inconsistent nr1, nr2, nr3 read',1)
do i=1,3
do j=1,3
do na=1,nat
do nb=1,nat
IF (ionode) read (ifn,*) ibid, jbid, nabid, nbbid
CALL mp_bcast(ibid, ionode_id, world_comm)
CALL mp_bcast(jbid, ionode_id, world_comm)
CALL mp_bcast(nabid, ionode_id, world_comm)
CALL mp_bcast(nbbid, ionode_id, world_comm)
if(i.ne.ibid.or.j.ne.jbid.or.na.ne.nabid.or.nb.ne.nbbid) then
write(stdout,*) i,j,na,nb,' <> ', ibid, jbid, nabid, nbbid
call errore ('readfG','error in reading',1)
else
IF (ionode) read (ifn,*) (((m1bid, m2bid, m3bid, &
frcg(m1,m2,m3,i,j,na,nb), &
m1=1,nr1),m2=1,nr2),m3=1,nr3)
endif
CALL mp_bcast(frcg(:,:,:,i,j,na,nb), ionode_id, world_comm)
end do
end do
end do
end do
!
IF (ionode) close(ifn)
!
return
end subroutine readfg
!
!
SUBROUTINE find_representations_mode_q ( nat, ntyp, xq, w2, u, tau, ityp, &
amass, num_rap_mode, nspin_mag )
USE kinds, ONLY : DP
USE cell_base, ONLY : at, bg
USE symm_base, ONLY : s, sr, ft, irt, nsym, nrot, t_rev, time_reversal,&
sname, copy_sym, s_axis_to_cart
IMPLICIT NONE
INTEGER, INTENT(IN) :: nat, ntyp, nspin_mag
REAL(DP), INTENT(IN) :: xq(3), amass(ntyp), tau(3,nat)
REAL(DP), INTENT(IN) :: w2(3*nat)
INTEGER, INTENT(IN) :: ityp(nat)
COMPLEX(DP), INTENT(IN) :: u(3*nat,3*nat)
INTEGER, INTENT(OUT) :: num_rap_mode(3*nat)
REAL(DP) :: gi (3, 48), gimq (3), sr_is(3,3,48), rtau(3,48,nat)
INTEGER :: irotmq, nsymq, nsym_is, isym, i, ierr
LOGICAL :: minus_q, search_sym, sym(48), magnetic_sym
!
! find the small group of q
!
time_reversal=.TRUE.
IF (.NOT.time_reversal) minus_q=.FALSE.
sym(1:nsym)=.true.
call smallg_q (xq, 0, at, bg, nsym, s, sym, minus_q)
nsymq=copy_sym(nsym,sym )
call s_axis_to_cart ()
CALL set_giq (xq,s,nsymq,nsym,irotmq,minus_q,gi,gimq)
!
! if the small group of q is non symmorphic,
! search the symmetries only if there are no G such that Sq -> q+G
!
search_sym=.TRUE.
IF ( ANY ( ABS(ft(:,1:nsymq)) > 1.0d-8 ) ) THEN
DO isym=1,nsymq
search_sym=( search_sym.and.(abs(gi(1,isym))<1.d-8).and. &
(abs(gi(2,isym))<1.d-8).and. &
(abs(gi(3,isym))<1.d-8) )
END DO
END IF
!
! Set the representations tables of the small group of q and
! find the mode symmetry
!
IF (search_sym) THEN
magnetic_sym=(nspin_mag==4)
CALL prepare_sym_analysis(nsymq,sr,t_rev,magnetic_sym)
sym (1:nsym) = .TRUE.
CALL sgam_lr (at, bg, nsym, s, irt, tau, rtau, nat)
CALL find_mode_sym_new (u, w2, tau, nat, nsymq, s, sr, irt, xq, &
rtau, amass, ntyp, ityp, 1, .FALSE., .FALSE., num_rap_mode, ierr)
ENDIF
RETURN
END SUBROUTINE find_representations_mode_q
Reference in New Issue View Git Blame Copy Permalink