LCPQ
/
quantum_package
mirror of https://github.com/LCPQ/quantum_package
10
0
Fork 0
Code Issues Releases Wiki Activity

Merge pull request #220 from scemama/master 

Avoid using Core v0.10.0 and stick to v0.9.1
This commit is contained in:
Thomas Applencourt 2017-12-20 18:10:08 -06:00 committed by GitHub
parent 279c225f3b 37c07671da
commit cfe0d25d03
No known key found for this signature in database
GPG Key ID: 4AEE18F83AFDEB23
7 changed files with 91 additions and 76 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
install/scripts/install_ocaml.sh
View File

@ -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}"

2
plugins/Symmetry/Symmetry.main.irp.f
View File

@ -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

28
plugins/Symmetry/aos.irp.f
View File

@ -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))

27
plugins/Symmetry/find_sym.irp.f
View File

@ -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)

78
plugins/Symmetry/nuclei.irp.f
View File

@ -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

25
plugins/Symmetry/sym_operation.irp.f
View File

@ -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)

5
src/Nuclei/inertia.irp.f
View File

@ -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