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

Ford-PW part 6

This commit is contained in:
fabrizio22 2019-08-26 16:35:24 +02:00 committed by Samuel Ponce
parent d7a1dcfb4e
commit be857a6690
9 changed files with 11852 additions and 10036 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

30
PW/src/eqvect.f90
View File

@ -6,19 +6,27 @@
! or http://www.gnu.org/copyleft/gpl.txt .
!
!-----------------------------------------------------------------------
logical function eqvect (x, y, f, accep )
LOGICAL FUNCTION eqvect( x, y, f, accep )
!-----------------------------------------------------------------------
!
! This function test if the difference x-y-f is an integer.
! x, y = 3d vectors in crystal axis, f = fractionary translation
!! This function tests if the difference x-y-f is an integer.
!
USE kinds
implicit none
real(DP), intent(in) :: x (3), y (3), f (3), accep
!
eqvect = abs( x(1)-y(1)-f(1) - nint(x(1)-y(1)-f(1)) ) < accep .and. &
abs( x(2)-y(2)-f(2) - nint(x(2)-y(2)-f(2)) ) < accep .and. &
abs( x(3)-y(3)-f(3) - nint(x(3)-y(3)-f(3)) ) < accep
IMPLICIT NONE
!
return
end function eqvect
REAL(DP), INTENT(IN) :: x(3)
!! first 3d vector in crystal axis
REAL(DP), INTENT(IN) :: y(3)
!! second 3d vector in crystal axis
REAL(DP), INTENT(IN) :: f(3)
!! fractionary translation
REAL(DP), INTENT(IN) :: accep
!! threshold of acceptability
!
eqvect = ABS( x(1)-y(1)-f(1) - NINT(x(1)-y(1)-f(1)) ) < accep .AND. &
ABS( x(2)-y(2)-f(2) - NINT(x(2)-y(2)-f(2)) ) < accep .AND. &
ABS( x(3)-y(3)-f(3) - NINT(x(3)-y(3)-f(3)) ) < accep
!
RETURN
!
END FUNCTION eqvect

15472
PW/src/esm.f90
View File

File diff suppressed because it is too large Load Diff

192
PW/src/ewald.f90
View File

@ -7,65 +7,72 @@
!
!
!-----------------------------------------------------------------------
function ewald (alat, nat, ntyp, ityp, zv, at, bg, tau, omega, g, &
gg, ngm, gcutm, gstart, gamma_only, strf)
FUNCTION ewald( alat, nat, ntyp, ityp, zv, at, bg, tau, omega, g, &
gg, ngm, gcutm, gstart, gamma_only, strf )
!-----------------------------------------------------------------------
!
! Calculates Ewald energy with both G- and R-space terms.
! Determines optimal alpha. Should hopefully work for any structure.
!
!! Calculates Ewald energy with both G- and R-space terms.
!! Determines optimal alpha. Should hopefully work for any structure.
!
USE kinds
USE constants, ONLY : tpi, e2
USE mp_bands, ONLY : intra_bgrp_comm, me_bgrp, nproc_bgrp
USE mp, ONLY : mp_sum
USE constants, ONLY : tpi, e2
USE mp_bands, ONLY : intra_bgrp_comm, me_bgrp, nproc_bgrp
USE mp, ONLY : mp_sum
USE martyna_tuckerman, ONLY : wg_corr_ewald, do_comp_mt
USE Coul_cut_2D, ONLY : do_cutoff_2D, cutoff_ewald
implicit none
USE Coul_cut_2D, ONLY : do_cutoff_2D, cutoff_ewald
!
! first the dummy variables
IMPLICIT NONE
!
integer, intent(in) :: nat, ntyp, ityp (nat), ngm, gstart
! input: number of atoms in the unit cell
! input: number of different types of atoms
! input: the type of each atom
! input: number of plane waves for G sum
! input: first non-zero G vector
logical, intent(in) :: gamma_only
real(DP), intent(in) :: tau (3, nat), g (3, ngm), gg (ngm), zv (ntyp), &
at (3, 3), bg (3, 3), omega, alat, gcutm
! input: the positions of the atoms in the cell
! input: the coordinates of G vectors
! input: the square moduli of G vectors
! input: the charge of each type of atoms
! input: the direct lattice vectors
! input: the reciprocal lattice vectors
! input: the volume of the unit cell
! input: lattice parameter
! input: cut-off of g vectors
complex(DP), intent(in) :: strf (ngm, ntyp)
! input: structure factor
real(DP) :: ewald
! output: the ewald energy
INTEGER, INTENT(IN) :: nat
!! number of atoms in the unit cell
INTEGER, INTENT(IN) :: ntyp
!! number of different types of atoms
INTEGER, INTENT(IN) :: ityp(nat)
!! the type of each atom
INTEGER, INTENT(IN) :: ngm
!! number of plane waves for G sum
INTEGER, INTENT(IN) :: gstart
!! first non-zero G vector
LOGICAL, INTENT(IN) :: gamma_only
!! gamma only
REAL(DP), INTENT(IN) :: tau(3,nat)
!! the positions of the atoms in the cell
REAL(DP), INTENT(IN) :: g(3,ngm)
!! the coordinates of G vectors
REAL(DP), INTENT(IN) :: gg(ngm)
!! the square moduli of G vectors
REAL(DP), INTENT(IN) :: zv(ntyp)
!! the charge of each type of atoms
REAL(DP), INTENT(IN) :: at(3,3)
!! the direct lattice vectors
REAL(DP), INTENT(IN) :: bg(3,3)
!! the reciprocal lattice vectors
REAL(DP), INTENT(IN) :: omega
!! the volume of the unit cell
REAL(DP), INTENT(IN) :: alat
!! lattice parameter
REAL(DP), INTENT(IN) :: gcutm
!! cut-off of g vectors
COMPLEX(DP), INTENT(IN) :: strf(ngm,ntyp)
!! structure factor
REAL(DP) :: ewald
!! output: the ewald energy
!
! here the local variables
! ... local variables
!
integer, parameter :: mxr = 50
INTEGER, PARAMETER :: mxr = 50
! the maximum number of R vectors included in r
integer :: ng, nr, na, nb, nt, nrm
INTEGER :: ng, nr, na, nb, nt, nrm
! counter over reciprocal G vectors
! counter over direct vectors
! counter on atoms
! counter on atoms
! counter on atomic types
! number of R vectors included in r sum
integer :: na_s, na_e, mykey
INTEGER :: na_s, na_e, mykey
! for parallelization of real-space sums
real(DP) :: charge, tpiba2, ewaldg, ewaldr, dtau (3), alpha, &
r (3, mxr), r2 (mxr), rmax, rr, upperbound, fact
!
REAL(DP) :: charge, tpiba2, ewaldg, ewaldr, dtau(3), alpha, &
r(3,mxr), r2(mxr), rmax, rr, upperbound, fact
! total ionic charge in the cell
! length in reciprocal space
! ewald energy computed in reciprocal space
@ -77,59 +84,59 @@ function ewald (alat, nat, ntyp, ityp, zv, at, bg, tau, omega, g, &
! the maximum radius to consider real space sum
! buffer variable
! used to optimize alpha
complex(DP) :: rhon
real(DP), external :: qe_erfc
COMPLEX(DP) :: rhon
REAL(DP), EXTERNAL :: qe_erfc
!
tpiba2 = (tpi / alat) **2
charge = 0.d0
do na = 1, nat
charge = charge+zv (ityp (na) )
enddo
DO na = 1, nat
charge = charge + zv( ityp(na) )
ENDDO
alpha = 2.9d0
100 alpha = alpha - 0.1d0
!
! choose alpha in order to have convergence in the sum over G
! upperbound is a safe upper bound for the error in the sum over G
!
if (alpha.le.0.d0) call errore ('ewald', 'optimal alpha not found', 1)
upperbound = 2.d0 * charge**2 * sqrt (2.d0 * alpha / tpi) * qe_erfc ( &
sqrt (tpiba2 * gcutm / 4.d0 / alpha) )
if (upperbound.gt.1.0d-7) goto 100
IF ( alpha <= 0.d0 ) CALL errore( 'ewald', 'optimal alpha not found', 1 )
upperbound = 2.d0 * charge**2 * SQRT(2.d0 * alpha / tpi) * qe_erfc( &
SQRT(tpiba2 * gcutm / 4.d0 / alpha) )
IF (upperbound > 1.0d-7) GOTO 100
!
! G-space sum here.
! Determine if this processor contains G=0 and set the constant term
!
IF (do_cutoff_2D) then ! cutoff ewald sums
CALL cutoff_ewald(alpha, ewaldg, omega)
IF ( do_cutoff_2D ) THEN ! cutoff ewald sums
CALL cutoff_ewald( alpha, ewaldg, omega )
ELSE
if (gstart==2) then
IF ( gstart==2 ) THEN
ewaldg = - charge**2 / alpha / 4.0d0
else
ELSE
ewaldg = 0.0d0
endif
if (gamma_only) then
ENDIF
IF ( gamma_only ) THEN
fact = 2.d0
else
ELSE
fact = 1.d0
end if
do ng = gstart, ngm
ENDIF
DO ng = gstart, ngm
rhon = (0.d0, 0.d0)
do nt = 1, ntyp
rhon = rhon + zv (nt) * CONJG(strf (ng, nt) )
enddo
ewaldg = ewaldg + fact * abs (rhon) **2 * exp ( - gg (ng) * tpiba2 / &
alpha / 4.d0) / gg (ng) / tpiba2
enddo
DO nt = 1, ntyp
rhon = rhon + zv(nt) * CONJG(strf(ng, nt))
ENDDO
ewaldg = ewaldg + fact * ABS(rhon)**2 * EXP( - gg(ng) * tpiba2 / &
alpha / 4.d0) / gg(ng) / tpiba2
ENDDO
ewaldg = 2.d0 * tpi / omega * ewaldg
!
! Here add the other constant term
!
if (gstart.eq.2) then
do na = 1, nat
ewaldg = ewaldg - zv (ityp (na) ) **2 * sqrt (8.d0 / tpi * &
IF (gstart==2) THEN
DO na = 1, nat
ewaldg = ewaldg - zv(ityp(na))**2 * SQRT(8.d0 / tpi * &
alpha)
enddo
endif
ENDDO
ENDIF
ENDIF
!
! R-space sum here
@ -143,37 +150,38 @@ function ewald (alat, nat, ntyp, ityp, zv, at, bg, tau, omega, g, &
CALL block_distribute( nat, me_bgrp, nproc_bgrp, na_s, na_e, mykey )
ewaldr = 0.d0
IF ( mykey == 0 ) THEN
rmax = 4.d0 / sqrt (alpha) / alat
rmax = 4.d0 / SQRT(alpha) / alat
!
! with this choice terms up to ZiZj*erfc(4) are counted (erfc(4)=2x10^-8
!
do na = na_s, na_e
do nb = 1, nat
dtau (:) = tau (:, na) - tau (:, nb)
DO na = na_s, na_e
DO nb = 1, nat
dtau(:) = tau(:,na) - tau(:,nb)
!
! generates nearest-neighbors shells
!
call rgen (dtau, rmax, mxr, at, bg, r, r2, nrm)
CALL rgen( dtau, rmax, mxr, at, bg, r, r2, nrm )
!
! and sum to the real space part
!
do nr = 1, nrm
rr = sqrt (r2 (nr) ) * alat
ewaldr = ewaldr + zv (ityp (na) ) * zv (ityp (nb) ) * qe_erfc ( &
sqrt (alpha) * rr) / rr
enddo
enddo
enddo
endif
DO nr = 1, nrm
rr = SQRT(r2(nr)) * alat
ewaldr = ewaldr + zv(ityp(na)) * zv(ityp(nb)) * &
qe_erfc(SQRT(alpha) * rr) / rr
ENDDO
ENDDO
ENDDO
ENDIF
ewald = 0.5d0 * e2 * (ewaldg + ewaldr)
if ( do_comp_mt ) ewald = ewald + wg_corr_ewald ( omega, ntyp, ngm, zv, strf )
IF ( do_comp_mt ) ewald = ewald + wg_corr_ewald( omega, ntyp, ngm, zv, strf )
!
call mp_sum( ewald, intra_bgrp_comm )
! call mp_sum( ewaldr, intra_bgrp_comm )
! call mp_sum( ewaldg, intra_bgrp_comm )
CALL mp_sum( ewald, intra_bgrp_comm )
! CALL mp_sum( ewaldr, intra_bgrp_comm )
! CALL mp_sum( ewaldg, intra_bgrp_comm )
! WRITE( stdout,'(/5x,"alpha used in ewald term: ",f4.2/
! + 5x,"R-space term: ",f12.7,5x,"G-space term: ",f12.7/)')
! + alpha, ewaldr, ewaldg
return
end function ewald
RETURN
!
END FUNCTION ewald

182
PW/src/ewald_dipole.f90
View File

@ -7,135 +7,135 @@
!
!
!-----------------------------------------------------------------------
subroutine ewald_dipole (tens,dipole)
SUBROUTINE ewald_dipole( tens, dipole )
!-----------------------------------------------------------------------
!! Calculates the ewald field on each atom due to the presence of dipole, or
!! the electic field on each atom due to the ionic charge of other atoms,
!! with both G- and R-space terms.
!! Determines optimal alpha. Should hopefully work for any structure.
!
! Calculates the ewald field on each atom due to the presence of dipole, or
! the electic field on each atom due to the ionic charge of other atoms,
! with both G- and R-space terms.
! Determines optimal alpha. Should hopefully work for any structure.
USE kinds, ONLY: DP
USE gvect, ONLY: gcutm, gstart, ngm, g, gg
USE constants, ONLY: tpi, e2, fpi, pi
USE cell_base, ONLY: tpiba2, omega, alat, at, bg
USE ions_base, ONLY: nat, ntyp => nsp, ityp, tau
USE vlocal, ONLY: strf
USE mp_bands, ONLY: intra_bgrp_comm
USE mp, ONLY: mp_sum
!
IMPLICIT NONE
!
USE kinds , ONLY : dp
USE gvect , ONLY : gcutm, gstart, ngm, g, gg
USE constants , ONLY : tpi, e2, fpi, pi
USE cell_base , ONLY : tpiba2, omega, alat, at, bg
USE ions_base, ONLY : nat, ntyp => nsp, ityp, tau
USE vlocal , ONLY : strf
USE mp_bands, ONLY : intra_bgrp_comm
USE mp, ONLY : mp_sum
COMPLEX(DP) :: tens(nat,3,3)
!! the ewald field/electric field
REAL(DP) :: dipole(ntyp)
!! the dipole
REAL(DP) :: charge, eta, arg, upperbound, temp
COMPLEX(DP) :: rhon
REAL(DP), EXTERNAL :: qe_erfc
COMPLEX(DP), ALLOCATABLE:: ewaldg(:,:,:), ewaldr(:,:,:)
INTEGER :: alpha, beta, na, ng, nt, ipol, nb, nrm, nr
!
implicit none
INTEGER, PARAMETER :: mxr = 50
REAL (DP) :: r(3,mxr), r2(mxr), rmax, rr, dtau(3)
REAL (DP) :: expcoeff
COMPLEX(DP) :: carg, recarg, recarg_dgg
!
real(DP) :: dipole(ntyp),charge, eta, arg, upperbound, temp
complex(DP) :: tens(nat,3,3)
complex(DP) :: rhon
real(DP), external :: qe_erfc
complex(DP), allocatable:: ewaldg(:,:,:), ewaldr(:,:,:)
integer :: alpha, beta, na, ng, nt, ipol, nb, nrm, nr
integer, parameter :: mxr = 50
real (DP) :: r(3,mxr), r2(mxr), rmax, rr, dtau(3)
real (DP) :: expcoeff
complex(DP) :: carg, recarg, recarg_dgg
allocate (ewaldg(nat,3,3))
allocate (ewaldr(nat,3,3))
ewaldg=(0.d0,0.d0)
ewaldr=(0.d0,0.d0)
! e2=1.d0 !hartree
ALLOCATE( ewaldg(nat,3,3) )
ALLOCATE( ewaldr(nat,3,3) )
!
ewaldg = (0.d0,0.d0)
ewaldr = (0.d0,0.d0)
!
! e2=1.d0 !hartree
charge = 0.d0
do na = 1, nat
charge = charge+dipole (ityp (na) )
enddo
DO na = 1, nat
charge = charge + dipole(ityp(na))
ENDDO
eta = 2.9d0
do
DO
eta = eta - 0.1d0
!
! choose alpha in order to have convergence in the sum over G
! upperbound is a safe upper bound for the error in the sum over G
!
if (eta.le.0.d0) call errore ('ewald_dipole', 'optimal eta not found', 1)
upperbound = 2.d0 * charge**2 * sqrt (2.d0 * eta / tpi) &
* qe_erfc ( sqrt (tpiba2 * gcutm / 4.d0 / eta) )
if (upperbound.le.1.0d-7) exit
enddo
IF ( eta <= 0.d0 ) CALL errore( 'ewald_dipole', 'optimal eta not found', 1 )
upperbound = 2.d0 * charge**2 * SQRT(2.d0 * eta / tpi) &
* qe_erfc ( SQRT(tpiba2 * gcutm / 4.d0 / eta) )
IF ( upperbound <= 1.0d-7 ) EXIT
ENDDO
!
! G-space sum here.
do ng = gstart, ngm
!
DO ng = gstart, ngm
rhon = (0.d0, 0.d0)
expcoeff = exp ( - gg (ng) * tpiba2 * 0.25d0 / eta )
do nt = 1, ntyp
rhon = rhon + dipole (nt) * CONJG(strf (ng, nt) )
enddo
do na=1, nat
arg = (g (1, ng) * tau (1, na) + g (2, ng) * tau (2, na) &
+ g (3, ng) * tau (3, na) ) * tpi
carg = CMPLX(cos(arg), -sin(arg),kind=DP)
expcoeff = EXP( - gg(ng) * tpiba2 * 0.25d0 / eta )
DO nt = 1, ntyp
rhon = rhon + dipole(nt) * CONJG(strf(ng, nt))
ENDDO
DO na=1, nat
arg = (g(1, ng) * tau(1, na) + g(2, ng) * tau(2, na) &
+ g(3, ng) * tau(3, na) ) * tpi
carg = CMPLX(COS(arg), -SIN(arg),KIND=DP)
recarg = rhon*expcoeff*carg
recarg_dgg = recarg / gg(ng)
do alpha = 1,3
do beta=1,3
ewaldg(na , alpha, beta) = ewaldg(na, alpha, beta) &
DO alpha = 1,3
DO beta=1,3
ewaldg(na, alpha, beta) = ewaldg(na, alpha, beta) &
- recarg_dgg * g(alpha,ng) * g(beta,ng)
enddo
ewaldg(na , alpha, alpha) = ewaldg(na, alpha, alpha) &
ENDDO
!
ewaldg(na, alpha, alpha) = ewaldg(na, alpha, alpha) &
+ 1.d0/3.d0 * recarg
enddo
enddo
enddo
ENDDO
ENDDO
ENDDO
ewaldg = e2 / 2.d0 * fpi / omega * ewaldg !Temp to compare with paratec
! ewaldg = e2 * fpi / omega * ewaldg
! ewaldg = e2 * fpi / omega * ewaldg
!
call mp_sum( ewaldg, intra_bgrp_comm )
CALL mp_sum( ewaldg, intra_bgrp_comm )
!
! R-space sum here (only for the processor that contains G=0)
!
ewaldr = 0.d0
if (gstart.eq.2) then
rmax = 4.d0 / sqrt (eta) / alat
IF ( gstart==2 ) THEN
rmax = 4.d0 / SQRT(eta) / alat
!
! with this choice terms up to ZiZj*erfc(4) are counted (erfc(4)=2x10^-8
!
do na = 1, nat
do nb = 1, nat
do ipol = 1, 3
dtau (ipol) = tau (ipol, na) - tau (ipol, nb)
enddo
DO na = 1, nat
DO nb = 1, nat
DO ipol = 1, 3
dtau(ipol) = tau(ipol, na) - tau(ipol, nb)
ENDDO
!
! generates nearest-neighbors shells
!
call rgen (dtau, rmax, mxr, at, bg, r, r2, nrm)
CALL rgen(dtau, rmax, mxr, at, bg, r, r2, nrm)
!
! and sum to the real space part
! and sum to the REAL space part
!
r = r * alat
do nr = 1, nrm
rr = sqrt (r2 (nr) ) * alat
temp= dipole (ityp (na)) * ( 3.d0 / rr**3 * qe_erfc ( sqrt (eta) * rr) &
+ (6.d0 * sqrt (eta/pi) * 1.d0 / rr*2 + 4.d0 * sqrt (eta**3/pi)) &
* exp(-eta* rr**2))
do alpha=1,3
do beta=1,3
DO nr = 1, nrm
rr = SQRT( r2(nr) ) * alat
temp= dipole(ityp(na)) * ( 3.d0 / rr**3 * qe_erfc( SQRT(eta) * rr) &
+ (6.d0 * SQRT(eta/pi) * 1.d0 / rr*2 + 4.d0 * SQRT(eta**3/pi)) &
* EXP(-eta* rr**2))
DO alpha=1,3
DO beta=1,3
ewaldr(na, alpha,beta) = ewaldr(na, alpha,beta) &
+ temp*r(alpha,nr)*r(beta,nr) / rr**2
enddo
ENDDO
ewaldr(na, alpha,alpha)= ewaldr(na, alpha,alpha) &
- 1.d0/3.d0 * temp
enddo
enddo
enddo
enddo
endif
ewaldr = e2 * ewaldr
ENDDO
ENDDO
ENDDO
ENDDO
ENDIF
ewaldr = e2 * ewaldr
!
call mp_sum( ewaldr, intra_bgrp_comm )
CALL mp_sum( ewaldr, intra_bgrp_comm )
!
tens=ewaldg+ewaldr
end subroutine ewald_dipole
tens = ewaldg + ewaldr
!
END SUBROUTINE ewald_dipole

145
PW/src/extfield.f90
View File

@ -8,89 +8,110 @@
!--------------------------------------------------------------------------
!
MODULE extfield
!
! ... The quantities needed in calculations with external field
!! The quantities needed in calculations with external field.
!
USE kinds, ONLY : DP
!
SAVE
!
LOGICAL :: &
tefield, &! if .TRUE. a finite electric field is added to the
! local potential
dipfield, &! if .TRUE. the dipole field is subtracted
! TB
relaxz, &! relax in z direction
block, &! add potential barrier
gate ! if .TRUE. and system is charged, charge is represented
! with charged plate (gate)
INTEGER :: &
edir ! direction of the field
REAL(DP) :: &
emaxpos, &! position of the maximum of the field (0<emaxpos<1)
eopreg, &! amplitude of the inverse region (0<eopreg<1)
eamp, &! field amplitude (in a.u.) (1 a.u. = 51.44 10^11 V/m)
etotefield, &! energy correction due to the field
el_dipole, &! electronic_dipole used when dipole correction is on
ion_dipole, &! ionic_dipole used when dipole correction is on
tot_dipole, &! total dipole used when dipole correction is on
! TB
zgate, &! position of charged plate
block_1, &! blocking potential
block_2, &
block_height, &
etotgatefield ! energy correction due to the gate
REAL(DP), ALLOCATABLE :: &
forcefield(:,:), &
forcegate(:,:) ! TB gate forces
LOGICAL :: tefield
!! if .TRUE. a finite electric field is added to the
!! local potential
LOGICAL :: dipfield
!! if .TRUE. the dipole field is subtracted
! TB
LOGICAL :: relaxz
!! relax in z direction
LOGICAL :: block
!! add potential barrier
LOGICAL :: gate
!! if .TRUE. and system is charged, charge is represented
!! with charged plate (gate)
!
INTEGER :: edir
!! direction of the field
REAL(DP) :: emaxpos
!! position of the maximum of the field (0<emaxpos<1)
REAL(DP) :: eopreg
!! amplitude of the inverse region (0<eopreg<1)
REAL(DP) :: eamp
!! field amplitude (in a.u.) (1 a.u. = 51.44 10^11 V/m)
REAL(DP) :: etotefield
!! energy correction due to the field
REAL(DP) :: el_dipole
!! electronic_dipole, used when dipole correction is on
REAL(DP) :: ion_dipole
!! ionic_dipole, used when dipole correction is on
REAL(DP) :: tot_dipole
!! total dipole, used when dipole correction is on
! TB
REAL(DP) :: zgate
!! position of charged plate
REAL(DP) :: block_1
!! ...blocking potential
REAL(DP) :: block_2
!! ...blocking potential
REAL(DP) :: block_height
!! ...blocking potential
REAL(DP) :: etotgatefield
!! energy correction due to the gate
REAL(DP), ALLOCATABLE :: forcefield(:,:)
!
REAL(DP), ALLOCATABLE :: forcegate(:,:)
!! TB gate forces
!
CONTAINS
!
!------------------------------------------------------------------------------!
!
FUNCTION saw(emaxpos,eopreg,x) RESULT (sawout)
FUNCTION saw( emaxpos, eopreg, x ) RESULT ( sawout )
!
IMPLICIT NONE
REAL(DP) :: emaxpos,eopreg,x
!
REAL(DP) :: emaxpos, eopreg, x
REAL(DP) :: y, sawout, z
!
z = x - emaxpos
y = z - floor(z)
if (y.le.eopreg) then
y = z - FLOOR(z)
!
IF (y <= eopreg) THEN
!
sawout = (0.5_DP - y/eopreg) * (1._DP-eopreg)
else
!
! I would use: sawout = y - 0.5_DP * ( 1.0_DP + eopreg )
!
!
ELSE
!
! I would use: sawout = y - 0.5_DP * ( 1.0_DP + eopreg )
!
sawout = (-0.5_DP + (y-eopreg)/(1._DP-eopreg)) * (1._DP-eopreg)
end if
!
ENDIF
!
END FUNCTION saw
!TB - start
!------------------------------------------------------------------------------!
!mopopla - add a potential of a charged plate (kflag = .true.)
! or the compensating background charge (kflag = .false.)
! I split those in order to plot both independently
! cite PRB 89, 245406 (2014)
!
FUNCTION mopopla(zgate,x,kflag) RESULT (mopoplaout)
FUNCTION mopopla( zgate, x, kflag ) RESULT ( mopoplaout )
!------------------------------------------------------
!! Add a potential of a charged plate (kflag = .true.)
!! or the compensating background charge (kflag = .false.).
!! Cite PRB 89, 245406 (2014).
! I split those in order to plot both independently
!
IMPLICIT NONE
!
REAL(DP) :: zgate,x
REAL(DP) :: mopoplaout, z
LOGICAL :: kflag
!
DO ! is x within the cell?
IF (x>1.0) x=x-1.0
IF (x<0.0) x=x+1.0
IF (x<=1.0.and.x>=0.0) EXIT
IF ( x>1.0 ) x=x-1.0
IF ( x<0.0 ) x=x+1.0
IF ( x<=1.0.AND.x>=0.0 ) EXIT
ENDDO
!
z = (x - zgate)
!
!Charged-plate
! if z < 0, we are below the plane
! z > 0, above
@ -98,22 +119,22 @@ CONTAINS
! in order to make it periodic
! z > 0.5, the same as for z < 0
!
IF (z<=-0.5) z=z+1
IF (z>=0.5) z=z-1
IF (z.LE.0) THEN
IF ( z<=-0.5 ) z=z+1
IF ( z>=0.5 ) z=z-1
IF ( z<=0 ) THEN
IF (kflag) THEN
mopoplaout = ( 1.0_DP*z )
ELSE
mopoplaout = ( z**2 )
ENDIF
ELSE
IF (kflag) THEN
IF ( kflag ) THEN
mopoplaout = ( -1.0_DP*z )
ELSE
mopoplaout = ( z**2 )
ENDIF
ENDIF
!
END FUNCTION mopopla
!
END MODULE extfield

4204
PW/src/exx.f90
View File

File diff suppressed because it is too large Load Diff

760
PW/src/exx_base.f90
View File

File diff suppressed because it is too large Load Diff

312
PW/src/fcp.f90
View File

@ -5,33 +5,34 @@
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
! Original version by Minoru Otani (AIST) and Nicephore Bonnet (AIST).
!
! This module controls the Fictitious Charge Particle (FCP) for constant-mu
! method developed by N. Bonnet, T. Morishita, O. Sugino, and M. Otani
! (see PRL 109, 266101 [2012]).
!
! Constant-mu scheme with the boundary condition 'bc2' and 'bc3' enables
! description of the system connected to a potentiostat which preserves
! the Fermi energy of the system as the target Fermi energy (mu).
!
! MDIIS algorithm is implemented for relax-calculation
! by Satomichi Nishihara (2016)
!
!----------------------------------------------------------------------------
MODULE fcp
!----------------------------------------------------------------------------
!! Original version by Minoru Otani (AIST) and Nicephore Bonnet (AIST).
!
USE kinds, ONLY : DP
USE io_global, ONLY : stdout
USE io_files, ONLY : seqopn
USE constants, ONLY : amu_ry, ry_to_kelvin, au_ps, rytoev, e2, fpi
USE control_flags, ONLY : tolp
USE dynamics_module, ONLY : dt, delta_t, nraise, &
control_temp, thermostat
USE fcp_variables, ONLY : fcp_mu, fcp_mass, fcp_temperature, &
fcp_relax, fcp_relax_step, fcp_relax_crit, fcp_mdiis_size, fcp_mdiis_step
USE mdiis, ONLY : mdiis_type, allocate_mdiis, deallocate_mdiis, update_by_mdiis
!! This module controls the Fictitious Charge Particle (FCP) for constant-mu
!! method developed by N. Bonnet, T. Morishita, O. Sugino, and M. Otani
!! (see PRL 109, 266101 [2012]).
!
!! Constant-mu scheme with the boundary condition 'bc2' and 'bc3' enables
!! description of the system connected to a potentiostat which preserves
!! the Fermi energy of the system as the target Fermi energy (mu).
!
!! MDIIS algorithm is implemented for relax-calculation by
!! Satomichi Nishihara (2016).
!
USE kinds, ONLY: DP
USE io_global, ONLY: stdout
USE io_files, ONLY: seqopn
USE constants, ONLY: amu_ry, ry_to_kelvin, au_ps, rytoev, e2, fpi
USE control_flags, ONLY: tolp
USE dynamics_module, ONLY: dt, delta_t, nraise, &
control_temp, thermostat
USE fcp_variables, ONLY: fcp_mu, fcp_mass, fcp_temperature, fcp_relax, &
fcp_relax_step, fcp_relax_crit, fcp_mdiis_size, &
fcp_mdiis_step
USE mdiis, ONLY: mdiis_type, allocate_mdiis, deallocate_mdiis, &
update_by_mdiis
!
IMPLICIT NONE
!
@ -39,45 +40,53 @@ MODULE fcp
!
PRIVATE
!
LOGICAL :: vel_defined
! if true, vel is used rather than tau_old to do the next step
REAL(DP) :: tau, tau_old, tau_new
REAL(DP) :: vel, acc
LOGICAL :: vel_defined
!! if .TRUE., vel is used rather than tau_old to do the next step
REAL(DP) :: tau
!! present position
REAL(DP) :: tau_old
!! previous step position
REAL(DP) :: tau_new
!! next step position
REAL(DP) :: vel
!! velocity
REAL(DP) :: acc
!! acceleration
!
! ... data of MDIIS
LOGICAL :: init_mdiis = .FALSE.
!!
TYPE(mdiis_type) :: mdiist
!!
!
PUBLIC :: fcp_line_minimisation, fcp_verlet, fcp_summary
PUBLIC :: fcp_mdiis_update, fcp_mdiis_end
!
CONTAINS
CONTAINS
!
! ... public methods
!
!------------------------------------------------------------------------
!------------------------------------------------------------------------------
SUBROUTINE fcp_verlet()
!------------------------------------------------------------------------
!-----------------------------------------------------------------------------
!! This routine performs one step of molecular dynamics evolution by
!! using the Verlet algorithm.
!
! ... This routine performs one step of molecular dynamics evolution
! ... using the Verlet algorithm.
!! Parameters (in \(\textrm{dynamics_module}\) and \(\textrm{fcp_variables}\)):
!
! ... Parameters:
! ... mass mass of the atoms
! ... dt time step
! ... temperature starting temperature
! ... The starting velocities of atoms are set accordingly
! ... to the starting temperature, in random directions.
! ... The initial velocity distribution is therefore a
! ... constant.
!! * mass: mass of the atoms;
!! * dt: time step;
!! * temperature: starting temperature. The starting velocities of atoms
!! are set accordingly to the starting temperature, in random directions.
!! Therefore the initial velocity distribution is a constant.
!
! ... Dario Alfe' 1997 and Carlo Sbraccia 2004-2006
!! Dario Alfe' 1997 and Carlo Sbraccia 2004-2006
!
USE cell_base, ONLY : alat
USE control_flags, ONLY : istep
USE ener, ONLY : ef
USE klist, ONLY : nelec, tot_charge
USE ions_base, ONLY : nat, ityp, zv
USE cell_base, ONLY: alat
USE control_flags, ONLY: istep
USE ener, ONLY: ef
USE klist, ONLY: nelec, tot_charge
USE ions_base, ONLY: nat, ityp, zv
!
IMPLICIT NONE
!
@ -85,24 +94,23 @@ CONTAINS
REAL(DP) :: temp_new, temp_av, elapsed_time, force
! istep counts all MD steps, including those of previous runs
LOGICAL :: file_exists
!
tau = nelec
vel_defined = .true.
temp_av = 0.D0
!
vel_defined = .TRUE.
temp_av = 0.0_DP
!
CALL seqopn( 4, 'fcp_md', 'FORMATTED', file_exists )
!
IF ( file_exists ) THEN
!
! ... the file is read : simulation is continuing
! ... the file is read: simulation is continuing
!
READ( UNIT = 4, FMT = * ) istep, tau_old
!
vel_defined = .false.
vel_defined = .FALSE.
!
READ( UNIT = 4, FMT = * ) &
temp_new, temp_av, fcp_mass, elapsed_time
READ( UNIT = 4, FMT = * ) temp_new, temp_av, fcp_mass, elapsed_time
!
CLOSE( UNIT = 4, STATUS = 'KEEP' )
!
@ -110,7 +118,7 @@ CONTAINS
!
CLOSE( UNIT = 4, STATUS = 'DELETE' )
!
! ... the file is absent : simulation is starting from scratch
! ... the file is absent: simulation is starting from scratch
!
CALL md_init()
!
@ -118,7 +126,7 @@ CONTAINS
!
! ... elapsed_time is in picoseconds
!
elapsed_time = elapsed_time + dt*2.D0*au_ps
elapsed_time = elapsed_time + dt*2.0_DP*au_ps
istep = istep + 1
!
IF ( control_temp ) CALL apply_thermostat()
@ -139,21 +147,21 @@ CONTAINS
!
ELSE
!
tau_new = 2.D0*tau - tau_old + acc * dt**2
tau_new = 2.0_DP*tau - tau_old + acc * dt**2
!
ENDIF
!
! ... the linear momentum and the kinetic energy are computed here
!
vel = ( tau_new - tau_old ) / ( 2.D0*dt )
vel = ( tau_new - tau_old ) / ( 2.0_DP*dt )
!
ekin = 0.5D0 * fcp_mass * vel**2
ekin = 0.5_DP * fcp_mass * vel**2
!
ekin = ekin*alat**2
!
! ... find the new temperature and update the average
!
temp_new = 2.D0 * ekin * ry_to_kelvin
temp_new = 2.0_DP * ekin * ry_to_kelvin
!
temp_av = temp_av + temp_new
!
@ -196,8 +204,7 @@ CONTAINS
'(/,5X,"Starting fcp temperature",T27," = ",F8.2," K")' ) &
fcp_temperature
!
SELECT CASE( trim( thermostat ) )
!
SELECT CASE( TRIM( thermostat ) )
CASE( 'andersen', 'Andersen' )
!
WRITE( UNIT = stdout, &
@ -223,17 +230,16 @@ CONTAINS
WRITE( UNIT = stdout, &
FMT = '(/,5X,"fcp temperature is controlled by ", &
& "velocity rescaling (",A,")"/)' )&
trim( thermostat )
TRIM( thermostat )
!
END SELECT
!
ENDIF
!
WRITE( UNIT = stdout, &
FMT = '(5X,"fcp_mass = ",F8.2)' ) fcp_mass
WRITE( UNIT = stdout, FMT = '(5X,"fcp_mass = ",F8.2)' ) fcp_mass
!
WRITE( UNIT = stdout, &
FMT = '(5X,"Time step",T27," = ",F8.2," a.u.,",F8.4, &
FMT = '(5X,"Time step",T27," = ",F8.2," a.u.,",F8.4, &
& " femto-seconds")' ) dt, dt*2.D+3*au_ps
!
! ... masses in rydberg atomic units
@ -243,65 +249,70 @@ CONTAINS
! ... initial thermalization. N.B. tau is in units of alat
!
CALL start_therm()
vel_defined = .true.
vel_defined = .TRUE.
!
temp_new = fcp_temperature
!
temp_av = 0.D0
temp_av = 0.0_DP
!
ELSE
!
vel = 0.0_DP
vel_defined = .true.
vel_defined = .TRUE.
!
ENDIF
!
elapsed_time = 0.D0
elapsed_time = 0.0_DP
!
END SUBROUTINE md_init
!
!
!--------------------------------------------------------------------
SUBROUTINE apply_thermostat()
!--------------------------------------------------------------------
!
USE random_numbers, ONLY : randy, gauss_dist
USE random_numbers, ONLY: randy, gauss_dist
!
IMPLICIT NONE
!
INTEGER :: nat_moved
REAL(DP) :: sigma, kt
!
IF(.not.vel_defined)THEN
IF (.NOT. vel_defined) THEN
vel = (tau - tau_old) / dt
ENDIF
!
SELECT CASE( trim( thermostat ) )
SELECT CASE( TRIM( thermostat ) )
CASE( 'rescaling' )
IF ( abs (temp_new-fcp_temperature) > tolp ) THEN
!
IF ( ABS(temp_new-fcp_temperature) > tolp ) THEN
!
WRITE( UNIT = stdout, &
FMT = '(/,5X,"FCP Velocity rescaling: T (",F6.1,"K) ", &
FMT = '(/,5X,"FCP Velocity rescaling: T (",F6.1,"K) ", &
& "out of range, reset to " ,F6.1)' ) &
temp_new, fcp_temperature
temp_new, fcp_temperature
CALL thermalize( 0, temp_new, fcp_temperature )
!
ENDIF
!
CASE( 'rescale-v', 'rescale-V', 'rescale_v', 'rescale_V' )
IF ( mod( istep, nraise ) == 0 ) THEN
!
IF ( MOD(istep, nraise) == 0 ) THEN
!
temp_av = temp_av / dble( nraise )
temp_av = temp_av / DBLE( nraise )
!
WRITE( UNIT = stdout, &
FMT = '(/,5X,"FCP Velocity rescaling: average T on ",i3, &
&" steps (",F6.1,"K) reset to ",F6.1)' ) &
nraise, temp_av, fcp_temperature
!
CALL thermalize( 0, temp_new, fcp_temperature )
!
temp_av = 0.D0
temp_av = 0.0_DP
!
ENDIF
!
CASE( 'rescale-T', 'rescale-t', 'rescale_T', 'rescale_t' )
!
IF ( delta_t > 0 ) THEN
!
fcp_temperature = temp_new*delta_t
@ -313,13 +324,15 @@ CONTAINS
CALL thermalize( 0, temp_new, fcp_temperature )
!
ENDIF
!
CASE( 'reduce-T', 'reduce-t', 'reduce_T', 'reduce_t' )
IF ( mod( istep, nraise ) == 0 .and. delta_t < 0 ) THEN
!
IF ( MOD( istep, nraise ) == 0 .AND. delta_t < 0 ) THEN
!
fcp_temperature = temp_new + delta_t
!
WRITE( UNIT = stdout, &
FMT = '(/,5X,"FCP Thermalization: T (",F6.1,"K) reduced ",&
FMT = '(/,5X,"FCP Thermalization: T (",F6.1,"K) reduced ",&
& "by ",F6.3)' ) temp_new, -delta_t
!
CALL thermalize( 0, temp_new, fcp_temperature )
@ -338,15 +351,16 @@ CONTAINS
kt = fcp_temperature / ry_to_kelvin
nat_moved = 0
!
IF ( randy() < 1.D0 / dble( nraise ) ) THEN
IF ( randy() < 1.D0 / DBLE( nraise ) ) THEN
!
nat_moved = nat_moved + 1
sigma = sqrt( kt / fcp_mass )
sigma = SQRT( kt / fcp_mass )
!
! ... N.B. velocities must in a.u. units of alat and are zero
! ... for fixed ions
!
vel = gauss_dist( 0.D0, sigma ) / alat
!
ENDIF
!
IF ( nat_moved > 0) WRITE( UNIT = stdout, &
@ -361,17 +375,17 @@ CONTAINS
!
! ... the old positions are updated to reflect the new velocities
!
IF(.not.vel_defined)THEN
IF (.NOT. vel_defined) THEN
tau_old = tau - vel * dt
ENDIF
!
END SUBROUTINE apply_thermostat
!
!
!-----------------------------------------------------------------------
SUBROUTINE start_therm()
!-----------------------------------------------------------------------
!
! ... Starting thermalization of the system
!! Starting thermalization of the system.
!
USE symm_base, ONLY : invsym, nsym, irt
USE cell_base, ONLY : alat
@ -384,12 +398,13 @@ CONTAINS
! ... next command prevents different MD runs to start
! ... with exactly the same "random" velocities
!
call set_random_seed ( )
CALL set_random_seed()
!
kt = fcp_temperature / ry_to_kelvin
!
! ... starting velocities have a Maxwell-Boltzmann distribution
!
sigma = sqrt( kt / fcp_mass )
sigma = SQRT( kt / fcp_mass )
!
! ... N.B. velocities must in a.u. units of alat
!
@ -417,8 +432,14 @@ CONTAINS
!
IMPLICIT NONE
!
REAL(DP), INTENT(in) :: system_temp, required_temp
INTEGER, INTENT(in) :: nraise
REAL(DP), INTENT(IN) :: system_temp
!! Temperature of the system
REAL(DP), INTENT(IN) :: required_temp
!! Required temperature
INTEGER, INTENT(IN) :: nraise
!! number of steps
!
! ... local variables
!
REAL(DP) :: aux
!
@ -428,10 +449,10 @@ CONTAINS
! ... the "rise time" is tau=nraise*dt so dt/tau=1/nraise
! ... Equivalent to traditional rescaling if nraise=1
!
IF ( system_temp > 0.D0 .and. required_temp > 0.D0 ) THEN
IF ( system_temp > 0.D0 .AND. required_temp > 0.D0 ) THEN
!
aux = sqrt( 1.d0 + (required_temp / system_temp - 1.d0) * &
(1.D0/dble (nraise) ) )
aux = SQRT( 1.d0 + (required_temp / system_temp - 1.d0) * &
( 1.D0/DBLE(nraise) ) )
!
ELSE
!
@ -443,9 +464,9 @@ CONTAINS
!
! ... rescale the velocities by a factor 3 / 2KT / Ek
!
IF ( system_temp > 0.D0 .and. required_temp > 0.D0 ) THEN
IF ( system_temp > 0.D0 .AND. required_temp > 0.D0 ) THEN
!
aux = sqrt( required_temp / system_temp )
aux = SQRT( required_temp / system_temp )
!
ELSE
!
@ -462,9 +483,8 @@ CONTAINS
!------------------------------------------------------------------------
SUBROUTINE fcp_line_minimisation( conv_fcp )
!------------------------------------------------------------------------
!
! ... This routine performs one step of relaxation
! ... using the steepest descent algorithm.
!! This routine performs one step of relaxation by using the
!! steepest descent algorithm.
!
USE control_flags, ONLY : iverbosity
USE ener, ONLY : ef
@ -473,9 +493,14 @@ CONTAINS
USE cell_base, ONLY : at, alat
!
IMPLICIT NONE
!
LOGICAL, INTENT(OUT) :: conv_fcp
REAL(DP) :: force, n_tmp, max_tot_charge, capacitance, &
ionic_charge
!!
!
! ... local variables
!
REAL(DP) :: force, n_tmp, max_tot_charge, capacitance, &
ionic_charge
!
LOGICAL , SAVE :: firstcall = .TRUE.
REAL(DP), SAVE :: force0 = 0.0_DP
@ -488,7 +513,9 @@ CONTAINS
!
capacitance = (at(1,1) * at(2,2) - at(2,1) * at(1,2)) * alat**2 &
/ (alat * at(3,3) / 2._DP ) / fpi
max_tot_charge = abs( capacitance * force / e2 )
!
max_tot_charge = ABS( capacitance * force / e2 )
!
IF ( firstcall .OR. ABS( force0 - force ) < 1.0D-20 ) THEN
firstcall = .FALSE.
nelec0 = nelec
@ -499,30 +526,31 @@ CONTAINS
nelec = ( nelec * force0 - nelec0 * force ) / ( force0 - force )
nelec0 = n_tmp
force0 = force
END IF
ENDIF
!
ionic_charge = SUM( zv(ityp(1:nat)) )
!
IF ( iverbosity > 1 ) THEN
write( stdout,'(/,5X,"Upper bound for tot_charge:",F12.6)') &
WRITE( stdout,'(/,5X,"Upper bound for tot_charge:",F12.6)') &
max_tot_charge
write( stdout,'(5X,"Original:",F12.6," Expected:",F12.6)') &
WRITE( stdout,'(5X,"Original:",F12.6," Expected:",F12.6)') &
ionic_charge - nelec0, ionic_charge - nelec
END IF
if( nelec-nelec0 < -max_tot_charge ) nelec= nelec0 - max_tot_charge
if( nelec-nelec0 > max_tot_charge ) nelec= nelec0 + max_tot_charge
ENDIF
!
IF ( nelec-nelec0 < -max_tot_charge ) nelec = nelec0 - max_tot_charge
IF ( nelec-nelec0 > max_tot_charge ) nelec = nelec0 + max_tot_charge
tot_charge = ionic_charge - nelec
IF ( iverbosity > 1 ) THEN
write( stdout,'(5X,"Next tot_charge:",F12.6)') tot_charge
END IF
WRITE( stdout,'(5X,"Next tot_charge:",F12.6)') tot_charge
ENDIF
!
conv_fcp = .FALSE.
!
IF( abs( force ) < fcp_relax_crit ) THEN
!
IF( ABS( force ) < fcp_relax_crit ) THEN
conv_fcp = .TRUE.
nelec = nelec0
tot_charge = ionic_charge - nelec
!
END IF
ENDIF
!
WRITE( stdout, FMT = 9001 ) force
!
@ -530,12 +558,12 @@ CONTAINS
!
END SUBROUTINE fcp_line_minimisation
!
!
!------------------------------------------------------------------------
SUBROUTINE fcp_mdiis_update( conv_fcp )
!------------------------------------------------------------------------
!
! ... This routine performs one step of relaxation
! ... using the MDIIS algorithm.
!! This routine performs one step of relaxation by
!! using the MDIIS algorithm.
!
USE control_flags, ONLY : iverbosity
USE ener, ONLY : ef
@ -543,35 +571,43 @@ CONTAINS
USE ions_base, ONLY : nat, ityp, zv
!
IMPLICIT NONE
!
LOGICAL, INTENT(OUT) :: conv_fcp
REAL(DP) :: ionic_charge
REAL(DP) :: nelec0
REAL(DP) :: nelec1(1)
REAL(DP) :: force
REAL(DP) :: force1(1)
!! .TRUE. if converged
!
! ... local variables
!
REAL(DP) :: ionic_charge
REAL(DP) :: nelec0
REAL(DP) :: nelec1(1)
REAL(DP) :: force
REAL(DP) :: force1(1)
!
IF ( .NOT. init_mdiis ) THEN
init_mdiis = .TRUE.
CALL allocate_mdiis(mdiist, fcp_mdiis_size, 1, fcp_mdiis_step, 1)
END IF
CALL allocate_mdiis( mdiist, fcp_mdiis_size, 1, fcp_mdiis_step, 1 )
ENDIF
!
nelec0 = nelec
nelec1(1) = nelec
force = fcp_mu - ef
force1(1) = force
CALL update_by_mdiis(mdiist, nelec1, force1)
!
CALL update_by_mdiis( mdiist, nelec1, force1 )
!
nelec = nelec1(1)
!
ionic_charge = SUM( zv(ityp(1:nat)) )
!
IF ( iverbosity > 1 ) THEN
write( stdout,'(5X,"Original:",F12.6," Expected:",F12.6)') &
ionic_charge - nelec0, ionic_charge - nelec
END IF
ENDIF
!
tot_charge = ionic_charge - nelec
IF ( iverbosity > 1 ) THEN
write( stdout,'(5X,"Next tot_charge:",F12.6)') tot_charge
END IF
ENDIF
!
conv_fcp = .FALSE.
!
@ -581,7 +617,7 @@ CONTAINS
nelec = nelec0
tot_charge = ionic_charge - nelec
!
END IF
ENDIF
!
WRITE( stdout, FMT = 9002 ) force
!
@ -596,12 +632,15 @@ CONTAINS
IMPLICIT NONE
!
IF ( init_mdiis ) THEN
CALL deallocate_mdiis(mdiist)
END IF
CALL deallocate_mdiis( mdiist )
ENDIF
!
END SUBROUTINE fcp_mdiis_end
!
SUBROUTINE fcp_summary ()
!
!------------------------------------------------------------------------
SUBROUTINE fcp_summary()
!------------------------------------------------------------------------
!
USE io_global, ONLY : stdout, ionode
USE constants, ONLY : rytoev, BOHR_RADIUS_ANGS
@ -619,18 +658,18 @@ CONTAINS
IF ( TRIM(fcp_relax) == 'lm' ) THEN
WRITE( UNIT = stdout, FMT = 9058 )
WRITE( UNIT = stdout, FMT = 9059 ) fcp_relax_step
ELSE IF ( TRIM(fcp_relax) == 'mdiis' ) THEN
ELSEIF ( TRIM(fcp_relax) == 'mdiis' ) THEN
WRITE( UNIT = stdout, FMT = 9060 )
WRITE( UNIT = stdout, FMT = 9061 ) fcp_mdiis_size
WRITE( UNIT = stdout, FMT = 9062 ) fcp_mdiis_step
END IF
ENDIF
WRITE( UNIT = stdout, FMT = 9063 ) fcp_relax_crit*rytoev, &
fcp_relax_crit
ELSE IF ( lfcpdyn ) THEN
fcp_relax_crit
ELSEIF ( lfcpdyn ) THEN
WRITE( UNIT = stdout, FMT = '(5x,"-->FCP optimiser activated<--")' )
WRITE( UNIT = stdout, FMT = 9064 ) fcp_mu*rytoev, fcp_mu
WRITE( UNIT = stdout, FMT = 9065 ) tot_charge
END IF
ENDIF
9056 FORMAT( ' Target Fermi energy = ', F9.4,' eV' &
/' = ', F9.4,' Ry')
9057 FORMAT( ' Initial tot_charge = ', F9.6)
@ -646,4 +685,5 @@ CONTAINS
!
END SUBROUTINE fcp_summary
!
!
END MODULE fcp

591
PW/src/find_group.f90
View File

@ -5,281 +5,342 @@
! in the root directory of the present distribution,
! or http://www.gnu.org/copyleft/gpl.txt .
!
SUBROUTINE find_group(nrot,smat,gname,code_group)
!
! Given a group of nrot rotation matrices smat (in cartesian coordinates)
! this routine finds the name of the point group. It assumes but does not
! check that:
! 1) The nrot matrices smat are actually a group.
! 2) The group is one of the thirty-two point groups.
!
USE kinds, ONLY : DP
IMPLICIT NONE
INTEGER :: nrot, code_group
REAL(DP) :: smat(3,3,nrot)
CHARACTER (LEN=11) :: gname, group_name
INTEGER :: noperation(6), irot, ts, tipo_sym
!
! For each possible group operation the function tipo_sym gives a code
! 1 identity,
! 2 inversion,
! 3 proper rotation <> 180,
! 4 proper rotation 180 degrees,
! 5 mirror,
! 6 improper rotation
! the variable noperation counts how many operations are present in the group.
!
noperation=0
DO irot=1,nrot
ts=tipo_sym(smat(1,1,irot))
noperation(ts)=noperation(ts)+1
END DO
IF (noperation(1).ne.1) call errore('find_group','the group has not identity',1)
code_group=0
IF (noperation(2)==0) THEN
!
! There is not inversion
!
IF (nrot==1) THEN
code_group=1 ! C_1
ELSEIF (nrot==2) THEN
IF (noperation(4)==1) code_group=4 ! C_2
IF (noperation(5)==1) code_group=3 ! C_s
ELSEIF (nrot==3) THEN
IF (noperation(3)==2) code_group=5 ! C_3
ELSEIF (nrot==4) THEN
IF (noperation(6)>0) code_group=26 ! S_4
IF (noperation(5)>0.and.code_group==0) code_group=12 ! C_2v
IF (noperation(3)>0.and.code_group==0) code_group=6 ! C_4
IF (noperation(4)>0.and.code_group==0) code_group=8 ! D_2
ELSEIF (nrot==6) THEN
IF (noperation(5)==3) code_group=13 ! C_3v
IF (noperation(5)==1) code_group=17 ! C_3h
IF (noperation(4)==3.and.code_group==0) code_group=9 ! D_3
IF (noperation(3)>0.and.code_group==0) code_group=7 ! C_6
ELSEIF (nrot==8) THEN
IF (noperation(5)==4) code_group=14 ! C_4v
IF (noperation(5)==2) code_group=24 ! D_2d
IF (noperation(3)>0.and.code_group==0) code_group=10 ! D_4
ELSEIF (nrot==12) THEN
IF (noperation(5)==6) code_group=15 ! C_6v
IF (noperation(5)==4) code_group=21 ! D_3h
IF (noperation(4)>6.and.code_group==0) code_group=11 ! D_6
IF (noperation(3)>0.and.code_group==0) code_group=28 ! T
ELSEIF (nrot==24) THEN
IF (noperation(5)>0) code_group=30 ! T_d
IF (noperation(5)==0) code_group=31 ! O
ELSE
CALL errore('find_group','wrong number of elements',1)
ENDIF
ELSEIF (noperation(2)==1) THEN
!
! There is inversion
!
IF (nrot==2) THEN
code_group=2 ! C_i
ELSEIF (nrot==4) THEN
code_group=16 ! C_2h
ELSEIF (nrot==6) THEN
code_group=27 ! S_6
ELSEIF (nrot==8) THEN
IF (noperation(5)==3) code_group=20 ! D_2h
IF (noperation(5)==1) code_group=18 ! C_4h
ELSEIF (nrot==12) THEN
IF (noperation(5)==3) code_group=25 ! D_3d
IF (noperation(5)==1) code_group=19 ! C_6h
ELSEIF (nrot==16) THEN
IF (noperation(5)==5) code_group=22 ! D_4h
ELSEIF (nrot==24) THEN
IF (noperation(5)>6) code_group=23 ! D_6h
IF (noperation(5)==3) code_group=29 ! T_h
ELSEIF (nrot==48) THEN
code_group=32 ! O_h
ELSE
CALL errore('find_group','wrong number of elements',1)
ENDIF
ELSE
CALL errore('find_group','too many inversions',1)
ENDIF
IF (code_group==0) call errore('find_group','incompatible operations',1)
gname=group_name(code_group)
RETURN
!---------------------------------------------------------------------------
SUBROUTINE find_group( nrot, smat, gname, code_group )
!--------------------------------------------------------------------------
!! Given a group of nrot rotation matrices smat (in cartesian coordinates)
!! this routine finds the name of the point group. It assumes but does not
!! check that:
!! * the nrot matrices smat are actually a group;
!! * the group is one of the thirty-two point groups.
!
USE kinds, ONLY : DP
!
IMPLICIT NONE
!
INTEGER :: nrot
!! number of rotation matrices
INTEGER :: code_group
!! code that identifies the group
REAL(DP) :: smat(3,3,nrot)
!! rotation matrices in cartesian coordinates
CHARACTER (LEN=11) :: gname
!! name of the group
!
! ... local variables
!
CHARACTER(LEN=11) :: group_name
INTEGER :: noperation(6), irot, ts, tipo_sym
!
! For each possible group operation the function tipo_sym gives a code
! 1 identity,
! 2 inversion,
! 3 proper rotation <> 180,
! 4 proper rotation 180 degrees,
! 5 mirror,
! 6 improper rotation
! the variable noperation counts how many operations are present in the group.
!
noperation=0
DO irot = 1, nrot
ts = tipo_sym(smat(1,1,irot))
noperation(ts) = noperation(ts)+1
ENDDO
!
IF (noperation(1) /= 1) CALL errore( 'find_group', 'the group has not identity', 1 )
!
code_group=0
!
IF (noperation(2)==0) THEN
!
! There is not inversion
!
SELECT CASE( nrot )
CASE( 1 )
code_group=1 ! C_1
CASE( 2 )
IF (noperation(4)==1) code_group=4 ! C_2
IF (noperation(5)==1) code_group=3 ! C_s
CASE( 3 )
IF (noperation(3)==2) code_group=5 ! C_3
CASE( 4 )
IF (noperation(6)>0) code_group=26 ! S_4
IF (noperation(5)>0.AND.code_group==0) code_group=12 ! C_2v
IF (noperation(3)>0.AND.code_group==0) code_group=6 ! C_4
IF (noperation(4)>0.AND.code_group==0) code_group=8 ! D_2
CASE( 6 )
IF (noperation(5)==3) code_group=13 ! C_3v
IF (noperation(5)==1) code_group=17 ! C_3h
IF (noperation(4)==3.AND.code_group==0) code_group=9 ! D_3
IF (noperation(3)> 0.AND.code_group==0) code_group=7 ! C_6
CASE( 8 )
IF (noperation(5)==4) code_group=14 ! C_4v
IF (noperation(5)==2) code_group=24 ! D_2d
IF (noperation(3)>0.AND.code_group==0) code_group=10 ! D_4
CASE( 12 )
IF (noperation(5)==6) code_group=15 ! C_6v
IF (noperation(5)==4) code_group=21 ! D_3h
IF (noperation(4)>6.AND.code_group==0) code_group=11 ! D_6
IF (noperation(3)>0.AND.code_group==0) code_group=28 ! T
CASE( 24 )
IF (noperation(5)>0) code_group=30 ! T_d
IF (noperation(5)==0) code_group=31 ! O
CASE DEFAULT
CALL errore( 'find_group','wrong number of elements', 1 )
END SELECT
!
ELSEIF (noperation(2)==1) THEN
!
! There is inversion
!
SELECT CASE( nrot )
CASE( 2 )
code_group=2 ! C_i
CASE( 4 )
code_group=16 ! C_2h
CASE( 6 )
code_group=27 ! S_6
CASE( 8 )
IF (noperation(5)==3) code_group=20 ! D_2h
IF (noperation(5)==1) code_group=18 ! C_4h
CASE( 12 )
IF (noperation(5)==3) code_group=25 ! D_3d
IF (noperation(5)==1) code_group=19 ! C_6h
CASE( 16 )
IF (noperation(5)==5) code_group=22 ! D_4h
CASE( 24 )
IF (noperation(5)>6) code_group=23 ! D_6h
IF (noperation(5)==3) code_group=29 ! T_h
CASE( 48 )
code_group=32 ! O_h
CASE DEFAULT
CALL errore( 'find_group', 'wrong number of elements', 1 )
END SELECT
ELSE
CALL errore( 'find_group', 'too many inversions', 1 )
ENDIF
!
IF (code_group==0) CALL errore( 'find_group', 'incompatible operations', 1 )
!
gname = group_name( code_group )
!
RETURN
!
END SUBROUTINE find_group
!--------------------------------------------------------------------------
FUNCTION group_name(code)
!--------------------------------------------------------------------------
! This function receives a code of the group and provides the name of the
! group. The order is the following:
!
! 1 "C_1 " 11 "D_6 " 21 "D_3h" 31 "O "
! 2 "C_i " 12 "C_2v" 22 "D_4h" 32 "O_h "
! 3 "C_s " 13 "C_3v" 23 "D_6h"
! 4 "C_2 " 14 "C_4v" 24 "D_2d"
! 5 "C_3 " 15 "C_6v" 25 "D_3d"
! 6 "C_4 " 16 "C_2h" 26 "S_4 "
! 7 "C_6 " 17 "C_3h" 27 "S_6 "
! 8 "D_2 " 18 "C_4h" 28 "T "
! 9 "D_3 " 19 "C_6h" 29 "T_h "
! 10 "D_4 " 20 "D_2h" 30 "T_d "
!
IMPLICIT NONE
INTEGER :: code
CHARACTER(LEN=11) :: group_name
CHARACTER(LEN=11) :: gname(32)
data gname / "C_1 (1) ", "C_i (-1) ", "C_s (m) ", "C_2 (2) ", &
"C_3 (3) ", "C_4 (4) ", "C_6 (6) ", "D_2 (222) ", &
"D_3 (32) ", "D_4 (422) ", "D_6 (622) ", "C_2v (mm2) ", &
"C_3v (3m) ", "C_4v (4mm) ", "C_6v (6mm) ", "C_2h (2/m) ", &
"C_3h (-6) ", "C_4h (4/m) ", "C_6h (6/m) ", "D_2h (mmm) ", &
"D_3h (-62m)", "D_4h(4/mmm)", "D_6h(6/mmm)", "D_2d (-42m)", &
"D_3d (-3m) ", "S_4 (-4) ", "S_6 (-3) ", "T (23) ", &
"T_h (m-3) ", "T_d (-43m) ", "O (432) ", "O_h (m-3m) " /
IF (code < 1 .OR. code > 32 ) CALL errore('group_name','code is out of range',1)
group_name=gname(code)
RETURN
!--------------------------------------------------------------------------
FUNCTION group_name( code )
!--------------------------------------------------------------------------
!! This function receives a code of the group and provides the name of the
!! group. The order is the following:
!
!! \begin{array}{|r|r|r|r|r|r|}
!! \hline
!! 1 & C_1 & 12 & C_2v & 23 & D_6h \\
!! \hline
!! 2 & C_i & 13 & C_3v & 24 & D_2d \\
!! \hline
!! 3 & C_s & 14 & C_4v & 25 & D_3d \\
!! \hline
!! 4 & C_2 & 15 & C_6v & 26 & S_4 \\
!! \hline
!! 5 & C_3 & 16 & C_2h & 27 & S_6 \\
!! \hline
!! 6 & C_4 & 17 & C_3h & 28 & T \\
!! \hline
!! 7 & C_6 & 18 & C_4h & 29 & T_h \\
!! \hline
!! 8 & D_2 & 19 & C_6h & 30 & T_d \\
!! \hline
!! 9 & D_3 & 20 & D_2h & 31 & O \\
!! \hline
!! 10 & D_4 & 21 & D_3h & 32 & O_h \\
!! \hline
!! 11 & D_6 & 22 & D_4h & & \\
!! \hline
!! \end{array}
!
IMPLICIT NONE
!
INTEGER :: code
!! input: code of the group
CHARACTER(LEN=11) :: group_name
!! output: name of the group
!
! ... local variables
!
CHARACTER(LEN=11) :: gname(32)
DATA gname / "C_1 (1) ", "C_i (-1) ", "C_s (m) ", "C_2 (2) ", &
"C_3 (3) ", "C_4 (4) ", "C_6 (6) ", "D_2 (222) ", &
"D_3 (32) ", "D_4 (422) ", "D_6 (622) ", "C_2v (mm2) ", &
"C_3v (3m) ", "C_4v (4mm) ", "C_6v (6mm) ", "C_2h (2/m) ", &
"C_3h (-6) ", "C_4h (4/m) ", "C_6h (6/m) ", "D_2h (mmm) ", &
"D_3h (-62m)", "D_4h(4/mmm)", "D_6h(6/mmm)", "D_2d (-42m)", &
"D_3d (-3m) ", "S_4 (-4) ", "S_6 (-3) ", "T (23) ", &
"T_h (m-3) ", "T_d (-43m) ", "O (432) ", "O_h (m-3m) " /
!
IF (code < 1 .OR. code > 32 ) CALL errore( 'group_name', 'code is out of range', 1 )
!
group_name = gname(code)
!
RETURN
!
END FUNCTION group_name
!
!
!--------------------------------------------------------------------------
FUNCTION tipo_sym(s)
!--------------------------------------------------------------------------
! This function receives a 3x3 orthogonal matrix which is a symmetry
! operation of the point group of the crystal written in cartesian
! coordinates and gives as output a code according to the following:
!
! 1 Identity
! 2 Inversion
! 3 Proper rotation of an angle <> 180 degrees
! 4 Proper rotation of 180 degrees
! 5 Mirror symmetry
! 6 Improper rotation
!
USE kinds, ONLY : DP
FUNCTION tipo_sym( s )
!--------------------------------------------------------------------------
!! This function receives a 3x3 orthogonal matrix which is a symmetry
!! operation of the point group of the crystal written in cartesian
!! coordinates and gives as output a code according to the following:
!! 1 - identity;
!! 2 - inversion;
!! 3 - proper rotation of an angle <> 180 degrees;
!! 4 - proper rotation of 180 degrees;
!! 5 - mirror symmetry;
!! 6 - improper rotation.
!
USE kinds, ONLY : DP
!
IMPLICIT NONE
!
REAL(DP) :: s(3,3)
!! input: see the function main comments
INTEGER :: tipo_sym
!! output: see the function main comments
!
! ... local variables
!
REAL(DP), PARAMETER :: eps=1.d-7
REAL(DP) :: det, det1
!
! Check for identity
!
IF ((ABS(s(1,1)-1.d0) < eps).AND. &
(ABS(s(2,2)-1.d0) < eps).AND. &
(ABS(s(3,3)-1.d0) < eps).AND. &
(ABS(s(1,2)) < eps).AND.(ABS(s(2,1)) < eps).AND.(ABS(s(2,3)) < eps).AND. &
(ABS(s(3,2)) < eps).AND.(ABS(s(1,3)) < eps).AND.(ABS(s(3,1)) < eps)) THEN
tipo_sym=1
RETURN
ENDIF
!
! Check for inversion
!
IF ((ABS(s(1,1)+1.d0) < eps).AND. &
(ABS(s(2,2)+1.d0) < eps).AND. &
(ABS(s(3,3)+1.d0) < eps).AND. &
(ABS(s(1,2)) < eps).AND.(ABS(s(2,1)) < eps).AND.(ABS(s(2,3)) < eps).AND. &
(ABS(s(3,2)) < eps).AND.(ABS(s(1,3)) < eps).AND.(ABS(s(3,1)) < eps)) THEN
tipo_sym=2
RETURN
ENDIF
!
! compute the determinant
!
det = s(1,1) * ( s(2,2) * s(3,3) - s(3,2) * s(2,3) ) - &
s(1,2) * ( s(2,1) * s(3,3) - s(3,1) * s(2,3) ) + &
s(1,3) * ( s(2,1) * s(3,2) - s(3,1) * s(2,2) )
!
! Determinant equal to 1: proper rotation
!
IF (ABS(det-1.d0) < eps) THEN
!
! check if an eigenvalue is equal to -1.d0 (180 rotation)
!
det1=(s(1,1)+1.d0)*((s(2,2)+1.d0)*(s(3,3)+1.d0)-s(3,2)*s(2,3)) - &
s(1,2)* (s(2,1)* (s(3,3)+1.d0)-s(3,1)*s(2,3)) + &
s(1,3)* (s(2,1)*s(3,2) -s(3,1)*(s(2,2)+1.d0))
IMPLICIT NONE
IF (ABS(det1) < eps) THEN
tipo_sym = 4 ! 180 proper rotation
ELSE
tipo_sym = 3 ! proper rotation <> 180
ENDIF
RETURN
ENDIF
!
! Determinant equal to -1: mirror symmetry or improper rotation
!
IF (ABS(det+1.d0) < eps) THEN
!
! check if an eigenvalue is equal to 1.d0 (mirror symmetry)
!
det1=(s(1,1)-1.d0)*((s(2,2)-1.d0)*(s(3,3)-1.d0)-s(3,2)*s(2,3)) - &
s(1,2)* (s(2,1)* (s(3,3)-1.d0)-s(3,1)*s(2,3)) + &
s(1,3)* (s(2,1)*s(3,2) -s(3,1)*(s(2,2)-1.d0))
REAL(DP), PARAMETER :: eps=1.d-7
REAL(DP) :: s(3,3), det, det1
INTEGER :: tipo_sym
!
! Check for identity
!
IF ((ABS(s(1,1)-1.d0) < eps).AND. &
(ABS(s(2,2)-1.d0) < eps).AND. &
(ABS(s(3,3)-1.d0) < eps).AND. &
(ABS(s(1,2)) < eps).AND.(ABS(s(2,1)) < eps).AND.(ABS(s(2,3)) < eps).AND. &
(ABS(s(3,2)) < eps).AND.(ABS(s(1,3)) < eps).AND.(ABS(s(3,1)) < eps)) THEN
tipo_sym=1
RETURN
ENDIF
!
! Check for inversion
!
IF ((ABS(s(1,1)+1.d0) < eps).AND. &
(ABS(s(2,2)+1.d0) < eps).AND. &
(ABS(s(3,3)+1.d0) < eps).AND. &
(ABS(s(1,2)) < eps).AND.(ABS(s(2,1)) < eps).AND.(ABS(s(2,3)) < eps).AND. &
(ABS(s(3,2)) < eps).AND.(ABS(s(1,3)) < eps).AND.(ABS(s(3,1)) < eps)) THEN
tipo_sym=2
RETURN
ENDIF
!
! compute the determinant
!
det = s(1,1) * ( s(2,2) * s(3,3) - s(3,2) * s(2,3) )- &
s(1,2) * ( s(2,1) * s(3,3) - s(3,1) * s(2,3) )+ &
s(1,3) * ( s(2,1) * s(3,2) - s(3,1) * s(2,2) )
!
! Determinant equal to 1: proper rotation
!
IF (abs(det-1.d0) < eps) THEN
!
! check if an eigenvalue is equal to -1.d0 (180 rotation)
!
det1=(s(1,1)+1.d0)*((s(2,2)+1.d0)*(s(3,3)+1.d0)-s(3,2)*s(2,3))- &
s(1,2)* (s(2,1)* (s(3,3)+1.d0)-s(3,1)*s(2,3))+ &
s(1,3)* (s(2,1)*s(3,2) -s(3,1)*(s(2,2)+1.d0))
IF (abs(det1) < eps) THEN
tipo_sym=4 ! 180 proper rotation
ELSE
tipo_sym=3 ! proper rotation <> 180
ENDIF
RETURN
ENDIF
!
! Determinant equal to -1: mirror symmetry or improper rotation
!
IF (abs(det+1.d0) < eps) THEN
!
! check if an eigenvalue is equal to 1.d0 (mirror symmetry)
!
det1=(s(1,1)-1.d0)*((s(2,2)-1.d0)*(s(3,3)-1.d0)-s(3,2)*s(2,3))- &
s(1,2)* (s(2,1)* (s(3,3)-1.d0)-s(3,1)*s(2,3))+ &
s(1,3)* (s(2,1)*s(3,2) -s(3,1)*(s(2,2)-1.d0))
IF (abs(det1) < eps) THEN
tipo_sym=5 ! mirror symmetry
ELSE
tipo_sym=6 ! improper rotation
ENDIF
RETURN
ELSE
call errore('tipo_sym','symmetry not recognized',1)
ENDIF
IF (ABS(det1) < eps) THEN
tipo_sym=5 ! mirror symmetry
ELSE
tipo_sym=6 ! improper rotation
ENDIF
!
RETURN
!
ELSE
CALL errore( 'tipo_sym', 'symmetry not recognized', 1 )
ENDIF
!
END FUNCTION tipo_sym
!
!--------------------------------------------------------------------------
FUNCTION laue_class(code)
FUNCTION laue_class( code )
!--------------------------------------------------------------------------
! This function receives a code of the point group and provides the
! code of the point group that defines the Laue class (that is the point
! group obtained by multipling by inversion).
! The order is the following:
!
! 1 "C_1 " -> 2 11 "D_6 " -> 23 21 "D_3h" -> 23 31 "O " -> 32
! 2 "C_i " -> 2 12 "C_2v" -> 20 22 "D_4h" -> 22 32 "O_h " -> 32
! 3 "C_s " -> 16 13 "C_3v" -> 25 23 "D_6h" -> 23
! 4 "C_2 " -> 16 14 "C_4v" -> 22 24 "D_2d" -> 22
! 5 "C_3 " -> 27 15 "C_6v" -> 23 25 "D_3d" -> 25
! 6 "C_4 " -> 18 16 "C_2h" -> 16 26 "S_4 " -> 18
! 7 "C_6 " -> 19 17 "C_3h" -> 19 27 "S_6 " -> 27
! 8 "D_2 " -> 20 18 "C_4h" -> 18 28 "T " -> 29
! 9 "D_3 " -> 25 19 "C_6h" -> 19 29 "T_h " -> 29
! 10 "D_4 " -> 22 20 "D_2h" -> 20 30 "T_d " -> 32
!
IMPLICIT NONE
INTEGER :: code
INTEGER :: laue_class
INTEGER :: laue(32)
DATA laue / 2, 2, 16, 16, 27, 18, 19, 20, 25, 22, &
23, 20, 25, 22, 23, 16, 19, 18, 19, 20, &
23, 22, 23, 22, 25, 18, 27, 29, 29, 32, &
32, 32 /
IF (code < 1 .OR. code > 32 ) CALL errore('laue_class','code is out of range',1)
laue_class=laue(code)
RETURN
!! This function receives a code of the point group and provides the
!! code of the point group that defines the Laue class (that is the point
!! group obtained by multipling by inversion).
!! The order is the following:
!
!! \begin{array}{|r|c|r|c|r|c|}
!! \hline
!! 1 & C_1\rightarrow 2 & 12 & C_2v\rightarrow 20 & 23 & D_6h\rightarrow 23 \\
!! \hline
!! 2 & C_i\rightarrow 2 & 13 & C_3v\rightarrow 25 & 24 & D_2d\rightarrow 22 \\
!! \hline
!! 3 & C_s\rightarrow 16 & 14 & C_4v\rightarrow 22 & 25 & D_3d\rightarrow 25 \\
!! \hline
!! 4 & C_2\rightarrow 16 & 15 & C_6v\rightarrow 23 & 26 & S_4\rightarrow 18 \\
!! \hline
!! 5 & C_3\rightarrow 27 & 16 & C_2h\rightarrow 16 & 27 & S_6\rightarrow 27 \\
!! \hline
!! 6 & C_4\rightarrow 18 & 17 & C_3h\rightarrow 19 & 28 & T\rightarrow 29 \\
!! \hline
!! 7 & C_6\rightarrow 19 & 18 & C_4h\rightarrow 18 & 29 & T_h\rightarrow 29 \\
!! \hline
!! 8 & D_2\rightarrow 20 & 19 & C_6h\rightarrow 19 & 30 & T_d\rightarrow 32 \\
!! \hline
!! 9 & D_3\rightarrow 25 & 20 & D_2h\rightarrow 20 & 31 & O\rightarrow 32 \\
!! \hline
!! 10 & D_4\rightarrow 22 & 21 & D_3h\rightarrow 23 & 32 & O_h\rightarrow 32 \\
!! \hline
!! 11 & D_6\rightarrow 23 & 22 & D_4h\rightarrow 22 & & \\
!! \hline
!! \end{array}
!
IMPLICIT NONE
!
INTEGER :: code
!! code of the point group
INTEGER :: laue_class
!! code of the point group that defines the Laue class
!
! ... local variables
!
INTEGER :: laue(32)
!
DATA laue / 2, 2, 16, 16, 27, 18, 19, 20, 25, 22, &
23, 20, 25, 22, 23, 16, 19, 18, 19, 20, &
23, 22, 23, 22, 25, 18, 27, 29, 29, 32, &
32, 32 /
!
IF (code<1 .OR. code>32) CALL errore( 'laue_class', 'code is out of range', 1 )
!
laue_class = laue(code)
!
RETURN
!
END FUNCTION laue_class