diff --git a/install/scripts/install_ocaml.sh b/install/scripts/install_ocaml.sh index 88d38845..f322bd0b 100755 --- a/install/scripts/install_ocaml.sh +++ b/install/scripts/install_ocaml.sh @@ -5,7 +5,7 @@ QP_ROOT=$PWD cd - # Normal installation -PACKAGES="core cryptokit.1.10 ocamlfind sexplib ZMQ ppx_sexp_conv ppx_deriving" +PACKAGES="core.v0.9.1 cryptokit.1.10 ocamlfind sexplib.v0.9.1 ZMQ ppx_sexp_conv ppx_deriving" # Needed for ZeroMQ export C_INCLUDE_PATH="${QP_ROOT}"/include:"${C_INCLUDE_PATH}" diff --git a/plugins/Symmetry/Symmetry.main.irp.f b/plugins/Symmetry/Symmetry.main.irp.f index 1d438af5..9540295f 100644 --- a/plugins/Symmetry/Symmetry.main.irp.f +++ b/plugins/Symmetry/Symmetry.main.irp.f @@ -15,7 +15,7 @@ program Symmetry print *, 'Symmetry operations : ', sym_operation(1:n_irrep) print *, 'Character table' do i=1,n_irrep - print *, character_table(i,:) + print *, i, real(character_table(i,:)) enddo PROVIDE mo_sym end diff --git a/plugins/Symmetry/aos.irp.f b/plugins/Symmetry/aos.irp.f index 5515e1b7..1ed567bc 100644 --- a/plugins/Symmetry/aos.irp.f +++ b/plugins/Symmetry/aos.irp.f @@ -45,38 +45,27 @@ subroutine compute_sym_ao_values(sym_points, n_sym_points, result) double precision, intent(in) :: sym_points(3,n_sym_points) double precision, intent(out) :: result(n_sym_points, ao_num) integer :: i, j + double precision :: point(3) double precision :: x, y, z double precision :: x2, y2, z2 integer :: k - + result (:,:) = 0.d0 -print *, sym_molecule_rotation_inv -print *, '' -print *, sym_molecule_rotation -stop do j=1,ao_num do i=1,n_sym_points - x2 = sym_points(1,i) - y2 = sym_points(2,i) - z2 = sym_points(3,i) - x = x2*sym_molecule_rotation_inv(1,1) + y2*sym_molecule_rotation_inv(2,1) + z2*sym_molecule_rotation_inv(3,1) - y = x2*sym_molecule_rotation_inv(1,2) + y2*sym_molecule_rotation_inv(2,2) + z2*sym_molecule_rotation_inv(3,2) - z = x2*sym_molecule_rotation_inv(1,3) + y2*sym_molecule_rotation_inv(2,3) + z2*sym_molecule_rotation_inv(3,3) - x = x - nucl_coord_transp(1,ao_nucl(j)) - y = y - nucl_coord_transp(2,ao_nucl(j)) - z = z - nucl_coord_transp(3,ao_nucl(j)) + call point_to_input_orientation(sym_points(:,i), point) + x = point(1) - nucl_coord_transp(1,ao_nucl(j)) + y = point(2) - nucl_coord_transp(2,ao_nucl(j)) + z = point(3) - nucl_coord_transp(3,ao_nucl(j)) x2 = x*x + y*y + z*z result(i,j) = 0.d0 do k=1,ao_prim_num(j) result(i,j) += ao_coef_normalized_ordered_transp(k,j)*exp(-ao_expo_ordered_transp(k,j)*x2) enddo -print *, real(x), ao_power(j,1), real(y), ao_power(j,2), real(z), ao_power(j,3) x = x**ao_power(j,1) y = y**ao_power(j,2) z = z**ao_power(j,3) -print *, result(i,j) result(i,j) = x*y*z*result(i,j) - print *, result(i,j) enddo enddo @@ -94,11 +83,6 @@ subroutine compute_sym_mo_values(sym_points, n_sym_points, result) double precision, allocatable :: tmp(:,:) allocate(tmp(n_sym_points,ao_num)) call compute_sym_ao_values(sym_points,n_sym_points,tmp) -integer :: i -do i=1,ao_num - print *, tmp(:,i) -enddo -stop call dgemm('N','N',n_sym_points,mo_tot_num,ao_num, & 1.d0, tmp,size(tmp,1), mo_coef, size(mo_coef,1), & 0.d0, result,size(result,1)) diff --git a/plugins/Symmetry/find_sym.irp.f b/plugins/Symmetry/find_sym.irp.f index c7e2ba37..38c0a6b7 100644 --- a/plugins/Symmetry/find_sym.irp.f +++ b/plugins/Symmetry/find_sym.irp.f @@ -84,7 +84,7 @@ END_PROVIDER END_DOC molecule_principal_axis = maxloc(sym_rotation_axis,1) if (molecule_principal_axis == 1) then - if (sym_rotation_axis(2) >= sym_rotation_axis(3)) then + if (sym_rotation_axis(2) > sym_rotation_axis(3)) then molecule_secondary_axis = 2 molecule_ternary_axis = 3 else @@ -92,7 +92,7 @@ END_PROVIDER molecule_ternary_axis = 2 endif else if (molecule_principal_axis == 2) then - if (sym_rotation_axis(1) >= sym_rotation_axis(3)) then + if (sym_rotation_axis(1) > sym_rotation_axis(3)) then molecule_secondary_axis = 1 molecule_ternary_axis = 3 else @@ -100,7 +100,7 @@ END_PROVIDER molecule_ternary_axis = 1 endif else if (molecule_principal_axis == 3) then - if (sym_rotation_axis(1) >= sym_rotation_axis(2)) then + if (sym_rotation_axis(1) > sym_rotation_axis(2)) then molecule_secondary_axis = 1 molecule_ternary_axis = 2 else @@ -375,11 +375,11 @@ BEGIN_PROVIDER [ integer, mo_sym, (mo_tot_num) ] call generate_sym_coord(n_sym_points,ref_points) call compute_sym_mo_values(ref_points,n_sym_points,val(1,1,2)) + possible_irrep = .True. do iop=1,n_irrep if (sym_operation(iop) == 'E') then cycle endif - possible_irrep = .True. if (sym_operation(iop) == 'i') then do ipoint=1,n_sym_points @@ -395,11 +395,12 @@ BEGIN_PROVIDER [ integer, mo_sym, (mo_tot_num) ] enddo else if (sym_operation(iop) == 'sv') then do ipoint=1,n_sym_points - call sym_apply_reflexion(molecule_secondary_axis,ref_points(1,ipoint),sym_points(1,ipoint)) + call sym_apply_reflexion(molecule_ternary_axis,ref_points(1,ipoint),sym_points(1,ipoint)) enddo else if (sym_operation(iop) == 'sd') then + angle = dble(maxval(sym_rotation_axis)) do ipoint=1,n_sym_points - call sym_apply_diagonal_reflexion(molecule_principal_axis,ref_points(1,ipoint),sym_points(1,ipoint)) + call sym_apply_diagonal_reflexion(angle,molecule_principal_axis,ref_points(1,ipoint),sym_points(1,ipoint)) enddo else if (sym_operation(iop) == 'C2''') then angle = 2.d0 @@ -452,17 +453,23 @@ BEGIN_PROVIDER [ integer, mo_sym, (mo_tot_num) ] sym_operations_on_mos(imo) += x enddo sym_operations_on_mos(imo) *= 1.d0/icount - print *, iop, imo, sym_operations_on_mos(imo) - print *, val(:,imo,1) + if (dabs(sym_operations_on_mos(imo) - 1.d0) < 1.d-2) then + sym_operations_on_mos(imo) = 1.d0 + else if (dabs(sym_operations_on_mos(imo) + 1.d0) < 1.d-2) then + sym_operations_on_mos(imo) = -1.d0 + else if (dabs(sym_operations_on_mos(imo)) < 1.d-2) then + sym_operations_on_mos(imo) = 0.d0 + endif + print *, imo, sym_operations_on_mos(imo) do i=1,n_irrep - if (character_table(i,iop) /= sym_operations_on_mos(imo)) then + if (dabs(character_table(i,iop) - sym_operations_on_mos(imo)) > 1.d-2) then possible_irrep(i,imo) = .False. endif enddo enddo enddo do imo=1,mo_tot_num - print *, imo + print *, 'MO ', imo do i=1,n_irrep if (possible_irrep(i,imo)) then print *, sym_irrep(i) diff --git a/plugins/Symmetry/nuclei.irp.f b/plugins/Symmetry/nuclei.irp.f index caa4781e..d680393d 100644 --- a/plugins/Symmetry/nuclei.irp.f +++ b/plugins/Symmetry/nuclei.irp.f @@ -1,3 +1,56 @@ +subroutine point_to_standard_orientation(point_in,point_out) + implicit none + double precision, intent(in) :: point_in(3) + double precision, intent(out) :: point_out(3) + BEGIN_DOC + ! Returns the coordinates of a point in the standard orientation + END_DOC + double precision :: point_tmp(3) + + point_tmp(1) = point_in(1) - center_of_mass(1) + point_tmp(2) = point_in(2) - center_of_mass(2) + point_tmp(3) = point_in(3) - center_of_mass(3) + + point_out(1) = point_tmp(1)*inertia_tensor_eigenvectors(1,1) + & + point_tmp(2)*inertia_tensor_eigenvectors(2,1) + & + point_tmp(3)*inertia_tensor_eigenvectors(3,1) + + point_out(2) = point_tmp(1)*inertia_tensor_eigenvectors(1,2) + & + point_tmp(2)*inertia_tensor_eigenvectors(2,2) + & + point_tmp(3)*inertia_tensor_eigenvectors(3,2) + + point_out(3) = point_tmp(1)*inertia_tensor_eigenvectors(1,3) + & + point_tmp(2)*inertia_tensor_eigenvectors(2,3) + & + point_tmp(3)*inertia_tensor_eigenvectors(3,3) + +end + +subroutine point_to_input_orientation(point_in,point_out) + implicit none + double precision, intent(in) :: point_in(3) + double precision, intent(out) :: point_out(3) + BEGIN_DOC + ! Returns the coordinates of a point in the input orientation + END_DOC + double precision :: point_tmp(3) + + point_tmp(1) = point_in(1)*inertia_tensor_eigenvectors(1,1) + & + point_in(2)*inertia_tensor_eigenvectors(1,2) + & + point_in(3)*inertia_tensor_eigenvectors(1,3) + + point_tmp(2) = point_in(1)*inertia_tensor_eigenvectors(2,1) + & + point_in(2)*inertia_tensor_eigenvectors(2,2) + & + point_in(3)*inertia_tensor_eigenvectors(2,3) + + point_tmp(3) = point_in(1)*inertia_tensor_eigenvectors(3,1) + & + point_in(2)*inertia_tensor_eigenvectors(3,2) + & + point_in(3)*inertia_tensor_eigenvectors(3,3) + + point_out(1) = point_tmp(1) + center_of_mass(1) + point_out(2) = point_tmp(2) + center_of_mass(2) + point_out(3) = point_tmp(3) + center_of_mass(3) + +end BEGIN_PROVIDER [ double precision, nucl_coord_sym, (nucl_num,3) ] implicit none @@ -9,15 +62,7 @@ BEGIN_PROVIDER [ double precision, nucl_coord_sym, (nucl_num,3) ] if (mpi_master) then integer :: i do i=1,nucl_num - nucl_coord_sym(i,1) = (nucl_coord(i,1) - center_of_mass(1))*inertia_tensor_eigenvectors(1,1) + & - (nucl_coord(i,2) - center_of_mass(2))*inertia_tensor_eigenvectors(2,1) + & - (nucl_coord(i,3) - center_of_mass(3))*inertia_tensor_eigenvectors(3,1) - nucl_coord_sym(i,2) = (nucl_coord(i,1) - center_of_mass(1))*inertia_tensor_eigenvectors(1,2) + & - (nucl_coord(i,2) - center_of_mass(2))*inertia_tensor_eigenvectors(2,2) + & - (nucl_coord(i,3) - center_of_mass(3))*inertia_tensor_eigenvectors(3,2) - nucl_coord_sym(i,3) = (nucl_coord(i,1) - center_of_mass(1))*inertia_tensor_eigenvectors(1,3) + & - (nucl_coord(i,2) - center_of_mass(2))*inertia_tensor_eigenvectors(2,3) + & - (nucl_coord(i,3) - center_of_mass(3))*inertia_tensor_eigenvectors(3,3) + call point_to_standard_orientation(nucl_coord(i,:), nucl_coord_sym(i,:)) enddo character*(64), parameter :: f = '(A16, 4(1X,F12.6))' @@ -74,19 +119,4 @@ BEGIN_PROVIDER [ double precision, nucl_coord_sym_transp, (3,nucl_num) ] enddo END_PROVIDER -BEGIN_PROVIDER [ double precision, sym_molecule_rotation, (3,3) ] - implicit none - BEGIN_DOC - ! Rotation of the molecule to go from input orientation to standard orientation - END_DOC - call find_rotation(nucl_coord, size(nucl_coord,1), nucl_coord_sym, 3, sym_molecule_rotation, 3) -END_PROVIDER - -BEGIN_PROVIDER [ double precision, sym_molecule_rotation_inv, (3,3) ] - implicit none - BEGIN_DOC - ! Rotation of the molecule to go from standard orientation to input orientation - END_DOC - call find_rotation(nucl_coord_sym, size(nucl_coord_sym,1), nucl_coord, 3, sym_molecule_rotation_inv, 3) -END_PROVIDER diff --git a/plugins/Symmetry/sym_operation.irp.f b/plugins/Symmetry/sym_operation.irp.f index 93e580e2..cfc86621 100644 --- a/plugins/Symmetry/sym_operation.irp.f +++ b/plugins/Symmetry/sym_operation.irp.f @@ -15,25 +15,18 @@ subroutine sym_apply_reflexion(iaxis,point_in,point_out) end -subroutine sym_apply_diagonal_reflexion(iaxis,point_in,point_out) +subroutine sym_apply_diagonal_reflexion(angle,iaxis,point_in,point_out) implicit none integer, intent(in) :: iaxis - double precision, intent(in) :: point_in(3) + double precision, intent(in) :: point_in(3), angle double precision, intent(out) :: point_out(3) - if (iaxis == 1) then - point_out(1) = point_in(1) - point_out(2) = point_in(3) - point_out(3) = point_in(2) - else if (iaxis == 2) then - point_out(1) = point_in(3) - point_out(2) = point_in(2) - point_out(3) = point_in(1) - else if (iaxis == 3) then - point_out(1) = point_in(2) - point_out(2) = point_in(1) - point_out(3) = point_in(3) - endif - + 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) diff --git a/src/Nuclei/inertia.irp.f b/src/Nuclei/inertia.irp.f index d2cb7cf8..c79fa243 100644 --- a/src/Nuclei/inertia.irp.f +++ b/src/Nuclei/inertia.irp.f @@ -25,8 +25,9 @@ END_PROVIDER ! Eigenvectors/eigenvalues of the inertia_tensor. Used to find normal orientation. END_DOC integer :: k - call lapack_diagd(inertia_tensor_eigenvalues,inertia_tensor_eigenvectors,inertia_tensor,3,3) + call lapack_diagd(inertia_tensor_eigenvalues,inertia_tensor_eigenvectors,-inertia_tensor,3,3) + inertia_tensor_eigenvalues = -inertia_tensor_eigenvalues print *, 'Rotational constants (GHZ):' - print *, (1805.65468542d0/(inertia_tensor_eigenvalues(k)+1.d-32), k=3,1,-1) + print *, (1805.65468542d0/(inertia_tensor_eigenvalues(k)+1.d-32), k=1,3) END_PROVIDER