mirror of https://gitlab.com/QEF/q-e.git
1168 lines
37 KiB
Fortran
1168 lines
37 KiB
Fortran
!
|
|
! Copyright (C) 2002-2011 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 cell_base
|
|
!------------------------------------------------------------------------------!
|
|
!! Cell parameters and initialization.
|
|
!
|
|
USE kinds, ONLY : DP
|
|
USE constants, ONLY : pi, bohr_radius_angs
|
|
USE io_global, ONLY : stdout
|
|
!
|
|
IMPLICIT NONE
|
|
SAVE
|
|
!
|
|
INTEGER :: ibrav
|
|
!! Index of the bravais lattice (see latgen.f90)
|
|
REAL(DP) :: celldm(6) = (/ 0.0_DP,0.0_DP,0.0_DP,0.0_DP,0.0_DP,0.0_DP /)
|
|
!! Old-style parameters of the simulation cell (se latgen.f90)
|
|
REAL(DP) :: a
|
|
!! Traditional crystallographic cell parameters (alpha=cosbc and so on)
|
|
REAL(DP) :: b, c, cosab, cosac, cosbc
|
|
!
|
|
CHARACTER(len=80) :: cell_units
|
|
!! Format of input cell parameters: 'alat','bohr','angstrom'
|
|
REAL(DP) :: alat = 0.0_DP
|
|
!! Lattice parameter - often used to scale quantities, or in
|
|
!! combination to other parameters/constants to define new units.
|
|
REAl(DP) :: omega = 0.0_DP
|
|
!! Volume of the simulation cell.
|
|
REAL(DP) :: tpiba = 0.0_DP
|
|
!! 2 PI/alat
|
|
REAL(DP) :: tpiba2 = 0.0_DP
|
|
!! tpiba2=tpiba^2
|
|
!
|
|
! direct and reciprocal lattice primitive vectors
|
|
REAL(DP) :: at(3,3) = RESHAPE( (/ 0.0_DP /), (/ 3, 3 /), (/ 0.0_DP /) )
|
|
!! The lattice vectors of the simulation cell, a_i, in alat units:
|
|
!! a_i(:) = at(:,i)/alat
|
|
REAL(DP) :: bg(3,3) = RESHAPE( (/ 0.0_DP /), (/ 3, 3 /), (/ 0.0_DP /) )
|
|
!! The reciprocal lattice vectors, b_i, in tpiba=2pi/alat units:
|
|
!! b_i(:) = bg(:,i)/tpiba
|
|
!
|
|
! parameters for reference cell
|
|
REAL(DP) :: ref_tpiba2 = 0.0_DP
|
|
!! Parameters for reference cell
|
|
REAL(DP) :: ref_at(3,3) = RESHAPE( (/ 0.0_DP /), (/ 3, 3 /), (/ 0.0_DP /) )
|
|
REAL(DP) :: ref_bg(3,3) = RESHAPE( (/ 0.0_DP /), (/ 3, 3 /), (/ 0.0_DP /) )
|
|
!
|
|
REAL(DP) :: init_tpiba2 = 0.0_DP
|
|
!! Parameter to store tpiba2 calculated from the input cell parameter
|
|
!! used in emass_preconditioning, required for restarting variable cell calculation
|
|
!! correctly in CP.
|
|
!
|
|
! -------------------------------------------------------------------------!
|
|
!
|
|
TYPE boxdimensions
|
|
!! Periodicity box: in the matrix "a" every row is the vector
|
|
!! of each side of the cell in the real space
|
|
REAL(DP) :: a(3,3)
|
|
!! direct lattice generators
|
|
REAL(DP) :: m1(3,3)
|
|
!! reciprocal lattice generators
|
|
REAL(DP) :: omega
|
|
!! cell volume = determinant of a
|
|
REAL(DP) :: g(3,3)
|
|
!! metric tensor
|
|
REAL(DP) :: gvel(3,3)
|
|
!! metric velocity
|
|
REAL(DP) :: pail(3,3)
|
|
!! stress tensor ( scaled coor. )
|
|
REAL(DP) :: paiu(3,3)
|
|
!! stress tensor ( cartesian coor. )
|
|
REAL(DP) :: hmat(3,3)
|
|
!! cell parameters ( transpose of "a" )
|
|
REAL(DP) :: hvel(3,3)
|
|
!! cell velocity
|
|
REAL(DP) :: hinv(3,3)
|
|
REAL(DP) :: deth
|
|
INTEGER :: perd(3)
|
|
END TYPE boxdimensions
|
|
!
|
|
! The following relations should always be kept valid:
|
|
! h = at*alat; ainv = h^(-1); ht=transpose(h)
|
|
!
|
|
REAL(DP) :: h(3,3) = 0.0_DP
|
|
!! simulation cell at time t
|
|
REAL(DP) :: ainv(3,3) = 0.0_DP
|
|
REAL(DP) :: hold(3,3) = 0.0_DP
|
|
!! simulation cell at time t-delt
|
|
REAL(DP) :: hnew(3,3) = 0.0_DP
|
|
!! simulation cell at time t+delt
|
|
REAL(DP) :: velh(3,3) = 0.0_DP
|
|
!! simulation cell velocity
|
|
REAL(DP) :: deth = 0.0_DP
|
|
!! determinant of h ( cell volume )
|
|
!
|
|
INTEGER :: iforceh(3,3) = 1
|
|
!! if iforceh( i, j ) = 0 then h( i, j ) is not allowed to move
|
|
LOGICAL :: enforce_ibrav = .FALSE.
|
|
!! True if ibrav representation is fix
|
|
LOGICAL :: fix_volume = .FALSE.
|
|
!! True if cell volume is kept fixed
|
|
LOGICAL :: fix_area = .FALSE.
|
|
!! True if area in xy plane is kept constant
|
|
LOGICAL :: isotropic = .FALSE.
|
|
!! True if volume option is chosen for cell_dofree
|
|
REAL(DP) :: wmass = 0.0_DP
|
|
!! cell fictitious mass
|
|
REAL(DP) :: press = 0.0_DP
|
|
!! external pressure
|
|
!
|
|
REAL(DP) :: frich = 0.0_DP
|
|
!! friction parameter for cell damped dynamics
|
|
REAL(DP) :: greash = 1.0_DP
|
|
!! greas parameter for damped dynamics
|
|
!
|
|
LOGICAL :: tcell_base_init = .FALSE.
|
|
!
|
|
INTERFACE cell_init
|
|
MODULE PROCEDURE cell_init_ht, cell_init_a
|
|
END INTERFACE
|
|
!
|
|
INTERFACE pbcs
|
|
MODULE PROCEDURE pbcs_components, pbcs_vectors
|
|
END INTERFACE
|
|
!
|
|
INTERFACE s_to_r
|
|
MODULE PROCEDURE s_to_r1, s_to_r1b, s_to_r3
|
|
END INTERFACE
|
|
!
|
|
INTERFACE r_to_s
|
|
MODULE PROCEDURE r_to_s1, r_to_s1b, r_to_s3
|
|
END INTERFACE
|
|
!
|
|
CONTAINS
|
|
!
|
|
!------------------------------------------------------------------------------!
|
|
SUBROUTINE cell_base_init( ibrav_, celldm_, a_, b_, c_, cosab_, cosac_, &
|
|
cosbc_, trd_ht, rd_ht, cell_units_ )
|
|
!-----------------------------------------------------------------------------!
|
|
!! Initialize cell_base module variables, set up crystal lattice.
|
|
!
|
|
IMPLICIT NONE
|
|
INTEGER, INTENT(IN) :: ibrav_
|
|
REAL(DP), INTENT(IN) :: celldm_ (6)
|
|
LOGICAL, INTENT(IN) :: trd_ht
|
|
REAL(DP), INTENT(IN) :: rd_ht (3,3)
|
|
CHARACTER(LEN=*), INTENT(IN) :: cell_units_
|
|
REAL(DP), INTENT(IN) :: a_ , b_ , c_ , cosab_, cosac_, cosbc_
|
|
!
|
|
REAL(DP) :: units
|
|
!
|
|
IF ( ibrav_ == 0 .and. .not. trd_ht ) THEN
|
|
CALL errore('cell_base_init', 'ibrav=0: must read cell parameters', 1)
|
|
ELSE IF ( ibrav_ /= 0 .and. trd_ht ) THEN
|
|
CALL errore('cell_base_init', 'redundant data for cell parameters', 2)
|
|
END IF
|
|
!
|
|
ibrav = ibrav_
|
|
celldm = celldm_
|
|
a = a_ ; b = b_ ; c = c_ ; cosab = cosab_ ; cosac = cosac_ ; cosbc = cosbc_
|
|
cell_units = cell_units_
|
|
units = 0.0_DP
|
|
!
|
|
IF ( trd_ht ) THEN
|
|
!
|
|
! ... crystal lattice vectors read from input: find units
|
|
!
|
|
SELECT CASE ( TRIM( cell_units ) )
|
|
CASE ( 'bohr' )
|
|
IF( celldm( 1 ) /= 0.0_DP .OR. a /= 0.0_dp ) CALL errore &
|
|
('cell_base_init','lattice parameter specified twice',1)
|
|
units = 1.0_DP
|
|
CASE ( 'angstrom' )
|
|
IF( celldm( 1 ) /= 0.0_DP .OR. a /= 0.0_dp ) CALL errore &
|
|
('cell_base_init','lattice parameter specified twice',2)
|
|
units = 1.0_DP / bohr_radius_angs
|
|
CASE ( 'alat' )
|
|
IF( celldm( 1 ) /= 0.0_DP ) THEN
|
|
units = celldm( 1 )
|
|
ELSE IF ( a /= 0.0_dp ) THEN
|
|
units = a / bohr_radius_angs
|
|
ELSE
|
|
CALL errore ('cell_base_init', &
|
|
'lattice parameter not specified',1)
|
|
END IF
|
|
! following case is deprecated and should be removed
|
|
CASE ( 'none' )
|
|
! cell_units is 'none' if nothing was specified
|
|
IF( celldm( 1 ) /= 0.0_DP ) THEN
|
|
units = celldm( 1 )
|
|
cell_units = 'alat'
|
|
ELSE IF ( a /= 0.0_dp ) THEN
|
|
units = a / bohr_radius_angs
|
|
cell_units = 'alat'
|
|
ELSE
|
|
units = 1.0_DP
|
|
cell_units = 'bohr'
|
|
END IF
|
|
!
|
|
CASE DEFAULT
|
|
CALL errore ('cell_base_init', &
|
|
'unexpected cell_units '//TRIM(cell_units),1)
|
|
END SELECT
|
|
!
|
|
! ... Beware the transpose operation between matrices ht and at!
|
|
!
|
|
at = TRANSPOSE( rd_ht ) * units
|
|
!
|
|
! ... at is in atomic units: find alat, bring at to alat units, find omega
|
|
!
|
|
IF( celldm( 1 ) /= 0.0_DP ) THEN
|
|
alat = celldm( 1 )
|
|
ELSE IF ( a /= 0.0_dp ) THEN
|
|
alat = a / bohr_radius_angs
|
|
ELSE
|
|
alat = SQRT ( at(1,1)**2+at(2,1)**2+at(3,1)**2 )
|
|
END IF
|
|
! for compatibility: celldm still used in phonon etc
|
|
celldm(1) = alat
|
|
!
|
|
at(:,:) = at(:,:) / alat
|
|
CALL volume( alat, at(1,1), at(1,2), at(1,3), omega )
|
|
!
|
|
ELSE
|
|
!
|
|
! ... crystal lattice vectors via ibrav + celldm parameters
|
|
!
|
|
IF ( celldm(1) == 0.D0 .and. a /= 0.D0 ) THEN
|
|
!
|
|
! ... convert crystallographic parameters into celldm parameters
|
|
!
|
|
CALL abc2celldm ( ibrav, a,b,c,cosab,cosac,cosbc, celldm )
|
|
!
|
|
ELSE IF ( celldm(1) /= 0.D0 .and. a /= 0.D0 ) THEN
|
|
!
|
|
CALL errore( 'input', 'do not specify both celldm and a,b,c!', 1 )
|
|
!
|
|
END IF
|
|
!
|
|
! ... generate at (in atomic units) from ibrav and celldm
|
|
!
|
|
CALL latgen( ibrav, celldm, at(1,1), at(1,2), at(1,3), omega )
|
|
!
|
|
! ... define lattice constants alat, divide at by alat
|
|
!
|
|
alat = celldm(1)
|
|
at(:,:) = at(:,:) / alat
|
|
!
|
|
END IF
|
|
IF ( alat < 1.9_dp ) CALL infomsg ('cell_base_init', &
|
|
'DEPRECATED: use true lattice parameter, not A to a.u. conversion factor')
|
|
!
|
|
! ... Generate the reciprocal lattice vectors
|
|
!
|
|
CALL recips( at(1,1), at(1,2), at(1,3), bg(1,1), bg(1,2), bg(1,3) )
|
|
!
|
|
tpiba = 2.0_DP * pi / alat
|
|
tpiba2 = tpiba * tpiba
|
|
init_tpiba2 = tpiba2 ! BS : this is used in CPV/src/init_run.f90
|
|
RETURN
|
|
!
|
|
END SUBROUTINE cell_base_init
|
|
!
|
|
!------------------------------------------------------------------------!
|
|
SUBROUTINE ref_cell_base_init( ref_alat, rd_ref_ht, ref_cell_units )
|
|
!--------------------------------------------------------------------!
|
|
!! Initialize cell_base module variables, set up crystal lattice
|
|
!
|
|
IMPLICIT NONE
|
|
REAL(DP), INTENT(IN) :: rd_ref_ht (3,3)
|
|
REAL(DP), INTENT(INOUT) :: ref_alat
|
|
CHARACTER(LEN=*), INTENT(IN) :: ref_cell_units
|
|
!
|
|
REAL(DP) :: units, ref_omega
|
|
!
|
|
! ... reference cell lattice vectors read from REF_CELL_PARAMETERS Card: find units
|
|
!
|
|
SELECT CASE ( TRIM( ref_cell_units ) )
|
|
!
|
|
CASE ( 'bohr' )
|
|
units = 1.0_DP
|
|
CASE ( 'angstrom' )
|
|
units = 1.0_DP / bohr_radius_angs
|
|
CASE DEFAULT
|
|
IF( ref_alat .GT. 0.0_DP ) THEN
|
|
units = ref_alat
|
|
ELSE
|
|
CALL errore('ref_cell_base_init', 'ref_alat must be set to a positive value (in A.U.) in SYSTEM namelist', 1)
|
|
END IF
|
|
!
|
|
END SELECT
|
|
!
|
|
! ... Beware the transpose operation between matrices ht and at!
|
|
!
|
|
ref_at = TRANSPOSE( rd_ref_ht ) * units
|
|
!
|
|
! ... ref_at is in atomic units: find ref_alat, bring ref_at to ref_alat units
|
|
!
|
|
ref_alat = SQRT ( ref_at(1,1)**2+ref_at(2,1)**2+ref_at(3,1)**2 )
|
|
!
|
|
ref_at(:,:) = ref_at(:,:) / ref_alat
|
|
!
|
|
! ... Generate the reciprocal lattice vectors from the reference cell
|
|
!
|
|
CALL recips( ref_at(1,1), ref_at(1,2), ref_at(1,3), ref_bg(1,1), ref_bg(1,2), ref_bg(1,3) )
|
|
!
|
|
ref_tpiba2 = (2.0_DP * pi / ref_alat)**2
|
|
!
|
|
CALL volume( ref_alat, ref_at(1,1), ref_at(1,2), ref_at(1,3), ref_omega )
|
|
!
|
|
WRITE( stdout, * )
|
|
WRITE( stdout, '(3X,"Reference Cell read from REF_CELL_PARAMETERS Card")' )
|
|
WRITE( stdout, '(3X,"Reference Cell alat =",F14.8,1X,"A.U.")' ) ref_alat
|
|
WRITE( stdout, '(3X,"ref_cell_a1 = ",1X,3f14.8)' ) ref_at(:,1)*ref_alat
|
|
WRITE( stdout, '(3X,"ref_cell_a2 = ",1X,3f14.8)' ) ref_at(:,2)*ref_alat
|
|
WRITE( stdout, '(3X,"ref_cell_a3 = ",1X,3f14.8)' ) ref_at(:,3)*ref_alat
|
|
WRITE( stdout, * )
|
|
WRITE( stdout, '(3X,"ref_cell_b1 = ",1X,3f14.8)' ) ref_bg(:,1)/ref_alat
|
|
WRITE( stdout, '(3X,"ref_cell_b2 = ",1X,3f14.8)' ) ref_bg(:,2)/ref_alat
|
|
WRITE( stdout, '(3X,"ref_cell_b3 = ",1X,3f14.8)' ) ref_bg(:,3)/ref_alat
|
|
WRITE( stdout, '(3X,"Reference Cell Volume",F16.8,1X,"A.U.")' ) ref_omega
|
|
!
|
|
RETURN
|
|
!
|
|
END SUBROUTINE ref_cell_base_init
|
|
!------------------------------------------------------------------------------!
|
|
! ... set box
|
|
! ... box%m1(i,1) == b1(i) COLUMN are B vectors
|
|
! ... box%a(1,i) == a1(i) ROW are A vector
|
|
! ... box%omega == volume
|
|
! ... box%g(i,j) == metric tensor G
|
|
!------------------------------------------------------------------------------!
|
|
|
|
SUBROUTINE cell_init_ht( what, box, hval )
|
|
TYPE (boxdimensions) :: box
|
|
REAL(DP), INTENT(IN) :: hval(3,3)
|
|
CHARACTER, INTENT(IN) :: what
|
|
IF( what == 't' .OR. what == 'T' ) THEN
|
|
! hval == ht
|
|
box%a = hval
|
|
box%hmat = TRANSPOSE( hval )
|
|
ELSE
|
|
! hval == hmat
|
|
box%hmat = hval
|
|
box%a = TRANSPOSE( hval )
|
|
END IF
|
|
CALL gethinv( box )
|
|
box%g = MATMUL( box%a(:,:), box%hmat(:,:) )
|
|
box%gvel = 0.0_DP
|
|
box%hvel = 0.0_DP
|
|
box%pail = 0.0_DP
|
|
box%paiu = 0.0_DP
|
|
RETURN
|
|
END SUBROUTINE cell_init_ht
|
|
|
|
!------------------------------------------------------------------------------!
|
|
|
|
SUBROUTINE cell_init_a( alat, at, box )
|
|
TYPE (boxdimensions) :: box
|
|
REAL(DP), INTENT(IN) :: alat, at(3,3)
|
|
INTEGER :: i
|
|
DO i=1,3
|
|
! this is HT: the rows are the lattice vectors
|
|
box%a(1,i) = at(i,1)*alat
|
|
box%a(2,i) = at(i,2)*alat
|
|
box%a(3,i) = at(i,3)*alat
|
|
! this is H : the column are the lattice vectors
|
|
box%hmat(i,1) = at(i,1)*alat
|
|
box%hmat(i,2) = at(i,2)*alat
|
|
box%hmat(i,3) = at(i,3)*alat
|
|
END DO
|
|
box%pail = 0.0_DP
|
|
box%paiu = 0.0_DP
|
|
box%hvel = 0.0_DP
|
|
CALL gethinv(box)
|
|
box%g = MATMUL( box%a(:,:), box%hmat(:,:) )
|
|
box%gvel = 0.0_DP
|
|
RETURN
|
|
END SUBROUTINE cell_init_a
|
|
|
|
!------------------------------------------------------------------------------!
|
|
|
|
SUBROUTINE r_to_s1 (r,s,box)
|
|
REAL(DP), intent(out) :: S(3)
|
|
REAL(DP), intent(in) :: R(3)
|
|
type (boxdimensions), intent(in) :: box
|
|
integer i,j
|
|
DO I=1,3
|
|
S(I) = 0.0_DP
|
|
DO J=1,3
|
|
S(I) = S(I) + R(J)*box%m1(J,I)
|
|
END DO
|
|
END DO
|
|
RETURN
|
|
END SUBROUTINE r_to_s1
|
|
|
|
!------------------------------------------------------------------------------!
|
|
|
|
SUBROUTINE r_to_s3 ( r, s, nat, hinv )
|
|
REAL(DP), intent(out) :: S(:,:)
|
|
INTEGER, intent(in) :: nat
|
|
REAL(DP), intent(in) :: R(:,:)
|
|
REAL(DP), intent(in) :: hinv(:,:) ! hinv = TRANSPOSE( box%m1 )
|
|
integer :: i, j, ia
|
|
DO ia = 1, nat
|
|
DO i=1,3
|
|
S(i,ia) = 0.0_DP
|
|
DO j=1,3
|
|
S(i,ia) = S(i,ia) + R(j,ia)*hinv(i,j)
|
|
END DO
|
|
END DO
|
|
END DO
|
|
RETURN
|
|
END SUBROUTINE r_to_s3
|
|
|
|
!------------------------------------------------------------------------------!
|
|
|
|
SUBROUTINE r_to_s1b ( r, s, hinv )
|
|
REAL(DP), intent(out) :: S(:)
|
|
REAL(DP), intent(in) :: R(:)
|
|
REAL(DP), intent(in) :: hinv(:,:) ! hinv = TRANSPOSE( box%m1 )
|
|
integer :: i, j
|
|
DO I=1,3
|
|
S(I) = 0.0_DP
|
|
DO J=1,3
|
|
S(I) = S(I) + R(J)*hinv(i,j)
|
|
END DO
|
|
END DO
|
|
RETURN
|
|
END SUBROUTINE r_to_s1b
|
|
|
|
|
|
!------------------------------------------------------------------------------!
|
|
|
|
SUBROUTINE s_to_r1 (S,R,box)
|
|
REAL(DP), intent(in) :: S(3)
|
|
REAL(DP), intent(out) :: R(3)
|
|
type (boxdimensions), intent(in) :: box
|
|
integer i,j
|
|
DO I=1,3
|
|
R(I) = 0.0_DP
|
|
DO J=1,3
|
|
R(I) = R(I) + S(J)*box%a(J,I)
|
|
END DO
|
|
END DO
|
|
RETURN
|
|
END SUBROUTINE s_to_r1
|
|
|
|
!------------------------------------------------------------------------------!
|
|
|
|
SUBROUTINE s_to_r1b (S,R,h)
|
|
REAL(DP), intent(in) :: S(3)
|
|
REAL(DP), intent(out) :: R(3)
|
|
REAL(DP), intent(in) :: h(:,:) ! h = TRANSPOSE( box%a )
|
|
integer i,j
|
|
DO I=1,3
|
|
R(I) = 0.0_DP
|
|
DO J=1,3
|
|
R(I) = R(I) + S(J)*h(I,j)
|
|
END DO
|
|
END DO
|
|
RETURN
|
|
END SUBROUTINE s_to_r1b
|
|
|
|
!------------------------------------------------------------------------------!
|
|
|
|
SUBROUTINE s_to_r3 ( S, R, nat, h )
|
|
REAL(DP), intent(in) :: S(:,:)
|
|
INTEGER, intent(in) :: nat
|
|
REAL(DP), intent(out) :: R(:,:)
|
|
REAL(DP), intent(in) :: h(:,:) ! h = TRANSPOSE( box%a )
|
|
integer :: i, j, ia
|
|
DO ia = 1, nat
|
|
DO I = 1, 3
|
|
R(I,ia) = 0.0_DP
|
|
DO J = 1, 3
|
|
R(I,ia) = R(I,ia) + S(J,ia) * h(I,j)
|
|
END DO
|
|
END DO
|
|
END DO
|
|
RETURN
|
|
END SUBROUTINE s_to_r3
|
|
|
|
|
|
!
|
|
!------------------------------------------------------------------------------!
|
|
!
|
|
|
|
SUBROUTINE gethinv(box)
|
|
USE matrix_inversion
|
|
IMPLICIT NONE
|
|
TYPE (boxdimensions), INTENT (INOUT) :: box
|
|
!
|
|
CALL invmat( 3, box%a, box%m1, box%omega )
|
|
box%deth = box%omega
|
|
box%hinv = TRANSPOSE( box%m1 )
|
|
!
|
|
RETURN
|
|
END SUBROUTINE gethinv
|
|
|
|
|
|
FUNCTION get_volume( hmat )
|
|
IMPLICIT NONE
|
|
REAL(DP) :: get_volume
|
|
REAL(DP) :: hmat( 3, 3 )
|
|
get_volume = hmat(1,1)*(hmat(2,2)*hmat(3,3)-hmat(2,3)*hmat(3,2)) + &
|
|
hmat(1,2)*(hmat(2,3)*hmat(3,1)-hmat(2,1)*hmat(3,3)) + &
|
|
hmat(1,3)*(hmat(2,1)*hmat(3,2)-hmat(2,2)*hmat(3,1))
|
|
RETURN
|
|
END FUNCTION get_volume
|
|
!
|
|
!------------------------------------------------------------------------------!
|
|
!
|
|
FUNCTION pbc(rin,box,nl) RESULT (rout)
|
|
IMPLICIT NONE
|
|
TYPE (boxdimensions) :: box
|
|
REAL (DP) :: rin(3)
|
|
REAL (DP) :: rout(3), s(3)
|
|
INTEGER, OPTIONAL :: nl(3)
|
|
|
|
s = matmul(box%hinv(:,:),rin)
|
|
s = s - box%perd*nint(s)
|
|
rout = matmul(box%hmat(:,:),s)
|
|
IF (present(nl)) THEN
|
|
s = REAL( nl, DP )
|
|
rout = rout + matmul(box%hmat(:,:),s)
|
|
END IF
|
|
END FUNCTION pbc
|
|
!
|
|
!------------------------------------------------------------------------------!
|
|
!
|
|
SUBROUTINE get_cell_param(box,cell,ang)
|
|
IMPLICIT NONE
|
|
TYPE(boxdimensions), INTENT(in) :: box
|
|
REAL(DP), INTENT(out), DIMENSION(3) :: cell
|
|
REAL(DP), INTENT(out), DIMENSION(3), OPTIONAL :: ang
|
|
! This code gets the cell parameters given the h-matrix:
|
|
! a
|
|
cell(1)=sqrt(box%hmat(1,1)*box%hmat(1,1)+box%hmat(2,1)*box%hmat(2,1) &
|
|
+box%hmat(3,1)*box%hmat(3,1))
|
|
! b
|
|
cell(2)=sqrt(box%hmat(1,2)*box%hmat(1,2)+box%hmat(2,2)*box%hmat(2,2) &
|
|
+box%hmat(3,2)*box%hmat(3,2))
|
|
! c
|
|
cell(3)=sqrt(box%hmat(1,3)*box%hmat(1,3)+box%hmat(2,3)*box%hmat(2,3) &
|
|
+box%hmat(3,3)*box%hmat(3,3))
|
|
IF (PRESENT(ang)) THEN
|
|
! gamma
|
|
ang(1)=acos((box%hmat(1,1)*box%hmat(1,2)+ &
|
|
box%hmat(2,1)*box%hmat(2,2) &
|
|
+box%hmat(3,1)*box%hmat(3,2))/(cell(1)*cell(2)))
|
|
! beta
|
|
ang(2)=acos((box%hmat(1,1)*box%hmat(1,3)+ &
|
|
box%hmat(2,1)*box%hmat(2,3) &
|
|
+box%hmat(3,1)*box%hmat(3,3))/(cell(1)*cell(3)))
|
|
! alpha
|
|
ang(3)=acos((box%hmat(1,2)*box%hmat(1,3)+ &
|
|
box%hmat(2,2)*box%hmat(2,3) &
|
|
+box%hmat(3,2)*box%hmat(3,3))/(cell(2)*cell(3)))
|
|
! ang=ang*180.0_DP/pi
|
|
|
|
ENDIF
|
|
END SUBROUTINE get_cell_param
|
|
|
|
!------------------------------------------------------------------------------!
|
|
|
|
SUBROUTINE pbcs_components(x1, y1, z1, x2, y2, z2, m)
|
|
!! This subroutine compute the periodic boundary conditions in the scaled
|
|
!! variables system.
|
|
USE kinds
|
|
INTEGER, INTENT(IN) :: M
|
|
REAL(DP), INTENT(IN) :: X1,Y1,Z1
|
|
REAL(DP), INTENT(OUT) :: X2,Y2,Z2
|
|
REAL(DP) MIC
|
|
MIC = REAL( M, DP )
|
|
X2 = X1 - DNINT(X1/MIC)*MIC
|
|
Y2 = Y1 - DNINT(Y1/MIC)*MIC
|
|
Z2 = Z1 - DNINT(Z1/MIC)*MIC
|
|
RETURN
|
|
END SUBROUTINE pbcs_components
|
|
|
|
!------------------------------------------------------------------------------!
|
|
|
|
SUBROUTINE pbcs_vectors(v, w, m)
|
|
!! This subroutine compute the periodic boundary conditions in the scaled
|
|
!! variables system.
|
|
USE kinds
|
|
INTEGER, INTENT(IN) :: m
|
|
REAL(DP), INTENT(IN) :: v(3)
|
|
REAL(DP), INTENT(OUT) :: w(3)
|
|
REAL(DP) :: MIC
|
|
MIC = REAL( M, DP )
|
|
w(1) = v(1) - DNINT(v(1)/MIC)*MIC
|
|
w(2) = v(2) - DNINT(v(2)/MIC)*MIC
|
|
w(3) = v(3) - DNINT(v(3)/MIC)*MIC
|
|
RETURN
|
|
END SUBROUTINE pbcs_vectors
|
|
|
|
!------------------------------------------------------------------------------!
|
|
|
|
SUBROUTINE set_h_ainv()
|
|
!! CP-PW compatibility: align CP arrays H and ainv to at and bg.
|
|
!
|
|
IMPLICIT NONE
|
|
!
|
|
!write(stdout,*) 'alat=',alat
|
|
!write(stdout,*) 'at=',at
|
|
!write(stdout,*) 'bg=',bg
|
|
!
|
|
h(:,:) = at(:,:)*alat
|
|
!
|
|
ainv(1,:) = bg(:,1)/alat
|
|
ainv(2,:) = bg(:,2)/alat
|
|
ainv(3,:) = bg(:,3)/alat
|
|
!
|
|
END SUBROUTINE set_h_ainv
|
|
|
|
!------------------------------------------------------------------------------!
|
|
|
|
SUBROUTINE cell_dyn_init( trd_ht, rd_ht, wc_ , total_ions_mass , press_ , &
|
|
frich_ , greash_ , cell_dofree )
|
|
|
|
USE constants, ONLY: au_gpa, amu_au
|
|
USE io_global, ONLY: stdout
|
|
|
|
IMPLICIT NONE
|
|
CHARACTER(LEN=*), INTENT(IN) :: cell_dofree
|
|
LOGICAL, INTENT(IN) :: trd_ht
|
|
REAL(DP), INTENT(IN) :: rd_ht (3,3)
|
|
REAL(DP), INTENT(IN) :: wc_ , frich_ , greash_ , total_ions_mass
|
|
REAL(DP), INTENT(IN) :: press_ ! external pressure from input
|
|
! ( in KBar = 0.1 GPa )
|
|
INTEGER :: j
|
|
!
|
|
press = press_ / 10.0_DP ! convert press in KBar to GPa
|
|
press = press / au_gpa ! convert to AU
|
|
! frich = frich_ ! for the time being this is set elsewhere
|
|
greash = greash_
|
|
|
|
WRITE( stdout, 105 )
|
|
WRITE( stdout, 110 ) press_
|
|
105 format(/,3X,'Simulation Cell Parameters (from input)')
|
|
110 format( 3X,'external pressure = ',f15.2,' [KBar]')
|
|
|
|
wmass = wc_
|
|
IF( wmass == 0.0_DP ) THEN
|
|
wmass = 3.0_DP / (4.0_DP * pi**2 ) * total_ions_mass
|
|
wmass = wmass * AMU_AU
|
|
WRITE( stdout,130) wmass
|
|
ELSE
|
|
WRITE( stdout,120) wmass
|
|
END IF
|
|
120 format(3X,'wmass (read from input) = ',f15.2,' [AU]')
|
|
130 format(3X,'wmass (calculated) = ',f15.2,' [AU]')
|
|
|
|
IF( wmass <= 0.0_DP ) &
|
|
CALL errore(' cell_dyn_init',' wmass out of range ',0)
|
|
|
|
IF ( trd_ht ) THEN
|
|
!
|
|
WRITE( stdout, 210 )
|
|
WRITE( stdout, 220 ) ( rd_ht( 1, j ), j = 1, 3 )
|
|
WRITE( stdout, 220 ) ( rd_ht( 2, j ), j = 1, 3 )
|
|
WRITE( stdout, 220 ) ( rd_ht( 3, j ), j = 1, 3 )
|
|
!
|
|
210 format(3X,'initial cell from CELL_PARAMETERS card')
|
|
220 format(3X,3F14.8)
|
|
!
|
|
END IF
|
|
!
|
|
ainv(1,:) = bg(:,1)/alat
|
|
ainv(2,:) = bg(:,2)/alat
|
|
ainv(3,:) = bg(:,3)/alat
|
|
!
|
|
CALL init_dofree ( cell_dofree )
|
|
!
|
|
tcell_base_init = .TRUE.
|
|
|
|
WRITE( stdout, 300 ) ibrav
|
|
WRITE( stdout, 305 ) alat
|
|
WRITE( stdout, 310 ) at(:,1)*alat
|
|
WRITE( stdout, 320 ) at(:,2)*alat
|
|
WRITE( stdout, 330 ) at(:,3)*alat
|
|
WRITE( stdout, * )
|
|
WRITE( stdout, 350 ) bg(:,1)/alat
|
|
WRITE( stdout, 360 ) bg(:,2)/alat
|
|
WRITE( stdout, 370 ) bg(:,3)/alat
|
|
WRITE( stdout, 340 ) omega
|
|
300 FORMAT( 3X, 'ibrav = ',I4)
|
|
305 FORMAT( 3X, 'alat = ',F14.8)
|
|
310 FORMAT( 3X, 'a1 = ',3F14.8)
|
|
320 FORMAT( 3X, 'a2 = ',3F14.8)
|
|
330 FORMAT( 3X, 'a3 = ',3F14.8)
|
|
350 FORMAT( 3X, 'b1 = ',3F14.8)
|
|
360 FORMAT( 3X, 'b2 = ',3F14.8)
|
|
370 FORMAT( 3X, 'b3 = ',3F14.8)
|
|
340 FORMAT( 3X, 'omega = ',F16.8)
|
|
|
|
|
|
RETURN
|
|
END SUBROUTINE cell_dyn_init
|
|
|
|
!------------------------------------------------------------------------------!
|
|
SUBROUTINE init_dofree ( cell_dofree )
|
|
!! Set constraints on cell dynamics/optimization
|
|
|
|
CHARACTER(LEN=*), INTENT(IN) :: cell_dofree
|
|
CHARACTER(LEN=80) :: cell_dofree_
|
|
|
|
IF(cell_dofree(1:5) == 'ibrav') THEN
|
|
iforceh = 1
|
|
enforce_ibrav = .true.
|
|
IF(cell_dofree(6:6)=="+")THEN
|
|
cell_dofree_ = cell_dofree(7:)
|
|
ELSE
|
|
cell_dofree_="default"
|
|
ENDIF
|
|
ELSE
|
|
cell_dofree_ = cell_dofree
|
|
ENDIF
|
|
|
|
SELECT CASE ( TRIM( cell_dofree_ ) )
|
|
|
|
CASE ( 'all', 'default', '' )
|
|
iforceh = 1
|
|
!CASE ( 'ibrav')
|
|
! iforceh = 1
|
|
! enforce_ibrav = .true.
|
|
CASE ( 'shape' )
|
|
iforceh = 1
|
|
fix_volume = .true.
|
|
! 2DSHAPE: CASE FOR SHAPE CHANGE IN xy PLANE WITH CONST AREA
|
|
! contribution from Richard Charles Andrew
|
|
! Physics Department, University of Pretoria
|
|
! South Africa, august 2012.
|
|
CASE ( '2Dshape' )
|
|
iforceh = 1
|
|
iforceh(3,3) = 0
|
|
iforceh(1,3) = 0
|
|
iforceh(3,1) = 0
|
|
iforceh(2,3) = 0
|
|
iforceh(3,2) = 0
|
|
fix_area = .true.
|
|
! 2DSHAPE
|
|
CASE ( 'volume' )
|
|
!CALL errore(' init_dofree ', &
|
|
! ' cell_dofree = '//TRIM(cell_dofree)//' not yet implemented ', 1 )
|
|
IF ( ibrav /= 1 ) THEN
|
|
CALL errore('cell_dofree', 'Isotropic expansion is only allowed for ibrav=1; i.e. for simple cubic', 1)
|
|
END IF
|
|
iforceh = 0
|
|
iforceh(1,1) = 1
|
|
iforceh(2,2) = 1
|
|
iforceh(3,3) = 1
|
|
isotropic = .TRUE.
|
|
CASE ('fixa')
|
|
iforceh = 1
|
|
iforceh(1,1) = 0
|
|
iforceh(2,1) = 0
|
|
iforceh(3,1) = 0
|
|
CASE ('fixb')
|
|
iforceh = 1
|
|
iforceh(1,2) = 0
|
|
iforceh(2,2) = 0
|
|
iforceh(3,2) = 0
|
|
CASE ('fixc')
|
|
iforceh = 1
|
|
iforceh(1,3) = 0
|
|
iforceh(2,3) = 0
|
|
iforceh(3,3) = 0
|
|
CASE ('a')
|
|
iforceh = 1
|
|
iforceh(1,1) = 0
|
|
CASE ('b')
|
|
iforceh = 1
|
|
iforceh(2,2) = 0
|
|
CASE ('c')
|
|
iforceh = 1
|
|
iforceh(3,3) = 0
|
|
CASE ('x')
|
|
iforceh = 0
|
|
iforceh(1,1) = 1
|
|
CASE ('y')
|
|
iforceh = 0
|
|
iforceh(2,2) = 1
|
|
CASE ('z')
|
|
iforceh = 0
|
|
iforceh(3,3) = 1
|
|
CASE ('xy')
|
|
iforceh = 0
|
|
iforceh(1,1) = 1
|
|
iforceh(2,2) = 1
|
|
! ... if you want the entire xy plane to be free, uncomment:
|
|
! iforceh(1,2) = 1
|
|
! iforceh(2,1) = 1
|
|
! 2DSHAPE THE ENTIRE xy PLANE IS FREE
|
|
CASE ('2Dxy')
|
|
iforceh = 0
|
|
iforceh(1,1) = 1
|
|
iforceh(2,2) = 1
|
|
iforceh(1,2) = 1
|
|
iforceh(2,1) = 1
|
|
! 2DSHAPE
|
|
CASE ('xz')
|
|
iforceh = 0
|
|
iforceh(1,1) = 1
|
|
iforceh(3,3) = 1
|
|
CASE ('yz')
|
|
iforceh = 0
|
|
iforceh(2,2) = 1
|
|
iforceh(3,3) = 1
|
|
CASE ('xyz')
|
|
iforceh = 0
|
|
iforceh(1,1) = 1
|
|
iforceh(2,2) = 1
|
|
iforceh(3,3) = 1
|
|
! epitaxial constraints (2 axes fixed, one free)
|
|
! added by ulrich.aschauer@dcb.unibe.ch on 2018-02-02
|
|
CASE ('epitaxial_ab')
|
|
!fix the a and b axis while allowing c to change
|
|
iforceh = 0
|
|
iforceh(1,3) = 1
|
|
iforceh(2,3) = 1
|
|
iforceh(3,3) = 1
|
|
CASE ('epitaxial_ac')
|
|
!fix the a and c axis while allowing b to change
|
|
iforceh = 0
|
|
iforceh(1,2) = 1
|
|
iforceh(2,2) = 1
|
|
iforceh(3,2) = 1
|
|
CASE ('epitaxial_bc')
|
|
!fix the b and c axis while allowing a to change
|
|
iforceh = 0
|
|
iforceh(1,1) = 1
|
|
iforceh(2,1) = 1
|
|
iforceh(3,1) = 1
|
|
CASE DEFAULT
|
|
CALL errore(' init_dofree ',' unknown cell_dofree '//TRIM(cell_dofree), 1 )
|
|
|
|
END SELECT
|
|
END SUBROUTINE init_dofree
|
|
|
|
!------------------------------------------------------------------------------!
|
|
|
|
SUBROUTINE cell_base_reinit( ht )
|
|
|
|
USE control_flags, ONLY: iverbosity
|
|
|
|
IMPLICIT NONE
|
|
REAL(DP), INTENT(IN) :: ht (3,3)
|
|
|
|
INTEGER :: j
|
|
|
|
alat = sqrt( ht(1,1)*ht(1,1) + ht(1,2)*ht(1,2) + ht(1,3)*ht(1,3) )
|
|
tpiba = 2.0_DP * pi / alat
|
|
tpiba2 = tpiba * tpiba
|
|
!
|
|
IF( iverbosity > 2 ) THEN
|
|
WRITE( stdout, 210 )
|
|
WRITE( stdout, 220 ) ( ht( 1, j ), j = 1, 3 )
|
|
WRITE( stdout, 220 ) ( ht( 2, j ), j = 1, 3 )
|
|
WRITE( stdout, 220 ) ( ht( 3, j ), j = 1, 3 )
|
|
END IF
|
|
210 format(3X,'Simulation cell parameters with the new cell:')
|
|
220 format(3X,3F14.8)
|
|
|
|
! matrix "ht" used in CP is the transpose of matrix "at"
|
|
! times the lattice parameter "alat"; matrix "ainv" is "bg" divided alat
|
|
!
|
|
at = TRANSPOSE( ht ) / alat
|
|
!
|
|
CALL recips( at(1,1), at(1,2), at(1,3), bg(1,1), bg(1,2), bg(1,3) )
|
|
CALL volume( alat, at(1,1), at(1,2), at(1,3), deth )
|
|
omega = deth
|
|
!
|
|
ainv(1,:) = bg(:,1)/alat
|
|
ainv(2,:) = bg(:,2)/alat
|
|
ainv(3,:) = bg(:,3)/alat
|
|
!
|
|
IF( iverbosity > 2 ) THEN
|
|
WRITE( stdout, 305 ) alat
|
|
WRITE( stdout, 310 ) at(:,1)*alat
|
|
WRITE( stdout, 320 ) at(:,2)*alat
|
|
WRITE( stdout, 330 ) at(:,3)*alat
|
|
WRITE( stdout, * )
|
|
WRITE( stdout, 350 ) bg(:,1)/alat
|
|
WRITE( stdout, 360 ) bg(:,2)/alat
|
|
WRITE( stdout, 370 ) bg(:,3)/alat
|
|
WRITE( stdout, 340 ) omega
|
|
END IF
|
|
|
|
305 FORMAT( 3X, 'alat = ',F14.8)
|
|
310 FORMAT( 3X, 'a1 = ',3F14.8)
|
|
320 FORMAT( 3X, 'a2 = ',3F14.8)
|
|
330 FORMAT( 3X, 'a3 = ',3F14.8)
|
|
350 FORMAT( 3X, 'b1 = ',3F14.8)
|
|
360 FORMAT( 3X, 'b2 = ',3F14.8)
|
|
370 FORMAT( 3X, 'b3 = ',3F14.8)
|
|
340 FORMAT( 3X, 'omega = ',F14.8)
|
|
|
|
|
|
RETURN
|
|
END SUBROUTINE cell_base_reinit
|
|
|
|
|
|
|
|
!------------------------------------------------------------------------------!
|
|
|
|
SUBROUTINE cell_steepest( hnew, h, delt, iforceh, fcell )
|
|
REAL(DP), INTENT(OUT) :: hnew(3,3)
|
|
REAL(DP), INTENT(IN) :: h(3,3), fcell(3,3)
|
|
INTEGER, INTENT(IN) :: iforceh(3,3)
|
|
REAL(DP), INTENT(IN) :: delt
|
|
INTEGER :: i, j
|
|
REAL(DP) :: dt2,fiso
|
|
dt2 = delt * delt
|
|
!
|
|
IF( isotropic ) THEN
|
|
!
|
|
! Isotropic force on the cell
|
|
!
|
|
fiso = (fcell(1,1)+fcell(2,2)+fcell(3,3))/3.0_DP
|
|
!
|
|
DO j=1,3
|
|
DO i=1,3
|
|
hnew(i,j) = h(i,j) + dt2 * fiso * REAL( iforceh(i,j), DP )
|
|
ENDDO
|
|
ENDDO
|
|
!
|
|
ELSE
|
|
!
|
|
DO j=1,3
|
|
DO i=1,3
|
|
hnew(i,j) = h(i,j) + dt2 * fcell(i,j) * REAL( iforceh(i,j), DP )
|
|
ENDDO
|
|
ENDDO
|
|
!
|
|
END IF
|
|
!
|
|
RETURN
|
|
END SUBROUTINE cell_steepest
|
|
|
|
|
|
!------------------------------------------------------------------------------!
|
|
|
|
SUBROUTINE cell_verlet( hnew, h, hold, delt, iforceh, fcell, frich, tnoseh, hnos )
|
|
REAL(DP), INTENT(OUT) :: hnew(3,3)
|
|
REAL(DP), INTENT(IN) :: h(3,3), hold(3,3), hnos(3,3), fcell(3,3)
|
|
INTEGER, INTENT(IN) :: iforceh(3,3)
|
|
REAL(DP), INTENT(IN) :: frich, delt
|
|
LOGICAL, INTENT(IN) :: tnoseh
|
|
|
|
REAL(DP) :: htmp(3,3)
|
|
REAL(DP) :: verl1, verl2, verl3, dt2, ftmp, v1, v2, v3, fiso
|
|
INTEGER :: i, j
|
|
|
|
dt2 = delt * delt
|
|
|
|
IF( tnoseh ) THEN
|
|
ftmp = 0.0_DP
|
|
htmp = hnos
|
|
ELSE
|
|
ftmp = frich
|
|
htmp = 0.0_DP
|
|
END IF
|
|
|
|
verl1 = 2.0_DP / ( 1.0_DP + ftmp )
|
|
verl2 = 1.0_DP - verl1
|
|
verl3 = dt2 / ( 1.0_DP + ftmp )
|
|
verl1 = verl1 - 1.0_DP
|
|
|
|
IF( isotropic ) THEN
|
|
!
|
|
fiso = (fcell(1,1)+fcell(2,2)+fcell(3,3))/3.0_DP
|
|
!
|
|
DO j=1,3
|
|
DO i=1,3
|
|
v1 = verl1 * h(i,j)
|
|
v2 = verl2 * hold(i,j)
|
|
v3 = verl3 * ( fiso - htmp(i,j) )
|
|
hnew(i,j) = h(i,j) + ( v1 + v2 + v3 ) * REAL( iforceh(i,j), DP )
|
|
ENDDO
|
|
ENDDO
|
|
!
|
|
ELSE
|
|
!
|
|
DO j=1,3
|
|
DO i=1,3
|
|
v1 = verl1 * h(i,j)
|
|
v2 = verl2 * hold(i,j)
|
|
v3 = verl3 * ( fcell(i,j) - htmp(i,j) )
|
|
hnew(i,j) = h(i,j) + ( v1 + v2 + v3 ) * REAL( iforceh(i,j), DP )
|
|
ENDDO
|
|
ENDDO
|
|
!
|
|
END IF
|
|
|
|
RETURN
|
|
END SUBROUTINE cell_verlet
|
|
|
|
!------------------------------------------------------------------------------!
|
|
|
|
subroutine cell_hmove( h, hold, delt, iforceh, fcell )
|
|
REAL(DP), intent(out) :: h(3,3)
|
|
REAL(DP), intent(in) :: hold(3,3), fcell(3,3)
|
|
REAL(DP), intent(in) :: delt
|
|
integer, intent(in) :: iforceh(3,3)
|
|
REAL(DP) :: dt2by2, fac
|
|
integer :: i, j
|
|
dt2by2 = 0.5_DP * delt * delt
|
|
fac = dt2by2
|
|
do i=1,3
|
|
do j=1,3
|
|
h(i,j) = hold(i,j) + fac * iforceh(i,j) * fcell(i,j)
|
|
end do
|
|
end do
|
|
return
|
|
end subroutine cell_hmove
|
|
|
|
!------------------------------------------------------------------------------!
|
|
|
|
subroutine cell_force( fcell, ainv, stress, omega, press, wmassIN )
|
|
USE constants, ONLY : eps8
|
|
REAL(DP), intent(out) :: fcell(3,3)
|
|
REAL(DP), intent(in) :: stress(3,3), ainv(3,3)
|
|
REAL(DP), intent(in) :: omega, press
|
|
REAL(DP), intent(in), optional :: wmassIN
|
|
integer :: i, j
|
|
REAL(DP) :: wmass, fiso
|
|
IF (.not. present(wmassIN)) THEN
|
|
wmass = 1.0
|
|
ELSE
|
|
wmass = wmassIN
|
|
END IF
|
|
do j=1,3
|
|
do i=1,3
|
|
fcell(i,j) = ainv(j,1)*stress(i,1) + ainv(j,2)*stress(i,2) + ainv(j,3)*stress(i,3)
|
|
end do
|
|
end do
|
|
do j=1,3
|
|
do i=1,3
|
|
fcell(i,j) = fcell(i,j) - ainv(j,i) * press
|
|
end do
|
|
end do
|
|
IF( wmass < eps8 ) &
|
|
CALL errore( ' movecell ',' cell mass is less than 0 ! ', 1 )
|
|
fcell = omega * fcell / wmass
|
|
! added this :
|
|
IF( isotropic ) THEN
|
|
!
|
|
! Isotropic force on the cell
|
|
!
|
|
fiso = (fcell(1,1)+fcell(2,2)+fcell(3,3))/3.0_DP
|
|
do i=1,3
|
|
fcell(i,i)=fiso
|
|
end do
|
|
END IF
|
|
!
|
|
return
|
|
end subroutine cell_force
|
|
|
|
!------------------------------------------------------------------------------!
|
|
|
|
subroutine cell_move( hnew, h, hold, delt, iforceh, fcell, frich, tnoseh, vnhh, velh, tsdc )
|
|
REAL(DP), intent(out) :: hnew(3,3)
|
|
REAL(DP), intent(in) :: h(3,3), hold(3,3), fcell(3,3)
|
|
REAL(DP), intent(in) :: vnhh(3,3), velh(3,3)
|
|
integer, intent(in) :: iforceh(3,3)
|
|
REAL(DP), intent(in) :: frich, delt
|
|
logical, intent(in) :: tnoseh, tsdc
|
|
|
|
REAL(DP) :: hnos(3,3)
|
|
|
|
hnew = 0.0
|
|
|
|
if( tnoseh ) then
|
|
hnos = vnhh * velh
|
|
else
|
|
hnos = 0.0_DP
|
|
end if
|
|
!
|
|
IF( tsdc ) THEN
|
|
call cell_steepest( hnew, h, delt, iforceh, fcell )
|
|
ELSE
|
|
call cell_verlet( hnew, h, hold, delt, iforceh, fcell, frich, tnoseh, hnos )
|
|
END IF
|
|
|
|
return
|
|
end subroutine cell_move
|
|
|
|
!------------------------------------------------------------------------------!
|
|
|
|
SUBROUTINE cell_gamma( hgamma, ainv, h, velh )
|
|
!! Compute hgamma = g^-1 * dg/dt that enters in the ions equation of motion.
|
|
!
|
|
IMPLICIT NONE
|
|
REAL(DP), INTENT(OUT) :: hgamma(3,3)
|
|
REAL(DP), INTENT(IN) :: ainv(3,3), h(3,3), velh(3,3)
|
|
REAL(DP) :: gm1(3,3), gdot(3,3)
|
|
!
|
|
! g^-1 inverse of metric tensor = (ht*h)^-1 = ht^-1 * h^-1
|
|
!
|
|
gm1 = MATMUL( ainv, TRANSPOSE( ainv ) )
|
|
!
|
|
! dg/dt = d(ht*h)/dt = dht/dt*h + ht*dh/dt ! derivative of metrix tensor
|
|
!
|
|
gdot = MATMUL( TRANSPOSE( velh ), h ) + MATMUL( TRANSPOSE( h ), velh )
|
|
!
|
|
hgamma = MATMUL( gm1, gdot )
|
|
!
|
|
RETURN
|
|
END SUBROUTINE cell_gamma
|
|
|
|
!------------------------------------------------------------------------------!
|
|
|
|
SUBROUTINE cell_update_vel( htp, ht0, htm, delt, velh )
|
|
!
|
|
IMPLICIT NONE
|
|
TYPE (boxdimensions) :: htp, ht0, htm
|
|
REAL(DP), INTENT(IN) :: delt
|
|
REAL(DP), INTENT(OUT) :: velh( 3, 3 )
|
|
|
|
velh(:,:) = ( htp%hmat(:,:) - htm%hmat(:,:) ) / ( 2.0d0 * delt )
|
|
htp%gvel = ( htp%g(:,:) - htm%g(:,:) ) / ( 2.0d0 * delt )
|
|
ht0%hvel = velh
|
|
|
|
RETURN
|
|
END SUBROUTINE cell_update_vel
|
|
|
|
|
|
!------------------------------------------------------------------------------!
|
|
|
|
subroutine cell_kinene( ekinh, temphh, velh )
|
|
use constants, only: k_boltzmann_au
|
|
implicit none
|
|
REAL(DP), intent(out) :: ekinh, temphh(3,3)
|
|
REAL(DP), intent(in) :: velh(3,3)
|
|
integer :: i,j
|
|
ekinh = 0.0_DP
|
|
do j=1,3
|
|
do i=1,3
|
|
ekinh = ekinh + 0.5_DP*wmass*velh(i,j)*velh(i,j)
|
|
temphh(i,j) = wmass*velh(i,j)*velh(i,j)/k_boltzmann_au
|
|
end do
|
|
end do
|
|
return
|
|
end subroutine cell_kinene
|
|
|
|
!------------------------------------------------------------------------------!
|
|
|
|
function cell_alat( )
|
|
real(DP) :: cell_alat
|
|
if( .NOT. tcell_base_init ) &
|
|
call errore( ' cell_alat ', ' alat has not been set ', 1 )
|
|
cell_alat = alat
|
|
return
|
|
end function cell_alat
|
|
!
|
|
!------------------------------------------------------------------------------!
|
|
END MODULE cell_base
|
|
!------------------------------------------------------------------------------!
|