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

888 lines
28 KiB
Fortran
Raw Permalink Blame History

!
! 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 ions_base
!------------------------------------------------------------------------------!
USE kinds, ONLY : DP
USE parameters, ONLY : ntypx
!
IMPLICIT NONE
SAVE
! nsp = number of species
! na(is) = number of atoms of species is
! nax = max number of atoms of a given species
! nat = total number of atoms of all species
INTEGER :: nsp = 0
INTEGER :: na(ntypx) = 0
INTEGER :: nax = 0
INTEGER :: nat = 0
! zv(is) = (pseudo-)atomic charge
! amass(is) = mass of ions, in atomic mass units
! rcmax(is) = Ewald radius (for ion-ion interactions)
REAL(DP) :: zv(ntypx) = 0.0_DP
REAL(DP) :: amass(ntypx) = 0.0_DP
REAL(DP) :: rcmax(ntypx) = 0.0_DP
! ityp( i ) = the type of i-th atom in stdin
! atm( j ) = name of the type of the j-th atomic specie
! tau( 1:3, i ) = position of the i-th atom
INTEGER, ALLOCATABLE :: ityp(:)
REAL(DP), ALLOCATABLE :: tau(:,:) ! initial positions read from stdin (in bohr)
REAL(DP), ALLOCATABLE :: vel(:,:) ! initial velocities read from stdin (in bohr)
REAL(DP), ALLOCATABLE :: tau_srt(:,:) ! tau sorted by specie in bohr
REAL(DP), ALLOCATABLE :: vel_srt(:,:) ! vel sorted by specie in bohr
INTEGER, ALLOCATABLE :: ind_srt(:) ! index of tau sorted by specie
INTEGER, ALLOCATABLE :: ind_bck(:) ! reverse of ind_srt
CHARACTER(LEN=3) :: atm( ntypx )
CHARACTER(LEN=3), ALLOCATABLE :: label_srt( : )
CHARACTER(LEN=80) :: tau_format ! format of input atomic positions:
! 'alat','crystal','bohr','angstrom'
! if_pos( x, i ) = 0 : x coordinate of i-th atom will be kept fixed
INTEGER, ALLOCATABLE :: if_pos(:,:) ! allowed values: 0 or 1 only
INTEGER, ALLOCATABLE :: iforce(:,:) ! if_pos sorted by specie
INTEGER :: fixatom = 0 ! number of frozen atoms
INTEGER :: ndofp =-1 ! ionic degree of freedom
INTEGER :: ndfrz = 0 ! frozen degrees of freedom
REAL(DP) :: fricp ! friction parameter for damped dynamics
REAL(DP) :: greasp ! friction parameter for damped dynamics
! ... taui = real ionic positions in the center of mass reference
! ... system at istep = 0
! ... this array is used to compute mean square displacements,
! ... it is initialized when NBEG = -1, NBEG = 0 and TAURDR = .TRUE.
! ... first index: x,y,z, second index: atom sorted by specie with respect input
! ... this array is saved in the restart file
REAL(DP), ALLOCATABLE :: taui(:,:)
! ... cdmi = center of mass reference system (related to the taui)
! ... this vector is computed when NBEG = -1, NBEG = 0 and TAURDR = .TRUE.
! ... this array is saved in the restart file
REAL(DP) :: cdmi(3), cdm(3)
! ... cdms = center of mass computed for scaled positions (taus)
REAL(DP) :: cdms(3)
!
REAL(DP), ALLOCATABLE :: extfor(:,:) ! external forces on atoms
LOGICAL :: tions_base_init = .FALSE.
LOGICAL, PRIVATE :: tdebug = .FALSE.
INTERFACE ions_vel
MODULE PROCEDURE ions_vel3, ions_vel2
END INTERFACE
!------------------------------------------------------------------------------!
CONTAINS
!------------------------------------------------------------------------------!
SUBROUTINE sort_tau( tausrt, isrt, tau, isp, nat, nsp )
IMPLICIT NONE
REAL(DP), INTENT(OUT) :: tausrt( :, : )
INTEGER, INTENT(OUT) :: isrt( : )
REAL(DP), INTENT(IN) :: tau( :, : )
INTEGER, INTENT(IN) :: nat, nsp, isp( : )
INTEGER :: ina( nsp ), na( nsp )
INTEGER :: is, ia
! ... count the atoms for each specie
na = 0
DO ia = 1, nat
is = isp( ia )
IF( is < 1 .OR. is > nsp ) &
CALL errore(' sorttau ', ' wrong species index for positions ', ia )
na( is ) = na( is ) + 1
END DO
IF ( ANY ( na(1:nsp) == 0 ) ) &
CALL errore ('sort_atoms', 'some atomic species have no atoms',1)
! ... compute the index of the first atom in each specie
ina( 1 ) = 0
DO is = 2, nsp
ina( is ) = ina( is - 1 ) + na( is - 1 )
END DO
! ... sort the position according to atomic specie
na = 0
DO ia = 1, nat
is = isp( ia )
na( is ) = na( is ) + 1
tausrt( :, na(is) + ina(is) ) = tau(:, ia )
isrt ( na(is) + ina(is) ) = ia
END DO
RETURN
END SUBROUTINE sort_tau
!------------------------------------------------------------------------------!
SUBROUTINE unsort_tau( tau, tausrt, isrt, nat )
IMPLICIT NONE
REAL(DP), INTENT(IN) :: tausrt( :, : )
INTEGER, INTENT(IN) :: isrt( : )
REAL(DP), INTENT(OUT) :: tau( :, : )
INTEGER, INTENT(IN) :: nat
INTEGER :: isa, ia
DO isa = 1, nat
ia = isrt( isa )
tau( :, ia ) = tausrt( :, isa )
END DO
RETURN
END SUBROUTINE unsort_tau
!-------------------------------------------------------------------------
SUBROUTINE ions_base_init( nsp_, nat_, na_, ityp_, tau_, vel_, amass_,&
atm_, if_pos_, tau_format_, alat_, at_, &
rcmax_ , extfor_ )
!-------------------------------------------------------------------------
!
USE constants, ONLY: bohr_radius_angs
USE io_global, ONLY: stdout
!
IMPLICIT NONE
!
INTEGER, INTENT(IN) :: nsp_, nat_, na_(:), ityp_(:)
REAL(DP), INTENT(IN) :: tau_(:,:)
REAL(DP), INTENT(IN) :: vel_(:,:)
REAL(DP), INTENT(IN) :: amass_(:)
CHARACTER(LEN=*), INTENT(IN) :: atm_(:)
CHARACTER(LEN=*), INTENT(IN) :: tau_format_
INTEGER, INTENT(IN) :: if_pos_(:,:)
REAL(DP), INTENT(IN) :: alat_, at_(3,3)
REAL(DP), INTENT(IN) :: rcmax_(:)
REAL(DP), INTENT(IN) :: extfor_(:,:)
!
INTEGER :: i, ia, is
!
!
nsp = nsp_
nat = nat_
!
IF ( nat < 1 ) &
CALL errore( 'ions_base_init ', 'nax out of range', 1 )
IF ( nsp < 1 ) &
CALL errore( 'ions_base_init ', 'nsp out of range', 1 )
IF ( nsp > SIZE( na ) ) &
CALL errore( 'ions_base_init ', &
& 'nsp too large, increase ntypx parameter ', 1 )
!
na(1:nsp) = na_(1:nsp)
nax = MAXVAL( na(1:nsp) )
!
atm(1:nsp) = atm_(1:nsp)
tau_format = TRIM( tau_format_ )
!
IF ( nat /= SUM( na(1:nsp) ) ) &
CALL errore( 'ions_base_init ','inconsistent nat and na ', 1 )
!
CALL deallocate_ions_base()
!
ALLOCATE( ityp( nat ) )
ALLOCATE( tau( 3, nat ) )
ALLOCATE( vel( 3, nat ) )
ALLOCATE( tau_srt( 3, nat ) )
ALLOCATE( vel_srt( 3, nat ) )
ALLOCATE( ind_srt( nat ) )
ALLOCATE( ind_bck( nat ) )
ALLOCATE( if_pos( 3, nat ) )
ALLOCATE( iforce( 3, nat ) )
ALLOCATE( taui( 3, nat ) )
ALLOCATE( label_srt( nat ) )
ALLOCATE( extfor( 3, nat ) )
!
ityp(1:nat) = ityp_(1:nat)
vel(:,1:nat) = vel_(:,1:nat)
if_pos(:,1:nat) = if_pos_(:,1:nat)
!
! ... radii, masses
!
DO is = 1, nsp_
!
rcmax(is) = rcmax_(is)
!
IF( rcmax(is) <= 0.0_DP ) &
CALL errore( 'ions_base_init ', 'invalid rcmax', is )
!
END DO
!
SELECT CASE ( TRIM( tau_format ) )
!
! ... convert input atomic positions to internally used format:
! ... tau in atomic units
!
CASE( 'alat' )
!
! ... input atomic positions are divided by a0
!
tau(:,1:nat) = tau_(:,1:nat) * alat_
vel(:,1:nat) = vel_(:,1:nat) * alat_
!
CASE( 'bohr' )
!
! ... input atomic positions are in a.u.: do nothing
!
tau(:,1:nat) = tau_(:,1:nat)
vel(:,1:nat) = vel_(:,1:nat)
!
CASE( 'crystal' )
!
! ... input atomic positions are in crystal axis ("scaled")
!
DO ia = 1, nat
!
DO i = 1, 3
!
tau(i,ia) = at_(i,1)*alat_ * tau_(1,ia) + &
at_(i,2)*alat_ * tau_(2,ia) + &
at_(i,3)*alat_ * tau_(3,ia)
!
vel(i,ia) = at_(i,1)*alat_ * vel_(1,ia) + &
at_(i,2)*alat_ * vel_(2,ia) + &
at_(i,3)*alat_ * vel_(3,ia)
END DO
!
END DO
!
CASE( 'angstrom' )
!
! ... atomic positions in A
!
tau(:,1:nat) = tau_(:,1:nat) / bohr_radius_angs
vel(:,1:nat) = vel_(:,1:nat) / bohr_radius_angs
!
CASE DEFAULT
!
CALL errore( 'ions_base_init',' tau_format = ' // &
& TRIM( tau_format ) // ' not implemented ', 1 )
!
END SELECT
!
! ... tau_srt : atomic species are ordered according to
! ... the ATOMIC_SPECIES input card. Within each specie atoms are ordered
! ... according to the ATOMIC_POSITIONS input card.
! ... ind_srt : can be used to restore the original position
!
CALL sort_tau( tau_srt, ind_srt, tau, ityp, nat, nsp )
!
vel_srt(:,:) = vel(:,ind_srt(:))
!
DO ia = 1, nat
!
label_srt( ia ) = atm( ityp( ind_srt( ia ) ) )
!
END DO
!
! ... generate ind_bck from ind_srt (reverse sort list)
!
DO ia = 1, nat
!
ind_bck(ind_srt(ia)) = ia
!
END DO
!
DO ia = 1, nat
!
extfor( :, ia ) = extfor_( :, ind_srt( ia ) )
!
END DO
!
IF( tdebug ) THEN
WRITE( stdout, * ) 'ions_base_init: unsorted position and velocities'
DO ia = 1, nat
WRITE( stdout, fmt="(A3,3D12.4,3X,3D12.4)") &
atm( ityp( ia ) ), tau(1:3, ia), vel(1:3,ia)
END DO
WRITE( stdout, * ) 'ions_base_init: sorted position and velocities'
DO ia = 1, nat
WRITE( stdout, fmt="(A3,3D12.4,3X,3D12.4)") &
atm( ityp( ind_srt( ia ) ) ), tau_srt(1:3, ia), vel_srt(1:3,ia)
END DO
END IF
!
! ... The constrain on fixed coordinates is implemented using the array
! ... if_pos whose value is 0 when the coordinate is to be kept fixed, 1
! ... otherwise.
!
if_pos = 1
if_pos(:,:) = if_pos_(:,1:nat)
!
iforce = 0
iforce(:,:) = if_pos(:,ind_srt(:))
!
fixatom=COUNT( if_pos(1,:)==0 .AND. if_pos(2,:)==0 .AND. if_pos(3,:)==0 )
ndofp = COUNT( iforce == 1 )
ndfrz = 3*nat - ndofp
!
amass(1:nsp) = amass_(1:nsp)
!
IF ( ANY( amass(1:nsp) <= 0.0_DP ) ) &
CALL errore( 'ions_base_init ', 'invalid mass', 1 )
!
CALL ions_cofmass( tau_srt, amass, na, nsp, cdmi )
!
DO ia = 1, nat
!
taui(1:3,ia) = tau_srt(1:3,ia) - cdmi(1:3)
!
END DO
!
tions_base_init = .TRUE.
!
RETURN
!
END SUBROUTINE ions_base_init
!
!-------------------------------------------------------------------------
SUBROUTINE deallocate_ions_base()
!-------------------------------------------------------------------------
!
IMPLICIT NONE
!
IF ( ALLOCATED( ityp ) ) DEALLOCATE( ityp )
IF ( ALLOCATED( tau ) ) DEALLOCATE( tau )
IF ( ALLOCATED( vel ) ) DEALLOCATE( vel )
IF ( ALLOCATED( tau_srt ) ) DEALLOCATE( tau_srt )
IF ( ALLOCATED( vel_srt ) ) DEALLOCATE( vel_srt )
IF ( ALLOCATED( ind_srt ) ) DEALLOCATE( ind_srt )
IF ( ALLOCATED( ind_bck ) ) DEALLOCATE( ind_bck )
IF ( ALLOCATED( if_pos ) ) DEALLOCATE( if_pos )
IF ( ALLOCATED( iforce ) ) DEALLOCATE( iforce )
IF ( ALLOCATED( taui ) ) DEALLOCATE( taui )
IF ( ALLOCATED( label_srt ) ) DEALLOCATE( label_srt )
IF ( ALLOCATED( extfor ) ) DEALLOCATE( extfor )
!
tions_base_init = .FALSE.
!
RETURN
!
END SUBROUTINE deallocate_ions_base
!
!-------------------------------------------------------------------------
SUBROUTINE ions_vel3( vel, taup, taum, na, nsp, dt )
!-------------------------------------------------------------------------
USE constants, ONLY : eps8
IMPLICIT NONE
REAL(DP) :: vel(:,:), taup(:,:), taum(:,:)
INTEGER :: na(:), nsp
REAL(DP) :: dt
INTEGER :: ia, is, i, isa
REAL(DP) :: fac
IF( dt < eps8 ) &
CALL errore( ' ions_vel3 ', ' dt <= 0 ', 1 )
fac = 1.0_DP / ( dt * 2.0_DP )
isa = 0
DO is = 1, nsp
DO ia = 1, na(is)
isa = isa + 1
DO i = 1, 3
vel(i,isa) = ( taup(i,isa) - taum(i,isa) ) * fac
END DO
END DO
END DO
RETURN
END SUBROUTINE ions_vel3
!------------------------------------------------------------------------------!
SUBROUTINE ions_vel2( vel, taup, taum, nat, dt )
USE constants, ONLY : eps8
IMPLICIT NONE
REAL(DP) :: vel(:,:), taup(:,:), taum(:,:)
INTEGER :: nat
REAL(DP) :: dt
INTEGER :: ia, i
REAL(DP) :: fac
IF( dt < eps8 ) &
CALL errore( ' ions_vel3 ', ' dt <= 0 ', 1 )
fac = 1.0_DP / ( dt * 2.0_DP )
DO ia = 1, nat
DO i = 1, 3
vel(i,ia) = ( taup(i,ia) - taum(i,ia) ) * fac
END DO
END DO
RETURN
END SUBROUTINE ions_vel2
!------------------------------------------------------------------------------!
SUBROUTINE ions_cofmass( tau, pmass, na, nsp, cdm )
USE constants, ONLY : eps8
IMPLICIT NONE
REAL(DP), INTENT(IN) :: tau(:,:), pmass(:)
REAL(DP), INTENT(OUT) :: cdm(3)
INTEGER, INTENT(IN) :: na(:), nsp
REAL(DP) :: tmas
INTEGER :: is, i, ia, isa
!
tmas=0.0_DP
do is=1,nsp
tmas=tmas+na(is)*pmass(is)
end do
if( tmas < eps8 ) &
call errore(' ions_cofmass ', ' total mass <= 0 ', 1 )
!
do i=1,3
cdm(i)=0.0_DP
isa = 0
do is=1,nsp
do ia=1,na(is)
isa = isa + 1
cdm(i)=cdm(i)+tau(i,isa)*pmass(is)
end do
end do
cdm(i)=cdm(i)/tmas
end do
!
RETURN
END SUBROUTINE ions_cofmass
!------------------------------------------------------------------------------!
SUBROUTINE randpos(tau, na, nsp, tranp, amprp, hinv, ifor )
USE cell_base, ONLY: r_to_s
USE io_global, ONLY: stdout
USE random_numbers, ONLY: randy
IMPLICIT NONE
REAL(DP) :: hinv(3,3)
REAL(DP) :: tau(:,:)
INTEGER, INTENT(IN) :: ifor(:,:), na(:), nsp
LOGICAL, INTENT(IN) :: tranp(:)
REAL(DP), INTENT(IN) :: amprp(:)
REAL(DP) :: oldp(3), rand_disp(3), rdisp(3)
INTEGER :: k, is, isa, isa_s, isa_e, isat
WRITE( stdout, 600 )
isat = 0
DO is = 1, nsp
isa_s = isat + 1
isa_e = isat + na(is)
IF( tranp(is) ) THEN
WRITE( stdout,610) is, na(is)
WRITE( stdout,615)
DO isa = isa_s, isa_e
oldp = tau(:,isa)
rand_disp(1) = randy ()
rand_disp(2) = randy ()
rand_disp(3) = randy ()
rand_disp = amprp(is) * ( rand_disp - 0.5_DP )
rdisp = rand_disp
CALL r_to_s( rdisp(:), rand_disp(:), hinv )
DO k = 1, 3
tau(k,isa) = tau(k,isa) + rand_disp(k) * ifor(k,isa)
END DO
WRITE( stdout,620) (oldp(k),k=1,3), (tau(k,isa),k=1,3)
END DO
END IF
isat = isat + na(is)
END DO
600 FORMAT(//,3X,'Randomization of SCALED ionic coordinates')
610 FORMAT( 3X,'Species ',I3,' atoms = ',I4)
615 FORMAT( 3X,' Old Positions New Positions')
620 FORMAT( 3X,3F10.6,2X,3F10.6)
RETURN
END SUBROUTINE randpos
!------------------------------------------------------------------------------!
SUBROUTINE ions_kinene( ekinp, vels, na, nsp, h, pmass )
IMPLICIT NONE
REAL(DP), intent(out) :: ekinp ! ionic kinetic energy
REAL(DP), intent(in) :: vels(:,:) ! scaled ionic velocities
REAL(DP), intent(in) :: pmass(:) ! ionic masses
REAL(DP), intent(in) :: h(:,:) ! simulation cell
integer, intent(in) :: na(:), nsp
integer :: i, j, is, ia, ii, isa
ekinp = 0.0_DP
isa = 0
do is=1,nsp
do ia=1,na(is)
isa = isa + 1
do j=1,3
do i=1,3
do ii=1,3
ekinp=ekinp+pmass(is)* h(j,i)*vels(i,isa)* h(j,ii)*vels(ii,isa)
end do
end do
end do
end do
end do
ekinp=0.5_DP*ekinp
return
END SUBROUTINE ions_kinene
!------------------------------------------------------------------------------!
subroutine ions_temp( tempp, temps, ekinpr, vels, na, nsp, h, pmass, ndega, nhpdim, atm2nhp, ekin2nhp )
!
use constants, only: k_boltzmann_au
!
implicit none
!
REAL(DP), intent(out) :: ekinpr, tempp
REAL(DP), intent(out) :: temps(:)
REAL(DP), intent(out) :: ekin2nhp(:)
REAL(DP), intent(in) :: vels(:,:)
REAL(DP), intent(in) :: pmass(:)
REAL(DP), intent(in) :: h(:,:)
integer, intent(in) :: na(:), nsp, ndega, nhpdim, atm2nhp(:)
!
integer :: nat, i, j, is, ia, ii, isa
REAL(DP) :: cdmvel(3), eks, eks1
!
call ions_cofmass( vels, pmass, na, nsp, cdmvel )
!
nat = SUM( na(1:nsp) )
!
ekinpr = 0.0_DP
temps( 1:nsp ) = 0.0_DP
ekin2nhp(1:nhpdim) = 0.0_DP
!
do i=1,3
do j=1,3
do ii=1,3
isa = 0
do is=1,nsp
eks = 0.0_DP
do ia=1,na(is)
isa = isa + 1
eks1 = pmass(is)*h(j,i)*(vels(i,isa)-cdmvel(i))*h(j,ii)*(vels(ii,isa)-cdmvel(ii))
eks=eks+eks1
ekin2nhp(atm2nhp(isa)) = ekin2nhp(atm2nhp(isa)) + eks1
end do
ekinpr = ekinpr + eks
temps(is) = temps(is) + eks
end do
end do
end do
end do
!
do is = 1, nhpdim
ekin2nhp(is) = ekin2nhp(is) * 0.5_DP
enddo
!
!
do is = 1, nsp
if( na(is) < 1 ) call errore(' ions_temp ', ' 0 number of atoms ', 1 )
temps( is ) = temps( is ) * 0.5_DP
temps( is ) = temps( is ) / k_boltzmann_au / ( 1.5_DP * na(is) )
end do
!
ekinpr = 0.5_DP * ekinpr
!
IF( ndega < 1 ) THEN
tempp = 0.0_DP
ELSE
tempp = ekinpr / k_boltzmann_au * 2.0_DP / DBLE( ndega )
END IF
!
return
end subroutine ions_temp
!------------------------------------------------------------------------------!
subroutine ions_thermal_stress( stress, nstress, pmass, omega, h, vels, nsp, na )
USE constants, ONLY : eps8
REAL(DP), intent(inout) :: stress(3,3),nstress(3,3)
REAL(DP), intent(in) :: pmass(:), omega, h(3,3), vels(:,:)
integer, intent(in) :: nsp, na(:)
integer :: i, j, is, ia, isa
nstress = 0.0_DP ! BS
isa = 0
if( omega < eps8 ) call errore(' ions_thermal_stress ', ' omega <= 0 ', 1 )
do is = 1, nsp
do ia = 1, na(is)
isa = isa + 1
do i = 1, 3
do j = 1, 3
stress(i,j) = stress(i,j) + pmass(is) / omega * &
& ( (h(i,1)*vels(1,isa)+h(i,2)*vels(2,isa)+h(i,3)*vels(3,isa)) * &
(h(j,1)*vels(1,isa)+h(j,2)*vels(2,isa)+h(j,3)*vels(3,isa)) )
! BS : to print Pressure of Nuclei at each step
nstress(i,j) = nstress(i,j) + pmass(is) / omega * &
& ( (h(i,1)*vels(1,isa)+h(i,2)*vels(2,isa)+h(i,3)*vels(3,isa)) * &
(h(j,1)*vels(1,isa)+h(j,2)*vels(2,isa)+h(j,3)*vels(3,isa)) )
! BS
enddo
enddo
enddo
enddo
return
end subroutine ions_thermal_stress
!------------------------------------------------------------------------------!
subroutine randvel( tempw, tau0, taum, na, nsp, iforce, &
amass, delt )
use constants, only: pi, k_boltzmann_au, amu_au
USE random_numbers, ONLY : randy
implicit none
REAL(DP), intent(out) :: taum(:,:)
REAL(DP), intent(in) :: tau0(:,:)
REAL(DP), intent(in) :: delt, amass(:), tempw
integer, intent(in) :: na(:), nsp
integer, intent(in) :: iforce(:,:)
REAL(DP) :: alfap, qr(3), alfar, gausp
REAL(DP) :: dt2by2
integer :: i, ia, is, nat, isa
dt2by2 = 0.5_DP * delt * delt
gausp = delt * sqrt( tempw * k_boltzmann_au )
nat = SUM( na( 1:nsp ) )
! randomize velocities, calculate center of mass vel.
do i=1,3
qr(i)=0.0_DP
isa = 0
do is=1,nsp
do ia=1,na(is)
isa = isa + 1
alfar=gausp/sqrt(amass(is)*amu_au)*cos(2.0_DP*pi*randy())*sqrt(-2.0_DP*log(randy()))
taum(i,isa)=alfar
qr(i)=qr(i)+alfar
end do
end do
qr(i)=qr(i)/nat
end do
!subtract center of mass velocity
isa = 0
do is=1,nsp
do ia=1,na(is)
isa = isa + 1
do i=1,3
alfar=taum(i,isa)-qr(i)
taum(i,isa)=tau0(i,isa)-iforce(i,isa)* &
& alfar
end do
end do
end do
return
end subroutine randvel
!------------------------------------------------------------------------------!
subroutine ions_vrescal( tcap, tempw, tempp, taup, tau0, taum, na, nsp, fion, iforce, &
pmass, delt )
use constants, only: pi, k_boltzmann_au, eps8
USE random_numbers, ONLY : randy
implicit none
logical, intent(in) :: tcap
REAL(DP), intent(inout) :: taup(:,:)
REAL(DP), intent(in) :: tau0(:,:), taum(:,:), fion(:,:)
REAL(DP), intent(in) :: delt, pmass(:), tempw, tempp
integer, intent(in) :: na(:), nsp
integer, intent(in) :: iforce(:,:)
REAL(DP) :: alfap, qr(3), alfar, gausp
REAL(DP) :: dt2by2
integer :: i, ia, is, nat, isa
dt2by2 = 0.5_DP * delt * delt
gausp = delt * sqrt( tempw * k_boltzmann_au )
nat = SUM( na( 1:nsp ) )
if(.not.tcap) then
if( tempp < eps8 ) call errore(' ions_vrescal ', ' tempp <= 0 ', 1 )
alfap = 0.5_DP * sqrt(tempw/tempp)
isa = 0
do is=1,nsp
do ia=1,na(is)
isa = isa + 1
do i=1,3
taup(i,isa) = tau0(i,isa) + &
& alfap*(taup(i,isa)-taum(i,isa)) + &
& dt2by2/pmass(is)*fion(i,isa)*iforce(i,isa)
end do
end do
end do
else
! randomize velocities, calculate center of mass vel.
do i=1,3
qr(i)=0.0_DP
isa = 0
do is=1,nsp
do ia=1,na(is)
isa = isa + 1
alfar=gausp/sqrt(pmass(is))*cos(2.0_DP*pi*randy())*sqrt(-2.0_DP*log(randy()))
taup(i,isa)=alfar
qr(i)=qr(i)+alfar
end do
end do
qr(i)=qr(i)/nat
end do
!subtract center of mass velocity
isa = 0
do is=1,nsp
do ia=1,na(is)
isa = isa + 1
do i=1,3
alfar=taup(i,isa)-qr(i)
taup(i,isa)=tau0(i,isa)+iforce(i,isa)* &
& (alfar+dt2by2/pmass(is)*fion(i,isa))
end do
end do
end do
end if
return
end subroutine ions_vrescal
!------------------------------------------------------------------------------!
subroutine ions_shiftvar( varp, var0, varm )
implicit none
REAL(DP), intent(in) :: varp(:,:)
REAL(DP), intent(out) :: varm(:,:), var0(:,:)
varm = var0
var0 = varp
return
end subroutine ions_shiftvar
!------------------------------------------------------------------------------!
SUBROUTINE ions_reference_positions( tau )
! Calculate the real position of atoms relative to the center of mass (cdm)
! and store them in taui
! cdmi: initial position of the center of mass (cdm) in cartesian coor.
IMPLICIT NONE
REAL(DP) :: tau( :, : )
INTEGER :: isa
CALL ions_cofmass( tau, amass, na, nsp, cdmi )
DO isa = 1, nat
taui(:,isa) = tau(:,isa) - cdmi(:)
END DO
RETURN
END SUBROUTINE ions_reference_positions
!------------------------------------------------------------------------------!
SUBROUTINE ions_displacement( dis, tau )
! Calculate the sum of the quadratic displacements of the atoms in the ref.
! of cdm respect to the initial positions.
! taui: initial positions in real units in the ref. of cdm
! ----------------------------------------------
! att! tau_ref: starting position in center-of-mass ref. in real units
! ----------------------------------------------
IMPLICIT NONE
REAL (DP), INTENT(OUT) :: dis(:)
REAL (DP), INTENT(IN) :: tau(:,:)
REAL(DP) :: rdist(3), r2, cdm(3)
INTEGER :: is, ia, isa
! ... Compute the current value of cdm "Center of Mass"
!
CALL ions_cofmass(tau, amass, na, nsp, cdm )
!
IF( SIZE( dis ) < nsp ) &
CALL errore(' displacement ',' size of dis too small ', 1)
isa = 0
DO is = 1, nsp
dis(is) = 0.0_DP
r2 = 0.0_DP
DO ia = 1, na(is)
isa = isa + 1
rdist = tau(:,isa) - cdm
r2 = r2 + SUM( ( rdist(:) - taui(:,isa) )**2 )
END DO
dis(is) = dis(is) + r2 / DBLE(na(is))
END DO
RETURN
END SUBROUTINE ions_displacement
!--------------------------------------------------------------------------
SUBROUTINE ions_cofmsub( tausp, iforce, nat, cdm, cdm0 )
!--------------------------------------------------------------------------
!
IMPLICIT NONE
!
REAL(DP), INTENT(INOUT) :: tausp(:,:)
INTEGER, INTENT(IN) :: iforce(:,:)
INTEGER, INTENT(IN) :: nat
REAL(DP), INTENT(IN) :: cdm(:), cdm0(:)
!
INTEGER :: i, ia
!
DO ia = 1, nat
!
DO i = 1, 3
!
tausp(i,ia) = tausp(i,ia) + DBLE( iforce(i,ia) ) * ( cdm0(i) - cdm(i) )
!
END DO
!
END DO
!
RETURN
!
END SUBROUTINE ions_cofmsub
REAL(DP) FUNCTION compute_eextfor( tau0 )
IMPLICIT NONE
REAL(DP), OPTIONAL, INTENT(IN) :: tau0(:,:)
INTEGER :: i
REAL(DP) :: e
compute_eextfor = 0.0d0
e = 0.0d0
IF( PRESENT( tau0 ) ) THEN
DO i = 1, SIZE( extfor,2 )
e = e + extfor( 3, i ) * tau0( 3, i ) &
+ extfor( 2, i ) * tau0( 2, i ) &
+ extfor( 1, i ) * tau0( 1, i )
END DO
ELSE
DO i = 1, SIZE( extfor,2 )
e = e + extfor( 3, i ) * tau( 3, i ) &
+ extfor( 2, i ) * tau( 2, i ) &
+ extfor( 1, i ) * tau( 1, i )
END DO
END IF
compute_eextfor = - e
RETURN
END FUNCTION compute_eextfor
!------------------------------------------------------------------------------!
END MODULE ions_base
!------------------------------------------------------------------------------!
Reference in New Issue View Git Blame Copy Permalink