When checkallsym finds that the symmetry is lower than the original one

the code symmetrizes the atomic (and iif needed the cell) configuration
before stopping, so as to allow to restart from a symmetric configuration
if desired.
SdG


git-svn-id: http://qeforge.qe-forge.org/svn/q-e/trunk/espresso@3425 c92efa57-630b-4861-b058-cf58834340f0
This commit is contained in:
degironc 2006-10-02 09:37:35 +00:00
parent 1148ec8480
commit c11b1c6adc
7 changed files with 412 additions and 245 deletions

View File

@ -221,6 +221,7 @@ sumkg.o \
sumkt.o \
summary.o \
swap.o \
symmetrize_at.o \
symrho.o \
symrho_mag.o \
symscalar.o \

View File

@ -7,7 +7,7 @@
!
!-----------------------------------------------------------------------
subroutine checkallsym (nsym, s, nat, tau, ityp, at, bg, nr1, nr2, &
nr3, irt, ftau)
nr3, irt, ftau, alat, omega)
!-----------------------------------------------------------------------
! given a crystal group this routine checks that the actual
! atomic positions and bravais lattice vectors are compatible with
@ -35,7 +35,7 @@ subroutine checkallsym (nsym, s, nat, tau, ityp, at, bg, nr1, nr2, &
! k-th symmetry operation, referred to the
! crystal axis and in units at/nr (0-nr-1)
real(DP) :: tau (3, nat), at (3, 3), bg (3, 3)
real(DP) :: tau (3, nat), at (3, 3), bg (3, 3), alat, omega
! input: cartesian coordinates of the atoms
! input: basis of the real-space lattice
! input: " " " reciprocal-space lattic
@ -126,14 +126,25 @@ subroutine checkallsym (nsym, s, nat, tau, ityp, at, bg, nr1, nr2, &
ft (3) = ftau (3, isym) / DBLE (nr3)
call checksym (isym, nat, ityp, xau, rau, ft, loksym, irt)
if (.not.loksym (isym) ) call errore ('checkallsym', &
' symmetry operation not satisfied ', isym)
enddo
!
! deallocate work space
!
deallocate(rau)
deallocate(xau)
!
!
do isym = 1,nsym
if (.not.loksym (isym) ) call errore ('checkallsym', &
'the following symmetry operation is not satisfied ', -isym)
end do
if (ANY (.not.loksym (1:nsym) ) ) then
call symmetrize_at(nsym, s, nat, tau, ityp, at, bg, nr1, nr2, &
nr3, irt, ftau, alat, omega)
call errore ('checkallsym', &
'some of the original symmetry operations not satisfied ',1)
end if
!
return
end subroutine checkallsym

File diff suppressed because it is too large Load Diff

View File

@ -240,7 +240,7 @@ SUBROUTINE move_ions()
! ... according to the symmetry of the system.
!
CALL checkallsym( nsym, s, nat, tau, ityp, &
at, bg, nr1, nr2, nr3, irt, ftau )
at, bg, nr1, nr2, nr3, irt, ftau, alat, omega )
!
END IF
!

View File

@ -144,7 +144,7 @@ SUBROUTINE read_file()
!
IF (nat>0) &
CALL checkallsym( nsym, s, nat, tau, &
ityp, at, bg, nr1, nr2, nr3, irt, ftau )
ityp, at, bg, nr1, nr2, nr3, irt, ftau, alat, omega )
!
! ... read pseudopotentials
!

View File

@ -44,7 +44,8 @@ SUBROUTINE setup()
USE io_global, ONLY : stdout
USE io_files, ONLY : tmp_dir, prefix, delete_if_present
USE constants, ONLY : pi, degspin
USE cell_base, ONLY : at, bg, alat, tpiba, tpiba2, ibrav, symm_type
USE cell_base, ONLY : at, bg, alat, tpiba, tpiba2, ibrav, symm_type,&
omega
USE ions_base, ONLY : nat, tau, ntyp => nsp, ityp, zv
USE basis, ONLY : startingpot, natomwfc
USE gvect, ONLY : gcutm, ecutwfc, dual, nr1, nr2, nr3
@ -576,7 +577,7 @@ SUBROUTINE setup()
!
IF (nat>0) &
CALL checkallsym( nsym, s, nat, tau, ityp, at, &
bg, nr1, nr2, nr3, irt, ftau )
bg, nr1, nr2, nr3, irt, ftau, alat, omega )
!
! ... if dynamics is done the system should have no symmetries
! ... (inversion symmetry alone is allowed)

150
PW/symmetrize_at.f90 Normal file
View File

@ -0,0 +1,150 @@
!
! Copyright (C) 2001 PWSCF 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 .
!
!
!-----------------------------------------------------------------------
subroutine symmetrize_at(nsym, s, nat, tau, ityp, at, bg, &
nr1, nr2, nr3, irt, ftau, alat, omega)
!-----------------------------------------------------------------------
!
! given a point group, this routine finds the subgroup which is
! the point group of the crystal under consideration
! non symmorphic groups non allowed, provided that fractional
! translations are commensurate with the FFT grid
!
! It sets the array sym, which for each operation of the original
! point group is true if this operation is also an operation of the
! total point group
!
#include "f_defs.h"
USE io_global, ONLY : stdout
USE cellmd, ONLY: at_old, lmovecell
USE kinds
implicit none
!
! input variables
!
integer :: nsym, s (3, 3, 48), nat, ityp (nat), nr1, nr2, nr3
real(DP) :: tau (3, nat), at (3, 3), bg (3, 3), alat, omega
! nsym : number of symmetry operation of the crystal
! s : symmetry operations of parent group
! nat : number of atoms in the unit cell
! ityp : species of each atom in the unit cell
! nr* : dimensions of the FFT mesh
! tau : cartesian coordinates of the atoms
! at : basis of the real-space lattice
! bg : " " " reciprocal-space lattice
!
! output variables
!
integer :: irt (48, nat), ftau (3, 48)
! irt(isym,na) : sym.op. isym sends atom na into atom irt(isym,na)
! ftau(:,isym) : fractional translation associated to sym.op. isym
! (in FFT coordinates: crystal axis, multiplied by nr*)
! sym(isym) : flag indicating if sym.op. isym in the parent group
! is a true symmetry operation of the crystal
!
! local variables
!
integer :: na, icar, ipol, jpol, kpol, lpol, nb, irot, i, j
! counters
real(DP) , allocatable :: xau (:,:)
! atomic coordinates in crystal axis
logical :: fractional_translations
real(DP) :: work(3), obnr(3), bg_old(3,3), sat(3,3), wrk(3,3), ba(3,3)
integer :: table(48,48), invs(3,3,48)
!
allocate(xau(3,nat))
!
! Compute the coordinates of each atom in the basis of
! the direct lattice vectors
!
xau = tau
tau = 0.d0
call cryst_to_cart( nat, xau, bg, -1 )
!
obnr(1) = 1.d0/dfloat(nr1)
obnr(2) = 1.d0/dfloat(nr2)
obnr(3) = 1.d0/dfloat(nr3)
do irot = 1, nsym
do na = 1, nat
do kpol = 1, 3
work (kpol) = &
s (1, kpol, irot) * xau (1, na) + &
s (2, kpol, irot) * xau (2, na) + &
s (3, kpol, irot) * xau (3, na) - &
ftau(kpol,irot)* obnr(kpol)
tau (kpol, irt(irot,na)) = tau (kpol, irt(irot,na)) + work(ipol) &
- nint(work(kpol)-xau(kpol,irt(irot,na)))
enddo
enddo
enddo
tau (:,:) = tau(:,:)/nsym
!
! If the cell is moving then the lattice vectors has to be
! symmetrized as well
!
if (lmovecell) then
CALL multable (nsym, s, table)
CALL inverse_s (nsym, s, table, invs)
CALL recips( at_old(1,1), at_old(1,2), at_old(1,3), &
bg_old(1,1), bg_old(1,2), bg_old(1,3) )
ba(ipol,jpol) = bg_old(1,ipol) * at(1,jpol) + &
bg_old(2,ipol) * at(2,jpol) + &
bg_old(3,ipol) * at(3,jpol)
at = 0.d0
!
! at(i) = 1/nsym sum_S at_old(m) S(l,m) <bg_old(l)|at(k)> invS(i,k)
!
do irot=1,nsym
do icar = 1, 3
do lpol = 1, 3
sat(icar,lpol) = at_old(icar,1) * s(lpol,1,irot) &
+ at_old(icar,2) * s(lpol,2,irot) &
+ at_old(icar,3) * s(lpol,3,irot)
end do
end do
do icar = 1, 3
do kpol =1, 3
wrk(icar,kpol) = sat(icar,1) * ba(1,kpol) &
+ sat(icar,2) * ba(2,kpol) &
+ sat(icar,3) * ba(3,kpol)
end do
end do
do icar = 1, 3
do ipol =1, 3
at(icar,ipol) = at(icar,ipol) &
+ wrk(icar,1) * invs(ipol,1,irot) &
+ wrk(icar,2) * invs(ipol,2,irot) &
+ wrk(icar,3) * invs(ipol,3,irot)
end do
end do
end do
at(:,:) = at(:,:) / nsym
CALL volume( alat, at(1,1), at(1,2), at(1,3), omega )
CALL recips( at(1,1), at(1,2), at(1,3), bg(1,1), bg(1,2), bg(1,3) )
end if
!
! deallocate work space
!
deallocate (xau)
!
!
!
call cryst_to_cart(nat, tau, at, 1)
write (stdout,*) " SYMMETRIZED ATOMIC COORDINATES) "
call output_tau(lmovecell)
!
return
end subroutine symmetrize_at