BEGIN_PROVIDER [ double precision, sym_transformation_matrices, (3,3,n_irrep) ] implicit none BEGIN_DOC ! Transformation matrices for the symmetry operations END_DOC integer :: iop, iangle, iaxis, iaxis2, l double precision :: angle, point_tmp1(3), point_tmp2(3), A(3,3) sym_transformation_matrices = 0.d0 do iop=1,n_irrep if (sym_operation(iop) == 'E') then sym_transformation_matrices(1,1,iop) = 1.d0 sym_transformation_matrices(2,2,iop) = 1.d0 sym_transformation_matrices(3,3,iop) = 1.d0 else if (sym_operation(iop) == 'i') then sym_transformation_matrices(1,1,iop) = -1.d0 sym_transformation_matrices(2,2,iop) = -1.d0 sym_transformation_matrices(3,3,iop) = -1.d0 else if (sym_operation(iop) == 'sh') then sym_transformation_matrices(1,1,iop) = 1.d0 sym_transformation_matrices(2,2,iop) = 1.d0 sym_transformation_matrices(3,3,iop) = 1.d0 iaxis = molecule_principal_axis sym_transformation_matrices(iaxis,iaxis,iop) = -1.d0 else if (sym_operation(iop) == 's') then sym_transformation_matrices(1,1,iop) = 1.d0 sym_transformation_matrices(2,2,iop) = 1.d0 sym_transformation_matrices(3,3,iop) = 1.d0 iaxis = molecule_principal_axis sym_transformation_matrices(iaxis,iaxis,iop) = -1.d0 else if (sym_operation(iop) == 'sv') then sym_transformation_matrices(1,1,iop) = 1.d0 sym_transformation_matrices(2,2,iop) = 1.d0 sym_transformation_matrices(3,3,iop) = 1.d0 iaxis = molecule_ternary_axis sym_transformation_matrices(iaxis,iaxis,iop) = -1.d0 else if (sym_operation(iop) == 'sd') then angle = dble(maxval(sym_rotation_axis)) sym_transformation_matrices(1,1,iop) = 1.d0 sym_transformation_matrices(2,2,iop) = 1.d0 sym_transformation_matrices(3,3,iop) = 1.d0 iaxis = molecule_principal_axis iaxis2 = molecule_secondary_axis call sym_apply_rotation(-angle,iaxis,sym_transformation_matrices(1,1,iop),point_tmp1) call sym_apply_reflexion(iaxis2,point_tmp1,point_tmp2) call sym_apply_rotation(angle,iaxis,point_tmp2,sym_transformation_matrices(1,1,iop)) call sym_apply_rotation(-angle,iaxis,sym_transformation_matrices(1,2,iop),point_tmp1) call sym_apply_reflexion(iaxis2,point_tmp1,point_tmp2) call sym_apply_rotation(angle,iaxis,point_tmp2,sym_transformation_matrices(1,2,iop)) call sym_apply_rotation(-angle,iaxis,sym_transformation_matrices(1,3,iop),point_tmp1) call sym_apply_reflexion(iaxis2,point_tmp1,point_tmp2) call sym_apply_rotation(angle,iaxis,point_tmp2,sym_transformation_matrices(1,3,iop)) else if (sym_operation(iop) == 'C2''') then angle = 2.d0 sym_transformation_matrices(1,1,iop) = 1.d0 sym_transformation_matrices(2,2,iop) = 1.d0 sym_transformation_matrices(3,3,iop) = 1.d0 iaxis = molecule_secondary_axis call sym_apply_rotation(angle,iaxis,sym_transformation_matrices(1,1,iop),point_tmp1) sym_transformation_matrices(1:3,1,iop) = point_tmp1(1:3) call sym_apply_rotation(angle,iaxis,sym_transformation_matrices(1,2,iop),point_tmp1) sym_transformation_matrices(1:3,2,iop) = point_tmp1(1:3) call sym_apply_rotation(angle,iaxis,sym_transformation_matrices(1,3,iop),point_tmp1) sym_transformation_matrices(1:3,3,iop) = point_tmp1(1:3) else if (sym_operation(iop) == 'C2"') then angle = 2.d0 sym_transformation_matrices(1,1,iop) = 1.d0 sym_transformation_matrices(2,2,iop) = 1.d0 sym_transformation_matrices(3,3,iop) = 1.d0 iaxis = molecule_ternary_axis call sym_apply_rotation(angle,iaxis,sym_transformation_matrices(1,1,iop),point_tmp1) sym_transformation_matrices(1:3,1,iop) = point_tmp1(1:3) call sym_apply_rotation(angle,iaxis,sym_transformation_matrices(1,2,iop),point_tmp1) sym_transformation_matrices(1:3,2,iop) = point_tmp1(1:3) call sym_apply_rotation(angle,iaxis,sym_transformation_matrices(1,3,iop),point_tmp1) sym_transformation_matrices(1:3,3,iop) = point_tmp1(1:3) else do l=2,len(sym_operation(iop)) if (sym_operation(iop)(l:l) == '^') exit enddo read(sym_operation(iop)(2:l-1), *) iangle if (l == len(sym_operation(iop))+1) then l=1 else read(sym_operation(iop)(l+1:), *, err=10, end=10) l 10 continue endif angle = dble(iangle)/(dble(l)) if (sym_operation(iop)(1:1) == 'C') then sym_transformation_matrices(1,1,iop) = 1.d0 sym_transformation_matrices(2,2,iop) = 1.d0 sym_transformation_matrices(3,3,iop) = 1.d0 iaxis = molecule_principal_axis call sym_apply_rotation(angle,iaxis,sym_transformation_matrices(1,1,iop),point_tmp1) sym_transformation_matrices(1:3,1,iop) = point_tmp1(1:3) call sym_apply_rotation(angle,iaxis,sym_transformation_matrices(1,2,iop),point_tmp1) sym_transformation_matrices(1:3,2,iop) = point_tmp1(1:3) call sym_apply_rotation(angle,iaxis,sym_transformation_matrices(1,3,iop),point_tmp1) sym_transformation_matrices(1:3,3,iop) = point_tmp1(1:3) else if (sym_operation(iop)(1:1) == 'S') then sym_transformation_matrices(1,1,iop) = 1.d0 sym_transformation_matrices(2,2,iop) = 1.d0 sym_transformation_matrices(3,3,iop) = 1.d0 iaxis = molecule_principal_axis call sym_apply_rotation(angle,iaxis,sym_transformation_matrices(1,1,iop),point_tmp1) call sym_apply_reflexion(iaxis,point_tmp1,point_tmp2) sym_transformation_matrices(1:3,1,iop) = point_tmp2(1:3) call sym_apply_rotation(angle,iaxis,sym_transformation_matrices(1,2,iop),point_tmp1) call sym_apply_reflexion(iaxis,point_tmp1,point_tmp2) sym_transformation_matrices(1:3,2,iop) = point_tmp2(1:3) call sym_apply_rotation(angle,iaxis,sym_transformation_matrices(1,3,iop),point_tmp1) call sym_apply_reflexion(iaxis,point_tmp1,point_tmp2) sym_transformation_matrices(1:3,3,iop) = point_tmp2(1:3) endif endif enddo END_PROVIDER BEGIN_PROVIDER [ double precision, ao_symm, (ao_num,ao_num,n_irrep) ] implicit none BEGIN_DOC ! where O is a symmetry operation END_DOC double precision, allocatable :: sym_points(:,:), ref_points(:,:) double precision, allocatable :: val(:,:,:) integer :: iangle, n_sym_points integer :: iop, i,j,k n_sym_points = 100000*n_irrep double precision :: f allocate(val(n_sym_points,ao_num,2), sym_points(3,n_sym_points), ref_points(3,n_sym_points)) call generate_sym_coord(n_sym_points,ref_points) ! call compute_sym_ao_values(ref_points, n_sym_points, val(1,1,1)) ! f = 1.d0/(dble(n_sym_points)) do iop=1,n_irrep call dgemm('N','N',3,n_sym_points,3,1.d0,sym_transformation_matrices(1,1,iop), & size(sym_transformation_matrices,1),& ref_points,size(ref_points,1),0.d0,sym_points,size(sym_points,1)) call compute_sym_ao_values(sym_points, n_sym_points, val(1,1,2)) do i=1,ao_num do j=1,i ao_symm(i,j,iop) = 0.d0 do k=1,n_sym_points ao_symm(i,j,iop) = ao_symm(i,j,iop) + val(k,i,1)*val(k,j,2) enddo ao_symm(i,j,iop) = ao_symm(i,j,iop) * f enddo enddo enddo do iop=2,n_irrep do i=1,ao_num do j=1,i f = 1.d0/dsqrt(ao_symm(i,i,1)*ao_symm(j,j,1)) ao_symm(i,j,iop) = ao_symm(i,j,iop)*f ao_symm(j,i,iop) = ao_symm(i,j,iop) enddo enddo enddo do i=1,ao_num do j=1,i-1 f = 1.d0/dsqrt(ao_symm(i,i,1)*ao_symm(j,j,1)) ao_symm(i,j,1) = ao_symm(i,j,1)*f ao_symm(j,i,1) = ao_symm(i,j,1) enddo enddo do i=1,ao_num ao_symm(i,i,1) = 1.d0 enddo END_PROVIDER BEGIN_PROVIDER [ double precision, mo_symm, (mo_tot_num,mo_tot_num,n_irrep) ] implicit none BEGIN_DOC ! where O is a symmetry operation END_DOC integer :: iop do iop=1,n_irrep call ao_to_mo(ao_symm(1,1,iop),size(ao_symm,1),mo_symm(1,1,iop),size(mo_symm,1)) enddo END_PROVIDER