LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-06-25 22:52:15 +02:00
Code Issues Releases Wiki Activity

CISD is cleaned. Perturbation/selection is broken

Browse Source
This commit is contained in:
Anthony Scemama 2014-05-25 01:18:41 +02:00
parent d5bbe35bff
commit 583e8859d0
24 changed files with 476 additions and 339 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

55
scripts/generate_h_apply.py
View File

@ -8,10 +8,14 @@ file.close()
keywords = """
subroutine
parameters
params_main
initialization
declarations
decls_main
keys_work
copy_buffer
finalization
generate_psi_guess
""".split()
class H_apply(object):
@ -21,9 +25,8 @@ class H_apply(object):
for k in keywords:
s[k] = ""
s["subroutine"] = "H_apply_%s"%(sub)
s["params_post"] = ""
self.openmp = openmp
if openmp:
s["subroutine"] += "_OpenMP"
self.selection_pt2 = None
self.perturbation = None
@ -47,6 +50,33 @@ class H_apply(object):
s["omp_do"] = "!$OMP DO SCHEDULE (static)"
s["omp_enddo"] = "!$OMP ENDDO NOWAIT"
s["keys_work"] += "call fill_H_apply_buffer_no_selection(key_idx,keys_out,N_int,iproc)"
s["generate_psi_guess"] = """
! Sort H_jj to find the N_states lowest states
integer :: i,k
integer, allocatable :: iorder(:)
double precision, allocatable :: H_jj(:)
double precision, external :: diag_h_mat_elem
allocate(H_jj(N_det),iorder(N_det))
!$OMP PARALLEL DEFAULT(NONE) &
!$OMP SHARED(psi_det,N_int,H_jj,iorder,N_det) &
!$OMP PRIVATE(i)
!$OMP DO
do i = 1, N_det
H_jj(i) = diag_h_mat_elem(psi_det(1,1,i),N_int)
iorder(i) = i
enddo
!$OMP END DO
!$OMP END PARALLEL
call dsort(H_jj,iorder,N_det)
do k=1,N_states
psi_coef(iorder(k),k) = 1.d0
enddo
deallocate(H_jj,iorder)
"""
if not openmp:
for k in s:
s[k] = ""
@ -54,6 +84,7 @@ class H_apply(object):
s["set_i_H_j_threshold"] = """
thresh = H_apply_threshold
"""
s["copy_buffer"] = "call copy_h_apply_buffer_to_wf"
self.data = s
def __setitem__(self,key,value):
@ -101,6 +132,21 @@ class H_apply(object):
sum_norm_pert_in = sum_norm_pert
sum_H_pert_diag_in = sum_H_pert_diag
"""
self.data["copy_buffer"] = ""
self.data["generate_psi_guess"] = ""
self.data["params_main"] = "pt2, norm_pert, H_pert_diag, N_st"
self.data["params_post"] = ","+self.data["params_main"] +", N_int"
self.data["decls_main"] = """ integer, intent(in) :: N_st
double precision, intent(inout):: pt2(N_st)
double precision, intent(inout):: norm_pert(N_st)
double precision, intent(inout):: H_pert_diag
PROVIDE reference_energy N_det_generators key_pattern_not_in_ref
pt2 = 0.d0
norm_pert = 0.d0
H_pert_diag = 0.d0
"""
if self.openmp:
self.data["omp_parallel"] += """&
!$OMP SHARED(N_st,Nint) PRIVATE(e_2_pert_buffer,coef_pert_buffer) &
@ -113,6 +159,11 @@ class H_apply(object):
self.selection_pt2 = pert
if pert is not None:
self.data["size_max"] = str(1024*128)
self.data["copy_buffer"] = """
call copy_h_apply_buffer_to_wf
selection_criterion_min = selection_criterion_min*0.1d0
selection_criterion = selection_criterion_min
"""
self.data["keys_work"] = """
e_2_pert_buffer = 0.d0
coef_pert_buffer = 0.d0

6
src/Bitmask/README.rst
View File

@ -69,9 +69,15 @@ Documentation
`n_int <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks.irp.f#L3>`_
Number of 64-bit integers needed to represent determinants as binary strings
`n_reference_bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks.irp.f#L112>`_
Number of bitmasks for reference
`ref_bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks.irp.f#L50>`_
Reference bit mask, used in Slater rules, chosen as Hartree-Fock bitmask
`reference_bitmask <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks.irp.f#L145>`_
Bitmasks for reference determinants. (N_int, alpha/beta, hole/particle, reference)
`bitstring_to_hexa <http://github.com/LCPQ/quantum_package/tree/master/src/Bitmask/bitmasks_routines.irp.f#L95>`_
Transform a bit string to a string in hexadecimal format for printing

4
src/Bitmask/bitmasks.ezfio_config
View File

@ -2,7 +2,7 @@ bitmasks
N_int integer
bit_kind integer
N_mask_gen integer
generators integer (bitmasks_N_int*bitmasks_bit_kind/4,2,2,bitmasks_N_mask_gen)
generators integer*8 (bitmasks_N_int*bitmasks_bit_kind/8,2,6,bitmasks_N_mask_gen)
N_mask_ref integer
reference integer (bitmasks_N_int*bitmasks_bit_kind/4,2,2,bitmasks_N_mask_ref)
reference integer*8 (bitmasks_N_int*bitmasks_bit_kind/8,2,6,bitmasks_N_mask_ref)

93
src/Bitmask/bitmasks.irp.f
View File

@ -88,10 +88,17 @@ BEGIN_PROVIDER [ integer, N_generators_bitmask ]
END_PROVIDER
BEGIN_PROVIDER [ integer(bit_kind), generators_bitmask, (N_int,2,2,N_generators_bitmask) ]
BEGIN_PROVIDER [ integer(bit_kind), generators_bitmask, (N_int,2,6,N_generators_bitmask) ]
implicit none
BEGIN_DOC
! Bitmasks for generator determinants. (N_int, alpha/beta, hole/particle, generator)
! Bitmasks for generator determinants. (N_int, alpha/beta, hole/particle, generator).
! 3rd index is :
! * 1 : hole for single exc
! * 1 : particle for single exc
! * 3 : hole for 1st exc of double
! * 4 : particle for 1st exc of double
! * 5 : hole for 2dn exc of double
! * 6 : particle for 2dn exc of double
END_DOC
logical :: exists
PROVIDE ezfio_filename
@ -100,12 +107,86 @@ BEGIN_PROVIDER [ integer(bit_kind), generators_bitmask, (N_int,2,2,N_generators_
if (exists) then
call ezfio_get_bitmasks_generators(generators_bitmask)
else
generators_bitmask(:,:,1,1) = HF_bitmask
generators_bitmask(:,1,2,1) = iand(not(HF_bitmask(:,1)),full_ijkl_bitmask(:,1))
generators_bitmask(:,2,2,1) = iand(not(HF_bitmask(:,2)),full_ijkl_bitmask(:,2))
generators_bitmask(:,:,s_hole ,1) = HF_bitmask
generators_bitmask(:,:,s_part ,1) = iand(not(HF_bitmask(:,:)),full_ijkl_bitmask(:,:))
generators_bitmask(:,:,d_hole1,1) = HF_bitmask
generators_bitmask(:,:,d_part1,1) = iand(not(HF_bitmask(:,:)),full_ijkl_bitmask(:,:))
generators_bitmask(:,:,d_hole2,1) = HF_bitmask
generators_bitmask(:,:,d_part2,1) = iand(not(HF_bitmask(:,:)),full_ijkl_bitmask(:,:))
call ezfio_set_bitmasks_generators(generators_bitmask)
endif
ASSERT (N_generators_bitmask > 0)
END_PROVIDER
BEGIN_PROVIDER [ integer, N_reference_bitmask ]
implicit none
BEGIN_DOC
! Number of bitmasks for reference
END_DOC
logical :: exists
PROVIDE ezfio_filename
call ezfio_has_bitmasks_N_mask_ref(exists)
if (exists) then
call ezfio_get_bitmasks_N_mask_ref(N_reference_bitmask)
integer :: N_int_check
integer :: bit_kind_check
call ezfio_get_bitmasks_bit_kind(bit_kind_check)
if (bit_kind_check /= bit_kind) then
print *, bit_kind_check, bit_kind
print *, 'Error: bit_kind is not correct in EZFIO file'
endif
call ezfio_get_bitmasks_N_int(N_int_check)
if (N_int_check /= N_int) then
print *, N_int_check, N_int
print *, 'Error: N_int is not correct in EZFIO file'
endif
else
N_reference_bitmask = 1
call ezfio_set_bitmasks_N_int(N_int)
call ezfio_set_bitmasks_bit_kind(bit_kind)
call ezfio_set_bitmasks_N_mask_ref(N_reference_bitmask)
endif
ASSERT (N_reference_bitmask > 0)
END_PROVIDER
BEGIN_PROVIDER [ integer(bit_kind), reference_bitmask, (N_int,2,2,N_reference_bitmask) ]
implicit none
BEGIN_DOC
! Bitmasks for reference determinants. (N_int, alpha/beta, hole/particle, reference)
END_DOC
logical :: exists
PROVIDE ezfio_filename
call ezfio_has_bitmasks_reference(exists)
if (exists) then
call ezfio_get_bitmasks_reference(reference_bitmask)
else
reference_bitmask(:,:,s_hole ,1) = HF_bitmask
reference_bitmask(:,:,s_part ,1) = iand(not(HF_bitmask(:,:)),full_ijkl_bitmask(:,:))
reference_bitmask(:,:,d_hole1,1) = HF_bitmask
reference_bitmask(:,:,d_part1,1) = iand(not(HF_bitmask(:,:)),full_ijkl_bitmask(:,:))
reference_bitmask(:,:,d_hole2,1) = HF_bitmask
reference_bitmask(:,:,d_part2,1) = iand(not(HF_bitmask(:,:)),full_ijkl_bitmask(:,:))
call ezfio_set_bitmasks_reference(reference_bitmask)
endif
END_PROVIDER
BEGIN_PROVIDER [ integer, i_bitmask_gen ]
implicit none
BEGIN_DOC
! Current bitmask for the generators
END_DOC
i_bitmask_gen = 1
END_PROVIDER
BEGIN_PROVIDER [ integer, i_bitmask_ref ]
implicit none
BEGIN_DOC
! Current bitmask for the reference
END_DOC
i_bitmask_ref = 1
END_PROVIDER

6
src/Bitmask/bitmasks_module.f90
View File

@ -2,4 +2,10 @@ module bitmasks
integer, parameter :: bit_kind_shift = 6 ! 5: 32 bits, 6: 64 bits
integer, parameter :: bit_kind_size = 64
integer, parameter :: bit_kind = 64/8
integer, parameter :: s_hole = 1
integer, parameter :: s_part = 2
integer, parameter :: d_hole1 = 3
integer, parameter :: d_part1 = 4
integer, parameter :: d_hole2 = 5
integer, parameter :: d_part2 = 6
end module

70
src/CISD/H_apply.irp.f
View File

@ -1,65 +1,9 @@
BEGIN_SHELL [ /bin/bash ]
./h_apply.py
! Generates subroutine H_apply_cisd
! ----------------------------------
BEGIN_SHELL [ /usr/bin/env python ]
from generate_h_apply import H_apply
H = H_apply("cisd",openmp=True)
print H
END_SHELL
subroutine fill_H_apply_buffer(n_selected,det_buffer,Nint,iproc)
use bitmasks
implicit none
BEGIN_DOC
! Fill the H_apply buffer with determiants for CISD
END_DOC
integer, intent(in) :: n_selected, Nint, iproc
integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected)
integer :: i,j,k
integer :: new_size
PROVIDE H_apply_buffer_allocated
new_size = H_apply_buffer(iproc)%N_det + n_selected
if (new_size > H_apply_buffer(iproc)%sze) then
call resize_h_apply_buffer(max(2*H_apply_buffer(iproc)%sze,new_size),iproc)
endif
do i=1,H_apply_buffer(iproc)%N_det
ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,1,i)) )== elec_alpha_num)
ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,2,i))) == elec_beta_num)
enddo
do i=1,n_selected
do j=1,N_int
H_apply_buffer(iproc)%det(j,1,i+H_apply_buffer(iproc)%N_det) = det_buffer(j,1,i)
H_apply_buffer(iproc)%det(j,2,i+H_apply_buffer(iproc)%N_det) = det_buffer(j,2,i)
enddo
ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,1,i+H_apply_buffer(iproc)%N_det)) )== elec_alpha_num)
ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,2,i+H_apply_buffer(iproc)%N_det))) == elec_beta_num)
H_apply_buffer(iproc)%coef(i,:) = 0.d0
enddo
H_apply_buffer(iproc)%N_det = new_size
do i=1,H_apply_buffer(iproc)%N_det
ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,1,i)) )== elec_alpha_num)
ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,2,i))) == elec_beta_num)
enddo
end
subroutine H_apply_cisd
implicit none
BEGIN_DOC
! Calls H_apply on the HF determinant and selects all connected single and double
! excitations (of the same symmetry).
END_DOC
integer(bit_kind) :: hole_mask(N_int,2)
integer(bit_kind) :: particle_mask(N_int,2)
ASSERT (N_det_generators == 1)
PROVIDE H_apply_buffer_allocated mo_bielec_integrals_in_map
call H_apply_cisd_OpenMP_monoexc(HF_bitmask, &
generators_bitmask(:,:,1,1), generators_bitmask(:,:,2,1))
call H_apply_cisd_OpenMP_diexc(HF_bitmask, &
generators_bitmask(:,:,1,1), generators_bitmask(:,:,2,1), &
generators_bitmask(:,:,1,1), generators_bitmask(:,:,2,1) )
call copy_h_apply_buffer_to_wf
end

21
src/CISD/README.rst
View File

@ -1,3 +1,14 @@
CISD
====
This is a test directory which builds a CISD by setting the follwoing rules:
* The only generator determinant is the Hartee-Fock (single-reference method)
* All generated determinants are included in the wave function (no perturbative
selection)
These rules are set in the ``H_apply.irp.f`` file.
Needed Modules
==============
@ -24,13 +35,6 @@ Documentation
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
`fill_h_apply_buffer <http://github.com/LCPQ/quantum_package/tree/master/src/CISD/H_apply.irp.f#L6>`_
Fill the H_apply buffer with determiants for CISD
`h_apply_cisd <http://github.com/LCPQ/quantum_package/tree/master/src/CISD/H_apply.irp.f#L43>`_
Calls H_apply on the HF determinant and selects all connected single and double
excitations (of the same symmetry).
`cisd_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/CISD/SC2.irp.f#L1>`_
CISD+SC2 method :: take off all the disconnected terms of a CISD (selected or not)
.br
@ -50,5 +54,8 @@ Documentation
`repeat_excitation <http://github.com/LCPQ/quantum_package/tree/master/src/CISD/SC2.irp.f#L143>`_
Undocumented
`cisd <http://github.com/LCPQ/quantum_package/tree/master/src/CISD/cisd_sc2.irp.f#L1>`_
Undocumented

42
src/CISD/cisd.irp.f
View File

@ -1,46 +1,18 @@
program cisd
implicit none
integer :: i,k
double precision, allocatable :: eigvalues(:),eigvectors(:,:)
PROVIDE ref_bitmask_energy H_apply_buffer_allocated mo_bielec_integrals_in_map
integer :: i
double precision :: pt2(10), norm_pert(10), H_pert_diag
double precision,allocatable :: H_jj(:)
double precision :: diag_h_mat_elem
integer,allocatable :: iorder(:)
! N_states = 3
! touch N_states
call H_apply_cisd
allocate(eigvalues(n_states),eigvectors(n_det,n_states))
print *, 'N_det = ', N_det
print *, 'HF = ', HF_energy
print *, 'N_states = ', N_states
psi_coef = - 1.d-4
allocate(H_jj(n_det),iorder(n_det))
do i = 1, N_det
H_jj(i) = diag_h_mat_elem(psi_det(1,1,i),N_int)
iorder(i) = i
enddo
call dsort(H_jj,iorder,n_det)
do k=1,N_states
psi_coef(iorder(k),k) = 1.d0
enddo
call davidson_diag(psi_det,psi_coef,eigvalues,size(psi_coef,1),N_det,N_states,N_int,output_CISD)
do i = 1, N_states
print*,'eigvalues(i) = ',eigvalues(i)
call H_apply_cisd
print *, 'N_det = ', N_det
do i = 1,N_states
print *, 'energy = ',CI_energy(i)
print *, 'E_corr = ',CI_electronic_energy(i) - ref_bitmask_energy
enddo
print *, '---'
print *, 'HF:', HF_energy
print *, '---'
do i = 1,1
print *, 'energy(i) = ',eigvalues(i) + nuclear_repulsion
print *, 'E_corr = ',eigvalues(i) - ref_bitmask_energy
enddo
! call CISD_SC2(psi_det,psi_coef,eigvalues,size(psi_coef,1),N_det,N_states,N_int)
! do i = 1, N_states
! print*,'eigvalues(i) = ',eigvalues(i)
! enddo
deallocate(eigvalues,eigvectors)
end

28
src/CISD/cisd_lapack.irp.f
View File

@ -1,21 +1,21 @@
program cisd
implicit none
integer :: i
PROVIDE ref_bitmask_energy H_apply_buffer_allocated mo_bielec_integrals_in_map
N_states = 3
diag_algorithm = "Lapack"
touch N_states diag_algorithm
print *, 'HF = ', HF_energy
print *, 'N_states = ', N_states
call H_apply_cisd
double precision, allocatable :: eigvalues(:),eigvectors(:,:)
allocate(eigvalues(n_det),eigvectors(n_det,n_det))
print *, 'N_det = ', N_det
call lapack_diag(eigvalues,eigvectors,H_matrix_all_dets,n_det,n_det)
! print *, H_matrix_all_dets
print *, '---'
print *, 'HF:', HF_energy
print *, '---'
do i = 1,20
print *, 'energy(i) = ',eigvalues(i) + nuclear_repulsion
do i = 1,N_states
print *, 'energy = ',CI_energy(i)
print *, 'E_corr = ',CI_electronic_energy(i) - ref_bitmask_energy
enddo
! print *, eigvectors(:,1)
deallocate(eigvalues,eigvectors)
end
! call CISD_SC2(psi_det,psi_coef,eigvalues,size(psi_coef,1),N_det,N_states,N_int)
! do i = 1, N_states
! print*,'eigvalues(i) = ',eigvalues(i)
! enddo
end

57
src/CISD/cisd_sc2.irp.f
View File

@ -1,46 +1,23 @@
program cisd
implicit none
integer :: i,k
double precision, allocatable :: eigvalues(:),eigvectors(:,:)
PROVIDE ref_bitmask_energy H_apply_buffer_allocated mo_bielec_integrals_in_map
integer :: i, j
double precision :: pt2(10), norm_pert(10), H_pert_diag
double precision,allocatable :: H_jj(:)
double precision :: diag_h_mat_elem
integer,allocatable :: iorder(:)
! N_states = 3
! touch N_states
call H_apply_cisd
allocate(eigvalues(n_states),eigvectors(n_det,n_states))
print *, 'N_det = ', N_det
print *, 'HF = ', HF_energy
print *, 'N_states = ', N_states
psi_coef = - 1.d-4
allocate(H_jj(n_det),iorder(n_det))
do i = 1, N_det
H_jj(i) = diag_h_mat_elem(psi_det(1,1,i),N_int)
iorder(i) = i
enddo
call dsort(H_jj,iorder,n_det)
call H_apply_cisd
print *, 'N_det = ', N_det
do i = 1,N_states
print *, 'energy = ',CI_energy(i)
print *, 'E_corr = ',CI_electronic_energy(i) - ref_bitmask_energy
do j=1,N_det
psi_coef(j,i) = CI_eigenvectors(j,i)
enddo
enddo
SOFT_TOUCH CI_electronic_energy CI_eigenvectors
call CISD_SC2(psi_det,psi_coef,CI_electronic_energy,size(psi_coef,1),N_det,N_states,N_int)
TOUCH CI_electronic_energy CI_eigenvectors
do k=1,N_states
psi_coef(iorder(k),k) = 1.d0
enddo
call davidson_diag(psi_det,psi_coef,eigvalues,size(psi_coef,1),N_det,N_states,N_int,output_CISD)
do i = 1, N_states
print*,'eigvalues(i) = ',eigvalues(i)
enddo
print *, '---'
print *, 'HF:', HF_energy
print *, '---'
do i = 1,1
print *, 'energy(i) = ',eigvalues(i) + nuclear_repulsion
print *, 'E_corr = ',eigvalues(i) - ref_bitmask_energy
enddo
call CISD_SC2(psi_det,psi_coef,eigvalues,size(psi_coef,1),N_det,N_states,N_int,output_CISD)
do i = 1, N_states
print*,'eigvalues(i) = ',eigvalues(i)
enddo
deallocate(eigvalues,eigvectors)
do i = 1, N_states
print*,'eigvalues(i) = ',CI_energy(i)
enddo
end

11
src/CISD/h_apply.py
View File

@ -1,11 +0,0 @@
#!/usr/bin/env python
from generate_h_apply import *
s = H_apply("cisd",openmp=True)
s["keys_work"] += """
call fill_H_apply_buffer(key_idx,keys_out,N_int,iproc)
"""
print s

41
src/Dets/H_apply.irp.f
View File

@ -194,8 +194,45 @@ subroutine copy_H_apply_buffer_to_wf
end
subroutine fill_H_apply_buffer_no_selection(n_selected,det_buffer,Nint,iproc)
use bitmasks
implicit none
BEGIN_DOC
! Fill the H_apply buffer with determiants for CISD
END_DOC
integer, intent(in) :: n_selected, Nint, iproc
integer(bit_kind), intent(in) :: det_buffer(Nint,2,n_selected)
integer :: i,j,k
integer :: new_size
PROVIDE H_apply_buffer_allocated
new_size = H_apply_buffer(iproc)%N_det + n_selected
if (new_size > H_apply_buffer(iproc)%sze) then
call resize_h_apply_buffer(max(2*H_apply_buffer(iproc)%sze,new_size),iproc)
endif
do i=1,H_apply_buffer(iproc)%N_det
ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,1,i)) )== elec_alpha_num)
ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,2,i))) == elec_beta_num)
enddo
do i=1,n_selected
do j=1,N_int
H_apply_buffer(iproc)%det(j,1,i+H_apply_buffer(iproc)%N_det) = det_buffer(j,1,i)
H_apply_buffer(iproc)%det(j,2,i+H_apply_buffer(iproc)%N_det) = det_buffer(j,2,i)
enddo
ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,1,i+H_apply_buffer(iproc)%N_det)) )== elec_alpha_num)
ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,2,i+H_apply_buffer(iproc)%N_det))) == elec_beta_num)
enddo
do j=1,N_states
do i=1,N_selected
H_apply_buffer(iproc)%coef(i,j) = 0.d0
enddo
enddo
H_apply_buffer(iproc)%N_det = new_size
do i=1,H_apply_buffer(iproc)%N_det
ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,1,i)) )== elec_alpha_num)
ASSERT (sum(popcnt(H_apply_buffer(iproc)%det(:,2,i))) == elec_beta_num)
enddo
end

99
src/Dets/H_apply_template.f
View File

@ -37,7 +37,7 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame
iproc = 0
$omp_parallel
!$ iproc = omp_get_thread_num()
!$ iproc = omp_get_thread_num()
allocate (keys_out(N_int,2,size_max),hij_tab(size_max))
!print*,'key_in !!'
@ -69,7 +69,7 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame
i_a = occ_hole(ii,ispin)
ASSERT (i_a > 0)
ASSERT (i_a <= mo_tot_num)
do jj=1,N_elec_in_key_part_1(ispin) !particle
j_a = occ_particle(jj,ispin)
ASSERT (j_a > 0)
@ -142,20 +142,20 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame
k = ishft(j_b-1,-bit_kind_shift)+1
l = j_b-ishft(k-1,bit_kind_shift)-1
key(k,other_spin) = ibset(key(k,other_spin),l)
! call i_H_j(key,key_in,N_int,hij_elec)
! if(dabs(hij_elec)>=thresh) then
key_idx += 1
do k=1,N_int
keys_out(k,1,key_idx) = key(k,1)
keys_out(k,2,key_idx) = key(k,2)
enddo
hij_tab(key_idx) = hij_elec
ASSERT (key_idx <= size_max)
if (key_idx == size_max) then
$keys_work
key_idx = 0
endif
! endif
! call i_H_j(key,key_in,N_int,hij_elec)
! if(dabs(hij_elec)>=thresh) then
key_idx += 1
do k=1,N_int
keys_out(k,1,key_idx) = key(k,1)
keys_out(k,2,key_idx) = key(k,2)
enddo
hij_tab(key_idx) = hij_elec
ASSERT (key_idx <= size_max)
if (key_idx == size_max) then
$keys_work
key_idx = 0
endif
! endif
enddo
enddo
endif
@ -181,20 +181,20 @@ subroutine $subroutine_diexc(key_in, hole_1,particl_1, hole_2, particl_2 $parame
key(k,ispin) = ibset(key(k,ispin),l)
!! a^((+)_j_b(ispin) a_i_b(ispin) : mono exc :: orb(i_b,ispin) --> orb(j_b,ispin)
! call i_H_j(key,key_in,N_int,hij_elec)
! if(dabs(hij_elec)>=thresh) then
key_idx += 1
do k=1,N_int
keys_out(k,1,key_idx) = key(k,1)
keys_out(k,2,key_idx) = key(k,2)
enddo
hij_tab(key_idx) = hij_elec
ASSERT (key_idx <= size_max)
if (key_idx == size_max) then
$keys_work
key_idx = 0
endif
! endif
! call i_H_j(key,key_in,N_int,hij_elec)
! if(dabs(hij_elec)>=thresh) then
key_idx += 1
do k=1,N_int
keys_out(k,1,key_idx) = key(k,1)
keys_out(k,2,key_idx) = key(k,2)
enddo
hij_tab(key_idx) = hij_elec
ASSERT (key_idx <= size_max)
if (key_idx == size_max) then
$keys_work
key_idx = 0
endif
! endif
enddo
enddo! kk
enddo ! ii
@ -246,7 +246,7 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1 $parameters )
iproc = 0
$omp_parallel
!$ iproc = omp_get_thread_num()
!$ iproc = omp_get_thread_num()
allocate (keys_out(N_int,2,size_max),hij_tab(size_max))
!!!! First couple hole particle
do j = 1, N_int
@ -255,13 +255,13 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1 $parameters )
particle(j,1) = iand(xor(particl_1(j,1),key_in(j,1)),particl_1(j,1))
particle(j,2) = iand(xor(particl_1(j,2),key_in(j,2)),particl_1(j,2))
enddo
call bitstring_to_list(particle(1,1),occ_particle(1,1),N_elec_in_key_part_1(1),N_int)
call bitstring_to_list(particle(1,2),occ_particle(1,2),N_elec_in_key_part_1(2),N_int)
call bitstring_to_list(hole (1,1),occ_hole (1,1),N_elec_in_key_hole_1(1),N_int)
call bitstring_to_list(hole (1,2),occ_hole (1,2),N_elec_in_key_hole_1(2),N_int)
allocate (ia_ja_pairs(2,0:(elec_alpha_num)*mo_tot_num,2))
do ispin=1,2
i=0
do ii=N_elec_in_key_hole_1(ispin),1,-1 ! hole
@ -315,3 +315,36 @@ subroutine $subroutine_monoexc(key_in, hole_1,particl_1 $parameters )
end
subroutine $subroutine($params_main)
implicit none
use bitmasks
BEGIN_DOC
! Calls H_apply on the HF determinant and selects all connected single and double
! excitations (of the same symmetry). Auto-generated by the ``generate_h_apply`` script.
END_DOC
$decls_main
PROVIDE H_apply_buffer_allocated mo_bielec_integrals_in_map N_det_reference psi_generators
integer :: imask
do imask=1,N_det_generators
call $subroutine_monoexc(psi_generators(1,1,imask), &
generators_bitmask(1,1,s_hole ,i_bitmask_gen), &
generators_bitmask(1,1,s_part ,i_bitmask_gen) &
$params_post)
call $subroutine_diexc(psi_generators(1,1,imask), &
generators_bitmask(1,1,d_hole1,i_bitmask_gen), &
generators_bitmask(1,1,d_part1,i_bitmask_gen), &
generators_bitmask(1,1,d_hole2,i_bitmask_gen), &
generators_bitmask(1,1,d_part2,i_bitmask_gen) &
$params_post)
enddo
$copy_buffer
$generate_psi_guess
SOFT_TOUCH psi_det psi_coef
end

3
src/Dets/README.rst
View File

@ -53,6 +53,9 @@ Documentation
`copy_h_apply_buffer_to_wf <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/H_apply.irp.f#L113>`_
Undocumented
`fill_h_apply_buffer_no_selection <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/H_apply.irp.f#L197>`_
Fill the H_apply buffer with determiants for CISD
`h_apply_buffer_allocated <http://github.com/LCPQ/quantum_package/tree/master/src/Dets/H_apply.irp.f#L14>`_
Buffer of determinants/coefficients/perturbative energy for H_apply.
Uninitialized. Filled by H_apply subroutines.

6
src/Dets/determinants.ezfio_config
View File

@ -1,9 +1,9 @@
determinants
N_int integer
bit_kind integer
n_dets integer
n_det integer
n_states integer
psi_coef double precision (determinants_n_dets,determinants_n_states)
psi_det integer (determinants_N_int*determinants_bit_kind/4,2,determinants_n_dets)
psi_coef double precision (determinants_n_det,determinants_n_states)
psi_det integer*8 (determinants_N_int*determinants_bit_kind/8,2,determinants_n_det)
H_apply_threshold double precision

37
src/Dets/determinants.irp.f
View File

@ -5,7 +5,16 @@ BEGIN_PROVIDER [ integer, N_states ]
BEGIN_DOC
! Number of states to consider
END_DOC
N_states = 1
logical :: exists
PROVIDE ezfio_filename
call ezfio_has_determinants_n_states(exists)
if (exists) then
call ezfio_get_determinants_n_states(N_states)
else
N_states = 1
call ezfio_set_determinants_n_states(N_states)
endif
ASSERT (N_states > 0)
END_PROVIDER
BEGIN_PROVIDER [ integer, N_det ]
@ -13,22 +22,33 @@ BEGIN_PROVIDER [ integer, N_det ]
BEGIN_DOC
! Number of determinants in the wave function
END_DOC
N_det = 1
logical :: exists
PROVIDE ezfio_filename
call ezfio_has_determinants_n_det(exists)
if (exists) then
call ezfio_get_determinants_n_det(N_det)
else
N_det = 1
call ezfio_set_determinants_n_det(N_det)
endif
ASSERT (N_det > 0)
END_PROVIDER
BEGIN_PROVIDER [ integer, psi_det_size ]
implicit none
BEGIN_DOC
! Size of the psi_det/psi_coef arrays
END_DOC
psi_det_size = 1000
psi_det_size = 1000*N_states
END_PROVIDER
BEGIN_PROVIDER [ integer(bit_kind), psi_det, (N_int,2,psi_det_size) ]
&BEGIN_PROVIDER [ double precision, psi_coef, (psi_det_size,N_states) ]
implicit none
BEGIN_DOC
! The wave function. Initialized with Hartree-Fock
! The wave function. Initialized with Hartree-Fock if the EZFIO file
! is empty
END_DOC
integer, save :: ifirst = 0
@ -52,3 +72,12 @@ END_PROVIDER
END_PROVIDER
BEGIN_PROVIDER [ integer, N_det_reference ]
implicit none
BEGIN_DOC
! Number of determinants in the reference wave function
END_DOC
N_det_reference = N_det
ASSERT (N_det_reference > 0)
END_PROVIDER

64
src/Dets/diagonalize_CI.irp.f Normal file
View File

@ -0,0 +1,64 @@
BEGIN_PROVIDER [ character*(64), diag_algorithm ]
implicit none
BEGIN_DOC
! Diagonalization algorithm (Davidson or Lapack)
END_DOC
if (N_det > 500) then
diag_algorithm = "Davidson"
else
diag_algorithm = "Lapack"
endif
END_PROVIDER
BEGIN_PROVIDER [ double precision, CI_energy, (N_states) ]
implicit none
BEGIN_DOC
! N_states lowest eigenvalues of the CI matrix
END_DOC
integer :: j
do j=1,min(N_states,N_det)
CI_energy(j) = CI_electronic_energy(j) + nuclear_repulsion
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, CI_electronic_energy, (N_states) ]
&BEGIN_PROVIDER [ double precision, CI_eigenvectors, (N_det,N_states) ]
implicit none
BEGIN_DOC
! Eigenvectors/values of the CI matrix
END_DOC
integer :: i,j
do j=1,N_states
do i=1,N_det
CI_eigenvectors(i,j) = psi_coef(i,j)
enddo
enddo
if (diag_algorithm == "Davidson") then
call davidson_diag(psi_det,CI_eigenvectors,CI_electronic_energy, &
size(CI_eigenvectors,1),N_det,N_states,N_int,output_Dets)
else if (diag_algorithm == "Lapack") then
double precision, allocatable :: eigenvectors(:,:), eigenvalues(:)
allocate (eigenvectors(size(H_matrix_all_dets,1),N_det))
allocate (eigenvalues(N_det))
call lapack_diag(eigenvalues,eigenvectors, &
H_matrix_all_dets,size(H_matrix_all_dets,1),N_det)
CI_electronic_energy(:) = 0.d0
do j=1,min(N_states,N_det)
do i=1,N_det
CI_eigenvectors(i,j) = eigenvectors(i,j)
enddo
CI_electronic_energy(j) = eigenvalues(j)
enddo
deallocate(eigenvectors,eigenvalues)
endif
END_PROVIDER

3
src/Dets/utils.irp.f
View File

@ -1,4 +1,4 @@
BEGIN_PROVIDER [ double precision, H_matrix_all_dets,(n_det,n_det) ]
BEGIN_PROVIDER [ double precision, H_matrix_all_dets,(N_det,N_det) ]
implicit none
BEGIN_DOC
! H matrix on the basis of the slater deter;inants defined by psi_det
@ -13,3 +13,4 @@ BEGIN_PROVIDER [ double precision, H_matrix_all_dets,(n_det,n_det) ]
enddo
enddo
END_PROVIDER

30
src/Perturbation/README.rst
View File

@ -82,8 +82,15 @@ Documentation
.. Do not edit this section. It was auto-generated from the
.. NEEDED_MODULES file.
`cisd <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/cisd_test.irp.f#L1>`_
Undocumented
`pt2_moller_plesset <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/Moller_plesset.irp.f#L1>`_
compute the standard Moller-Plesset perturbative first order coefficient and second order energetic contribution
.br
for the various n_st states.
.br
c_pert(i) = <psi(i)|H|det_pert>/(difference of orbital energies)
.br
e_2_pert(i) = <psi(i)|H|det_pert>^2/(difference of orbital energies)
.br
`pt2_epstein_nesbet <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/epstein_nesbet.irp.f#L1>`_
compute the standard Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
@ -95,7 +102,7 @@ Documentation
e_2_pert(i) = <psi(i)|H|det_pert>^2/( E(i) - <det_pert|H|det_pert> )
.br
`pt2_epstein_nesbet_2x2 <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/epstein_nesbet.irp.f#L34>`_
`pt2_epstein_nesbet_2x2 <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/epstein_nesbet.irp.f#L35>`_
compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution
.br
for the various n_st states.
@ -105,7 +112,7 @@ Documentation
c_pert(i) = e_2_pert(i)/ <psi(i)|H|det_pert>
.br
`pt2_epstein_nesbet_2x2_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/epstein_nesbet.irp.f#L118>`_
`pt2_epstein_nesbet_2x2_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/epstein_nesbet.irp.f#L119>`_
compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution
.br
for the various n_st states.
@ -127,7 +134,7 @@ Documentation
c_pert(i) = e_2_pert(i)/ <psi(i)|H|det_pert>
.br
`pt2_epstein_nesbet_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/epstein_nesbet.irp.f#L68>`_
`pt2_epstein_nesbet_sc2 <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/epstein_nesbet.irp.f#L69>`_
compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
.br
for the various n_st states,
@ -149,7 +156,7 @@ Documentation
e_2_pert(i) = <psi(i)|H|det_pert>^2/( E(i) - (<det_pert|H|det_pert> ) )
.br
`pt2_epstein_nesbet_sc2_projected <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/epstein_nesbet.irp.f#L170>`_
`pt2_epstein_nesbet_sc2_projected <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/epstein_nesbet.irp.f#L171>`_
compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
.br
for the various n_st states,
@ -175,6 +182,9 @@ Documentation
NOTE :::: if you satisfy Brillouin Theorem, the singles don't contribute !!
.br
`cisd <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/perturbation_test.irp.f#L1>`_
Undocumented
`fill_h_apply_buffer_selection <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/selection.irp.f#L1>`_
Fill the H_apply buffer with determiants for the selection
@ -190,14 +200,6 @@ Documentation
`diagonalize <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/temporary_stuff.irp.f#L18>`_
Undocumented
`h_apply_cisd_pt2 <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/temporary_stuff.irp.f#L63>`_
Calls H_apply on the HF determinant and selects all connected single and double
excitations (of the same symmetry).
`h_apply_cisd_selection <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/temporary_stuff.irp.f#L97>`_
Calls H_apply on the HF determinant and selects all connected single and double
excitations (of the same symmetry).
`n_det_ref <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation/temporary_stuff.irp.f#L36>`_
Undocumented

7
src/Perturbation/epstein_nesbet.irp.f
View File

@ -24,6 +24,7 @@ subroutine pt2_epstein_nesbet(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet,n_s
ASSERT (Nint > 0)
call i_H_psi(det_pert,psi_ref,psi_ref_coef,Nint,ndet,psi_ref_size,n_st,i_H_psi_array)
H_pert_diag = diag_H_mat_elem(det_pert,Nint)
print *, H_pert_diag
do i =1,n_st
c_pert(i) = i_H_psi_array(i) / (reference_energy(i) - H_pert_diag + eps)
e_2_pert(i) = c_pert(i) * i_H_psi_array(i)
@ -72,7 +73,7 @@ subroutine pt2_epstein_nesbet_SC2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,ndet
integer(bit_kind), intent(in) :: det_pert(Nint,2)
double precision , intent(out) :: c_pert(n_st),e_2_pert(n_st),H_pert_diag
double precision :: i_H_psi_array(N_st)
integer :: idx_repeat(ndet)
integer :: idx_repeat(0:ndet)
BEGIN_DOC
! compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution
@ -122,7 +123,7 @@ subroutine pt2_epstein_nesbet_2x2_SC2(det_pert,c_pert,e_2_pert,H_pert_diag,Nint,
integer(bit_kind), intent(in) :: det_pert(Nint,2)
double precision , intent(out) :: c_pert(n_st),e_2_pert(n_st),H_pert_diag
double precision :: i_H_psi_array(N_st)
integer :: idx_repeat(ndet)
integer :: idx_repeat(0:ndet)
BEGIN_DOC
! compute the Epstein-Nesbet 2x2 diagonalization coefficient and energetic contribution
@ -174,7 +175,7 @@ subroutine pt2_epstein_nesbet_SC2_projected(det_pert,c_pert,e_2_pert,H_pert_diag
integer(bit_kind), intent(in) :: det_pert(Nint,2)
double precision , intent(out) :: c_pert(n_st),e_2_pert(n_st),H_pert_diag
double precision :: i_H_psi_array(N_st)
integer :: idx_repeat(ndet)
integer :: idx_repeat(0:ndet)
BEGIN_DOC
! compute the Epstein-Nesbet perturbative first order coefficient and second order energetic contribution

2
src/Perturbation/perturbation_template.f
View File

@ -24,7 +24,7 @@ subroutine perturb_buffer_$PERT(buffer,buffer_size,e_2_pert_buffer,coef_pert_buf
ASSERT (N_st > 0)
do i = 1,buffer_size
c_ref = connected_to_ref(buffer(1,1,i),psi_det,Nint,N_det_generators,N_det,h_apply_threshold)
c_ref = connected_to_ref(buffer(1,1,i),psi_det,Nint,N_det_generators,N_det_reference,h_apply_threshold)
if (c_ref /= 0) then
cycle

6
src/Perturbation/perturbation_test.irp.f
View File

@ -2,11 +2,6 @@ program cisd
implicit none
integer :: i,k
double precision, allocatable :: eigvalues(:),eigvectors(:,:)
PROVIDE ref_bitmask_energy H_apply_buffer_allocated mo_bielec_integrals_in_map
! N_states = 8
! TOUCH N_states
double precision, allocatable :: pt2(:), norm_pert(:)
double precision :: H_pert_diag
integer :: N_st
@ -14,7 +9,6 @@ program cisd
allocate (pt2(N_st), norm_pert(N_st))
call H_apply_cisd_pt2(pt2, norm_pert, H_pert_diag, N_st)
allocate(eigvalues(n_states),eigvectors(n_det,n_states))
print *, 'N_det = ', N_det
print *, 'N_states = ', N_states
print *, 'pt2 = ', pt2

72
src/Perturbation/temporary_stuff.irp.f
View File

@ -60,75 +60,3 @@ END_PROVIDER
END_PROVIDER
subroutine H_apply_cisd_pt2(pt2, norm_pert, H_pert_diag, N_st)
implicit none
BEGIN_DOC
! Calls H_apply on the HF determinant and selects all connected single and double
! excitations (of the same symmetry).
END_DOC
integer, intent(in) :: N_st
double precision, intent(inout):: pt2(N_st)
double precision, intent(inout):: norm_pert(N_st)
double precision, intent(inout):: H_pert_diag
integer(bit_kind) :: hole_mask(N_int,2)
integer(bit_kind) :: particle_mask(N_int,2)
hole_mask(:,1) = HF_bitmask(:,1)
hole_mask(:,2) = HF_bitmask(:,2)
particle_mask(:,1) = iand(not(HF_bitmask(:,1)),full_ijkl_bitmask(:,1))
particle_mask(:,2) = iand(not(HF_bitmask(:,2)),full_ijkl_bitmask(:,2))
PROVIDE reference_energy N_det_generators key_pattern_not_in_ref
PROVIDE mo_bielec_integrals_in_map
PROVIDE H_apply_buffer_allocated
pt2 = 0.d0
call H_apply_cisd_pt2_OpenMP_monoexc(HF_bitmask, &
hole_mask, particle_mask, &
pt2,norm_pert,H_pert_diag,N_st,N_int )
call H_apply_cisd_pt2_OpenMP_diexc(HF_bitmask, &
hole_mask, particle_mask, &
hole_mask, particle_mask, &
pt2,norm_pert,H_pert_diag,N_st,N_int )
end
subroutine H_apply_cisd_selection(pt2, norm_pert, H_pert_diag, N_st)
implicit none
BEGIN_DOC
! Calls H_apply on the HF determinant and selects all connected single and double
! excitations (of the same symmetry).
END_DOC
integer, intent(in) :: N_st
double precision, intent(inout):: pt2(N_st)
double precision, intent(inout):: norm_pert(N_st)
double precision, intent(inout):: H_pert_diag
integer(bit_kind) :: hole_mask(N_int,2)
integer(bit_kind) :: particle_mask(N_int,2)
hole_mask(:,1) = HF_bitmask(:,1)
hole_mask(:,2) = HF_bitmask(:,2)
particle_mask(:,1) = iand(not(HF_bitmask(:,1)),full_ijkl_bitmask(:,1))
particle_mask(:,2) = iand(not(HF_bitmask(:,2)),full_ijkl_bitmask(:,2))
PROVIDE reference_energy N_det_generators key_pattern_not_in_ref
PROVIDE mo_bielec_integrals_in_map selection_criterion
PROVIDE H_apply_buffer_allocated
pt2 = 0.d0
norm_pert = 0.d0
H_pert_diag = 0.d0
call H_apply_cisd_selection_OpenMP_monoexc(HF_bitmask, &
hole_mask, particle_mask, &
pt2,norm_pert,H_pert_diag,N_st,N_int )
call H_apply_cisd_selection_OpenMP_diexc(HF_bitmask, &
hole_mask, particle_mask, &
hole_mask, particle_mask, &
pt2,norm_pert,H_pert_diag,N_st,N_int )
call copy_h_apply_buffer_to_wf
selection_criterion_min = selection_criterion_min*0.1d0
selection_criterion = selection_criterion_min
end

52
src/Utils/LinearAlgebra.irp.f
View File

@ -196,31 +196,42 @@ subroutine lapack_diag(eigvalues,eigvectors,H,nmax,n)
double precision, intent(in) :: H(nmax,n)
double precision,allocatable :: eigenvalues(:)
double precision,allocatable :: work(:)
integer ,allocatable :: iwork(:)
double precision,allocatable :: A(:,:)
allocate(A(nmax,n),eigenvalues(nmax),work(4*nmax))
integer :: LWORK, info, i,j,l,k
integer :: lwork, info, i,j,l,k, liwork
allocate(A(nmax,n),eigenvalues(n))
A=H
lwork = 2*n*n + 6*n+ 1
liwork = 5*n + 3
allocate (work(lwork),iwork(liwork))
! if (n<30) then
! do i=1,n
! do j=1,n
! print *, j,i, H(j,i)
! enddo
! print *, '---'
! enddo
! print *, '---'
! endif
LWORK = 4*nmax
call dsyev( 'V', 'U', n, A, nmax, eigenvalues, work, LWORK, info )
lwork = -1
liwork = -1
call DSYEVD( 'V', 'U', n, A, nmax, eigenvalues, work, lwork, &
iwork, liwork, info )
if (info < 0) then
print *, irp_here, ': the ',-info,'-th argument had an illegal value'
stop 1
else if (info > 0) then
print *, irp_here, ': the algorithm failed to converge; ',info,' off-diagonal'
print *, 'elements of an intermediate tridiagonal form did not converge to zero.'
stop 1
print *, irp_here, ': DSYEVD: the ',-info,'-th argument had an illegal value'
stop 2
endif
lwork = int( work( 1 ) )
liwork = iwork(1)
deallocate (work,iwork)
allocate (work(lwork),iwork(liwork))
call DSYEVD( 'V', 'U', n, A, nmax, eigenvalues, work, lwork, &
iwork, liwork, info )
deallocate(work,iwork)
if (info < 0) then
print *, irp_here, ': DSYEVD: the ',-info,'-th argument had an illegal value'
stop 2
else if( info > 0 ) then
write(*,*)'DSYEVD Failed'
stop 1
end if
eigvectors = 0.d0
eigvalues = 0.d0
do j = 1, n
@ -229,4 +240,5 @@ subroutine lapack_diag(eigvalues,eigvectors,H,nmax,n)
eigvectors(i,j) = A(i,j)
enddo
enddo
deallocate(A,eigenvalues)
end