quantum_package/plugins/Symmetry/sym_operation.irp.f

71 lines
2.4 KiB
Fortran

subroutine sym_apply_inversion(point_in, point_out)
implicit none
double precision, intent(in) :: point_in(3)
double precision, intent(out) :: point_out(3)
point_out(:) = -point_in(:)
end
subroutine sym_apply_reflexion(iaxis,point_in,point_out)
implicit none
integer, intent(in) :: iaxis
double precision, intent(in) :: point_in(3)
double precision, intent(out) :: point_out(3)
point_out(:) = point_in(:)
point_out(iaxis) = -point_in(iaxis)
end
subroutine sym_apply_diagonal_reflexion(angle,iaxis,point_in,point_out)
implicit none
integer, intent(in) :: iaxis
double precision, intent(in) :: point_in(3), angle
double precision, intent(out) :: point_out(3)
double precision :: point_tmp1(3), point_tmp2(3)
integer :: iaxis2
iaxis2 = mod(iaxis,3)+1
! iaxis2 = mod(iaxis2,3)+1
call sym_apply_rotation(-angle,iaxis,point_in,point_tmp1)
call sym_apply_reflexion(iaxis2,point_tmp1,point_tmp2)
call sym_apply_rotation(angle,iaxis,point_tmp2,point_out)
end
subroutine sym_apply_rotation(angle,iaxis,point_in,point_out)
implicit none
double precision, intent(in) :: angle
integer, intent(in) :: iaxis
double precision, intent(in) :: point_in(3)
double precision, intent(out) :: point_out(3)
double precision :: theta, sin_t, cos_t
double precision, parameter :: pi = dacos(-1.d0)
theta = 2.d0*pi/angle
sin_t = dsin(theta)
cos_t = dcos(theta)
if (iaxis==1) then
point_out(1) = point_in(1)
point_out(2) = point_in(2)*cos_t - point_in(3)*sin_t
point_out(3) = point_in(2)*sin_t + point_in(3)*cos_t
else if (iaxis==2) then
point_out(1) = point_in(1)*cos_t - point_in(3)*sin_t
point_out(2) = point_in(2)
point_out(3) = point_in(1)*sin_t + point_in(3)*cos_t
endif
if (iaxis==3) then
point_out(1) = point_in(1)*cos_t - point_in(2)*sin_t
point_out(2) = point_in(1)*sin_t + point_in(2)*cos_t
point_out(3) = point_in(3)
endif
end
subroutine sym_apply_improper_rotation(angle,iaxis,point_in,point_out)
implicit none
double precision, intent(in) :: angle
integer, intent(in) :: iaxis
double precision, intent(in) :: point_in(3)
double precision, intent(out) :: point_out(3)
double precision :: point_tmp(3)
call sym_apply_rotation(angle,iaxis,point_in,point_tmp)
call sym_apply_reflexion(iaxis,point_tmp,point_out)
end