LCPQ/quantum_package
10
0
Fork 0
mirror of https://github.com/LCPQ/quantum_package synced 2024-07-22 18:57:31 +02:00
Code Issues Releases Wiki Activity

Merge pull request #18 from scemama/master

Browse Source 
Improved parallel Davidson
This commit is contained in:
garniron 2016-10-12 09:43:46 +02:00 committed by GitHub
commit e1e896dab3
40 changed files with 842 additions and 871 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

2
.travis.yml
View File

@ -24,7 +24,7 @@ python:
script:
- ./configure --production ./config/gfortran.cfg
- source ./quantum_package.rc ; qp_module.py install Full_CI Hartree_Fock CAS_SD MRCC_CASSD All_singles
- source ./quantum_package.rc ; qp_module.py install Full_CI Full_CI_ZMQ Hartree_Fock CAS_SD mrcepa0 All_singles
- source ./quantum_package.rc ; ninja
- source ./quantum_package.rc ; cd ocaml ; make ; cd -
- source ./quantum_package.rc ; cd tests ; ./run_tests.sh #-v

5
README.md
View File

@ -7,11 +7,14 @@ Set of quantum chemistry programs and libraries.
For more information, you can visit the [wiki of the project](http://github.com/LCPQ/quantum_package/wiki>), or below for the installation instructions.
Demo
====
[![Full-CI energy of C2 in 2 minutes](https://i.vimeocdn.com/video/555047954_295x166.jpg)](https://vimeo.com/scemama/quantum_package_demo "Quantum Package Demo")
[![Frozen-core Full-CI energy of Ti](https://raw.githubusercontent.com/LCPQ/quantum_package/master/data/Titanium.png)](https://raw.githubusercontent.com/LCPQ/quantum_package/master/data/Titanium.png "Convergence of Ti in cc-pv{DTQ}Z")
# Installation
@ -155,7 +158,7 @@ Program exited with code 139.
#### Why ?
It's caused when we call the DGEM routine of LAPACK.
It's caused when we call the DGEMM routine of LAPACK.
##### Fix

BIN
data/Titanium.png Normal file
View File

Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 54 KiB

4
plugins/FOBOCI/H_apply_dressed_autonom.irp.f
View File

@ -273,7 +273,7 @@ subroutine H_apply_dressed_pert_monoexc(key_in, hole_1,particl_1,i_generator,ipr
integer,parameter :: size_max = 3072
integer, intent(in) :: Ndet_generators
double precision, intent(in) :: E_ref
double precision, intent(inout) :: E_ref
double precision, intent(inout) :: delta_ij_generators_(Ndet_generators,Ndet_generators)
integer(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,Ndet_generators)
@ -438,7 +438,7 @@ subroutine H_apply_dressed_pert(delta_ij_generators_, Ndet_generators,psi_det_g
integer, intent(in) :: Ndet_generators
double precision, intent(in) :: E_ref
double precision, intent(inout) :: E_ref
double precision, intent(inout) :: delta_ij_generators_(Ndet_generators,Ndet_generators)
integer(bit_kind), intent(in) :: psi_det_generators_input(N_int,2,Ndet_generators)

2
plugins/FOBOCI/NEEDED_CHILDREN_MODULES
View File

@ -1 +1 @@
Perturbation Selectors_no_sorted Hartree_Fock Davidson
Perturbation Selectors_no_sorted Hartree_Fock Davidson CISD

7
plugins/FOBOCI/hcc_1h1p.irp.f
View File

@ -15,11 +15,10 @@ subroutine routine
call diagonalize_CI
call test_hcc
call test_mulliken
! call SC2_1h1p(psi_det,psi_coef,energies, &
! diag_H_elements,size(psi_coef,1),N_det,N_states_diag,N_int,threshold_convergence_SC2)
allocate(H_matrix(N_det,N_det))
call SC2_1h1p_full(psi_det,psi_coef,energies, &
H_matrix,size(psi_coef,1),N_det,N_states_diag,N_int,threshold_convergence_SC2)
stop 'SC2_1h1p_full is not in the git!'
! call SC2_1h1p_full(psi_det,psi_coef,energies, &
! H_matrix,size(psi_coef,1),N_det,N_states_diag,N_int,threshold_convergence_SC2)
deallocate(H_matrix)
integer :: i,j
double precision :: accu,coef_hf

29
plugins/Full_CI_ZMQ/fci_zmq.irp.f
View File

@ -46,14 +46,16 @@ program fci_zmq
PROVIDE psi_det_sorted
call diagonalize_CI
call save_wavefunction
if (N_det > N_det_max) then
psi_det = psi_det_sorted
psi_coef = psi_coef_sorted
N_det = N_det_max
soft_touch N_det psi_det psi_coef
call diagonalize_CI
call save_wavefunction
endif
call save_wavefunction
print *, 'N_det = ', N_det
print *, 'N_states = ', N_states
@ -79,17 +81,13 @@ program fci_zmq
E_CI_before(1:N_states) = CI_energy(1:N_states)
call ezfio_set_full_ci_energy(CI_energy)
enddo
if (N_det > N_det_max) then
N_det = N_det_max
touch N_det psi_det psi_coef
call diagonalize_CI
endif
if(do_pt2_end)then
print*,'Last iteration only to compute the PT2'
threshold_selectors = 1.d0
threshold_generators = 0.9999d0
E_CI_before(1:N_states) = CI_energy(1:N_states)
call ZMQ_selection(1, pt2)
call ZMQ_selection(0, pt2)
print *, 'Final step'
print *, 'N_det = ', N_det
print *, 'N_states = ', N_states
@ -108,7 +106,7 @@ end
subroutine ZMQ_selection(N, pt2)
subroutine ZMQ_selection(N_in, pt2)
use f77_zmq
use selection_types
@ -116,13 +114,14 @@ subroutine ZMQ_selection(N, pt2)
character*(512) :: task
integer(ZMQ_PTR) :: zmq_to_qp_run_socket
integer, intent(in) :: N
integer, intent(in) :: N_in
type(selection_buffer) :: b
integer :: i
integer :: i, N
integer, external :: omp_get_thread_num
double precision, intent(out) :: pt2(N_states)
N = max(N_in,1)
provide nproc
provide ci_electronic_energy
call new_parallel_job(zmq_to_qp_run_socket,"selection")
@ -147,16 +146,18 @@ subroutine ZMQ_selection(N, pt2)
if (i==0) then
call selection_collector(b, pt2)
else
call selection_dressing_slave_inproc(i)
call selection_slave_inproc(i)
endif
!$OMP END PARALLEL
call end_parallel_job(zmq_to_qp_run_socket, 'selection')
call fill_H_apply_buffer_no_selection(b%cur,b%det,N_int,0) !!! PAS DE ROBIN
call copy_H_apply_buffer_to_wf()
if (N_in > 0) then
call fill_H_apply_buffer_no_selection(b%cur,b%det,N_int,0) !!! PAS DE ROBIN
call copy_H_apply_buffer_to_wf()
endif
end subroutine
subroutine selection_dressing_slave_inproc(i)
subroutine selection_slave_inproc(i)
implicit none
integer, intent(in) :: i

107
plugins/Full_CI_ZMQ/selection_davidson_slave.irp.f Normal file
View File

@ -0,0 +1,107 @@
program selection_slave
implicit none
BEGIN_DOC
! Helper program to compute the PT2 in distributed mode.
END_DOC
read_wf = .False.
SOFT_TOUCH read_wf
call provide_everything
call switch_qp_run_to_master
call run_wf
end
subroutine provide_everything
PROVIDE H_apply_buffer_allocated mo_bielec_integrals_in_map psi_det_generators psi_coef_generators psi_det_sorted_bit psi_selectors n_det_generators n_states generators_bitmask zmq_context mo_mono_elec_integral
! PROVIDE ci_electronic_energy mo_tot_num N_int
end
subroutine run_wf
use f77_zmq
implicit none
integer(ZMQ_PTR), external :: new_zmq_to_qp_run_socket
integer(ZMQ_PTR) :: zmq_to_qp_run_socket
double precision :: energy(N_states_diag)
character*(64) :: states(2)
integer :: rc, i
call provide_everything
zmq_context = f77_zmq_ctx_new ()
states(1) = 'selection'
states(2) = 'davidson'
zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
do
call wait_for_states(states,zmq_state,2)
if(trim(zmq_state) == 'Stopped') then
exit
else if (trim(zmq_state) == 'selection') then
! Selection
! ---------
print *, 'Selection'
call zmq_get_psi(zmq_to_qp_run_socket,1,energy,N_states_diag)
!$OMP PARALLEL PRIVATE(i)
i = omp_get_thread_num()
call selection_slave_tcp(i, energy)
!$OMP END PARALLEL
print *, 'Selection done'
else if (trim(zmq_state) == 'davidson') then
! Davidson
! --------
print *, 'Davidson'
call davidson_miniserver_get()
!$OMP PARALLEL PRIVATE(i)
i = omp_get_thread_num()
call davidson_slave_tcp(i)
!$OMP END PARALLEL
print *, 'Davidson done'
endif
end do
end
subroutine update_energy(energy)
implicit none
double precision, intent(in) :: energy(N_states_diag)
BEGIN_DOC
! Update energy when it is received from ZMQ
END_DOC
integer :: j,k
do j=1,N_states
do k=1,N_det
CI_eigenvectors(k,j) = psi_coef(k,j)
enddo
enddo
call u_0_S2_u_0(CI_eigenvectors_s2,CI_eigenvectors,N_det,psi_det,N_int)
if (.True.) then
do k=1,size(ci_electronic_energy)
ci_electronic_energy(k) = energy(k)
enddo
TOUCH ci_electronic_energy CI_eigenvectors_s2 CI_eigenvectors
endif
call write_double(6,ci_energy,'Energy')
end
subroutine selection_slave_tcp(i,energy)
implicit none
double precision, intent(in) :: energy(N_states_diag)
integer, intent(in) :: i
call run_selection_slave(0,i,energy)
end

23
plugins/Full_CI_ZMQ/selection_slave.irp.f
View File

@ -23,20 +23,19 @@ subroutine run_wf
integer(ZMQ_PTR), external :: new_zmq_to_qp_run_socket
integer(ZMQ_PTR) :: zmq_to_qp_run_socket
double precision :: energy(N_states_diag)
character*(64) :: states(2)
character*(64) :: states(1)
integer :: rc, i
call provide_everything
zmq_context = f77_zmq_ctx_new ()
states(1) = 'selection'
states(2) = 'davidson'
zmq_to_qp_run_socket = new_zmq_to_qp_run_socket()
do
call wait_for_states(states,zmq_state,2)
call wait_for_states(states,zmq_state,1)
if(trim(zmq_state) == 'Stopped') then
@ -52,24 +51,10 @@ subroutine run_wf
!$OMP PARALLEL PRIVATE(i)
i = omp_get_thread_num()
call selection_dressing_slave_tcp(i, energy)
call selection_slave_tcp(i, energy)
!$OMP END PARALLEL
print *, 'Selection done'
else if (trim(zmq_state) == 'davidson') then
! Davidson
! --------
print *, 'Davidson'
call davidson_miniserver_get()
!!$OMP PARALLEL PRIVATE(i)
!i = omp_get_thread_num()
call davidson_slave_tcp(0)
!!$OMP END PARALLEL
print *, 'Davidson done'
endif
end do
@ -98,7 +83,7 @@ subroutine update_energy(energy)
call write_double(6,ci_energy,'Energy')
end
subroutine selection_dressing_slave_tcp(i,energy)
subroutine selection_slave_tcp(i,energy)
implicit none
double precision, intent(in) :: energy(N_states_diag)
integer, intent(in) :: i

34
plugins/MRCC_CASSD/.gitignore vendored
View File

@ -1,34 +0,0 @@
# Automatically created by $QP_ROOT/scripts/module/module_handler.py
.ninja_deps
.ninja_log
AO_Basis
Bitmask
Determinants
Electrons
Ezfio_files
Generators_full
Hartree_Fock
IRPF90_man
IRPF90_temp
Integrals_Bielec
Integrals_Monoelec
MOGuess
MO_Basis
MRCC_Utils
Makefile
Makefile.depend
Nuclei
Perturbation
Properties
Pseudo
Psiref_CAS
Psiref_Utils
Selectors_full
Utils
ZMQ
ezfio_interface.irp.f
irpf90.make
irpf90_entities
mrcc_cassd
mrcc_noiter
tags

17
plugins/MRCC_CASSD/EZFIO.cfg
View File

@ -1,17 +0,0 @@
[energy]
type: double precision
doc: Calculated energy
interface: ezfio
[thresh_mrcc]
type: Threshold
doc: Threshold on the convergence of the MRCC energy
interface: ezfio,provider,ocaml
default: 1.e-5
[n_it_mrcc_max]
type: Strictly_positive_int
doc: Maximum number of MRCC iterations
interface: ezfio,provider,ocaml
default: 10

1
plugins/MRCC_CASSD/NEEDED_CHILDREN_MODULES
View File

@ -1 +0,0 @@
Perturbation Selectors_full Generators_full Psiref_CAS MRCC_Utils

80
plugins/MRCC_CASSD/README.rst
View File

@ -1,80 +0,0 @@
===========
MRCC Module
===========
MRCC as a coupled cluster on a CAS+SD wave function.
Needed Modules
==============
.. Do not edit this section. It was auto-generated from the
.. by the `update_README.py` script.
.. image:: tree_dependency.png
* `Perturbation <http://github.com/LCPQ/quantum_package/tree/master/src/Perturbation>`_
* `Selectors_full <http://github.com/LCPQ/quantum_package/tree/master/src/Selectors_full>`_
* `Generators_full <http://github.com/LCPQ/quantum_package/tree/master/src/Generators_full>`_
* `Psiref_CAS <http://github.com/LCPQ/quantum_package/tree/master/src/Psiref_CAS>`_
* `MRCC_Utils <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_Utils>`_
Documentation
=============
.. Do not edit this section. It was auto-generated from the
.. by the `update_README.py` script.
`mrcc <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_CASSD/mrcc_cassd.irp.f#L1>`_
Undocumented
`print_cas_coefs <http://github.com/LCPQ/quantum_package/tree/master/src/MRCC_CASSD/mrcc_cassd.irp.f#L11>`_
Undocumented
Needed Modules
==============
.. Do not edit this section It was auto-generated
.. by the `update_README.py` script.
.. image:: tree_dependency.png
* `Perturbation <http://github.com/LCPQ/quantum_package/tree/master/plugins/Perturbation>`_
* `Selectors_full <http://github.com/LCPQ/quantum_package/tree/master/plugins/Selectors_full>`_
* `Generators_full <http://github.com/LCPQ/quantum_package/tree/master/plugins/Generators_full>`_
* `Psiref_CAS <http://github.com/LCPQ/quantum_package/tree/master/plugins/Psiref_CAS>`_
* `MRCC_Utils <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_Utils>`_
Documentation
=============
.. Do not edit this section It was auto-generated
.. by the `update_README.py` script.
`mrcc <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_CASSD/mrcc_cassd.irp.f#L1>`_
Undocumented
`mrcc_noiter <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_CASSD/mrcc_noiter.irp.f#L1>`_
Undocumented
`n_it_mrcc_max <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_CASSD/ezfio_interface.irp.f#L6>`_
Maximum number of MRCC iterations
`print_cas_coefs <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_CASSD/mrcc_noiter.irp.f#L78>`_
Undocumented
`run <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_CASSD/mrcc_noiter.irp.f#L17>`_
Undocumented
`run_pt2 <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_CASSD/mrcc_noiter.irp.f#L36>`_
Undocumented
`thresh_mrcc <http://github.com/LCPQ/quantum_package/tree/master/plugins/MRCC_CASSD/ezfio_interface.irp.f#L25>`_
Threshold on the convergence of the MRCC energy

120
plugins/MRCC_CASSD/mrcc_cassd.irp.f
View File

@ -1,120 +0,0 @@
program mrcc
implicit none
double precision, allocatable :: energy(:)
allocate (energy(N_states))
read_wf = .True.
SOFT_TOUCH read_wf
call print_cas_coefs
call set_generators_bitmasks_as_holes_and_particles
call run(N_states,energy)
if(do_pt2_end)then
call run_pt2(N_states,energy)
endif
deallocate(energy)
end
subroutine run(N_st,energy)
implicit none
integer, intent(in) :: N_st
double precision, intent(out) :: energy(N_st)
integer :: i
double precision :: E_new, E_old, delta_e
integer :: iteration
double precision :: E_past(4), lambda
E_new = 0.d0
delta_E = 1.d0
iteration = 0
lambda = 1.d0
do while (delta_E > thresh_mrcc)
iteration += 1
print *, '==========================='
print *, 'MRCC Iteration', iteration
print *, '==========================='
print *, ''
E_old = sum(ci_energy_dressed(1:N_st))
call write_double(6,ci_energy_dressed(1),"MRCC energy")
call diagonalize_ci_dressed(lambda)
E_new = sum(ci_energy_dressed(1:N_st))
delta_E = dabs(E_new - E_old)
call save_wavefunction
call ezfio_set_mrcc_cassd_energy(ci_energy_dressed(1))
if (iteration > n_it_mrcc_max) then
exit
endif
enddo
call write_double(6,ci_energy_dressed(1),"Final MRCC energy")
energy(1:N_st) = ci_energy_dressed(1:N_st)
end
subroutine run_pt2(N_st,energy)
implicit none
integer :: i,j,k
double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:)
integer, intent(in) :: N_st
double precision, intent(in) :: energy(N_st)
allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st))
pt2 = 0.d0
print*,'Last iteration only to compute the PT2'
threshold_selectors = 1.d0
threshold_generators = 0.999d0
N_det_generators = lambda_mrcc_pt2(0) + N_det_cas
do i=1,N_det_cas
do k=1,N_int
psi_det_generators(k,1,i) = psi_ref(k,1,i)
psi_det_generators(k,2,i) = psi_ref(k,2,i)
enddo
do k=1,N_st
psi_coef_generators(i,k) = psi_ref_coef(i,k)
enddo
enddo
do i=N_det_cas+1,N_det_generators
j = lambda_mrcc_pt2(i)
do k=1,N_int
psi_det_generators(k,1,i) = psi_non_ref(k,1,j)
psi_det_generators(k,2,i) = psi_non_ref(k,2,j)
enddo
do k=1,N_st
psi_coef_generators(i,k) = psi_non_ref_coef(j,k)
enddo
enddo
SOFT_TOUCH N_det_generators psi_det_generators psi_coef_generators ci_eigenvectors_dressed ci_eigenvectors_s2_dressed ci_electronic_energy_dressed
call H_apply_mrcc_PT2(pt2, norm_pert, H_pert_diag, N_st)
print *, 'Final step'
print *, 'N_det = ', N_det
print *, 'N_states = ', N_states
print *, 'PT2 = ', pt2
print *, 'E = ', energy
print *, 'E+PT2 = ', energy+pt2
print *, '-----'
call ezfio_set_full_ci_energy_pt2(energy+pt2)
deallocate(pt2,norm_pert)
end
subroutine print_cas_coefs
implicit none
integer :: i,j
print *, 'CAS'
print *, '==='
do i=1,N_det_cas
print *, psi_cas_coef(i,:)
call debug_det(psi_cas(1,1,i),N_int)
enddo
call write_double(6,ci_energy(1),"Initial CI energy")
end

91
plugins/MRCC_CASSD/mrcc_noiter.irp.f
View File

@ -1,91 +0,0 @@
program mrcc_noiter
implicit none
double precision, allocatable :: energy(:)
allocate (energy(N_states))
read_wf = .True.
threshold_generators = .9999d0
SOFT_TOUCH read_wf threshold_generators
call print_cas_coefs
call set_generators_bitmasks_as_holes_and_particles
call run(N_states,energy)
if(do_pt2_end)then
call run_pt2(N_states,energy)
endif
deallocate(energy)
end
subroutine run(N_st,energy)
implicit none
integer, intent(in) :: N_st
double precision, intent(out) :: energy(N_st)
integer :: i,j
do j=1,N_states_diag
do i=1,N_det
psi_coef(i,j) = CI_eigenvectors_dressed(i,j)
enddo
enddo
SOFT_TOUCH psi_coef ci_energy_dressed
call write_double(6,ci_energy_dressed(1),"Final MRCC energy")
call ezfio_set_mrcc_cassd_energy(ci_energy_dressed(1))
call save_wavefunction
energy(:) = ci_energy_dressed(:)
end
subroutine run_pt2(N_st,energy)
implicit none
integer :: i,j,k
double precision, allocatable :: pt2(:), norm_pert(:), H_pert_diag(:)
integer, intent(in) :: N_st
double precision, intent(in) :: energy(N_st)
allocate (pt2(N_st), norm_pert(N_st),H_pert_diag(N_st))
pt2 = 0.d0
print*,'Last iteration only to compute the PT2'
threshold_selectors = 1.d0
threshold_generators = 0.999d0
N_det_generators = lambda_mrcc_pt2(0)
do i=1,N_det_generators
j = lambda_mrcc_pt2(i)
do k=1,N_int
psi_det_generators(k,1,i) = psi_non_ref(k,1,j)
psi_det_generators(k,2,i) = psi_non_ref(k,2,j)
enddo
do k=1,N_st
psi_coef_generators(i,k) = psi_non_ref_coef(j,k)
enddo
enddo
SOFT_TOUCH N_det_generators psi_det_generators psi_coef_generators ci_eigenvectors_dressed ci_eigenvectors_s2_dressed ci_electronic_energy_dressed
call H_apply_mrcc_PT2(pt2, norm_pert, H_pert_diag, N_st)
print *, 'Final step'
print *, 'N_det = ', N_det
print *, 'N_states = ', N_states
print *, 'PT2 = ', pt2
print *, 'E = ', energy
print *, 'E+PT2 = ', energy+pt2
print *, '-----'
call ezfio_set_full_ci_energy_pt2(energy+pt2)
deallocate(pt2,norm_pert)
end
subroutine print_cas_coefs
implicit none
integer :: i,j
print *, 'CAS'
print *, '==='
do i=1,N_det_cas
print *, psi_cas_coef(i,:)
call debug_det(psi_cas(1,1,i),N_int)
enddo
call write_double(6,ci_energy(1),"Initial CI energy")
end

BIN
plugins/MRCC_CASSD/tree_dependency.png
View File

Binary file not shown.
Side by Side

Before

Width:  |  Height:  |  Size: 144 KiB

120
plugins/MRCC_Utils/mrcc_utils.irp.f
View File

@ -253,7 +253,7 @@ BEGIN_PROVIDER [ double precision, CI_energy_dressed, (N_states_diag) ]
integer :: j
character*(8) :: st
call write_time(output_determinants)
do j=1,min(N_det,N_states_diag)
do j=1,min(N_det,N_states)
write(st,'(I4)') j
CI_energy_dressed(j) = CI_electronic_energy_dressed(j) + nuclear_repulsion
call write_double(output_determinants,CI_energy_dressed(j),'Energy of state '//trim(st))
@ -349,8 +349,9 @@ integer function searchDet(dets, det, n, Nint)
do while(.true.)
searchDet = (l+h)/2
c = detCmp(dets(1,1,searchDet), det(1,1), Nint)
if(c == 0) return
if(c == 1) then
if(c == 0) then
return
else if(c == 1) then
h = searchDet-1
else
l = searchDet+1
@ -498,12 +499,12 @@ subroutine tamise_exc(key, no, n, N_key)
BEGIN_DOC
! Uncodumented : TODO
END_DOC
integer,intent(in) :: no, n, N_key
integer,intent(in) :: no, n, N_key
integer*2,intent(inout) :: key(4, N_key)
integer :: k,j
integer*2 :: tmp(4)
logical :: exc_inf
integer :: ni
integer :: k,j
integer*2 :: tmp(4)
logical :: exc_inf
integer :: ni
k = no
j = 2*k
@ -633,7 +634,6 @@ END_PROVIDER
allocate(A_ind(0:N_det_ref+1, nex), A_val(N_det_ref+1, nex))
allocate(AtA_ind(N_det_ref * nex), AtA_val(N_det_ref * nex)) !!!!! MAY BE TOO SMALL ? !!!!!!!!
allocate(x(nex), AtB(nex))
allocate(A_val_mwen(nex), A_ind_mwen(nex))
allocate(N_col(nex), col_shortcut(nex))
allocate(x_new(nex))
@ -645,7 +645,6 @@ END_PROVIDER
AtA_val = 0d0
x = 0d0
A_val_mwen = 0d0
A_ind_mwen = 0
N_col = 0
col_shortcut = 0
@ -689,7 +688,8 @@ END_PROVIDER
end do
!$OMP END DO
deallocate(lref)
!$OMP END PARALLEL
!$OMP END PARALLEL
print *, 'Done building A_val, A_ind'
AtB = 0d0
AtA_size = 0
@ -698,6 +698,8 @@ END_PROVIDER
!$OMP PARALLEL default(none) shared(k, psi_non_ref_coef, A_ind, A_val, x, N_det_ref, nex, N_det_non_ref)&
!$OMP private(at_row, a_col, t, i, r1, r2, wk, A_ind_mwen, A_val_mwen)&
!$OMP shared(col_shortcut, N_col, AtB, AtA_size, AtA_val, AtA_ind, s)
allocate(A_val_mwen(nex), A_ind_mwen(nex))
A_ind_mwen = 0
!$OMP DO schedule(dynamic, 100)
do at_row = 1, nex
wk = 0
@ -711,10 +713,10 @@ END_PROVIDER
r1 = 1
r2 = 1
do while ((A_ind(r1, at_row) /= 0).and.(A_ind(r2, a_col) /= 0))
if(A_ind(r1, at_row) < A_ind(r2, a_col)) then
r1 = r1+1
else if(A_ind(r1, at_row) > A_ind(r2, a_col)) then
if(A_ind(r1, at_row) > A_ind(r2, a_col)) then
r2 = r2+1
else if(A_ind(r1, at_row) < A_ind(r2, a_col)) then
r1 = r1+1
else
t = t - A_val(r1, at_row) * A_val(r2, a_col)
r1 = r1+1
@ -748,7 +750,8 @@ END_PROVIDER
!$OMP END CRITICAL
end if
end do
!$OMP END DO
!$OMP END DO NOWAIT
deallocate (A_ind_mwen, A_val_mwen)
!$OMP END PARALLEL
if(AtA_size > size(AtA_val)) stop "SIZA"
@ -807,14 +810,18 @@ END_PROVIDER
if(res < 1d-8) exit
end do
! rho_mrcc now contains A.X
norm = 0.d0
do i=1,N_det_non_ref
norm = norm + rho_mrcc(i,s)*rho_mrcc(i,s)
enddo
! Norm now contains the norm of A.X
do i=1,N_det_ref
norm = norm + psi_ref_coef(i,s)*psi_ref_coef(i,s)
enddo
! Norm now contains the norm of Psi + A.X
print *, k, "res : ", res, "norm : ", sqrt(norm)
@ -826,30 +833,46 @@ END_PROVIDER
if (rho_mrcc(i,s) == 0.d0) then
rho_mrcc(i,s) = 1.d-32
endif
! f is such that f.\tilde{c_i} = c_i
f = psi_non_ref_coef(i,s) / rho_mrcc(i,s)
! Avoid numerical instabilities
f = min(f,2.d0)
f = max(f,-2.d0)
norm = norm + f*f *rho_mrcc(i,s)*rho_mrcc(i,s)
rho_mrcc(i,s) = f
enddo
! norm now contains the norm of |T.Psi_0>
! rho_mrcc now contains the f factors
f = 1.d0/norm
norm = 0.d0
do i=1,N_det_non_ref
norm = norm + psi_non_ref_coef(i,s)*psi_non_ref_coef(i,s)
enddo
f = dsqrt(f*norm)
! f now contains 1/ <T.Psi_0|T.Psi_0>
print *, 'norm of |T Psi_0> = ', norm*f*f
norm = 1.d0
do i=1,N_det_ref
norm = norm - psi_ref_coef(i,s)*psi_ref_coef(i,s)
enddo
! norm now contains <Psi_SD|Psi_SD>
f = dsqrt(f*norm)
! f normalises T.Psi_0 such that (1+T)|Psi> is normalized
norm = norm*f
print *, 'norm of |T Psi_0> = ', dsqrt(norm)
do i=1,N_det_ref
norm = norm + psi_ref_coef(i,s)*psi_ref_coef(i,s)
enddo
do i=1,N_det_non_ref
rho_mrcc(i,s) = rho_mrcc(i,s) * f
enddo
! rho_mrcc now contains the product of the scaling factors and the
! normalization constant
end do
print *, "done"
END_PROVIDER
@ -857,13 +880,17 @@ BEGIN_PROVIDER [ double precision, dij, (N_det_ref, N_det_non_ref, N_states) ]
integer :: s,i,j
double precision, external :: get_dij_index
print *, "computing amplitudes..."
!$OMP PARALLEL DEFAULT(shared) PRIVATE(s,i,j)
do s=1, N_states
!$OMP DO
do i=1, N_det_non_ref
do j=1, N_det_ref
dij(j, i, s) = get_dij_index(j, i, s, N_int)
end do
end do
!$OMP END DO
end do
!$OMP END PARALLEL
print *, "done computing amplitudes"
END_PROVIDER
@ -876,10 +903,9 @@ double precision function get_dij_index(II, i, s, Nint)
double precision :: HIi, phase
if(lambda_type == 0) then
get_dij_index = get_dij(psi_ref(1,1,II), psi_non_ref(1,1,i), s, Nint)
! get_dij_index = get_dij_index * rho_mrcc(i,s)
call get_phase(psi_ref(1,1,II), psi_non_ref(1,1,i), phase, N_int)
get_dij_index = get_dij_index * rho_mrcc(i,s) * phase
get_dij_index = get_dij(psi_ref(1,1,II), psi_non_ref(1,1,i), s, Nint) * phase
get_dij_index = get_dij_index * rho_mrcc(i,s)
else
call i_h_j(psi_ref(1,1,II), psi_non_ref(1,1,i), Nint, HIi)
get_dij_index = HIi * lambda_mrcc(s, i)
@ -1059,18 +1085,20 @@ subroutine apply_hole_local(det, exc, res, ok, Nint)
res = det
if(h1 /= 0) then
ii = (h1-1)/bit_kind_size + 1
pos = mod(h1-1, 64)!iand(h1-1,bit_kind_size-1) ! mod 64
if(iand(det(ii, s1), ishft(1_bit_kind, pos)) == 0_8) return
res(ii, s1) = ibclr(res(ii, s1), pos)
ii = (h1-1)/bit_kind_size + 1
pos = iand(h1-1,bit_kind_size-1) ! mod 64
if(iand(det(ii, s1), ishft(1_bit_kind, pos)) == 0_8) then
return
endif
res(ii, s1) = ibclr(res(ii, s1), pos)
end if
ii = (h2-1)/bit_kind_size + 1
pos = mod(h2-1, 64)!iand(h2-1,bit_kind_size-1)
if(iand(det(ii, s2), ishft(1_bit_kind, pos)) == 0_8) return
res(ii, s2) = ibclr(res(ii, s2), pos)
ii = (h2-1)/bit_kind_size + 1
pos = iand(h2-1,bit_kind_size-1) ! mod 64
if(iand(det(ii, s2), ishft(1_bit_kind, pos)) == 0_8) then
return
endif
res(ii, s2) = ibclr(res(ii, s2), pos)
ok = .true.
end subroutine
@ -1094,16 +1122,20 @@ subroutine apply_particle_local(det, exc, res, ok, Nint)
res = det
if(p1 /= 0) then
ii = (p1-1)/bit_kind_size + 1
pos = mod(p1-1, 64)!iand(p1-1,bit_kind_size-1)
if(iand(det(ii, s1), ishft(1_bit_kind, pos)) /= 0_8) return
res(ii, s1) = ibset(res(ii, s1), pos)
ii = (p1-1)/bit_kind_size + 1
pos = iand(p1-1,bit_kind_size-1)
if(iand(det(ii, s1), ishft(1_bit_kind, pos)) /= 0_8) then
return
endif
res(ii, s1) = ibset(res(ii, s1), pos)
end if
ii = (p2-1)/bit_kind_size + 1
pos = mod(p2-1, 64)!iand(p2-1,bit_kind_size-1)
if(iand(det(ii, s2), ishft(1_bit_kind, pos)) /= 0_8) return
res(ii, s2) = ibset(res(ii, s2), pos)
ii = (p2-1)/bit_kind_size + 1
pos = iand(p2-1,bit_kind_size-1)
if(iand(det(ii, s2), ishft(1_bit_kind, pos)) /= 0_8) then
return
endif
res(ii, s2) = ibset(res(ii, s2), pos)
ok = .true.

2
plugins/Psiref_CAS/overwrite_with_cas.irp.f
View File

@ -1,3 +1,5 @@
program overwrite_w_cas
read_wf = .True.
TOUCH read_wf
call extract_ref
end

2
plugins/QmcChem/dressed_dmc.irp.f
View File

@ -57,7 +57,7 @@ program dressed_dmc
enddo
call davidson_diag_hjj(psi_det,psi_coef,H_jj,energies,size(psi_coef,1),N_det,N_states,N_states_diag,,N_int,6)
call davidson_diag_hjj(psi_det,psi_coef,H_jj,energies,size(psi_coef,1),N_det,N_states,N_states_diag,N_int,6)
call save_wavefunction
call write_spindeterminants

4
plugins/Selectors_full/zmq.irp.f
View File

@ -101,6 +101,7 @@ subroutine zmq_get_psi(zmq_to_qp_run_socket, worker_id, energy, size_energy)
print *, '77_zmq_recv(zmq_to_qp_run_socket,psi_coef,psi_det_size*N_states*8,ZMQ_SNDMORE)'
stop 'error'
endif
TOUCH psi_det psi_coef
rc = f77_zmq_recv(zmq_to_qp_run_socket,energy,size_energy*8,0)
if (rc /= size_energy*8) then
@ -110,11 +111,12 @@ subroutine zmq_get_psi(zmq_to_qp_run_socket, worker_id, energy, size_energy)
if (N_det_generators_read > 0) then
N_det_generators = N_det_generators_read
TOUCH N_det_generators
endif
if (N_det_selectors_read > 0) then
N_det_selectors = N_det_selectors_read
TOUCH N_det_selectors
endif
SOFT_TOUCH psi_det psi_coef N_det_selectors N_det_generators
end

28
plugins/mrcepa0/EZFIO.cfg
View File

@ -3,3 +3,31 @@ type: Positive_int
doc: lambda type
interface: ezfio,provider,ocaml
default: 0
[energy]
type: double precision
doc: Calculated energy
interface: ezfio
[energy_pt2]
type: double precision
doc: Calculated energy with PT2 contribution
interface: ezfio
[energy]
type: double precision
doc: Calculated energy
interface: ezfio
[thresh_dressed_ci]
type: Threshold
doc: Threshold on the convergence of the dressed CI energy
interface: ezfio,provider,ocaml
default: 1.e-4
[n_it_max_dressed_ci]
type: Strictly_positive_int
doc: Maximum number of dressed CI iterations
interface: ezfio,provider,ocaml
default: 10

12
plugins/mrcepa0/mrcepa0_general.irp.f
View File

@ -17,8 +17,8 @@ subroutine run(N_st,energy)
thresh_mrcc = 1d-7
n_it_mrcc_max = 10
thresh_mrcc = thresh_dressed_ci
n_it_mrcc_max = n_it_max_dressed_ci
if(n_it_mrcc_max == 1) then
do j=1,N_states_diag
@ -48,8 +48,8 @@ subroutine run(N_st,energy)
E_new = sum(ci_energy_dressed)
delta_E = dabs(E_new - E_old)
call save_wavefunction
call ezfio_set_mrcc_cassd_energy(ci_energy_dressed(1))
if (iteration > n_it_mrcc_max) then
call ezfio_set_mrcepa0_energy(ci_energy_dressed(1))
if (iteration >= n_it_mrcc_max) then
exit
endif
enddo
@ -184,7 +184,7 @@ subroutine run_pt2_old(N_st,energy)
print *, '-----'
! call ezfio_set_full_ci_energy_pt2(energy+pt2)
call ezfio_set_mrcepa0_energy_pt2(energy(1)+pt2(1))
end
@ -238,5 +238,7 @@ subroutine run_pt2(N_st,energy)
print *, 'E+PT2 = ', energy+pt2
print *, '-----'
call ezfio_set_mrcepa0_energy_pt2(energy(1)+pt2(1))
end

227
src/Davidson/davidson_parallel.irp.f
View File

@ -4,12 +4,12 @@
use bitmasks
use f77_zmq
subroutine davidson_process(blockb, blocke, N, idx, vt, st, bs)
subroutine davidson_process(blockb, blockb2, N, idx, vt, st, bs, istep)
implicit none
integer , intent(in) :: blockb, blocke, bs
integer , intent(in) :: blockb, bs, blockb2, istep
integer , intent(inout) :: N
integer , intent(inout) :: idx(bs)
double precision , intent(inout) :: vt(N_states_diag, bs)
@ -20,11 +20,12 @@ subroutine davidson_process(blockb, blocke, N, idx, vt, st, bs)
double precision :: s2, hij
logical, allocatable :: wrotten(:)
allocate(wrotten(bs))
wrotten = .false.
PROVIDE dav_det
do sh = blockb, blocke
ii=0
sh = blockb
do sh2=1,shortcut_(0,1)
exa = 0
do ni=1,N_int
@ -32,7 +33,7 @@ subroutine davidson_process(blockb, blocke, N, idx, vt, st, bs)
end do
if(exa > 2) cycle
do i=shortcut_(sh,1),shortcut_(sh+1,1)-1
do i=blockb2+shortcut_(sh,1),shortcut_(sh+1,1)-1, istep
ii = i - shortcut_(blockb,1) + 1
org_i = sort_idx_(i,1)
@ -48,8 +49,8 @@ subroutine davidson_process(blockb, blocke, N, idx, vt, st, bs)
ext = ext + popcnt(xor(sorted_i(ni), sorted_(ni,j,1)))
end do
if(ext <= 4) then
call i_h_j (dav_det(1,1,org_j),dav_det(1,1,org_i),n_int,hij)
call get_s2(dav_det(1,1,org_j),dav_det(1,1,org_i),n_int,s2)
call i_h_j (dav_det(1,1,org_j),dav_det(1,1,org_i),n_int,hij)
if(.not. wrotten(ii)) then
wrotten(ii) = .true.
idx(ii) = org_i
@ -64,39 +65,39 @@ subroutine davidson_process(blockb, blocke, N, idx, vt, st, bs)
enddo
enddo
enddo
enddo
do sh=blockb,min(blocke, shortcut_(0,2))
do sh2=sh, shortcut_(0,2), shortcut_(0,1)*51
do i=shortcut_(sh2,2),shortcut_(sh2+1,2)-1
ii += 1
org_i = sort_idx_(i,2)
do j=shortcut_(sh2,2),shortcut_(sh2+1,2)-1
if(i == j) cycle
org_j = sort_idx_(j,2)
ext = 0
do ni=1,N_int
ext = ext + popcnt(xor(sorted_(ni,i,2), sorted_(ni,j,2)))
if (blockb <= shortcut_(0,2)) then
sh=blockb
do sh2=sh, shortcut_(0,2), shortcut_(0,1)
do i=blockb2+shortcut_(sh2,2),shortcut_(sh2+1,2)-1, istep
ii += 1
org_i = sort_idx_(i,2)
do j=shortcut_(sh2,2),shortcut_(sh2+1,2)-1
if(i == j) cycle
org_j = sort_idx_(j,2)
ext = 0
do ni=1,N_int
ext = ext + popcnt(xor(sorted_(ni,i,2), sorted_(ni,j,2)))
end do
if(ext == 4) then
call i_h_j (dav_det(1,1,org_j),dav_det(1,1,org_i),n_int,hij)
call get_s2(dav_det(1,1,org_j),dav_det(1,1,org_i),n_int,s2)
if(.not. wrotten(ii)) then
wrotten(ii) = .true.
idx(ii) = org_i
vt (:,ii) = 0d0
st (:,ii) = 0d0
end if
do istate=1,N_states_diag
vt (istate,ii) += hij*dav_ut(istate,org_j)
st (istate,ii) += s2*dav_ut(istate,org_j)
enddo
end if
end do
if(ext == 4) then
call i_h_j (dav_det(1,1,org_j),dav_det(1,1,org_i),n_int,hij)
call get_s2(dav_det(1,1,org_j),dav_det(1,1,org_i),n_int,s2)
if(.not. wrotten(ii)) then
wrotten(ii) = .true.
idx(ii) = org_i
vt (:,ii) = 0d0
st (:,ii) = 0d0
end if
do istate=1,N_states_diag
vt (istate,ii) += hij*dav_ut(istate,org_j)
st (istate,ii) += s2*dav_ut(istate,org_j)
enddo
end if
end do
end do
enddo
enddo
enddo
endif
N=0
do i=1,bs
@ -107,16 +108,17 @@ subroutine davidson_process(blockb, blocke, N, idx, vt, st, bs)
st(:,N) = st(:,i)
end if
end do
end subroutine
subroutine davidson_collect(blockb, blocke, N, idx, vt, st , v0t, s0t)
subroutine davidson_collect(N, idx, vt, st , v0t, s0t)
implicit none
integer , intent(in) :: blockb, blocke
integer , intent(in) :: N
integer , intent(in) :: idx(N)
double precision , intent(in) :: vt(N_states_diag, N)
@ -138,28 +140,48 @@ subroutine davidson_collect(blockb, blocke, N, idx, vt, st , v0t, s0t)
end subroutine
subroutine davidson_init(zmq_to_qp_run_socket)
subroutine davidson_init(zmq_to_qp_run_socket,n,n_st_8,ut)
use f77_zmq
implicit none
integer(ZMQ_PTR), intent(out) :: zmq_to_qp_run_socket
integer(ZMQ_PTR), intent(out) :: zmq_to_qp_run_socket
integer, intent(in) :: n, n_st_8
double precision, intent(in) :: ut(n_st_8,n)
integer :: i,k
touch dav_size dav_det dav_ut
dav_size = n
touch dav_size
do i=1,n
do k=1,N_int
dav_det(k,1,i) = psi_det(k,1,i)
dav_det(k,2,i) = psi_det(k,2,i)
enddo
enddo
do i=1,n
do k=1,N_states_diag
dav_ut(k,i) = ut(k,i)
enddo
enddo
touch dav_det dav_ut
call new_parallel_job(zmq_to_qp_run_socket,"davidson")
end subroutine
subroutine davidson_add_task(zmq_to_qp_run_socket, blockb, blocke)
subroutine davidson_add_task(zmq_to_qp_run_socket, blockb, blockb2, istep)
use f77_zmq
implicit none
integer(ZMQ_PTR) ,intent(in) :: zmq_to_qp_run_socket
integer ,intent(in) :: blockb, blocke
integer ,intent(in) :: blockb, blockb2, istep
character*(512) :: task
write(task,*) blockb, blocke
write(task,*) blockb, blockb2, istep
call add_task_to_taskserver(zmq_to_qp_run_socket, task)
end subroutine
@ -188,7 +210,7 @@ subroutine davidson_run_slave(thread,iproc)
integer, intent(in) :: thread, iproc
integer :: worker_id, task_id, blockb, blocke
integer :: worker_id, task_id, blockb
character*(512) :: task
integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket
@ -228,7 +250,7 @@ subroutine davidson_slave_work(zmq_to_qp_run_socket, zmq_socket_push, worker_id)
character*(512) :: task
integer :: blockb, blocke
integer :: blockb, blockb2, istep
integer :: N
integer , allocatable :: idx(:)
double precision , allocatable :: vt(:,:)
@ -241,10 +263,10 @@ subroutine davidson_slave_work(zmq_to_qp_run_socket, zmq_socket_push, worker_id)
do
call get_task_from_taskserver(zmq_to_qp_run_socket,worker_id, task_id, task)
if(task_id == 0) exit
read (task,*) blockb, blocke
bs = shortcut_(blocke+1,1) - shortcut_(blockb, 1)
read (task,*) blockb, blockb2, istep
bs = shortcut_(blockb+1,1) - shortcut_(blockb, 1)
do i=blockb, shortcut_(0,2), shortcut_(0,1)
do j=i, min(i+blocke-blockb, shortcut_(0,2))
do j=i, min(i, shortcut_(0,2))
bs += shortcut_(j+1,2) - shortcut_(j, 2)
end do
end do
@ -255,11 +277,11 @@ subroutine davidson_slave_work(zmq_to_qp_run_socket, zmq_socket_push, worker_id)
allocate(st(N_states_diag, bs))
end if
call davidson_process(blockb, blocke, N, idx, vt, st, bs)
call davidson_process(blockb, blockb2, N, idx, vt, st, bs, istep)
call task_done_to_taskserver(zmq_to_qp_run_socket,worker_id,task_id)
call davidson_push_results(zmq_socket_push, blockb, blocke, N, idx, vt, st, task_id)
call davidson_push_results(zmq_socket_push, blockb, blockb2, N, idx, vt, st, task_id)
end do
end subroutine
@ -340,28 +362,28 @@ end subroutine
subroutine davidson_collector(zmq_to_qp_run_socket, zmq_socket_pull , v0, s0)
subroutine davidson_collector(zmq_to_qp_run_socket, zmq_socket_pull , v0, s0, LDA)
use f77_zmq
implicit none
integer(ZMQ_PTR), intent(in) :: zmq_to_qp_run_socket
integer(ZMQ_PTR), intent(in) :: zmq_socket_pull
double precision ,intent(inout) :: v0(dav_size, N_states_diag)
double precision ,intent(inout) :: s0(dav_size, N_states_diag)
integer :: more, task_id
integer :: LDA
integer(ZMQ_PTR), intent(in) :: zmq_to_qp_run_socket
integer(ZMQ_PTR), intent(in) :: zmq_socket_pull
integer :: blockb, blocke
integer :: N
integer , allocatable :: idx(:)
double precision , allocatable :: vt(:,:), v0t(:,:), s0t(:,:)
double precision , allocatable :: st(:,:)
double precision ,intent(inout) :: v0(LDA, N_states_diag)
double precision ,intent(inout) :: s0(LDA, N_states_diag)
integer :: more, task_id, taskn
integer :: blockb, blocke
integer :: N
integer , allocatable :: idx(:)
double precision , allocatable :: vt(:,:), v0t(:,:), s0t(:,:)
double precision , allocatable :: st(:,:)
integer :: msize
msize = (max_workload + max_blocksize)*2
msize = (1 + max_blocksize)*2
allocate(idx(msize))
allocate(vt(N_states_diag, msize))
allocate(st(N_states_diag, msize))
@ -376,32 +398,41 @@ subroutine davidson_collector(zmq_to_qp_run_socket, zmq_socket_pull , v0, s0)
do while (more == 1)
call davidson_pull_results(zmq_socket_pull, blockb, blocke, N, idx, vt, st, task_id)
!DIR$ FORCEINLINE
call davidson_collect(blockb, blocke, N, idx, vt, st , v0t, s0t)
call davidson_collect(N, idx, vt, st , v0t, s0t)
call zmq_delete_task(zmq_to_qp_run_socket,zmq_socket_pull,task_id,more)
end do
deallocate(idx,vt,st)
call dtranspose(v0t,size(v0t,1), v0, size(v0,1), N_states_diag, dav_size)
call dtranspose(s0t,size(s0t,1), s0, size(s0,1), N_states_diag, dav_size)
integer :: i,j
!DIR$ IVDEP
do j=1,N_states_diag
!DIR$ IVDEP
do i=1,dav_size
v0(i,j) = v0t(j,i)
s0(i,j) = s0t(j,i)
enddo
enddo
deallocate(v0t,s0t)
end subroutine
subroutine davidson_run(zmq_to_qp_run_socket , v0, s0)
subroutine davidson_run(zmq_to_qp_run_socket , v0, s0, LDA)
use f77_zmq
implicit none
integer(ZMQ_PTR), intent(in) :: zmq_to_qp_run_socket
integer :: LDA
integer(ZMQ_PTR), intent(in) :: zmq_to_qp_run_socket
integer(ZMQ_PTR),external :: new_zmq_to_qp_run_socket
integer(ZMQ_PTR) :: zmq_collector
integer(ZMQ_PTR) :: zmq_collector
integer(ZMQ_PTR), external :: new_zmq_pull_socket
integer(ZMQ_PTR) :: zmq_socket_pull
integer :: i
integer, external :: omp_get_thread_num
double precision , intent(inout) :: v0(dav_size, N_states_diag)
double precision , intent(inout) :: s0(dav_size, N_states_diag)
double precision , intent(inout) :: v0(LDA, N_states_diag)
double precision , intent(inout) :: s0(LDA, N_states_diag)
call zmq_set_running(zmq_to_qp_run_socket)
@ -412,19 +443,20 @@ subroutine davidson_run(zmq_to_qp_run_socket , v0, s0)
PROVIDE nproc
!$OMP PARALLEL DEFAULT(shared) private(i) num_threads(nproc+2)
i = omp_get_thread_num()
if (i==0) then
call davidson_collector(zmq_collector, zmq_socket_pull , v0, s0)
call end_zmq_to_qp_run_socket(zmq_collector)
call end_zmq_pull_socket(zmq_socket_pull)
call davidson_miniserver_end()
else if(i==1) then
call davidson_miniserver_run()
else
call davidson_slave_inproc(i)
endif
!$OMP PARALLEL NUM_THREADS(nproc+2) PRIVATE(i)
i = omp_get_thread_num()
if (i == 0 ) then
call davidson_collector(zmq_collector, zmq_socket_pull , v0, s0, LDA)
call end_zmq_to_qp_run_socket(zmq_collector)
call end_zmq_pull_socket(zmq_socket_pull)
call davidson_miniserver_end()
else if (i == 1 ) then
call davidson_miniserver_run ()
else
call davidson_slave_inproc(i)
endif
!$OMP END PARALLEL
call end_parallel_job(zmq_to_qp_run_socket, 'davidson')
end subroutine
@ -505,15 +537,28 @@ end subroutine
BEGIN_PROVIDER [ integer(bit_kind), dav_det, (N_int, 2, dav_size) ]
END_PROVIDER
BEGIN_PROVIDER [ double precision, dav_ut, (N_states_diag, dav_size) ]
BEGIN_PROVIDER [ integer(bit_kind), dav_det, (N_int, 2, dav_size) ]
&BEGIN_PROVIDER [ double precision, dav_ut, (N_states_diag, dav_size) ]
use bitmasks
implicit none
BEGIN_DOC
! Temporary arrays for parallel davidson
!
! Touched in davidson_miniserver_get
END_DOC
dav_det = 0_bit_kind
dav_ut = -huge(1.d0)
END_PROVIDER
BEGIN_PROVIDER [ integer, dav_size ]
implicit none
BEGIN_DOC
! Size of the arrays for Davidson
!
! Touched in davidson_miniserver_get
END_DOC
dav_size = 1
END_PROVIDER

8
src/Davidson/davidson_slave.irp.f
View File

@ -26,10 +26,10 @@ program davidson_slave
print *, 'Davidson slave running'
! !$OMP PARALLEL PRIVATE(i)
!i = omp_get_thread_num()
call davidson_slave_tcp(0)
!!$OMP END PARALLEL
!$OMP PARALLEL PRIVATE(i)
i = omp_get_thread_num()
call davidson_slave_tcp(i)
!$OMP END PARALLEL
end do
end

46
src/Davidson/diagonalization_hs2.irp.f
View File

@ -1,4 +1,4 @@
subroutine davidson_diag_hs2(dets_in,u_in,dim_in,energies,sze,N_st,N_st_diag,Nint,iunit)
subroutine davidson_diag_hs2(dets_in,u_in,s2_out,dim_in,energies,sze,N_st,N_st_diag,Nint,iunit)
use bitmasks
implicit none
BEGIN_DOC
@ -22,7 +22,7 @@ subroutine davidson_diag_hs2(dets_in,u_in,dim_in,energies,sze,N_st,N_st_diag,Nin
integer, intent(in) :: dim_in, sze, N_st, N_st_diag, Nint, iunit
integer(bit_kind), intent(in) :: dets_in(Nint,2,sze)
double precision, intent(inout) :: u_in(dim_in,N_st_diag)
double precision, intent(out) :: energies(N_st)
double precision, intent(out) :: energies(N_st), s2_out(N_st)
double precision, allocatable :: H_jj(:), S2_jj(:)
double precision :: diag_h_mat_elem
@ -46,6 +46,9 @@ subroutine davidson_diag_hs2(dets_in,u_in,dim_in,energies,sze,N_st,N_st_diag,Nin
!$OMP END PARALLEL
call davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_st,N_st_diag,Nint,iunit)
do i=1,N_st_diag
s2_out(i) = S2_jj(i)
enddo
deallocate (H_jj,S2_jj)
end
@ -79,7 +82,8 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s
END_DOC
integer, intent(in) :: dim_in, sze, N_st, N_st_diag, Nint
integer(bit_kind), intent(in) :: dets_in(Nint,2,sze)
double precision, intent(in) :: H_jj(sze), S2_jj(sze)
double precision, intent(in) :: H_jj(sze)
double precision, intent(inout) :: S2_jj(sze)
integer, intent(in) :: iunit
double precision, intent(inout) :: u_in(dim_in,N_st_diag)
double precision, intent(out) :: energies(N_st_diag)
@ -200,7 +204,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s
call normalize(u_in(1,k),sze)
enddo
do while (.not.converged)
do k=1,N_st_diag
@ -239,12 +243,12 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s
! enddo
! enddo
call dgemm('T','N', shift2, N_st_diag, sze, &
1.d0, U, size(U,1), W(1,shift+1), size(W,1), &
call dgemm('T','N', shift2, N_st_diag, sze, &
1.d0, U, size(U,1), W(1,shift+1), size(W,1), &
0.d0, h(1,shift+1), size(h,1))
call dgemm('T','N', shift2, N_st_diag, sze, &
1.d0, U, size(U,1), S(1,shift+1), size(S,1), &
call dgemm('T','N', shift2, N_st_diag, sze, &
1.d0, U, size(U,1), S(1,shift+1), size(S,1), &
0.d0, s_(1,shift+1), size(s_,1))
! Diagonalize h
@ -254,14 +258,14 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s
! Compute S2 for each eigenvector
! -------------------------------
call dgemm('N','N',shift2,shift2,shift2, &
1.d0, s_, size(s_,1), y, size(y,1), &
call dgemm('N','N',shift2,shift2,shift2, &
1.d0, s_, size(s_,1), y, size(y,1), &
0.d0, s_tmp, size(s_tmp,1))
call dgemm('T','N',shift2,shift2,shift2, &
1.d0, y, size(y,1), s_tmp, size(s_tmp,1), &
call dgemm('T','N',shift2,shift2,shift2, &
1.d0, y, size(y,1), s_tmp, size(s_tmp,1), &
0.d0, s_, size(s_,1))
do k=1,shift2
s2(k) = s_(k,k) + S_z2_Sz
enddo
@ -326,8 +330,8 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s
! enddo
do k=1,N_st_diag
do i=1,sze
R(i,k) = (lambda(k) * U(i,shift2+k) - W(i,shift2+k) ) &
* (1.d0 + s2(k) * U(i,shift2+k) - S(i,shift2+k) - S_z2_Sz)
R(i,k) = (lambda(k) * U(i,shift2+k) - W(i,shift2+k) ) &
* (1.d0 + s2(k) * U(i,shift2+k) - S(i,shift2+k) - S_z2_Sz)
enddo
if (k <= N_st) then
residual_norm(k) = u_dot_u(R(1,k),sze)
@ -367,10 +371,10 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s
! enddo
! enddo
!
call dgemv('T',sze,N_st_diag*iter,1.d0,U,size(U,1), &
U(1,shift2+k),1,0.d0,c,1)
call dgemv('N',sze,N_st_diag*iter,-1.d0,U,size(U,1), &
c,1,1.d0,U(1,shift2+k),1)
call dgemv('T',sze,N_st_diag*iter,1.d0,U,size(U,1), &
U(1,shift2+k),1,0.d0,c,1)
call dgemv('N',sze,N_st_diag*iter,-1.d0,U,size(U,1), &
c,1,1.d0,U(1,shift2+k),1)
!
! do l=1,k-1
! c(1) = u_dot_v(U(1,shift2+k),U(1,shift2+l),sze)
@ -379,7 +383,7 @@ subroutine davidson_diag_hjj_sjj(dets_in,u_in,H_jj,S2_jj,energies,dim_in,sze,N_s
! enddo
! enddo
!
call dgemv('T',sze,k-1,1.d0,U(1,shift2+1),size(U,1), &
call dgemv('T',sze,k-1,1.d0,U(1,shift2+1),size(U,1), &
U(1,shift2+k),1,0.d0,c,1)
call dgemv('N',sze,k-1,-1.d0,U(1,shift2+1),size(U,1), &
c,1,1.d0,U(1,shift2+k),1)

8
src/Davidson/diagonalize_CI.irp.f
View File

@ -59,13 +59,13 @@ END_PROVIDER
! size(CI_eigenvectors,1), &
! N_det,min(N_det,N_states),min(N_det,N_states_diag),N_int,output_determinants)
!
call davidson_diag_HS2(psi_det,CI_eigenvectors, &
! call u_0_S2_u_0(CI_eigenvectors_s2,CI_eigenvectors,N_det,psi_det,N_int,&
! min(N_det,N_states_diag),size(CI_eigenvectors,1))
call davidson_diag_HS2(psi_det,CI_eigenvectors, CI_eigenvectors_s2, &
size(CI_eigenvectors,1),CI_electronic_energy, &
N_det,min(N_det,N_states),min(N_det,N_states_diag),N_int,output_determinants)
call u_0_S2_u_0(CI_eigenvectors_s2,CI_eigenvectors,N_det,psi_det,N_int,&
min(N_det,N_states_diag),size(CI_eigenvectors,1))
else if (diag_algorithm == "Lapack") then

62
src/Davidson/diagonalize_CI_SC2.irp.f
View File

@ -1,62 +0,0 @@
BEGIN_PROVIDER [ double precision, CI_SC2_energy, (N_states_diag) ]
implicit none
BEGIN_DOC
! N_states_diag lowest eigenvalues of the CI matrix
END_DOC
integer :: j
character*(8) :: st
call write_time(output_determinants)
do j=1,N_states_diag
CI_SC2_energy(j) = CI_SC2_electronic_energy(j) + nuclear_repulsion
write(st,'(I4)') j
call write_double(output_determinants,CI_SC2_energy(j),'Energy of state '//trim(st))
enddo
END_PROVIDER
BEGIN_PROVIDER [ double precision, threshold_convergence_SC2]
implicit none
BEGIN_DOC
! convergence of the correlation energy of SC2 iterations
END_DOC
threshold_convergence_SC2 = 1.d-10
END_PROVIDER
BEGIN_PROVIDER [ double precision, CI_SC2_electronic_energy, (N_states_diag) ]
&BEGIN_PROVIDER [ double precision, CI_SC2_eigenvectors, (N_det,N_states_diag) ]
&BEGIN_PROVIDER [ double precision, Diag_H_elements_SC2, (N_det) ]
implicit none
BEGIN_DOC
! Eigenvectors/values of the CI matrix
END_DOC
integer :: i,j
do j=1,N_states_diag
do i=1,N_det
CI_SC2_eigenvectors(i,j) = psi_coef(i,j)
enddo
CI_SC2_electronic_energy(j) = CI_electronic_energy(j)
enddo
call CISD_SC2(psi_det,CI_SC2_eigenvectors,CI_SC2_electronic_energy, &
! size(CI_SC2_eigenvectors,1),N_det,N_states_diag,N_int,threshold_convergence_SC2)
diag_H_elements_SC2,size(CI_SC2_eigenvectors,1),N_det,N_states_diag,N_int,threshold_convergence_SC2)
END_PROVIDER
subroutine diagonalize_CI_SC2
implicit none
BEGIN_DOC
! Replace the coefficients of the CI states_diag by the coefficients of the
! eigenstates of the CI matrix
END_DOC
integer :: i,j
do j=1,N_states_diag
do i=1,N_det
psi_coef(i,j) = CI_SC2_eigenvectors(i,j)
enddo
enddo
SOFT_TOUCH psi_coef CI_SC2_electronic_energy CI_SC2_energy CI_SC2_eigenvectors diag_h_elements_sc2
! SOFT_TOUCH psi_coef CI_SC2_electronic_energy CI_SC2_energy CI_SC2_eigenvectors
end

76
src/Davidson/u0Hu0.irp.f
View File

@ -24,6 +24,7 @@ subroutine u_0_H_u_0(e_0,u_0,n,keys_tmp,Nint,N_st,sze_8)
do i=1,N_st
e_0(i) = u_dot_v(v_0(1,i),u_0(1,i),n)/u_dot_u(u_0(1,i),n)
enddo
deallocate (H_jj, v_0)
end
@ -56,10 +57,9 @@ subroutine H_u_0_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,N_st,sze_8)
integer :: N_st_8
integer, external :: align_double
!!!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: vt, ut
!DIR$ ATTRIBUTES ALIGN : $IRP_ALIGN :: vt, ut
if(N_st /= N_states_diag) stop "H_u_0_nstates N_st /= N_states_diag"
N_st_8 = N_states_diag ! align_double(N_st)
N_st_8 = align_double(N_st)
ASSERT (Nint > 0)
ASSERT (Nint == N_int)
@ -165,7 +165,7 @@ subroutine H_u_0_nstates(v_0,u_0,H_jj,n,keys_tmp,Nint,N_st,sze_8)
v_0(i,istate) += H_jj(i) * u_0(i,istate)
enddo
enddo
!deallocate (shortcut, sort_idx, sorted, version, ut)
deallocate (shortcut, sort_idx, sorted, version, ut)
end
BEGIN_PROVIDER [ double precision, psi_energy, (N_states) ]
@ -200,76 +200,82 @@ subroutine H_S2_u_0_nstates(v_0,s_0,u_0,H_jj,S2_jj,n,keys_tmp,Nint,N_st,sze_8)
integer :: i,j,k,l, jj,ii
integer :: i0, j0
integer, allocatable :: shortcut(:,:), sort_idx(:,:)
integer(bit_kind), allocatable :: sorted(:,:,:), version(:,:,:)
integer, allocatable :: shortcut(:,:), sort_idx(:)
integer(bit_kind), allocatable :: sorted(:,:), version(:,:)
integer(bit_kind) :: sorted_i(Nint)
integer :: sh, sh2, ni, exa, ext, org_i, org_j, endi, istate
integer :: N_st_8
integer, external :: align_double
integer :: workload, blockb, blocke
integer :: blockb, blockb2, istep
double precision :: ave_workload, workload
integer(ZMQ_PTR) :: handler
if(N_st /= N_states_diag .or. sze_8 < N_det) stop "assert fail in H_S2_u_0_nstates"
N_st_8 = N_st !! align_double(N_st)
N_st_8 = N_st ! align_double(N_st)
ASSERT (Nint > 0)
ASSERT (Nint == N_int)
ASSERT (n>0)
PROVIDE ref_bitmask_energy
allocate (shortcut(0:n+1,2), sort_idx(n,2), sorted(Nint,n,2), version(Nint,n,2))
allocate (shortcut(0:n+1,2), sort_idx(n), sorted(Nint,n), version(Nint,n))
allocate(ut(N_st_8,n))
v_0 = 0.d0
s_0 = 0.d0
if(n /= N_det) stop "n /= N_det"
do i=1,n
do istate=1,N_st
ut(istate,i) = u_0(i,istate)
enddo
enddo
call sort_dets_ab_v(keys_tmp, sorted(1,1,1), sort_idx(1,1), shortcut(0,1), version(1,1,1), n, Nint)
call sort_dets_ba_v(keys_tmp, sorted(1,1,2), sort_idx(1,2), shortcut(0,2), version(1,1,2), n, Nint)
call sort_dets_ab_v(keys_tmp, sorted, sort_idx, shortcut(0,1), version, n, Nint)
call sort_dets_ba_v(keys_tmp, sorted, sort_idx, shortcut(0,2), version, n, Nint)
dav_size = n
touch dav_size
dav_det = psi_det
dav_ut = ut
workload = 0
blockb = shortcut(0,1)
blocke = blockb
call davidson_init(handler)
call davidson_init(handler,n,N_st_8,ut)
ave_workload = 0.d0
do sh=1,shortcut(0,1)
ave_workload += shortcut(0,1)
ave_workload += (shortcut(sh+1,1) - shortcut(sh,1))**2
do i=sh, shortcut(0,2), shortcut(0,1)
do j=i, min(i, shortcut(0,2))
ave_workload += (shortcut(j+1,2) - shortcut(j, 2))**2
end do
end do
enddo
ave_workload = ave_workload/dble(shortcut(0,1))
do sh=shortcut(0,1),1,-1
workload += (shortcut(sh+1,1) - shortcut(sh,1))**2
if(workload > max_workload) then
blocke = sh
call davidson_add_task(handler, blocke, blockb)
blockb = sh-1
workload = 0
end if
workload = shortcut(0,1)+dble(shortcut(sh+1,1) - shortcut(sh,1))**2
do i=sh, shortcut(0,2), shortcut(0,1)
do j=i, min(i, shortcut(0,2))
workload += (shortcut(j+1,2) - shortcut(j, 2))**2
end do
end do
istep = 1+ int(0.5d0*workload/ave_workload)
do blockb2=0, istep-1
call davidson_add_task(handler, sh, blockb2, istep)
enddo
enddo
if(blockb > 0) call davidson_add_task(handler, 1, blockb)
call davidson_run(handler, v_0, s_0)
call davidson_run(handler, v_0, s_0, size(v_0,1))
do istate=1,N_st
do i=1,n
v_0(i,istate) = v_0(i,istate) + H_jj(i) * u_0(i,istate)
s_0(i,istate) = s_0(i,istate) + s2_jj(i)* u_0(i,istate)
enddo
enddo
deallocate(shortcut, sort_idx, sorted, version)
deallocate(ut)
end
BEGIN_PROVIDER [ integer, max_workload ]
max_workload = 1000
END_PROVIDER

6
src/Integrals_Bielec/map_integrals.irp.f
View File

@ -364,7 +364,7 @@ double precision function get_mo_bielec_integral(i,j,k,l,map)
integer(key_kind) :: idx
type(map_type), intent(inout) :: map
real(integral_kind) :: tmp
PROVIDE mo_bielec_integrals_in_map
PROVIDE mo_bielec_integrals_in_map mo_integrals_cache
if ( (i >= mo_integrals_cache_min) .and. &
(j >= mo_integrals_cache_min) .and. &
(k >= mo_integrals_cache_min) .and. &
@ -393,7 +393,7 @@ double precision function get_mo_bielec_integral_schwartz(i,j,k,l,map)
integer(key_kind) :: idx
type(map_type), intent(inout) :: map
real(integral_kind) :: tmp
PROVIDE mo_bielec_integrals_in_map
PROVIDE mo_bielec_integrals_in_map mo_integrals_cache
if (mo_bielec_integral_schwartz(i,k)*mo_bielec_integral_schwartz(j,l) > mo_integrals_threshold) then
double precision, external :: get_mo_bielec_integral
!DIR$ FORCEINLINE
@ -411,7 +411,7 @@ double precision function mo_bielec_integral(i,j,k,l)
END_DOC
integer, intent(in) :: i,j,k,l
double precision :: get_mo_bielec_integral
PROVIDE mo_bielec_integrals_in_map
PROVIDE mo_bielec_integrals_in_map mo_integrals_cache
!DIR$ FORCEINLINE
mo_bielec_integral = get_mo_bielec_integral(i,j,k,l,mo_integrals_map)
return

17
tests/bats/cassd.bats Normal file
View File

@ -0,0 +1,17 @@
#!/usr/bin/env bats
source $QP_ROOT/tests/bats/common.bats.sh
@test "CAS_SD H2O cc-pVDZ" {
test_exe cas_sd_selected || skip
INPUT=h2o.ezfio
qp_edit -c $INPUT
ezfio set_file $INPUT
ezfio set perturbation do_pt2_end False
ezfio set determinants n_det_max 1000
qp_set_mo_class $INPUT -core "[1]" -inact "[2,5]" -act "[3,4,6,7]" -virt "[8-24]"
qp_run cas_sd_selected $INPUT
energy="$(ezfio get cas_sd energy)"
eq $energy -76.22213389282479 1.E-5
}

44
tests/bats/common.bats.sh Normal file
View File

@ -0,0 +1,44 @@
#!/usr/bin/env bats
# floating point number comparison
# Compare two numbers ($1, $2) with a given precision ($3)
# If the numbers are not equal, the exit code is 1 else it is 0
# So we strip the "-", is the abs value of the poor
function eq() {
declare -a diff
diff=($(awk -v d1=$1 -v d2=$2 -v n1=${1#-} -v n2=${2#-} -v p=$3 'BEGIN{ if ((n1-n2)^2 < p^2) print 0; print 1 " " (d1-d2) " " d1 " " d2 }'))
if [[ "${diff[0]}" == "0" ]]
then
return 0
else
echo "Test : " ${BATS_TEST_DESCRIPTION}
echo "Error : " ${diff[1]}
echo "Reference : " ${diff[3]}
echo "Computed : " ${diff[2]}
exit 1
fi
}
# ___
# | ._ o _|_
# _|_ | | | |_
#
source ${QP_ROOT}/install/EZFIO/Bash/ezfio.sh
TEST_DIR=${QP_ROOT}/tests/work/
mkdir -p "${TEST_DIR}"
cd "${TEST_DIR}" || exit 1
function test_exe() {
l_EXE=$(awk "/^$1 / { print \$2 }" < "${QP_ROOT}"/data/executables)
l_EXE=$(echo $l_EXE | sed "s|\$QP_ROOT|$QP_ROOT|")
if [[ -x "$l_EXE" ]]
then
return 0
else
return 127
fi
}

27
tests/bats/convert.bats Normal file
View File

@ -0,0 +1,27 @@
#!/usr/bin/env bats
source $QP_ROOT/tests/bats/common.bats.sh
#=== Convert
@test "gamess convert HBO.out" {
cp ${QP_ROOT}/tests/input/HBO.out .
qp_convert_output_to_ezfio.py HBO.out
qp_edit -c HBO.out.ezfio
ezfio set_file HBO.out.ezfio
qp_run SCF HBO.out.ezfio
# Check energy
energy="$(ezfio get hartree_fock energy)"
eq $energy -100.0185822590964 1.e-10
}
@test "g09 convert H2O.log" {
cp ${QP_ROOT}/tests/input/h2o.log .
qp_convert_output_to_ezfio.py h2o.log
qp_edit -c h2o.log.ezfio
ezfio set_file h2o.log.ezfio
qp_run SCF h2o.log.ezfio
# Check energy
energy="$(ezfio get hartree_fock energy)"
eq $energy -76.0270218704265 1E-10
}

52
tests/bats/fci.bats Normal file
View File

@ -0,0 +1,52 @@
#!/usr/bin/env bats
source $QP_ROOT/tests/bats/common.bats.sh
function run_FCI() {
thresh=5.e-5
test_exe full_ci || skip
qp_edit -c $1
ezfio set_file $1
ezfio set perturbation do_pt2_end True
ezfio set determinants n_det_max $2
ezfio set davidson threshold_davidson 1.e-10
qp_run full_ci $1
energy="$(ezfio get full_ci energy)"
eq $energy $3 $thresh
energy_pt2="$(ezfio get full_ci energy_pt2)"
eq $energy_pt2 $4 $thresh
}
function run_FCI_ZMQ() {
thresh=5.e-5
test_exe full_ci || skip
qp_edit -c $1
ezfio set_file $1
ezfio set perturbation do_pt2_end True
ezfio set determinants n_det_max $2
ezfio set davidson threshold_davidson 1.e-10
qp_run fci_zmq $1
energy="$(ezfio get full_ci energy)"
eq $energy $3 $thresh
energy_pt2="$(ezfio get full_ci energy_pt2)"
eq $energy_pt2 $4 $thresh
}
#=== H2O
@test "qp_set_mo_class H2O cc-pVDZ" {
qp_set_mo_class h2o.ezfio -core "[1]" -act "[2-12]" -del "[13-24]"
}
@test "FCI H2O cc-pVDZ" {
run_FCI h2o.ezfio 2000 -0.761255633582109E+02 -0.761258377850042E+02
}
@test "FCI-ZMQ H2O cc-pVDZ" {
run_FCI_ZMQ h2o.ezfio 2000 -0.761255633582109E+02 -0.761258377850042E+02
}

27
tests/bats/foboci.bats Normal file
View File

@ -0,0 +1,27 @@
#!/usr/bin/env bats
source $QP_ROOT/tests/bats/common.bats.sh
function run_all_1h_1p() {
thresh=1.e-6
test_exe all_1h_1p || skip
qp_edit -c $1
ezfio set_file $1
ezfio set determinants n_det_max $2
ezfio set perturbation pt2_max $3
ezfio set davidson threshold_davidson 1.e-10
qp_run all_1h_1p $1 | tee $1.F1h1p.out
energy="$(ezfio get all_singles energy)"
eq $energy $4 $thresh
}
#=== DHNO
@test "all_1h_1p DHNO chipman-dzp" {
qp_set_mo_class -inact "[1-8]" -act "[9]" -virt "[10-64]" dhno.ezfio
run_all_1h_1p dhno.ezfio 10000 0.0000000001 -130.4466283766202
}

52
tests/bats/hf.bats Normal file
View File

@ -0,0 +1,52 @@
#!/usr/bin/env bats
source $QP_ROOT/tests/bats/common.bats.sh
function run_init() {
cp "${QP_ROOT}/tests/input/$1" .
qp_create_ezfio_from_xyz $1 -o $3 $2
qp_edit -c $3
}
function run_HF() {
thresh=1.e-7
test_exe SCF || skip
qp_edit -c $1
ezfio set_file $1
ezfio set hartree_fock thresh_scf 1.e-11
qp_run SCF $1
energy="$(ezfio get hartree_fock energy)"
eq $energy $2 $thresh
}
#=== DHNO
@test "init DHNO chipman-dzp" {
run_init dhno.xyz "-b chipman-dzp -m 2" dhno.ezfio
}
@test "SCF DHNO chipman-dzp" {
run_HF dhno.ezfio -130.4278777822
}
#=== HBO
@test "init HBO STO-3G" {
run_init HBO.xyz "-b STO-3G" hbo.ezfio
}
@test "SCF HBO STO-3G" {
run_HF hbo.ezfio -98.8251985678084
}
#=== H2O
@test "init H2O cc-pVDZ" {
run_init h2o.xyz "-b cc-pvdz" h2o.ezfio
}
@test "SCF H2O cc-pVDZ" {
run_HF h2o.ezfio -0.760270218692179E+02
}

70
tests/bats/mrcepa0.bats Normal file
View File

@ -0,0 +1,70 @@
#!/usr/bin/env bats
source $QP_ROOT/tests/bats/common.bats.sh
#=== H2O
@test "MRCC-lambda H2O cc-pVDZ" {
INPUT=h2o.ezfio
EXE=mrcc
test_exe $EXE || skip
qp_edit -c $INPUT
ezfio set_file $INPUT
ezfio set determinants threshold_generators 1.
ezfio set determinants threshold_selectors 1.
ezfio set determinants read_wf True
ezfio set mrcepa0 lambda_type 1
ezfio set mrcepa0 n_it_max_dressed_ci 3
qp_run $EXE $INPUT
energy="$(ezfio get mrcepa0 energy)"
eq $energy -76.22903276183061 1.e-4
}
@test "MRCC H2O cc-pVDZ" {
INPUT=h2o.ezfio
EXE=mrcc
test_exe $EXE || skip
qp_edit -c $INPUT
ezfio set_file $INPUT
ezfio set determinants threshold_generators 1.
ezfio set determinants threshold_selectors 1.
ezfio set determinants read_wf True
ezfio set determinants read_wf True
ezfio set mrcepa0 lambda_type 0
ezfio set mrcepa0 n_it_max_dressed_ci 3
qp_run $EXE $INPUT
energy="$(ezfio get mrcepa0 energy)"
eq $energy -76.22899302846875 1.e-4
}
@test "MRSC2 H2O cc-pVDZ" {
INPUT=h2o.ezfio
EXE=mrsc2
test_exe $EXE || skip
qp_edit -c $INPUT
ezfio set_file $INPUT
ezfio set determinants threshold_generators 1.
ezfio set determinants threshold_selectors 1.
ezfio set determinants read_wf True
ezfio set mrcepa0 lambda_type 0
ezfio set mrcepa0 n_it_max_dressed_ci 3
qp_run $EXE $INPUT
energy="$(ezfio get mrcepa0 energy)"
eq $energy -76.22647345292708 1.e-4
}
@test "MRCEPA0 H2O cc-pVDZ" {
INPUT=h2o.ezfio
EXE=mrcepa0
test_exe $EXE || skip
qp_edit -c $INPUT
ezfio set_file $INPUT
ezfio set determinants threshold_generators 1.
ezfio set determinants threshold_selectors 1.
ezfio set determinants read_wf True
ezfio set mrcepa0 lambda_type 0
ezfio set mrcepa0 n_it_max_dressed_ci 3
qp_run $EXE $INPUT
energy="$(ezfio get mrcepa0 energy)"
eq $energy -76.23199784430074 1.e-4
}

53
tests/bats/pseudo.bats Normal file
View File

@ -0,0 +1,53 @@
#!/usr/bin/env bats
source $QP_ROOT/tests/bats/common.bats.sh
function run_init() {
cp "${QP_ROOT}/tests/input/$1" .
qp_create_ezfio_from_xyz $1 -o $3 $2
qp_edit -c $3
}
function run_HF() {
thresh=1.e-7
test_exe SCF || skip
qp_edit -c $1
ezfio set_file $1
ezfio set hartree_fock thresh_scf 1.e-11
qp_run SCF $1
energy="$(ezfio get hartree_fock energy)"
eq $energy $2 $thresh
}
function run_FCI_ZMQ() {
thresh=5.e-5
test_exe full_ci || skip
qp_edit -c $1
ezfio set_file $1
ezfio set perturbation do_pt2_end True
ezfio set determinants n_det_max $2
ezfio set davidson threshold_davidson 1.e-10
qp_run fci_zmq $1
energy="$(ezfio get full_ci energy)"
eq $energy $3 $thresh
energy_pt2="$(ezfio get full_ci energy_pt2)"
eq $energy_pt2 $4 $thresh
}
#=== H2O Pseudo
@test "init H2O VDZ pseudo" {
run_init h2o.xyz "-p bfd -b vdz-bfd" h2o_pseudo.ezfio
}
@test "SCF H2O VDZ pseudo" {
run_HF h2o_pseudo.ezfio -16.9483703905461
}
@test "FCI H2O VDZ pseudo" {
qp_set_mo_class h2o_pseudo.ezfio -core "[1]" -act "[2-12]" -del "[13-23]"
run_FCI_ZMQ h2o_pseudo.ezfio 2000 -0.170399597228904E+02 -0.170400168816800E+02
}

205
tests/bats/qp.bats
View File

@ -1,205 +0,0 @@
#!/usr/bin/env bats
#
# |\/| o _ _
# | | | _> (_
#
# floating point number comparison
# Compare two numbers ($1, $2) with a given precision ($3)
# If the numbers are not equal, the exit code is 1 else it is 0
# So we strip the "-", is the abs value of the poor
function eq() {
declare -a diff
diff=($(awk -v d1=$1 -v d2=$2 -v n1=${1#-} -v n2=${2#-} -v p=$3 'BEGIN{ if ((n1-n2)^2 < p^2) print 0; print 1 " " (d1-d2) " " d1 " " d2 }'))
if [[ "${diff[0]}" == "0" ]]
then
return 0
else
echo "Test : " ${BATS_TEST_DESCRIPTION}
echo "Error : " ${diff[1]}
echo "Reference : " ${diff[3]}
echo "Computed : " ${diff[2]}
exit 1
fi
}
# ___
# | ._ o _|_
# _|_ | | | |_
#
source ${QP_ROOT}/install/EZFIO/Bash/ezfio.sh
TEST_DIR=${QP_ROOT}/tests/work/
mkdir -p "${TEST_DIR}"
cd "${TEST_DIR}" || exit 1
function run_init() {
cp "${QP_ROOT}/tests/input/$1" .
qp_create_ezfio_from_xyz $1 -o $3 $2
qp_edit -c $3
}
function test_exe() {
EXE=$(awk "/^$1 / { print \$2 }" < "${QP_ROOT}"/data/executables)
EXE=$(echo $EXE | sed "s|\$QP_ROOT|$QP_ROOT|")
if [[ -x "$EXE" ]]
then
return 0
else
return 127
fi
}
function run_HF() {
thresh=1.e-7
test_exe SCF || skip
ezfio set_file $1
ezfio set hartree_fock thresh_scf 1.e-11
qp_run SCF $1
energy="$(ezfio get hartree_fock energy)"
eq $energy $2 $thresh
}
function run_FCI() {
thresh=5.e-5
test_exe full_ci || skip
ezfio set_file $1
ezfio set perturbation do_pt2_end True
ezfio set determinants n_det_max $2
ezfio set davidson threshold_davidson 1.e-10
qp_run full_ci $1
energy="$(ezfio get full_ci energy)"
eq $energy $3 $thresh
energy_pt2="$(ezfio get full_ci energy_pt2)"
eq $energy_pt2 $4 $thresh
}
function run_all_1h_1p() {
thresh=1.e-6
test_exe all_1h_1p || skip
ezfio set_file $1
ezfio set determinants n_det_max $2
ezfio set perturbation pt2_max $3
ezfio set davidson threshold_davidson 1.e-10
qp_run all_1h_1p $1 | tee $1.F1h1p.out
energy="$(ezfio get all_singles energy)"
eq $energy $4 $thresh
}
# ___
# | _ _ _|_
# | (/_ _> |_
#
#=== DHNO
@test "init DHNO chipman-dzp" {
run_init dhno.xyz "-b chipman-dzp -m 2" dhno.ezfio
}
@test "SCF DHNO chipman-dzp" {
run_HF dhno.ezfio -130.4278777822
}
@test "all_1h_1p DHNO chipman-dzp" {
qp_set_mo_class -inact "[1-8]" -act "[9]" -virt "[10-64]" dhno.ezfio
run_all_1h_1p dhno.ezfio 10000 0.0000000001 -130.4466283766202
}
#=== HBO
@test "init HBO STO-3G" {
run_init HBO.xyz "-b STO-3G" hbo.ezfio
}
@test "SCF HBO STO-3G" {
run_HF hbo.ezfio -98.8251985678084
}
#=== H2O
@test "init H2O cc-pVDZ" {
run_init h2o.xyz "-b cc-pvdz" h2o.ezfio
}
@test "SCF H2O cc-pVDZ" {
run_HF h2o.ezfio -0.760270218692179E+02
}
@test "FCI H2O cc-pVDZ" {
qp_set_mo_class h2o.ezfio -core "[1]" -act "[2-12]" -del "[13-24]"
run_FCI h2o.ezfio 2000 -0.761255633582109E+02 -0.761258377850042E+02
}
@test "CAS_SD H2O cc-pVDZ" {
test_exe cas_sd_selected || skip
INPUT=h2o.ezfio
ezfio set_file $INPUT
ezfio set perturbation do_pt2_end False
ezfio set determinants n_det_max 1000
qp_set_mo_class $INPUT -core "[1]" -inact "[2,5]" -act "[3,4,6,7]" -virt "[8-24]"
qp_run cas_sd_selected $INPUT
energy="$(ezfio get cas_sd energy)"
eq $energy -0.762219854008117E+02 1.E-5
}
@test "MRCC H2O cc-pVDZ" {
test_exe mrcc_cassd || skip
INPUT=h2o.ezfio
ezfio set_file $INPUT
ezfio set determinants threshold_generators 1.
ezfio set determinants threshold_selectors 1.
ezfio set determinants read_wf True
qp_run mrcc_cassd $INPUT
energy="$(ezfio get mrcc_cassd energy)"
eq $energy -76.2288648023833 1.e-4
}
#=== H2O Pseudo
@test "init H2O VDZ pseudo" {
run_init h2o.xyz "-p bfd -b vdz-bfd" h2o_pseudo.ezfio
}
@test "SCF H2O VDZ pseudo" {
run_HF h2o_pseudo.ezfio -16.9483703905461
}
@test "FCI H2O VDZ pseudo" {
qp_set_mo_class h2o_pseudo.ezfio -core "[1]" -act "[2-12]" -del "[13-23]"
run_FCI h2o_pseudo.ezfio 2000 -0.170399597228904E+02 -0.170400168816800E+02
}
#=== Convert
@test "gamess convert HBO.out" {
cp ${QP_ROOT}/tests/input/HBO.out .
qp_convert_output_to_ezfio.py HBO.out
ezfio set_file HBO.out.ezfio
qp_run SCF HBO.out.ezfio
# Check energy
energy="$(ezfio get hartree_fock energy)"
eq $energy -100.0185822590964 1.e-10
}
@test "g09 convert H2O.log" {
cp ${QP_ROOT}/tests/input/h2o.log .
qp_convert_output_to_ezfio.py h2o.log
ezfio set_file h2o.log.ezfio
qp_run SCF h2o.log.ezfio
# Check energy
energy="$(ezfio get hartree_fock energy)"
eq $energy -76.0270218704265 1E-10
}
# TODO N_int = 1,2,3,4,5
# TODO mod(64) MOs
# TODO All G2 SCF energies
# TODO Long and short tests
# TODO MP2
# TODO CISD_selected

4
tests/bats_to_sh.py
View File

@ -1,6 +1,8 @@
#!/usr/bin/env python
with open('bats/qp.bats','r') as f:
import sys
with open(sys.argv[1],'r') as f:
raw_data = f.read()
output = []

39
tests/run_tests.sh
View File

@ -1,18 +1,39 @@
#!/bin/bash
LIST="
convert.bats
hf.bats
foboci.bats
pseudo.bats
fci.bats
cassd.bats
mrcepa0.bats
"
export QP_PREFIX="timeout -s 9 300"
export QP_TASK_DEBUG=1
BATS_FILE=bats/qp.bats
rm -rf work output
if [[ "$1" == "-v" ]]
then
echo "Verbose mode"
./bats_to_sh.py $BATS_FILE | bash
else
bats $BATS_FILE
fi
for BATS_FILE in $LIST
do
echo
echo "-~-~-~-~-~-~"
echo
echo "Running tests for ${BATS_FILE%.bats}"
echo
BATS_FILE=bats/$BATS_FILE
if [[ "$1" == "-v" ]]
then
echo "Verbose mode"
./bats_to_sh.py $BATS_FILE | bash
else
bats $BATS_FILE
fi
done